# Impact of input forcings variability on the global sensitivity analysis of a hydrological model

Katarina Radišić, Claire Lauvernet, Arthur Vidard

► To cite this version:

Katarina Radišić, Claire Lauvernet, Arthur Vidard. Impact of input forcings variability on the global sensitivity analysis of a hydrological model. Environmental Modelling and Software, 2025, 192, pp.106522. 10.1016/j.envsoft.2025.106522. hal-04849856v2

HAL Id: hal-04849856

https://hal.inrae.fr/hal-04849856v2

Submitted on 15 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Impact of input forcings variability on the global sensitivity analysis of a hydrological model

Katarina Radišić<sup>a,b ID,*</sup>, Claire Lauvernet<sup>a ID</sup>, Arthur Vidard<sup>b ID</sup>

## ARTICLE INFO

Dataset link: doi.org/10.57745/C455R6, https://www.uqlab.com, doi.org/10.57745/4MDYN U, doi.org/10.15454/2HAU8R

**Keywords:**

Global sensitivity analysis

Hydrology

Sobol' indices

Polynomial chaos expansion

Forcing variability

Uncertainty

## ABSTRACT

Hydrological models are subject to natural variability in forcing conditions, yet traditional global sensitivity analysis (GSA) often overlooks these uncertainties, potentially limiting the validity of results to specific conditions. To address this, we propose an approach that accounts for the uncontrollable nature of rain forcing variability. We treat the model parameters' Sobol' indices as random variables dependent on forcing variability and use polynomial chaos expansion metamodels to reduce the computational cost of their estimation. We apply this methodology to study soil moisture sensitivity in the physically-based distributed hydrological model PESHMELBA. Our results show that parameter rankings vary with forcing conditions. To consolidate these diverging GSA results, we propose a unique ranking based on aggregated sensitivity indices that accounts for parameter contributions across the entire domain of forcing conditions. This approach enhances the robustness of GSA to natural variability in forcings, thereby improving the reliability of subsequent GSA-based decisions.

## 1. Introduction

Global sensitivity analysis (GSA) provides information on the importance of model parameters for the simulation of model outputs. As GSA identifies the parameters having a larger effect on model outputs, it is a useful tool in all stages of model development and validation (Tarantola et al., 2024). GSA is increasingly used in environmental sciences (Tang et al., 2007; Nossent et al., 2011; Garcia et al., 2019; Perrin et al., 2021; Heredia et al., 2021; Alipour et al., 2022; Jamous et al., 2023; Colombi et al., 2024), hydrology (Song et al., 2015; Nagel et al., 2020) and water quality (Fox et al., 2010; Lauvernet and Muñoz-Carpena, 2018; Hong and Purucker, 2018; Gatel et al., 2019; D'Andrea et al., 2020; Faúndez Urbina et al., 2020).

For models with large input sets, as it often is the case in hydrology, GSA can help in choosing a smaller, and thus more manageable size of parameters to be considered at different stages of model development (Wagener and Pianosi, 2019). This can be done by means of screening, that is classifying the parameters in two groups: influential and non-influential, or by ranking the parameters in decreasing order of their influence, then choosing to focus on the parameters appearing at the top of the ranking list. Screening is often performed with the Morris method (Morris, 1991), or more recently with the Hilbert-Schmidt Independence Criterion HSIC independence tests (De Lozzo and Marrel, 2016), ranking is typically performed with variance-based decomposition (Sobol') methods. Additionally, a good practice for removing

non-identifiable parameters at the calibration stage, is to perform a sensitivity analysis on the model's performance metrics (Mai, 2023).

The inputs of most deterministic hydrological models can be separated in two sets:

1. *operational parameters (x):* parameters whose values are fitted to obtain model outputs closer to observations, such as hydraulic conductivity, Manning coefficient, Van Genuchten parameters, etc.

2. *external forcings (z):* forcing terms such as rainfall, evapotranspiration, boundary conditions (piezometric levels or upstream/downstream river stages), or the dates of pesticide usage on the cultivated lands.

The uncertainty of the output of the (deterministic) hydrological model is thus due to the uncertainties of both types of inputs. However, the uncertainty of *x* and that of *z* differ in their nature. Indeed, while the uncertainty of *z* represents the natural fluctuations of the different forcings under which the hydrological model is expected to operate, the uncertainty of *x* is due to a misspecification of the parameters.

In a context of model development, the uncertainty of *x* should be reduced, for example, by obtaining more direct measurements when possible, or, when this is not the case, estimating parameter values through calibration or data assimilation processes. On the other hand, the uncertainty of *z* cannot be reduced, as it represents the uncertainty

<footnote><sup>a</sup> INRAE, RiverLy, 5 rue de la Doua CS 20244, 69625 Villeurbanne CEDEX, France</footnote>

<footnote><sup>b</sup> Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France</footnote>

<footnote>* Corresponding author at: INRAE, RiverLy, 5 rue de la Doua CS 20244, 69625 Villeurbanne CEDEX, France. E-mail address: katarina.radisic@inrae.fr (K. Radišić).</footnote>
of the natural conditions under which the hydrological model operates.

One of the goals of GSA is to help choose the inputs whose uncertainty should be reduced with highest priority. This is usually done by performing a GSA for a fixed set of external forcings z. However, ignoring the stochastic nature of z amounts to the GSA results of parameters x being valid only locally in the forcing conditions used. Indeed, due to the interactions of hydrological model inputs, there is no guarantee that the results of a GSA obtained under one set of external forcings can be extrapolated to other values of the external forcings, even if performed on the same catchment and similar conditions.

Furthermore, model forcings such as rainfall are typically multidimensional quantities, having a temporal or spatio-temporal dimension. To avoid making hypotheses on the forcing quantities' structure, there is a need for a non-intrusive approach.

To make the stochastic nature of z explicit, without making intrusive assumptions on its structure, it can be represented as depending on an abstract random event $\omega \subset \Omega$ defined on a probability space $(\Omega, \mathcal{F}, P)$. The probability space represents the variability of forcings under which the model is expected to operate, with $\Omega$ being the sample space, $\mathcal{F}$ the event space and $P$ the probability measure. The deterministic value of z is thus considered as a realization of $Z(\omega)$, and the (deterministic) hydrological model can be rewritten as a stochastic model $\mathcal{M}_s$ such that:

$$ \mathcal{M}_s : D \times \Omega \rightarrow \mathbb{R}^q, \\ (x, \omega) \mapsto \mathcal{M}_s(x, Z(\omega)), \qquad (1) $$

where the operational parameters x are represented with a p-dimensional vector with values in $D \subset \mathbb{R}^p$.

While many methods have been successfully developed for the sensitivity analysis of deterministic models (Da Veiga et al., 2021; Saltelli et al., 2008), not as much attention has been given to their stochastic counterparts. One of the most popular approaches in GSA is the variance-based sensitivity analysis with its Sobol' indices (Sobol, 2001). Zhu and Sudret (2021) list three extensions of Sobol' indices to stochastic simulators. Firstly, Iooss and Ribatet (2009) propose a direct extension of the input x with the stochastic inputs Z, thus treating them as additional input variables, thus not distinguishing the different natures of the two sources of uncertainty. A second choice, is that of studying the sensitivity of a summarizing statistical quantity of interest (QoI) of the model output distribution. Typical QoIs are: the output's mean value, its variance, higher moments (Dell'Oca et al., 2017), a quantile or differential entropy (Azzi, 2020). Hence, this second extension amounts to the classical definition of Sobol' index being applied on the QoI.

A third extension of Sobol' indices to stochastic models comes from considering the Sobol' indices as random variables. In this case, the randomness of the Sobol' indices comes from seeing them as functions of the stochastic inputs Z. Such approach was considered when studying the moments of Sobol' indices (Hart et al., 2017; Jimenez et al., 2017). Recently, random Sobol' indices and their robustness was studied in an application to stochastic agent-based models, (Carmona-Cabrero et al., 2024). This extension embodies the operational approach to GSA, introduced in Dell'Oca (2023), which argues for a distinction between the operational parameters (considered controllable by the modeler) and the (uncontrollable) stochastic inputs.

In this work, we perform a global sensitivity analysis of a hydrological model output by distinguishing the uncertainties of an uncontrollable external forcing (rain) and the uncertainties of controllable operational parameters (soil hydrodynamic properties). The outputs of interest are soil moisture profiles simulated by the PESHMELBA hydrological model (Rouzies et al., 2019). We obtain the distributions of the Sobol' indices, which are considered random variables, and study the robustness of the ranking under different realizations of the rain in a same event-based context. We then propose a unique ranking based on the aggregated sensitivity indices ASI, which describe the

part of the total output uncertainty (coming from both controllable and uncontrollable sources).

Section 2 presents briefly the Sobol' indices, their efficient calculation with polynomial chaos expansion, the notion of Sobol' indices as random variables (trajectory-based Sobol' indices) and the aggregated sensitivity indices (ASI). Section 3 presents the case study: the setup of the PESHMELBA model, the six operational parameters whose sensitivity indices are of interest and the rain measurements considered as realizations of the stochastic external forcing. Section 4 presents the results: the obtained distributions of the Sobol' indices for each operational parameter and shows the effect of the rain on the parameter ranking.

Section 5 discusses the type of model in which one can expect the stochastic input to impact the sensitivity analysis results, the computational cost of such an approach and other potential sources of stochastic inputs in hydrological models.

## 2. Methodology

In Section 2.1 the definition of Sobol' indices are recalled, then in Section 2.2 their estimation with polynomial chaos expansion method is explained. Section 2.3 introduces the trajectory-based Sobol' indices, that is, the notion of Sobol' indices being random variables dependent on the realization of the external forcing. Section 2.4 presents the aggregated sensitivity indices, used in this work to synthesize the information obtained from Sobol' indices in different external forcings.

### 2.1. Sobol' indices

We recall in this section the definition and basic notions of classical Sobol' indices. Let the function of interest $f_d$ be a deterministic function with a scalar output:

$$ f_d : D \rightarrow \mathbb{R}, \\ X \mapsto Y = f_d(X), \qquad (2) $$

where $X = (X_1, X_2, ..., X_p) \in D \subset \mathbb{R}^p$ are p independent scalar inputs. Sobol' indices are a variance-based GSA approach, meaning that they consider variance as a metric to quantify the contribution of each uncertain parameter to the uncertainty of the output. The classical Sobol' indices quantify the part of the output variance explained by the variations of each input parameter, or their interaction. The first and total order Sobol' indices of the input $X_i$, $i \in \{1, ..., p\}$, are defined respectively as:

$$ S_i = \frac{\mathrm{Var}[\mathbb{E}[Y | X_i]]}{\mathrm{Var}[Y]}, \qquad S_{T_i} = 1 - \frac{\mathrm{Var}[\mathbb{E}[Y | X_{i-1}]]}{\mathrm{Var}[Y]}, \qquad (3) $$

where $\sim i$ denotes the set of all inputs, but the input $X_i$, (Saltelli et al., 2008; Homma and Saltelli, 1996). The Sobol' indices lie in the interval $[0, 1]$, and a larger Sobol' value indicates a greater importance of the input variable or the group of input variables related to this index. The input variables can then be ranked by their order of importance with respect to their first or total Sobol' index. The first order index $S_i$ gives the information about the impact that the input $X_i$ has on the output when varied alone, while the total index $S_{T_i}$ gives the importance of $X_i$ and all of its interactions with other variables.

Theoretically, the Sobol' indices can be obtained for any function $f_d$. Practically, the estimation of these indices can be obtained by Monte Carlo simulations, on a large experimental design (minimum $M(p+1)$, with $M > 1000$ and $p$ the number of inputs), that can be very expensive in terms of evaluations of the function $f_d$, (Da Veiga et al., 2021; Iooss and Lemaître, 2015; Saltelli et al., 2004; Sobol, 1993). The number of necessary evaluations can be reduced either by sampling more intricate experimental designs, (Tissot and Prieur, 2015; Saltelli et al., 1999, 2008), or through the use of metamodels, such as polynomial chaos expansion (PCE), which is the approach used in this study.
## 2.2. Polynomial chaos expansion

Let $\mathbf{X}$ be a random vector with values in $\mathcal{D} \subset \mathbb{R}^p$, with finite variance and joint probability density function $g_X$ satisfying certain conditions (Xiu and Karniadakis, 2002; Ernst et al., 2012). Consider a deterministic function $f_d$ mapping a set of input parameters $\mathbf{X} \in \mathcal{D} \subset \mathbb{R}^p$ to the output $\mathbf{y} \in \mathbb{R}$. If $\mathbf{Y} = f_d(\mathbf{X})$ has finite variance, the function $f_d$ admits a polynomial chaos expansion (PCE):

$$ \mathbf{Y} = f_d(\mathbf{X}) = \sum_{\alpha \in \mathcal{N}_p} c_{\alpha} \psi_{\alpha}(\mathbf{X}), \quad (4) $$

where $p$ is the number of input parameters, $\{\psi_{\alpha}\}_{\alpha \in \mathcal{N}_p}$ is a basis of multivariate orthonormal polynomials (chosen according to the marginal probability density functions $g_X$ of the input parameters), $\alpha \in \mathcal{N}_p$ are multi-indices, i.e. $\alpha = (\alpha_1, \alpha_2, \dots, \alpha_p)$ where each $\alpha_i \in \mathbb{N}$ corresponds to the partial degree with which the $i$th input parameter $X_i$ is represented in the basis component $\psi_{\alpha}$ and $c_{\alpha} \in \mathbb{R}^{\mathcal{N}_p}$ are the PCE coefficients.

Once the function $f_d$ is rewritten via its polynomial chaos expansion, the Sobol' indices are obtained analytically from the coefficients $c_{\alpha}$ (Sudret, 2008). For example, to obtain the first order Sobol' index of the input parameter $X_i$, we consider only the coefficients $c_{\alpha}$ that correspond to basis polynomials $\psi_{\alpha}$ which only depend on the parameter $X_i$. In other words, consider only the coefficients $c_{\alpha}$ corresponding to basis polynomials $\psi_{\alpha}$ where the partial degrees of all input parameters other than $X_i$ are equal to zero, i.e. $\alpha_{j \neq i} = 0$. The first order Sobol' index of parameter $X_i$ is then given by the sum of squares of the concerned coefficients, normalized by the total variance of the output:

$$ S_i = \sum_{\alpha \in \mathcal{N}_i} c_{\alpha}^2 / D, \quad I_i = \{\alpha \in \mathcal{N}_p : \alpha_i > 0, \alpha_{j \neq i} = 0\}, \quad (5) $$

$$ D = \mathrm{Var} \left[ \sum_{\alpha \in \mathcal{N}_p} c_{\alpha} \psi_{\alpha}(\mathbf{X}) \right] = \sum_{\substack{\alpha \in \mathcal{N}_p \\ \alpha \neq \emptyset}} c_{\alpha}^2 \quad (6) $$

The total Sobol' index of input parameter $X_i$ is obtained similarly. However, the filtering criteria on $c_{\alpha}$ is changed, since here all interactions with other parameters are of interest. Thus, all appearances of a non-null partial degree $\alpha_i > 0$ are of interest:

$$ S_{T_i} = \sum_{\alpha \in \mathcal{N}_i} c_{\alpha}^2 / D, \quad I_{T_i} = \{\alpha \in \mathcal{N}_p : \alpha_i > 0\}, \quad (7) $$

where $D$ is as in Eq. (6)

As the expansion cannot use infinite terms, a truncated basis $\mathcal{A}$ must be chosen. This can be done efficiently by taking advantage of the sparsity-of-effects, a heuristic stating that most models describing physical phenomena are dominated by main effects and interactions of low order. This translates to PCE having coefficients $c_{\alpha}$ that decay quickly as the partial degree increases, and most information being contained in $c_{\alpha}$ with low values of $\alpha_i$, see Lüthen et al. (2021) for an overview. Thus, higher order interactions can be left out in the PCE. With the basis truncation, the previous equations become approximations:

$$ \begin{align} Y &= f_d(\mathbf{X}) \approx f_{\mathrm{PCE}}(\mathbf{X}) = \sum_{\alpha \in \mathcal{A}} c_{\alpha} \psi_{\alpha}(\mathbf{X}), \nonumber \\ \hat{S}_i &= \sum_{\alpha \in \mathcal{A}_i} c_{\alpha}^2 / \hat{D}, \quad \mathcal{A}_i = \{\alpha \in \mathcal{A} : \alpha_i > 0, \alpha_{j \neq i} = 0\}, \nonumber \\ \hat{S}_{T_i} &= \sum_{\alpha \in \mathcal{A}_{T_i}} c_{\alpha}^2 / \hat{D}, \quad \mathcal{A}_{T_i} = \{\alpha \in \mathcal{A} : \alpha_i > 0\}, \tag{8} \end{align} $$

$$ \hat{D} = \mathrm{Var} \left[ \sum_{\alpha \in \mathcal{A}} c_{\alpha} \psi_{\alpha}(\mathbf{X}) \right] = \sum_{\substack{\alpha \in \mathcal{A} \\ \alpha \neq \emptyset}} c_{\alpha}^2 $$

The quality of the metamodel $f_{\mathrm{PCE}}$ is evaluated by comparing it with a test set. This is quantified through the coefficient $Q^2$ which measures the quality of the prediction of a linear regression:

$$ Q^2 = 1 - \frac{\sum_{n=1}^{N_{\mathrm{test}}} (f(\mathbf{x}_n) - f_{\mathrm{PCE}}(\mathbf{x}_n))^2}{\sum_{n=1}^{N_{\mathrm{test}}} (f(\mathbf{x}_n) - \bar{f})^2}, \quad (9) $$

where $f(\mathbf{x}_n)$ and $f_{\mathrm{PCE}}(\mathbf{x}_n)$ are, respectively, the evaluations of the original function and of the PCE metamodel in the points $\{\mathbf{x}_n\}_{n=1}^{N_{\mathrm{test}}}$, and $\bar{f}$ is the mean of the evaluations of the original function, $\bar{f} = \sum_{n=1}^{N_{\mathrm{test}}} f(\mathbf{x}_n)$. The closer the $Q^2$ is to 1, the better is $f_{\mathrm{PCE}}$ in approximating $f$.

With sparse PCE, the number of model runs for Sobol' index calculation can be drastically reduced, to 10$\rho$ or less model runs, depending on the model complexity. The precision of the estimated Sobol' indices can be evaluated through bootstrap confidence intervals (Dubreuil et al., 2014). They are obtained by creating multiple resamples with replacement of the original experimental design and calculating the Sobol' indices from each resample.

## 2.3. Trajectory-based Sobol' indices

Let $f_s$ be a stochastic model with a scalar output:

$$ \begin{align} f_s : D \times \Omega &\rightarrow \mathbb{R}, \nonumber \\ (x, \omega) &\mapsto f_s(x, Z(\omega)), \tag{10} \end{align} $$

where $x$ is the input vector belonging to the input space $D \subset \mathbb{R}^p$ and $\Omega$ denotes the probability space representing the intrinsic stochasticity of $\omega$. Note that, for a fixed value of the parameters $x = x_0 \in \mathbb{R}^p$, the model reduces to $f_s(x_0, \cdot) : \Omega \rightarrow \mathbb{R}$. Thus, for a fixed $x_0$, the output is a random variable, whose randomness stems from that of $\omega$.

On the other hand, by fixing a random realization of the abstract event $\omega = \omega_0$, the output is not stochastic anymore, but reduces to a deterministic function in the input parameters $f_s(\cdot, Z(\omega_0)) : D \rightarrow \mathbb{R}$. In other terms, $f_s$ can be considered a random field, indexed by the input parameters $x \in D$. The random field, for a fixed realization of $\omega$, becomes a deterministic function of $x$. These deterministic functions, each corresponding to a realization of $\omega$, are called trajectories in some stochastic metamodeling literature, (Zhu and Sudret, 2021; Lüthen et al., 2023) and the same term is used in this paper.

By observing one particular trajectory, the stochasticity of the model is removed, and classical Sobol' indices (Eq. (3)) can be evaluated on the fixed trajectory. By iterating this approach on multiple trajectories, the classical Sobol' indices become a function of $\omega$, and thus random variables. This extension of Sobol' indices for stochastic models is called trajectory-based Sobol' indices:

$$ S_i(\omega) = \frac{\mathrm{Var}_{X_i} [\mathbb{E}[Y | X_i, \omega]]}{\mathrm{Var}[Y | \omega]}, \quad S_{T_i}(\omega) = 1 - \frac{\mathrm{Var}_{X_{-i}} [\mathbb{E}[Y | X_{-i}, \omega]]}{\mathrm{Var}[Y | \omega]} \quad (11) $$

The statistical properties of trajectory-based indices can be studied at the cost of repeating a standard Sobol' analysis for different realizations of the stochastic inputs. This assumes that it is possible for the modeler to perform replications with the same internal randomness, i.e. it is possible to obtain multiple simulations of the stochastic model with the same realization of the abstract event $\omega = \omega_0$ and varying values of $x$.

However, this also means that the number of model simulation needed for the estimation of trajectory-based Sobol' indices is multiplied by the number of realization of the abstract event. The need for an efficient means of computing the Sobol' indices is accentuated. The general framework for an efficient approximation is presented in Algorithm 1.

The algorithm returns $R$ realizations of Sobol' indices of $f_{\mathrm{PCE}}^{(r)}(\mathbf{X})$, denoted with $\hat{S}_i^{(r)}$ and $\hat{S}_{T_i}^{(r)}$ for the first order and total order indices, $r \in \{1..R\}$, where $R$ is the number of samples over $\Omega$. The $m$th moment of the Sobol' indices can be estimated with the following estimator:

$$ \hat{\mu}_{S_i}^{[m]}(\omega) = \frac{1}{R} \sum_{r=1}^{R} (\hat{S}_i^{(r)})^m, \quad (12) $$

The estimator bias and variance are derived by Hart et al. (2017). For the special case where $m = 1$, an estimate of the error can be
for r = 1 to R do

Take the r^th realization of ω, leading to z_r = Z(ω_r)

Generate an experimental design {X_j^{(r)}}_{j=1}^{N_{train}}
in the parameter space

Simulate $f_s(\mathbf{X}_j^{(r)}, \mathbf{z}_r)$, for each $j = 1, \dots, n$

Construct metamodel $f_{PCE}^{(r)}(\mathbf{X}) \approx$

$f_s(\mathbf{X}, \mathbf{z}_r)$ using data set $\{(\mathbf{X}_j^{(r)}, f_s(\mathbf{X}_j^{(r)}, \mathbf{z}_r))\}_{j=1}^{N_{train}}$

Compute $\hat{S}_i^{(r)}$ and $\hat{S}_i^{(r)}$ using the metamodel

$f_{PCE}^{(r)}(\mathbf{X})$ for $i \in \{1, 2, \dots, p\}$

end

**Algorithm 1:** Estimation of the trajectory-based Sobol' indices. For each realization $\omega_r$ an approximation of the Sobol' indices is performed with a PCE metamodel. The possibility of having different experimental designs for each trajectory $r \in \{1..R\}$ is indicated by the dependence of the experimental design $\{X_j^{(r)}\}_{j=1}^{N_{train}}$ on $r$. In this study, the same experimental design is used for all trajectories, thus the index is dropped in the following.

obtained for the point estimator of the expected value of $S_i(\omega)$ given by the sample mean:

$$ \mathbb{E}_{\omega} \left\{ \left( \hat{\mu}_{i}^{[1]} - \mathbb{E}_{\omega} \{ S_{i} \} \right)^{2} \right\} \leq \mathbb{E}_{\omega} \{ \hat{S}_{i} - S_{i} \}^{2} + \frac{1}{4R} \quad (13) $$

Thus, the estimator bias comes exclusively from the approximation error due to the metamodel, while the estimator variance is bounded by $\frac{1}{4R}$.

## 2.4. Aggregated sensitivity indices

For models with multiple outputs, (for example spatialized or dynamic outputs in distributed hydrological models), the information provided by Sobol' indices on individual outputs can be aggregated into scalar values by weighting the indices with the corresponding estimated variances of each output (Radišić et al., 2023). Such aggregated Sobol indices (ASI) are proposed in Gamboa et al. (2014) and are defined as:

$$ ASI_i = \frac{\sum_{r=1}^{R} \mathrm{Var}(Y^{(r)}) S_i^{(r)}}{\sum_{r=1}^{R} \mathrm{Var}(Y^{(r)})}, \quad (14) $$

where $Y^{(r)}$ is the rth component of the R-dimensional output and $S_i^{(r)}$ is the sensitivity index of the set of input parameters $i$ on the rth component of the output. Thus, by considering each trajectory $f_d^{(r)} = f_s^{(r)}(\cdot, \omega_r)$ as a different model output, the ASI indices become:

$$ ASI_i = \frac{\sum_{r=1}^{R} \mathrm{Var}(f_d^{(r)}(\mathbf{X})) S_i^{(r)}}{\sum_{r=1}^{R} \mathrm{Var}(f_d^{(r)}(\mathbf{X}))}, \quad (15) $$

In the special case where the variances conditioned on these realizations are equal, $\mathrm{Var}(f_d^{(r)}(\mathbf{X})) = \mathrm{const.}$ for each $r \in \{1..R\}$, the estimator of the expected value of the trajectory-based Sobol' indices (Eq. (12)) and the ASI (Eq. (15)) coincide. Other ways for considering multivariate outputs are proposed in Li et al. (2016).

## 3. Case study

In this section, the PESHMELBA model is introduced along with its representation of the soil moisture profiles, which are the output of interest in this study. The uncertainty in soil moisture profiles is considered to come from two sources: the model parameters describing the soil hydrodynamic properties and the natural variability of the rainfall events under which the model is expected to operate. Here, the soil hydrodynamic properties are considered controllable, whereas the natural variability of the rain events is considered as the uncontrollable source of uncertainty.

### 3.1. PESHMELBA model setup

PESHMELBA (Rouzies et al., 2019) is a process-oriented, physically-based model representing the spatial organization of an agricultural landscape. It simulates the fate of pesticides in small agricultural catchments, transferred by water and contaminant surface and subsurface fluxes. PESHMELBA has a modular structure, simulating the physical processes of each soil compartment, then coupling them via the Open-PALM coupler (Buis et al., 2006). The water transfers occur via three types of pathways: infiltration via the Richards' equation solved by Ross (2003), surface runoff (kinematic wave equation) and lateral exchanges (Darcy's law) occurring in subsurface saturated zones. The model setup in this case study is based on a virtual version of the Morcille catchment in Beaujoulais, France (Gouy et al., 2021) and is available on Rouzies and Lauvernnet (2022).

The water and pesticide transfers between the spatial elements are driven by the changes of capillary pressure is the soil, which is directly tied to the soil moisture observed in the columns through Van Genuchten equations. The moisture profile, i.e. the distribution of the soil volumetric water content along a vertical soil section, depends (at a given time) on previous weather conditions (rainfall, evapotranspiration), on the initial state, and on the soil hydrodynamic properties.

The rain is considered as an inherently stochastic input. This can be formalized by considering a random event $\omega$ (Eq. (1)) to represent the state of nature (humidity, temperature) leading to a specific rain fall timeseries realization $Z(\omega)$. The probability space $\Omega$ is unknown, except indirectly, through previous measurements of rain fall timeseries $\{z_1, z_2, \dots, z_R\}$. Here, we choose to limit the probability space to a specific type of event, that is six-hour long summer rains measured locally between 1997 and 2015 (Lagouy et al., 2024) and following Cat-alogne et al. (2016). On the other hand, the uncertainty on the initial state of the soil moisture is not considered, and is fixed to the same initial conditions in all simulations. The chosen initial conditions are characterized by setting a hydrostatic equilibrium.

The probability distributions of the soil parameters represent the variability of the parameters across the catchment. The distributions are set to be representative of the available data on the soil, vegetation and topological properties of the Morcille catchment. They are assigned based on experimental measurements from the catchment of application, available scientific literature and expert knowledge. For more details, we refer to Rouzies et al. (2023). To account for the vertical variability of the soil characteristics, the vertical soil column is decomposed in three soil horizons, each considered internally homogeneous when it comes to soil hydrodynamical properties.

As only one soil column is studied in this work (Figure, 1), only six parameters whose sensibility indices are of interest are chosen out of the original 145 PESHMELBA parameters by performing a screening with the Hilbert-Schmidt Independence Criterion (HSIC) criterion using 3000 model simulations (Iooss et al., 2023). The input parameters representing soil properties used in this study, along with their distributions, are listed in Table 1. The catchment soil is mainly sandy.

### 3.2. Representation of soil moisture profiles in PESHMELBA

In this study case, each spatial element is represented by a four-meter deep vertical column composed of 25 cells of varying thickness. The cell thickness varies from 5 mm up to 1 m, with the finest discretization being at the soil surface. The 25 cells are considered to come from three soil horizons. The soil horizons are called: surface — surf. (depth 1 cm, composed of two cells), intermediate — inter. (depth 10 cm, seven cells) and deep (depth 4 m, 16 cells), see Fig. 2.

Fig. 2(a) illustrates the rains used in the study. The cumulative amount of the rainfall events amounts to 15–75 mm in 6 h, which presents a significant range of rain events. The same Fig. 2(b) highlights the differences in moisture profiles coming from both rains and soil
<caption>Table 1</caption>

Distributions of the Parameters Considered for the Global Sensitivity Analysis. The suffixes *surf*, *inter*, *deep* denote the soil horizon. The parameter $\theta_{s,deep}$ is readjusted through a truncation on the gaussian to the interval [0, 1].

<table><thead><tr><th>Name</th><th>Definition</th><th>Unit</th><th>Distribution</th></tr></thead><tbody><tr><td>$\theta_{s,surf}$</td><td>water content at saturation</td><td>$[L^3 L^{-3}]$</td><td>$N^*(0.3375, 0.0338^2)$</td></tr><tr><td>$\theta_{s,inter}$</td><td>water content at saturation</td><td>$[L^3 L^{-3}]$</td><td>$N^*(0.3322, 0.0332^2)$</td></tr><tr><td>$\theta_{s,deep}$</td><td>water content at saturation</td><td>$[L^3 L^{-3}]$</td><td>$N^*(0.316, 0.0316^2)$</td></tr><tr><td>$\theta_{r,deep}$</td><td>residual water content</td><td>$[L^3 L^{-3}]$</td><td>$N^*(0.0612, 0.0153^2)$</td></tr><tr><td>$mn_{deep}$</td><td>Van Genuchten retention curve parameter</td><td>$[-]$</td><td>$N^*(0.1791, 0.0179^2)$</td></tr><tr><td>$hg_{deep}$</td><td>Van Genuchten retention curve parameter</td><td>$[-]$</td><td>$N^*(-9.69, 0.969^2)$</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 1. The decomposition of the virtual catchment in 14 plots belonging to one of the three possible soil units, depicted by different fillings, and two vegetation types, (vegetative filter strips and vineyard plots). The separation of the soil column under each plot into three soil horizons (surface, intermediary and deep) is illustrated. The fourth parcel is put into evidence, as its soil moisture profile is observed. Image adapted from Rouzies et al. (2023).</figcaption></figure>

hydrodynamic properties. It can be seen that the heavier rain, $z_B$, leads to a moisture profile saturated up to the 19th cell for both sets of soil hydrodynamic properties. The rain $z_A$ leads to moisture profiles with a shallower saturation front (14th and 15th cell, depending on the set of hydrodynamic properties). On the other hand, the most visible difference for the example shown between two moisture profiles coming from the same rain, but different hydrodynamic properties, is the difference between the water contents of the deep horizon, more precisely the values of cells 8 to 19 for $M_s(\cdot, z_B)$, resp. cells 8 to 14 for $M_s(\cdot, z_A)$.

In this study, the model PESHMELBA is denoted with $M_s$, while the operator $\mathcal{H}$ is used to denote the mapping of the model state to the moisture profiles, as only the moisture profiles are our output of interest. The (scalar) function of interest whose sensitivity to input parameters is studied, is defined as the $L^2$ norm of the difference between the simulated moisture profile and an observed one $y_{obs}$, which is considered fixed in all simulations:

$$f_s: D \times \Omega \to \mathbb{R}, \\ (X, \omega) \mapsto |\mathcal{H}[M_s(X, Z(\omega))] - y_{obs}|^2, \qquad (16)$$

where $\omega$ is an abstract random event, whose realization leads to one particular rain measurement $z_i = Z(\omega_i)$ taking place and $X$ are the input parameters representing soil hydrodynamic properties whose trajectory-based Sobol' indices we are interested in.

## 4. Results

In this section, we first apply classical global sensitivity analysis to two separate rainfall events, to see how the most influential parameters differ in these two cases. Then, the Sobol' indices and the associated rankings are obtained for all the available measured rainfall events described in Section 3. As all Sobol' indices are obtained through PCE metamodels, the precision of the metamodels is also commented and evaluated through the $Q^2$ coefficient. Finally, to obtain a single ranking

over the $\Omega$ space, the aggregated sensitivity indices are calculated and compared with the ranking of the Sobol' indices means. Replication data and code are available on the public repositories Radišić (2024a,b).

### 4.1. Global sensitivity analysis in two forcings

To see whether the GSA results differ in two cases of forcings, the Sobol' indices are obtained for the six input parameters under two rain events, denoted $z_A$ and $z_B$ (introduced in Fig. 2). For the calculation of the Sobol' indices a space-filling experimental design of $N_{train} = 50$ points is generated in the six-dimensional parameter space as a Latin Hypercube Sample, that is a space-filling design which fills optimally the six-dimensional parameter space, (Carnell, 2024). A PCE metamodel is built from the 50 PESHMELBA simulations, with least-angle regression (Marelli and Sudret, 2014). Sobol' indices are obtained from the PCE following Eq. (5).

Fig. 3 illustrates the first order and total Sobol' indices of the six input parameters for the rain events $z_A$ and $z_B$. The GSA results under $z_B$ show that the most influential parameter is $\theta_{s,deep}$, accounting for more than 70% of the output variance, whereas the last ranked parameter is $h_{g,deep}$, presenting a virtually null Sobol' index, both for the first order and interactions. The parameters $\theta_{s,inter}$ and $\theta_{s,surf}$ account for around 10% of the output variance, and the closeness of the values of the first and total Sobol' indices shows the absence of interactions with other input parameters. Inversely, $mn_{deep}$ and $hg_{deep}$ show that most of their contribution to output variability is through indirect effects as interactions with other parameters.

The GSA results obtained under rain $z_B$ show that to reduce the output uncertainty conditioned to this rain, the efforts should primarily be put into a more precise knowledge of the parameter $\theta_{s,deep}$. However, when the rain event $z_A$ is used as the model forcing, the Sobol' index of the parameter $\theta_{s,deep}$ is significantly smaller, leading to this parameter not being ranked first anymore. Indeed, under the rain $z_A$, it is the parameter $\theta_{s,inter}$ that accounts for the largest part of the output variance. The parameter sensitivities differ significantly in the two situations, as the first ranked parameter is not the same under the two rain events. In the next section, the sensitivity analysis is explored for all available rain events.

### 4.2. Trajectory-based Sobol' indices

In this section, we explore the values of Sobol' indices under all available rain realizations. First, the quality of the metamodels used for the estimation of Sobol' indices is evaluated. Then, the distributions of the Sobol' indices under the available rain realizations are studied. Finally, two summarizing quantities, the mean of the Sobol' indices and the aggregated sensitivity index (ASI, Eq. (15)) are compared and commented.

To obtain Sobol' indices under all available rain measurements, a PCE metamodel is built for each rain fall timeseries realization $z_r$. Each metamodel is fitted on a training set of the original model evaluations $\{(x_j, f_s(x_j, z_r))\}_{j=1}^{N_{train}}$. The points $\{x_j\}_{j=1}^{N_{train}}$ are chosen as a Latin Hypercube Sample. To reduce the number of training points needed, least-angle regression is used (Marelli and Sudret, 2014). To ensure consistency between the original model and the metamodels,
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** (a) Cumulative rainfall events as a function of time, two rain realizations (z<sub>A</sub> and z<sub>B</sub>) are put in evidence by differing colors and line types. (b) Moisture profiles at the end of the six hours. Two sets of soil hydrodynamic properties are used in two different simulations for each of the rains z<sub>A</sub> and z<sub>B</sub>. The heavier rain, z<sub>B</sub> (orange dotted line) leads to moisture profiles with a deeper saturation front than the lighter rain z<sub>A</sub> (green two-dashed line). Horizontal dashed lines show the soil horizons. (c) Soil description of the vertical column and its three soil horizons.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** First and total order Sobol' indices with corresponding 95% bootstrap confidence intervals for the six input parameters, when the rain events are z<sub>B</sub> (heavier rain) and z<sub>A</sub> (lighter rain), introduced in Fig. 2.</figcaption></figure>

the samples size of the experimental design {$x_j$}$_{j=1}^{N_{\text{rain}}}$ is increased until reaching a coefficient $Q^2 > 0.95$ for each metamodel $f_{\text{PCE}}^{(r)}(x)$. The $Q^2$ (Eq. (9)) is evaluated on an independent test set of 100 PESHMELBA simulations for each rain. For this case study, this meant using a LHS of $N_{\text{train}} = 50$ points for the training set, which is very low compared to the classical Sobol' sampling, Fig. 4.

The total Sobol' indices obtained from the PCE metamodels are shown on Fig. 5 as histograms. The total Sobol' indices are presented for five input parameters, as the Sobol' indices of parameter $hg_{\text{deep}}$ are found to be virtually null in all rain realizations. It can be seen that the histogram of $\theta_{\text{s,deep}}$ covers a wide range of values, reaching from a total Sobol' index from 0.20 up to 0.80. The parameter $\theta_{\text{s,inter}}$ covers ranges from 0.05 to 0.45, overlapping with the histogram of $\theta_{\text{s,deep}}$ and having a larger Sobol' index for certain rain realizations. The histograms of $mn_{\text{deep}}$ and $\theta_{\text{r,deep}}$ cover smaller ranges, showing a more consistent influence on the output over different rain realizations.

For an easier interpretation of the results, the parameters are sorted in decreasing order from 1 to 6 for each rain realization. The number of appearances at each rank is counted. Fig. 5 lists the six input parameters along with their number of appearances at each rank. A color coding

is introduced to distinguish between two groups of ranking results: rankings listing $\theta_{\text{s,deep}}$ first and rankings that do not. When the color coding is applied to the cumulative rains leading to such rankings, it can be seen that the rains with high cumulative precipitation (more than 50 mm in 6 h) all lead to $\theta_{\text{s,deep}}$ being ranked first. This seems reasonable, as heavier rains lead to deeper soil saturation fronts, thus the major changes of moisture profile after such event come from the change of the water content at saturation.

However, among the rains with less than 50 mm accumulated in 6 h, the ranking of the first parameter is not consistent. While $\theta_{\text{s,deep}}$ remains the first ranked parameter for most rain realizations, there is a significant proportion of rain events where variations in $\theta_{\text{s,inter}}$ contribute the most to output variability. The lack of a clear distinction between these two groups of rain events shows that it is not straightforward to predict the results of a GSA solely based on the similarity of cumulative rain series. This highlights the importance of conducting a GSA under an ensemble of forcing scenarios to accurately assess the impact of different parameters.
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Left: the increase of the $Q^2$ coefficient with the number of training points used for evaluating the PCE metamodel, for all rain events. Each line corresponds to a metamodel constructed for one trajectory, i.e. one fixed rain event. When the number of training points reaches $N_{train} = 50$, the coefficient $Q^2$ is above 0.95 for all rain realizations. Right: comparison between outputs produced with the original model $f^{(r)}(x)$ and the outputs produced with the PCE metamodels $f_{PCE}^{(r)}(x)$, for all rain events, $r \in \{1..R\}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Top left: histograms of the total Sobol' indices, obtained for the input parameters under different rain realizations. The parameter $hg_{deep}$ is not listed, as its Sobol' index is virtually null in all rain realizations. Top right: the cumulative rain over six hours of the 192 rainfall events. Bottom: the number of appearances of the input parameters at each rank (from 1st to 6th), for a ranking based on the total Sobol' index. The two colors correspond to the rank of the parameter $\theta_{s,deep}$, distinguishing between the rankings that list $\theta_{s,deep}$ as 1st and the rankings that do not.</figcaption></figure>

Furthermore, when GSA is used to identify non-influential parameters, disregarding the full range of rain variability can result in overlooking rain events that highlight its contribution to output variability. For example, in the case of the histogram representing the Sobol' parameters of $\theta_{s,surf}$ (Fig. 5), certain realizations are very close to zero. This suggests that ignoring the rain variability could lead to not observing rain events that highlight its contribution, and could lead to a mistaken conclusion that $\theta_{s,surf}$ has no impact on output variability. In this case, the only parameter found to have no influence under any of the rain events examined is $hg_{deep}$.

### 4.3. Aggregation over the uncontrollable space

The results in Fig. 5 show how the rankings of the parameters can differ with respect to the rain realization used for the forcing of the model. The emphasis is put on the disagreement among the parameters appearing at the top of the list, as it can potentially lead to diverging choices in further data acquisition and model calibration.

Indeed, in a context where the aim of the GSA is to determine the order of parameters in which the reduction of their uncertainty contributes the most to the reduction of the output uncertainty, the results show that different conclusions can be made depending on the rain realization. Thus, while reducing the uncertainty of the parameter $\theta_{s,deep}$ contributes greatly to the reduction of output uncertainty under rain $z_B$, it would contribute less to the reduction of output uncertainty under $z_A$ (Fig. 3).

When using GSA to select the parameters whose uncertainty reduction is a top priority, a unique ranking of the parameters should consider the contributions of each parameter to the total output variance over the entire probability space $\Omega$ representing the uncontrollable uncertainties. To do this, the contributions of each parameter should be weighted by the total variance of the output conditioned on each rain event. In other words, if the output uncertainty conditioned on $z_A$ is larger than the one conditioned on $z_B$, prioritizing the reduction of uncertainties associated with parameters that have an impact on
<figure><img src="image_9.png" /><figcaption>Fig. 6. The sample mean (Eq. (12)) and the aggregated sensitivity indices ASI indices (Eq. (15)) of the total Sobol' indices for the six input parameters, with corresponding confidence intervals.</figcaption></figure>

the (already small) uncertainty under $z_B$ may not contribute to overall uncertainty reduction across the whole domain of rain realizations.

The Sobol' indices previously obtained are used to calculate the aggregated sensitivity indices (ASI, Eq. (15)). The values of ASI and the mean coincide in cases where the conditional variances of the output are the same in all rain realizations. Fig. 6 reports the ASI$_{T_i}$ indices, along with the sample mean of the Sobol' indices (Eq. (12)). The two quantities are very close, reflecting the steadiness of the output variance conditioned to different rain realizations.

The difference between the ASI and the Sobol' means in our case study is not very significant, as both the mean and the ASI give essentially the same rankings. The parameter $\theta_{s,deep}$ is ranked first, both when looking the total Sobol' index mean (in mean, $\theta_{s,deep}$ accounts for 50% of the output variance) and the corresponding ASI (50% of the total output variance). It is only for the parameter $\theta_{s,inter}$ that the ASI is slightly lower than the $\mu_{T_i}$, meaning that the output variances conditioned on rains leading to high ranks of $\theta_{s,inter}$ are lower than the output variances conditioned on rains ranking $\theta_{s,inter}$ lower. The parameter $hg_{deep}$ shows a virtually null ASI index, confirming that its influence over the whole domain of forcing variations is of no importance.

## 5. Discussion and conclusion

In this work, we perform a global sensitivity analysis of a model output by disentangling the uncertainty of an uncontrollable external forcing (rain), and the uncertainty of the controllable operational parameters (soil hydrodynamic properties). To do so, we consider the Sobol' indices of the parameters as random variables. The randomness of the Sobol' indices comes from the stochastic nature of the rain. We show that the ranking of the parameters differs with respect to the rain realization used. Indeed, while the most influential parameter in most cases is the water content at saturation of the deep horizon, there is also a significant amount of rain realizations under which the lead is taken by the water content at saturation of the intermediate horizon. Finally, we aggregate the Sobol' indices by taking into consideration the contributions of each parameter to overall uncertainty reduction across the whole domain of rain realizations.

### 5.1. On the source of Sobol' index variability

Variable Sobol' indices of operational parameters can be expected in models where interactions are present between the model parameters and the uncontrollable external forcing. Indeed, consider a stochastic simulator as the one used in this work (Eq. (10)), and consider its

Hoeffding-Sobol' (ANOVA) decomposition. The summands of this decomposition can then be grouped according to whether they involve operational parameters, stochastic inputs or both:

$$f_s(\mathbf{X}, \mathbf{Z}) = f_0 + f_{s1}(\mathbf{X}) + f_{s2}(\mathbf{Z}) + f_{s12}(\mathbf{X}, \mathbf{Z}) \quad (17)$$

Consider the interaction term between the stochastic and operational variables to be negligible. In this case, the target function is essentially additive with respect to the two groups of variables:

$$f_s(\mathbf{X}, \mathbf{Z}) \approx f_0 + f_{s1}(\mathbf{X}) + f_{s2}(\mathbf{Z}) \quad (18)$$

Note that, for models where the interaction term is negligible, the Sobol' indices of the operational variables $X$ do not depend on the realization of the stochastic input. Indeed, the variations of the operational parameters influences the output only through the term $f_{s1}$. Hence, the trajectory-based Sobol' indices (Eq. (11)) for a realization $z_r$ can be rewritten:

$$\begin{aligned} S_i^{(r)} &= \frac{\mathrm{Var}[\mathbb{E}[f_s(\mathbf{X}, z_r)|\mathbf{X}_i]]}{\mathrm{Var}[f_s(\mathbf{X}, z_r)]} \\ &= \frac{\mathrm{Var}[\mathbb{E}[f_{s1}(\mathbf{X})|\mathbf{X}_i] + f_0 + f_{s2}(z_r)]}{\mathrm{Var}[f_0 + f_{s1}(\mathbf{X}) + f_{s2}(z_r)]} \\ &= \frac{\mathrm{Var}[\mathbb{E}[f_{s1}(\mathbf{X})|\mathbf{X}_i]]}{\mathrm{Var}[f_{s1}(\mathbf{X})]}, \end{aligned}$$

showing that the Sobol' indices are independent of the realization $z_r$. Thus, it is only in the presence of a non-negligible interaction term $f_{s12}(\mathbf{X}, \mathbf{Z})$ the Sobol' indices of the operational parameters can vary with respect to the realization of the uncontrollable forcing.

This shows that, in presence of interactions between the parameters and the external forcing, the GSA-informed decisions about uncertainty reduction of the controllable parameters that would benefit the model under one realization of the external forcing, do not benefit the model over the whole uncontrollable domain $\Omega$.

### 5.2. on the implications of social policies

To support policy-making and the formulation of emergency measures, global sensitivity analysis can be used to identify the most impactful operational parameters on outputs of interest. This can help policy-makers focus on the crucial quantities to be defined. For example, in the serious game CAUSERIE (Gouy Boussada et al., 2024), management scenarios based on real catchments are discussed among watershed management, river syndicates and farmers. Some of the possible management decisions are choosing the widths of the vegetative filter zones (for example grass strips). GSA can help identify the grass strips whose width is the most important to define, thus focusing the managerial decisions on the most important operational parameters.

However, due to the limitations of traditional GSA when facing uncertain external forcing, the most influential operational parameter can vary with the external forcing. To integrate the uncertainty of external forcings in the decision-making process, Sobol' index distributions should be observed. However, given the possible complexity of such information, a more synthetic view is proposed with Aggregated Sensitivity Indices (ASI). The ranking based on ASI weights the contributions of each decision with its impact on the total output uncertainty, integrating it over the whole domain of forcing uncertainty.

### 5.3. On the aggregation

To enhance interdisciplinary cooperation and decision-making, it is crucial to transform complex technical model outputs into easily understandable management information. While the distributions of Sobol' indices provide detailed insights, they can be challenging to interpret. We propose using ASI to simplify this information through a weighted average, making it more accessible for decision-makers. ASI reduces the complexity of Sobol' indices distributions to a single index per parameter, facilitating translation into management information.
In this work, the aggregation considers all $ω$ realizations equally likely. Still, if one was interested in having different probabilities for the events, this could be accounted for by introducing corresponding weights in Eq. (15). Weights can be also be introduced if certain realizations of $Ω$ are more important than others (for example, realizations exceeding a certain threshold and for which the model should be particularly precise). For example, in flood management, instead of giving equal weights to all realizations $ω$, rainfall events leading to floods could be weighted more, leading to the parameters influencing the variability of the output under these conditions being considered a more important source of uncertainty.

## 5.4. On the computational cost

From a computational aspect, estimating the distributions of the sensitivity indices adds an outer loop in the sensitivity analysis, thus needing even more model simulations than ordinary GSA. In this work, we have resorted to PCE metamodels as means of limiting the number of simulations needed in the inner loop. Optimizing the experimental design in the $Ω$ space is not discussed, as the number of rain measurements in our case study is rather limited. If a larger amount of realizations of the stochastic input is available, either from measurements or from stochastic simulators, (Leblois and Creutin, 2013), different ways of choosing a space-filling experimental design in $Ω$ could be compared. However, the optimal choice of a space-filling design greatly depends on the nature of $Ω$ and the way in which its randomness is propagated through the model and to the target function. If the distribution of the Sobol' indices is not the primary interest, cheaper methods can be used. Furthermore, if one only needs screening results, variance-based indices can be replaced with HSIC-based sensitivity indices. Fellmann et al. (2024) propose an HSIC adaptation for set-based outputs. This approach could be used to separate the uncertainty coming from operational parameters and external forcings by replacing the maps with the domain of forcing inputs.

## 5.5. On the sources of stochastic uncertainty

In our study, the stochastic uncertainty comes from the natural variability of rain realizations, and is thus specific to the rainfall characteristics of the studied region. It is important to have in mind that this variability, represented with the probability space $Ω$, is specific to the rainfall characteristics of the case study scenario and would be different for urban catchments, agricultural areas, or mountainous regions.

The proposed approach is generalizable to any other stochastic uncertainty, as it does not call for a specific structure of the abstract event $ω$. In other words, the stochastic input may not be parametrized, but a sufficient number of realizations of the event should be available to correctly represent the studied domain.

Other sources of stochastic uncertainty in hydrological models are the initial saturation in the soil columns or the choice of the vertical discretization of the soil column, but also temperature or evapotranspiration. Another common uncontrollable forcing in hydrological modeling is that of the future climate projections, which comes with a new form of uncertainty that may drastically change the model sensitivity results.

In our study, the size of the parameter space is relatively small. However, it already illustrates the possible differences of parameter importance under varying uncontrollable sources of uncertainty. The benefit coming from the distinction of the two sources of uncertainty is yet to be seen on larger, more complex case studies, where we could expect to find more differences between the sensitivity analyses. In the case of the PESHMELBA model, a particularly interesting source of uncertainty to be studied in future works is the uncertainty of the pesticide application dates, which have a significant impact on the model simulations. We believe that this approach can contribute to making more informed choices in the further steps of model uncertainty reduction, in particular for robust calibration and model design.

# CRediT authorship contribution statement

Katarina Radišić: Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Claire Lauverneta: Writing – review & editing, Supervision, Resources, Methodology, Funding acquisition, Conceptualization. Arthur Vidard: Writing – review & editing, Supervision, Resources, Methodology, Funding acquisition, Conceptualization.

# Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

# Acknowledgments

This research was supported by the INRAE-Inria thesis funding program and the LEFE-MANU MARQUISE project http://programmes.insu.cnrs.fr/lefe.

# Software and data availability

The data supporting all the figures and the code for recreating them are available on Zenodo doi.org/10.57745/C455R6, (Radišić, 2024a). The scripts used to obtain Sobol' indices use the open-source MATLAB software for uncertainty quantification UQLab Version 2.0.0 https://www.uqlab.com, (Marelli and Sudret, 2014). The simulations used for this work are available on doi.org/10.57745/4MDYNU, (Radišić, 2024b). They are obtained with the PESHMELBA version archived on Zenodo doi.org/10.15454/2HAU8R, (Rouzies and Lauverneta, 2022).

# References

Alipour, A., Jafarzadegan, K., Moradkhani, H., 2022. Global sensitivity analysis in hydrodynamic modeling and flood inundation mapping. Environ. Model. Softw. 152, 105398. http://dx.doi.org/10.1016/j.envsoft.2022.105398.

Azzi, S., 2020. Surrogate Modeling of Stochastic Simulators (Ph.D. thesis). Institut Polytechnique de Paris, France.

Buis, S., Piacentini, A., Déclat, D., the PALM Group, 2006. PALM: a computational framework for assembling high-performance computing applications. Concurr. Comput.: Pr. Exp. 18 (2), 231–245. http://dx.doi.org/10.1002/cpe.914.

Carmona-Cabrero, A., Munoz-Carpena, R., Oh, W.S., Muneepeerakul, R., 2024. Decomposing variance decomposition for stochastic models: Application to a proof-of-concept human migration agent-based model. J. Artif. Soc. Soc. Simul. 27 (1), 16. http://dx.doi.org/10.18564/jass.5174.

Carnell, R., 2024. lhs: Latin hypercube samples. R package version 1.1.7.

Catalogne, C., Lauverneta, C., Leblois, E., 2016. Analyse Des Structures Temporelles De Pluies Pour La Définition De Hyetogrammes En Entrée De La Chaine De Dimensionnement Des Bandes Tampons Végétalisées BUVARD. Technical Report, INRAE, Pollutions Diffuses, Villeurbanne, France.

Colombi, A., Bancheri, M., Acutis, M., Basile, A., Botta, M., Perego, A., 2024. A sound understanding of a cropping system model with the global sensitivity analysis. Environ. Model. Softw. 173, 105932. http://dx.doi.org/10.1016/j.envsoft.2023.105932.

Da Veiga, S., Gamboa, F., Iooss, B., Prieur, C., 2021. Basics and Trends in Sensitivity Analysis. Society for Industrial and Applied Mathematics, Philadelphia, PA, http://dx.doi.org/10.1137/1.9781611976694.

D'Andrea, M.F., Letourneau, G., Rousseau, A.N., Brodeur, J.C., 2020. Sensitivity analysis of the Pesticide in Water Calculator model for applications in the Pampa region of Argentina. Sci. Total Environ. 698, 134232. http://dx.doi.org/10.1016/j.scitotenv.2019.134232.

De Lozzo, M., Marrel, A., 2016. New improvements in the use of dependence measures for sensitivity analysis and screening. J. Stat. Comput. Simul. 86 (15), 3038–3058. http://dx.doi.org/10.1080/00949655.2016.1149854.

Dell'Oca, A., 2023. Sensitivity analysis: An operational picture. Water Resour. Res. 59, http://dx.doi.org/10.1029/2022WR033780, e2022WR033780.

Dell'Oca, A., Riva, M., Guadagnini, A., 2017. Moment-based metrics for global sensitivity analysis of hydrological systems. Hydrol. Earth Syst. Sci. 21 (12), 6219–6234. http://dx.doi.org/10.5194/hess-21-6219-2017.
Dubreuil, S., Berveiller, M., Petitjean, F., Salaün, M., 2014. Construction of bootstrap confidence intervals on sensitivity indices computed by polynomial chaos expansion. Reliab. Eng. Syst. Saf. 121, 263-275. http://dx.doi.org/10.1016/j.ress.2013.09.011.

Ernst, O.G., Mugler, A., Starkloff, H.-J., Ullmann, E., 2012. On the convergence of generalized polynomial chaos expansions. ESAIM Math. Model. Numer. Anal. 46 (2), 317-339. http://dx.doi.org/10.1051/m2an/2011045.

Faúndez Urbina, C.A., Van Den Berg, F., Van Dam, J.C., Tang, D.W.S., Ritsema, C.J., 2020. Parameter sensitivity of SWAP-PEARL models for pesticide leaching in macroporous soils. Vadose Zone J. 19 (1), e20075. http://dx.doi.org/10.1002/vzj2.20075.

Fellmann, N., Blanchet-Scalliet, C., Helbert, C., Spagnol, A., Sinoquet, D., 2024. Kernel-based sensitivity analysis for (excursion) sets. Technometrics http://dx.doi.org/10.1080/00401706.2024.2336537.

Fox, G.A., Muñoz-Carpena, R., Sabbagh, G.J., 2010. Influence of flow concentration on parameter importance and prediction uncertainty of pesticide trapping by vegetative filter strips. J. Hydrol. 384 (1), 164-173. http://dx.doi.org/10.1016/j.jhydrol.2010.01.020.

Gamboa, F., Janon, A., Klein, T., Lagnoux, A., 2014. Sensitivity analysis for multidimensional and functional outputs. Electron. J. Stat. 8 (1), http://dx.doi.org/10.1214/14-EJS895.

García, D., Arostegui, I., Prellezo, R., 2019. Robust combination of the Morris and Sobol methods in complex multidimensional models. Environ. Model. Softw. 122, 104517. http://dx.doi.org/10.1016/j.envsoft.2019.104517.

Gatel, L., Lauvernnet, C., Carluer, N., Weill, S., Paniconi, C., 2019. Sobol global sensitivity analysis of a coupled surface/subsurface water flow and reactive solute transfer model on a real hillslope. Water 12 (1), 121. http://dx.doi.org/10.3390/w12010121.

Gouy, V., Liger, L., Ahrouch, S., Bonnineau, C., Carluer, N., Chaumot, A., Coquery, M., Dabrin, A., Margoum, C., Pesce, S., 2021. Ardières-Morcillein the Beaujolais, France: A research catchment dedicated to study of the transport and impacts of diffuse agricultural pollution in rivers. Hydrol. Process. 35 (10), http://dx.doi.org/10.1002/hyp.14384.

Gouy Boussada, V., Seguin, L., Malingrey, S., Armani, G., Barreteau, O., Leteurtre, E., Adoir, E., Luzi, H., Carluer, N., Grillot, J., Rabotin, M., Abrami, G., Lauvernnet, C., 2024. Can a serious hybrid game centred on the landscape encourage the emergence of levers for collective action to reconcile agriculture and water quality? Agronomic Innov. 96, 15-24. http://dx.doi.org/10.17180/CIAG-2024-VOL96-ART02-GB.

Hart, J.L., Alexanderian, A., Gremaud, P.A., 2017. Efficient computation of Sobol' indices for stochastic models. SIAM J. Sci. Comput. 39 (4), A1514-A1530. http://dx.doi.org/10.1137/16M106193X.

Heredia, M.B., Prieur, C., Eckert, N., 2021. Nonparametric estimation of aggregated Sobol' indices: Application to a depth averaged snow avalanche model. Reliab. Eng. Syst. Saf. 212, 107422. http://dx.doi.org/10.1016/j.ress.2020.107422.

Homma, T., Saltelli, A., 1996. Importance measures in global sensitivity analysis of nonlinear models. Reliab. Eng. Syst. Saf. 52 (1), 1-17. http://dx.doi.org/10.1016/0951-8320(96)00002-6.

Hong, T., Purucker, S.T., 2018. Spatiotemporal sensitivity analysis of vertical transport of pesticides in soil. Environ. Model. Softw. 105, 24-38. http://dx.doi.org/10.1016/j.envsoft.2018.03.018.

Iooss, B., Da Veiga, S., Janon, A., Pujol, G., Broto, w.c.f.B., Boumhaout, K., Delage, T., Amri, R.E., Fruth, J., Gilquin, L., Guillaume, J., Herin, M., Idrissi, M.I., Gratiet, L.L., Lemaître, P., Marrel, A., Meynaoui, A., Nelson, B.L., Monari, F., Oomen, R., Rakovec, O., Ramos, B., Roustant, O., Sarazin, G., Song, E., Staum, J., Sueur, R., Touati, T., Verges, V., Weber, F., 2023. sensitivity R package: Global sensitivity analysis of model outputs. http://dx.doi.org/10.32614/CRAN/package.sensitivity.

Iooss, B., Lemaître, P., 2015. A review on global sensitivity analysis methods. In: Dellino, G., Meloni, C. (Eds.), Uncertainty Management in Simulation-Optimization of Complex Systems. Operations Research/Computer Science Interfaces Series. Springer, Boston, MA, pp. 101-122. http://dx.doi.org/10.1007/978-1-4899-7547-8_5.

Iooss, B., Ribatet, M., 2009. Global sensitivity analysis of computer models with functional inputs. Reliab. Eng. Syst. Saf. 94 (7), 1194-1204. http://dx.doi.org/10.1016/j.ress.2008.09.010.

Jamous, M., Marsooli, R., Ayyad, M., 2023. Global sensitivity and uncertainty analysis of a coastal morphodynamic model using Polynomial Chaos Expansions. Environ. Model. Softw. 160, 105611. http://dx.doi.org/10.1016/j.envsoft.2022.105611.

Jimenez, M.N., Le Maître, O.P., Knio, O.M., 2017. Nonintrusive polynomial chaos expansions for sensitivity analysis in stochastic differential equations. SIAM/ ASA J. Uncertain. Quantif. 5 (1), 378-402. http://dx.doi.org/10.1137/16M1061989.

Lagouy, M., Branger, F., Thollet, F., Breil, P., Dramais, G., 2024. Suivi hydrologique du bassin versant périurbain de l'Yzeron. http://dx.doi.org/10.17180/OBS.YZERON.

Lauvernnet, C., Muñoz-Carpena, R., 2018. Shallow water table effects on water, sediment, and pesticide transport in vegetative filter strips Part 2: model coupling, application, factor importance, and uncertainty. Hydrol. Earth Syst. Sci. 22 (1), 71-87. http://dx.doi.org/10.5194/hess-22-71-2018.

Leblois, E., Creutin, J.-D., 2013. Space-time simulation of intermittent rainfall with prescribed advection field: Adaptation of the turning band method: Simulation of Rainfall with Advection Field. Water Resour. Res. 49 (6), 3375-3387. http://dx.doi.org/10.1002/wrcr.20190.

Li, L., Lu, Z., Wu, D., 2016. A new kind of sensitivity index for multivariate output. Reliab. Eng. Syst. Saf. 147, 123-131. http://dx.doi.org/10.1016/j.ress.2015.11.006.

Lüthen, N., Marelli, S., Sudret, B., 2021. Sparse polynomial chaos expansions: Literature survey and benchmark. SIAM/ ASA J. Uncertain. Quantif. 9 (2), 593-649. http://dx.doi.org/10.1137/20M1315774.

Lüthen, N., Marelli, S., Sudret, B., 2023. A spectral surrogate model for stochastic simulators computed from trajectory samples. Comput. Methods Appl. Mech. Engrg. 406, 115875. http://dx.doi.org/10.1016/j.cma.2022.115875.

Mai, J., 2023. Ten strategies towards successful calibration of environmental models. J. Hydrol. 620, 129414. http://dx.doi.org/10.1016/j.jhydrol.2023.129414.

Marelli, S., Sudret, B., 2014. UQLab: A framework for uncertainty quantification in matlab. In: Vulnerability, Uncertainty, and Risk. American Society of Civil Engineers, University of Liverpool, United Kingdom, pp. 2554-2563. http://dx.doi.org/10.1061/9780784413609.257.

Morris, M.D., 1991. Factorial sampling plans for preliminary computational experiments. Technometrics 33 (2), 161-174. http://dx.doi.org/10.1080/00401706.1991.10484804.

Nagel, J.B., Rieckermann, J., Sudret, B., 2020. Principal component analysis and sparse polynomial chaos expansions for global sensitivity analysis and model calibration: Application to urban drainage simulation. Reliab. Eng. Syst. Saf. 195, 106737. http://dx.doi.org/10.1016/j.ress.2019.106737.

Nossent, J., Elsen, P., Bauwens, W., 2011. Sobol' sensitivity analysis of a complex environmental model. Environ. Model. Softw. 26 (12), 1515-1525. http://dx.doi.org/10.1016/j.envsoft.2011.08.010.

Perrin, T., Roustant, O., Rohmer, J., Alata, O., Naulin, J., Idier, D., Pedreros, R., Moncoulon, D., Tinard, P., 2021. Functional principal component analysis for global sensitivity analysis of model with spatial output. Reliab. Eng. Syst. Saf. 211, 107522. http://dx.doi.org/10.1016/j.ress.2021.107522.

Radišić, K., 2024a. Replication Data and Code for: "Impact of input forcing variability on the global sensitivity analysis of a hydrological model". http://dx.doi.org/10.57745/C455R6, Recherche Data Gouv.

Radišić, K., 2024b. Replication Data for: Impact of input forcings variability on the global sensitivity analysis of a hydrological model. http://dx.doi.org/10.57745/4MDYNU, Recherche Data Gouv. [Dataset].

Radišić, K., Rouzies, E., Lauvernnet, C., Vidard, A., 2023. Global sensitivity analysis of the dynamics of a distributed hydrological model at the catchment scale. Socio-Environ. Syst. Model. 5, 18570. http://dx.doi.org/10.18174/sesmo.18570.

Ross, P.J., 2003. Modeling soil water and solute transport—Fast, simplified numerical solutions. Agron. J. 95 (6), 1352-1361. http://dx.doi.org/10.2134/agronj2003.1352.

Rouzies, E., Lauvernnet, C., 2022. Peshmelba v2.0 code for global sensitivity analysis. http://dx.doi.org/10.15454/2HAU8R, Recherche Data Gouv.

Rouzies, E., Lauvernnet, C., Barachet, C., Morel, T., Branger, F., Braud, I., Carluer, N., 2019. From agricultural catchment to management scenarios: A modular tool to assess effects of landscape features on water and pesticide behavior. Sci. Total Environ. 671, 1144-1160. http://dx.doi.org/10.1016/j.scitotenv.2019.03.060.

Rouzies, E., Lauvernnet, C., Sudret, B., Vidard, A., 2023. How is a global sensitivity analysis of a catchment-scale, distributed pesticide transfer model performed? Application to the PESHMELBA model. Geosci. Model. Dev. 16 (11), 3137-3163. http://dx.doi.org/10.5194/gmd-16-3137-2023.

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S., 2008. Global Sensitivity Analysis: the Primer. John Wiley & Sons, http://dx.doi.org/10.1002/9780470725184.

Saltelli, A., Tarantola, S., Campolongo, F., Ratto, M., 2004. Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models. John Wiley & Sons, ISBN: 9780470870952, http://dx.doi.org/10.1002/0470870958.

Saltelli, A., Tarantola, S., Chan, K.-S., 1999. A quantitative model-independent method for global sensitivity analysis of model output. Technometrics 41 (1), 39-56. http://dx.doi.org/10.1080/00401706.1999.10485594.

Sobol, I.M., 1993. Sensitivity estimates for nonlinear mathematical models. Math. Model. Comput. Exp. 1, 407-414.

Sobol, I., 2001. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math. Comput. Simulation 55 (1-3), 271-280. http://dx.doi.org/10.1016/S0378-4754(00)00270-6.

Song, X., Zhang, J., Zhan, C., Xuan, Y., Ye, M., Xu, C., 2015. Global sensitivity analysis in hydrological modeling: Review of concepts, methods, theoretical framework, and applications. J. Hydrol. 523, 739-757. http://dx.doi.org/10.1016/j.jhydrol.2015.02.013.

Sudret, B., 2008. Global sensitivity analysis using polynomial chaos expansions. Reliab. Eng. Syst. Saf. 93 (7), 964-979. http://dx.doi.org/10.1016/j.ress.2007.04.002.
Tang, Y., Reed, P., Wagener, T., Van Werkhoven, K., 2007. Comparing sensitivity analysis methods to advance lumped watershed model identification and evaluation. Hydrol. Earth Syst. Sci. 11 (2), 793–817. http://dx.doi.org/10.5194/hess-11-793-2007.

Tarantola, S., Ferretti, F., Lo Piano, S., Kozlova, M., Lachi, A., Rosati, R., Puy, A., Roy, P., Vannucci, G., Kuc-Czarnecka, M., Saltelli, A., 2024. An annotated timeline of sensitivity analysis. Environ. Model. Softw. 174, 105977. http://dx.doi.org/10.1016/j.envsoft.2024.105977.

Tissot, J.-Y., Prieur, C., 2015. A randomized orthogonal array-based procedure for the estimation of first- and second-order Sobol' indices. J. Stat. Comput. Simul. 85 (7), 1358–1381. http://dx.doi.org/10.1080/00949655.2014.971799.

Wagener, T., Pianosi, F., 2019. What has Global Sensitivity Analysis ever done for us? A systematic review to support scientific advancement and to inform policy-making in earth system modelling. Earth-Sci. Rev. 194, 1–18. http://dx.doi.org/10.1016/j.earscirev.2019.04.006.

Xiu, D., Karniadakis, G.E., 2002. The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM J. Sci. Comput. 24 (2), 619–644. http://dx.doi.org/10.1137/S1064827501387826.

Zhu, X., Sudret, B., 2021. Global sensitivity analysis for stochastic simulators based on generalized lambda surrogate models. Reliab. Eng. Syst. Saf. 214, 107815. http://dx.doi.org/10.1016/j.ress.2021.107815.