Cécile, Erwan Manon Guillaume

Aussi

Cynthia Nino

Keywords: Snow avalanche numerical models, Bayesian inference, aggregated Sobol' indices, aggregated

Snow avalanche is a natural hazard defined as a snow mass in fast motion. Since the thirties, scientists have been designing snow avalanche models to describe snow avalanches. However, these models depend on some poorly known input parameters that cannot be measured. To understand better model input parameters and model outputs, the aims of this thesis are (i) to propose a framework to calibrate input parameters and (ii) to develop methods to rank input parameters according to their importance in the model taking into account the functional nature of outputs. Within these two purposes, we develop statistical methods based on Bayesian inference and global sensitivity analyses. All the developments are illustrated on test cases and real snow avalanche data.

First, we propose a Bayesian inference method to retrieve input parameter distribution from avalanche velocity time series having been collected on experimental test sites. Our results show that it is important to include the error structure (in our case the autocorrelation) in the statistical modeling in order to avoid bias for the estimation of friction parameters.

Second, to identify important input parameters, we develop two methods based on variance based measures. For the first method, we suppose that we have a given data sample and we want to estimate sensitivity measures with this sample. Within this purpose, we develop a nonparametric estimation procedure based on the Nadaraya-Watson kernel smoother to estimate aggregated Sobol' indices. For the second method, we consider the setting where the sample is obtained from acceptance/rejection rules corresponding to physical constraints. The set of input parameters become dependent due to the acceptance-rejection sampling, thus we propose to estimate aggregated Shapley effects (extension of Shapley effects to multivariate or functional outputs). We also propose an algorithm to construct bootstrap confidence intervals. For the snow avalanche model application, we consider different uncertainty scenarios to model the input parameters. Under our scenarios, the release avalanche position and volume are the most crucial inputs.

Our contributions should help avalanche scientists to (i) account for the error structure in model calibration and (ii) rank input parameters according to their importance in the models using statistical methods.

Résumé

Une avalanche de neige est un danger naturel défini comme une masse de neige en mouvement rapide. Depuis les années 30, scientifiques conçoivent des modèles d'avalanche de neige pour décrire ce phénomène. Cependant, ces modèles dépendent de certains paramètres d'entrée mal connus qui ne peuvent pas être mesurés. Pour mieux comprendre les paramètres d'entrée du modèle et les sorties du modèle, les objectifs de cette thèse sont (i) de proposer un cadre pour calibrer les paramètres d'entrée et (ii) de développer des méthodes pour classer les paramètres d'entrée en fonction de leur importance dans le modèle en tenant compte la nature fonctionnelle des sorties. Dans ce cadre, nous développons des méthodes statistiques basées sur l'inférence bayésienne et les analyses de sensibilité globale. Nos développements sont illustrés sur des cas de test et des données réelles des avalanches de neige. D'abord, nous proposons une méthode d'inférence bayésienne pour récupérer la distribution des paramètres d'entrée à partir de séries chronologiques de vitesse d'avalanche ayant été collectées sur des sites de test expérimentaux. Nos résultats montrent qu'il est important d'inclure la structure d'erreur (dans notre cas l'autocorrélation) dans la modélisation statistique afin d'éviter les biais dans l'estimation des paramètres de frottement.

Deuxièmement, pour identifier les paramètres d'entrée importants, nous développons deux méthodes basées sur des mesures de sensibilité basées sur la variance. Pour la première méthode, nous supposons que nous avons un échantillon de données et nous voulons estimer les mesures de sensibilité avec cet échantillon. Dans ce but, nous développons une procédure d'estimation non paramétrique basée sur l'estimateur de Nadaraya-Watson pour estimer les indices agrégés de Sobol. Pour la deuxième méthode, nous considérons le cadre où l'échantillon est obtenu à partir de règles d'acceptation/rejet correspondant à des contraintes physiques. L'ensemble des paramètres d'entrée devient dépendant du fait de l'échantillonnage d'acceptation-rejet, nous proposons donc d'estimer les effets de Shapley agrégés (extension des effets de Shapley à des sorties multivariées ou fonctionnelles). Nous proposons également un algorithme pour construire des intervalles de confiance bootstrap. Pour l'application du modèle d'avalanche de neige, nous considérons différents scénarios d'incertitude pour modéliser les paramètres d'entrée. Dans nos scénarios, la position et le volume de départ de l'avalanche sont les entrées les plus importantes.

Nos contributions peuvent aider les spécialistes des avalanches à (i) prendre en compte la structure d'erreur dans la calibration du modèle et (ii) proposer un classement des paramètres d'entrée en fonction de leur importance dans les modèles en utilisant des approches statistiques.

Mot-clés : Modèles numériques d'avalanche de neige, inférence bayésienne, indices de Sobol' agrégés, effets de Shapley agrégés, réduction de la dimension v 1

Introduction

In this chapter, we give a general overview of the thesis, of the scientific questions to be addressed in this work and of the way we handled them. The concepts and tools employed are detailed in Chapter 2.

Context

Snow covered mountains lie in apparent equilibrium till an avalanche suddenly occurs. Snow avalanche is a natural hazard defined as a snow mass in fast motion on a slope [START_REF] Ancey | Snow Avalanches[END_REF]. In fact, avalanche velocities can reach 60 ms -1 [START_REF] Schaer | Particle densities, velocities and size distributions in large avalanches from impact-sensor measurements[END_REF]. There exist many factors that interact to explain the release of an avalanche for example, the average inclination of starting zones (between 27 • to 50 • ), the ground surface roughness, snowpack structure or temperature changes [START_REF] Ancey | Snow Avalanches[END_REF]. But more important than its causes are its consequences: from building damage or destruction to casualties. Around 100 people die each year in the European Alps due to snow avalanches [START_REF] Techel | Avalanche fatalities in the european alps: long-term trends and statistics[END_REF]. To try in someway to prevent these accidents, a better comprehension of avalanche phenomena might be crucial.

Since the thirties, scientists have been designing and improving snow avalanche models (e.g., physical models, probabilistic models) to describe and understand better snow avalanches. Models consists in a simplification of natural phenomena and they often depend on some poorly known input parameters that cannot be necessarily measured. To cite an example, friction parameters (representing the behavior of snow in motion) cannot be measured from in situ experiments. That's why friction parameter meaning and their link to physical properties of snow is still debated in the literature of the field [START_REF] Naaim | Dense avalanche friction coefficients: influence of physical properties of snow[END_REF].

To go further, the following questions may be addressed. First, if we are capable to measure some other attributes of an avalanche as e.g., its velocity, how to calibrate the unobserved parameters? Within this purpose, avalanche data (e.g, from photogrammetric techniques) with high temporal resolution (1 s) have been collected in experimental test sites to study snow avalanches [e.g., the French test site at Lautaret Pass [START_REF] Thibert | The full-scale avalanche test-site at Lautaret Pass (French Alps)[END_REF]. One solution to calibrate model parameters given a set of observations is to solve an inverse problem, which consists into finding the parameter configuration(s) that approach the best reality, where reality is known through the observations. Second, which parameters should be calibrated with high accuracy and which others might not? The purpose of this thesis is to answer these two questions using a set of advanced statistical techniques.

To answer the first question, we use Bayesian inference. Bayesian inference is a powerful tool to solve inverse problems. Its main advantage is its application even if sample size is small because prior knowledge, given by scientists, can be provided to improve parameter inference. Many scientists from different fields (e.g., social and environmental sciences) are turning to Bayesian inference to solve problems where samples are scarce. In particular, avalanche observations are scarce since avalanches are rare events. Thus, Bayesian inference is the ideal framework for avalanche models and indeed, it is becoming more and more popular in the avalanche community. In particular, the velocity avalanche profile we disposed has 21 observations. For this reason, we propose to answer the first question using a Bayesian framework to calibrate friction parameters for snow in motion. For the snow avalanche field, the novelty of our approach is that we consider the error (defined as observationsmodel simulations) structure in the calibration. Indeed, our observed velocities are a time series and the errors are autocorrelated. Therefore within our results, we demonstrate that it is important to consider the correct structure of the errors in the calibration framework because if errors are not properly modeled, estimation of friction parameters might be biased.

To answer the second question: how to determine which input parameters contribute the most to model output uncertainty, we use global sensitivity analysis. Global sensitivity analysis (GSA) is a probabilistic technique which consists in modeling unknown input parameters by a probability distribution which propagates through the model to the outputs. Input parameter distributions are often guided by practitioner's belief. Then, input parameters are ordered according to their contribution on the model outputs by computing sensitivity measures. By answering this second question, we gain (i) a better comprehension of avalanche model and, (ii) we might put more efforts to calibrate important input parameters and, fix the other inputs to literature nominal values without losing model output accuracy. In fact, a major drawback of Bayesian inference is the high number of model simulations required to get algorithms convergence. Indeed, at least 1000 simulations might be required to get convergence. Moreover, algorithm convergence could depend on the number of input parameters. Thus, the number of model simulations required could be reduced by identifying influential inputs.

From a GSA scientist's perspective, avalanche model outputs offer an interesting application framework because (i) they are functional (e.g., the functional flow velocity) and, (ii) their input/output samples could be selected by acceptance-rejection (AR) algorithms. For example, we could use only avalanche simulations with certain particular characteristics e.g., having a volume superior to a prescribed threshold (e.g., 7000m 3 ) and, runout distance superior to a prescribed position. In that context, independence of input parameters is no longer true and the computation and interpretation of sensitivity measures is a difficult task. Thus, the functional output nature and dependence between input parameters were our guidelines for proposing a framework for GSA in our context of snow avalanches.

In GSA literature, there exist different sensitivity measures. We focus in two variancebased sensitivity measures: Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] and Shapley effects [START_REF] Shapley | A Value for n-Person Games[END_REF]. Sobol' indices are sensitivity measures commonly used. They can be applied if inputs are independent. This sensitivity measure is easy to interpret because, (i) all the indices are normalized between 0 and 1 and (ii) sum to 1, (iii) a high index shows the input is relevant and a value close to zero shows that it is not. When inputs are dependent, distributions are chosen as the posterior of the Bayesian inference. This work was submitted to Reliability Engineering & System Safety journal in March 2020 and presented in Chapter 4:

Heredia MB, Prieur C, Eckert N. Nonparametric estimation of aggregated Sobol' indices: application to a depth averaged snow avalanche model. Submitted in March 2020.

Eventually, we consider the setting where the sample is obtained from an acceptance/rejection rule to physical constraints. The corresponding set of input parameters present dependencies thus, we propose to estimate Shapley effects extended to multivariate or functional model outputs. We call these new sensitivity measures aggregated Shapley effects. More precisely, we propose an algorithm to estimate aggregated Shapley effects in a given data framework based in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF]. We also build bootstrap confidence intervals, adapting the ideas of [START_REF] Benoumechiara | Shapley effects for sensitivity analysis with dependent inputs: bootstrap and kriging-based algorithms[END_REF]. This work has been submitted to the SIAM/ASA Journal on Uncertainty Quantification in July 2020 and presented in Chapter 5:

Heredia MB, Prieur C, Eckert N., Given data inference of aggregated Shapley effects, application to avalanche modeling. Submitted in July 2020.

Organization of the manuscript

The manuscript is organized as follows. In Chapter 2, we propose a short introduction to (i) 1-D avalanche modeling, (ii) Bayesian inference applied to parameter model calibration and (iii) global sensitivity analysis with a particular emphasis when inputs are independent or dependent, and multivariate or functional. This chapter can be seen as a toolbox to understand the developments we have made. Then, the manuscript is divided in two parts:

Part I: Bayesian Inference. Part I is composed by Chapter 3. In Chapter 3, we propose a Bayesian framework to calibrate an avalanche model using velocities extracted from images.

Part II: Global Sensitivity Analysis

In Part II, we develop global sensitivity methods applied to dynamic models, more particular to a 1-D avalanche model. Part II is composed by Chapter 4 and Chapter 5. In Chapter 4, we proposed a non parametric method to estimate aggregated Sobol' indices from an i.i.d. sample and, in Chapter 5, we extend Shapley effects to multivariate and functional outputs, we also propose algorithms for the estimation and the computation of bootstrap confidence intervals.

Finally, conclusions and perspectives for this thesis are exposed in Chapter 6.

2

Avalanche models, Bayesian inference and global sensitivity analysis Overview: In this chapter, we introduce the principal tools and concepts employed in the core chapters of this manuscript (Chapters 3 to 4). A short introduction to snow avalanche dynamic modeling is presented in Section 2.1. We put particular attention to previous studies of model calibration and sensitivity analysis in the avalanche field. Then, model calibration using Bayesian inference is explained in Section 2.2 and global sensitivity analysis methods are described in Section 2.3. We put emphasis on global sensitivity methods that can be applied to models with dependent input parameters and functional model outputs.

A short introduction to 1-D snow avalanche dynamic modeling

Avalanche numerical models provide a simplification of avalanche flow on the basis of physical conservation laws. To situate avalanche models in human history, one of the first avalanche models was probably proposed in the 1930s [START_REF] Salm | A short and personal history of snow avalanche dynamics[END_REF]. This first model had a dry friction coefficient (Coulomb friction) [START_REF] Salm | A short and personal history of snow avalanche dynamics[END_REF]. Then, [START_REF] Voellmy | On the destructive force of avalanches[END_REF] proposed a based on hydrodynamics model, that had two friction coefficients: a turbulent friction coefficient ξ and a Coulomb-like dry coefficient µ. Since then, different analytical and numerical models have been proposed based in Coulomb or Voellmy friction laws [e.g., [START_REF] Perla | A Two-Parameter Model of Snow-Avalanche Motion[END_REF][START_REF] Bruno Salm | Berenchnung von Fliesslanwinen, eine Anleitung für Paktiker mit Beispielen[END_REF]. Analytical (e.g., sliding block models) and numerical models (e.g., Saint-Venant proposed by [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF]) have advantages and disadvantages as for example their accuracy in test site experiments, computational cost, number of parameters, etc. Also, some of them may be only capable to output flow velocities and runout distances1 while others can also compute flow depths. Nevertheless, avalanche modeling is still challenging, due to lack of measurements and due to avalanche complexity (e.g., flow regime transitions, snow mass balance) [START_REF] Harbitz | Conclusions from a recent survey of avalanche computational models[END_REF]]. However, depending on the application, analytical or numerical models (or both) might be used. A more detailed review of these models can be found in e.g., [START_REF] Harbitz | Conclusions from a recent survey of avalanche computational models[END_REF], [START_REF] Salm | A short and personal history of snow avalanche dynamics[END_REF] and [START_REF] Jamieson | Application and limitations of dynamic models for snow avalanche hazard mapping[END_REF]. For illustration purpose, Figure 2.1 shows a graphical simplification of a sliding block model. These kind of models represent an avalanche as a sliding rigid body experimenting frictional forces. The model is based in the motion equation:

du dt = g sin θ - F m , (2.1)
where u = || u|| is the velocity norm, g is the gravity constant, θ is the local slope, h is the mean flow depth, m the mass and F = || F || the friction force. For example, Coulomb friction law F = µmg cos θ or Voellmy friction law F = µmg cos θ + mg ξh u2 could be considered in Equation (2.1).

Model input parameters

All of the aforementioned models, independently of their complexity, depend on unknown parameters that cannot be measured. More specifically, friction parameters cannot be measured from experiments. To guide practitioners, some authors have created parameter guidelines. For example, [START_REF] Bruno Salm | Berenchnung von Fliesslanwinen, eine Anleitung für Paktiker mit Beispielen[END_REF] proposed reference values for µ according to the avalanche characteristics (e.g., wet snow, large avalanche) and for ξ according to the avalanche path characteristics (e.g., path local slope). Author guideline is called the Swiss guideline. For example, Christen et al. [2010b] applied the guideline recommendations for back calculation of the In den Arlen 2 avalanche using the 2-D RAMMS model. However, parameter in situ calibration is still necessary for a better comprehension of avalanche models. Deterministic inversion, which consists in writing friction coefficients as a function of observed quantities, has been employed to calibrate parameters from in situ experiments. In the work of [START_REF] Dent | Density, velocity and friction measurements in a dry-snow avalanche[END_REF], friction coefficients were computed by measuring normal forces. [START_REF] Ancey | Inverse problem in avalanche dynamics models[END_REF] proposed a deterministic inversion framework to retrieve friction parameters from runout distances and the sum of snow fall amount during the three days previous to the avalanche event. Indeed, the snow fall during the days preceding the avalanche event is a key meteorological parameter in avalanche phenomena [START_REF] Ancey | Inverse problem in avalanche dynamics models[END_REF]. In [START_REF] Ancey | Estimating bulk rheological properties of flowing snow avalanches from field data[END_REF] a back analysis was performed to infer friction coefficients from velocities. [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF] performed a back analysis to model validation, also the authors found an empirical relationship between the dry friction coefficient µ and runout distances. Recently, [START_REF] Oda | Inverse simulation for extracting the flow characteristics of artificial snow avalanches based on computation fluid dynamics[END_REF] performed small scale 2.1. A short introduction to 1-D snow avalanche dynamic modeling experiments using a rotating drum device to measure flow velocities and flow depths in order to perform a deterministic inversion to retrieve friction parameters. Nevertheless, parameter point estimation does not take into account parameter uncertainty which is important for uncertainty propagation.

One of the principal difficulties of parameter calibration in avalanche problems is the low number of observations recorded. Therefore, more complex schemes of calibration based on Bayesian inference have been proposed. For example, [START_REF] Ancey | Monte Carlo Calibration of Avalanches Described as Coulomb Fluid Flows[END_REF], [START_REF] Straub | A Bayesian probabilistic framework for avalanche modelling based on observations[END_REF] and [START_REF] Eckert | Long-term avalanche hazard assessment with a Bayesian depth averaged propagation model[END_REF] used runout distances to model calibration. [START_REF] Gauer | Optimization and probabilistic calibration of avalanche block models[END_REF] used the velocity at the lower avalanche track and runout distances to model calibration. Also, procedures based on optimization have been proposed. For example, [START_REF] Naaim | Dense avalanche friction coefficients: influence of physical properties of snow[END_REF] performed the calibration of friction parameters using runout distances from historical avalanche events. Regarding 3-D avalanche model calibration, in [START_REF] Fischer | Evaluation of probabilistic snow avalanche simulation ensembles with Doppler radar observations[END_REF], the 3-D SamosAT (Snow Avalanche Modeling and Simulation Advanced technology) model is calibrated using flow depths and velocities deduced from Doppler radar measurements. [START_REF] Fischer | Multivariate parameter optimization for computational snow avalanche simulation[END_REF] proposed a multivariate parameter optimization to calibrate SamosAT.

In the last two decades, the interest in Bayesian approaches is increasing in the avalanche community. To demonstrate this assertion, we cite some examples: the already described works of [START_REF] Ancey | Monte Carlo Calibration of Avalanches Described as Coulomb Fluid Flows[END_REF], [START_REF] Straub | A Bayesian probabilistic framework for avalanche modelling based on observations[END_REF], [START_REF] Gauer | Optimization and probabilistic calibration of avalanche block models[END_REF], [START_REF] Eckert | Long-term avalanche hazard assessment with a Bayesian depth averaged propagation model[END_REF], [START_REF] María | Bayesian calibration of an avalanche model from autocorrelated measurements along the flow: application to velocities extracted from photogrammetric images[END_REF] and other recent works, [START_REF] Schläppy | Validation of extreme snow avalanches and related return periods derived from a statistical-dynamical model using tree-ring techniques[END_REF] reconstructed runout distances from dendropomorphic historical data, then authors applied a Bayesian approach based on [START_REF] Eckert | Long-term avalanche hazard assessment with a Bayesian depth averaged propagation model[END_REF] to model calibration and to estimate return periods. Recently, [START_REF] Fischer | Bayesian inference in snow avalanche simulation with r.avaflow[END_REF] provided a comprehensive tool for snow avalanche simulation using Bayesian inference which is implemented in the open source software r.avaflow. The software r.avaflow was developed by [START_REF] Mergili | Pudasaini. r.avaflow v1, an advanced open-source computational framework for the propagation and interaction of two-phase mass flows[END_REF].

Sensitivity analyses for avalanche models

As mentioned before, avalanche models depend on poorly known parameters (e.g., friction parameters and initial conditions corresponding to the avalanche released). Authors have studied the parameter sensitivity, for example, [START_REF] Barbolini | Estimate of uncertainties in avalanche hazard mapping[END_REF] used a Monte Carlo approach to analyze the sensitivity for the runout distances and impact pressures outputs of the VARA model. Authors found two interesting results: (i) for relatively frequent events, friction parameters and release depth and length are important for runout distances and, (ii) for extreme rare events, the friction coefficient has a higher importance than the two other parameters. [START_REF] Borstad | Sensitivity analyses in snow avalanche dynamics modeling and implications when modeling extreme events[END_REF] developed a sensitivity analysis of an avalanche model with Coulomb friction law. Their results demonstrated that the friction coefficient is important for the runout positions. [START_REF] Bühler | Sensitivity of snow avalanche simulations to digital elevation model quality and resolution[END_REF] studied the influence of the digital elevation model resolution (DEM) on the outputs of the RAMMS model. They found that DEM resolution and quality is critical for modeling and a spatial resolution of around 25m is sufficient for large-scale modeling. [START_REF] Fischer | Evaluation of probabilistic snow avalanche simulation ensembles with Doppler radar observations[END_REF] performed a sensitivity analysis for RAMMS model by estimating the Spearman rank correlation between parameters and outputs. [START_REF] Buhler | Sensitivity Analysis of the RAMMS Avalanche Dynamics Model in a Canadian Transitional Snow Climate[END_REF] made a sensitivity analysis of runout distance, volume and avalanche velocity in the runout zone with respect to the initial volume in the RAMMS model. They found that the released volume has a large effect in runout distance and avalanche velocity. However, even if all the previous Chapter 2. Avalanche models, Bayesian inference and global sensitivity analysis works provide useful information about avalanche models, none of them has used variancebased global sensitivity measures to quantify input parameter importance.

Few words about land-use planning

Finally, why avalanche input parameters and thus, model outputs are important? The reason is that model outputs (e.g., avalanche speed and runout distance) might be used in land use planning which typically consists into dividing a region in zones according to their potential risk [START_REF] Salm | A short and personal history of snow avalanche dynamics[END_REF]. In some countries (e.g., Austria, France, Switzerland), this assignation is done by colors: green (or white in France [Eckert et al., 2018b]), blue and red. The colors assignation are sometimes determined by runout return periods. A runout return period is the mean time in which a given runout distance is reached or exceeded at a given path's position [START_REF] Eckert | Revisiting statistical-topographical methods for avalanche predetermination: Bayesian modelling for runout distance predictive distribution[END_REF][START_REF] Eckert | Bayesian stochastic modelling for avalanche predetermination: from a general system framework to return period computations[END_REF][START_REF] Schläppy | Validation of extreme snow avalanches and related return periods derived from a statistical-dynamical model using tree-ring techniques[END_REF]. For example in Switzerland, red zones have return periods lower than 300 years and avalanche pressures can exceed 30kNm -2 [WSL]. Constructions are forbidden in red zone [WSL]. Therefore, parameter estimation and associated uncertainty might be useful for computing more accurate avalanche simulations and as a consequence, for a better land use planning. The importance of avalanche model calibration and current issues for hazard mapping are nicely described in [START_REF] Jamieson | Application and limitations of dynamic models for snow avalanche hazard mapping[END_REF].

To summarize, the two main axes of this thesis are model calibration using Bayesian inference and global sensitivity analysis applied to avalanche models. These two subjects are explained in Section 2.2 and Section 2.3.

Model calibration using Bayesian inference

Model calibration consists in finding the most likely combination(s) of parameters given some observations. We are not only interested in determining an optimal single combination of parameters, but also in the uncertainty associated to it [START_REF] Kuczera | Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm[END_REF]. Bayesian calibration offers a solution to this problem, because it allows to assess parameter uncertainty by deriving the posterior distribution of parameters from observations [START_REF] Kennedy | Bayesian calibration of computer models[END_REF]. Moreover, advantages of Bayesian calibration are twofold: (i) if model is used for prediction, all the sources of uncertainty are included, and (ii) since no model is perfect, the approach cope for any inadequacy between observations and model [START_REF] Kennedy | Bayesian calibration of computer models[END_REF].

Let us denote f a prediction model for a given phenomenon (e.g., an avalanche model which predicts the avalanche runout distance), f depends on parameters Θ ∈ R p and some forcing variables x ∈ R d . Let us denote by y obs a vector of n noisy observations:

y obs i = f (Θ, x) + i , ∀i ∈ {1, . . . , n} (2.2)
where models observation errors. For example, y obs could be a set of observed runout distances to calibrate an avalanche model. Note that within this formulation, we do not make distinction between the observation errors and the model inadequacy. The distinction has been studied in, e.g., [START_REF] Kennedy | Bayesian calibration of computer models[END_REF], [START_REF] Carmassi | Bayesian calibration of a numerical code for prediction[END_REF]. In a statistical framework to solve the inverse problem of retrieving the model parameters Θ, assumptions of the distribution of the residuals or errors might be required. For example, they could be assumed i.i.d. centered Gaussian random variables i ∼ N (0, σ 2 ) ∀i ∈ {1, . . . , n}. Let us denote γ the residual parameters. In our example, γ = {σ 2 }. An advantage of Bayesian inference is that practitioner experience or beliefs could be plugged into the inference by the prior distribution denoted π(Θ, γ). Then, the posterior distribution is: π(Θ, γ|y obs , x) = L(y obs |Θ, γ, x)π(Θ, γ) L(y obs |Θ, γ, x)π(Θ, γ)dΘdγ

(2.3)

where L(y obs |Θ, γ, x) is the likelihood of the observations which can be obtained from the error distribution. It is usually difficult to compute L(y obs |Θ, γ, x)π(Θ, γ)dΘdγ, except when there is an analytical expression of the posterior distribution, which is achieved by using conjugate priors. When there is no explicit posterior formulation, Monte Carlo Markov chain (MCMC) methods can be used to sample from the posterior distribution. MCMC methods consist into constructing a stationary and ergodic Markov chain that converges, under mild conditions, to the posterior distribution π(Θ, γ|y obs , x).

Among the MCMC methods, there is the Metropolis-Hasting algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo Sampling Methods Using Markov Chains and Their Applications[END_REF] which constructs sequentially a Markov chain by applying acceptance-rejection rules defined by a proposal distribution q. The convergence of the MH algorithm depends on q [START_REF] Robert | The Metropolis-Hastings algorithm[END_REF] and q is usually assumed as multinormal distribution.

If there is an analytical expression of all the parameter conditional distributions π(θ |y obs , x, θ 1 , . . . θ -1 , θ +1 , . . . θ d ), Gibbs sampling, a particular instance of MH algorithm, could be applied. The Hamiltonian Monte Carlo [START_REF] Simon Duane | Hybrid monte carlo[END_REF] algorithm, also known as hybrid Monte Carlo, has shown good results in applications. Nevertheless, the gradient of the distribution functions are required and their computations might be expensive.

When observations are time series, errors could be autocorrelated or/and suffer from heteroscedasticity. In the calibration of hydrological models, authors have studied the importance of considering the correlations and/or heteroscedasticity in residuals [e.g., [START_REF] Kuczera | Improved parameter inference in catchment models: 1. Evaluating parameter uncertainty[END_REF][START_REF] Kuczera | Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm[END_REF][START_REF] Evin | Comparison of joint versus postprocessor approaches for hydrological uncertainty estimation accounting for error autocorrelation and heteroscedasticity[END_REF]. Indeed, they found that if autocorrelation is not included parameter estimation could be biased. In Chapter 3, given an avalanche velocity measured time series, we analyze the impact of considering or not the error autocorrelation in a Bayesian inference calibration framework. A block model based on the Voellmy friction law is considered. [START_REF] Saltelli | Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models[END_REF] defined precisely the sensitivity analysis (SA) as "the study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input"3 . Under a SA framework, we consider a deterministic model f (computational code, black box model, etc) that depends on d inputs X = (X 1 , . . . , X d ). First for simplicity, we consider the output Y as scalar Y = f (X). The inputs are modeled by random variables on a probability space (Ω, F, P). Thus, the output Y is also a random variable (Ω, F, P) because it depends on the random vector X. Let us denote (x, y) a random sample from (X, Y ). Figure 2.2 illustrates the SA framework. The following exposition is inspired by the works of [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF] and [START_REF] Iooss | Introduction to sensitivity analysis[END_REF].

Global sensitivity analysis

SA is divided in two classes: local and global methods. Local methods consist in studying the impact of small input perturbations around certain nominal values on the output. An example of local methods are partial derivatives, they consist into calculating ∂f ∂x i at a given point x 0 = (x 0 1 , . . . , x 0 d ) in the domain of X. Local methods are sometimes used when model computational cost or number of parameters is high. For example, [START_REF] Castaings | Sensitivity analysis and parameter estimation for distributed hydrological modeling: potential of variational methods[END_REF] proposed the adjoint state method to compute efficiently the derivatives of the output of a distributed hydrological model with respect to parameter inputs to perform local SA. [START_REF] Castillo | A general method for local sensitivity analysis with application to regression models and other optimization problems[END_REF] proposed a local SA method based on the duality property of mathematical programming to compute partial derivatives. In contrast to local methods, global methods explore all the input domain and not only variations around nominal values.

Global methods provide a more convincing tool for highly nonlinear models often encountered in environmental sciences. Indeed, global methods provide an exploration of the whole input space [Saltelli et al., 2008, Iooss and[START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF]. Nonetheless, there is a price to pay to use global sensitivity methods in terms of number of model simulations. But there are some ways to reduce computational costs which will be explained in Subsection 2.3.5. [START_REF] Tang | Comparing sensitivity analysis methods to advance lumped watershed model identification and evaluation[END_REF] compared local and global methods applied to a watershed model and they concluded that global methods yield robust results. [Saltelli et al., 2008, page 11] recommended strongly to use global methods rather than local ones. Hereafter, we focus on global sensitivity analysis methods.

Global sensitivity analysis (GSA) consists into identifying the input parameters that contribute the most (in a sense to be precised) to a given quantity of interest defined from the output of the model [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF]. The objectives of GSA are numerous: they are to provide a better comprehension of the numerical model, to identify non influential parameters to set them to nominal values, among others [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF]. To guide practitioners, [START_REF] Saltelli | Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models[END_REF] proposed four settings of GSA objectives: factor prioritization, factor fixing, variance cutting and factor mapping. Before performing a GSA, it is important to have clear the objectives of the study [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF] and also, to select the adequate GSA method to avoid misinterpretations [START_REF] Saltelli | Why so many published sensitivity analyses are false: A systematic review of sensitivity analysis practices[END_REF]. Then, parameter inputs might be ordered according to their importance or contributions by computing sensitivity measures.

GSA has been widely applied in many fields as for example, in environmental sciences and engineering. We cite some recent works in different application domains, [START_REF] Sarrazin | Global sensitivity analysis of environmental models: Convergence and validation[END_REF] proposed a guide to apply GSA in environmental models, more particularly applied to hydrological models. [START_REF] Mavromatidis | Uncertainty and global sensitivity analysis for the optimal design of distributed energy systems[END_REF] developed a GSA for distributed energy systems. These models are important for an efficient and sustainable energy future. [START_REF] Prieur | A global sensitivity analysis approach for marine biogeochemical modeling[END_REF] performed a GSA for a marine biogeochemical model with a high number of inputs (74 input parameters).

The GSA measures may be divided into: screening approaches, variance-based measures, moment free measures and, other measures. In the manuscript, we only focus on variance based measures which are explained in Subsection 2.3.2. Even if we do not use the other three GSA measures, they are briefly explained in Subsections 2. 3.1, 2.3.3 and, 2.3.4 to provide a general overview of GSA measures.

Screening approaches

Screening approaches consist in the discretization of inputs X in levels [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF]. If the model has a large number of inputs d 1, these sensitivity measures could help the user to identify non influential input parameters. Therefore, screening methods could be a preliminary step before computing other, in general more expensive, GSA measures when input dimension is high. One of the most popular screening methods is the OAT (one at a time) Morris method [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF].

Here we describe briefly the Morris method. The d random inputs X are supposed independent and are modeled with uniform distributions on [0, 1], i.e. X i ∼ U[0, 1] for all i ∈ {1, . . . , d}. The assumption of uniformity on [0, 1] is not constraining as for X i with cumulative distribution function (cdf) F , we can apply the inverse transform sampling to come back to that setting.

The rectangular domain [0, 1] of each input is divided in a p-level grid, then each x i may take values from {0, 1/(p -1), 2/(p -1), . . . , 1}. Let Ω = {0, 1/(p -1), 2/(p -1), . . . , 1} d . Let ∆ be a predetermined multiple of 1/(p -1). For a given value x ∈ ω, the elementary effect of the input X i is defined according to:

d i (x) = f (x 1 , x 2 , . . . , x i-1 , x i + ∆, x i+1 , . . . , x d ) -f (x 1 , . . . , x d ) ∆
x must be such that x i + ∆ remains in {0, 1/(p -1), 2/(p -1), . . . , 1}. We sample n points x 1 , . . . , x n from X on Ω such that x i + ∆ is in {0, 1/(p -1), 2/(p -1), . . . , 1}. Then, we compute d i (x ) for each ∈ {1, . . . , n}. Finally, we can estimate the mean absolute value μi and standard deviation s i of the elementary effects. By calculating the indices μi and s i , the Morris method allows us to identify three kind of inputs (i) non influential inputs (μ i will have a low value), (ii) influential inputs with linear effects (μ i will be large and s i will be small) and, (iii) influential inputs with non-linear and/or interactions effects (μ i and s i will both be large). Graphical 2-D representation of (μ i , s i ) points might help the user to classify the inputs [see [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF]]. Here we will not explain the details of the Morris measure estimation (e.g., the structure of the design of experiment) and parameter tuning (e.g., ∆, p, the threshold for defining small or large ûi and s i ) since we do not use them in the core chapters of the manuscript, but more information can be found in the works of [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF] or [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF], the last work proposed an improved sampling strategy.

Example We illustrate the method in the Morris function [modified example from Morris, 1991]. The Morris function has d = 20 independent inputs X i ∼ U[0, 1] ∀i ∈ {1, . . . , 20} and is constructed as follows: (2.4) where w i = 2(X i -1/2) except for i = 3, 5 and 7 where w i = 2(1.1X i /(X i + 0.1) -1/2). The β coefficients are defined as:

y = β 0 + 20 i=1 β i w i + 20 i<j β ij w i w j + 20 i<j< β ij w i w j w + 20 i<j< <s β ij s w i w j w w s ,
β i = 20, i = 1, . . . , 7, β 8 = 25, β 9 = 30, β 10 = 40, β ij = -15, i, j = 1, . . . , 6, β ij = -10, i, j, = 1, . . . , 5, β ij s = 5,
i, j, , s = 1, . . . , 4, the remainder coefficients are set to 0. Figure 2.3 displays (μ i , s i ) with i ∈ {1, . . . , 20}. The standard error of the mean which is estimated as s i / √ n (result from the central limit theorem) allow us to display confidence intervals for the mean corresponding to µ ± 2s/ √ n (black line). Within this line, we can deduce the following results: inputs X 11 to X 20 are negligible because their means and standard deviations are 0. Inputs X 1 to X 10 are influential. More particularly, inputs X 8 , X 9 and X 10 have large means and small standard deviations which suggest that they are influential with linear effects. And finally, the remaining inputs (X 3 to X 10 ) have non linear and/or interaction effects. As we have seen in this example, the Morris method gives useful insights about input importance. However, to rank unambiguously the inputs according to their importance, we will need other sensitivity measures and they are explained in the following section.

Variance-based measures

Variance-based measures analyze how uncertainty in inputs impacts the variance of the output(s). The two main variance-based measures are Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] and Shapley effects [START_REF] Owen | Sobol' Indices and Shapley Value[END_REF]. Sobol' indices are based on a functional ANOVA decomposition [START_REF] Efron | The jackknife estimate of variance[END_REF], which is unique if inputs are independent, and Shapley effects are based on the cooperative game theory concept of Shapley value [START_REF] Shapley | A Value for n-Person Games[END_REF]. The main difference between both measures is that Shapley effects are easier to interpret if the inputs are not independent.

Sobol' indices

If the output of a model is scalar, its total variance can be split into partial variances by using the Hoeffding decomposition [START_REF] Hoeffding | A Class of Statistics with Asymptotically Normal Distribution[END_REF]. If the inputs are independent, each of the partial variances is associated to an input or to an interaction between inputs through the functional ANOVA decomposition [START_REF] Efron | The jackknife estimate of variance[END_REF]. Then, the Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] are calculated as the ratio of each partial variance to the total variance. In mathematical terms, this is expressed as follows: let us suppose the random variables X 1 , . . . , X d are independent from each other and f ∈ L 2 (P X ) where P X = d i=1 P X i . Thus, f can be decomposed as [START_REF] Hoeffding | A Class of Statistics with Asymptotically Normal Distribution[END_REF]:

f (X) = f ∅ + d i=1 f i (X i ) + 1≤i< ≤d f i (X i , X ) + . . . + f {1,...,d} (X 1 , . . . , X d ), (2.5)
where

f ∅ = E(Y ) and f u (X u ) = E(Y |X u ) - v u f v (X v ) = v u (-1) |u|-|v| E(f (X)|X i , i ∈ v).
This decomposition is the so-called functional ANOVA decomposition. From (2.5) and from the fundamental independence assumption we get [START_REF] Efron | The jackknife estimate of variance[END_REF]:

Var(Y ) = d i=1 Var(f i (X i )) + 1≤i< ≤d Var(f i (X i , X )) + . . . + Var(f {1,...,d} (X 1 , . . . , X d )).
(2.6) Then, Sobol' indices are defined:

• First-order Sobol' index of input X i :

S i = Var(f i (X i )) Var(f (X)) = Var(E(Y |X i )) Var(f (X)) .
(2.7)

• Second-order Sobol' index of input X i and X j , to quantify interactions:

S i = Var(E(Y |X i , X )) Var(f (X)) -S i -S .
(2.8)

• Total Sobol' index of input X i to quantify the effect of X i and all its interaction with the other inputs:

S T i = 1 - Var(E(Y |X 1 , . . . , X i-1 , X i+1 , . . . , X d )) Var(f (X)) .
(2.9) q q q q q q q q q q q q q q q q q q q q 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 X 2 X 1 a) q q q q q q q q q q q q q q q q q q q q 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8

X 2 X 1 b) Figure 2.4: n = 20 samples of X 1 , X 2 random variables uniformly distributed on [0, 1] a) struc-
tured sample design required to estimate Sobol' indices using [START_REF] Tissot | A randomized Orthogonal Array-based procedure for the estimation of first-and second-order Sobol' indices[END_REF] method. b) A random sample.

Higher order indices might be also defined (e.g., S i 1 ,i 2 ,i 3 with i 1 , i 2 , i 3 ∈ {1, . . . , d}) but in general, only the aforementioned indices are computed in applications. Sobol' indices interpretation is easy (if inputs are independent, otherwise it becomes challenging) because they take values in [0, 1], they sum to one and, a high value means the input is influential.

Estimation methods The estimation of Sobol' indices is a challenging issue because, except for rare cases, a high number of simulations is required. Indeed, to estimate a single Sobol' index with an uncertainty of 10%, it could be required to perform 10 4 model runs [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF]. For explanation purpose, it is convenient to divide the estimation techniques according to the sampling structure: structured sample design and given data methods. The structured sample design techniques can be applied if we have access to the model f but, in some applications we only dispose of a random i.i.d. n sample from (X, Y ). Figure 2.4 shows a structured sample and a non structured sample designs.

Among the structured sample methods, we can mention the Fast Amplitude Sensitivity Test (FAST) introduced by [START_REF] Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF] [see, e.g., Saltelli et al., 1999, and references therein], Sobol' pick-freeze schemes [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF], FAST mixed with Random Balance Designs [START_REF] Tarantola | Random balance designs for the estimation of first order global sensitivity indices[END_REF]. [START_REF] Tissot | A randomized Orthogonal Array-based procedure for the estimation of first-and second-order Sobol' indices[END_REF] proposed a pick-freeze procedure based on replicated sampling which reduce the estimation cost. In [START_REF] Prieur | Variance-Based Sensitivity Analysis: Theory and Estimation Algorithms[END_REF], the most efficient estimation methods are well presented. Moreover, depending on the estimation method, confidence intervals or significance test could be deduced from asymptotic results of the central limit theorem [see more details in Janon et al., 2014a]. However, bootstrap intervals can be also built, under some conditions, to construct error bars.

To overcome the drawback of structured sample methods (computational cost, access to the model), authors propose given data methods. Among the given data methods, we can cite the EASI spectral method of [START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF], [START_REF] Plischke | Global sensitivity measures from given data[END_REF] method relies on the notion of class-conditional densities, where a class is a sub-sample stemming from a suitable partition of the dataset, the nonparametric estimation methods of Da Veiga and Gamboa [2013] or [START_REF] Solís | Non-parametric estimation of the first-order Sobol indices with bootstrap bandwidth[END_REF], the fully Bayesian given data procedure proposed by [START_REF] Antoniano-Villalobos | Nonparametric estimation of probabilistic sensitivity measures[END_REF]. And more recently, the [START_REF] Gamboa | Global sensitivity analysis: a new generation of mighty estimators based on rank statistics[END_REF] estimators based on rank statistics. It is worth to mention that the EASI method gives good results in test cases even with a small bunch of observations (e.g., 200 samples). The two precursory given data methods, the EASI method and NSD by [START_REF] Plischke | Global sensitivity measures from given data[END_REF] method, are explained in more detail and compared to the structured sample method of [START_REF] Tissot | A randomized Orthogonal Array-based procedure for the estimation of first-and second-order Sobol' indices[END_REF] in Appendix A. In Chapter 4, we propose a given data method to estimate first-order Sobol' indices, based on a nonparametric Nadaraya-Watson [START_REF] Nadaraya | On estimating regression[END_REF][START_REF] Watson | Smooth regression analysis[END_REF] kernel regression with bias corrected [START_REF] Racine | Bias-Corrected Kernel Regression[END_REF] to quantify the importance of the inputs of an avalanche model.

Dependent inputs Independence assumption on input parameters is crucial although unrealistic in many applications for a comprehensive interpretation of Sobol' indices. In fact, if input parameters are dependent, the variance decomposition in Equation (2.6) used for the interpretation of Sobol' indices is not true anymore. Some authors have proposed strategies to estimate Sobol' indices if input parameters are dependent. [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF] propose to separate the dependent inputs in independent groups and to compute Sobol' indices for each group. In other words, they propose to organize the inputs in s independent groups X = (X u 1 , . . . , X us ) with u ⊂ {1, . . . , d} and then compute, the Sobol' index S u = Var(E(Y |X u ))/Var(Y ) for all ∈ {1, . . . , s}. This method has two main limitations: it assumes all variables are not correlated and, we cannot recognize the influence of single inputs. [START_REF] Da Veiga | Local polynomial estimation for sensitivity analysis on models with correlated inputs[END_REF] proposed to estimate Sobol' indices using local polynomials. [START_REF] Broto | Sensitivity indices for independent groups of variables[END_REF] studied the Sobol' indices for independent groups in the particular case of block-additive models.

Other works focused on alternative ANOVA based indices when inputs are dependent. [START_REF] Xu | Uncertainty and sensitivity analysis for models with correlated parameters[END_REF] [see also [START_REF] Zhang | A new framework of variance based global sensitivity analysis for models with correlated inputs[END_REF] proposed to decompose the contribution in two parts: the correlated and uncorrelated contributions, the method assumes that the linearity of both the model and dependence structure which limits its application. [START_REF] Li | Global sensitivity analysis for systems with independent and/or correlated inputs[END_REF] proposed three sensitivity indices to reflect the total, structural (linked to the system structure) and correlative contributions (linked to the correlated input distribution) extending the results of [START_REF] Xu | Uncertainty and sensitivity analysis for models with correlated parameters[END_REF] to more general models using the concept of High Dimensional Model representation (HDMR), Mara and Tarantola [2012] proposed a first step consisting into decorrelating the inputs using the Gram-Schmidt procedure and then, to perform the ANOVA-HDMR decomposition and compute the sensitivity measures. However, the decorrelation procedure depends on the variable order which makes intricate the index interpretation [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: numerical methods[END_REF]. [START_REF] Chastaing | Generalized Hoeffding-Sobol decomposition for dependent variables -application to sensitivity analysis[END_REF] define new sensitivity measures based on the Generalized Hoeffding decomposition proposed by [START_REF] Stone | The use of polynomial splines and their tensor products in multivariate function estimation[END_REF], these new indices are called the generalized Sobol' indices. However, the index construction requires specific assumptions on the joint distribution function of the inputs [START_REF] Chastaing | Generalized Hoeffding-Sobol decomposition for dependent variables -application to sensitivity analysis[END_REF]. Also, their interpretation is not trivial because even if they sum to one, they could be negative [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: numerical methods[END_REF]. [START_REF] Kucherenko | Estimation of global sensitivity indices for models with dependent variables[END_REF] proposed a copula-based approach to sample multivariate distributions but, it requires a priori knowledge of the conditional densities. Moreover, index interpretation may not be clear because if input correlation is high, first order index are higher than the total one. [START_REF] Thierry | Non-parametric methods for global sensitivity analysis of model output with dependent inputs[END_REF] established a link between [START_REF] Kucherenko | Estimation of global sensitivity indices for models with dependent variables[END_REF] and Mara and Tarantola [2012] works proposing the estimation of four indices, namely full and independent first-order and total indices. Recently, [START_REF] Hart | An approximation theoretic perspective of sobol' indices with dependent variables[END_REF] reformulated total Sobol' indices in terms of approximation errors rather than variance analysis. More precisely, they proposed to see total Sobol' indices as a measure of discrepancy between the model and a surrogate using less variables, i.e. the error to approximate f (X) by f (g(X -u ), X -u ) where g(X -u ) is an approximation of X u . However, all the afore-mentioned works have drawbacks and specially they do not provide an univocal way of partitioning the influence of input parameters on the output. For these reasons, Shapley effects seems more adapted in a dependent framework. Sensitivity measures when inputs are dependent is a still ongoing research subject which has been recently studied in the Master thesis of [START_REF] Wiederkehr | Global sensitivity analysis with dependent inputs[END_REF] and Ph. D. thesis of [START_REF] Broto | Sensitivity analysis with dependent variables: Estimation of the Shapley effects for unknown input distribution and linear Gaussian models[END_REF].

Shapley effects

In contrast to Sobol' indices, Shapley effects [START_REF] Owen | Sobol' Indices and Shapley Value[END_REF] are meaningful in the framework of dependent input parameters [START_REF] Owen | On Shapley value for measuring importance of dependent inputs[END_REF]. This measure is based on the Shapley value which is a cooperative game theory concept. Briefly speaking, Shapley value ensures a fair distribution of a gain among team players according to their individual contributions. As a sensitivity measure, [START_REF] Owen | Sobol' Indices and Shapley Value[END_REF] adapted the Shapley value [START_REF] Shapley | A Value for n-Person Games[END_REF] into the Shapley effects by considering model input parameters as players and the gain function as the output variance. The main advantage of such an approach is that it is possible to attribute a non negative sensitivity index to each parameter, and the sum of the indices is equal to one [Broto et al., 2020, Iooss andPrieur, 2019].

The definition of Shapley effects is based on Shapley value concept which is next explained. Given a set of d players in a coalitional game and a characteristic function val : S → R, val(∅) = 0, where S is the set of all subsets of {1, . . . , d}, the Shapley value (φ 1 , . . . , φ d ) is the only distribution of the total gains val({1, . . . , d}) to the players satisfying the desirable properties listed below: It is proved in [START_REF] Shapley | A Value for n-Person Games[END_REF] that according to the Shapley value, the amount that player i gets given a coalitional game (val, d) is:

φ i = 1 d u⊆-{i} d -1 |u| -1 (val(u ∪ {i}) -val(u)) ∀i ∈ {1, . . . , d}.
(2.10)

The Shapley effects are defined by considering the characteristic function of the game as:

val(u) = Var(E(Y |X u )) Var(Y ) , u ⊆ {1, . . . , d}
in Equation (2.10). Thus, the scalar Shapley effect of input i is defined as:

Sh i = 1 d Var(Y ) u⊆-{i} d -1 |u| -1 (Var(E(Y |X u∪i )) -Var(E(Y |X u ))).
(2.11)

Estimation methods [START_REF] Song | Shapley effects for global sensitivity analysis: Theory and computation[END_REF] mentioned that Shapley effects estimation is costly because one has to consider all possible subsets of players (inputs) which is 2 d . In their work, two estimators are proposed: a first estimator constructed by enumerating all possible players permutations, the permutations set is denoted S d and a second estimator, which is a modified and improved version of [START_REF] Castro | Polynomial calculation of the shapley value based on sampling[END_REF] algorithm, which selects randomly m permutations from S d . Note that m should be chosen as large as possible. [START_REF] Benoumechiara | Shapley effects for sensitivity analysis with dependent inputs: bootstrap and kriging-based algorithms[END_REF] proposed bootstrap confidence intervals for [START_REF] Song | Shapley effects for global sensitivity analysis: Theory and computation[END_REF] estimators. [START_REF] Song | Shapley effects for global sensitivity analysis: Theory and computation[END_REF] algorithms assume that it is possible to sample from the distribution of all subsets of the input parameters conditionally to the complementary set of input parameters which might be complicated in some applications. [START_REF] Plischke | Computing Shapley Effects for Sensitivity Analysis[END_REF] proposed an estimation algorithm based on [START_REF] Song | Shapley effects for global sensitivity analysis: Theory and computation[END_REF] algorithm and the Möbious inverse to reduce estimation computational cost. Their algorithm also requires the explicit knowledge of the conditional distributions. Finally, [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF] proposed given data estimators based on nearest-neighbor, i.e. estimators that can be computed from an i.i.d. sample, which is in general more convenient for real applications. In Chapter 5, we extend Broto et al. [2020] work to functional outputs (see Subsection 2.3.6) and we construct bootstrap confidence intervals for these estimators.

Example To illustrate the difficulty of ranking dependent inputs using Sobol' indices, we consider a linear model with three Gaussian inputs when two of them are highly correlated [example from Section 3.4 [START_REF] Iooss | Shapley effects for sensitivity analysis with correlated inputs: Comparisons with Sobol' indices, numerical estimation and applications[END_REF]. Let us consider the linear model:

Y = f (X) = β 0 + β T X with β = (β 1 , β 2 , β 3 ) T = (1, 1, 1) T , X = (X 1 , X 2 , X 3 ) T , X ∼ N (µ, Σ
) with µ = (0, 0, 0) T and Σ ∈ R 3×3 a positive-definite matrix. X 1 is independent from X 2 and X 3 , and X 2 and X 3 are highly correlated (correlation ρ = 0.9). We set σ 1 = σ 2 = 1 and σ 3 = 2. The covariance matrix is then:

Σ =    σ 2 1 0 0 0 σ 2 2 ρσ 2 σ 3 0 ρσ 2 σ 3 σ 2 3    =    1 0 0 0 1 1.8 0 1.8 4   
The variance of the output is

σ 2 = Var[f (X)] = 3 j=1 β 2 j σ 2 j + 2ρβ 2 β 3 σ 2 σ 3 .
There exists analytical expressions for Sobol' and total indices, and Shapley effects for this test case and they are the following:

• Shapley effects [analytical expression from [START_REF] Iooss | Shapley effects for sensitivity analysis with correlated inputs: Comparisons with Sobol' indices, numerical estimation and applications[END_REF]:

Sh 1 = β 2 1 σ 2 1 /σ 2 , Sh 2 = [β 2 2 σ 2 2 + ρβ 2 β 3 σ 2 σ 3 + ρ 2 2 (β 2 3 σ 2 3 -β 2 2 σ 2 2 )]/σ 2 , Sh 3 = [β 2 3 σ 2 3 + ρβ 2 β 3 σ 2 σ 3 + ρ 2 2 (β 2 2 σ 2 2 -β 2 3 σ 2 3 )]/σ 2 .
• Sobol' indices [analytical expression from [START_REF] Kucherenko | Estimation of global sensitivity indices for models with dependent variables[END_REF]:

0.0 0.2 0.4 0.6 0.8 1.0 Input Sensitivity measure • • • X 1 X 2 X 3 •

Shapley effects

First-order Sobol' Total Sobol'

Figure 2.5: Comparison between Shapley effects and Sobol' indices when inputs are correlated on a linear model with three Gaussian inputs. X 1 is independent from X 2 and X 3 , and X 2 and X 3 are highly correlated (correlation ρ = 0.9).

S 1 = 1 2 + σ 2 + 2ρσ , S 2 = (1 + ρσ) 2 2 + σ 2 + 2ρσ S 3 = (σ + ρ) 2 2 + σ 2 + 2ρσ
.

• Total Sobol' indices [analytical expression from [START_REF] Kucherenko | Estimation of global sensitivity indices for models with dependent variables[END_REF]:

S T 1 = 1 2 + σ 2 + 2ρσ , S T 2 = 1 -ρ 2 2 + σ 2 + 2ρσ S T 3 = σ 2 (1 -ρ 2 ) 2 + σ 2 + 2ρσ
.

Figure 2.5 shows the analytical sensitivity measures for the linear model example. This example illustrates how difficult it is to interpret Sobol' indices if inputs are correlated since first-order indices are higher than total indices (see the indices of X 2 and X 3 ). On the contrary to Sobol' indices, Shapley effects are easy to interpret.

Moment-free measures

In this section, we describe briefly other useful global sensitivity measures that could be applied when inputs are correlated. These measures are not used in the manuscript core chapters but they offer an interesting approach to GSA when inputs are dependent. Moment-free sensitivity measures were early proposed by [START_REF] Borgonovo | A new uncertainty importance measure[END_REF] [see also [START_REF] Borgonovo | Moment independent importance measures: New results and analytical test cases[END_REF]. These sensitivity measures are useful because they analyze the influence of input uncertainty on the whole output distribution, this is why they are called moment-free, and they can be defined also in the presence of correlated inputs [START_REF] Borgonovo | A new uncertainty importance measure[END_REF]. Indeed, [START_REF] Borgonovo | A new uncertainty importance measure[END_REF] affirms that the influence of an input on the model output is represented by the entire uncertainty distribution and not only its variance.

The Borgonovo's δ sensitivity measure for input X i is defined according to:

δ i = 1 2 E X i (s(X i )),
where s(

X i ) = |f Y (y)-f Y |X i (y)|dy, f Y is the probability density function of Y , f Y |X i is the conditional density function of Y given X i .
Borgonovo' δ quantifies the shift area between f Y and f Y |X i . Borgonovo's δ has some nice properties: (i) δ i lies between 0 and 1 for all i ∈ {1, . . . , p}, and (ii) it is equal to 0, if Y is independent of X i . Borgonovo's δ can also be applied to group of inputs [see Equation 21 of [START_REF] Borgonovo | A new uncertainty importance measure[END_REF].

Da [START_REF] Da | Global sensitivity analysis with dependence measures[END_REF] extended the work of [START_REF] Borgonovo | A new uncertainty importance measure[END_REF] by proposing a general GSA framework based on the concept of dissimilarity measures (e.g. the Csiszàr f-divergence, integral probability metrics) and dependence measures (e.g., Hilbert-Schmidt independence criterion HSIC). More precisely, the impact of X i on Y could be given by [see also [START_REF] Baucells | Invariant probabilistic sensitivity analysis[END_REF]:

D i = E X i (d(P Y , P Y |X i )),
where P Y and P Y |X i are the probability measures of Y and Y |X i , respectively and, d is a dissimilarity measure between two probability measures. Borgonovo's δ is a particular case where the dissimilarity measure is the Csiszàr f-divergence. The first-order Sobol' index is also a particular case if

d(P Y , P Y |X i ) = (E(Y ) -E(Y |X i )) 2 .
It is worth to mention that HSIC measure developed for GSA by [START_REF] Da | Global sensitivity analysis with dependence measures[END_REF] [see also [START_REF] Meynaoui | New statistical methodology for second level global sensitivity analysis[END_REF] can be applied to perform target sensitivity analysis (TSA) which consists into measuring the input importance on a restricted or critical domain of the model output [see more details in [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF]. TSA has a wild spectrum of applications e.g., in reliability and risk management where the critical output domain is sometimes associated to a low probability [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF]. Also, HSIC measure can be applied to optimization problems with variable selection. The variable selection is performed through GSA. Indeed, in optimization problems with large number of parameters, identifying the influential inputs might be crucial to reduce the problem dimension and therefore, the computational cost of such complex applications [see more detail in [START_REF] Spagnol | Global sensitivity analysis for optimization with variable selection[END_REF].

Other measures

There are other global sensitivity measures, to mention some of them we can cite the derivative based global sensitivity measures (DGSM) proposed by [START_REF] Sobol | Derivative based global sensitivity measures and their link with global sensitivity indices[END_REF]. DGSM extend local sensitivity measures based on derivatives by computing integrals of the square derivatives on the whole input domain i.e. E ∂f ∂X i 2 (see also [START_REF] Lamboni | Derivative-based global sensitivity measures: General links with Sobol' indices and numerical tests[END_REF]). DGSM requires less model evaluations in comparison to other GSA measures [START_REF] Lamboni | Derivative-based global sensitivity measures: General links with Sobol' indices and numerical tests[END_REF]. However to apply DGSM, regularity assumptions might be required. Also, they are not recommended for highly nonlinear functions because the ranking obtained may suggest false conclusions [START_REF] Sobol | Derivative based global sensitivity measures and their link with global sensitivity indices[END_REF]. Global sensitivity metrics from active subspaces, called activity scores, were recently proposed by [START_REF] Constantine | Global sensitivity metrics from active subspaces[END_REF]. Active subspaces are emerging tools based on spectral decomposition for dimension reduction [START_REF] Constantine | Global sensitivity metrics from active subspaces[END_REF]. Finally, one could fit a linear model to obtain some information about the model, and then compute Pearson correlation or partial correlation coefficients, these methods are called the sampling-based GSA methods [START_REF] Helton | Survey of samplingbased methods for uncertainty and sensitivity analysis[END_REF]. However, these kind of methods require linear and/or monotony hypotheses which limit their application [see Iooss and Lemaître, 2015, Section 3.1.]. Nevertheless, scatterplots will always give some preliminary information about the input-output dependence structure (e.g., linear, non linear) thus, they should not be neglected if the number of parameters is moderate4 . This exposition was not exhaustive but it provides good insights of the current works developed in the GSA community. Other surveys of GSA can be found in e.g., [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF] and Borgonovo et al. [2016].

A little bit about metamodels

As already mentioned, the estimation of GSA measures, require large number of model simulations. If a model is computationally expensive, the construction of a metamodel or surrogate could be necessary to compute GSA measures. A metamodel is a mathematical construction to approximate f with a low computational cost. There exists a lot of metamodels (e.g., linear regression model, polynomials, splines). In particular, some metamodels provide Sobol' indices analytical expressions which means, they could be directly deduced without any additional cost. For example, Sobol' indices can be deduced from Gaussian process metamodels [START_REF] Marrel | Calculations of sobol indices for the gaussian process metamodel[END_REF] [see also [START_REF] Durrande | Anova kernels and rkhs of zero mean functions for model-based sensitivity analysis[END_REF], for a special case of ANOVA kernels], polynomial of chaos [START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF].

Multivariate and/or functional model outputs

Until here, we supposed Y scalar but some avalanche model outputs are functional. For example, Figure 2.6 shows the functional high density region (FHDR) bloxplots [see FHDR in [START_REF] Hyndman | Rainbow plots, bagplots, and boxplots for functional data[END_REF] of the functional outputs (panel a flow depth and panel b velocity) and scalar output (panel c runout distance) of n = 5000 random samples from the Saint-Venant [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF] avalanche model. If the output is multivariate f (X) = (Y 1 , . . . , Y p ) or functional Y = f s (X) with s ∈ τ ⊂ R where τ is the temporal or spatial domain of the phenomenon, it is possible to apply sensitivity analysis to each component of the multivariate output but this could lead to redundancies in the results, particularly in the setting of discretized functional outputs. [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] and Gamboa et al. [2013] extended Sobol' indices to multivariate or functional outputs, see also [START_REF] Alexanderian | Variance-based sensitivity analysis for time-dependent processes[END_REF] whose work extended Sobol' indices to time-dependent outputs.

The aggregated Sobol' indices are defined as:

GS i = p j=1 Var(Y j )S j i p j=1 Var(Y j )
, where S j i is the first-order Sobol' index of the output Y j with respect to the input X i given in Equation (2.7). x start

Dimension reduction techniques for functional outputs

x runout

• • • • • • • • • • • •
topo start and stop for each avalanche sample Shapley effects can be naturally extended to multivariate output models by following the ideas presented by [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF], Gamboa et al. [2013]. Thus, aggregated Shapley effect of an input i is then defined as:

GSh i = p j=1 Var(Y j )Sh j i p j=1 Var(Y j ) ,
(2.12)

where Sh j i is the scalar Shapley effect of input X i on output Y j given in Equation (2.11). In Chapter 5, we introduce more precisely the extension of Shapley effects to multivariate or functional outputs. Note that aggregated HSIC indices proposed by [START_REF] Da | Global sensitivity analysis with dependence measures[END_REF] could be also used for multivariate outputs.

If the output dimension is high (as it is the case, e.g., when considering the discretization of a functional output), a dimension reduction should be applied as a preliminary step to estimate efficiently aggregated Sobol' or Shapley effects. Some examples of dimension reduction techniques are presented in Section 2.4.

Dimension reduction techniques for functional outputs

Campbell et al. [2006] suggested that if model output is functional, sensitivity analysis should be carried out in two steps: first, to expand the functional output in terms of a set of basis functions and then, to perform sensitivity analysis on the coefficients of the expansion. In general, dimension reduction consists into projecting the output matrix Y ∈ R n×Ns , composed by n discretized output curves in a grid of size N s , in a basis of lower dimension Φ q×Ns with q ≤ N s :

Y -Ȳ ≈ HΦ T ,
where Ȳ = n -1 11 T Y with 1 a n-vector of ones, the columns of Φ codes the basis functions and, H the coefficients of the expansion. Thus, each curve ∈ {1, . . . , n} can be approximated as:

y (s) -ȳ(s) ≈ q k=1 h k φ k (s) for all s ∈ τ
Principal component analysis (PCA) is one of the most common techniques employed for dimension reduction. PCA consists in the spectral decomposition, i.e. eigenvalues and vector decomposition, of the variance-covariance matrix of the output Y. Note that if the output is functional, functional PCA (fPCA) might be preferred (see [START_REF] Ramsay | Functional Data Analysis[END_REF] and [START_REF] Yao | Functional data analysis for sparse longitudinal data[END_REF]). Legendre polynomials, Haar functions or wavelet bases are other bases that could be used for dimension reduction [START_REF] Campbell | Sensitivity analysis when model outputs are functions[END_REF]. We choose to use PCA in Chapter 4 and fPCA in Chapter 5, whose use in GSA is common (e.g., [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF][START_REF] Lamboni | Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models[END_REF], [START_REF] Alexanderian | Variance-based sensitivity analysis for time-dependent processes[END_REF]).

Part I Model calibration using Bayesian inference

3

Bayesian calibration of an avalanche model from autocorrelated measurements along the flow: application to velocities extracted from photogrammetric images Abstract: Physically-based avalanche propagation models must still be locally calibrated to provide robust predictions, e.g. in long-term forecasting and subsequent risk assessment. Friction parameters cannot be measured directly and need to be estimated from observations. Rich and diverse data are now increasingly available from test-sites, but for measurements made along flow propagation, potential autocorrelation should be explicitly accounted for. To this aim, this work proposes a comprehensive Bayesian calibration and statistical model selection framework. As a proof of concept, the framework was applied to an avalanche sliding block model with the standard Voellmy friction law and high rate photogrammetric images. An avalanche released at the Lautaret test-site and a synthetic data set based on the avalanche are used to test the approach and to illustrate its benefits. Results demonstrate (i) the efficiency of the proposed calibration scheme, and (ii) that including autocorrelation in the statistical modelling definitely improves the accuracy of both parameter estimation and velocity predictions. Our approach could be extended without loss of generality to the calibration of any avalanche dynamics model from any type of measurement stemming from the same avalanche flow.

Introduction

The complexity of snow avalanche physics is related to the variability and changing nature of snow [START_REF] Schweizer | Snow avalanche formation[END_REF][START_REF] Ancey | Dynamique des avalanches[END_REF][START_REF] Castebrunet | Snow and weather climatic control on snow avalanche occurrence fluctuations over 50 yr in the French Alps[END_REF][START_REF] Steinkogler | Influence of snow cover properties on avalanche dynamics[END_REF]. Evidence obtained at full-scale experimental slopes [START_REF] Sovilla | Impact pressures and flow regimes in dense snow avalanches observed at the vallée de la sionne test site[END_REF][START_REF] Vriend | High-resolution radar measurements of snow avalanches[END_REF], Prokop et al., 2015[START_REF] Faug | Looking beyond the powder/dense flow avalanche dichotomy[END_REF] shows a myriad of avalanche propagation and stopping regimes [START_REF] Köhler | The dynamics of surges in the 3 february 2015 avalanches in vallée de la sionne[END_REF], and numerical propagation models can reproduce these observations with increasing realism [START_REF] Bartelt | Configurational energy and the formation of mixed flowing/powder snow and ice avalanches[END_REF][START_REF] Gaume | Dynamic anticrack propagation in snow[END_REF]. However, knowledge concerning the mechanical behaviour of snow during motion and associated processes (granulation, erosion/deposition, etc.) remains incomplete [START_REF] Steinkogler | Granulation of snow: From tumbler experiments to discrete element simulations[END_REF][START_REF] Ancey | Snow avalanches[END_REF][START_REF] Truong | Refining the processing of paired time series data to improve velocity estimation in snow flows[END_REF]. From a macroscopic point of view, experimental approaches [START_REF] Casassa | Measurements of friction coefficients of snow blocks[END_REF][START_REF] Rognon | Rheology of dense snow flows: Inferences from steady state chuteflow experiments[END_REF][START_REF] Kern | Measured shear rates in large dry and wet snow avalanches[END_REF] and the proposal of [START_REF] Voellmy | On the destructive force of avalanches[END_REF] suggest rheological behaviours which remain ad-hoc. This renders on-site calibration on the basis of local data unavoidable [START_REF] Ancey | Estimating bulk rheological properties of flowing snow avalanches from field data[END_REF][START_REF] Salm | A short and personal history of snow avalanche dynamics[END_REF], Eckert et al., 2012] to, e.g. predict high-return-period avalanches in land use planning and assess the related risk [START_REF] Keylock | Avalanche risk mapping by simulation[END_REF], Meunier et al., 2004, Favier et al., 2014b,a]. This is all the more true given that studies have shown that avalanche propagation models are highly sensitive to their friction parameter values [see e.g., [START_REF] Borstad | Sensitivity analyses in snow avalanche dynamics modeling and implications when modeling extreme events[END_REF]McClung, 2009, Fischer, 2013].

After deterministic inversion methods had shown their limits [START_REF] Dent | Modeling of snow flow[END_REF][START_REF] Dent | Density, velocity and friction measurements in a dry-snow avalanche[END_REF][START_REF] Ancey | Inverse problem in avalanche dynamics models[END_REF], and following progress made in many fields where accurate numerical model calibration is now recognised as a crucial issue [e.g., [START_REF] Oakley | Bayesian inference for the uncertainty distribution of computer model outputs[END_REF]O'Hagan, 2002, Carmassi et al., 2018], the Bayesian framework has become an appealing avenue in snow science over the past years, especially in the frequent case of small data samples [e.g., [START_REF] Ancey | Monte Carlo Calibration of Avalanches Described as Coulomb Fluid Flows[END_REF][START_REF] Straub | A Bayesian probabilistic framework for avalanche modelling based on observations[END_REF][START_REF] Eckert | Revisiting statistical-topographical methods for avalanche predetermination: Bayesian modelling for runout distance predictive distribution[END_REF][START_REF] Eckert | Bayesian stochastic modelling for avalanche predetermination: from a general system framework to return period computations[END_REF], 2009[START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF][START_REF] Schläppy | Validation of extreme snow avalanches and related return periods derived from a statistical-dynamical model using tree-ring techniques[END_REF]. Specifically, [START_REF] Gauer | Optimization and probabilistic calibration of avalanche block models[END_REF] used a Bayesian framework to calibrate the friction parameters of three avalanche sliding block models and [START_REF] Fischer | Evaluation of probabilistic snow avalanche simulation ensembles with Doppler radar observations[END_REF] proposed a method to evaluate simulations compared to Doppler radar observations. Ultimately, [START_REF] Naaim | Dense avalanche friction coefficients: influence of physical properties of snow[END_REF] could establish empirical links between friction parameters and physical properties of snow.

However, most of these existing approaches remain limited to rather coarse field data (e.g., samples of runout distances supplemented by input conditions), and when more comprehensive data sets have been considered, improper likelihood formulations have been used more often than not [START_REF] Fischer | Multivariate parameter optimization for computational snow avalanche simulation[END_REF]. For instance, little attention has been given so far to the specific difficulty induced by potential autocorrelation between the data used for calibration. Whereas this is not a matter for example, a sample of runout corresponding to distinct avalanche events, the assumption of independent observations is much more questionable in the current context of the increasingly diverse and rich measurements made on test-sites within the same avalanche. Many environmental applications have indeed demonstrated that neglecting potential autocorrelation between different measurements used in a calibration scheme can lead to biases in parameter esti-mation and/or to lower predictive performances [see e.g., [START_REF] Kuczera | Improved parameter inference in catchment models: 1. Evaluating parameter uncertainty[END_REF][START_REF] Mcinerney | Improving probabilistic prediction of daily streamflow by identifying Pareto optimal approaches for modeling heteroscedastic residual errors[END_REF][START_REF] Schaefli | Bayesian spectral likelihood for hydrological parameter inference[END_REF][START_REF] Sun | Effect of heteroscedasticity treatment in residual error models on model calibration and prediction uncertainty estimation[END_REF]. If this happens, both the physical interpretation of parameter estimates and the operational use of model predictions can be questioned.

On this basis, in this paper, we propose a Bayesian approach to calibrate the friction parameters of an avalanche propagation model using data of high temporal resolution (of the order of 1 s). Inspired by studies done for the calibration of hydrological models [e.g., [START_REF] Kuczera | Improved parameter inference in catchment models: 1. Evaluating parameter uncertainty[END_REF]Parent, 1998, Evin et al., 2014], the main novelty of our work is to explicitly account for potential autocorrelation between measurements made along the avalanche flow within the calibration framework. In what follows, we demonstrate that doing so results in a model better supported by the data that improves the accuracy of friction parameter estimates and of velocity predictions. Application is made on a well documented avalanche event from the Lautaret full-scale test-site [START_REF] Thibert | The full-scale avalanche test-site at Lautaret Pass (French Alps)[END_REF]. The used velocity data were obtained from high rate positioning from photogrammetric images. Such data were already applied by various authors to avalanche simulations performed with reference friction parameter values [e.g., [START_REF] Dreier | Comparison of simulated powder snow avalanches with photogrammetric measurements[END_REF][START_REF] Gauer | Comparison of avalanche front velocity measurements and implications for avalanche models[END_REF][START_REF] Turnbull | Mass and Momentum Balance Model of a Mixed Flowing/Powder Snow Avalanche[END_REF], but without including them within an explicit calibration scheme. That taking into account potential autocorrelation between measurements within the calibration is unavoidable to get unbiased estimates is further demonstrated with synthetic data analogue to the case-study.

In this work, as a proof of concept, we use the sliding block avalanche model, also known as 1-D Voellmy model instead of state of the art depth-averaged models [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF], Christen et al., 2010a[START_REF] Bartelt | Modeling mass-dependent flow regime transitions to predict the stopping and depositional behavior of snow avalanches[END_REF]. We are aware of its limitations, notably that it makes it impossible to depict, e.g., flow depth variations in space and time. Also, some authors have shown that it may underestimate avalanche velocities [see e.g., Ancey andMeunier, 2004, Gauer, 2014]. However, for hazard mitigation, simple models with few parameters remain useful [START_REF] Salm | A short and personal history of snow avalanche dynamics[END_REF] and have the advantage to allow fast computation in comparison to more complex ones. A very simple avalanche propagation model is therefore a good choice for developing a calibration approach which could be in the future be applied to any other, more advanced, avalanche propagation model as soon as autocorrelation in measurements series is suspected.

Avalanche model calibration principle

Sliding block propagation model

Our model considers the avalanche as a rigid body sliding over a bidimensional curvilinear profile starting from the top of the path. The mass m and body shape variations of the avalanche are neglected. Under these assumptions, the motion equation of the avalanche mass center is:

du dt = g sin θ - F m , (3.1)
where u = || u||, du dt is the acceleration, g is the gravity constant, θ is the local slope angle and F = || F || is the frictional force. In this study, we consider the classical Voellmy friction law [START_REF] Voellmy | On the destructive force of avalanches[END_REF]. This means that the friction force is:

F = µmg cos θ + mg ξh u 2 , (3.2)
depending on two friction parameters Θ = {µ, ξ}. Often it is assumed that µ evolves with the physical properties of snow whereas ξ may correspond to the geometry of the avalanche path and to terrain roughness [START_REF] Ancey | Inverse problem in avalanche dynamics models[END_REF]. However, whereas this interpretation is sound or not is not our debate here. We do not make any further assumption regarding the linkages between (µ, ξ) and snow and topographical variables and simply search for the best couple on the basis of the data. The propagation model also depends on three forcing variables x = {T, h, x start } where T defines the topography of the terrain, h the mean flow depth of the avalanche and x start the release abscissa (the protected runout length) of the avalanche mass center [Ancey andMeunier, 2004, Eckert et al., 2007].

Statistical model formulation

Let us denote f our avalanche propagation model. The model f predicts the avalanche speed (m s -1 ) and position along the slope (m) at a time t (s). The model depends on parameters Θ ∈ R p and forcing variables x ∈ R d . The observed velocity, collected on the field, at time t is noted by v t and we denote v obs = {v 1 , . . . v n } the set of observations where n denotes the number of observations. The aim of model calibration is i) to find the optimal combination of parameters Θ that minimizes the discrepancy between the observations v obs and the model simulation f (Θ, x) and ii) to rigorously quantify the associated uncertainty. To this end, we use the generic statistical model: (3.3) where f t (Θ, x) denotes the simulation of the avalanche velocity at time t and t is the model error. With this classical additive formulation, propagation model errors and observation errors are modelled altogether in the residuals t .

M : v t = f t (x, Θ) + t , ∀t ∈ {1, . . . n},
In nearly all existing avalanche model calibration approaches, model errors are implicitly or explicitly assumed as independent and identically normally distributed (i.i.d.). In other words, t i.i.d.

∼ N (0, σ 2 ), ∀t ∈ {1, . . . n}, where N denotes the Normal distribution and σ 2 is the common variance of the errors. However, this may be a too strong assumption for different measurements made along the same avalanche. For instance, the errors t are likely to present a non-negligible correlation between two consecutive time steps (in other words, t and t-1 may be correlated). To include the errors' autocorrelation in the calibration we propose to model them as an autoregressive (AR) process. Specifically, we consider an AR model of order 1 only (AR1). AR models of higher order could be used but this would imply the estimation of additional parameters, a non-trivial task in our case of a limited velocity sample. In addition, results obtained for the application indeed suggest that an AR(1) is sufficient to accurately account for the data variability (see Section 3.4). Thus, the errors t can be expressed as: (3.4) where φ ∈ R and η t are the coefficient and innovations of the AR(1) process, respectively. If |φ| < 1, ( t ) t∈N is defined as the unique stationary solution of equation (3.4). In

t = φ t-1 + η t , η t i.i.d. ∼ N (0, σ 2 ) ∀t ∈ {1, . . . n},

Model Errors

Parameters M 0 Independent and identically distributed N (0, σ 2 ) Θ, σ 2 M 1 AR(1) Θ, φ, σ 2 Table 3.1: Considered statistical models and corresponding parameters.

our study, |φ| < 1 (see Section 3.4) that guarantees the stationarity of ( t ) t∈Z . Note also that with this model the innovations η t are normally independent distributed but not the errors t .

Hereafter, the model with the assumption of normally distributed and independent errors is denoted M 0 and the model with AR(1) errors is denoted M 1 . Depending on the M i model with i ∈ {0, 1}, the whole set of parameters to be estimated is different. Table 3.1 summarizes the parametrization of the competing statistical models considered. Hence, model M 1 with φ = 0 corresponds to M 0 .

Bayesian framework

The probability of the data can be maximised with respect to the model parameters [START_REF] Fisher | On the mathematical foundations of theoretical statistics[END_REF]. Bayesian statistics is an useful framework to estimate model parameters from scarce data. Within this approach, the quantification of uncertainty in parameter estimation is straightforward. In fact, the advantage of the Bayesian approach is that the uncertainty on parameters is assessed through credibility intervals by contrast to traditional methods (confidence intervals) [Bayes, 1763, Bernardo and[START_REF] José | Bayesian theory[END_REF]. Hence, Bayesian statistics is now widely accepted as a reasonable option in environmental sciences [START_REF] Berliner | Physical-statistical modeling in geophysics[END_REF][START_REF] Clark | Why environmental scientists are becoming Bayesians[END_REF] and we use this framework in what follows.

For simplicity, let us denote γ, the set of parameters related to the errors, it means γ = {σ 2 } for M 0 and γ = {φ, σ 2 } for M 1 , and x = {T, h, x start } the forcing variables. Under a Bayesian framework, the joint posterior distribution of the parameters is the following one:

π(Θ, γ|v obs , x) ∝ L(v obs |Θ, γ, x)π(Θ, γ), (3.5)
where L(v obs |Θ, γ, x) is the likelihood of the observations and π(Θ, γ) is the joint prior distribution. Thus, the likelihood L(v obs |Θ, γ, x) is required. Under the M 0 model, the observations follow a Gaussian distribution and the likelihood L(v obs |Θ, γ, x, M 0 ) writes as:

1 (2πσ 2 ) n/2 exp - 1 2σ 2 n t=1 (v t -f t (Θ, x)) 2 .
(3.6)

Under the M 1 model, the likelihood L(v obs |Θ, γ, x, M 1 ) writes as (for more detail, see Appendix A):

1 -φ 2 2πσ 2 exp - 1 -φ 2 2σ 2 (v 1 -f 1 (Θ, x)) 2 × (3.7) 1 (2πσ 2 ) (n-1)/2 exp - 1 2σ 2 n t=2 η 2 t .

Avalanche model calibration principle

Note that, in this study, we do not calibrate the quantities x start , h, T because they were inferred from the data to put the effort on the friction law calibration. The sensitivity to these quantities in avalanche models has been studied in the works of e.g., [START_REF] Barbolini | Estimate of uncertainties in avalanche hazard mapping[END_REF], [START_REF] Borstad | Sensitivity analyses in snow avalanche dynamics modeling and implications when modeling extreme events[END_REF], [START_REF] Buhler | Sensitivity Analysis of the RAMMS Avalanche Dynamics Model in a Canadian Transitional Snow Climate[END_REF]. The authors found that changes in avalanche volume have a larger effect on both runout distance and avalanche velocity and that the friction coefficient µ (in a Coulomb model) is of high importance. However, this question should be analyzed deeper in a formal statistical framework and it is out of the scope of this work.

Metropolis-Hastings algorithm

The main practical difficulty in Bayesian inference is how to compute the normalizing constant in Bayes theorem. We overcome it by implementing a sequential Metropolis-Hastings algorithm, hereafter denoted MH [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo Sampling Methods Using Markov Chains and Their Applications[END_REF]. The MH algorithm proposes a generic way to construct a stationary and ergodic Markov chain that converges, under mild conditions, to the posterior distribution π(Θ, γ|v obs , x). The Markov chain returned by the algorithm can be considered as a sample from π(Θ, γ|v obs , x).

The following description of the MH algorithm is obtained from [START_REF] Robert | The Metropolis-Hastings algorithm[END_REF]. Let us denote, for simplicity, ψ = {Θ, γ} the set of the error and model parameters. For the application of the MH algorithm, it is needed an initial value ψ (0) and a proposal distribution q. Each iteration k of the algorithm consists in:

1. Generating ψ ∼ q(.|ψ (k-1) ).

2. Calculating u ∼ U(0, 1).

Taking

ψ (k) = ψ if u ≤ α, ψ (k-1) otherwise (3.8) with α = min L(v obs |ψ , x)π(ψ ) L(v obs |ψ (k-1) , x)π(ψ (k-1) ) q(ψ (k-1) |ψ ) q(ψ |ψ (k-1) ) , 1 (3.9)
where we recall that L(v obs |ψ, x) and π(ψ) = π(Θ, γ) stand for the likelihood of velocity observations and the joint prior distribution, respectively. As mentioned in [START_REF] Robert | The Metropolis-Hastings algorithm[END_REF], the performance of the algorithm depends on the choice of q. For example, a standard choice for the proposal distribution q is a multinormal distribution centered in ψ (k-1) and with a given covariance Σ q , which defines a random walk. For our application, Σ q was tuned according to the optimal acceptance rates that grant fast convergence of the MH algorithm [START_REF] Gelman | Inference from Iterative Simulation Using Multiple Sequences[END_REF].

Model selection using Bayes Factor

To compare the accuracy of the two competing models M 0 and M 1 , we use the Bayes factor. The Bayes factor is a criterion of the evidence provided by the data (in our case v obs ) to reject a model (M 0 ) compared to another one (M 1 ). The Bayes factor is defined as the ratio of marginal probabilities [see more details in [START_REF] Kass | Bayes Factors[END_REF]:

B 10 = p(v obs |M 1 , x) p(v obs |M 0 , x) , (3.10) where p(v obs |M i , x) = L(v obs |ψ i , M i , x)π(ψ i |M i , x)dψ i , i ∈ {0, 1}. (3.11)
The conditioning by M i highlights that likelihood, prior and posterior distributions depend on the considered model. The Bayes factor can be estimated by applying importance sampling and using the MH sample drawn from the posterior density π(ψ i |v obs , x, M i ) (see more details in Appendix B Section 3.7). An interpretation of the numerical value obtained to determine if there is evidence provided by the data to reject the model M 0 compared to M 1 according to the value of log 10 B 10 is proposed in the work of [START_REF] Kass | Bayes Factors[END_REF], where log 10 denotes the common logarithm. Notably, a log 10 B 10 value between 1 and 2 indicates strong evidence in favour of rejecting M 0 compared to M 1 , and log 10 B 10 > 2 a decisive evidence in favour of M 1 .

Application data

Avalanche data from the Lautaret test-site

The case study is an avalanche released at the Lautaret full-scale test-site [START_REF] Thibert | The full-scale avalanche test-site at Lautaret Pass (French Alps)[END_REF]. This test-site, located in the French Alps, holds a succession of avalanche paths.

Here, the path referred as "path number 2" was used to artificially trigger an avalanche on 13 February 2013 (Figure 3.1). This path is 450 m long, dropping from 2400 to 2100 m a.s.l. (Figure 3.2). Its upper part part, where the acceleration of the flow occurs, is steep with an average inclination of 37 • and a maximal slope around 45 • in the starting area. This part of the path is steep-sided and around 10 m-wide, so that the flow is channelized. The lower part of the path is mainly the run-out zone, a large and naturally open slope. At the transition between the two parts, a road (Col du Galibier road, open in summer only) crosses the avalanche path, which generates a local slope rupture in the avalanche track.

The avalanche released on February 13 2013 was composed of a dense part with a limited saltation layer on top. Properties of the snow involved in the flow were characterized with a density, temperature and hardness profile of the released snow layer. Snow grain types and dimensions were also characterized. The avalanche was released at 11h58 AM when air temperature was about -10 • Celsius. A 0.25 m thick layer of fragmented and decomposing snow particles was released. The mean density was 250 kg.m -3 , ranging between 270 and 225 kg.m -3 from the bottom to the surface of the released snow cover. In the released layer, particles diameter is less than 0.5 mm. Snow temperature was between -4.7 and -5 • Celcius. Hardness was "fist" (hand index) and measured as 20 N in Ram Resistance Equivalents [START_REF] Crla Fierz | The International Classification for Seasonal Snow on the Ground[END_REF].

Avalanche front positions were determined using a high rate photogrammetric system that was specifically developed [START_REF] Soruco | Measurement of avalanche front velocity from high-speed terrestrial digital photogrammetry[END_REF], Thibert et al., 2015]. We used a lowcost non-metric imagery system (numerical reflex Nikon D2Xs cameras in DX format, CMOS sensor with Nikon 85 mm f/1.4 AF fixed focal lenses). An advanced ad hoc tuning was performed to account for the radial distortion of the lenses, the decentration 3.3. Application data of the principal point (principal point shift) and the exact focal length of the lenses required for the correct scaling of the images [START_REF] Faig | A non-metric stereo system and its calibration[END_REF]. The resulting error in positioning avalanche fronts was estimated to be less than 25 cm after image orientation on ground control points and a comparison of direct photogrammetric measurements to laser scanning on a test area. Synchronization between the two cameras was achieved within a precision of 6.10 -6 s, therefore the error associated with the time sequence is also negligible. Eventually, terrestrial laser scanning was used to retrieve terrain elevation before and after avalanche triggering [START_REF] Prokop | Assessing the applicability of terrestrial laser scanning for spatial snow depth measurements[END_REF], Prokop et al., 2015], so as to quantify the snow mass transfer by the avalanche (3.1). In Appendix D, Figure 3.10 shows some examples of images used in the photogrammetric process.

The sliding block model is a one dimensional model representing the successive positions of the center of mass of the avalanche, thus we used the following procedure to determine the velocity vector v obs . At each time t, we computed the center of mass position of the avalanche using a set of front points D t ∈ R 3 from the delineated front. D t excludes the lateral front points because they are less active in the avalanche flow [START_REF] Pulfer | Retrieving avalanche basal friction law from high rate positioning of avalanches[END_REF]:

c t = pt∈Dt p t n t , (3.12)
where n t is the number of elements in D t . Note that, from the data available, it is not possible to calculate with certainty the position of the mass center of each avalanche front. In this study, we estimated each of them using equation (3.12), arguably a rough estimation but sufficient for obtaining an approximation of the avalanche velocity consistent with the modelling framework we are using. Figure 3.1 shows the front of the avalanche and the center of mass at each five seconds. Then, the velocity at time t is estimated as the mean velocity of the center of mass between two successive images:

v t = c t -c t-1 ∆t . (3.13)
In our application, ∆t = 1s and the data set is composed of 21 observations. We consider the first 21 seconds before the avalanche splits in sub-avalanches to ensure the validity of the application of a sliding block model. The mean flow depth h value was calculated as the mean of the difference between the snow depth of the digital elevation model taken after the release of the avalanche and the snow depth registered at the mass center locations during the avalanche motion. The value calculated was h = 2.19m with a standard deviation of 0.78m. This estimation is rough but it is a reasonable value given the characteristics of the avalanche studied.

The topography T was constructed as follows: the site under study has a simple geometry. It is a rather straight avalanche path starting with a cornice and of limited and rather constant width. The chosen 2-D topography is the main flow path starting from the cornice in the middle of the path (see more details in [START_REF] Thibert | The full-scale avalanche test-site at Lautaret Pass (French Alps)[END_REF]). From a fine-scale Digital Elevation Model, it was extracted as a grid of horizontal resolution of 1.4m. This is largely small enough to make the impact of the numerical approximation on the computed velocities negligible. Indeed, [START_REF] Bühler | Sensitivity of snow avalanche simulations to digital elevation model quality and resolution[END_REF] showed that a spatial resolution of 25m is sufficient for avalanche modelling.

Figure 3.2 shows the topography T and the points where the fronts were recorded.

Parameter MAP estimate

µ MAP 0.32 ξ MAP [m.s -2 ] 495.73 σ 2 MAP [m 2 .s -2 ] 8.07 φ MAP 0.65
Table 3.2: Parameter values used for the generation of the synthetic data set. They correspond to the maximum (MAP) obtained with model M 1 for the avalanche of the 13 February 2013.

Synthetic data generation

In order to further evaluate the accuracy of parameter estimation with the two models, we used synthetic data. This standard approach in statistical developments allows validation because "truth" (i.e. the parameter values used for the synthetic data generation) is known. Any parameter values could have been used at this stage. To resemble the most possible to the measured avalanche velocities, we chose the combination of parameter estimates corresponding to the measured case study. Also, to highlight that taking into account autocorrelation between measurements within the calibration is crucial to get unbiased estimates when such autocorrelation actually exists, our synthetic data were generated with model M 1 . Specifically, synthetic velocity time series were generated from the maximum a posteriori (MAP) resulting from the application of our calibration approach to the avalanche described above with model M 1 . The maximum a posteriori is the mode of the posterior distribution, namely the most plausible value given the data (Table 3.2). In detail, after the MAP estimators of the µ, ξ, σ 2 and φ parameters of M 1 model were determined for the avalanche, we proceed as follows:

1. An avalanche model simulation was conducted using the µ MAP and ξ MAP parameters.

The result of this simulation is a velocity time series denoted by f (Θ MAP , x).

2. An AR(1) error of parameters σ 2 MAP and φ MAP was simulated. 3. The AR(1) error simulated was added to f (Θ MAP , x). If negative velocity was obtained, the AR(1) error was resampled.

We generated 100 samples of synthetic data following this procedure and, for each of these, the 21 velocity values which correspond to the same location of the measurement were kept (to stay with the same data size as for the avalanche). Figure 3.3 shows some examples of the velocity samples generated versus the avalanche data.

Application, results and discussion

The avalanche Prior distributions

Our Bayesian framework was used to calibrate the parameters of the M 0 and M 1 models (see Table 3.1). Marginal prior distributions used are described in Table 3.3. Parameters were assumed marginally independent a priori. This assumption is not strong because the q q q q q q q q q q q q q q q q q q q q q q 0 10 20 µ ∼ Γ(9, 0.03) dependence of the parameters is reflected in their joint posterior distribution [START_REF] Walter R Gilks | Markov chain Monte Carlo in practice[END_REF], Eckert et al., 2010]. Informative marginal priors were used for all parameters but one, φ, for which a vague (poorly informative) prior was used instead, namely an uniform distribution U(-1, 1). This latter choice was done because prior knowledge for φ was unavailable. Instead, informative marginal prior distributions for other parameters could be determined on the basis of expert knowledge and well known reference values from [START_REF] Bruno Salm | Berenchnung von Fliesslanwinen, eine Anleitung für Paktiker mit Beispielen[END_REF]. In addition, a comprehensive prior sensitivity analysis was conducted (see Appendix C in Section 3.8). It demonstrates that our results are highly robust to the prior choice so that this question is no longer further considered in what follows.

0
σ 2 ∼ InvGamma(3, 20) ξ ∼ Γ(25, 40) φ ∼ U(-1, 1)

Convergence of the MH algorithm

For both models, we generated 50,000 MH samples and the last 25,000 were kept. The first 25,000 iterations were discarded as a burn-in period. This step avoids dependence on initial values. To further assess the convergence of the MH algorithm, we computed standard diagnoses [START_REF] Kuczera | Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm[END_REF][START_REF] Torre | The Metropolis-Hastings algorithm, a handy tool for the practice of environmental model estimation: illustration with biochemical oxygen demand data[END_REF][START_REF] Robert | Introducing Monte Carlo Methods with R (Use R)[END_REF], Eckert et al., 2010]. The [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] package coda created by [START_REF] Plummer | CODA: Convergence Diagnosis and Output Analysis for MCMC[END_REF] was used to calculate some of these. Specifically, three parallel chains were generated starting from different initial points and it was graphically checked that they were mixing well enough (not displayed). Also, the Kolmogorov-Smirnov test (H 0 : the two samples where drawn from the same continuous distribution) was applied to the samples from the different chains and results showed that these were indeed drawn from the same distribution. In addition, for both models and all parameters, it was checked that acceptance rates (defined as the number of candidates accepted over the total number of iterations of the MH algorithm) were close enough to 0.25, the optimal value according to [START_REF] Robert | The Metropolis-Hastings algorithm[END_REF] to reach fast convergence in a random walk. Finally, the Gelman-Rubin convergence criterion was computed [START_REF] Gelman | Inference from Iterative Simulation Using Multiple Sequences[END_REF]. It is based on the difference between a weighted estimator of the variance and the variance of the estimators from the different chains [START_REF] Robert | Introducing Monte Carlo Methods with R (Use R)[END_REF]. In our case, it equals 1 for all parameters, a perfect result. This all demonstrates that convergence was clearly reached for both models applied to the case study, leading to meaningful numerical approximations of the joint posterior distributions and, hence, reliable posterior estimates of the parameters of both models and of the Bayes factor.

Posterior distributions and posterior estimates and3.6 for M 0 and M 1 , respectively. For the two models, the marginal posterior distributions have lower variance than their priors, meaning that the observations have conveyed information into the analysis (see Figure 3.4), an expected result.

The variances of the posterior distribution of the friction parameters µ and ξ in the M 0 model are lower than in M 1 model. This result (arguably the sole undesirable one with M 1 compared to M 0 ) could be a consequence of the interactions between the friction coefficients and error parameters in the more parametrized M 1 model (four free parameters instead of three with M 0 ). Indeed, there is strong correlation between the parameters ξ and φ (0.49) for the model M 1 , which may preclude reaching sharp estimates of friction parameters with model M 1 . Conversely, even if there is a high correlation between the parameters µ and ξ with both models, switching from M 0 to M 1 reduces it from 0.85 to 0.50. The high correlation between the two parameters of the Voellmy friction law is known since the first calibration approaches of, e.g. [START_REF] Dent | Modeling of snow flow[END_REF]. This usually limits robust interpretation of obtained estimates, so that reducing it thanks to M 1 should be seen as advantageous. This correlation reduction may be a collateral effect of having one more free parameters, allowing a bit more flexibility to fit the data. Similar effect is reflected in the lower correlation between σ 2 and the friction parameters with M 1 than with M 0 . Note by the way that fairly assessing such correlations is a real strong point of our formal Bayesian calibration approach.

Remarkably, posterior estimates of friction parameters are very contrasted under both models. Even if posterior densities are not significantly different at the 95% credibility level (notably because of the higher a posteriori variance with model M 1 ), differences in posterior estimates reaches 13% for µ (relative difference between 0.34 with M 1 instead of 0.3 with M 0 ), and 70% for ξ (relative difference between 696 m 2 .s -1 with model M 1 and 410 m 2 .s -1 with model M 0 ). This indicates that using either M 0 or M 1 makes a huge difference if one aims at interpreting the value of posterior estimates, for, e.g., relating avalanche friction characteristics to snow and topographical conditions.

Another remarkable result is the diminution of the error variance σ 2 from model M 0 to model M 1 (Table 3.4). Indeed, the errors variance σ 2 is higher if the autocorrelation φ is not included. Specifically, the standard deviation of the errors decreases from 3.6 to 3.1 m.s -1 indicating that with model M 1 velocity predictions may be seen as 13% more accurate (relative difference between both estimates), another desirable property for practical use in snow science. Eventually, the autocorrelation parameter φ is largely positive, with a posterior mean of 0.71. Also, its posterior 95% credibility interval whose lover bound is 0.33 firmly excludes the zero value corresponding to M 0 . This pleads for a high and significant autocorrelation between velocities along flow propagation.

Model M 0 M 1 µ ξ[m.s -2 ] σ 2 [m 2 .s -2 ] µ ξ[m.s -2 ] σ 2 [m 2 .s -2 ] φ Mean 0.

Model posterior estimates

We calculated for model M 0 the posterior estimate of Voellmy model simulations f (Θ, x) from the posterior estimates of parameters µ and ξ (Table 3.4). Resulting posterior estimates of model errors t were analyzed by applying the Ljung-Box test (H 0 : the data are independently distributed). We found a p-value lower than 0.05 indicating that there is a significant autocorrelation of model errors. In other words, cor( t , t-1 ) for all t > 1 is significant. Hence, one of the assumptions underlying statistical model M 0 is not fulfilled by the data. On the other hand, applying standard statistical tests to posterior estimates of model errors shows no evidence of non-Gaussian errors (Shapiro test H 0 : the sample is normally distributed, p-value > 0.05) or heteroscedasticity. The assumptions of Gaussian errors of common variance underlying statistical model M 0 are thus fulfilled.

We proceeded similarly for statistical model M 1 , applying the same statistical tests to the posterior estimates of the innovations η t of the AR(1) model. Remember that, in the M 1 model, the independence assumption is on the innovations η t and not on the errors t . Test results show that innovations are indeed independent (p-value of the Ljung-Box > 0.05), and normally distributed (p-value of Shapiro test > 0.05) with common variance. These results show that M 1 represents correctly the autocorrelation of the errors, and, more widely, that contrary to model M 1 all of its underlying assumption are fulfilled, which promotes its use from a strict statistical point of view.

Predictive velocity distributions

To analyze how the statistical model choice affects velocity estimation, we propagated parameter uncertainty up to model predictions. Two sets of posterior predictive simulations were performed, the first integrating over the posterior distribution of friction parameters only, leading the posterior predictive distribution p(v x |v obs , x, M i ) of the avalanche propagation model, and the second integrating over the distribution of both friction parameters and model error parameters leading the full posterior predictive distribution of avalanche velocities for the case study p (v x |v obs , x, M i ). More precisely, the first writes:

p(v x |v obs , x, M i ) = f (Θ, x)p(Θ|v obs , x, M i )dΘ, (3.14)
and, the second one :

p (v x |v obs , x, M i ) = (f (Θ, x) + )p(Θ, γ|v obs , x, M i )dΘdγ. (3.15)
Eventually, to quantify the added value of parameter inference, the prior distribution of the velocity was also calculated by integrating over the prior distribution on the friction and model error parameters:

π(v x |x, M i ) = (f (Θ, x) + )π(Θ, γ|M i )dΘdγ.
(3.16)

By using a Bayesian approach, we obtained a sample of parameters (µ, ξ, σ 2 ) for M 0 (and also of φ for M 1 ). Then, we computed avalanche model simulations using these two samples to obtain several curves of velocity. Finally, we drawn the percentile of each curve in Figure 3.7. Also, this figure shows the resulting 90% credibility intervals. For both models, comparison between prior and posterior credible intervals shows that predicted velocities are logically shifted towards observations. Uncertainty reduction (change in the width of the 90% credible intervals) is stronger with M 0 , but also exists with M 1 . This simply reflects the larger a posterior variance of model M 1 parameters highlighted above.

In detail, almost all the elements of v obs , except two with model M 0 and one with model M 1 , are inside the 90% posterior credibility intervals drawn on p (v x |v obs , x, M i ) (blue dotted lines Figure 3.7), so that from this perspective the added value of M 1 is not obvious. However, why M 1 model is advantageous clearly appears if one focuses on the deterministic Voellmy model component (dotted orange lines in Figure 3.7): predictions come closer to observations and notably the velocity underestimation generally attributed to the Voellmy sliding block model [see e.g., Ancey andMeunier, 2004, Gauer, 2014] is reduced. Indeed, only 33% of the observations are inside the 90% posterior credibility intervals of the sliding block model simulations p(v x |v obs , x, M 0 ), compared to 62% with M 1 model.

Another quantification is provided by the systematic difference between both models that we evaluated through the mean difference and the mean quadratic difference between the medians of predicted velocities under M 0 (Figure 3.7 panel a), black line) and M 1 (Figure 3.7 panel b), black line). If we denote q50 vx,M i the median (50% percentile) of the posterior predictive distribution of velocities under model M i , the first writes:

∆ 1 = 1 240 240 0 (q50 vx,M 1 -q50 vx,M 0 )dx
and the second:

∆ 2 = 1 240 240 0 (q50 vx,M 1 -q50 vx,M 0 ) 2 dx
Obtained values are ∆ 1 = 1.31 m.s -1 and ∆ 2 = 1.72 m.s -1 . The positive value of the mean difference and the fact that it is not that much lower than the mean quadratic difference clearly shows that the inclusion of autocorrelation into the modelling truly leads to systematically higher velocity predictions. Predicted velocities are rather constant between x = 50m and x = 250m (plateau phase), around 15m.s -1 with M 0 . The underestimation with M 0 with regards to M 1 can be estimated by the ratio between ∆ 2 and this velocity to 11.5%.

All in all, the analysis of predictive velocity distributions further confirms that model M 1 should be preferred to get sharper, less underestimated, point estimates of avalanche velocities and to assess the related uncertainty fairly (credibility intervals lager but more likely to realistically describe the range of possible results).

Model selection

Finally, the logarithm of the Bayes Factor log 10 B 10 was calculated to compare the M 0 model to the M 1 model and we obtained a value of 2.02. According to the interpretation suggested in [Kass and Raftery, 1995, page 777] evidence, given by the data, to reject the M 0 model compared to the M 1 model is thus decisive. Indeed, log 10 B 10 > 2 simply implies that given the data model M 1 is more than 100 times more likely than model model M 0 .

Synthetic data

We then applied our Bayesian framework to calibrate the parameters of the 100 synthetic velocity profiles generated using the MAP estimators of model M 1 . We used the same marginal prior distributions as before (see Table 3.3). For each synthetic avalanche, we generated 20,000 MH samples and the last 10,000 were kept. Convergence was verified with the same diagnoses as for the measured data. We could then determine the ability of the M 0 and M 1 models to retrieve the true model parameter values used for the synthetic data generation. For this, we calculated the 90% coverage rates for the different parameters, this means, the number of 90% credibility intervals recovering the true value of Table 3.2. According to the results of Table 3.5, M 1 model has a much better ability to determine the true model parameter values compared to the M 0 model. This especially q q q q q q q q q q q q q q q q q q q q q q 0 10 20 a) q q q q q q q q q q q q q q q q q q q q q q 0 10 20 M 0 M 1 Parameter µ ξ σ 2 µ ξ σ 2 φ Coverage rate 64 84 0 98 82 98 93 Table 3.5: The 90% coverage rates for the synthetic data sample. For each parameter and both models, the 90% coverage rate corresponds to the number of times over the sample of 100 synthetic avalanches for which the 90% posterior credibility interval includes the true value used for data generation. In other words, 90% is the perfect validation score. Parameter φ applies to model M 1 only.

applies to σ 2 for which M 0 model is fully unable to identify the true value and, to a lower extent, to µ, for which the true value is in the posterior 90% credibility interval less than two times over three. Conversely, 90% coverage rates with M 1 model are rather fair, varying between 82% and 98% for all the four parameters.

To further compare the precision of the estimates led by the two models, the MAP estimators corresponding to the full synthetic sample are presented in Figure 3.8 in panels from a) to d). In mean, both models result in an underestimation of the µ parameter, but the underestimation with M 1 is much smaller than with M 0 . In addition, the true value is well within the range of variability of the different estimates corresponding to the 100 synthetic avalanches with M 1 whereas it is clearly outside with M 0 . Similarly, parameter σ 2 is much better estimated by M 1 model. Parameter ξ is slightly underestimated by M 0 model and slightly overestimated by M 1 model, but both models perform reasonably. Finally, M 1 model estimates correctly the autocorrelation parameter φ. Overall, this analysis confirms that only M 1 model is able to retrieve the true parameter values as soon as autocorrelation between measurements actually exists. In other words, not accounting 3.5. Conclusion and outlook q q q q q q 0.2 0.3 0.4 0.5

M0 M1

Model µ a) q q q q q q q q 200 400 600 800 1000 M0 M1

Model

ξ ( ms -2 ) b) q q q q q q 0 50 100 M0 M1 Model σ 2 ( m 2 s -2 ) c) q q q 0.0 0.5 1.0 M1 Model φ d)
q q q q q q 10 20 30

Num. Iterations for autocorrelation within model calibration when autocorrelation actually exists carries high risk to lead to biased estimates. This provides another very strong argument in favour of M 1 since, as evidenced by the case study, significant autocorrelation may indeed exist between measurements made along the same avalanche flow.

Eventually, the common logarithm of the Bayes factor calculated for each synthetic avalanche (see Figure 3.8 panel e) indicates for all synthetic velocity series decisive evidence to reject the M 0 model compared to the M 1 model. This result is all the more logical given that the Bayes factor asymptotically selects the true model when it is included within a sample of competing models. However, this can be seen as a last strong point to advise using model M 1 as soon as autocorrelation is suspected.

Conclusion and outlook

In this work, a Bayesian calibration of an avalanche flow dynamics model from data of high temporal resolution [1s] was developed. The objective of this work was to show how potential autocorrelation between measurements made along an avalanche flow can be considered within the calibration of an avalanche model, and to demonstrate that this improves the accuracy of friction parameters estimation and velocity predictions.

Two statistical models representing the discrepancies between observations and simulations were proposed: the first one, M 0 , classically considered the errors as independent and identically normally distributed, and the second one, M 1 , modelled the errors as an autoregressive process of order 1. The latter accounts for potential autocorrelation between measures made along the same avalanche flow, a question which has been poorly addressed in the snow science literature so far. Our objective was to determine the accuracy of parameter estimation when the autocorrelation is included or not in the modelling. More generally, we wanted to develop a framework able to link in a more mathematically consistent way, the rich and diverse data now increasingly available from avalanche testsites with numerical propagation models.

Application was made on a well documented avalanche event and on synthetic data. The avalanche was released at the Lautaret full-scale test-site on February 13 2013 and the corresponding velocity time series was obtained from high rate positioning photogrammetry. A synthetic data set of 100 velocities time series mimicking the avalanche was created in order to further test our approach and to illustrate its benefits. Results for the avalanche showed strong and significant autocorrelation between predicted velocities. It also appeared that the velocity time series was correctly modelled by the statistical model M 1 only (Ljung-Box test p-values > 0.05) so that there was decisive evidence to reject M 0 compared to M 1 . In addition, resulting posterior estimates for friction parameters µ and ξ and velocities along the path were shown to be different and arguably more realistic (for the estimated velocities) with M 1 than those obtained with M 0 . Eventually, our synthetic data confirmed that in presence of autocorrelation between measurements, not accounting for it may lead to biased estimates. This all demonstrates the necessity and usefulness to explicitly account for autocorrelation within the calibration, both to realistically predict avalanche characteristics and to get friction parameter estimates that can be related to snow and topographical characteristics in a meaningful way, two points of crucial importance in snow science. Hence, since avalanche velocities are time-series, as a good practice of modelling, potential autocorrelation of the errors should always be envisaged, since we clearly demonstrate that neglecting this effect may lead to undesirable consequences.

To develop and illustrate our approach, we chose the simple sliding block model with the Voellmy friction law because it is faster to run than more computationally intensive state of the art avalanche models. As it was mentioned in the Introduction, we used this simple model as a proof of concept but our approach could be generalized to other friction laws (e.g., the Coulomb one which has the advantage of having a single parameter (µ)), and to other ways of describing avalanche flows (e.g., depth-averaged equations or a sliding block with a PCM formulation, which would avoid the estimation of the flow depth) after suitable adaptation of the algorithm. Also, on the existing basis, it should be possible to include within the calibration other quantities which can be measured simultaneously as for example, the snow depth in the release area, the runout distance, the deposited volume, etc. This would better constrain the calibration, leading to potentially less uncertain posterior estimates. In such a case, an appropriate likelihood would need to be proposed and a higher number of observations would be required to calibrate the new likelihood parameters. To overcome this shortcoming using several avalanches altogether could be an option which was out of the scope of this work.

One of the main drawbacks of a Bayesian approach is the high number of simulations that are needed (at least 1000 simulations) to reach convergence in a Markov Chain Monte Carlo setting. Also, a compromise between the number of parameters that can be estimated and the size of the data must be found. Hence, whether or not our approach can be practically implemented with the most computationally intensive avalanche model and how many unknown parameters can be inferred as function of the available data sets remain questions to be investigated.

Finally, it is worth to mention that the avalanche size and velocities we found are typical of a medium size avalanche (class 3 on the CAA international avalanche scale). Specifically, friction parameter values correspond to the range of values that can be found in the literature. In future work, it could therefore be informative to apply our approach to a much large sample of avalanches to exploit its potential for inferring relevant physics (e.g., to find the relationship between friction parameters and snow conditions). Ultimately, this all should i) lead to a more accurate evaluation of the highest-return-period events required for avalanche risk assessment, and ii) simultaneously improve our knowledge of the relevant physics by providing sharper quantification of underlying processes.

Appendix A: Likelihood of an AR(1) process

Since |φ| < 1, the t error has a representation in terms of the innovations η t (also known as MA(∞) representation, see [START_REF] Gouriéroux | Séries temporelles et modèles dynamiques[END_REF] for more details):

t = φ t 0 + t-1 i=0 φ i η t-i , (3.17)
where we recall that φ and η t-i are the coefficient and the innovations of the AR(1) process, respectively (see equation 3.4).

From the last equation, the unique stationary solution is obtained for 0 with zero mean and variance equal to σ 2 1-φ 2 , independent from (η 1 , . . . , η t ), and we get that:

E( t ) = E( 0 ) + t-1 i=0 φ i E(η t-i ) = 0, (3.18) Var( t ) = φ 2t Var( 0 ) + t-1 i=0 φ 2i Var(η t-i ) = σ 2 (1 -φ 2 ) . (3.19)
The marginal distribution of t is Gaussian because it is the sum of independent and identically normally distributed η t . In particular the first error term 1 is Gaussian:

1 ∼ N 0, σ 2 (1 -φ 2 ) . (3.20) Considering that 2 = φ 1 + η 2 , the conditional distribution of 2 given 1 is 2 | 1 ∼ N (φ 1 , σ 2 ).
In a more general way:

t | t-1 ∼ N (φ t-1 , σ 2 ).
(3.21)

Then, the joint distribution of errors is:

p( 1 , • • • , n |φ, σ 2 ) = p( 1 ) n i=2 p( i | i-1 ). (3.22)
From this equation, the likelihood expression of equation (3.7) is obtained.

Appendix B: Numerical evaluation of the Bayes Factor

Thanks to the Bayes theorem, we can write:

1 p(v obs |M i , x) = π(ψ i |v obs , M i , x) L(v obs |ψ i , M i , x)π(ψ i |M i , x) , (3.23)
where i ∈ {0, 1}. If g is a density function defined on an ensemble Ω:

1 p(v obs |M i , x) = 1 p(v obs |M i , x) ψ i ∈Ω g(ψ i )dψ i = ψ i ∈Ω g(ψ i )π(ψ i |v obs , M i , x) L(v obs |ψ i , M i , x)π(ψ i |M i , x) dψ i .
From Monte Carlo simulations, x) can thus be estimated as:

1 p(v obs |M i ,
1 p(v obs |M i , x) = 1 N N j=1 g(ψ (j) i ) L(v obs |ψ (j) i , M i , x)π(ψ (j) i |M i , x)
.

(3.24)

where {ψ [START_REF] Kass | Bayes factors and model uncertainty[END_REF] page 19 equation 12).

(j) i ; j = 1, • • • , N } is a sample from the posterior distribution of M i (see
For this application, the g function was chosen as a multinormal distribution with mean equal to the empirical mean and covariance matrix equal to the estimated covariance matrix of our Metropolis-Hastings sample, respectively.

Appendix C: Prior sensitivity analysis

To study the robustness of our results, we conducted a prior sensitivity analysis by varying prior information for model parameters (µ, ξ, σ 2 ). We explored the range of priors from the marginal informative priors used in the paper core up to vague (poorly informative) priors as much as possible. In the case of M 0 model, marginal vague priors could be tested for all the three parameters instead of informative ones. For M 1 model, however, the informative prior should be always kept for ξ, because, if for this parameter a vague prior was used, the MH algorithm convergence could not be achieved. This is explained by the high correlation between the ξ and φ parameters (close to 0.5, see Figure 3.6). For the φ parameter, a uniform distribution U(-1, 1) was used in all the analysis. The different combinations of informative and non-informative priors tested are shown in Table 3.6. Corresponding posterior distributions are shown in Figure 3.9.

Even if slight differences between the different posteriors can be noted as one replaces one prior by another, overall, for both models, results remain quite similar, leading in all cases to posterior estimates very close to the ones obtained with the informative priors used in the paper core (see Table 3.4). Also, the common logarithm of the Bayes Factors was calculated between all combinations of M 1 and M 0 posterior samples and, in nearly all cases, there is strong to decisive evidence to reject M 0 compared to M 1 model, as demonstrated in the paper core with informative priors. Indeed, the obtained log 10 (B 10 ) values are between 0.57 and 2.49 with a mean of 1.38 and a standard deviation of 0.61. Hence, all in all, posterior inference and model selection show little sensitivity to prior specification. Since it is advised to use priors with a finite domain as much as possible to avoid spurious results [START_REF] Kuczera | Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm[END_REF]] we kept the informative priors in the paper core. 3.9 Appendix D: Images used in the photogrammetric process

In this appendix, some examples of images used in the photogrammetric process are shown.

Part II

Global sensitivity analysis Abstract: Avalanche models are increasingly employed for elaborating land-use maps and designing defense structures, but they rely on poorly known parameters. Careful uncertainty assessment is thus required but difficulty arises from the nature of the outputs of these models, which are commonly both functional and scalar. Hence, so far in the avalanche field, few sensitivity analyses have been performed. In this work, we propose to determine the most influential inputs of an avalanche model by estimating aggregated first-order Sobol' indices. We propose a nonparametric estimation procedure based on the Nadaraya-Watson kernel smoother, which allows to estimate the aggregated Sobol' indices from a given random sample of small to moderate size. Due to the limited size of the sample, the kernel estimation is biased. Therefore, we propose a bootstrap based bias correction before selecting the bandwidth by cross-validation. To estimate the aggregated Sobol' indices, we reduce the dimension of the output using principal components analysis. After different test cases showing the efficiency of our approach, it is applied to a real avalanche case. Results show that the friction parameters and the snow depth in the release zone are the most influential parameters determining the avalanche characteristics.

Introduction

In mountain environments, snow avalanches are a constant threat for settlements and their inhabitants [McClung andSchaerer, 2006, Ancey, 2006]. This creates conflict between development and, safety and land-use planning is an efficient way to reduce death tolls [START_REF] Gruber | Winter 1999: a valuable test of the avalanche-hazard mapping procedure in Switzerland[END_REF]Margreth, 2001, Eckert et al., 2018a]. To this aim, avalanche flow dynamic models are increasingly employed for elaborating land-use maps and designing defense structures [START_REF] Naaim | Return period calculation and passive structure design at the taconnaz avalanche path, france[END_REF], Favier et al., 2014b]. However, all existing physics-based snow avalanche models remain based on some parameters which are poorly known [Ancey andMeunier, 2004, Eckert et al., 2008b]. This applies, for instance, to friction parameters representing the behavior of snow in motion and to initial conditions corresponding to the avalanche release [START_REF] Jamieson | Application and limitations of dynamic models for snow avalanche hazard mapping[END_REF], Eckert et al., 2010, Naaim et al., 2013]. Thus, a careful assessment of the impact of the uncertainty of the inputs on the outputs should be carried out [START_REF] Fischer | A novel approach to evaluate and compare computational snow avalanche simulation[END_REF]. Difficulty arises from the nature of the outputs of these models, which are commonly both functional (e.g., the velocity and flow depth of the avalanche as function of space and time) and scalar (e.g., the runout distance, aka the point of further reach on the 2D topography of the avalanche flow path).

One of the main purposes of sensitivity analysis is to determine the inputs which are the most influential on the output or outputs of a model where the model can be of any nature, for example a black-box model or a complex computational code [START_REF] Saltelli | Sensitivity Analysis: Gauging the Worth of Scientific Models[END_REF]. More specifically, [START_REF] Saltelli | Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models[END_REF] proposed four settings as a guide of objectives of a sensitivity analysis: factor prioritization, factor fixing, variance cutting and factor mapping [Saltelli et al., 2004, Iooss andPrieur, 2017]. It depends on the analyst to define the objectives of its sensitivity analysis. In the global sensitivity (GSA) framework, the input parameters are modeled by random variables, which will be assumed independent in this paper. The probability distribution chosen to model the input vector is often guided by practitioner's belief. The output is then random as it depends on the inputs through the model. In this work, our objective is factor prioritization which consists in identifying which inputs or factors once fixed would reduce the variance of the output [START_REF] Saltelli | Global sensitivity analysis. The primer[END_REF] at most.

There exist different sensitivity measures in the literature for quantifying the influence of each input on the output: variance based indices, also known as Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF], density based measures [Borgonovo et al., 2016], entropy-based sensitivity measures [START_REF] Auder | Global sensitivity analysis based on entropy[END_REF], etc. A detailed review of sensitivity measures can be found, e.g., in [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF] and Borgonovo and Plischke [2016].

We focus our study on variance based first-order Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF]] commonly used for factor prioritization. In short, Sobol' indices are constructed as follows: if the output of a model is scalar, its total variance can be split into partial variances by using the Hoeffding decomposition [START_REF] Hoeffding | A Class of Statistics with Asymptotically Normal Distribution[END_REF]. If the inputs are independent, each of the partial variances is associated to an input or to an interaction between inputs through the ANOVA decomposition [START_REF] Efron | The jackknife estimate of variance[END_REF]. Then, the Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] are calculated as the ratio of each partial variance to the total variance. Moreover, Sobol' indices generalize the linear correlation coefficient R 2 for any kind of models.

Sobol' indices are easy to interpret: a high index value shows the input is relevant and a value close to zero shows that it is not, all the indices are normalized between 0 and 1 and sum to 1. Moreover, Sobol' indices have a natural extension to multivariate or functional outputs [see, e.g., [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF], Gamboa et al., 2013] which is the framework of the application that motivated this study.

If the output is multivariate, it is possible to apply sensitivity analysis to each component of the multivariate output but this could lead to redundancies in the results, particularly in the setting of discretized functional outputs. It seems therefore interesting to turn to the so-called aggregated indices first introduced in [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] which summarize the information. A preliminary step of output dimension reduction such as, e.g., principal component analysis (PCA) or partial least-squares may also be applied first [see, e.g., [START_REF] Campbell | Sensitivity analysis when model outputs are functions[END_REF]. For our real application in the avalanche field, dimension reduction was performed by simultaneous PCA introduced in Ramsay and Silverman [2005] [see also Nanty et al., 2017, and references therein].

We now discuss the estimation of aggregated Sobol' indices. There are many methods to estimate Sobol' indices: the Fast Amplitude Sensitivity Test (FAST) [see, e.g., Saltelli et al., 1999, and references therein], Random Balance Design [START_REF] Tarantola | Random balance designs for the estimation of first order global sensitivity indices[END_REF]], Sobol' pick-freeze schemes [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF]. The main drawback of all the aforementioned methods [START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF] is that they are based on sampling designs of particular type. Nevertheless, for many memory and time consuming real applications, the cost (in terms of number of model evaluations) of these approaches is prohibitive. For example, to estimate a single Sobol' index with an uncertainty of 10%, it could be required to perform 10 4 model runs [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF].

To overcome the drawback of the aforementioned methods, authors use approaches not based on structured sampling designs to estimate the indices. These methods are known as given data or one sample methods. They also correspond to green sensitivity analysis because available data from previous model runs can be reused. Among these approaches, we can cite the effective algorithm for computing global sensitivity indices (EASI) method proposed by [START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF], which is a spectral method based on the Fast Fourier Transform. We also mention the work in [START_REF] Plischke | Global sensitivity measures from given data[END_REF], which relies on the notion of class-conditional densities, where a class is a sub-sample stemming from a suitable partition of the dataset. Much more recently, Antoniano-Villalobos et al. [2019] proposed a fully Bayesian given data procedure.

In the present paper, we propose a new given data method to estimate first-order Sobol' indices, based on a nonparametric Nadaraya-Watson [START_REF] Nadaraya | On estimating regression[END_REF][START_REF] Watson | Smooth regression analysis[END_REF] bias corrected kernel regression [START_REF] Racine | Bias-Corrected Kernel Regression[END_REF]. Our approach is close to the one introduced in [START_REF] Solís | Non-parametric estimation of the first-order Sobol indices with bootstrap bandwidth[END_REF]. The main difference is in the formula we apply to correct the bias, which seems to be more efficient, at least in the results of our simulation study. We then extend this estimation procedure to the estimation of the aggregated Sobol' indices [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF]. To our knowledge, it is the first time a nonparametric given data procedure is proposed in this framework. It must be noted that the estimation of higher-order and total Sobol' indices may not be possible using kernel estimation with moderate sample sizes due to the curse of dimensionality Da Veiga et al. [2009], Hence, the estimation method proposed in this work focuses on aggregated first-order Sobol' indices. Kernel regression with bandwidth selection via cross-validation was already used in [START_REF] Doksum | Nonparametric estimation of global functionals and a measure of the explanatory power of covariates in regression[END_REF] to estimate the relative measure of the importance of a subset of input variables relative to the full set of input variables. It is also of interest to cite [START_REF] Sparkman | Importance Sampling-based Post-Processing Method for Global Sensitivity Analysis[END_REF], in which the authors proposed a sample weighting scheme similar to kernel regression for estimating the moments of the conditional distributions in the calculation of first-order and total Sobol' indices. Da Veiga et al. [2009] also proposed a methodology based on local linear regression for the estimation of first-order Sobol' indices. The Nadaraya-Watson estimator can be seen as a particular case of the wider class of nonparametric estimators, called local polynomial estimators used in [START_REF] Da Veiga | Local polynomial estimation for sensitivity analysis on models with correlated inputs[END_REF]. Specifically, it corresponds to performing a local constant fit. However the procedure in [START_REF] Da Veiga | Local polynomial estimation for sensitivity analysis on models with correlated inputs[END_REF] is clearly different from ours, as it consists in splitting the sample in two disjoint sub-samples, the first one used to emulate the conditional expectation, the second one used to estimate each first-order Sobol' index, replacing the conditional expectation with its corresponding emulator. Note that metamodels have been widely used to estimate Sobol' indices (e.g., linear regression models, polynomials, splines). Most of the time, the metamodel is built to emulate f rather than the conditional expectation. Some of these metamodels, such as polynomial chaos (see, e.g., [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF], Le [START_REF] Gratiet | Metamodel-based sensitivity analysis: Polynomial chaos expansions and Gaussian processes[END_REF]), lead to an analytical expression of Sobol' indices (using, e.g., Parseval equality if polynomial chaos are considered). In [START_REF] Marrel | Calculations of sobol indices for the gaussian process metamodel[END_REF], the authors derive an analytical expression of Sobol' indices based on Gaussian process regression. In [START_REF] Durrande | Anova kernels and rkhs of zero mean functions for model-based sensitivity analysis[END_REF], the authors prone the use of ANOVA kernels. In the framework of sensitivity analysis with a nonparametric estimator based on nearest neighbors, [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF] proposed to estimate Shapley effects introduced by [START_REF] Owen | Sobol' Indices and Shapley Value[END_REF].

Within the avalanche field, only a few studies have analyzed the influence of inputs on the outputs of dynamic propagation models. For example, [START_REF] Barbolini | Estimate of uncertainties in avalanche hazard mapping[END_REF] used a Monte Carlo approach to analyze the sensitivity of runout distances and impact pressures in the VARA model. [START_REF] Jamieson | Application and limitations of dynamic models for snow avalanche hazard mapping[END_REF] discussed the main sources of uncertainty in the inputs of avalanche models. [START_REF] Borstad | Sensitivity analyses in snow avalanche dynamics modeling and implications when modeling extreme events[END_REF] developed a sensitivity analysis of an avalanche model with a Coulomb-type friction law. [START_REF] Bühler | Sensitivity of snow avalanche simulations to digital elevation model quality and resolution[END_REF] explored the influence of the digital elevation models resolution on the outputs of the RAMMS avalanche model. Eventually, [START_REF] Buhler | Sensitivity Analysis of the RAMMS Avalanche Dynamics Model in a Canadian Transitional Snow Climate[END_REF] developed a sensitivity analysis of released volumes, runout distances and avalanche velocities in the runout zone with respect to the initial released volume with the RAMMS avalanche model. However, the previous studies did not apply formal statistical methods to quantify the respective importance of the inputs. Moreover, authors considered only scalar outputs in their approaches. By contrast, in our paper, we focus on the avalanche model proposed by [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF] and we quantify the importance of its inputs on its outputs by estimating the aggregated Sobol' indices with our new nonparametric procedure. The uncertainty on the input parameters was defined based on data obtained from an avalanche released on 13 February 2013 at the Lautaret full-scale test-site [START_REF] Thibert | The full-scale avalanche test-site at Lautaret Pass (French Alps)[END_REF].

In summary, the aims of this study are: (i) to propose a nonparametric estimation method for the aggregated Sobol' indices in a given data framework when sample size is small (ii) to quantify the input importance in avalanche models, having complex outputs (e.g., a mix of functional and scalar outputs). The approach can be easily adapted to other avalanche models and more widely in environmental sciences in the frequent case of complex models with outputs which are both functional and scalar. This work is organized as follows: the aggregated Sobol' indices and the estimation method are described in Section 4.2. In Section 4.3, we test the estimation method on toy functions. Then, in Section 4.4, the avalanche model is described and the results are presented. Finally, in Section 4.5, the conclusions and perspectives are discussed.

Aggregated Sobol' indices

Let us denote by f the model which takes as inputs the vector X = (X 1 , . . . , X d ). The inputs X 1 , . . . , X d are modeled by random variables defined on a probability space (Ω, F, P) and valued in a measurable space E = E 1 × E 2 . . . × E d . The output of the model f is the p multivariate vector Y = (Y 1 , . . . , Y p ) T . It means, we have:

f (X) = f (X 1 , . . . , X d ) = Y = (Y 1 , . . . , Y p ) T .
In the following, we assume that E(Y 2 1 + . . . + Y 2 p ) < ∞ and that the random variables X 1 , . . . , X d are independent from each other. This condition guarantees the uniqueness of the ANOVA decomposition [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF]. To summarize the importance of each input X i on the multivariate output Y, we aim at computing the aggregated Sobol' index GS i introduced in Lamboni et al. [2009] [see also Gamboa et al., 2013] defined as:

GS i = p j=1 Var(Y j )S j i p j=1 Var(Y j )
, where S j i is the first-order Sobol' index of the output Y j with respect to the input X i namely

S j i = Var(E(Y j |X i )) Var(Y j ) (4.1) Notice that GS i ∈ [0, 1].
As in the case of the scalar indices, the main advantage of the aggregated indices is their easy interpretation: a high value means the input is important, a value close to zero means it is not, and d i=1 GS i = 1. In the following, we propose to estimate nonparametrically GS i from given data, proposing a new procedure based on the Nadaraya-Watson kernel smoother [START_REF] Nadaraya | On estimating regression[END_REF][START_REF] Watson | Smooth regression analysis[END_REF].

Nonparametric estimation procedure

Let (X i , Y j ) 1≤ ≤n be an independent identically distributed random sample of (X i , Y j ). The Nadaraya-Watson kernel smoother consists in estimating the conditional mean of E(Y j |X i ) at a point x in the domain of X i by:

m h (x) = n =1 Y j W ,h (x), (4.2)
where

W ,h (x) = K x-X i h n =1 K x-X i h , (4.3)
with K a kernel function and h a bandwidth. In Equation (4.3), different weights W ,h can be used, as for example, local linear regression proposed in [START_REF] Da Veiga | Local polynomial estimation for sensitivity analysis on models with correlated inputs[END_REF] and nearest neighbor smoothers as proposed by [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF].

When the sample size is small, the selection of the bandwidth h is critical. Indeed, to avoid undersmoothing or oversmoothing, an optimal value of h which balances bias and variance must be found [Tsybakov, 2008, page 17]. An option widely used to select h is cross-validation but if the sample size is small, the estimator will suffer of finitesample bias [see [START_REF] Racine | Bias-Corrected Kernel Regression[END_REF][START_REF] Solís | Non-parametric estimation of the first-order Sobol indices with bootstrap bandwidth[END_REF]. Therefore, if the sample size is small, it is recommended to perform a bias correction to the kernel smoother estimator preliminary to the bandwidth selection.

It is important to underline that the criterion we introduce here to select the bandwidth is different from the one proposed in [START_REF] Solís | Non-parametric estimation of the first-order Sobol indices with bootstrap bandwidth[END_REF], which was based on bagging minimization. The selection of the bandwidth we propose to apply in our paper is described in the next section.

Once the bandwidth h has been chosen, the estimation of the aggregated Sobol' index is straightforward: given a n sample (X k i , Y k j ) 1≤k≤n of (X i , Y j ),

GS i = p j=1 Var (Y k j ) 1≤k≤n Ŝj i p j=1 Var (Y k j ) 1≤k≤n , (4.4)
where Ŝj i is the estimation of the first-order Sobol' index:

Ŝj i = Var m h,-k (X k i ) 1≤k≤n Var (Y k j ) 1≤k≤n , (4.5)
with Var denoting empirical variance and

m h,-k (X k i ) = =k y j K X k i -X i h =k K X k i -X i h • (4.6)
Note that the same sample (X k i , Y k j ) 1≤k≤n is used for both the estimation of E(Y j |X i ) and the outer variance. This is the reason why we use the so-called leave-one-out estimate m h,-k and not the estimate m h defined in (4.2).

Bias correction and bandwidth selection

As mentioned previously, cross-validation is one option to select the bandwidth when using the Nadaraya-Watson kernel smoother defined in (4.2). In more details, given a n sample X = (X 1 , . . . , X n ) of X i and the corresponding evaluations Y = (Y 1 , . . . , Y n ) of Y j , the bandwidth can be selected as (see, e.g., Takezawa [2005]):

h cv = argmin h>0 1 n n k=1 Y k -m h (X k ) 2 1 -W k,h (X k ) (4.7)
with W k,h (X k ) defined by (4.3). For notation simplicity, K h (.) denotes K . h and we dropped all the i and j subindices, but keep in mind that we refer to a sample of X i and Y j .

Even if the kernel smoother with cross-validation is asymptotically unbiased, it is well known that it has a relatively large finite-sample bias. To tackle this issue, [START_REF] Racine | Bias-Corrected Kernel Regression[END_REF] proposed a bootstrap bias correction as a preliminary step to bandwidth selection. In this paper, we adapt the procedure introduced in [START_REF] Racine | Bias-Corrected Kernel Regression[END_REF] to the framework of Sobol' index estimation. The steps of our procedure are described hereafter:

1. Given a n sample X = (X 1 , . . . , X n ) of X i , the corresponding evaluations Y j = (Y 1 , . . . , Y n ) of Y j and a bandwidth h > 0 (for example, an initial bandwidth value could be h = h cv given by Equation (4.7)), we calculate the leave-one-out Nadaraya-Watson kernel regression estimate m h,-k (X k ) defined in (4.6). If

=k K h (X k - X ) = 0, we set m h,-k (X k ) = 0.
2. We compute the errors:

ˆ k = Y k -m h,-k (X k ), k = 1, . . . , n.
These errors may be heteroscedastic because Y j depends also on X 1 , . . . , X i-1 , X i+1 , . . . , X d , therefore, a standardization of the errors must be applied before bootstraping.

3. The errors {ˆ k } 1≤k≤n are standardized as:

υk = ˆ k -¯ r(X k ) , k = 1, . . . , n
where ¯ is the empirical mean of {ˆ k } 1≤k≤n and r(X k ) is the estimate of the conditional standard deviation of the random variable |X i at X k . Note that r(X k ) is estimated using nonparametric kernel smoother estimation. More precisely, we estimate r 2 (x) by applying the Nadaraya-Watson kernel smoother as defined in (4.2) with Y = 2 , = 1, . . . , n and selecting the bandwidth via cross-validation defined in (4.7). h,-k as in (4.6). 6. The bias of m h,-k (X k ) is estimated as:

A bootstrap sample Y

(b) of Y is created as: Y (b)k = m h,-k (X k ) + r(X k )υ (b)k , k =
1 B B b=1 m (b) h,-k (X k ) -m h,-k (X k ).
Then, the corrected kernel smoother m c h,-k (X k ) is computed as:

m c h,-k (X k ) = m h,-k (X k ) - 1 B B b=1 m (b) h,-k (X k ) -m h,-k (X k ) (4.8) = 2 m h,-k (X k ) - 1 B B b=1 m (b) h,-k (X k ).
(4.9)

Note that the last two equations correspond to the standard bootstraping bias correction proposed by [START_REF] Efron | A leisurely look at the bootstrap, the jackknife, and cross-validation[END_REF].

7. Finally, the bandwidth h is selected as

h boot = argmin h>0 1 n n k=1 Y k -m c h,-k (X k ) 2 .
We then estimate the aggregated Sobol' index GS i with (4.4) by replacing m h,-k (X k i ) by m c h boot ,-k (X k i ) in (4.5). [START_REF] Racine | Bias-Corrected Kernel Regression[END_REF] also proposed an iterative version of the bootstrap bias corrections. In this work, we used the non iterated version because a major drawback of the iterative one is its computational cost. A simplified version of the iterative procedure was run on a simple test case (see Figure 4.13 in the Appendix) The improvements were not that convincing with regards to the additional cost.

Following [START_REF] Da Veiga | Local polynomial estimation for sensitivity analysis on models with correlated inputs[END_REF], another approach could have been (i) to estimate Var(Y ) as far as to define the kernel smoother m h defined in ( 2), (ii) to perform the bias correction and selection of the bandwidth, and (iii) then to sample from the joint distribution of input vector X to estimate Var(E(Y j |X i )) using the bias corrected kernel smoother as a metamodel. By contrast, the method proposed is metamodel-free and belongs to the class of nonparametric given data (one sample) methods. In the following, our estimator with bias correction is denoted by cnp and the same estimator without bias correction, for which the bandwidth is selected by solving argmin h>0

1 n n k=1 Y k -m h,-k (X k ) 2 , is denoted by np.

Dimension reduction based on principal component analysis

If the output is high-dimensional or even functional, it may be computationally interesting to reduce the output dimension in a preliminary step. There are different reduction techniques. Principal component analysis, also known as Karhunen-Loève decomposition in the functional framework, consists in projecting the output on a new basis so that most information is concentrated in the first few components [see, e.g., [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF][START_REF] Loève | Probability Theory[END_REF]. Note that if the outputs are the discretization of more than one functional random variable (e.g., velocity and flow depth), simultaneous PCA may provide efficient reduction. For more details, see e.g., the work of [START_REF] Nanty | Uncertainty quantification for functional dependent random variables[END_REF] and reference therein.

In brief, simultaneous PCA consists in applying PCA to the matrix composed by the concatenation of the discretized functions divided by a normalization factor. [START_REF] Nanty | Uncertainty quantification for functional dependent random variables[END_REF] show that using the maximum of each functional as normalization factor provides the best results. Thus, in the avalanche application, we follow their procedure.

There exist dimension reduction tools specific to functional data, such as the regularized functional PCA described in Chapter 9 of [START_REF] Ramsay | Functional Data Analysis[END_REF] [see also the work of [START_REF] Yao | Functional data analysis for sparse longitudinal data[END_REF] which has been applied in different fields [e.g., [START_REF] Locantore | Robust principal component analysis for functional data[END_REF][START_REF] Antoniadis | Spatiotemporal metamodeling for West African monsoon[END_REF]. Also, it has been shown that in some applications, functional PCA provides better results than usual PCA applied after a discretization of the functional data [e.g., [START_REF] Viviani | Functional principal component analysis of fMRI data[END_REF]. On our test cases, usual PCA applied on the discretized data set and functional PCA provided similar results. The results we show in this paper are the ones obtained by applying usual PCA on discretized data.

The PCA consists in decomposing the variance-covariance matrix of the output. More precisely, let us denote by Σ the variance-covariance matrix of the output vector Y. The principal component decomposition of Y is based on the expansion:

Σ = p j=1 λ j v j v T j ,
where λ 1 ≥ . . . ≥ λ p ≥ 0 are the ordered eigenvalues and v 1 , . . . v p are the orthonormal eigenvectors of Σ. Therefore, Y can be decomposed as:

Y = E(Y) + p j=1 (Y -E(Y)) T v j v j = E(Y) + p j=1 H j v j
where H j is the jth principal component of Y. Note that Var(H j ) = λ j . The output can be approximated using the q ≤ p first components which capture the major part of the output variance:

Y ≈ E(Y) + q j=1 H j v j .
Then, the aggregated Sobol' indices can be computed in the first q principal components. For example, for input X i :

GS i ≈ q j=1 λ j S i (H j ) q j=1 λ j (4.10)
where S i (H j ) denotes the Sobol' index of H j with respect to input X i . The estimation of the aggregated Sobol' indices will be done for the reduced output H = (H 1 , . . . , H q ) T .

Test cases

The cnp method is tested on three toy functions and then, it is applied to an avalanche dynamic model. The cnp accuracy is compared with np (bandwidth selected with crossvalidation and without bias correction), the method of [START_REF] Solís | Non-parametric estimation of the first-order Sobol indices with bootstrap bandwidth[END_REF] denoted by Solis2019 and implemented in R in the library sobolnp, and the method of [START_REF] Tissot | Bias correction for the estimation of sensitivity indices based on random balance designs[END_REF], denoted by Tissot2012. [START_REF] Tissot | Bias correction for the estimation of sensitivity indices based on random balance designs[END_REF] proposed a bias correction for the index estimation based on random balance designs. To perform step 3 of the bias correction algorithm, we used the function npreg of the R package np developed by [START_REF] Hayfield | Nonparametric econometrics: The np package[END_REF].

To compare our method to the one proposed in [START_REF] Solís | Non-parametric estimation of the first-order Sobol indices with bootstrap bandwidth[END_REF], we use the Epanechnikov kernel of order 2:

K(u) = 3 4 (1 -u 2 )I (|u|≤1) ,
where I denotes the indicator function. We compared the estimation results with the ones obtained with a Gaussian kernel and the results were quite similar (see Figure 4.14 in the Appendix). All the computations were performed in [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF].

Toy functions

For each toy function, to test the accuracy of the method, first-order Sobol' indices were estimated with N = 100 independent samples of sizes n ∈ {100, 200, 300}. For cnp and Solis2019, the number of bootstrap samples was set to B = 100.

Scalar g-Sobol function

The scalar g-Sobol function is defined as:

Y = f (X 1 , . . . , X d ) = d i=1 |4X i -2| + a i 1 + a i , a i ∈ R, i ∈ {1, . . . d} where (X 1 , . . . , X d ) ∼ U([0, 1] d ).
For this model, an analytical expression of the scalar Sobol' indices is available:

S i = 1 3(1+a i ) 2 d i=1 3(a i ) 2 +6a i +4 3(1+a i ) 2 -1 , i ∈ {1, . . . , d}.
For our experiments, we chose d = 8 and a = (0, 1, 4. 5,9,99,99,99,99). First-order indices were estimated in N = 100 independent samples of sizes n ∈ {100, 200, 300} using the four methods cnp, np, Tissot2012 and Solis2019. The results can be seen in Figures 4.1 and 4.2.

We remark that the cnp procedure we have proposed in this work behaves well even for small data sets. The bias correction seems to outperform the usual np procedure, except maybe for n = 100.

The variances of both cnp and np are improved with respect to the nonparametric bagging approach in [START_REF] Solís | Non-parametric estimation of the first-order Sobol indices with bootstrap bandwidth[END_REF]. [START_REF] Tissot | Bias correction for the estimation of sensitivity indices based on random balance designs[END_REF] and cnp methods provide good results for samples of size n = 300. Contrary to Tissot2012, the advantage of cnp is that it does not require a structured sample design. At least, we remark a decrease of the efficiency of our approach for very small Sobol' indices.

Multivariate g-Sobol function

We introduce now a multivariate version of the g-Sobol function. It is defined as:

(Y 1 , . . . , Y p ) T = f (X 1 , . . . , X d ) with (X 1 , . . . , X d ) ∼ U([0, 1] d )
where

Y j = d i=1 |4X i -2| + a j i 1 + a j i , a j i , with j ∈ {1, . . . , p}, i ∈ {1, . . . d}.
An analytical expression of aggregated Sobol' indices is available for this test case:

GS i = p j=1 1 3(1+a j i ) 2 p j=1 d i=1 3(a j i ) 2 +6a j i +4 3(1+a j i ) 2
-1 • q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0

0.2 0.4 0.6 0.8 1.0 X 4 n Si a)
q q q q cnp np Solis2019 Tissot2012 100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0

0.2 0.4 0.6 0.8 1.0 X 2 n Si b)
100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0

0.2 0.4 0.6 0.8 1.0 X 5 n Si a)
q q q q cnp np Solis2019 Tissot2012 100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0

0.2 0.4 0.6 0.8 X 6 n Si b)
100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 X 7 n Si c) 100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 We chose p = 2, d = 6 and the coefficients a j i , i ∈ {1, . . . , 6}, j ∈ {1, 2} were coded by the matrix A 2 ∈ R 6×2 :

A T 2 = 0 0.5 3 9 99 99 1 1 1 1 1 1 Figure 4.3 shows the estimation of first-order Sobol' indices using cnp, np, Solis2019 and Tissot2012 in N = 100 independent samples of sizes n ∈ {100, 200, 300}. The bias accuracy of cnp compared to the other three methods is the best (see Figure 4.3). Furthermore, the variance of cnp estimation is lower than the other three methods (see Figure 4.3). Additionally, for this test case, the bias correction performance is more evident than in the previous section (see Figures 4.1

A functional example: the mass-spring model

Before applying our procedure to the avalanche application, we illustrate its performance on a test case with discretized functional output: the functional mass-spring model proposed by Gamboa et al. [2013], where the displacement of a mass connected to a spring is considered:

mx (t) + cx (t) + kx(t) = 0, (4.11)
with initial conditions x(0) = l, x (0) = 0, and t ∈ [1, 40]. There exists an analytical solution of Equation (4.11). This model has four inputs (see more details in Table 4.1). The model output is the vector:

Y = f (X) = (x(t 1 ), . . . , x(t 800 )), t i = 0.05i with i ∈ {1, . . . , 800}.
The discretized output is high-dimensional (p = 800). Therefore, we first reduce the dimension of the output by using PCA before estimating the aggregated Sobol' indices in the first q principal components with q p. The indices were estimated in N = 100 independent samples of size n ∈ {100, 200, 300}. Each sample was also used for PCA. We tested if the results were more accurate if we used two independent samples, one for PCA and the other one for the estimation of Sobol' indices. At a constant cost in terms of model evaluations, the results were even better using only one sample (see Figure 4.12 in the Appendix Section 4.6).

The explained variance as a function of the number of PCs q is shown in Figure 4.4. Notice that 70% of the variance is explained with a small number of PCs (q = 3) and almost all the variance (95%) is explained by 6 PCs. Thus, the aggregated Sobol' indices were estimated using the first 6 PCs. q q q q q q q q q q q q q 0.2 0.4 0.6 0.8 X 1 n GSi a) q q q q cnp np Solis2019 Tissot2012 100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q 0.1

0.3 0.5 0.7 X 2 n GSi b)
100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 X 3 n GSi c) 100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0

0.2 0.4 0.6 0.8 X 4 n GSi d)
100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 X 5

n GSi e) 100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 5 10 15 20 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Basis size Explained variance q q q q q q q q q q q q q q q q q q q q Figure 4.4: Mass-spring model: Explained variance as a function of the decomposition basis size (number of PCs). The gray line is displayed at the 95% of the variance explained which corresponds to 6 PCs.

The true aggregated Sobol' indices for this instance are unknown but they are considered close to the values estimated with a sample of large size n = 100 000 using the pick-and-freeze method of [START_REF] Tissot | A randomized Orthogonal Array-based procedure for the estimation of first-and second-order Sobol' indices[END_REF].

The evolution of the estimation of the aggregated Sobol' indices using the four methods in function of the sample size n ∈ {100, 200, 300} are shown in Figure 4.5. The cnp estimation procedure presents a much smaller variance than the np and Solis2019 ones. The bias correction allows to correct the bias of the np estimation procedure, even if a residual bias seems to persist. The bias decreases as n increases, as expected. As it was mentioned before, our aim is factor prioritization and this objective is achieved with the cnp method in the sense that the most influential factor X 3 is clearly identified.

In general, the cnp method has a better accuracy compared to np, and Solis2019 in all the test cases we considered. Even if Tissot2012 shows better results than cnp for some of the first-order indices, it suffers from two drawbacks: (i) it requires a structured sampling design which makes impossible its application in a given data framework, (ii) it requires the choice of the number of harmonics M , for which no adaptive procedure is available yet (in our experiments M was fixed to 6, as recommended in the literature). For these reasons we apply in the following section the cnp method to estimate the aggregated Sobol' indices of an avalanche model.

Application: the avalanche model

The avalanche model used in this study represents the avalanche motion as a fluid using depth-averaged Saint-Venant equations [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF]. In short, the model considers only the dense layer of the avalanche. The depth of the flow is then small compared to its length. Under this assumption, shallow-water approximations of the mass and momentum equations can be used. Also, the model assumes the avalanche is flowing on a curvilinear q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0

0.1 0.2 0.3 0.4 X 1 n GSi a)
q q q q cnp np Solis2019 Tissot2012 100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0

0.1 0.2 0.3 0.4 0.5 X 2 n GSi b)
100 200 300 q q q q q q q q q q q q q q q q q q q q 0.2 0.3 0.4 0.5 0.6 X 3 n GSi c) 100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 profile z = l(x), where z is the elevation and x is the projected runout length distance measured from the starting abscissa point of the avalanche. The Voellmy friction law is considered. We denote hereafter the avalanche model by f . The avalanche model depends on six inputs: two friction parameters µ and ξ, the length l start of the avalanche release zone, the snow depth h start within the release zone, the release position of the avalanche denoted by x start and the discretized topography of the flow path, denoted by D = (x, z) ∈ R T ×2 where x ∈ R T is the vector of projected runout length from the starting point of the avalanche and z = l(x) ∈ R T is the elevation vector. T is the number of points of the discretized path. In addition to these parameters, we included the term σ to code the error of the digital elevation model z on the path's topography. To do so, a Gaussian error term is added to each element of z to construct a new topography D = (x, z ) for each new model run, whose elevation is calculated as follows:

z = z + N (0, I T σ 2 )
where I T is the T × T identity matrix. This error codes the imprecision of the digital elevation model which allows analysing how it affects the simulation results.

In this study, we analyzed precisely the sensitivity of the model for the simulations of a single avalanche event. The input uncertainty intervals are those corresponding to an avalanche released at the Lautaret full-scale test-site on 13 February 2013. More details about this avalanche event can be found in [START_REF] Pulfer | Retrieving avalanche basal friction law from high rate positioning of avalanches[END_REF][START_REF] María | Bayesian calibration of an avalanche model from autocorrelated measurements along the flow: application to velocities extracted from photogrammetric images[END_REF], thus parameter (x,z) is not uncertain anymore.

For this avalanche event, the release position of the avalanche x start is precisely known because it was fixed by the experimental team (artificial avalanche release). Thanks to the sophisticated equipment of the test-site, the uncertainty of this input can therefore be neglected. The l start and h start input uncertainty intervals could be determined thanks to measures taken during the experiment. By contrast, the inputs µ and ξ cannot be measured. In another study, we developed a Bayesian calibration approach to infer these inputs from the measurements made along the flow. Thus, the uncertainty intervals of µ and ξ considered here correspond to the 95% credibility intervals of the posterior distribution we obtained [START_REF] María | Bayesian calibration of an avalanche model from autocorrelated measurements along the flow: application to velocities extracted from photogrammetric images[END_REF]. However, to avoid privileging some values in the sensitivity analysis, uniform distributions were chosen to model the d = 5 inputs, including the friction parameters (i.e. we did not directly use their posterior distribution). The uncertain inputs of the model and their uncertainty intervals are summarized in Table 4.2.

The outputs of the avalanche model are the functional maximal velocity v and the functional maximal flow depth h of the avalanche on the discretized grid corresponding to

Inputs

Model Outputs the topography D, and the runout distance of the avalanche denoted by x runout . In other words, the model has two vectorial and one scalar outputs.

µ ξ σ 2 topo l start h start f Y =    v : D ⊂ R → R h : D ⊂ R → R x runout ∈ R +
The vectors whose components are the evaluations on the discretized topography D of the velocity and flow depth functional outputs are denoted v = (v 1 , . . . , v T ) and h = (h 1 , . . . , h T ), respectively. Note that we have p = 2 × T + 1 outputs. A brief summary of our framework is drawn in Figure 4.6.

The global sensitivity analysis is conducted only on a section of the path where all the avalanche simulations are flowing. More precisely, considering e.g., the velocity, our sensitivity analysis is performed on a subset of length n, v = {v 1 , . . . , v T 1 } ∈{1,...n} , with T 1 ≤ T and such that v k > 0 for all k ∈ {1, . . . , T 1 }.

Sobol' indices

To estimate the sensitivity indices, we developed n = 300 model simulations using the random input distributions shown in Table 4.2. Figure 4.8 shows the topography of the avalanche path and the gray box shows the region where all the avalanche simulations are flowing. The output values in this subset are used in the sensitivity analysis. Also, in Figure 4.8, we show the percentage of avalanches flowing in the path as function of the position along the topography.

Scalar Sobol' indices

The scalar Sobol' indices calculated using n = 300 samples are shown in Figure 4.9. For the velocity output, in the first 50 m of the path, the input µ is the most important, then it is followed by ξ. For the rest of the path, the input ξ is the most important. For the flow depth output, the h start input is the most important. The input σ 2 is the least important parameter for both outputs. The scalar sensitivity indices for the runout distance are shown in Figure 4.10. For this output, the h start input is the most important. These three figures give us valuable information about the model. However, it could also be interesting to get a summarized information by computing aggregated sensitivity indices and this is done in the following section.

Aggregated Sobol' indices

To estimate the aggregated indices, we used two ways of aggregating the output information. In the first one, we performed two PCAs separately: one for the velocity output and the other one for the flow depth output. And in the second one, we used simultaneous PCA to estimate the indices of the three outputs. The results of the aggregated indices using two PCAs are shown in the panel a) of Figure 4.11 and using simultaneous PCA in the panel b) of Figure 4.11. This summarizes the importance of the inputs.

Let us analyze the results of the aggregated indices calculated when performing a PCA on each functional output: for the velocity output, the ξ and µ inputs are the most relevant and for the flow depth output, the h start input is by far the most important. These results are consistent with findings of other studies [e.g., [START_REF] Barbolini | Estimate of uncertainties in avalanche hazard mapping[END_REF][START_REF] Borstad | Sensitivity analyses in snow avalanche dynamics modeling and implications when modeling extreme events[END_REF][START_REF] Fischer | Multivariate parameter optimization for computational snow avalanche simulation[END_REF].

Finally, let us analyze the results of the aggregated indices calculated when performing a simultaneous PCA for the whole set of outputs (the two functional ones and the scalar runout distance): the most important inputs are ξ, h start and µ.

Conclusions and perspectives

In this work, we proposed a nonparametric method to estimate the aggregated Sobol' indices from a given random sample of small size, and called this method cnp. The method is based on the Nadaraya-Watson kernel smoother. Due to the small size of the sample at hand, the kernel estimation is biased. Therefore, to remove the bias of the estimation, we proposed a bias correction using bootstrapping samples based on the works of [START_REF] Racine | Bias-Corrected Kernel Regression[END_REF] and [START_REF] Solís | Non-parametric estimation of the first-order Sobol indices with bootstrap bandwidth[END_REF] before bandwidth selection based on cross-validation.

We tested the accuracy of the method on a scalar and two multivariate test cases. In general, the method cnp is accurate, even if the sample size is low (n = 300). The method was developed to estimate the aggregated indices for an avalanche model which has three outputs: the functional velocity and flow depth discretized on a given topography and 4.6. Appendix the runout distance of the avalanche. In this work, we developed a sensitivity analysis of a single avalanche event. The event corresponds to an avalanche released at the Lautaret test-site the 13 February 2013. The results have shown that for this particular avalanche event: for the velocity output, the µ and ξ inputs are the most relevant and for the snow depth output, h start is the most important. Finally, for the runout distance output, h start is the most important. We also showed how this information can be aggregated in one single set of indices summing up the model sensitivity which could be a very useful information for avalanche practitioners.

Note that in this work we proposed a method to quantify the importance of the inputs of a particular avalanche model but this method can be widely applied to other avalanche models. Moreover, the method could be applied to other avalanche events with the same model to generalize the results. Eventually, it could be adapted to various problems for which complex models with outputs which are both functional and scalar are employed.

Eventually, in this paper, we did not provide confidence intervals (CI) associated to the estimation of Sobol' or aggregated Sobol' indices. This would be an interesting task for a future work. However, bootstrap based CI are costly in that framework as they would involve two boostrap stages: one for the bias correction and one for the computation of CI.

Appendix

In this appendix, we compare the accuracy of Sobol' index estimation for the mass-spring test case when using a single sample of size n = 200 for both PCA and Sobol' index estimation or two independent samples of size n = 100, one for PCA, the other one for Sobol' index estimation. Boxplots presented in Figure 4.12 were obtained with N = 100 independent replications. We also compare the estimation of the indices using an Epanechnikov of second order kernel (left panel) or a Gaussian kernel (right panel) in the g-Sobol test case. The results are shown in Figure 4.14. q q q q q q q 0.00 0.05 0.10 0.15 0.20 X 1 GSI a) q q q cnp np Tissot2012 one sample two samples q q q q q q q q q q q q q q q q q q q q q q 0.00 0.05 0.10

0.15 X 2 GSI b)
one sample two samples q q q q q 0.15 0.20 0.25 0.30 X 3

GSI c)

one sample two samples q q q q q q q q q 0.00 Abstract: Dynamic models are simplified representations of some real-world entities that change over time. They are essential analytical tools with significant applications, e.g., in environmental and social sciences. The outputs produced by dynamic models are typically time and/or space dependent. Due to physical constraints, their parameters are confined to a non-rectangular domain. Also, they can be significantly sensitive to variations of input parameters. A global sensitivity analysis (GSA) consists in modeling input parameters by a probability distribution which propagates through the model to the outputs. Then, input parameters are ordered according to their contribution on the model outputs by computing sensitivity measures. In this paper, we extend Shapley effects, a sensitivity measure well suited for dependent input parameters, to the framework of dynamic models. Note that if input parameters are constrained to a non-rectangular domain, they cannot be modeled by independent random variables. We also propose an algorithm to estimate the so-called aggregated Shapley effects and to construct bootstrap confidence intervals for the estimation of scalar and aggregated Shapley effects. We measure the performances of the estimation procedure and the accuracy of the probability of coverage of the bootstrap confidence intervals on toy models. Finally, our procedure is applied to perform a GSA of an avalanche flow dynamic model, for which the sample is obtained by restricting the initial input/output sample to keep only combinations that are physically meaningful. More precisely, we analyze the sensitivity in two different settings: (i) little knowledge on the input parameter probability distribution, and (ii) well-calibrated input parameter distributions. Probative linkages between local slope and sensitivity indices demonstrate the usefulness of our approach for practical problems.

Introduction

Dynamic models are simplified representations of some real-world entity that change over time, in equations or computer code. These models are useful for the analysis of realworld phenomena, e.g., in environmental or social sciences. For a better understanding of a phenomenon or for forecasting purposes, it might be important to identify which input parameters entering in the formulation of such dynamic models are influential on the outputs of interest. Determining these influential parameters is one aim of global sensitivity analysis (GSA). A global sensitivity analysis (GSA) consists in modeling unknown input parameters by a probability distribution which propagates through the model to the outputs. Then, input parameters are ordered according to their contribution on the model outputs by computing sensitivity measures. In the literature, there exist different global sensitivity measures, e.g., variance based measures such as Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] or Shapley effects [START_REF] Owen | Sobol' Indices and Shapley Value[END_REF], density based measures [Borgonovo et al., 2016[START_REF] Da | Global sensitivity analysis with dependence measures[END_REF], etc. A review of global sensitivity measures can be found in, e.g., Borgonovo and Plischke [2016] or [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF].

A particularity of dynamic models considered in this paper is that the outputs they produce are typically time and/or space dependent (see e.g., [START_REF] Alexanderian | Variance-based sensitivity analysis for time-dependent processes[END_REF], [START_REF] Lamboni | Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models[END_REF]). More specifically, the application that motivated our study is an avalanche flow dynamic model (proposed by [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF]). The input parameters in our application cannot be considered as stochastically independent. Indeed, even if they are first sampled independently, we only analyze in the GSA the ones leading to snow avalanche simulations that are both realistic and of interest for risk purposes. For example, an avalanche simulation is rejected if its volume is inferior to 7000 m 3 as such small events do not represent a significant threat for settlements downslope. It is quite common that such physical constraints lead to dependent inputs by excluding infeasible sets of input values. Let us cite, e.g., [START_REF] Majdi | Shapley effect application for variance-based sensitivity analysis of the few-group cross-sections[END_REF] dealing with the response of a nuclear reactor or [START_REF] López | A case study on global sensitivity analysis with dependent inputs: The natural gas transmission model[END_REF] interested in a natural gas transmission model (see also [START_REF] Kucherenko | Applications of Probabilistic Graphical Models in Dependability, Diagnosis and Prognosis[END_REF]). For that reason, we develop a GSA methodology capable to handle dependent inputs, with unknown dependence structure, as far as dynamic models through functional outputs (or high-dimensional multivariate outputs as functional outputs are most of the time represented on a discrete grid).

Although the independence assumption on input parameters is unrealistic in many applications, it is traditionally required to interpret or to compute sensitivity measures.

Many alternatives to Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] based on functional ANOVA decomposition have been proposed in the literature to handle GSA in the framework of dependent input parameters. We refer the interested reader to the introduction in [START_REF] Iooss | Shapley effects for sensitivity analysis with correlated inputs: Comparisons with Sobol' indices, numerical estimation and applications[END_REF] and to references therein for a complete overview of these alternatives.

In this paper, we focus on Shapley effects introduced in [START_REF] Owen | Sobol' Indices and Shapley Value[END_REF], consisting in variance based sensitivity measures, still meaningful in the framework of dependent input parameters [START_REF] Owen | On Shapley value for measuring importance of dependent inputs[END_REF]. This approach is based on the Shapley value which is a cooperative game theory concept. Briefly speaking, Shapley value ensures a fair distribution of a gain among team players according to their individual contributions. As a sensitivity measure, [START_REF] Owen | Sobol' Indices and Shapley Value[END_REF] adapted the Shapley value into the Shapley effects by considering model input parameters as players and the gain function as the output variance. The main advantage of such an approach is that it is possible to attribute a non negative sensitivity index to each parameter, and the sum of the indices is equal to one [Broto et al., 2020, Iooss andPrieur, 2019].

A first algorithm to estimate Shapley effects was proposed in [START_REF] Song | Shapley effects for global sensitivity analysis: Theory and computation[END_REF]. A block bootstrap procedure was proposed in [START_REF] Benoumechiara | Shapley effects for sensitivity analysis with dependent inputs: bootstrap and kriging-based algorithms[END_REF] to add the construction of confidence intervals to the aforementioned algorithm. More recently, [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF] proposed a subset aggregation procedure which leads to a significant reduction of the variance of Shapley effect estimation. Let us also cite the algorithm proposed in [START_REF] Plischke | Computing Shapley Effects for Sensitivity Analysis[END_REF] based on the Möbious inverse, which offers a computationally efficient alternative for the estimation of Shapley effects. The main advantage of the subset aggregation procedure introduced in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF] is that it has a version, based on nearest-neighbors, which does not require the ability to sample from the exact conditional distributions of the input parameters.

In the present paper, we propose an extension of Shapley effects to multivariate or functional outputs. Although it is possible to compute a Shapley effect for each component of a multivariate output, it leads to results difficult to interpret as the balance of sensitivities may change from one component to the other [START_REF] Alexanderian | Variance-based sensitivity analysis for time-dependent processes[END_REF] and to redundancies in the information, particularly when considering the discretization of a time and/or space dependent output [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF]. For that reason, [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] (see also Gamboa et al. [2013]) extended Sobol' indices to multivariate or functional outputs. Let us also cite [START_REF] Alexanderian | Variance-based sensitivity analysis for time-dependent processes[END_REF] focused on time-dependent outputs. However, the so-called aggregated Sobol' indices cannot easily be interpreted if the input parameters are dependent. Following these papers, we introduce and study the properties of what we call aggregated Shapley effects, which can be interpreted in both the framework of independent and dependent input parameters. If the output dimension is high (as it is the case, e.g., when considering a fine discretization of a functional output), a dimension reduction can be applied as a preliminary step to the computation of aggregated Shapley effects. We use the Karhunen-Loève (KL) expansion as in [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF][START_REF] Alexanderian | Variance-based sensitivity analysis for time-dependent processes[END_REF]. More precisely to perform KL expansion, we use the functional principal component analysis proposed by [START_REF] Yao | Functional data analysis for sparse longitudinal data[END_REF]. Note that the extension of Shapley effects to multivariate outputs has been early used for nuclear safety application in [START_REF] Delipei | Développement d'une méthodologie de Quantification d'Incertitudes pour une analyse Multi-Physique Best Estimate et application sur un Accident d' Éjection de Grappe dans un Réacteur à Eau Pressurisée[END_REF]. In this paper, we propose a proper definition of aggregated Shapley indices, we study their invariance properties and we propose a given data estimation procedure, extending the nearest neighbor approach introduced in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF] and a block bootstrap procedure to build confidence intervals. By given data procedure, we mean that a single sample of size n is needed for the estimation of all Shapley effects. Eventually, we demonstrate the usefulness of this methodology in a complex avalanche application where the set of output simulations is constrained in order to keep only meaningful snow avalanche events. More precisely, we perform the GSA of our avalanche flow dynamic model in two different settings: (i) little knowledge on the input parameter probability distribution, and (ii) well-calibrated input parameter distribution [START_REF] Eckert | Long-term avalanche hazard assessment with a Bayesian depth averaged propagation model[END_REF]). Notably, this methodology which could be applied to many other applications where samples are constrained to satisfy physical or user constraints opens a new field of applications to GSA.

The paper is organized as follows. In Section 5.2, aggregated Shapley effects and their main properties are described. In Section 5.3, we propose an estimator for aggregated Shapley effects in a given data framework by extending the Monte-Carlo nearest-neighbor estimator of scalar Shapley effects introduced in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF]. At the end of the section, we describe the functional principal components analysis algorithm to perform model dimension reduction proposed by [START_REF] Yao | Functional data analysis for sparse longitudinal data[END_REF]. In Section 5.4, we propose a bootstrap algorithm to construct confidence intervals of the scalar and aggregated Shapley effect estimations based on [START_REF] Benoumechiara | Shapley effects for sensitivity analysis with dependent inputs: bootstrap and kriging-based algorithms[END_REF]. In Section 5.5, we test our estimation procedure on two toy models: a multivariate linear Gaussian model and the mass-spring model. Finally in Section 5.6, our GSA procedure is applied to our avalanche dynamic model. We discuss our conclusions and perspectives in Section 5.7.

Aggregated Shapley effects

Shapley effects are sensitivity measures to quantify input importance proposed by [START_REF] Owen | Sobol' Indices and Shapley Value[END_REF]. These measures are particularly useful when inputs are dependent. Shapley effects are based in the concept of Shapley value, introduced in the framework of game theory [START_REF] Shapley | A Value for n-Person Games[END_REF], which consists into dividing a game gain among a group of players in an equitable way. As sensitivity measures, Shapley effects consider model inputs as players and output variance as game function. Shapley effects can be naturally extended to multivariate output models by following the ideas presented in [Gamboa et al., 2013] and [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] to generalize Sobol' indices to multivariate output models (see also [START_REF] Alexanderian | Variance-based sensitivity analysis for time-dependent processes[END_REF] for time-dependent models). We call these new sensitivity measures aggregated Shapley effects.

Definition

Let us define Y = (Y 1 , . . . , Y j , . . . , Y p ) = f (X) the p multivariate output of a model f that depends on d random inputs X = (X 1 , . . . , X d ). The inputs are defined on some probability space (Ω, F, P X ) and f ∈ L 2 (P X ). For any u ⊆ {1, . . . , d}, let us define -u = {1, . . . , d} \ u its complement. We set X u = (X i ) i∈u . Note that the inputs are not necessary independent. In the framework of our application to avalanche modeling, the model produces outputs of the form Y = (Y 1 = f (s 1 , X), . . . , Y p = f (s p , X)), with s 1 , . . . , s p ∈ R the p discretization points along the avalanche corridor.

In this section we recall the definition and main properties of the Shapley value, on which the definition of Shapley effects is based. Given a set of d players in a coalitional game and a charateristic function val : 2 d → R, val(∅) = 0, the Shapley value (φ 1 , . . . , φ d ) is the only distribution of the total gains val({1, . . . , d}) to the players satisfying the It is proved in [START_REF] Shapley | A Value for n-Person Games[END_REF] that according to the Shapley value, the amount that player i gets given a coalitional game (val, d) is:

φ i = 1 d u⊆-{i} d -1 |u| -1
(val(u ∪ {i})val(u)) ∀i ∈ {1, . . . , d}.

(5.1)

The Shapley value also satisfies the linearity property:

5. (Linearity) Let λ ∈ R, if λval and val have Shapley values φ and φ, then φ i = λφ i for all i ∈ {1, . . . , d}.

The linearity property is a direct consequence from ( 5.1) and it is used to prove some of the nice properties of aggregated Shapley effects (see Propositions 1 and 2 further).

The Shapley effects are defined by considering the characteristic function of the game as: 5.2) in (5.1). Thus, the scalar Shapley effect of input i in output j is defined as:

val j (u) = Var(E(Y j |X u )) Var(Y j ) , u ⊆ {1, . . . , d} (
Sh j i = 1 d Var(Y j ) u⊆-{i} d -1 |u| -1 (Var(E(Y j |X u∪i )) -Var(E(Y j |X u ))).
(5.3) Shapley effects can be naturally extended to models with multivariate outputs following ideas from [Gamboa et al., 2013] and [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] where authors proposed to extend Sobol' indices to multivariate outputs. Aggregated Shapley effect of an input i is then defined as: .4) where Sh j i is the scalar Shapley effect of input X i in output Y j . This sensitivity measure is a weighted sum of the scalar Shapley effects where weights correspond to the proportion of the variance of each output over the sum of all individual variances.

GSh i = p j=1 Var(Y j )Sh j i p j=1 Var(Y j ) , ( 5

Properties

In this section, we prove some nice properties of aggregated Shapley effects.

val(u) = p j=1 Var(Y j )val j (u) p j=1 Var(Y j ) = p j=1 Var(E(Y j |X u )) p j=1 Var(Y j ) = tr(Σ u ) tr(Σ) (5.6)
where Σ u is the covariance matrix of E(Y|X u ) and Σ is the covariance matrix of Y. Note that the characteristic function val of aggregated Shapley effects corresponds to the definition of the aggregated Sobol' indices introduced in [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF], Gamboa et al., 2013]. In the next proposition, we prove that aggregated Shapley effects accomplish the natural requirements for a sensitivity measure mentioned in Proposition 3.1 in [Gamboa et al., 2013].

Proposition 2. Let i ∈ {1, . . . d}. The following items hold true.

i. 0 ≤ GSh i ≤ 1.
ii. GSh i is invariant by left-composition by any nonzero scaling of f , which means, for any λ ∈ R, the aggregated Shapley effect GSh i of λf (X) is GSh i .

iii. GSh i is invariant by left-composition of f by any isometry of R p , which means, for any O ∈ R p×p such that O t O = I, the aggregated Shapley effect GSh i of Of (X) is GSh i for all i ∈ {1, . . . , d}.

Proof. i. As for all j ∈ {1, . . . , p} 0 ≤ Sh j i ≤ 1 and as the sum of the non negative weights Var(Y j )/ p =1 Var(Y ) is one, we deduce that 0 ≤ GSh i ≤ 1. ii. Note that GSh i can be written as GSh i = p j=1 Var(λY j )Sh j i / p j=1 Var(λY j ), where Sh j i is the Shapley effect associated to the characteristic function val j . Note that val j (u) = Var(E(λY j |X u ))/Var(λY j ) = val j (u) for all u ⊆ {1, . . . , d}. Thus, Sh j i = Sh j i from where GSh i = GSh i which means the aggregated Shapley effect is invariant by any nonzero scaling of f . iii. Let us write g(X) = Of (X) = OY = U. The characteristic function associated to the aggregated Shapley effect GSh i of U is then (see Equation (5.6)) val'(u) = tr(Σ U u )/tr(Σ U ) where Σ U u is the covariance matrix of E(U|X u ) and Σ U is the covariance matrix of U. Then,

val'(u) = tr(Σ U u ) tr(Σ U ) = tr(OΣ Y u O t ) tr(OΣ Y O t ) = tr(Σ Y u ) tr(Σ Y )
= val(u) for all u ⊆ {1, . . . , d}.

As val(u) has a unique Shapley value GSh i , val (u) has Shapley value GSh i which proves that GSh i = GSh i for all i ∈ {1, . . . , d}.

In this section, we have proven that aggregated Shapley effects are sensitivity measures. In the next section, we describe the estimation procedure we propose for aggregated Shapley effects, based the estimation procedure of scalar Shapley effects proposed in [Broto et al., 2020, Section 6] when observing an i.i.d. sample of (X, Y). Such a procedure, which does not require a specific form for the design of experiments is also called given data procedure.

Estimation procedure for scalar and aggregated Shapley effects

The aggregated Shapley effect estimation procedure we propose in this section is based on the given data estimation procedure of the scalar Shapley effects introduced in [Broto et al., 2020, Section 6.1.1.]. In the application we consider in Section 5.6, the standard pick-freeze estimation procedure (see, e.g., Janon et al. [2014b]) cannot be used as it is based on a specific pick-freeze design of experiments. We recall indeed that input parameters in our application are first sampled independently, then some combinations are rejected as they do not lead to snow avalanche avalanche simulations that are both physically meaningful and of interest for risk purposes. It is the reason why we focus our attention to the given data estimation procedure of scalar Shapley effects introduced in [Broto et al., 2020, Section 6.1.1.]. For sake of clarity, we first present the estimation procedure for scalar Shapley effects in Subsection 5.3.1 before extending it to the estimation of aggregated Shapley effects in Subsection 5.3.2.

5.

3.1 Given data estimator of scalar Shapley effects proposed by [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF] Let u ⊂ {1, . . . d}. Let -u denote the complementary set of u in {1, . . . , d}. As noticed in [Song et al., 2016, Theorem 1], replacing the characteristic function

u → cj (u) = Var(E(Y j |X u )) by u → c j (u) = E(Var(Y j |X -u
)) in Equation ( 5.3) does not change the definition of Shapley effects. The estimation of the Shapley effects Sh j i , 1 ≤ i ≤ d, relies on the estimation of the c j (u) (resp. cj (u)) for all u ⊆ {1, . . . , d}. To estimate c j (u) = E(Var(Y j |X -u )) (resp. cj (u)), one can use a double Monte Carlo estimator: an inner loop to estimate the variance (resp. the expectation) and an outer loop to estimate the expectation (resp. the variance). However, [START_REF] Song | Shapley effects for global sensitivity analysis: Theory and computation[END_REF] pointed out that the double Monte Carlo estimator of cj (u) can suffer from a non neglectable bias if the inner loop sample is small, while in contrast the double Monte Carlo estimator of c j (u) is unbiased for any sample size. For that reason, we turned to the double Monte Carlo estimator of scalar Shapley effects introduced in [Broto et al., 2020, Section 6.1.1.] based on the double Monte Carlo estimation of c j (u) with a subset aggregation procedure. They also introduce an adaptation of this estimator to the given data framework, well suited to our application as it is based on a given and unstructured sample.

The estimator works as follows: assume we are given a n independent random sample (X (k) , Y (k) ) 1≤k≤n of (X, Y), then for 1 ≤ j ≤ p, ∅ u {1, . . . , d}, let 1 ≤ N u ≤ n be an integer, let (s ) 1≤ ≤Nu be a sample of uniformly distributed integers in {1, . . . , n} (with or without replacement). The inner Monte Carlo loop estimates the conditional variance Var(Y j |X -u ) using N I nearest neighbors of each X s -u , 1 ≤ ≤ N u . If all the N u were chosen equal to n, the computation cost for the search of the nearest neighbors would be in O n(2 d -2) . The integers (N u ) ∅ u {1,...,d} are related to the estimation accuracy of the conditional variances (Var(Y j |X -u )) ∅ u {1,...,d} and are referred as the accuracies hereafter.

To reduce the estimation cost to a prescribed computational budget N tot , the authors in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF] proposed to choose the set of accuracies that minimize the estimation variance for the prescribed total cost N tot = ∅ u {1,...,d} N u . In the framework where it is possible to sample from the conditional distributions of the input vector, it can be proven under a strong assumption (see [Broto et al., 2020, Proposition 4.2., Remark 4.3.]) that there exists an optimal allocation of the accuracies, namely

N * u = N tot d |u| -1 (d -1) -1 ,
for all ∅ u {1, . . . , d}. Even if this strong assumption, that all the conditional elements E(Var(Y |X -u )), u ⊆ {1, . . . , d}, u = ∅, u = {1, . . . , d} can be estimated with the same variance, does not hold true in the given data framework, the authors present numerical experiments in [Broto et al., 2020, Section 6.3], showing that if one chooses the accuracies N * u for all ∅ u {1, . . . , d}, one gets reasonable results, at least on the test case they consider. Note that the estimation cost in terms of number of model evaluations is n while the cost in terms of nearest-neighbor search is N tot . In Section 5.8, we compare the choice of the accuracy N * u for all ∅ u {1, . . . , d} to a uniform allocation N u = N tot /(2 d -2) for all ∅ u {1, . . . , d} for a test case in dimension d = 10. The allocation N * u suggested in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF] leads to better results, at least on that experiment. Regarding the number of neighbors, we analyzed the impact of using different numbers of neighbors N I on a linear model with d = 2. The results are presented in Section 5.9. We see in Figure 5.14 that there is no clear evidence if N I ≥ 3 is always the best choice. However in the following we choose N I = 3 following the recommendations in [START_REF] Song | Shapley effects for global sensitivity analysis: Theory and computation[END_REF], where the authors prove in the non given data framework that choosing N I = 3 is better than choosing N I = 2 and that the gain from increasing N I from 3 decreases as N I increases.

Let us write more precisely the inner Monte Carlo loop to estimate the conditional variance Var(Y j |X -u ):

E j u,s = 1 N I -1 N I i=1   υ:X υ -u ∈B -u, Y υ j - 1 N I Ȳs   2 with Ȳs = 1 N I υ:X υ -u ∈B -u, Y υ j , (5.7)
with B -u, the set of N I closest neighbors (chosen two by two distinct) of X s -u . The closest neighbors are computed using the Euclidean distance in R ku with k u denoting the cardinal of -u. Then, we compute the outer Monte Carlo loop to estimate the expectation:

c j (u) = 1 N * u N * u =1 E j u,s .
(5.8)

Note that c j (∅) = 0 and that c j ({1, . . . , d}) = Var(Y j ), which is assumed to be known in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF], and that is estimated by the empirical variance in the present paper. In [Broto et al., 2020, Theorem 6.6.], it is proved that under mild assumptions, c j (u) converges in probability to c j (u) as n and N * u go to ∞. Finally, the subset aggregation procedure consists in estimating scalar Shapley effects by plugging (5.8) in Equation ( 5.3):

Sh j i = 1 d σ 2 j u⊆-i d -1 |u| -1 ( c j (u ∪ {i}) -c j (u)) (5.9)
where σ2 j is the empirical estimator of Var(Y j ).

Estimator of the aggregated Shapley effects

Given scalar Shapley effect estimators whose definition is recalled in the previous section, we propose to estimate the aggregated Shapley effects by:

GSh i = p j=1 σ 2 j Sh j i p j=1 σ 2 j = 1 d p j=1 σ 2 j p j=1 u⊆-i d -1 |u| -1 ( c j (u ∪ {i}) -c j (u)),
(5.10) with σ2 j the empirical estimator of Var(Y j ) and c j (u) defined by (5.8).

Dimension reduction: functional principal component analysis

If model f is space or time-dependent, inspired by [START_REF] Alexanderian | Variance-based sensitivity analysis for time-dependent processes[END_REF] and [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF], we perform a functional principal component analysis (fPCA) based on the Karhunen-Loève (KL) expansion to obtain a low-rank model representation. Then we compute the aggregated Shapley effects for the vectorial output composed with the first principal components. This parsimonious representation of the functional output usually allows easier interpretation for real applications. To perform KL expansion, we use the principal component analysis through conditional expectation (PACE) method introduced in [START_REF] Yao | Functional data analysis for sparse longitudinal data[END_REF]. More precisely, we have a collection of n independent trajectories of a smooth random function f (., X) with unknown mean µ(s) = E(f (s, X)), s ∈ τ , where τ ⊆ R is a bounded and closed interval, and covariance function G(s 1 , s 2 ) = Cov(f (s 1 , X), f (s 2 , X)), s 1 , s 2 ∈ τ . We assume that G has a L 2 orthogonal expansion in terms of eigenfunctions ξ k and non increasing eigenvalues λ k such that:

G(s 1 , s 2 ) = k≥1 λ k ξ k (s 1 , X)ξ k (s 2 , X), s 1 , s 2 ∈ τ.
The KL orthogonal expansion of f (s, X) is: (5.11) where α k (X) = τ f (s, X)ξ k (s)ds is the k-th functional principal component (fPC) and q is a truncation level. Following the authors in [START_REF] Yao | Functional data analysis for sparse longitudinal data[END_REF], we estimate first μ(s) using local linear smoothers and G(s 1 , s 2 ) using local linear surface smoothers [START_REF] Fan | Local polynomial modelling and its applications[END_REF]), then we estimate eigenfunctions and eigenvalues by solving the following integral equations:

f (s, X) = µ(s) + k≥1 α k (X)ξ k (s) ≈ µ(s) + q k=1 α k (X)ξ k (s), s ∈ τ,
τ G(s 1 , s) ξ k (s 1 )ds 1 = λ k ξ k (s), s ∈ τ,
with τ ξ(s)ds = 1 and τ ξ k (s) ξ m (s) = 0 for all m = k ≤ q. The problem is solved by using a discretization of the smoothed covariance (see further details in [START_REF] Rice | Estimating the mean and covariance structure nonparametrically when the data are curves[END_REF] and [START_REF] Capra | An accelerated-time model for response curves[END_REF]). Finally, fPCs α k (X) = τ f (s, X) ξ k (s)ds are solved by numerical integration. Aggregated Shapley effects are finally estimated using the estimated low rank KL model representation with truncation level q, in other words, they are computed with only the q first estimated fPCs:

GSh i = 1 d q k=1 λ k q k=1 u⊆-i d -1 |u| -1 E(Var( α k (X)|X u∪{i} )) -E(Var( α k (X)|X u )) .
(5.12)

Remark. (5.12) can be estimated using (5.10).

Remark. To avoid splitting our sample, we use the same sample to estimate the KL representation and to estimate aggregated Shapley effects. The comparison with the alternative procedure consisting in splitting the initial sample in two subsamples, one to estimate the KL representation and the other one to estimate aggregated Shapley effects is presented and discussed for the mass-spring model in Section 5.10. Based on these numerical results, we decided to apply the procedure without splitting in the following and in particular on our snow avalanche application.

Bootstrap confidence intervals with percentile bias correction

Confidence intervals are a valuable tool to quantify uncertainty in estimation. We consider non parametric bootstrap confidence intervals with bias percentile correction (see, e.g., [START_REF] Efron | Nonparametric standard errors and confidence intervals[END_REF], [START_REF] Efron | Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy[END_REF]). More precisely, we propose to construct confidence intervals, with a block bootstrap procedure, following ideas in [START_REF] Benoumechiara | Shapley effects for sensitivity analysis with dependent inputs: bootstrap and kriging-based algorithms[END_REF]. Indeed, bootstrap by blocks is necessary to preserve the nearest-neighbor structure in Equation ( 5.7) and to avoid potential equalities in distance (see Assumption 6.3 in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF]). We describe in Algorithm 1 how to create B bootstrap samples for scalar Shapley effect estimators Sh j i and aggregated Shapley effect estimators GSh i , and then we describe the percentile bias correction method.

If model output is scalar, only Steps 1 to 3 of Algorithm 1 should be used. The block bootstrap procedure is described by Steps 3.1 to 3.3. Also, the same sample (X, Y) is used to estimate the variance of the outputs Y j , 1 ≤ j ≤ p, and the Shapley effects.

For 1 ≤ i ≤ d, 1 ≤ j ≤ p, let R i = { GSh (1) i , . . . , GSh (B) i } and R j i = { Sh j,(1) i , . . . , Sh j,(B) i
}, the bias-corrected percentile method presented in [START_REF] Efron | Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy[END_REF]] is applied. Let us denote by Φ the standard normal cumulative distribution function and by Φ -1 its inverse. A bias correction constant z 0 , estimated as 5.4. Bootstrap confidence intervals with percentile bias correction Algorithm 1 B bootstrap samples for Sh j i and GSh i Inputs:

(i) A n i.i.d. random sample (X k , Y k ) k∈{1,...,n} with X k ∈ R d and Y k ∈ R p .
(ii) For each ∅ u {1, . . . , d}, a N u random sample (s ) 1≤ ≤Nu from {1, . . . , n}.

Outputs: B bootstrap samples for Sh

j u⊆-i d -1 |u| -1 c j (b) (u ∪ {i}) -c j (b) (u) .
4. Compute the b bootstrap sample of GSh i using (5.10):

GSh (b) i = 1 d p j=1 σ 2,(b) j p j=1 u⊆-i d -1 |u| -1 c j (b) (u ∪ {i}) -c j (b) (u) .
end for

ẑ0 = Φ -1 (#{ GSh (b) i ∈ R i s. t. GSh (b) i ≤ GSh i }/B) is computed (similar for Sh j i ).
Then, the corrected quantile estimate q(β) for β ∈]0, 1[ is defined as qi (β) = Φ(2ẑ 0 + z β ), where z β satisfies Φ(z β ) = β. Corrected bootstrap confidence interval of level 1α is estimated by the interval whose endpoints are qi (α/2) and qi (1α/2).

To guarantee the validity of the previous BC corrected confidence interval, there must exist an increasing transformation g, z 0 ∈ R and τ > 0 such that g( GSh i ) ∼ N (GSh iτ z 0 , τ 2 ) and g( GSh * i ) ∼ N ( GSh iτ z 0 , τ 2 ) where GSh * i is the bootstrapped GSh i for fixed sample (see [START_REF] Efron | Nonparametric standard errors and confidence intervals[END_REF]). Normality hypothesis can be tested using traditional normality tests as Shapiro test or using graphical methods as empirical normal quantile-quantile plots. In our application and test cases, we observed that g can be chosen as the identity. The author in [START_REF] Efron | Better bootstrap confidence intervals[END_REF] advises the use of the bias correction to improve the accuracy of bootstrap intervals. In our application and test cases, there was no much difference between the percentile and BC intervals. However, as a good practice we performed the bias correction, especially since the algorithm is simple to implement.

To validate the coverage accuracy of the bootstrap intervals whose construction is described in Algorithm 1, we compute the empirical probability of coverage (POC) of simultaneous intervals using Bonferroni correction for the test cases where the true value of the effects is known. The POC with Bonferroni correction is the probability that the interval

[ Ĝ-1 • qi (α/(2d)), Ĝ-1 • qi (1 -α/(2d))],
where Ĝ is the cumulative distribution function for GSh i , contains GSh i for all i ∈ {1, . . . , d} simultaneously. To be more precise, if the confidence intervals are computed in N independent samples of size n of (X, Y), the POC is estimated as 2d)) for all i, and 0 otherwise.

P OC = N k=1 w k /N, where w k is equal to 1 if Ĝ-1 • qi (α/(2d)) ≤ GSh i ≤ Ĝ-1 • qi (1 -α/(

Test cases

In this section, we numerically study the performance of the estimation procedure and the probability coverage of the boostrap confidence intervals we introduced in the previous section. We consider two test cases: a multivariate linear Gaussian model and the functional mass-spring model proposed in the work of Gamboa et al. [2013]. To estimate the scalar Shapley effects, we use the function shapleySubsetMc of the R package sensitivity corresponding to the estimation procedure defined by (5.8), (5.7) and (5.9). Functional PCA is performed using the R package FPCA [START_REF] Chen | fdapace: Functional Data Analysis and Empirical Dynamics[END_REF].

Multivariate linear Gaussian model in dimension d = 2

We consider a multivariate linear model with two Gaussian inputs based on the example from [START_REF] Iooss | Shapley effects for sensitivity analysis with correlated inputs: Comparisons with Sobol' indices, numerical estimation and applications[END_REF]. For this toy function, there is an analytical expression of the scalar and aggregated Shapley effects (see [START_REF] Iooss | Shapley effects for sensitivity analysis with correlated inputs: Comparisons with Sobol' indices, numerical estimation and applications[END_REF]).

The model f is defined as Y = f (X) = B T X with X ∼ N (µ, Γ), Γ ∈ R d×d a positivedefinite matrix and B ∈ R d×p . In this example, we consider d = 2 and p = 3 which means Y = (Y 1 , Y 2 , Y 3 ). The variance of the centered random variables X 1 and X 2 are equal to σ 2 1 = 1 and σ 2 2 = 3, respectively and their correlation ρ = 0.4. Thus the covariance matrix of X is given by:

Γ = σ 2 1 ρσ 1 σ 2 ρσ 1 σ 2 σ 2 2 = 1 0.69 0.69 3 ,
and the coefficients of B = (β ij ) ∈ R 2×3 are chosen as:

B = 1 4 0.1 1 3 0.9 .
The variance of the output Y j with j ∈ {1, 2, 3} is σ 2

Y j = β 2 1j σ 2 1 + 2ρβ 1j β 2j σ 1 σ 2 + β 2 2j σ 2 2 .
The scalar Shapley effects are:

σ 2 Y j φ j 1 = β 2 1j σ 2 1 1 - ρ 2 2 + ρβ 1j β 2j σ 1 σ 2 + β 2 2 σ 2 2 ρ 2 2 , σ 2 Y j φ j 2 = β 2 2j σ 2 2 1 - ρ 2 2 + ρβ 1j β 2j σ 1 σ 2 + β 2 1 σ 2 1 ρ 2 2 .
Then, the aggregated Shapley effects for i ∈ {1, 2} are computed according to (5.10).

Let us first focus on scalar Shapley effect estimation and the associated confidence intervals, for example scalar Shapley effects for the first component of the output, Y 1 . For Y 1 , the most important input is X 2 , with a Shapley effect around 0.66. In Figure 5.1, we analyze the evolution of the estimation accuracy through the mean absolute error (MAE) and the POC as a function of n and N tot . The value of N tot varies from 100 to n(2 d -2) = 2n. For each combination of n and N tot , N = 300 independent random samples are used for estimating the MAE and the POC. The bootstrap sample size B is fixed to 500. The evolution of the 5% and 95% quantiles of the MAE is also drawn. As expected, the estimation accuracy increases with n and N tot . Also, for fixed n, it increases with N tot . As N tot reaches n, which means N u = n/2 for all u as d = 2, the decrease of the MAE slows down from exponential to linear. This behavior seems coherent with the result stated in [Broto et al., 2020, Corollary 6.8] which asserts that there exits 0 < C < +∞ such that, for N u ≥ Cn 1/(d-|u|) and δ > 0:

| c 1 (u) -c 1 (u)| = o p 1 n 1/(2(d-|u|)-δ .
(5.13)

The bound in (5.13) does not depend on N tot anymore. Finally, the behavior of the POC is as expected, around 0.9 whatever n and N tot . The evolution of the MAE and POC for the estimation of aggregated Shapley effects as function of n and N tot is displayed in Figure 5.2. Similarly as for the estimation of scalar effects, the accuracy increases with n and N tot . However for fixed n, the POC deteriorates for large N tot (POC is around 0.8 for N u = n). The reason for this deterioration is not well understood but is probably due to the use of bootstrap with nearest neighbors. In the following, we keep all the N u strictly below n to guarantee the POC.

We estimate Shapley effects if input correlation is higher (ρ = 0.9). POC and bias results are also satisfactory (not shown). In fact, POC values vary also around 0.9 and bias decreases and goes to 0 when n and N tot increases. For this simple test case, we have shown that confidence intervals using Algorithm 1 reach accurate coverage probability and that bias reduces when n and N tot increase. Nevertheless in this test case, estimation is effortless because d = 2. The empirical coverage probability of the 90% bootstrap simultaneous confidence intervals is displayed with dotted lines. The theoretical coverage probability 0.9 is also shown with a gray plain line. The bootstrap sample size is fixed to B = 500.

Input Description Distribution

m mass (kg) U[10, 12] c damping constant (Nm -1 s) U[0.4, 0.8] k spring constant (Nm -1 ) U[70, 90] l initial elongation (m) U[-1, -0.25]
Table 5.1: Mass-spring model: Input description and uncertainty intervals. U denotes the uniform distribution. The inputs are independent from each other.

Mass-spring model

The method is illustrated on a test case with discretized functional output: the functional mass-spring model proposed by [Gamboa et al., 2013], where the displacement of a mass connected to a spring is considered:

m (t) + c (t) + k (t) = 0, (5.14)
with initial conditions (0) = l, (0) = 0, and t ∈ [0, 40]. There exists an analytical solution to Equation (5.14). This model has four inputs (see more details in Table 5.1). The model output is the vector Y = f (X) = ( (t 1 ), . . . , (t 800 )), t i = 0.05i with i ∈ {1, . . . , 800}.

Inputs are considered independent. The true aggregated Shapley effects are unknown but they are approximated using a high sample size n = 25 000 and N tot = 10 000. Then, the Shapley effects estimated are GS m = 0.38, GS c = 0.01, GS k = 0.51 and, GS l = 0.09. Given these results, inputs ranking is: k, m, l and c which corresponds to the same ranking obtained using Sobol' indices (see Table 3 of [Gamboa et al., 2013]).

The discretized output is high-dimensional (p = 800). We perform fPCA (see Subsection 5.3.3) to estimate the effects using the first q p fPCs. Figure 5.3 shows the POC and bias evolution if different values for n and N tot are used for the aggregated effect estimation. We use the first 6 fPCs which explain 95% of the output variance (see Figure 5.3 a). For each n and N tot combination, the aggregated Shapley effects are estimated for N = 100 independent samples and confidence intervals are estimated with B = 500 bootstrap samples. The accuracy is not very good if sample size is small n = 1000 (see Figure 5.3 b). However, it reduces drastically when the sample size n increases as expected. On our experiments, the smallest MAE is achieved for n = 5000 and N tot = 2002 (see Figure 5.3 d). The confidence interval reaches the expected POC 0.9 as soon as n and N tot are large enough (see Figure 5.3 b).

Snow avalanche modeling

Model

Our avalanche model is based on depth-averaged Saint-Venant equations and considers the avalanche as a fluid in motion. In more detail, it considers only the dense layer of the avalanche. The flow depth is then small compared to its length. The model assumes the avalanche is flowing on a curvilinear profile z = l(x), where z is the elevation and x is the projected distance measured from the top of the avalanche path. Under these Explained variance q q q q q q q q q q a) 0.02 0.04 0.06 0.08 0.10

n = 1000

Ntot Mean absolute error The probability of coverage of the 90% bootstrap simultaneous intervals is displayed with a dotted line. The theoretical POC 0.9 is highlighted with a plain gray line. The bootstrap sample size is fixed to B = 500.

assumptions, shallow-water approximations of the mass and momentum equations can be used:

∂h ∂t + ∂hv ∂x = 0 ∂hv ∂t + ∂ ∂x hv 2 + h 2 2 = h(g sin θ -F)
where v = v is the flow velocity, h is the flow depth, θ is the local angle, t is the time, g is the gravity constant and F = F is a frictional force. The model uses the Voellmy frictional force F = µgcosθ + g/(ξh)v 2 , where µ and ξ are friction parameters. The partial differential equation (PDE) system is solved numerically using a finite volume scheme on the path curvilinear profile with a 5m resolution (see full details about the PDE scheme in [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF]). Hence, even if the PDEs are expressed in a two-dimensional frame, they in fact represent a one-dimensional flow on the curvilinear profile ( [START_REF] Greve | Unconfined flow of granular avalanches along a partly curved surface. i. theory[END_REF]).

The numerical model depends on six inputs: the friction parameters µ and ξ, the length l start of the avalanche release zone, the snow depth h start within the release zone, the abscissa corresponding to the beginning of the release zone denoted by x start and the discretized topography of the flow path, denoted by D = (x, z) ∈ R Ns×2 where x ∈ R N S is the vector of projected abscissa positions and z = l(x) ∈ R N S is the elevation vector. N s is the number of points of the discretized path. Subsection 5.6.1 shows an illustration of the avalanche model. We use for D the topography of a path located in Bessans, France. This particular path already considered in other works [Eckert et al., 2008a[START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF], Favier et al., 2014b, Eckert et al., 2018b] and one scalar outputs and these three outputs are the objects of the GSA study. We develop our GSA in two contexts or scenarios: (i) little knowledge on the input parameter probability distribution, and (ii) well-calibrated input parameter distributions.

Scenario 1 Principle

We first determine the most influential input parameters sampled from uniform distributions. We thus expect from the GSA a better understanding of the numerical model. Inputs µ and ξ vary in their physical value ranges. Inputs l start and h start vary in their spectrum of reasonable values given the characteristics of the avalanche path. The x start input distribution is determined by calculating the abscissa interval where the release zone average slope is superior to 30 • . Indeed, the slope remains above 30 • during the first 1600m of the path. A good approximation of avalanche release zones is commonly obtained this way. Since different studies [Bartelt et al., 2012, Brian Dade and[START_REF] Dade | Long-runout rockfalls[END_REF]] suggest that the volume of snow is a critical quantity that controls flow dynamics, we consider vol start as input of the GSA instead of h start and l start . The latter is evaluated as vol start = l start ×h start ×72.3/ cos(35 • ). The mean width and slope of the release zone equal to 72.3m and 35 • , respectively. All uncertainty intervals are summed-up in Table 5.2. The input correlations are close to 0 since we assume they are a priori independent. For a given avalanche simulation, its functional velocity and flow depth outputs have a high number of zeros because they are null above the beginning of the release zone and after the runout position. Also, there might be some avalanche simulations that are meaningless in physical terms and/or not useful to assess the related risk. Therefore to perform GSA, we select simulations that accomplish the following rules: (i) avalanche simulation is flowing in the interval [1600m, 2412m],

(ii) its volume is superior to 7000 m 3 and, (iii) avalanche runout distance is inferior to 2500m which corresponds to the end of the path.

In the following, we use the terminology constrained sample to deal with the sample obtained after applying rules (i) to (iii).

The return period of avalanches in the interval [1600m, 2412m] varies from 1 to 10 000 years according to the work of [START_REF] Eckert | Long-term avalanche hazard assessment with a Bayesian depth averaged propagation model[END_REF]. Roughly speaking, a return period is the mean time in which a given runout distance is reached or exceeded at a given path's position. Also, we focus on medium, large and very large avalanches which have a high potential damage. correlation original / AR 0/0.53 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 2500 5000 7500 10000 0.2 0.4 0.6 µ ξ 0e+00 5e-05 1e-04 0 2500 5000 7500 10000 ξ density correlation original / AR 0/-0.12 correlation original / AR 0/-0.43 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 500 1000 1500 0.2 0.4 0.6 µ xstart q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 500 1000 xstart density correlation original / AR 0/-0.13 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20000 40000 60000 0.2 0.4 0.6 µ volstart q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20000 40000 60000 0 2500 5000 7500 10000 ξ volstart q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20000 40000

Global sensitivity analysis results

We first ran n 0 = 100 000 avalanche simulations from an i.i.d. sample of input distributions described in Table 5.2. Then, by applying (i) to (iii) the size of the constrained sample was reduced to n 1 = 6152 ( Subsection 5.6.2). Even if the initial sample size is high (n 0 = 100 000) and if the corresponding input parameter sample does not present any significant correlation structure, the size of the constrained sample is low and we can observe a correlation structure. For example, inputs µ and ξ were independent for the initial sample but the correlation computed after applying the constraints is 0.31. Note that the input parameter correlations induced by the constraints were the main motivation to compute Shapley effects and not Sobol' indices in this first scenario. On Subsection 5.6.2 are plotted highest density region (HDR) boxplots for the velocity and flow depth, obtained by using the R package rainbow developed by [START_REF] Hyndman | Rainbow plots, bagplots, and boxplots for functional data[END_REF]. The HDR boxplot is a vizualization tool for functional data based on kernel density estimation of the scores associated to the two first principal components of the functional data (see [START_REF] Hyndman | Computing and graphing highest density regions[END_REF] for further details). In the data we consider in our study, the avalanche velocity ranges from 0.1ms -1 to 71.56ms -1 and avalanches are decelerating (see Subsection 5.6.2 a). Flow depths vary from 0.03m to 7.52m. The flow depth curves exhibit high fluctuations in the interval [2100m, 2300m] (see Subsection 5.6.2 b) which corresponds to a region where path's topography is mostly convex. Runout distances vary from 2409m to 2484m (see Subsection 5.6.2 c). On Subsection 5.6.2 panels a and b, ubiquitous (pointwise) Shapley effects of velocity and flow depth curves are shown, respectively. Depending on the output, results are quite different. For velocity, x start is the most relevant during a large part of the path but its importance decreases along the path and, conversely, the importance of the other inputs increases. For the flow depth output, the most important input is vol start since the corresponding Shapley effects vary from 0.4 to 0.2 along the path. Nevertheless, other inputs are not completely negligible. Input importance also varies according to the topography. In fact, the ubiquitous effect variation corresponds to local slope changes (see Subsection 5.6.2 a and b). Correlations between ubiquitous effects and local slope have been computed and are rather high. For example, for the velocity, the absolute value of the correlation is higher than 0.51 for all input parameters. This implies that local slope changes play an important role on the input contribution to output variations, a nice results showing the relevance of the GSA analysis to understand the dynamical properties of the flow. Eventually, for runout distance, the four inputs are relevant. Our results are in accordance with previous works documenting the sensitivity of snow avalanche models to their input conditions using less formalised methods [START_REF] Borstad | Sensitivity analyses in snow avalanche dynamics modeling and implications when modeling extreme events[END_REF]. They also relate to more general results regarding snow avalanche dynamics. Especially, it is well known that in the transit area of the path, an avalanche accelerates up to a maximal velocity and this limit velocity can be explicitly computed on a constant slope under the assumption of a Voellmy friction law [START_REF] Bruno Salm | Berenchnung von Fliesslanwinen, eine Anleitung für Paktiker mit Beispielen[END_REF][START_REF] Gauer | Comparison of avalanche front velocity measurements and implications for avalanche models[END_REF][START_REF] Mcclung | Maximum frontal speeds, alpha angles and deposit volumes of flowing snow avalanches[END_REF]. It is therefore all the more logical that the influence of x start on velocity decreases as the avalanche flows downslope. More generally, with depth-averaged flow equations, it is known that initial conditions corresponding to the avalanche release are important only in the early phases of the flow and are forgotten later, especially in the runout area where velocity and flow depths become primarily constrained by local slope. The fPCs convey useful information about the avalanche model and also, its sensitivity. Subsection 5.6.2 (resp., Subsection 5.6.2) shows the correlations between the fPCs and the velocity (resp., flow depth) output for scenario 1. The first fPC has a positive correlation with the velocity, always higher than 0.57, however the correlation decreases drastically from the top to the bottom of the path. By contrast, the second fPC is less correlated to the velocity, the correlation is positive at the top of the path, then the correlation decays and even becomes highly negative for x ≥ 2300 m. Hence, these two components identify well two regions of the path where velocity behaves differently: the first fPC corresponds to the release and transit area of the path, where avalanches accelerate and maintain their velocity while the second fPC identifies the runout area where avalanches are decelerating and eventually stop. The behavior is similar for the flow depth, with a sharp peak in the second fPC corresponding to the location in the runout area where slope becomes null and then negative. To strengthen our physical interpretation, we estimated the scalar Shapley effects for each fPC with n = 6152 and N tot = 2002. As far as velocity is concerned, for the first fPC, x start is the most important input and for the second fPC, vol start , x start and µ are the most relevant. For the flow depth, for the first fPC, vol start is the most relevant. For the second fPC, the inputs x start , vol start and µ are the most important inputs. These results are in accordance with previous works documenting the sensitivity of snow avalanche models to their input conditions, and notably the high importance of the Coulomb friction coefficient µ and snow volume in the runout area [START_REF] Dent | Modeling of snow flow[END_REF][START_REF] Borstad | Sensitivity analyses in snow avalanche dynamics modeling and implications when modeling extreme events[END_REF][START_REF] Bartelt | Modeling mass-dependent flow regime transitions to predict the stopping and depositional behavior of snow avalanches[END_REF]. Similar results are obtained for scenario 2 (see Figures 5.16 and 5.17 Subsection 5.6.2 shows aggregated Shapley effects and 90% confidence intervals computed over space intervals [x, 2412] where x ∈ {1600, 1700, . . . , 2412}. The aggregated effects are computed in the first fPCs explaining more than 95% of the output variance. Aggregated effects seem more robust than ubiquitous effects, specially in zones where local slope shows high variations (see Subsection 5.6.2 compared to Subsection 5.6.2). For explaining more than 95% of the velocity output variance, 2 fPCs are required, while, for explaining more than 95% of the flow depth output variance, at most 4 fPCs are required, depending on x. For the velocity output, the most important input is x start in the interval [1600m, 2100m] but its importance decreases along the path. In the interval [2017m, 2412m] where return periods are non trivial, x start and vol start are the most important followed by µ and ξ. For the flow depth output, vol start is the most relevant but its importance decreases along the path. At the end of the path from 2300m to 2412m where return periods are high (between 100 to 10 000 years), confidence intervals intersect. It seems thus difficult to deduce a clear ranking of the inputs for these last portions of the path. Nevertheless, it seems that none of the inputs is negligible, even at the very end of the path. In summary, to estimate velocities with accuracy, the release zone and volume are the most important parameters and, for the flow depth, a good approximation of the volume released is essential. q q q q q q q q q 0.1 0.3 0.5 0.7 velocity abscissa position [m] GSh q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q µ ξ x a) q q q q q q q q q 0.1 0.2 0.3 0.4 0.5 0.6

flow depth

abscissa position [m] GSh q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1600 1800 2000 2200 2415 Aggregated effects (see Subsection 5.6.3) present less fluctuations and are easier to interpret (see Subsection 5.6.3). In summary, under this second scenario, it is fundamental to have a good approximation of the released volume and abscissa for velocity forecasting, while for flow depth forecasting, a good approximation of released volume is desirable. Nevertheless, none of the other inputs are negligible. Note that the uncertainty associated to the estimation of Shapley effects at 2204m is high (see the width of the corresponding confidence intervals on Subsection 5.6.3). To outperform the estimation accuracy at the end of the path generating a larger initial sample of avalanches is possible, but the computational burden is prohibitive.

Conclusions and perspectives

In this work, we extended Shapley effects to models with multivariate or functional outputs. We proved that aggregated Shapley effects accomplish the natural requirements for a GSA measure. For the estimation, we proposed to extend the subset aggregation procedure with double Monte Carlo given data estimator of [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF]. We proposed an analysis of the impact of the total budget N tot , the number of neighbors N I and the optimal allocation N * u on the accuracy of Shapley effect estimation. Also, we proposed an algorithm to construct bootstrap confidence intervals for scalar and aggregated Shapley effects based on the ideas of [START_REF] Benoumechiara | Shapley effects for sensitivity analysis with dependent inputs: bootstrap and kriging-based algorithms[END_REF]. In test cases, the convergence of our estimator was empirically studied. Also, we empirically demonstrated that the bootstrap confidence intervals we proposed have accurate coverage probability. Estimation and bootstrap confidence interval algorithms well behave. Nevertheless, high sample sizes (n = 5000 and N tot = 2002) are required to guarantee accurate results. Remark that it is well known that Shapley effects estimation is costly. It would be interesting to study theoretically the asymptotic properties of our estimator, but this q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 0.6 velocity abscissa position [m] GSh q q q q q q q q q q q q q q q q q q q q q 1600 1800 2000 2205 10 20 30 40

Local slope [degrees] a) q q q q q q q 0.0 0.2 0.4 0.6 flow depth abscissa position [m] GSh q q q q q q q q q q q q q q q q q q q q q 1600 1800 2000 2205 However, we did not use this function in this is work because there is no paper that analyzes the convergence of the estimator and, this subject was out of the scope of this work. We rather used the shapleySubsetMc function which corresponds to the estimator introduced in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF] on which our estimator for aggregated Shapley effects is based.

We applied our GSA methods to an avalanche propagation model under two different settings. Due to physical modeling constraints input parameters were not confined in a rectangular region. For these reasons, it was not possible to consider independence of input parameters. Results showed probative linkages between local slope and sensitivity indexes. Notably, aggregated Shapley effects were more stable and easier to interpret than ubiquitous effects, as already observed by [START_REF] Alexanderian | Variance-based sensitivity analysis for time-dependent processes[END_REF] in the case of aggregated Sobol' indices. This demonstrates the usefulness of our approach for many practical problems. Especially, it could be applied to other avalanche paths to generalize the results obtained in terms of respective weights of the inputs and interpret the sensitivity indices more deeply in terms of physical properties of avalanche flows. Eventually, application was challenging because constrained samples were of moderate size, for example, from the 100 000 initial sample, the physical constraints produced a 6000 to 1200 sample, depending on the setting. In a future work, it would be useful to construct a surrogate of the avalanche model to generate larger constrained samples, improve the accuracy of aggregated Shapley effect estimation and thus reduce confidence intervals width. q q q q q q q q q 0.25 0.35 0.45 m Number of PCs GSh 1 q q q q q q q q q 1 2 3 4 5 6 7 8 800 q q one sample two samples a) q q q q q q q q q 0.00 0.04 0.08 c Number of PCs GSh 2 q q q q q q q q q 1 2 3 4 5 6 7 8 800 b) q q q q q q q q q 0.50 0.60 0.70 k Number of PCs GSh 3 q q q q q q q q q 1 2 3 4 5 6 7 8 800 c) q q q q q q q q q 0.05 0.10 0.15 l Number of PCs GSh 4 q q q q q q q q q 1 2 3 4 5 6 7 8 800 d)

Figure 5.15: Mass-spring: estimation and associated confidence intervals for aggregated Shapley effects using a single sample of size n = 5000 to perform the dimension reduction and to estimate the effects or two samples (one to perform the dimension reduction of size n 1 = 1000 and another one to estimate the effects n 2 = 4000).

The estimation as a function of the basis size (number of principal components) is displayed. All the pannels were produced with N tot = 2000 and B = 500 (bootstrap sample size for confidence intervals).

Estimation of Shapley effects using one or two samples

We compare the accuracy of Shapley effect estimation for the mass-spring test case when using a single sample of size n = 5000 or two independent samples of size n 1 = 1000 and n 2 = 4000 for estimating the KL representation and the aggregated Shapley effects. In Figure 5.15 we see that the procedure without splitting produces better results for three of the four inputs. For the last input, results with two samples are slightly better. Based on these results, we decided to apply the procedure without splitting on our snow avalanche application.

5.11 Functional principal components and Shapley effects for scenario 2 5.16 (resp.,Figure 5.17) shows the correlations between the fPCs and the velocity (resp., flow depth) output for scenario 2. 6

Conclusions and perspectives

In this thesis, we addressed the calibration and global sensitivity analysis of avalanche dynamic models. In general, we introduced and proposed advanced statistical methods for understanding better input parameters and outputs of snow avalanche models. All the techniques were illustrated with synthetic data and with real snow avalanche data to demonstrate their efficiency and usefulness for real applications. Our contributions may be useful for avalanche scientists in order to (i) account for the error structure in model calibration and (ii) to rank input parameters according to their importance in the output model.

• A framework based on Bayesian inference to calibrate avalanche dynamic models using high temporal resolution data.

• A nonparametric method based on Nadaraya-Watson kernel with bootstrap bias correction to estimate aggregated Sobol' indices from a given sample.

• A given data estimation method and associated confidence intervals for aggregated Shapley effects based on nearest-neighbor estimation.

The last two methods were applied to rank input parameters relatively to their impact on the variability of functional outputs of an avalanche dynamic model. In the next sections, we revisit in detail the contributions and conclusions of this manuscript and finally, we present the perspectives and open research questions.

Conclusions

Accounting for error structure in model calibration

In Chapter 3, we proposed a Bayesian inference framework to calibrate an avalanche dynamic model from an avalanche velocity time series with ∆t = 1s. The time series was obtained from high rate positioning photogrammetry images of an avalanche released at the Lautaret test site on February 13 2013. We considered a sliding block model with Voellmy friction law. To perform the calibration under a Bayesian framework, assumptions about the model error distributions (defined as observationsmodel simulation) are required. Traditionally, the errors are assumed i.i.d. and normally distributed. However, as our data is a time series, the errors are autocorrelated. Thus, the objective of our work was to show how taking into account the autocorrelation might improve the friction parameter estimation and the prediction of velocities. Within this purpose, we compared two statistical models: a first one that considers i.i.d. Gaussian centered errors and a second one that considers the errors modeled by an autoregressive process of order 1. We also generated synthetic data to illustrate how omitting the autocorrelation generates a bias for the estimation of friction parameters. Results showed that it is important to take into account the error structure (in our case, the autocorrelation) for model calibration to avoid bias.

Non parametric estimation method for aggregated Sobol' indices

In Chapter 4, our objective was to perform a GSA of an avalanche model. We considered the particular context of independent input parameters. The other constraint was that we had to propose a given data estimation procedure, accurate even for a sample of moderate size. Within this purpose, we proposed a nonparametric method to estimate the aggregated Sobol' indices from a given data sample. The nonparametric method was based on the Nadaraya-Watson kernel smoother. Since the sample size was moderate, the kernel estimation was biased and it required a bias correction. To remove such a bias, we proposed a bias correction based on the works of [START_REF] Racine | Bias-Corrected Kernel Regression[END_REF] and [START_REF] Solís | Non-parametric estimation of the first-order Sobol indices with bootstrap bandwidth[END_REF]. We tested the accuracy of the method on toy functions, and results were good even if sample sizes were moderate (n = 300). Then, we applied the method to estimate aggregated Sobol' indices of the depth averaged avalanche model proposed by [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF]. In this study, we analyzed the sensitivity of the model for the avalanche event used in Chapter 3. More precisely, the friction parameter uncertainty intervals were given by the posterior Bayesian inference distributions and the uncertainty for other model input quantities were obtained from measures taken during the avalanche release. In this chapter, the input parameters were assumed to be independent. For this particular avalanche event, the input ξ, followed by µ, were the most important for the velocity output. For the flow depth and runout distance outputs, the snow depth within the released zone was fundamental.

Aggregated Shapley effects with bootstrap confidence intervals

In Chapter 5, our objective was to perform GSA of avalanche dynamic models in a more general framework than the one of Chapter 4. Indeed, in Chapter 4 we considered the input distribution associated to a specific avalanche event to guarantee the input independence. To perform GSA in a more general framework, we had to guarantee that all the samples were physically meaningful. Non physical sets of input parameters were then removed from the experimental design. For example, some combinations of friction parameters and released volume were preventing from an avalanche release or were producing small avalanches (volume below a threshold) which were not interesting for our study. Within this purpose, we developed a GSA method that can be applied to samples obtained from acceptance/rejection (AR) rules (corresponding to the physical constraints of the study). AR rules were leading to dependence among input parameters. Therefore we proposed to estimate Shapley effects extended to multivariate or functional outputs, which are more meaningful in the framework of dependent inputs. We called these sensitivity measures aggregated Shapley effects. For their estimation, we extended the given data estimator of [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF]. We also built confidence intervals by adapting ideas of [START_REF] Benoumechiara | Shapley effects for sensitivity analysis with dependent inputs: bootstrap and kriging-based algorithms[END_REF]. Then, we applied our GSA method to the avalanche dynamic model proposed by [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF]. We found that the release position and volume were the most important inputs. Also, by looking at the evolution of ubiquitous Shapley effects along the path, we noticed that it might exist a relationship between the effects and the local slope.

Perspectives

For calibration using Bayesian inference

To reduce computational cost A drawback of the framework proposed is that the Metropolis-Hasting algorithm, a MCMC algorithm used for sampling from the parameter posterior distribution, might require a high number of model runs to converge. In our application, we required at least 1000 model runs and computational cost was not a problem since, model runs were fast. However, for using Bayesian inference to calibrate a more complex model (with high number of parameters and high computational cost) some solutions could be: (i) to perform a GSA of the model with the aim of identifying the most important input parameters. Then, to only calibrate those ones and, fix the others to reference values. (ii) To construct a metamodel for the avalanche model and to perform the calibration with the metamodel. Users must also take into account the error induced by approximating the model with the metamodel. (iii) To apply other MCMC algorithms that might require less model runs to converge as for example the adaptive Metropolis-Hasting algorithm introduced by [START_REF] Haario | An adaptive metropolis algorithm[END_REF]. In the adaptive MH, the proposal distribution q is updated in each iteration which makes convergence faster. Nevertheless, adaptive MH is non Markovian and some other theoretical nice properties of the MH algorithm might be lost.

To include other quantities within the calibration

In the study, we only used the avalanche velocity for the calibration. Thus, it could be interesting to add other quantities as for example, flow depth time series and runout distance. To do so, an appropriate likelihood function should be proposed (e.g., Gaussian process or Dirichlet process). However, a higher quantity of data would also be needed to optimize the new likelihood.

To apply the approach to several snow avalanches It could be useful to apply the same approach to several avalanches to investigate, e.g., the relationship between friction parameters and snow conditions. Such findings could improve deeply our knowledge about the physics of the avalanche phenomenon. However, to do so, it would be necessary to acquire high-resolution avalanche data and the snow conditions of the snowpack for the same the avalanche events.

For the Global sensitivity analysis

To construct confidence intervals for non parametric indices estimation It would be important to build confidence intervals for nonparametric bias corrected estimation of aggregated Sobol' indices. This task is not trivial because we cannot use bootstrap based confidence intervals since we already use bootstrap for correcting the bias and a double bootstrap is far too costly. Some interesting clues to explore could be (i) to use Bayesian inference to deduce credible intervals for the Sobol' indices [see Antoniano-Villalobos et al., 2019] and (ii) to prove asymptotic convergence properties of the aggregated Sobol' estimator. Nevertheless, the proof of such properties is a challenging task.

To apply GSA to other avalanche events and models It could be useful to apply our GSA methods to other snow avalanche models and snow avalanche events to see if results could be generalized. The methods should be tested for example for 3-D spatio-temporal avalanche dynamic models [as for example the r.avalflow model of [START_REF] Mergili | Pudasaini. r.avaflow v1, an advanced open-source computational framework for the propagation and interaction of two-phase mass flows[END_REF]. The output of such models are functional and dependent, thus the methods should be adjusted to this kind of outputs [see for example [START_REF] Marrel | Global sensitivity analysis for models with spatially dependent outputs[END_REF].

To include the topography as input parameter

In our studies, the topography was fixed but it could also be included in the GSA as a functional input. In such a case, our GSA method should be adapted to functional inputs [see for example, Iooss andRibatet, 2009, Fruth et al., 2015]. Also, the fact of including the topography as input could provide more insights about the potential connection between the input parameters and the topographical characteristics of avalanche paths.

To compare our ranking with other input rankings

The GSA framework applied in Chapter 5 could be seen as a particular case of conditional sensitivity analysis (CSA) which aims at measuring the influence of the inputs on the output within a restricted domain [see [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF]. In our application, the output domain was restricted by physical constraints. Thus it could be interesting to compare our rank to the input rank obtained by using the aggregated HSIC indices of [START_REF] Da | Global sensitivity analysis with dependence measures[END_REF]. Moreover, Shapley effects might be used within the CSA framework, but this should be first study in test cases.

A.2 The spectral approach EASI

The effective algorithm for computing global sensitivity indices (EASI) proposed by [START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF] is a spectral method to estimate first-order Sobol' indices using the Fast Fourier Transform. Spectral approaches (as FAST and EASI) should be preferred if the model f has some regularity [START_REF] Prieur | Variance-Based Sensitivity Analysis: Theory and Estimation Algorithms[END_REF]. Let (x 1 , . . . , x d , y ) ∈{1,...,n} be a n sample of (X, Y ). The sample vector x i = (x i ) ∈{1,...,n} of input X i is ordered increasingly to obtain a new vector x ( ) i which elements satisfy x

(1) i ≤ . . . ≤ x (n) i . Then, all odd indices from (x ( ) i ) ∈{1,...,n} are ordered increasingly and followed by all even indices in decreasing order to obtain the vector (x [ ] i ):

x

[ ] i = x (2 -1) i ≤ n+1 2 , x 2(n+1-) i > n+1 2
The elements of (x [ ] i ) satisfy a zig-zag relation:

x

[ ] i ≤ x [ +1] i if ≤ n + 1 2 and x [ ] i ≥ x [ +1] i if ≥ n + 1 2 .
Let us denote by π the permutation of {1, . . . , n} such that x π( ) i

= x

[ ] i

The first-order sensitivity index is then estimated by:

Si = 2 M m=1 |c m | 2 m =0 |c m | 2 (A.1)
where c m is the m-complex coefficient of the discrete Fourier transform of (y π( ) ) ∈{1,... , with m = 0, ±1, ±2, . . . , n 2 with . denoting the integer part. The maximum harmonic M is usually taken as 4 or 6 [START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF]. The estimation introduced in (A.1) is biased with respect to the maximal number of harmonics M , thus [START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF] proposed a bias correction which we use in the test cases:

Ši = 1 n -2M
(n Si -2M ) (A.2)

A.3 Test cases

We study the accuracy of NSD and EASI on two test cases. The methods are also compared to the pick-and-freeze method based on replicated orthogonal array proposed by [START_REF] Thierry | Comparison of some efficient methods to evaluate the main effect of computer model factors[END_REF] and further studied in [START_REF] Tissot | A randomized Orthogonal Array-based procedure for the estimation of first-and second-order Sobol' indices[END_REF], we denoted it by sobol-roalhs. This last method requires specific designs thus is not a given data method. It is used as a benchmark approach. The R package sensitivity is used for its estimation.

q q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 S 1 100 200 500 1000 q q q sobol. NSD EASI a) q q q q q q q q q q q q q q q q q q q q q q q q q 0.2 0.3 0.4 0.5 0.6 0.7 S 2 100 200 500 1000 b) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -0.6 -0.4 -0.

A.3.1 Ishigami function

The Ishigami function [START_REF] Ishigami | An importance quantification technique in uncertainty analysis for computer models[END_REF]] is a classical test function which is nonlinear and non monotonic:

f (X 1 , X 2 , X 3 ) = sin(X 1 ) + 7 sin 2 (X 2 ) + 0.1X 4 3 sin(X 1 ), the inputs X 1 , X 2 , X 3 are uniformly distributed on [-π, π]. The number of harmonics of EASI was set to 6 following the recommendation of [START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF]. When sample sizes are small, the three methods show small bias. However, as the sample size increases the bias reduces. EASI shows accurate results even with small sample sizes for S 1 and S 3 (n = 200). NSD does not estimate accurately indices close to 0 (S 3 ).

A.3.2 Bartley et al. function

The second test function is the Bartley et al. function: f (X 1 , . . . , X 6 ) = 6 i=1 i j=1 X j , (A.3) X = (X 1 , X 2 , . . . , X 6 ) is assumed to be uniformly distributed on [0, 1] 6 . In conclusion, within these two example we showed that EASI is a nice option to estimate Sobol' indices with small given data samples. The choice of M does not seem to q q q q q q q q q q q q q q q q q q q q q q q q q 0.4 sobol. NSD EASI a) q q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 S 2 100 200 500 1000 b) q q q q q q q q q q q q q q q -0.3 -0.1 0.0 0.1 0.2 0.3 S 3 100 200 500 1000 c) q q q q q q q q q q q q q q q -0.2 0.0 0.2 0.4 S 4 100 200 500 1000 d) q q q q q q q q q q q q q q q q q q q -0.3 -0.1 0.1 0.2 0.3 0.4 S 5 100 200 500 1000 e) q q q q q q q q q q q q q -0. impact that much the results. However, for highly non regular models we could expect that K = 6 is too small, even if the reordering of the sample regularizes the functions. NSD is biased, and to our knowledge, there does not exist any selection procedure for the number of classes for small to moderate sample sizes. q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 10 20 30 40 0.15 0.20 0.25 0.30 0.35 0.40 0.45

Ubiquitous effects Aggregated effects Flow velocity [m s

Local slope [degrees] Sh q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q a) q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 10 20 30 40 0.0 0.2 0.4 0.6 0.8

Local slope [degrees] Sh q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q µ x start vol start

Ubiquitous effects Aggregated effects Flow velocity [m s

  Contents 2.1 A short introduction to 1-D snow avalanche dynamic modeling . . . . 5 2.2 Model calibration using Bayesian inference . . . . . . . . . . . . . . . 8 2.3 Global sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Dimension reduction techniques for functional outputs . . . . . . . . . 21

Figure 2 . 1 :

 21 Figure 2.1: Avalanche block model.

Figure 2 . 2 :

 22 Figure 2.2: SA framework.

Figure 2 . 3 :

 23 Figure 2.3: Morris function. μi and s i estimated with n = 30, p = 8, ∆ = 4/7 (a total of 630 model runs). Black line corresponds to µ = 2s/ √ 30.

  1. (Efficiency) d i=1 φ i = val({1, . . . , d}). 2. (Symmetry) If val(u∪{i}) = val(u∪{ }) for all u ⊆ {1, . . . , d}-{i, j}, then φ i = φ . 3. (Dummy) If val(u ∪ {i}) = val(u) for all u ⊆ {1, . . . , d}, then φ i = 0. 4. (Additivity) If val and val' have Shapley values φ and φ respectively, then the game with characteristic function val + val' has Shapley value φ i + φ i for i ∈ {1, . . . , d}.

Figure 2 . 6 :

 26 Figure 2.6: FHDR boxplots of a random n = 5000 sample of avalanche runs a) flow depth, b) flow velocity. c) Some examples (5 selected runs) of released and runnout positions. The topography (black line) corresponds to an avalanche path at the Lautaret test site.

Figure 3 . 1 :

 31 Figure 3.1: 3D-topography of path number 2 of the Lautaret test-site, some front positions (black lines) of the avalanche released on 13 February 2013 as determined from photogrammetric measurements and changes in snow depth before/after the avalanche as inferred from terrestrial laser scanning.

)Figure 3 . 3 :

 33 Figure 3.3: Samples of the generated synthetic longitudinal velocity profiles (colour lines) and observations corresponding to the studied avalanche (dots).

Figure 3 .

 3 Figure3.4 shows the parameter prior and posterior densities according to the MH samples. In Figure3.4 panels a) to c) (resp., Figure3.4 panels d to g), the M 0 densities are shown (resp., M 1 ) and Table3.4 sums up the corresponding descriptive statistics. Eventually, posterior correlations between model parameters are shown in Figures 3.5 and 3.6 for M 0 and M 1 , respectively. For the two models, the marginal posterior distributions have lower variance than their priors, meaning that the observations have conveyed information into the analysis (see Figure3.4), an expected result.The variances of the posterior distribution of the friction parameters µ and ξ in the M 0 model are lower than in M 1 model. This result (arguably the sole undesirable one with M 1 compared to M 0 ) could be a consequence of the interactions between the friction coefficients and error parameters in the more parametrized M 1 model (four free parameters instead of three with M 0 ). Indeed, there is strong correlation between the parameters ξ and φ (0.49) for the model M 1 , which may preclude reaching sharp estimates of friction parameters with model M 1 . Conversely, even if there is a high correlation between the parameters µ and ξ with both models, switching from M 0 to M 1 reduces it from 0.85 to 0.50. The high correlation between the two parameters of the Voellmy friction law is known since the first calibration approaches of, e.g.[START_REF] Dent | Modeling of snow flow[END_REF]. This usually limits robust interpretation of obtained estimates, so that reducing it thanks to M 1 should be seen as advantageous. This correlation reduction may be a collateral effect of having one more free parameters, allowing a bit more flexibility to fit the data. Similar effect is reflected in the lower correlation between σ 2 and the friction parameters with M 1 than with M 0 . Note by the way that fairly assessing such correlations is a real strong point of our formal Bayesian calibration approach.Remarkably, posterior estimates of friction parameters are very contrasted under both models. Even if posterior densities are not significantly different at the 95% credibility level (notably because of the higher a posteriori variance with model M 1 ), differences in posterior estimates reaches 13% for µ (relative difference between 0.34 with M 1 instead of 0.3 with M 0 ), and 70% for ξ (relative difference between 696 m 2 .s -1 with model M 1 and 410 m 2 .s -1 with model M 0 ). This indicates that using either M 0 or M 1 makes a huge difference if one aims at interpreting the value of posterior estimates, for, e.g., relating avalanche friction characteristics to snow and topographical conditions.Another remarkable result is the diminution of the error variance σ 2 from model M 0 to model M 1 (Table3.4). Indeed, the errors variance σ 2 is higher if the autocorrelation φ is not included. Specifically, the standard deviation of the errors decreases from 3.6 to 3.1 m.s -1 indicating that with model M 1 velocity predictions may be seen as 13% more accurate (relative difference between both estimates), another desirable property for practical use in snow science. Eventually, the autocorrelation parameter φ is largely positive, with a posterior mean of 0.71. Also, its posterior 95% credibility interval whose lover bound is 0.33 firmly excludes the zero value corresponding to M 0 . This pleads for a high and significant autocorrelation between velocities along flow propagation.

Figure 3 . 4 :

 34 Figure 3.4: Prior and posterior densities of model parameters. Panels a) to c) model M 0 and panels d) to g) model M 1 .

Figure 3 . 5 :

 35 Figure 3.5: Joint posterior distribution of parameters of M 0 model highlighting interparameter correlations.

Figure 3 . 6 :

 36 Figure 3.6: Joint posterior distribution of parameters of M 1 model highlighting interparameter correlations.

Figure 3

 3 Figure 3.7: Predictive velocity distributions versus data: a) model M 0 and b) model M 1 . 90% posterior credibility intervals (CI) are computed according to Equations (3.14) and (3.15) for the Voellmy propagation model (orange dotted lines), and the complete propagation and error model (blue dotted lines), respectively. Black plain line denotes the posterior median of the complete model. 90% prior credibility intervals are computed according to Eq. (3.16) (green dotted lines). Observations used for calibration appear as red points.

Figure 3 . 8 :

 38 Figure 3.8: Panels a) to d): Boxplots of the MAP estimators under models M 0 (red colour)and M 1 (blue colour) corresponding to the 100 synthetic avalanches. The true parameter values used for synthetic data generation are shown with a green colour line. Panel e) Boxplot of log 10 B 10 obtained for the 100 synthetic avalanches. The value of 2 in green corresponds to the reference value of[START_REF] Kass | Bayes Factors[END_REF] above which evidence in favour of M 1 is decisive.

Figure 3 . 9 :

 39 Figure 3.9: Prior sensitivity analysis. Postk, where k ranges for 1 to 6, are the posterior distributions obtained with the priors 1 to 6 of Table 3.6. Prior1 is the informative prior used in the paper core. Panels a) to c): M 0 model; panels d) to g): M 1 model.

Figure 3 .

 3 Figure 3.10: Images used in the photogrammetric process. a) and b) are the left and right images for the avalanche released on 13 February 2013 for t=20s after triggering. This couple is used for restitution and to map the location of the avalanche head. Markers plot the permanent ground control points used for image orientation. To illustrate the spatial extension, the coordinates and elevations (in meters) are indicated for the upper and lower control points (m3, t6), and for control point m5 at location of Col du Galibier road. The right image c) has been taken after the avalanche stops. Temporary control points t1 to t6 are setup after the avalanche in the runout area to improve image orientation in this area according to the avalanche deposit.

  1, . . . , n, where υ(b)k is a bootstrap sample of {υ k } 1≤k≤n . 5. B bootstrap samples of Y are created using Step 4. For each bootstrap sample (X, Y (b) ), we estimate the leave-one-out kernel smoother m (b)

Figure 4 . 1 :

 41 Figure 4.1: Scalar g-Sobol function: estimation of first-order Sobol' indices for X 1 , X 2 , X 3 and X 4 using four estimation methods: cnp, np, Solis2019 and Tissot2012. To draw the boxplots, the indices were computed in N = 100 independent samples of sizes n ∈ {100, 200, 300}. For each iteration of cnp and Solis2019, we used B = 100 bootstrap samples.

Figure 4 . 2 :

 42 Figure 4.2: Scalar g-Sobol function: estimation of first-order Sobol' indices for X 5 , X 6 , X 7and X 8 using four estimation methods: cnp, np, Solis2019 and Tissot2012. To draw the boxplots, the indices were computed in N = 100 independent samples of sizes n ∈ {100, 200, 300}. For each iteration of cnp and Solis2019, we used B = 100 bootstrap samples.

  Figure4.3 shows the estimation of first-order Sobol' indices using cnp, np, Solis2019 and Tissot2012 in N = 100 independent samples of sizes n ∈ {100, 200, 300}. The bias accuracy of cnp compared to the other three methods is the best (see Figure4.3). Furthermore, the variance of cnp estimation is lower than the other three methods (see Figure4.3). Additionally, for this test case, the bias correction performance is more evident than in the previous section (seeFigures 4.1 and 4.2 compared to Figure 4.3).

Figure 4 . 3 :

 43 Figure 4.3: Multivariate g-Sobol function: estimation of the aggregated Sobol' indices using cnp, np, Solis2019 and Tissot2012. Boxplots were drawn from N=100 independent samples of sizes n ∈ {100, 200, 300}. To use cnp and Solis2019, we set B = 100.

Figure 4 . 5 :

 45 Figure 4.5: Mass-spring model: Estimation of the aggregated Sobol' indices using np, Solis2019 and Tissot2012. Boxplots were drawn from N=100 independent samples of different sizes n ∈ {100, 200, 300}. To use cnp and Solis2019, we set B = 100.

Figure 4 . 6 :Figure 4 . 7 :

 4647 Figure 4.6: Avalanche model: Inputs and outputs. x start is fixed. All other quantities are taken as random and considered in the sensitivity analysis.

Figure 4 .

 4 Figure 4.7 shows the functional high density region (HDR) boxplots of 300 velocity (a) and snow depth curves (b) and the boxplot for the runout distance (c). The functional HDR plots are a tool for visualizing large amounts of functional data based on the estimation of the bivariate kernel density function of the two first components of the decomposition of the functional data [see Hyndman and Shang, 2010, for more details].The HDR boxplots show the 50% HDR and the 100% HDR in light and dark gray, respectively. The modal curve is plotted with a solid line. These Figures were obtained using the R package rainbow developed by[START_REF] Hyndman | Rainbow plots, bagplots, and boxplots for functional data[END_REF].

Figure 4 . 8 :

 48 Figure 4.8: Avalanche model: path's topography (black line) with T = 84 discretization points and the subset of the path in gray box (T 1 = 33 discretization points) where the sensitivity analysis is performed. The panel highlights the percentage of avalanche flowing as function of the position within the path.

Figure 4 . 9 :

 49 Figure 4.9: Avalanche model: a) scalar Sobol' indices for the velocity and b) flow depth outputs. The white line corresponds to the path's topography.

Figure 4 . 10 :

 410 Figure 4.10: Avalanche model: scalar Sobol' indices for the runout distance output.

Figure 4 . 11 :

 411 Figure 4.11: Avalanche model: aggregated indices estimated with two ways of aggregating the information. a) Two PCAs performed separately for the two functional outputs: one for the velocity output and the other one for the flow depth output. b) Simultaneous PCA to estimate the indices of the three outputs (the two functional ones and the scalar runout distance) altogether.

Figure 4 .

 4 13 shows the comparison between the non iterative version of[START_REF] Racine | Bias-Corrected Kernel Regression[END_REF] algorithm (cnp) and the iterative version (ite. cnp) in the g-Sobol test case.

Figure 4 Figure 4 . 13 :Figure 4 . 14 :

 4413414 Figure 4.12: Mass-spring model: Estimation of the aggregated Sobol' indices using cnp, np and Tissot2012. Boxplots were drawn from N=100 independent samples of sizes 200. The estimation results using the same sample (one sample) are two different samples (two samples) to estimate Sobol' indices and PCs are shown. We set B = 100 for cnp.

  ) d i=1 φ i = val({1, . . . , d}). 2. (Symmetry) If val(u∪{i}) = val(u∪{ }) for all u ⊆ {1, . . . , d}-{i, j}, then φ i = φ . 3. (Dummy) If val(u ∪ {i}) = val(u) for all u ⊆ {1, . . . , d}, then φ i = 0. 4. (Additivity) If val and val' have Shapley values φ and φ respectively, then the game with characteristic function val + val' has Shapley value φ i + φ i for i ∈ {1, . . . , d}.

Proposition 1 .

 1 The aggregated Shapley effects GSh i , i ∈ {1, . . . , d}, correspond to the Shapley value with characteristic function defined as: val(u) = p j=1 Var(Y j )val j (u) p j=1 Var(Y j ) , u ⊆ {1, . . . , d}. (5.5) Proof. The proof is straightforward. It is a direct consequence of the linearity and additivity properties of the Shapley value. Let i ∈ {1, . . . , d} and j ∈ {1, . . . , p}. The characteristic function val j (see (5.2)) has Shapley value Sh j i , i ∈ {1, . . . , d}. Thanks to the linearity and additivity properties (see properties 4. and 5. of the Shapley value), the characteristic function u ⊆ {1, . . . , d} → p j=1 Var(Y j )val j (u) p i=1 Var(Y j ) leads to the Shapley value p j=1 Var(Y j )Sh j i p i=1 Var(Y j ) • The characteristic function (5.5) can be written in matricial form:

ji

  and GSh i . for b = 1 to b = B do 1. Create a n bootstrap sample Y (b) by sampling with replacement from the rows of Y. 2. Compute, for 1 ≤ j ≤ p, σ 2,(b) j the empirical variance of Y (b) j . 3. For each j ∈ {1, . . . , p}:3.1. For all u and for all (s ) 1≤ ≤Nu compute E j u,s using (5.7). 3.2. For all u, create a N u bootstrap sample E j,(b) u,s by sampling with replacement from E j u,s 1≤ ≤Nu computed in Step 3

Figure 5 . 1 :Figure 5 . 2 :

 5152 Figure 5.1: Linear Gaussian model in dimension d = 2: MAE for the estimation of scalar Shapley effects of the output Y 1 estimated from N=300 i.i.d. samples as a function of N tot using different sample sizes a) n = 1000, b) n = 2000 and c) n = 5000.The 0.05 and 0.95 pointwise quantiles of the absolute error are drawn with gray polygons. The empirical coverage probability of the 90% bootstrap simultaneous confidence intervals is displayed with dotted lines. The theoretical coverage probability 0.9 is also shown with a plain gray line. The bootstrap sample size is fixed to B = 500.

Figure 5 . 3 :

 53 Figure 5.3: Mass-spring model: a) Explained variance as a function of the number of principal components. The gray line is displayed at 95% of the variance explained which corresponds to 6 components. The MAE computed from N = 100 independent realizations of the aggregated Shapley effect estimates, using the first 6 principal components, is drawn as a function of N tot with b) n = 1000, c) n = 2000 and d) n = 5000. The 0.05 and 0.95 pointwise quantiles of the MAE are drawn with gray polygons.The probability of coverage of the 90% bootstrap simultaneous intervals is displayed with a dotted line. The theoretical POC 0.9 is highlighted with a plain gray line. The bootstrap sample size is fixed to B = 500.

  in well documented in the French avalanche database[START_REF] Bourova | A new web-based system to improve the monitoring of snow avalanche hazard in France[END_REF]. The model outputs are the flow velocity, flow depth trajectories in the path D and runout distance of an avalanche. Note that the model has two functional Input Description Distribution µ Static friction coefficient U[0.05, 0.65] ξ Turbulent friction [m/s 2 ] U[400, 10000] l start Length of the release zone [m] U[5, 300] h start Mean snow depth in the release zone [m] U[0.05, 3] x start Release abscissa [m] U[0, 1600]

Figure 5 . 4 :

 54 Figure 5.4: The avalanche model.

Figure 5 . 5 :

 55 Figure 5.5: Avalanche model, scenario 1: scatter-plots of initial (black points) and constrained (gray points) samples. In the figure's diagonal, the density function of the initial (gray color) and constrained (transparent) samples are displayed. Input correlations of the original and constrained samples are shown. 1000 subsamples of original and constrained samples are used for illustration purpose.

Figure 5 . 6 :

 56 Figure 5.6: Avalanche model, scenario 1: a) and b) functional HDR boxplots of velocity and flow depth curves, resp. It is shown 50% HDR (dark gray), 100% HDR (light gray) and modal curve, the curve in the sample with the highest density (black line). c) runout distance boxplot. The constrained sample size is n 1 = 6152.

Figure 5 . 7 :

 57 Figure 5.7: Avalanche model, scenario 1: a) and b) ubiquitous Shapley effects of velocity and flow depth curves, resp. and, c) runout distance Shapley effects. Shapley effects are estimated with a sample of size 6152 and Ntot=2002. The local slope is displayed with a white line. A gray dotted rectangle box is displayed at interval [2017, 2412] where snow avalanche return periods vary from 10 to 10 000 years. The bootstrap sample size is fixed to B = 500.

Figure 5 . 8 :

 58 Figure 5.8: Avalanche model, scenario 1: a) and c) correlations between the first and second fPCs and the velocity. Scalar Shapley effects for b) fPC 1 and d) fPC 2 are estimated with n = 6152 and N tot = 2002. For the confidence intervals, B = 500 is used. The local slope is displayed with a gray line.

Figure 5 . 9 :

 59 Figure 5.9: Avalanche model, scenario 1: a) and c) correlations between the first and second fPCs and the flow depth output. Scalar Shapley effects for b) fPC 1 and d) fPC 2 are estimated with n = 6152 and N tot = 2002. For the confidence intervals, B = 500 is used. The local slope is displayed with a gray line.

Figure 5 . 10 :Figure 5 .

 5105 Figure 5.10: Avalanche model, scenario 1: a) and b) aggregated Shapley effects of velocity and flow depth curves calculated over space intervals [x, 2412m] wherex ∈ {1600m, 1700m, . . . , 2412m}. Shapley effects are estimated with samples of size 6152 and Ntot=2002. Effects are estimated using the first fPCs explaining more than 95% of the output variance. The local slope is displayed with a gray line. A gray dotted rectangle is displayed at[2017m, 2412m] where snow avalanche return periods vary from 10 to 10 000 years. The bootstrap sample size is fixed to B = 500.

Figure 5 . 12 :

 512 Figure 5.12: Avalanche model, scenario 2: a) and b) aggregated Shapley effects of velocity and flow depth curves calculated over space intervals [x, 2204] where x ∈ {1600, 1700, . . . , 2204} and using the first fPCs which have 95% of output variance. Shapley effects are estimated with samples of size 1284 and Ntot=800. The local slope is displayed with a gray line. A gray dotted rectangle is displayed at[2017m, 2204m] where return periods vary from 10 to 300 years. The bootstrap sample size is fixed to B = 500.

Figure 5 .

 5 Figure 5.13: Linear Gaussian model with d = 10 inputs: comparison between N u = N tot /(2 d -2) for all u (called N u uniform) and N u = N * u for all u (called N u optimal) with n = 10 000 and N tot = 54 000. The true scalar Shapley effects are displayed with black lines.

Figure 5 . 14 :

 514 Figure 5.14: Linear Gaussian model in dimension d = 2: boxplots of scalar Shapley effects and evolution of the estimation by varying the number of closest neighbors N I = {2, 3, 5, 10, 20, 30, 100, 300}. Boxplots are drawn with N = 100 independent random samples with n = 5000, N tot = 5000 and N * u = N tot d |u| -1 (d -1) -1 ,

Figure

  Figure 5.16 (resp., Figure 5.17) shows the correlations between the fPCs and the velocity (resp., flow depth) output for scenario 2.

Figure 5 .Figure 5 . 17 :

 5517 Figure 5.16: Avalanche model, scenario 2: a) and c) correlations between the first and second fPCs and the velocity output. Scalar Shapley effects for b) fPC 1 and d) fPC 2 are estimated with n = 1284 and N tot = 800. For the confidence intervals, B = 500 is used. The local slope is displayed with a gray line.

Figure A. 1 :

 1 Figure A.1: Ishigami function. Estimation of Sobol' indices using sobolroalhs, NSD and EASI. For each sample size n ∈ {100, 200, 500, 1000}, N = 100 samples are used. The true sensitivity indices values are shown with black lines.

  Figure A.1 shows the estimation of the Sobol' indices using NSD, EASI and sobolroalhs. The true sensitivity indices of the Ishigami function are known and shown with black lines. In order to see the evolution of the accuracy of estimation, the number of model evaluations n varies from 100 to 1000 . For each sample size n, N = 100 random samples are used. The number of classes required by NSD was set to √ n.

  Figure A.2 shows the estimation of Sobol' indices using the three methods for the Bartley et al. functions. Independently of the sample size, NSD estimation shows always a bias. EASI shows good results when the sample size increases.

Figure A. 2 :

 2 Figure A.2: Bartley et al. function. Estimation of Sobol' indices using sobolroalhs, NSD and EASI. For each sample size n ∈ {100, 200, 500, 1000}, N = 100 samples are used. The true sensitivity indices values are shown with black lines.

  Figure B.2: Scenario 1. Scatterplots of aggregated Shapley effects and local slope for a) flow velocity and b) flow depth.

  Figure B.3: Scenario 2. Scatterplots of ubiquitous Shapley effects and local slope for a) flow velocity and b) flow depth.

Figure B. 4 :

 4 Figure B.4: Scenario 2. Scatterplots of aggregated Shapley effects and local slope for a) flow velocity and b) flow depth.

Table 3

 3

.3: Parameters marginal prior distributions. The prior distribution of φ applies only to M 1 model. Γ represents the Gamma distribution and U the uniform distribution.

Table 3 .

 3 6: Marginal prior distributions used for the sensitivity analysis. Priors 3 and 5 are used with M 0 model only, and prior 6 with M 1 model only. Γ denotes the Gamma distribution and InvΓ denotes the inverse Gamma distribution. Prior1 is the one used in the paper core. "-" denotes vague marginal priors.

	Parameter	µ	ξ	σ 2
	Prior1	Γ(9, 0.03)	Γ(25, 40)	InvΓ(3, 20)
	Prior2	-	Γ(25, 40)	InvΓ(3, 20)
	Prior3	Γ(9, 0.03)	-	InvΓ(3, 20)
	Prior4	Γ(9, 0.03)	Γ(25, 40)	-
	Prior5	-	-	-
	Prior6	-	InvΓ(3, 20)	-

Table 4 .

 4 1: Mass spring model: Input description and uncertainty intervals.

Table 5 .

 5 2: Avalanche model, scenario 1: Input description and uncertainty intervals. In the the GSA, we consider vol start = l start × h start × 72.3/ cos(35 • ) instead of h start and l start .

  Shapley effects of velocity and flow depth curves calculated over space intervals[x, 2204] where x ∈ {1600, 1700, . . . , 2204} and using the first fPCs which have 95% of output variance. Shapley effects are estimated with samples of size 1284 and Ntot=800. The local slope is displayed with a gray line. A gray dotted rectangle is displayed at[2017m, 2204m] where return periods vary from 10 to 300 years. The bootstrap sample size is fixed to B = 500.

study is out of the scope of this paper. Recently, in the R package sensitivity the function sobolshap knn to estimate Shapley effects with n and N tot from a given data sample has been implemented. This function uses a tree based technique to approximate nearest-neighbor search which reduces drastically computation times. The function is particularly attractive if n and N tot are high, we could even use N tot = (2 d -2) × n.

  -1 ] Flow depth [m] Flow velocity [m s -1 ] Flow depth [m]

	µ	-0.79	-0.28	-0.93	-0.70
	ξ	-0.50	0.54	-0.20	-0.75
	x start	0.92	-0.76	0.89	-0.92
	vol start -0.92	0.49	-0.87	0.82

Table B.1: Scenario 1. Correlations between Shapley effects and local slope.

  -1 ] Flow depth [m] Flow velocity [m s -1 ] Flow depth [m]Table B.2: Scenario 2. Correlations between Shapley effects and local slope.

	µ	-0.89	-0.67	-0.96	-0.93
	x start	0.64	-0.76	0.95	0.59
	vol start 0.87	-0.42	0.83	0.54

The runout distance is the position where the avalanche stops flowing.

Avalanche occurred in Switzerland the 27 January 1968. The event destroyed a farmhouse and killed four people.

This definition is, in general, widely accepted in the SA community.

If the number of inputs is high (e.g., GPT-3 autoregressive model which has 175 billion of parameters), it might be better to use Morris method first.

Remerciements

31.25+87.5h start µ|h start , x nstart N (0.449 -0.013x nstart + 0.025h start , 0.11 2 ) Table 5.3: Avalanche model: Scenario 2. Input description and uncertainty intervals. x nstart is a normalization of x start . A linear relationship between h start and l start inferred from the local data is used [START_REF] Eckert | Long-term avalanche hazard assessment with a Bayesian depth averaged propagation model[END_REF]. Gamma distribution parameters are its shape and rate. In the the GSA, we consider vol start = l start × h start × 72.3/ cos(35 • ) instead of h start and l start .

Scenario 2 Principle

The aim is now to determine the most influential inputs in the context of strong knowledge regarding input distributions. In [START_REF] Eckert | Long-term avalanche hazard assessment with a Bayesian depth averaged propagation model[END_REF], the authors developed a Bayesian framework to estimate input distributions from available avalanche observations. The objective is long-term avalanche hazard assessment in order to assess the related risk for buildings and people inside. In the avalanche literature, it is assumed that ξ depends on the path topography, so it is a parameter and not a variable varying from one avalanche to another in [START_REF] Eckert | Long-term avalanche hazard assessment with a Bayesian depth averaged propagation model[END_REF]]'s model. The input ξ is therefore fixed to its posterior estimate, 1300. Other input variables in this scenario are dependent. The dependence between h start and l start is modeled with a linear function l start = 31.25 + 87.5h start , and similarly as in scenario 1, we consider vol start as input of the GSA instead of h start and l start .

The complete input distribution resulting from the Bayesian inference on the studied path is described in Table 5.3. Input correlations have been computed. As an example, the correlation between µ and vol start is 0.8. As a preliminary step to GSA for this scenario, we apply the following constraints:

(i) avalanche is flowing in the interval [1600m, 2204m] where snow avalanche return periods vary from 10 to 300 years, (ii) avalanche volume is superior to 7 000 m 3 and, (iii) µ coefficient is inferior to 0.39 as we focus on dry snow avalanches.

Under these conditions, we sample the full set of dry snow avalanches that could cause strong material or human damages on the studied site.

Global sensitivity analysis results

We first ran n 0 = 100 000 avalanches from an i.i.d. sample of input distribution following Table 5.3. After applying the constraints, the sample size was reduced to n 2 = 1284 and the input distribution was modified. For example, µ and vol start correlation changes from 0.8 to 0.2 which is still non negligible. Ubiquitous Shapley effects are displayed on Subsection 5.6.3 panels a and b. For the velocity, the three inputs have a similar importance till 1900m, then vol start importance decreases and µ and x start importance increases (see Subsection 5.6.3 a). Similarly, as in scenario 1, the effects show fluctuations 5.4: Input dimension d = 10: N u values for both strategies (i) and (ii) for n = 10 000 (number of simulations) and N tot = 54 000 (total cost).

Comparison between the uniform allocation

) , and the one introduced in [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF]]

, d}

The total cost N tot to estimate all the c(u), ∅ u {1, . . . , d} could be uniformly divided as

To test if there is a difference in the estimation by using N u uniform or N * u optimal, we consider the Gaussian linear model example from [START_REF] Owen | On Shapley value for measuring importance of dependent inputs[END_REF]:

where X i are independent centered Gaussian variables with standard deviation σ i = i 2 and, β i = 1 for all i ∈ {1, . . . , d}. The scalar Shapley effects can be computed analytically (see e.g., [START_REF] Owen | On Shapley value for measuring importance of dependent inputs[END_REF]):

We used N = 100 independent samples of size n = 10 000 and N tot = 54 000 to draw the boxplots shown in Figure 5.13. Table 5.4 summarizes the sample size used to estimate each conditional element E(Var(Y |X -u )) in function of |u| in dimension d = 10. Figure 5.13 shows that in high dimension, at least for that toy model, the optimal allocation N * u leads to a better accuracy than the uniform allocation. When dimension is low (for example d = 4 as in the mass-spring model), taking N * u or N u give similar results because the two options are quite similar.

Comparison between different values of N I

In Figure 5.14, we analyzed the impact of using different numbers of neighbors N I in the Shapley effect estimation for the linear Gaussian model in dimension d = 2 test case described in Subsection 5.5.1. For the estimation we used n = 5000, N tot = 5000 and

Parameter N I is varied from 2 to 300.

We see in Figure 5.14 that there is no clear evidence if N I ≥ 3 is always the best choice. However in the following we choose N I = 3 following the recommendations in [START_REF] Song | Shapley effects for global sensitivity analysis: Theory and computation[END_REF], where the authors prove in the non given data framework that choosing N I = 3 is better than choosing N I = 2 and that the gain from increasing N I from 3 decreases as N I increases.

Appendices

A

Given data estimation methods for Sobol' indices

A.1 Given data method [START_REF] Plischke | Global sensitivity measures from given data[END_REF] [see also Borgonovo et al., 2016] proposed a given data method which relies on the notion of class-conditional densities, where a class is a sub-sample stemming from a suitable partition of the data. As its name suggests, the method does not require any specific sampling design to estimate the Sobol' indices (or other sensitivity measures).

Hereafter, it is denoted by NSD (non sampling design).

Let us consider a scalar output Y = f (X). Let (x j 1 , . . . , x j d , y j ) j∈{1,...,n} be a n sample of (X, Y ). The support of the input X i , denoted by

For notation simplicity, we remove i from the classes C but they depend on X i . The estimation of S i , the first-order scalar Sobol' index, proposed by [START_REF] Plischke | Global sensitivity measures from given data[END_REF] is:

where n counts the number of realizations of X i in C , ȳ and s 2 y denotes the empirical mean and variance of Y , resp., ȳ represents the within sample class mean and it is used to estimate the conditional mean of

The main drawback of this estimation method is the bias induced by the choice of the partition [START_REF] Antoniano-Villalobos | Nonparametric estimation of probabilistic sensitivity measures[END_REF]. Moreover, in [Antoniano-Villalobos et al., 2019, Appendix A.2] different values of L were analyzed but till this moment, there is no universal valid rule to choose L. [Theorem 1 of Borgonovo et al., 2016] states that the estimator is consistent if the number of classes can be written as L = g(n) where g is a monotonically increasing function such that:

For the test cases shown in Section A.3, we take L = √ n.

B

Correlations between Shapley effects and local slope

In Chapter 5, we found that the correlations between Shapley effects and local slope were significant. The local slope term assigns the angle between the abscissa and the corridor. This fact might suggest that local slope changes play a role in the input importance. B.1 and Table B.2 for scenarios 1 and 2, respectively. An interesting finding is that the correlations between µ and the local slope are always negative, this could suggest that if the local slope is high the importance of µ is low. However, to generalize the results, Shapley effects should be estimated in other corridors. It would also be interesting to include the topography in the global sensitivity analysis.

q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q -10 0 10 20 30 40 0.0 0.2 0.4 0.6 0.8

Local slope [degrees] Sh q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q µ ξ x start vol start a) q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -10 0 10 20 30 40 Local slope [degrees] Sh q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q b)