Jeremy Rohmer
email: j.rohmer@brgm.fr

Olivier Roustant

Sophie Lecacheux

Jean-Charles Manceau

Revealing the dependence structure of scenario-like inputs in numerical environmental simulations using Gaussian Process regression

Keywords: Categorical variables, Computationally intensive simulator, Metamodel, Kriging, Model selection

come

Revealing the dependence structure of scenario-like inputs in numerical environmental simulations using Gaussian Process regression

Introduction

High-resolution numerical simulators are key ingredients of environmental science in order to get deeper insights into the natural systems' behavior. Some examples are [START_REF] Veeck | Scale dynamics of the HIDROPIXEL high-resolution DEM-based distributed hydrologic modeling approach[END_REF] for hydrologic modelling; [START_REF] Zhao | Large-scale, high-resolution agricultural systems modeling using a hybrid approach combining grid computing and parallel processing[END_REF] for agricultural modelling; [START_REF] Vandromme | ALICE (Assessment of Landslides Induced by Climatic Events): A single tool to integrate shallow and deep landslides for susceptibility and hazard assessment[END_REF] for landslide modelling; [START_REF] Abily | Spatial Global Sensitivity Analysis of High Resolution classified topographic data use in 2D urban flood modelling[END_REF] for urban flooding; [START_REF] Idier | Coastal flood: a composite method for past events characterisation providing insights in past, present and future hazards-joining historical, statistical and modelling approaches[END_REF] for marine flooding, etc. To model the natural system, these simulators all have in common to involve a large spectrum of assumptions related to the system's geometry, to the loading/forcing conditions acting on the system, to the system's intrinsic physical law, to the properties' values, etc. While most of these modelling assumptions can mathematically be represented by means of continuous variables (like geotechnical properties of a given soil formation, or time series of rainfall conditions at a given location, etc.), some of them involve scenario-like variables.

This type of variable is assigned to different modelling scenarios like the choice in the structure/form of the model (e.g. 1D versus 2D modelling, [START_REF] Leandro | Comparison of 1D/1D and 1D/2D coupled (sewer/surface) hydraulic models for urban flood simulation[END_REF], the selection of the physical processes regarded as "relevant and prominent" (e.g. account for spatial heterogeneity, [START_REF] Liu | Impact of surface-heterogeneity on atmosphere and land-surface interactions[END_REF], the use of alternative physical laws (e.g. different soil water retention curves, Silva [START_REF] Silva Ursulino | Modelling soil water dynamics from soil hydraulic parameters estimated by an alternative method in a tropical experimental basin[END_REF], the system's future evolution (e.g., future gas emission scenarios, Le [START_REF] Cozannet | Evaluating uncertainties of future marine flooding occurrence as sea-level rises[END_REF] or land use change, [START_REF] Mishra | Assessment of future flood inundations under climate and land use change scenarios in the Ciliwung River Basin, Jakarta[END_REF], etc.

Depending on the modelling scenarios, the simulation results can differ, hence resulting in uncertainty. This category of uncertainty can be termed as structural since it is associated to the structure/form of the model or to the unambiguous choice of the "best" model to be used: this type of uncertainty should be understood with respect to parametric uncertainties, which are related to the difficulties in estimating the model input parameters (in a broad sense) due to the limited number, poor representativeness (caused by time, space and financial limitations), and imprecision of observations/data (e.g., [START_REF] Hill | Risk assessment and uncertainty in natural hazards[END_REF].

A pragmatic approach for mathematically representing a scenario-like variable is via a categorical variable: a multi-level indicator that takes up a finite number of discrete values; each discrete level being associated to a different scenario (e.g. level a is associated to scenario a); see some real case applications in the domain of: safety analysis of radioactive waste disposal by [START_REF] Storlie | Analysis of computationally demanding models with continuous and categorical inputs[END_REF]; earthquake risk assessments by [START_REF] Rohmer | Weighing the importance of model uncertainty against parameter uncertainty in earthquake loss assessments[END_REF]; marine flooding induced by sea level rise by [START_REF] Cozannet | Evaluating uncertainties of future marine flooding occurrence as sea-level rises[END_REF]; reservoir engineering for CO2 geological storage by [START_REF] Manceau | Post-injection trapping of mobile CO 2 in deep aquifers: Assessing the importance of model and parameter uncertainties[END_REF]; pollution risk analysis and management by [START_REF] Lauvernet | Metamodeling methods that incorporate qualitative variables for improved design of vegetative filter strips[END_REF], etc. Characterizing and quantifying structural uncertainty by means of categorical multi-level variables raises however several practical questions regarding: (Q1) the dependence structure of the modelling scenarios: should each level (i.e. modelling scenario) of the considered categorical variable be treated equivalently with respect to the numerically simulated variable of interest? is there any dependence among the levels?, and can several levels be grouped?;

(Q2) the prediction using categorical variables: how to derive a predictive statistical model when only a few model runs (of the order of 100s) are available because of too high computation time cost? What is the performance of the derived model? (e.g. does it explain well the observations? what is its predictive capability?).

The objective of the present work is to explore how Gaussian Processes [START_REF] Williams | Gaussian processes for machine learning[END_REF], denoted GP, with mixed continuous/categorical inputs (e.g., [START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF][START_REF] Qian | Gaussian process models for computer experiments with qualitative and quantitative factors[END_REF][START_REF] Zhang | A latent variable approach to Gaussian process modeling with qualitative and quantitative factors[END_REF] can bring valuable insights into the interdependencies of the modelling assumptions, while fulfilling key characteristics like predictability (i.e. capability of the derived model of predicting "yet-unseen" input configurations) and explainability (i.e. capability of the derived model of representing the data). These characteristics are further detailed in Sect. 2.4. The advantage of the GP-based method lies in the flexibility brought by the use of covariance functions (also known as kernels) which specify how similar (i.e. correlated) two instances of the variable of interest, e.g., y and y0, are expected to be at two input values u and u0 (i.e., in our case, at two levels of the scenario-like input). This "similarity" function can be encoded in different manners depending on the assumptions regarding the categorical variable (nominal/ ordinal, inter-dependence between the levels, interactions between given levels, etc.); see for instance [START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF]; [START_REF] Lauvernet | Metamodeling methods that incorporate qualitative variables for improved design of vegetative filter strips[END_REF]. Once fitted, the resulting correlation matrix provides the mapping of the dependence structure that can be used to support (or contradict) some physical assumptions regarding the scenario-like input's influence.

To illustrate the type of results that can be derived, Figure 1(a) depicts an unknown relation between a continuous and a categorical input variable with 5 levels (i.e. five scenarios). Figure 1(b) depicts the correlation matrix derived from the GP-based analysis given 25 model runs: this summarizes the interplay between the levels; a group of highly correlated levels are identified for u1-4 as well as the decreasing correlation of u5 with the others (from 25 to 15% considering u4 to u1). These observations are consistent with the test function. If the functional relation in Figure 1(a) had been known, this result would have been straightforward, but here the structure is unknown and can be learnt only with a limited number of numerical results (here with only 25 model runs). This result (that is further discussed in Sect. 3.1) depends on how the GP kernel is defined, which raises the question of covariance kernel model selection that is addressed in the present study by relying on a multi-criterion analysis.

The present paper is organized as follows. Section 2 describes the different steps of the proposed procedure as well as the statistical methods. In this section, a multi-criterion approach for selecting the categorical covariance kernel model is detailed. Section 3 presents the application to the synthetic case (described in Figure 1) and to two real cases, namely for cyclone-induced wave numerical modelling [START_REF] Rohmer | Dynamic parameter sensitivity in numerical modelling of cyclone-induced waves: a multi-look approach using advanced meta-modelling techniques[END_REF], and for reservoir modelling of CO2 storage [START_REF] Manceau | Post-injection trapping of mobile CO 2 in deep aquifers: Assessing the importance of model and parameter uncertainties[END_REF]. The results are then discussed in Section 4 based on which practical recommendations are defined.

Methods

Description of the procedure

The proposed procedure holds as follows:

-Step 1: a series of random computer experiments are performed by running the expensive-to-evaluate environmental simulator considering a limited number of randomly selected input variables' configurations; -Step 2: using the set of random computer experiments (training dataset), different hypotheses regarding the structure of the considered input categorical variable are tested and modelled by means of different GP kernel (covariance) formulations (see further details in Sect. 2.2 and 2.3); -Step 3: the question of selecting the most appropriate GP kernel is examined by analysing different aspects, i.e. by considering different criteria as described in Sect.

2.4. The objective is to select the resulting GP model which can achieve a trade-off between the different criteria; -Step 4: since the practitioner is preferably interested in the dependence structure between the scenarios, the correlation matrix derived from the covariance matrix is analysed (see further explanation in Sect. 2.2): this summarizes the inter-dependencies between the levels (scenarios). This result can be confronted to some a priori physicallybased interpretation of the scenario-like variable's influence that the practitioner may have before analysing the computer experiments.

Gaussian Process for mixed continuous and categorical inputs

Let us consider the set of d continuous input variables x=(x1,…,xd), and the set of J categorical inputs u=(u1,…uJ) with L1,…, LJ levels that represent the scenario-like inputs. The output y is then computed using the numerical environmental simulator f(.) as 𝑦 = f(𝐱, 𝐮) = f(𝐰).

In the context of Gaussian Process (GP) modelling (also named as kriging, [START_REF] Williams | Gaussian processes for machine learning[END_REF], the function f(.) is assumed to be a realization of a GP (Y(w)) with a constant mean m and a covariance function k(.,.), named kernel, that can be written as follows:

∀𝐰, 𝐰 ′ , k(𝐰, 𝐰 ′ ) = cov(𝑌(𝐰), 𝑌(𝐰 ′ )) (1)
Let denote (𝐰 1 , … , 𝐰 𝑛 ) the training samples and 𝒚 = (y 1 = f(𝒘 1 ), … , y 𝑛 = f(𝒘 𝑛 )) the corresponding results. The prediction at a new observation 𝐰 * is given by the kriging mean 𝑌 ̂(𝐰 * ) as follows:

𝑌 ̂(𝐰 * ) = E(𝑌(𝐰 * )|𝑌(𝐰 1 ) = 𝑦 1 , … , 𝑌(𝐰 𝑛 ) = 𝑦 𝑛 ) = 𝑚 + 𝒄 𝐰 * 𝑇 . 𝐂 -1 . (𝒚 -𝑚𝐈) (2)
where C is the covariance matrix between the points 𝑌(𝐰 1 ),…,𝑌(𝐰 𝑛 ) whose element is 𝐶[i, j] = 𝑘(𝐰 𝑖 , 𝐰 𝑗 ); 𝒄 𝐰 * is the vector composed of the covariance between Y(𝐰 * ) and the points 𝑌(𝐰 1 ),…, 𝑌(𝐰 𝑛 ), and I is the vector of ones of length n.

The prediction at 𝐰 * can be associated to an error estimate provided by the kriging variance 𝜎 ̂2 given by: 𝜎 ̂2(𝐰 * ) = Var(𝑌(𝐰 * )|𝑌(𝐰 1 ) = 𝑦 1 , … , 𝑌(𝐰 𝑛 ) = 𝑦 𝑛 ) = C(𝐰 . Hence, the covariance function 𝑘 cont can be described by kernel models that are commonly-used in the computer experiment community. In the present study, we restrict the analysis to the stationary two differentiable Matérn 5/2 model [START_REF] Santner | The design and analysis of computer experiments[END_REF]. The categorical covariance functions 𝑘 cat 𝑖 (𝑖 = 1, … , J) can be described in different manners depending on the assumption related to the scenario-like (categorical) input, as described in Sect. 2.3.

In practices, 𝑘 cat 𝑖 can be interpreted, under the homoscedastic assumption, as the kernel (up to a multiplicative constant) of the 1D section 𝑢 𝑖 → 𝑌(𝒙, 𝑢 𝑖 , 𝒖 -𝑖 ) where x, and 𝒖 -𝑖 = (𝑢 1 , … , 𝑢 𝑖-1 , 𝑢 𝑖+1 , … 𝑢 𝐽 ) are fixed. In practices, we preferably use the scaled form of the covariance, i.e. the correlation to ease the interpretation on the dependencies between the levels of the categorical variable (i.e. of the scenarios). With the same notations, the correlation kernel (derived from 𝑘 cat 𝑖 ) is interpreted as the correlation of the 1D section 𝑢 𝑖 → 𝑌(𝒙, 𝑢 𝑖 , 𝒖 -𝑖 ), whatever 𝒙 and 𝒖 -𝑖 . Thus the inspection of 𝑘 cat 𝑖 reveals the correlation of the simulator output explained by the i th categorical input, the others being fixed. Note that such correlation does not depend on the other inputs, which is a result of constructing k(.,.) by tensor product.

Covariance kernel models for categorical inputs

The different options for defining a categorical covariance function are summarized in Table 1.

For sake of clarity, we restrict the presentation to the case of a single input with L levels, which is hereafter denoted by u, because the kernel is constructed by tensor product of 1D categorical inputs (see Sect. 2.2). When the practitioner assumes that no preference can be given to the L levels (i.e. all considered scenarios are considered as having the same influence), 𝑘 cat can be described by an exchangeable covariance [START_REF] Qian | Gaussian process models for computer experiments with qualitative and quantitative factors[END_REF] -also named compound symmetry (denoted CS)function as follows:

𝑘 cat CS (𝑢, 𝑢 ′ ) = { 𝜎 2 if 𝑢 = 𝑢′ 𝜌. 𝜎 2 if 𝑢 ≠ 𝑢′ (4)
where 𝜌 is a unique correlation coefficient satisfying

𝜌 ∈ ]- 1 𝐿-1 , 1[.
When the practitioner assumes that the variable of interest will act differently depending on the considered scenario, but without excluding some dependencies between these different responses, 𝑘 cat can then be described by the most general (and complex) dependence structure where each pairwise coefficient can take a different value depending on the considered levels 𝑢, 𝑢 ′ . The covariance function reads as follows:

𝑘 cat Gen (𝑢, 𝑢′) = { 𝑐 𝑢,𝑢 ′ if 𝑢 ≠ 𝑢 ′ 𝑣 𝑢 if 𝑢 = 𝑢 ′ (5)
The latter structure can be simplified by adding a priori information on the dependence between the levels; for instance by relying on expert-based information. A possible option is to assume that some scenarios perform similarly and that they can be grouped. Assume that the L levels of u are partitioned in G groups, and denote g(u) the group number corresponding to a given level u. Then, the covariance function can be written as [START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF]:

𝑘 cat E (𝑢, 𝑢 ′ ) = 𝑘 cat Gen (𝑔(𝑢), 𝑔(𝑢 ′ )) = { 𝑐 𝑔(𝑢),𝑔(𝑢 ′ ) if g(𝑢) ≠ 𝑔(𝑢 ′ ) 𝑣 𝑔(𝑢) if 𝑔(𝑢) = 𝑔(𝑢 ′ ) (6)
where for all i, j ∈ {1, … , G}, the terms between-group correlations. The structure can be simplified by assuming that the correlation value for each pair of groups is unique by means of a compound symmetry matrix [START_REF] Pinheiro | Mixed-effects models in S and S-PLUS[END_REF].

Instead of deriving the groups based on expert information, a possible option is to derive the groups from the data using a low rank approximation [START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF]) so that the covariance matrix can be defined as:

𝑘 cat LR (𝑢, 𝑢 ′ ) = 𝐅. 𝐅 T (7)
where the matrix F is of size Lq where q is low, with typical values of 2 or 3.

When the practitioner assumes that the levels can be ordered, this means that the categorical variable can be described by an ordinal continuous variable that is not directly observed (i.e. it is said to be "latent") and the levels are seen as discretized values of this ordinal variable [START_REF] Qian | Gaussian process models for computer experiments with qualitative and quantitative factors[END_REF]. The corresponding covariance function can be defined by taking advantage of the tools available for the continuous variables as follows:

𝑘 cat O (𝑢, 𝑢 ′ ) = 𝑘 ̃cont (𝐹(𝑢), 𝐹(u ′ )) (8)
where F( .) is a one-dimensional non-decreasing function (also called warping) so that 𝐹: {1, … , 𝐿} → ℝ, and 𝑘 ̃cont is a one-dimensional continuous kernel.

Model selection

As aforementioned, different kernel modelling choices can be made to represent the categorical variable associated to the scenario-like variable. This raises the question of selecting the most appropriate kernel covariance model. Depending on the modelling objective, different approaches exist for selecting an optimal model with respect to a specific criterion [START_REF] Burnham | Model Selection and Inference A Practical Information Theoretic Approach[END_REF]; for instance, a model that satisfactorily represents (i.e. explains) the relationships between inputs and outputs might not necessarily perform as well as for prediction. Therefore, we propose to examine different viewpoints on the problem of kernel selection by examining different criteria. This multi-criterion approach shares similarities with the data science framework proposed by [START_REF] Yu | Veridical data science[END_REF], who advocate analysing three core principles, predictability, computability, and stability. To select the most appropriate kernel model (given the training dataset), we investigate whether the considered GP model is capable of:

Predictability. It is related to whether the GP model is capable of predicting "yet-unseen" input configurations, i.e. samples that have not been used for training. This can be examined by using cross-validation approaches (e.g. [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]. Two indicators are estimated.

The first indicator, denoted Q², measures the deviation from the true output value. Given a test set T, Q² is defined as follows:

𝑄 2 = 1 - ∑ (𝑦 𝑖 -𝑦 ̂𝑖 ) 2 𝑖∈𝑻 ∑ (𝑦 𝑖 -𝑦 ̅ ) 2 𝑖∈𝑻 (9)
where 𝑦 ̂𝑖 is the i th GP-based prediction of the model output 𝑦 𝑖 , and 𝑦 ̅ = The second criterion is the coverage of the prediction intervals (denoted PI) at the given confidence level  calculated on the test set T and holds as follows:

𝐶𝐴 = 1 |𝑻| ∑ 𝟏 (𝑦 𝑖 𝑃𝐼 𝛼 (𝒘 𝑖 )) 𝑖∈𝑻 ( 10
)
where 1A is the indicator function. CA evaluates whether the model output 𝑦 𝑖 is within the bounds of the prediction interval computed at the input variables 𝒘 𝑖 using the GP mean and variance (Eqs. 2&3). The GP-derived PI is "optimal" when CA reaches the theoretical value of .

Explainability and simplicity. The former concept relates to whether the considered GP model is capable of representing the data, for instance by analysing the likelihood l. However, adding more model parameters results in increasing the explainability. To counterbalance this tendency (related to overfitting), a penalty term is generally introduced (see e.g., [START_REF] Höge | A primer for model selection: The decisive role of model complexity[END_REF] to select simpler model. Here, simplicity is understood with respect to the number of GP model parameters. By assuming that the true GP model exists and that is among the set of candidate GP models, we propose to rely on the Bayesian Information Crierion BIC [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF] described as follows:

𝐵𝐼𝐶 = 2 log(𝑙) + 𝑘. log(𝑛) (11)
where k is the number of parameters, n is the number of observations.

Stability. We explore to which extend the kernel correlation matrix (derived from the covariance matrix) is stable to the perturbations in the training dataset. By following a cross-validation procedure, we evaluate the following error term:

𝑒𝑟𝑟 = 1 𝑛 𝑝 ∑ ( 1 𝑛 cv ∑ (𝑐î ,0 -𝑐î ,𝑗 )² 𝑛 cv 𝑗=1 ) 𝑛 𝑝 𝑖=1 (12)
where 𝑛 𝑝 is the number of terms in the correlation matrix (by restricting to the non-diagonal elements in the upper triangular part of the matrix), 𝑛 cv is the number of folds of the crossvalidation procedure, 𝑐î ,0 is the i th coefficient of the correlation matrix (read in the column order for instance) derived from the GP model fitted using the whole training dataset, 𝑐î ,𝑗 is the i th corresponding coefficient derived from the GP model fitted using the training dataset at the j th iteration of the cross validation procedure (i.e. for a 5-fold cross validation, this training dataset corresponds to the whole training dataset, from which 20% of the observations have been randomly removed).

Application

In this section, we apply the GP-based procedure described in Sect. 2 to a synthetic test function (Sect. 3.1) and to two real cases, in the domain of cyclone-induced wave modelling (Sect. 3.2) and reservoir engineering (Sect. 3.3). For all GP models, we consider a Matérn 5/2 covariance matrix for the continuous variables. A GP model with constant trend is fitted using the R package kergp [START_REF] Deville | kergp: Gaussian process laboratory[END_REF] by applying a pre-scaling and centering of the continuous input variables and of the variable of interest. The covariance parameters are estimated via a maximum likelihood approach using the derivative-free constrained optimiser by linear approximations named COBYLA developed by [START_REF] Powell | A direct search optimization method that models the objective and constraint functions by linear interpolation[END_REF] with 250 randomly selected initial starts. The predictability and the stability are assessment via a 5-fold cross-validation procedure repeated 25 times.

Synthetic case

We first consider a synthetic test function using a modified version of the two-dimensional

Branin function1 where one continuous variable u is assumed to reach only discrete values as follows:

𝑦 = { ℎ(𝑥, -20), 𝑖𝑓 𝑢 = 1 ℎ(𝑥, -10),
𝑖𝑓 𝑢 = 2 ℎ(𝑥, -7.5), 𝑖𝑓 𝑢 = 3 ℎ(𝑥, -5.0), 𝑖𝑓 𝑢 = 4 -5. ℎ(𝑥, 20), 𝑖𝑓 𝑢 = 5 (13) where ℎ(𝑥, 𝑧) = (𝑧 -

5 4 *  2 𝑥 2 + 5  𝑥 -6) 2 + 10. (1 - 1 8
) . cos(𝑥) + 10, with 𝑥 ∈ [-5,10]. By construction, levels 1-4 are highly correlated (see Fig. 1a). Different categorical kernels are defined as follows:

-Compound Symmetry (denoted CS): no preference is given to the scenarios i.e. to the levels of the categorical variable u;

-General (denoted Gen): this is the most generic structure where the pairwise correlation coefficient differs from one level to another; -Expert-based groups: levels are clustered based on expert information. We assume that two experts have given their opinions: the first one (denoted E) indicates a 'realistic' grouping of levels (i.e. consistent with the true function), and the second one (denoted EW) indicates a unrealistic grouping (i.e. (u1,u3), and (u2,u4,u5)); -Low rank approximation: groups of levels are identified through the data-driven low rank approach. A matrix rank of 2 is here tested (denoted LR2); -Ordinal variable: the levels are ordered by following the level index, and a continuous kernel is defined via a spline-based warping (kernel denoted O). As for the expert-based grouping, we assume that another expert does not know the intrinsic ordering and assumes an unrealistic ordering (denoted OW), i.e. u4<u2<u5<u3<u1.

The training dataset is defined through random sampling by considering m points per level with m=4, 5 and 6, i.e. with different training dataset sizes of 20, 25 and 30 respectively. Fig. 2 depicts the GP-derived correlation matrices for each of the afore-described kernel assumptions using the training dataset of intermediate size (m=5). The application of the expertbased and of the ordinal kernel assumptions (Fig. 2c,e) shows some consistent structures among the levels, namely the high-correlated group for u1 to u4 and the particular behaviour of u5. The magnitude of the inter-dependencies between the identified group and u5 slightly differs: assumption E indicates a low-to-moderate correlation (~10%), whereas assumption O indicates a decreasing correlation from 25 to 15% for u4-1. The structure of inter-level dependencies is here richer for the general (Fig. 2a) and the LR2-based GP models (Fig. 2b). Though the former assumption is hardly exploitable, the latter still allows identifying some interesting features, namely the particular behavior of u5, and a possible grouping of u2-4. In this case, the exchangeable assumption (i.e. compound symmetry CS) leads to a low-to-moderate inter-correlation coefficient of ~40% (not shown). The four criteria for kernel model selection are examined in Fig. 3. For the sake of comparability between the different assumptions for m, we preferably plot the difference of BIC with the minimum one (named BIC). Several observations can be made:

-The criterion Q² is commonly used in the computer experiment community to rank different models with respect to their predictive capability. In our case, basing the analysis on this unique criterion is difficult: at low size of the training dataset (m=4), models E, CS and O all minimize 1-Q², and show similar performance (see median values in Fig. 3b): selecting one of them is here hardly achievable. Besides, the model associated to the unrealistic expert-based grouping EW turns to have a satisfactory predictive capability at low m value (though it should be noted that the width of the confidence interval is larger than the others); -For large enough m value -here m=6 -(i.e. with the largest size of the training dataset), the Q² criterion allows to identify model O as the most appropriate one with respect to the predictive capability (log10(1-Q²) is the lowest in Fig. 3b); -The other facet of predictability related to the coverage of the prediction intervals turns to be informative to discard kernel assumptions Gen and LR2 (because CA is here far larger than the level of the prediction interval), but hardly allows differentiating the other assumptions whatever m (Fig. 3c); -Explainability measured by BIC (Fig. 3a) appears here efficient to exclude the unrealistic assumptions OW and EW: they present very large BIC differences whatever m. For instance, [START_REF] Burnham | Model Selection and Inference A Practical Information Theoretic Approach[END_REF] suggested a difference of the considered information criterion (relative to the minimum value) of at least 10 to support the ranking between model candidates with confidence;

-BIC appears to be very informative to discriminate the kernel assumptions, and clearly selects the ordinal assumption as the most appropriate one.

-The stability criterion (Fig. 3d) tends to suffer from the same sensitivity than Q² to the size of the training dataset, but has a higher discriminative power: at low m value (m=4), both kernel models CS and O are selected as very stable (because of low log10(err) values in Fig. 3d). It should however be noted that a high stability is also reached for EW. At a high enough m value, the ordinal assumption is then unambiguously selected as the one leading to the most stable correlation matrix. From this analysis, we can conclude that the ordinal assumption allows to successfully fulfill three of the criteria (explainability, predictability, and stability) especially at high enough m value (m5). For coverage, the ordinal assumption is not ranked first, but CA appears of reasonable order of magnitude (median value of ~85%), i.e. with moderate deviation from the level of the 95% predictive interval. For low size of the training dataset (m=4), selecting unambiguously EW and CS as valid assumptions turns to be difficult if the explainability criterion (BIC criterion) is not taken into account. This can partly be explained by the analysis of the EW-based correlation matrix (Fig. 2d), which reveals a quasi-homogeneous structure with correlation coefficients ranging from 44 to 56%, i.e. of the same order of magnitude than the CS assumption (of ~40%), hence indicating that both GP models should perform similarly.

This suggests that at a low number of training points, the CS assumption remains the most reasonable assumption (when the goal is the joint maximization of predictability, stability and explainability).

Real case application 1: cyclone-induced waves

The first real case is based on [START_REF] Rohmer | Dynamic parameter sensitivity in numerical modelling of cyclone-induced waves: a multi-look approach using advanced meta-modelling techniques[END_REF] and deals with the modelling of waves induced by cyclones at Sainte Suzanne city located in the North East of Reunion Island (Fig. 4b). The aim is to analyse the evolution of the significant wave height Hs (maximum value over time) as a function of the cyclone characteristics. These are modelled by means of five scalar continuous input parameters, namely the maximum wind speed, the radius of maximum winds (i.e. the distance from the cyclone eye at which the maximum wind intensity is reached); the shift around the central pressure; the forward speed defined as the translation speed of the cyclone eye, and the landfall position, that both characterize the minimum distance and the relative position of the track to the studied site. A set of seven historical cyclone tracks are considered: these are randomly shifted (via the continuous input variable modelling the relative position of the track) from their original track so that they cross the centre of Reunion Island (see Fig. 4b): a categorical variable is here defined; each level corresponding to a given track.

A series of 100 computer experiments were performed by randomly sampling the five scalar inputs using a Latin Hypercube Sampling approach combined with a maximin criterion [START_REF] Johnson | Minimax and maximin distance designs[END_REF]. The sampling of the cyclone tracks is done by sampling with replacement. An a priori physical interpretation of the track influence speculates that Hs is strongly related to the angle of approach of the cyclone in the vicinity of the studied site, which increases from 0° (East direction) to 180° (West direction); 90° being the North. This was confirmed by the sensitivity analysis of [START_REF] Rohmer | Dynamic parameter sensitivity in numerical modelling of cyclone-induced waves: a multi-look approach using advanced meta-modelling techniques[END_REF]. Despite the spatial variability of the track (as illustrated in Fig. 4b), the effect of the track scenarios could be summarized by a scalar continuous input, i.e. in relation to the angle of approach. The analysis of the boxplots in Fig. 4a seems to support this hypothesis; in particular, the median value of the maximum Hs appears to increase as the angle of approach increases from 0 to 90° (from Gael to Dumile cyclone).

Yet, the tendency from 90° to 180° (from Dumile to Banzi cyclone) is less clear, especially for Banzi; the difficulty in the interpretation may here be related to the complexity of this cyclone track compared to the quasi linear shape of the others (Fig. 4b). To support the evidence of the angle of approach's influence, a more rigorous analysis is here needed: the validity of this hypothesis is further investigated using a GP model with different categorical kernels as follows:

-Compound Symmetry (denoted CS): no preference is given to the tracks; -General (denoted Gen): this is the most generic structure where the pairwise correlation coefficient differs from one level to another; -Expert-based groups (denoted E): levels are clustered based on expert information. A first assumption relies on the selection of two groups: one composed of the 4 cyclones (Gael, Giovanna, Hollanda, and Dumile) coming from the north-eastern (NE) quadrant and another one composed of 3 tracks (Bejisa, Haliba, and Banzi) coming from the north-western (NW) quadrant (Fig. 4). A second assumption based on three groups is also tested by differentiating the track whose angle is quasi at 90° (North), namely Dumile cyclone (Fig. 4). In addition, two assumptions are made regarding the link between the groups by specifying a general (assumption E2 and E3), or a compound symmetry covariance (assumption E2cs and E3cs); -Low rank approximation: groups of levels are identified through the low rank approach with a assumed rank of 2 or of 3 (assumptions LR2 and LR3); -Ordinal variable: the levels are ordered following the angle of approach, and a continuous kernel is defined via a spline-based warping (kernel denoted O).

Fig. 5 depicts the GP-derived correlation matrices for each of the afore-described kernel assumptions. Fig. 5a reflects the hypothesis of a single group without any preference among the tracks. The derived correlation appears to be high (~80%). The application of alternative kernel assumptions show some consistent structures among the cyclone tracks. Two groups of cyclones appear to be highly correlated (coefficient >75%), namely the ones coming from NE, and the ones coming from NW as shown by Fig. 5b (General formulation), and Fig. 5e,f (E2 andE2cs). The high correlation among these groups is also indicated by the other assumptions:

(1) the low rank approximation, LR2 and LR3 (Fig. 5c,d) -but these assumptions lead to a richer inter-dependency structure; (2) to a lesser extent,the expert-based assumption E3cs (Fig. 5h).

However, we can note that there is an ambiguity for Dumile cyclone (see in particular, Fig. 5g), which is either highly correlated with the NE cyclones (E2, E2cs) or with the NW cyclones (LR3, E3cs). The different assumptions all suggest a moderate correlation (of the order of 40-60%) between groups of cyclones coming from NE and NW. These observations are consistent with the ordinal assumption (Fig. 5h), which indicates a decreasing correlation from Gael (the track with the lowest angle of approach) to Banzi cyclone (the track with the largest angle of approach), and a central role of Dumile, with a decreasing correlation either with the NE or with the NW cyclones. The four criteria of kernel model selection are examined in Fig. 6. Regarding explainability and simplicity, the ordinal assumption O appears to be the most appropriate: the derived GP model presents the minimum BIC value and the differences with the alternative models is large; here larger than 20. Regarding predictability, the ordinal assumption also leads to the best model regarding this criterion. Yet, it should be noted that the predictability of the expert-based model candidates (E2, E2cs, E3, E3cs) remains of moderate-to-high degree (with median value of Q² around 90%). Regarding CA, the use of E3cs allows to reach the level of the 95% prediction interval, but alternative assumptions (E2, E2cs, E3cs, CS and O) lead to a satisfactory coverage of the prediction intervals as well. Finally, Fig. 6d indicates the poor stability of Gen and LR3 kernel, i.e. the high sensitivity of the correlation coefficients' estimates, which may be related to the high number of coefficients to estimate (respectively 21 and 11). This also shows the satisfactory stability of two alternative assumptions, i.e. the expert-based ones and the ordinal one. On this basis, it can be concluded that the ordinal assumption allows to reach a satisfactory trade-off between the four criteria. This appears to be consistent with the afore-described physical intuition; the "latent" continuous variable being here related to the angle of approach. cyclones Gael-Giovanna (and to a lesser extent, for Hollanda as well), and a quasi-linear increasing influence to Banzi. This warping is the basis for computing the correlation matrix (Fig. 5i). To ease the interpretation, let us focus on a single row of Fig. 5i, i.e. the pairwise correlation between Gael track with the others. Fig. 7b provides a clear indication of a highly correlated group of cyclones coming from NE (within the red-coloured envelope), with correlation coefficient exceeding 80%, with decreasing correlation with those coming from NW. The analysis of the correlations derived from the cross-validation iterations (for each repetitions) -grey lines in Fig. 7b -confirms this result: more than 75% of the cross-validationderived results show high correlation (>75%) of Dumile with the NE cyclones, hence in agreement with E2 assumption (Fig. 5e).

Compared to O, the analysis of E2 performance criteria shows that this assumption can also be considered a reasonable one with very satisfactory stability of the correlation coefficients (Fig. 6d), though the criterion values appear to be higher. The stability criterion appears to be lower than for O, which may be related to the lower number of correlation coefficients to be estimated (of 2 for E2 assumption, and of 7 for O assumption).

Real case application 2: CO2 geological storage

The second real case application corresponds to the modelling of the long term fate of stored CO2 in a deep aquifer on a potential project in the Paris basin (France) as described by [START_REF] Manceau | Post-injection trapping of mobile CO 2 in deep aquifers: Assessing the importance of model and parameter uncertainties[END_REF]. The injection of 30 Mt of CO2 during 30 years in the lower Triassic sandstone formation at approximately 1,000m depth was numerically simulated and the evolution of the quantity of mobile CO2 for a time period of 150 years after the injection stops was investigated as a function of:

-two continuous input variables, namely the porosity and the intrinsic permeability of the aquifer rock formation;

-four categorical input variables related to: the scenarios of permeability anisotropy ("minor", "medium", and "large"), the scenarios of regional hydraulic gradient ("absence" and "activated"), the scenarios related to the capillary effect ("absence" and "activated"), and the physical laws used to model the relative permeability as a function of CO2 saturation (ten choices) as depicted in Fig. 8. Due to the importance of the latter parameter (as shown by [START_REF] Manceau | Post-injection trapping of mobile CO 2 in deep aquifers: Assessing the importance of model and parameter uncertainties[END_REF], the following analysis focuses on this variable.

A series of 100 computer experiments were performed by randomly sampling both scalar inputs using a Latin Hypercube Sampling approach combined with a maximin criterion. The sampling of the categorical variables is done by sampling with replacement. Contrarily to the case described in Sect. 3.2, the physical intuition on an a priori influence of scenario-like variable related to the relative permeability laws is harder to give: this is related to the richness of the information associated to the process of residual trapping that is related to different aspects: (1) the capacity of the porous medium to allow the flow of the gaseous phase in the presence of another phase: this is represented as a function of the saturation (see examples in Fig. 8c). This capacity is associated to a potential hysteretic effect resulting in a non-unique dependence over time;

(2) the capacity of the alternative phase (water) that is represented by another function of the saturation (Fig. 8d); (3) the considered phase (gaseous or aqueous) progressively becomes isolated when its saturation decreases in a porous medium.

This leads to a saturation that cannot be reduced: these specific saturations are called residual situation for the gaseous phase (Fig. 8b) and irreducible saturation for the aqueous phase.

To help in the formulation of a physically-based assumption about the inter-dependencies, boxplots in Fig. 8a allow to identify some specific law behaviours especially for law 1; some tendency can also be noticed when ordering the laws in a specific order with respect to the median values of the variable of interest. To get a clearer picture, we test the validity of these observations via the proposed GP-based approach by defining different categorical kernels as follows:

-Compound Symmetry (denoted CS): no preference is given to the physical laws;

-General (denoted Gen): the most generic dependence structure;

-Expert-based groups: the grouping should account for the three facets of the CO2 flow in porous media, i.e. by integrating the three pieces of information depicted in Fig. 8bd: dissimilarities in both the imbibition, the drainage curve's shape and in the residual trapping model. On this basis, the following grouping is proposed: (law 2; law 10); (law 6-9), (law 1); (law 3); (law 4); (law 5). In addition, an assumption is made regarding the link between the groups by specifying a general or a compound symmetry covariance

(assumption respectively denoted E and Ecs);

-Low rank approximation: groups of levels that are identified through the low rank approach. A rank of 2 and of 6 (i.e. of the same number of the expert-based groups) are tested (kernel respectively denoted LR2 and LR6);

-Ordinal variable: the levels are ordered with respect to the value of the maximum CO2 residual saturation (Fig. 8b), and a continuous kernel is defined via a spline-based warping (assumption denoted O).

The other categorical variables are respectively assigned a Compound Symmetric kernel for the regional hydraulic gradient and for the capillary effect. Owing to the "natural" ordering of the levels of the categorical variable for the permeability anisotropy, an ordinal kernel with spline-based warping is defined. A logarithm (base 10) transformation of the variable of interest is applied due to high asymmetry in the variable histogram.

Fig. 9 depicts the GP-derived correlation matrices for each of the afore-described kernel assumptions. Some consistent structures can be noticed regarding law 1, which appears to be anti-correlated with the others as shown in Fig. 9a (general assumption), and in Fig. 9b,c (low rank approximation LR2 and LR6). The specificity of law 1 is also outlined by the expert-based grouping in Fig. 9d, which indicates here a moderate positive correlation of 35-40%. This is in agreement with the ordinal assumption (Fig. 9f) which indicates that the pairwise correlation coefficients of law 1 with the others (see last row of Fig. 9f) rapidly decreases. Though disagreeing on the correlation magnitude, LR2, Gen, and E models all suggest a moderate correlation among all laws except for law 1. The assumption LR6 also suggests the particular behaviour of law 5, which goes in the same direction as the expert-based clustering of considering it as a belonging to a single group. The assumption Ecs (i.e. using a CS betweengroup covariance assumption) leads to a less complex correlation structure and clearly highlights the grouping of law 6-9 (as suggested by the experts, see also the purple-coloured laws in Fig. 8), which is in agreement with the group of highly correlated laws as outlined by Fig. 9f (though the size of the group is larger and includes law 2 and 4 as well).

To summarize, the inspection of the correlation matrices is here more difficult than for the cyclonic application case (Sect. 3.2), where all assumptions more or less agree regarding the information supplied by the modellers. Still, this inspection highlights the specificities of law 1 and law 5 (LR6 assumption), that both strongly differ from the other laws: this is suggested by the irreducible water saturation (red curve in Fig. 8c,d): law 1 is even more different with an irreducible water saturation associated to a large maximum gas residual saturation, which seems to explain why the model behaves so specifically when this law is accounted for. The four criteria for kernel model selection are examined in Fig. 10. For the considered case, we show that the expert-based (Ecs, i.e. with the simplified correlation structure between the groups) and the CS assumption both lead to GP models that satisfactorily fulfill the four criteria, with a slightly higher performance for the simpler structure of CS. This means that the expertbased grouping of laws, in particular of laws 6-9 (in purple in Fig. 8), is informative (in the sense that it leads to a competitive GP model), but only makes a slight difference with the simpler structure especially regarding explainability (with BIC difference ~10) and regarding stability (due to the lowest number of CS correlation coefficients, namely of 1). From this analysis, it can be concluded that, given the 100 simulation results, there is only a mild evidence supporting the assumption of the structure associated to the permeability laws that was intuited from the analysis of the boxplots in Fig. 8a and of the curve similarities (Fig. 8b-d). Similarly as for the synthetic case, additional simulation results should here be performed in order to discriminate unambiguously the most appropriate kernel assumption. Yet, the procedure does not ensure that a unique model is selected: multiple assumptions may eventually turn to be valid (with respect to the four selection criteria) or a trade-off may be difficult to find. In the reservoir case, the application (Sect. 3.3) only moderately supports the evidence of some dependence structure; the compound symmetric and the expert-based kernel model perform similarly with respect to the four selection criteria. Though this result is informative per se, in particular in situations where the practitioner is preferably interested in explaining the numerical results, additional investigations are necessary to confirm this conclusion. This is shown in the synthetic case (Sect. 3.1), where the ordinal assumption was also successfully identified provided that a minimum number of training samples are available.

On the one hand, if additional model runs are computationally affordable, a possible option is to rely on an adaptive sampling strategy. This question deserves however further investigation in the presence of a mixture of continuous and categorical variables, and could be based on recent advances in the context of optimization by [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF] and Munoz Zuniga and Sinoquet (2020). On the other hand, if additional model runs are not possible, an option is to aggregate the information provided by the "plausible" GP models (i.e. the ones that satisfactorily fulfil the criteria) while accounting for some weight reflecting their "plausibility"

(with respect to the selection criteria). This option can take advantage of adequate averaging techniques developed, for instance, within the Bayesian framework as proposed by [START_REF] Zhang | Bayesian model averaging for Kriging regression structure selection[END_REF] for uncertainty quantification, and by [START_REF] Ginsbourger | Discrete Mixtures of Kernels for Kriging-based optimization[END_REF] for optimization problems.

Comparison to a popular alternative

It should be acknowledged that GPs are not the first statistical modelling option that comes to mind when addressing the problem of scenario-like variables. A popular approach relies on tree-based methods like regression decision trees, denoted DT [START_REF] Breiman | Classification and Regression Trees[END_REF], and random forest regression, denoted RF [START_REF] Breiman | Random forests[END_REF], which both natively handle categorical predictors without having to first transform them (e.g., by using feature engineering techniques) because they are based on binary recursive partitioning. Examples of real case applications are provided by Jaxa-Rozen and Kwakkel (2018) and Rohmer et al. (2018).

From a practical viewpoint, the advantage of DT is to provide the structure of dependence (as well as the interactions) with a graphical presentation of the results in the form of a tree, which greatly eases the interpretation. For instance, in the cyclone real case application, Fig. 11a gives the tree structure derived from the analysis of the cyclone real case. chosen at that node: this is confirmed by Table 2. Yet, the high predictability comes at the expense of losing some interpretability, i.e. the ability to represent the structure via the easily understandable tree representation (RF being an ensemble of randomized DT models) though some developments are available to extract some meaningful rules from RF (see e.g. [START_REF] Fokkema | Fitting prediction rule ensembles with R package pre[END_REF]. We acknowledge that there is room for improving the comparison exercise; in particular, we have used commonly-used parametrisations of the tested tree-based methods: further works should include more advanced developments like the enhancement of RF for smooth non-linear relations [START_REF] Friedberg | Local linear forests[END_REF], and adequate splitting rules for categorical inputs (see an extended discussion by [START_REF] Wright | Splitting on categorical predictors in random forests[END_REF].

To conclude, GP models can be seen as a good compromise considering the results on the test cases, because:

(1) They clearly achieve the higher predictability with Q² value >90% given a moderate size of the training dataset (typically 100-200), whatever the considered case (Table 2), while achieving features of high interest (explainability, simplicity, stability) provided that the kernel function is appropriately selecting;

(2) The correlation matrices (Figs. 2,5,9) derived from the kernel function provide a concise and graphical way to get insights into the inter-dependencies among the scenarios (using the formal interpretation provided in Sect. 2.2).

These results highlight the added values of our GP-based approach and show that it can be considered a valuable tool to complement the toolbox of any geo-scientist in order to both explore and characterize model uncertainties related to scenario-like inputs.

Figure 1

 1 Figure 1 (a). Synthetic test function with continuous input variable x. Each color indicates a different level of the scenario-like (categorical) variable u. Each dot corresponds to a model run. (b) Correlation matrix for u derived from the GP-based analysis using 25 model runs.

  the test set. A coefficient Q² close to 1.0 indicates that the GP model is successful in matching the new observations that have not been used for the training.

Figure 2 .

 2 Figure 2. Correlation matrix for the synthetic test case (with m=5) considering different assumptions regarding the kernel associated to the categorical input variable denoted u. For the expert-based assumptions (c and d), the ordering of matrix coefficients follows the expert-based clustering, i.e. (u1-u4), and (u5).

Figure 3 .

 3 Figure 3. Selection criterion for the synthetic test case considering different numbers of points per level m: (a) BIC difference with respect to the minimum value; (b) Predictability measured by 1-Q² (log10 scale); (c) Coverage measured by 1-CA. The horizontal dashed line corresponds to 5%, i.e. the threshold consistent with the level of the 95% prediction interval; (d) Stability error err (log10 scale). The three latter criteria are derived from the 5-fold cross validation repeated 25 times: the height of the barplot is the median value and the lower and upper bounds are defined using the 25 th and 75 th percentiles.

Figure 4 .

 4 Figure 4. (a) Boxplot of the maximum significant wave height Hs (m) considering each cyclone track ordered according to the angle of approach from 0° to 180° (from the Eastleftmost part, to the Westrightmost part); (b) Cyclone tracks used for modelling the waves at Saint Suzanne city (Reunion island).

Figure 5 .

 5 Figure 5. Correlation matrix for the cyclonic test case considering different assumptions regarding the kernel associated to the categorical input variable.

Figure 6 .

 6 Figure 6. Selection criterion for the application case 1: (a) BIC; (b) Predictability measured by 1-Q²; (c) Coverage measured by 1-CA. The horizontal dashed line corresponds to 5%, i.e. the threshold consistent with the level of the 95% prediction interval; (d) Stability error err. Criteria (b-d) are derived from a 5-fold cross validation repeated 25 times: the height of the barplot is the median value and the lower and upper bounds are defined using the 25 th and 75 th percentiles.

Figure

  Figure 7a summarizes this dependence via the spline-based warping function F(.) used to setup the ordinal covariance kernel model 𝑘 cat O (see Eq. 8): this shows a strong link between NE

Figure 7 .

 7 Figure 7. (a) Spline-based (un-normalized) warping function F (see Eq. 8) used in the ordinal assumption for the cyclone case. (b) Pairwise correlation between Gael track with the other ones. The red-coloured envelope indicates cyclones coming from North-East (denoted NE). The grey-coloured curves are the correlations derived at each iteration of the 5-fold crossvalidation procedure (repeated 25 times). The red dashed line indicates the median value.

Figure 8 .

 8 Figure 8. (a) Boxplot of the amount of mobile CO2 considering each physical law index (1-10). Colors indicate the expert-based grouping of the laws. (b) Maximum residual saturation of CO2 for each law. Relative permeability law used in the reservoir test case: (c) water; (d) CO2.

Figure 9 .

 9 Figure 9. Correlation matrix for the reservoir test case considering different assumptions regarding the kernel associated to the categorical input variable. For the expert-based assumptions (d and e), the ordering of matrix coefficients follows the expert-based clustering, i.e. (law 2; law 10); (law 6-9), (law 1); (law 3); (law 4); (law 5). For the ordinal assumption (f), ordering of matrix coefficients follows the ordering of the maximum residual saturation of CO2.

Figure 10 .

 10 Figure 10. Selection criterion for the application case 2: (a) BIC; (b) Predictability measured by 1-Q²; (c) Coverage measured by 1-CA. The horizontal dashed line corresponds to 5%, i.e. the threshold consistent with the level of the 95% prediction interval; (d) Stability error err.

Figure 11 .

 11 Figure 11. (a) Tree structures derived from the analysis of the cyclone real case; (b-d) Tree structures considering three iterations of the 5-fold cross validation procedure. The decision rule is provided on each respective branch. The leaf (bottom node) are coloured according to the mean of the variable of interest (scaled value of HS); the number in percentage provides the number of samples in each leaf. The track name associated to the 'Track' node corresponds to the two first letters of the cyclone names provided in Fig. 4. The green rectangle outlines a particular grouping of tracks.

Table 1 .

 1 Different options for representing the scenario-like input variable using a kernel

	(covariance) model
	Assumption			Interpretation	of	Type	of	Symbol Equation
					categorical variable		covariance
							matrix
	No preference among the	Each level acts similarly Compound	CS	4
	scenarios					Symmetry
							(also named
							Exchangeable)
	The variable of interest	Each level has its own	General	Gen	5
	will	act	differently	variance coefficient and a
	depending	on	the	between-level structure
	considered scenario, but	exists
	without excluding inter-
	level dependencies
	The variable of interest	Some scenarios can be	Group	E	6
	may	act	similarly	grouped based on expert
	depending	on	some	information
	subsets of scenarios
	A grouping is a realistic	The grouping can be	Low	rank	LR	7
	option, but cannot be	learnt from the data		approximation
	unambiguously defined
	The scenarios can be	The levels are discretized	Ordinal	O	8
	ordered				values of a latent ordinal
					variable

http://www.sfu.ca/~ssurjano/branin.html

Acknowledgements

This research was conducted within the frame of the Chair in Applied Mathematics OQUAIDO, gathering partners in technological research (BRGM, CEA, IFPEN, IRSN, Safran, Storengy) and academia (CNRS, Ecole Centrale de Lyon, Mines Saint-Etienne, University of Grenoble, University of Nice, University of Toulouse) around advanced methods for Computer Experiments.

Criteria (b-d) are derived from a 5-fold cross validation repeated 25 times: the height of the barplot is the median value and the lower and upper bounds are defined using the 25 th and 75 th percentiles.

Summary, Discussion and Further works

Model (related to the structure/form of the model or to the unambiguous choice of "appropriate" physical laws) are generally integrated in environmental models via scenariolike variables, i.e. categorical variables. In the present study, we have explored the applicability of GP meta-modelling in order to:

-Inform on the structure of dependence for scenario-like inputs (research question Q1) using different formulations of the categorical covariance function (exchangeable, ordinal, group, etc.); -Derive a predictive statistical model using a limited number of computer experiments, i.e. 100-200 (research question Q2) that satisfactory fulfil predictive capability, explainability and stability of parameters' estimates.

Regarding Q1, we discuss the added value of the approach in Sect. 4.1, as well as its limitations.

Regarding Q2, we compare in Sect. 4.2 the GP-based procedure with popular modelling alternatives, namely tree-based methods, with respect to predictability as well as from a practical viewpoint. Different lines for further research works are also outlined.

Selecting the dependence structure

We have proposed to rely on a multiple-criterion selection approach to provide indications on which kernel model is the most appropriate to represent the dependence structure of the scenario-like inputs. By construction, the proposed framework should provide a level of sufficient flexibility to critically analyse the different physically-based assumptions: the advantage of testing and confronting each assumption to one another is to keep the procedure transparent to the practitioner. This flexibility is well shown in the cyclone application case (Sect. 3.2), where the proposed procedure allows us to confront an a priori physically-based assumption (i.e. the cyclone track effect can be summarized by a scalar ordinal variable) to alternative views on the dependence structure, and to support the evidence of the a priori assumption.

However, we note, from the analysis of Fig. 11, that several differences between the structure constructed using the DT trained using the whole dataset (Fig. 11a) and the ones at each iteration of the cross-validation procedure (i.e. using three DT models setup with "perturbed" training database, Fig. 11b-d); in particular for the leftmost part related to the track variable (outlined by a green rectangle), the grouping of tracks are similar for Fig. 11a and Fig. 11c, but differs for Fig. 11b and is even absent for Fig. 11d. This high sensitivity of the derived structure (to the changes in the training dataset) has already been identified in the literature [START_REF] Breiman | Random forests[END_REF]: in addition to bringing some confusion regarding the scenario dependencies, the drawback is also a poorer predictability: this is shown by the low Q² values for each test case in Table 2.

Table 2. Predictability measured by the Q² indicator for different statistical models. For the synthetic case with different points per level m, Q² is derived from the leave-one-out crossvalidation procedure. For the application cases 1 and 2, Q² is derived from the 5fold crossvalidation procedure (repeated 25 times): the median values are given together with and the 25 th and 75 th indicated in brackets.

Model Regression


Lucian Florin

Constantin
email: fl.constantin@univ-lille1.fr

Characterization of a Static Magnetic Field with Two-Photon Rotational Spectroscopy of Cold Trapped HD +

Keywords: two-photon spectroscopy, trapped ions, hydrogen molecular ions, Zeeman shift, ab-initio calculation, vector magnetometry

Two-photon rotational spectroscopy of cold trapped HD + ions may be exploited for characterization of the magnitude and orientation of a static magnetic field in the ion trap. The experimental setup and the approach for calibration of the magnetic field vector are described. A sensitivity at the 10 -11 T level may be reached with Zeeman spectroscopy of a hyperfine component of the (v,L)=(0,0)→(0,2) transition. The orientation of a magnetic field with a magnitude at the 1 µT level may be characterized with an uncertainty better than 30 mrad.

I. INTRODUCTION

Precise sensing of the magnetic fields has applications in geophysics [START_REF] Stuart | Earth's field magnetometry[END_REF], medical imaging [START_REF] Boto | On the potential of a new generation of magnetometers for MEG: A beamformer simulation study[END_REF][START_REF] Alem | Fetal magnetocardiography measurements with an array of microfabricated optically pumped magnetometers[END_REF], and precision physics measurements [START_REF] Budker | Optical Magnetometry[END_REF]. The atomic magnetometers exploiting alkali atoms vapours are extremely sensitive [START_REF] Kominis | A subfemtotesla multichannel atomic magnetometer[END_REF], but their response, based on Larmor precession, is intrinsically scalar. This work evaluates the potential for SI-traceable full characterization of a static magnetic field from precision measurements of transitions of trapped molecular ions compared with the predictions of the state-of-the-art ab-initio energy level calculations. Precisely, the energy levels of the hydrogen molecular ions (HMI), the simplest molecular quantum systems, were predicted with quantum electrodynamics (QED) ab-initio calculations reaching a fractional accuracy at the 10 -12 level [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF]. The shifts of the energy levels due to various external fields were also predicted by ab-initio calculations [START_REF] Bakalov | Magnetic field effects in the transitions of the HD + molecular ion and precision spectroscopy[END_REF][START_REF] Schiller | Static and dynamic polarizability and the Stark and blackbody-radiation frequency shifts of the molecular hydrogen ions H2 + , HD + , and D2 +[END_REF]. Particularly, the effects of a static magnetic field on the HD + energy levels were evaluated in the nonrelativistic approximation at the order α 2 . The ab-initio calculations allowed identifying and exploiting for precision measurements transitions which are insensitive to the magnetic field and to its orientation with respect to the polarization of the radiation field. The rovibrational and rotational spectroscopy of trapped and sympathetically cooled HD + ions was performed by resonantenhanced multiphoton dissociation (REMPD). The Dopplerlimited infrared spectroscopy allowed a fractional uncertainty at the 10 -9 level [START_REF] Biesheuvel | Probing QED and fundamental constants through laser spectroscopy of vibrational transitions in HD +[END_REF]. A significant increase in the accuracy and resolution can be provided by the Doppler-free spectroscopy. The rotational spectroscopy of HD + in the Lamb-Dicke regime allowed kHz-level resolution and a fractional accuracy at the 10 -10 level [START_REF] Alighanbari | Rotational spectroscopy of cold and trapped molecular ions in the Lamb-Dicke regime[END_REF]. A double-resonance two-photon spectroscopy scheme of HD + may provide fractional accuracy and resolution at the 10 -12 level or better [START_REF] Constantin | Double-resonance two-photon spectroscopy of hydrogen molecular ions for improved determination of fundamental constants[END_REF].

In precision measurements, a static magnetic field, generated with a set of coils, allows to define the quantization axis of the molecular ions and to split their energy levels. This contribution proposes to exploit the scalar response obtained from the HD + ions by two-photon Zeeman spectroscopy to characterize the magnitude and orientation of the magnetic field.

II. THEORY AND SPECTROSCOPY OF HD +

The rovibrational energy levels of HD + in its ground electronic state are calculated as the nonrelativistic Schrödinger energy plus a series expansion of QED correction terms [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF]. At high resolution, the energy levels reveal a hyperfine structure described with the coupling of the proton p

         + = + = + =
. The hyperfine eigenstates |vLFSJ>, labelled with the quantum numbers for the vibration v and for the angular momentum coupling scheme, are calculated using a Breit-type spin Hamiltonian and accurate variational Coulomb wavefunctions [START_REF] Bakalov | High-precision calculation of the hyperfine structure of the HD + ion[END_REF]. A small magnetic field is used to define the quantization axis and to split the hyperfine states in magnetic subcomponents, labelled with the quantum number Jz of the projection of J  on the quantization axis. The Zeeman shift of an energy level in an external static magnetic field may be described with a quadratic dependence in function of the magnetic field magnitude [START_REF] Bakalov | Magnetic field effects in the transitions of the HD + molecular ion and precision spectroscopy[END_REF]. The total energy of a magnetic level can be expressed as :

( ) ( ) ( ) ( ) ;B z J v,L,F,S,J, z E v,L,F,S,J hf E v,L rv E z J v,L,F,S,J, E Δ + + = ( 1
)
in function of the rovibrational energy

( ), v,L rv E the hyperfine energy ( ) v,L,F,S,J hf E
, and the Zeeman shift ( )

;B z J v,L,F,S,J, z E Δ , respectively.
Accurate measurements of the HD + transitions can be performed with an experimental setup based on cold trapped ions [START_REF] Biesheuvel | Probing QED and fundamental constants through laser spectroscopy of vibrational transitions in HD +[END_REF][START_REF] Alighanbari | Rotational spectroscopy of cold and trapped molecular ions in the Lamb-Dicke regime[END_REF]. Precisely, ~10 2 HD + ions and ~10 3 Be + ions are trapped together in a radiofrequency trap. The Be + ions are laser cooled with a 313 nm laser. The electrostatic interactions between two species allow sympathetical cooling of the HD + ions. When the secular motion of the HD + ions is excited in the trap, the increase in fluorescence of the laser-cooled Be + ions is used to monitor the number of the trapped HD + ions. A static magnetic field, controlled with three coil pairs, driven with three independent current sources, is applied to the ion trap. The measurements proposed in this contribution are based on the double resonance spectroscopy scheme discussed in [START_REF] Constantin | Double-resonance two-photon spectroscopy of hydrogen molecular ions for improved determination of fundamental constants[END_REF], that addresses the two-photon rotational transition in the vibrational ground state (v,L)=(0,0)→(0,2) and the twophoton rovibrational transition (v,L)=(0,0)→(2,0). The detection is performed by the photodissociation of the (v,L)=(2,0) level with a 175 nm laser. The Doppler-free spectroscopy is performed with two counterpropagating waves tuned around the rovibrational transition at 55.909 THz, and the rotational transition at 1.968 THz, respectively (Fig. 1). The two-photon transition rates between the magnetic levels are derived using the two-photon operator formalism and the second-order time-dependent perturbation theory. The interaction with the blackbody radiation at room temperature recycles continuously the population among the rotational levels of the ground vibrational state. Moreover, the population in an excited energy level decays radiatively. The change of the population in the rovibrational levels is calculated with a set of rate equations in order to derive the lineshapes of the two-photon resonances [START_REF] Constantin | Double-resonance two-photon spectroscopy of hydrogen molecular ions for improved determination of fundamental constants[END_REF]. The full-width half measured linewidth is 12.1 Hz for the rovibrational line and 0.447 Hz for the rotational line, by assuming a twophoton rotational transition rate of 2000 s -1 , a two-photon rovibrational transition rate of 10 s -1 , a dissociation rate of 200 s -1 and a REMPD time of 10 s. The fractional frequency uncertainty is estimated with the Allan variance for the quantum projection noise limit :

( ) τ π 1 τ σ c ion y T N Q × = (2)
expressed in function of the quality factor of the two-photon transition Q = f2ph/ΔfHWHM in terms of the half-linewidth determined previously. The cycle time Tc is associated to a single measurement with Nion ions at the two-photon resonance. Successive measurements are averaged during an interrogation time τ. Assuming the case of single ion spectroscopy experiment with Tc = τ, the frequency uncertainty for (v,L)=(0,0)→(2,0) line is estimated at 1.9 Hz and for (v,L)=(0,0)→(0,2) line at 71 mHz.

The Zeeman shift of a transition |v,L,F,S,J,Jz>=(n,Jz)→(n',J'z) may be calculated approximatively with a quadratic dependence :

( ) ( ) ( ) ( ) 2 2 0 B B B ' J , J , r , q , t , r , q , t B ' J , J , r , q , t , r , q , t f f z z ' n ' n ' n n n n B z z ' n ' n ' n n n n B ' n , n ' n , n η + η = = - ( 3
)
in function of the magnitude of the magnetic field B and a set

of six theoretical parameters { } { } ' n ' n ' n n n n th r , q , t , r , q , t U =
that were calculated ab-initio for the relevant hyperfine energy levels of HD + in [START_REF] Bakalov | Magnetic field effects in the transitions of the HD + molecular ion and precision spectroscopy[END_REF]. The favoured (intense) twophoton electric dipole (E1-E1) transitions obey the selection rules ΔF=ΔS=0, ΔJ=0,±2. The linear and quadratic Zeeman shift coefficients are calculated and displayed in Table I for a selection of magnetic subcomponents of the rotational and rovibrational two-photon transitions. The subcomponents of the rotational transition, in general more sensitive than the subcomponents of the rovibrational transition, are suitable for magnetic field characterization. When the THz source used for two-photon spectroscopy is referenced to a frequency standard, this approach allows absolute calibration of the magnetic field in the ion trap, using Zeeman spectroscopy measurements, ab-initio calculations of the HD + energy levels and of their Zeeman shifts. v,L,F,S,J;Jz v',L',F',S',J';J'z ΔηB(kHz/T) ΔηB2(MHz/T 2 ) 0, 0, 1, 0, 0; 0 2, 0, 1, 0, 0; 0 0 40900 0, 0, 0, 1, 1; ± 1 2, 0, 0, 1, 1; ± 1  155 4700 0, 0, 1, 2, 2;  2 0, 2, 1, 2, 4; 0 ± 6992320 1855850 0, 0, 1, 1, 1;  1 0, 2, 1, 1, 3; ± 1 ± 6123100 397050

III. CHARACTERIZATION OF A STATIC MAGNETIC FIELD

This part will describe a procedure for measuring the magnitude and orientation of a static magnetic field in the Cartesian Laboratory Coordinate Frame 3). A choice between the two values of the magnetic field derived from the second-order equation may be made by comparison with the value obtained with an accurate magnetometer. The uncertainty of the magnetic field may be calculated with the error propagation formula, in function of the contribution from the frequency shift measurement uncertainty Hz, 0.1 σ = rot exp, estimated with the quantum projection noise limit, and the contributions from the theoretical parameters, which are assumed correlated and having the same uncertainty

kHz/T. 5 σ , MHz/T 0 5 σ 2 = = t , th r , q , th
Here, the Zeeman spectroscopy of the subcomponent |v,L,F,S,J,Jz>=|0,0,1,2,2,-2>→|0,2,1,2,4,0>, having the highest Zeeman shift of 717.791 kHz at B = 10 -4 T, is exploited for measuring the magnetic field. The dependence of the uncertainty with the magnitude of the magnetic field in the ion trap is plotted in Figure 2. The sensitivity of the twophoton Zeeman spectroscopy method is given by the experimental uncertainty of 1.4×10 -11 T. Moreover, a magnetic field smaller than 10 -4 T may be calibrated with a precision better than 7×10 -11 T. Fig. 2. Total uncertainty (magenta line), the contribution from the frequency measurement (red line), and the contribution from the theoretical calculation (blue line) of the magnitude of the magnetic field measured by Zeeman spectroscopy of the |v,L,F,S,J,Jz>=|0,0,1,2,2,-2>→|0,2,1,2,4,0> transition of HD + . The frequency uncertainty is estimated with the quantum projection noise-limited instability for a single-ion experiment with integration time 1 s and cycle time 1 s.

The directional measurements of the magnetic field are determined by a set of experimental parameters: three currents ( )

3 2 1 I , I , I
that drive independently the coil pairs in Helmholtz configuration, the current-to-field parameters of the coil pairs ( )

3 2 1 k , k , k
, three angles that define the relative orientation of the coil pairs ( )

z y x , , α α α
, and the offset magnetic field components ( )

03 02 01 B , B , B
, due for example to the Earth's magnetic field. As discussed previously, the magnitude of the magnetic field can be expressed in function of the experimental frequency of a two-photon transition

: { } ( ) i i , th i exp, i , 3 i , 2 i , 1 i k f f I , I , I ; V B η - =  (4)
f f f + =
of the rotational frequency and the hyperfine frequency. A nonlinear least-squares adjustment of N Zeeman subcomponents can be exploited for absolute calibration of the experimental parameters { } k V (at least N=9 measurements are required). In order to estimate the uncertainties, the magnetic field dependence of eq. ( 4) is linearized around the design values of the parameters { }

0 ; k V : { } ( ) { } ( ) { } ( ) { } { } ( ) 0 3 2 1 3 2 1 0 3 2 1 0 ; k k k V V k i , i , i , i k i , i , i , i ; k i , i , i , i k V V V I , I , I ; V B I , I , I ; V B I , I , I ; V B ; k k - ∂ ∂ + =  =    (5)
The 9×9 covariance matrix G of the adjusted parameters may be estimated in function of the N×N covariance matrix X of the input data

( ) i i , th i exp, i f f b η - =
, and the N×9 Jacobian matrix ( )

{ } { }       ∂ ∂ = = 0 ; k k V V k , i k i V B A  : ( ) 1 1 - - = A X A G T ( 6
δ         + +         +         +         = i i , hf i , rot i i i , hf i i , rot i i exp, i f f f f f b (7)
The nondiagonal elements of the { } ij x X = matrix arise from the covariances between the theoretical parameters, expressed as :

( ) ( )         +         + + + + ×         =

Fig. 1 .

 1 Fig. 1. The experimental setup and the reference coordinate frames.

  two-photon Zeeman spectroscopy of HD + (Fig. 1). The three coil pairs define a Coil Coordinate Frame necessary orthogonal. The relative orientations of the CCF axis relative to LCF are defined with the Euler angles : ( ) the Zeeman shift δf of a magnetic subcomponent of the two photon rotational transition |vLFSJJz>=|n,Jz>→|n ' ,J'z> allows to calculate the magnitude of the magnetic field in function of the relevant values for the theoretical parameters { } th U , by inverting eq. (

  Zeeman effect from the theoretical value i , th f calculated at zero-field in proportion with the linear Zeeman shift coefficient i η

)

  The uncertainties of the input data are estimated as the quadratic sum of the contributions from the experimental uncertainty

  calibration is performed by Zeeman spectroscopy of the |v,L,F,S,J,Jz>=|0,0,1,2,2,-2>→|0,2,1,2,4,0> HD + transition. The Zeeman two-photon frequencies are measured for 27 sets ( ) intensities, where I0=1 A, Ioff1=3.82 mA, Ioff2=208

TABLE I .

 I

ZEEMAN SHIFT COEFFICIENTS FOR SELECTED COMPONENTS OF THE TWO-PHOTON ROTATIONAL AND ROVIBRATIONAL TRANSITIONS.

mA, Ioff3=422 mA, and n1,2,3=0,1,2. The uncertainties of the values of the experimental parameters, that can be determined by the nonlinear adjustment, are estimated with eq. (4-8) :

T/A, 10 4 5 Cancellation of the Earth's magnetic field in the ion trap is realized using the values of the current intensities that minimize eq. ( 4) : Inull1 = -6.34 mA, Inull2 = 230 mA, Inull3 = 434 mA. The uncertainty in setting the magnitude of the magnetic field is calculated with the root sum of squares of the contributions from the uncertainties of setting the current intensities (assumed at a fractional value of 10 -3 of each current intensity), and from the uncertainties of the experimental parameters determined previously by Zeeman spectroscopy. Particularly, the uncertainty in setting the null magnetic field in the ion trap is estimated at 3.9×10 -8 T, if the theoretical uncertainties are neglected. In addition, a conservative estimation of the orientation uncertainty for a magnetic field with a given magnitude may be derived by assuming that the total magnetic field uncertainty is perpendicular to the magnetic field direction. If the three pairs of coils are driven by (-6.33 mA, 240 mA, 444 mA) currents, the magnitude of the magnetic field is 1.7×10 -6 T and the orientation uncertainty is 30 mrad.

IV. CONCLUSION

This contribution describes the experimental setup and the calibration technique for the magnitude and the orientation of the magnetic field in an ion trap using precision Zeeman spectroscopy of HD + ions. The Zeeman shift for a selected hyperfine component of a two-photon rotational transition is exploited for the determination of the magnitude of a static magnetic field with an uncertainty at the 10 -11 T level. The approach allows coil calibration that enables to predict the magnitude and orientation of the magnetic field in function of the coil currents. The approach allows cancellation of the Earth's magnetic field at the 39 nT level. The orientation uncertainty of a magnetic field of 1.7 µT is better than 30 mrad.


Hao Xu

Max Langer

Françoise Peyrin

SEGMENTATION OF BONE VESSELS IN 3D MICRO-CT IMAGES USING THE MONOGENIC SIGNAL PHASE AND WATERSHED

Keywords: synchrotron radiation, micro-tomography, watershed, monogenic signal, vessel segmentation

We propose an algorithm based on marker-controlled watershed and the monogenic signal phase asymmetry for the segmentation of bone and micro-vessels in mouse bone. The images are acquired using synchrotron radiation micro-computed tomography (SR-µCT). The marker image is generated with hysteresis thresholding and morphological filters. The control surface is generated using the phase asymmetry of the monogenic signal in order to detect edgelike structures only, as well as improving detection in low contrast areas, such as bone-vessel interfaces. The quality of segmentation is evaluated by comparing to manually segmented images using the Dice coefficient. The proposed method shows substantial improvement compared to a previously proposed method based on hysteresis thresholding, as well as compared to using the gradient image as control surface. The algorithm was applied to images of healthy and metastatic bone, permitting quantification of both bone and vessel structures.

INTRODUCTION

Breast cancer is one of the most common cancer in women, with approximately 508,000 deaths annually worldwide [START_REF] Tulotta | Endogenous Production of IL1B by Breast Cancer Cells Drives Metastasis and Colonization of the Bone Microenvironment[END_REF] [START_REF] Coleman | 6. Box plots of parameters extracted from the imaged mouse bone[END_REF]. Around 70% of patients with breast cancer develop bone metastases, leading to osteolytic and osteoblastic lesions [START_REF] Akhtari | Biology of breast cancer bone metastasis[END_REF]. Bone metastasis is currently incurable, and median overall survival times have been reported in the range of 24-65 months [START_REF] Pulido | Bone metastasis risk factors in breast cancer[END_REF]- [START_REF] Lote | Bone Metastasis Prognosis, Diagnosis and Treatment[END_REF]. The vessel networks in bone is thought to play an important role in cancer and metastasis growth. In breast cancer bone metastasis, there is an abnormal growth of blood vessels in the bone [START_REF] Herve Nyangoga | Three-Dimensional Characterization of the Vascular Bed in Bone Metastasis of the Rat by Microcomputed Tomography (MicroCT)[END_REF].

Previously, three-dimensional (3D) X-ray µCT with a contrast agent has been used to visualize the bone vascularization [START_REF] Zhang | Periosteal Progenitor Cell Fate in Segmental Cortical Bone Graft Transplantations: Implications for Functional Tissue Engineering[END_REF]. This required demineralization of the bone however, precluding simultaneous quantitative analysis of bone and vessels. SR X-ray phase contrast CT has been used for the simultaneous visualization of calcified bone tissue and vasculature in C57BL/6 mice post mortem without any staining or contrast agent [START_REF] Núñez | Simultaneous visualisation of calcified bone microstructure and intracortical vasculature using synchrotron Xray phase contrast-enhanced tomography[END_REF]. However, due to the high resolution necessary for this imaging protocol, the field of view is limited. Further, the phase contrast in vessels was relatively weak, so that vessels had to be segmented manually, making the applicability to large data sets limited.

X-ray µCT has been used to examine alterations in cortical bone micro-porosity, mineralization, and cancellous bone architecture due to estrogen deficiency in an ovariectomized rat model of postmenopausal osteoporosis [START_REF] Sharma | The effects of estrogen deficiency on cortical bone microporosity and mineralization[END_REF]. In this study, the vessels were not directly imaged, but only cortical vascular canals. We previously developed a protocol to simultaneously analyze in 3D the micro-vascularization and bone structure in rat bone using SR-µCT images [START_REF] Langer | Simultaneous 3D imaging of bone and vessel microstructure in a rat model[END_REF]. This protocol was based on hysteresis thresholding and morphological filtering, and permitted the quantification of large data sets.

The aim of this work is to quantify bone and vascularization in mouse bone based on SR-µCT images. Transition from rats to mice makes models of many pathologies available, but presents a number of challenges. In mice bone, as opposed to rat, at the chosen imaging resolution (3.5 µm voxel size), vessels appear to be in contact with the bone surface, precluding the use of hysteresis thresholding for the separation of bone and vessels. In addition, the contrast between bone and vessels contrast agent may be weaker.

To this aim, we present a 3D segmentation method based on the marker-controlled watershed algorithm [START_REF] Meyer | Morphological segmentation[END_REF] to separate bone and vessel compartments. An important step for this algorithm is the generation of markers. Here, we assume three segments: bone, vessels and background. Ideally, the marker image should contain no false positives, while including a marker in each connected component of each compartment. For this purpose, we used hysteresis thresholding and morphological filters. The other essential step is the generation of a control surface image from the image to be segmented, where the "flooding" of the watershed algorithm happens. Classically, the gradient of the original image is used. Due to the relatively low contrast between bone and vessel compared to vessel or bone to background, the gradient image does not yield a satisfying segmentation. Therefore, we propose to use the local phase asymmetry derived from the monogenic signal to generate the control surface, which quantifies whether the image is locally odd (edge-like) or even (line-like). This measure has the interesting property of invariance to signal intensity [13][14].

The algorithm was evaluated using manually segmented volumes extracted from the original data. Segmentation quality was measured using the Dice coefficient. The method was compared to our previous method based on hysteresis thresholding as well as marker-controlled watershed on the gradient. The proposed method shows substantial improvement. The algorithm was used to extract parameters to characterize bone and vessel structure in healthy and metastatic bone.

MATERIALS AND METHODS

Sample preparation and image acquisition

8 week-old female Balb/c nude mice were inoculated with luciferase-expressing human B02 breast cancer cells in the tail artery [START_REF] Bachelier | Combination of anti-angiogenic therapies reduces osteolysis and tumor burden in experimental breast cancer bone metastasis[END_REF]. 22 days after the inoculation of tumor cells, the mice were sacrificed. The mice were perfused with contrast agent in the vascular system for imaging. We consider here two groups of samples: healthy and metastatic, with 5 mice per group.

Imaging was performed using SR-μCT at the European Synchrotron Radiation Facility (ESRF) on the ID19 beamline. A sketch of the experimental setup is shown in Fig. 1. The X-ray energy was set to 26 keV using the first harmonic of undulator radiation. Exposure time was set to 0.15 s per radiograph for approximately 8 minutes acquisition for one sample. 2000 radiographs of the sample were recorded using a 2048×2048 pixel CCD-based detector at evenly spaced angles of view over a 360° rotation. This gave a field of view (FOV) as a cylinder of diameter 7 mm. At the end, 2000×2000×1200 voxels 3D images with an isotropic voxel size of 3.5 μm were reconstructed from the projections using filtered back-projection.

Marker based watershed

The marker-controlled watershed algorithm permits to initialize the watershed from seeds in the image, which are known by some means [START_REF] Meyer | Morphological segmentation[END_REF]. This can reduce the number of iterations of the watershed algorithm, and simplifies the merging the resulting segments. The marker image should ideally contain at least one segment in each connected component of each class, while having a false positive rate close to zero. Here, we propose to use hysteresis thresholding and morphological thinning to generate the marker image.

Markers are generated for three segments in the image: bone, vessels and background. The procedure for the segmentation of the markers are given in Fig. 2. Background markers are given by blue arrows. Hysteresis thresholding is used to make an initial segmentation of the background. In order to ensure a low falsepositive rate, one iteration of morphological thinning is applied. Vessels markers are given by red arrows (Fig. 2). As with the background markers, an initial segmentation is performed using hysteresis thresholding which is followed by one iteration of morphological thinning. Bone markers are given by green arrows (Fig. 2). First, both bone and vessels are segmented together using hysteresis thresholding. Then, the vessel markers are subtracted and morphological opening is used to remove isolated small particles.

Finally, one iteration of morphological thinning is applied to minimize false positives. Background markers, vessel markers, and bone markers are shown in blue, red and green respectively in the final marker image (Fig. 2).

Monogenic signal phase

To steer the watershed algorithm, a control surface is generated from the original image. This can be the image itself, but is often based on the magnitude of the gradient. Here, the magnitude of the gradient does not perform satisfactory, due to the relatively low contrast at interfaces between bone and vessel, compared to vessel or bone to background. To alleviate this problem, we propose to use the local phase asymmetry of the 3D monogenic signal as the control surface for the watershed algorithm. The monogenic signal phase has the interesting property of being invariant to signal intensity: it only depends on the local shape of the signal. Here, we want to detect edge like structures only, so that a line yields two edge responses. For this purpose, we use the phase asymmetry measure, which quantifies whether the signal is locally odd (edge like) or even (line like).

The monogenic signal is the isotropic and multidimensional extension of 1-D analytic signal [START_REF] Felsberg | The monogenic signal[END_REF]. The 3D monogenic signal is formed by the combination of the original signal and and three components that are Riesz transforms of the original signal [START_REF] Chenouard | 3D Steerable Wavelets in Practice[END_REF]: 𝑓 𝑚 (𝒙) = (𝑓(𝒙), 𝑅 1 𝑓(𝒙), 𝑅 2 𝑓(𝒙), 𝑅 3 𝑓(𝒙))

(1)

𝑹𝑓(𝒙) = (𝑅 1 𝑓(𝒙), 𝑅 2 𝑓(𝒙), 𝑅 3 𝑓(𝒙)) (2)
where 𝑓(𝒙) is the original image, 𝑹 is the Riesz transform and 𝒙 = (𝑥 1 , 𝑥 2 , 𝑥 3 ) ∈ ℝ 3 are the spatial coordinates. The Riesz operator component 𝑅 𝛼 : ℝ 3 ⟼ ℝ is characterized by its frequency response

where 𝝎 = (𝜔 1 , 𝜔 2 , 𝜔 3 ) ∈ ℝ 3 is the angular frequency variable conjugate to 𝒙. 𝛼 gives the 3D direction of the Riesz transformed components, corresponding to the directions of the basis vectors. ℱ is the Fourier transform operator. Since structure may be scale dependent, features at a useful scale have to be isolated, usually by the use of bandpass filters [START_REF] Bridge | An Introduction To The Monogenic Signal[END_REF]. Here we use log-Gabor filters. Thus, the 3D monogenic signal 𝑓 𝑚𝑔 (𝒙) is constructed as: where 𝑔(𝒙) is the spatial domain representation of the log-Gabor filter, and * denotes convolution operation. The 3D monogenic signal can also be expressed in terms of its even and odd components, which yields [START_REF] Kashif Rajpoot | Local phase based 3D boundary detection using monogenic signal and its application to real-time 3D echocardiography images[END_REF]:

𝑓 𝑚𝑔 (𝒙) = 𝑓 𝑚𝑔𝑒 (𝒙) + 𝑖 * 𝑓 𝑚𝑔𝑜 (𝒙) (5)
where the even component is

𝑓 𝑚𝑔𝑒 (𝒙) = 𝑓(𝒙) * 𝑔(𝒙), (6)
and the odd component

𝑓 𝑚𝑔𝑜 (𝒙) = (∑ |𝑅 𝛼 𝑓(𝒙) * 𝑔(𝒙)| 2 3 𝛼=1 ) 1/2 (7)
For the detection of edges in the image, we use the multiscale phase asymmetry measure, which is given by [START_REF] Bridge | An Introduction To The Monogenic Signal[END_REF]:

𝐴(𝒙) = ∑ ⌊|𝑓 𝑚𝑔𝑜,𝜆 𝑖 (𝒙)|-|𝑓 𝑚𝑔𝑒,𝜆 𝑖 (𝒙)|-𝑇⌋
√(𝑓𝑚𝑔𝑜,𝜆 𝑖 (𝒙)) 2 +(𝑓 𝑚𝑔𝑒,𝜆 𝑖 (𝒙)) 2 +𝜖 𝑖 [START_REF] Zhang | Periosteal Progenitor Cell Fate in Segmental Cortical Bone Graft Transplantations: Implications for Functional Tissue Engineering[END_REF] where T is a threshold for setting the sensitivity and suppressing noisy signals, ⌊. ⌋ is an operator which replaces negative values with zero, 𝜖 is small number to avoid division by zero, and {𝜆 𝑖 } are a set of center-wavelengths (𝜆 𝑖 =2𝜋 𝜔 𝑖 ⁄ ) of the log-Gabor filters.

RESULTS

Qualitative results of the 3D edge detection are shown in Fig. 3: (a) shows the original image, which is a zoom on a 2D slice from the 3D volume. The red arrows show low contrast interfaces between bone and vessel. Fig. 3 (b) shows the magnitude of the gradient. The boundary between bone and vessel can not be detected due to the low contrast. Fig. 3 (c) shows the multiscale 3D phase asymmetry of the original image using (λ=6,7,8). Lower scales contain too much edge details so that useful structures can not be recognized. Conversely, higher scales miss too many details. Edge detection at the low contrast interfaces between bone and vessel is substantially improved using the phase asymmetry, compared to the magnitude of the gradient. There are some spurious edges detected in homogeneous areas in the vessels. The contrast is relatively weak, however, and these edges, being inside solid areas, are always covered by a marker. Hence they do not influence the segmentation. Fig. 4 illustrates the relationship between markers and control surfaces, and shows segmentation using hysteresis and watershed using a gradient based control surface. overlayed on the monogenic phase asymmetry based control surface. Fig. 4 (d) shows the segmentation result using hysteresis thresholding, (e) shows the segmentation using the gradient based control surface watershed, and (f) shows the segmentation using the monogenic phase asymmetry based control surface watershed.

To validate the segmentation results, a representative volume (64×64×64 voxels) was manually segmented using VGStudio Max (Volume Graphics GmbH, Heidelberg, Germany).

Segmentation quality was measured using the Dice coefficient. The proposed method was compared to our previous protocol based on hysteresis thresholding as well as to marker-controlled watershed using a control surface based on the magnitude of the gradient. A comparison is shown in Tab. 1. The proposed method substantially improved segmentation in terms of the Dice coefficient. 3D renderings of a healthy sample and a sample with metastases segmented using the proposed method are shown in Fig. 5. The metastatic sample shows evident large bones lesions. Comparing (c) and (d), vessels in the metastatic sample appear thicker than in the healthy sample, which indicates an increased and abnormal vascularization, as expected.

Bone volume (BV) and vessel volume (VV) were measured by counting voxels in the corresponding compartments. Total volume (TV) of the samples was considered as the volume spanned by the outer contour of the bone. In order to compare the different samples, the normalized ratios: BV/TV, VV/TV were calculated. Mean vessel thickness (V.Th) was also extracted. Box plots of these parameters are shown in Fig. 6. Normal distribution of the parameters was verified using the Lilliefors test (p > 0.05). BV/TV was significantly lower in the metastatic group (ANOVA F-test, p < 0.05), which is consistent with the apparent large metastatic lesions visible in Fig. 5 (b). Although there was no significant difference in VV/TV (p = 0.30) and V.Th (p = 0.09), possibly due to the limited sample size, the box plots in Fig. 6 still indicate that VV/TV and V.Th is increased in the metastatic bone, meaning there are both more and thicker vessels in the metastatic bone.

CONCLUSION

We presented an algorithm to segment bone and vessels in 3D SR-μCT images of mouse bone with contrast agent, using the markercontrolled watershed on the monogenic signal phase asymmetry. The marker-controlled watershed permitted to address the problem of having vessels in contact to the bone surface, which precludes the correct use of hysteresis thresholding. In addition, the introduction of the monogenic signal phase asymmetry as control surface brought a substantial improvement of the segmentation of the relatively weakly contrasted bone and vessel interfaces. Our evaluation using manually segmented images showed substantial improvement in segmentation quality using the proposed method, compared to using hysteresis thresholding as well as using a gradient based control surface in the marker-controlled watershed. The algorithm was used to analyze healthy and metastatic bone. Several quantitative parameters were extracted to characterize bone and vessel structure. Statistical analysis revealed a decrease in bone volume, and indicated an increase in vessel volume and thickness in metastatic bone, as expected. The proposed protocol will be used to study the effects of anti-angiogenic drugs in bone metastasis.

Fig. 1 .

 1 Fig. 1. Schematic of the imaging setup.

Fig. 2 .

 2 Fig. 2. Generation of the marker image for the marker-controlled watershed.

Fig. 3 .

 3 Fig. 3. Illustration of 3D edge detection. (a) zoom in an original gray level image, (b) gradient magnitude image. (c) multiscale 3D asymmetry computed from the monogenic signal;

  𝑓 𝑚𝑔 (𝒙) = (𝑓(𝒙) * 𝑔(𝒙), 𝑅 1 𝑓(𝒙) * 𝑔(𝒙), 𝑅 2 𝑓(𝒙) * 𝑔(𝒙), 𝑅 3 𝑓(𝒙) * 𝑔(𝒙)) (4)

Fig. 4 .

 4 Fig. 4. Illustration of 3D segmentation results. (a) original image, (b) vessels marker (red), bone marker (green), background marker (blue) and the gradient based control surface. (c) markers and the monogenic signal phase based control surface. (d) segmentation result using hysteresis thresholding. (e) segmentation using gradient control surface based marker watershed. (f) segmentation with the proposed method using marker-controlled watershed and monogenic phase asymmetry.

Fig. 5 .

 5 Fig. 5. 3D rendering volumes. (a) bone and vessel rendering in healthy group. (b) bone and vessel rendering in metastases group. (c) vessel rendering in healthy group. (d) vessel rendering in metastases group.

TABLE 1 .

 1 Dice coefficient for the different methods

	Dice	Vessel	Bone	Background
	Hysteresis	0.8434	0.8055	0.8517
	thresholding
	Gradient +	0.9234	0.9127	0.9300
	watershed
	Monogenic +	0.9695	0.9240	0.9711
	watershed