# Global sensitivity analysis of the dynamics of a distributed hydrological model at the catchment scale

Katarina Radišić, Emilie Rouzies, Claire Lauvernet, Arthur Vidard

► To cite this version:

Katarina Radišić, Emilie Rouzies, Claire Lauvernet, Arthur Vidard. Global sensitivity analysis of the dynamics of a distributed hydrological model at the catchment scale. Socio-Environmental Systems Modelling, 2023, 5, pp.18570:1-17. 10.18174/sesmo.18570. hal-03884443v2

HAL Id: hal-03884443

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

Submitted on 26 Oct 2023

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
# Global sensitivity analysis of the dynamics of a distributed hydrological model at the catchment scale

Radišić Katarina<fnref n="1" />,<fnref n="2" />, Rouzies Emilie<fnref n="1" />, Lauvernet Claire<fnref n="1" />, and Vidard Arthur<fnref n="2" />

## Abstract

The PESHMELBA model simulates water and pesticide transfers at the catchment scale. Its objective is to help the process of decision-making in the common management of long-term water quality. Performing the global sensitivity analysis (GSA) of this type of model is necessary to trace the output variability to the input parameters. The goal of the present work is to perform a GSA, while considering the spatio-temporal nature and the high dimensionality of the model. The output considered is the surface moisture simulated over a two-month period on a catchment of assorted mesh elements (plots). The GSA is performed on the dynamical outputs, rewritten through their functional principal components. Sobol' indices are then estimated through polynomial chaos expansion on each principal component. The analysis differs between the two types of behavior observed in the surface moisture outputs. The hydrodynamic properties of the surface soil have a dominant influence on the average surface moisture. Nonetheless, the parameters describing deeper soil layers influence the output dynamics of those plots where the surface moisture is saturated. We obtain Sobol' indices with high precision while using a limited number of model estimations and considering the model's spatio-temporal nature. The physical interpretation of the GSA confirms and augments our knowledge on the model.

## Keywords

Global sensitivity analysis; hydrology; functional principal components; polynomial chaos expansion; distributed model

## 1 Introduction

Pesticide use is a major issue in sustainable agriculture and water quality. Therefore, it is important to have the knowledge and the tools to best estimate the risks associated with their use and propose appropriate corrective actions. The PESHMELBA model (Pesticides and Hydrology: modeling at the catchment scale), Rouzies et al. (2019) simulates processes involved in water and pesticide transfers at the catchment scale, in order to compare different scenarios of landscape management and their impact on water quality.

An important step in the journey to PESHMELBA's operational use is to understand the model variability. Global sensitivity analysis (GSA) is useful in this case, as it can trace the output variability back to its input parameters. It can also verify the model consistency with respect to the physical processes and enhance the understanding of the modeled behavior. However, PESHMELBA outputs are both space and time dependent and the application of global sensitivity analysis to spatio-temporal environmental models can be challenging and dependent on the particular case studied.

Traditionally, variance-based sensitivity analysis Sobol' (2001) was developed for scalar outputs (Saltelli et al., 2008). For models with multivariate output (spatial or temporal), extensions of the Sobol' indices were introduced in Lamboni et al. (2011); Gamboa et al. (2014); Xu et al. (2019) or Roux et al. (2021). Other extensions explore dependence measures that are not variance-based. Such measures include the HSIC (Da Veiga, 2015; Marrel et al., 2015) or random forests (Antoniadis et al., 2021).

<footnote n="1">¹INRAE, RiverLy, 69625, Villeurbanne, France</footnote>

<footnote n="2">²Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France</footnote>
40

45

50

However, by obtaining the sensitivity indices for the whole multivariate output, one can overlook the sensitivities of locally-based behaviors. Thereby, Saint-Geours et al. (2014) distinct between GSA at multiple scales. The smallest scale would be to perform GSA on each component of the output, while larger scale GSA studies indices with respect to several components included in a domain of interest. As a matter of fact, if Sobol' indices are available for each component, it is possible to aggregate them through a weighted sum to obtain the generalized Sobol' indices from Gamboa et al. (2014). In such a manner, the non-aggregated indices bring insight into local behaviors' sensitivities, while the aggregated indices synthesize this information over the whole domain of interest (spatial or temporal). Still, for models whose outputs have high temporal (or spatial) resolution, performing local sensitivity analysis on each component (time step or mesh cell) is not always possible. Furthermore, Campbell et al. (2006) argue that this can lead to redundant results, and propose to project high dimensional outputs to polynomial or adaptive bases, then perform GSA of the components in the reduced-dimension space. This approach is particularly interesting if the extracted components can also be interpreted as physical signatures of the model behavior. The case of functional principal component basis and the corresponding aggregation of the Sobol' indices obtained on each component is studied in (Lamboni et al., 2011; Xiao & Li, 2016).

55

In environmental science modeling, the extensions were used in various areas. For instance, a sensitivity analysis on the spatialized maximum water depth of a coastal flooding risk model in Perrin et al. (2021); both a temporal and a spatial sensitivity analysis of a model describing nitrogen transfers in Ferrer Savall et al. (2019); and a comparison between various sensitivity indices on a spatio-temporal radionuclide atmospheric dispersion model in De Lozzo & Marrel (2017). A sensitivity analysis was also performed for spatially distributed PESH-MELBA outputs in Rouzies et al. (2023) by calculating site specific sensitivity indices, then aggregating them to blocks via the definition provided in Gamboa et al. (2014), but without studying the sensitivity of the dynamics of PESHMELBA.

60

As the previous examples show, there is no universal methodology for performing a sensitivity analysis on spatio-temporal models. Both the choice of the sensitivity index and its estimator can be guided by various model properties such as the space resolution, time resolution or the nature of the input parameters (scalar or multivariate). Other important aspects are the number of input parameters and the cost of a model run, since the number of model runs needed for sensitivity analysis increases with the number of input parameters.

65

70

Metamodels, in particular polynomial chaos expansion (PCE), can be used for sensitivity index calculation at a reasonable cost and high efficiency. They are especially useful when the model studied is computationally expensive or has a large input space. Indeed, building a PCE via Least Angle Regression (Blatman & Sudret, 2011) needs a significantly lower number of model simulations than other Monte Carlo based methods. Furthermore, once the metamodel is built, the sensitivity index calculation is analytical (Sudret, 2008). Applications of PCE for Sobol index calculation can be found in (Fajraoui et al., 2011) or in (Rouzies et al., 2023). For multidimensional outputs, principal component analysis (PCA) can be used prior to PCE to reduce the dimensionality, as done for an urban drainage model in Nagel et al. (2020). Furthermore, PCA can be used to find particular geometric features in the output. In that case, sensitivity indices can bring insight in the underlying physical processes of the model, as proposed in Garcia-Cabrejo & Valocchi (2014) and tested on a simple problem of reactive transport in porous media. To the best of our knowledge, such methods were not yet deployed on hydrological or pesticide transfer models.

75

In the present work, we perform a global sensitivity analysis of the spatio-temporal, highly dimensional, process-based water and pesticide transfer model PESHMELBA. We are particularly interested in detecting different behaviors in the model dynamics, and analyzing which inputs drive their occurrence. Another objective is to deploy this type of method for spatio-temporal outputs, such as those of distributed hydrological models (whether physically-based or semi-conceptual), as this has never been done on spatio-temporal variables in hydrology, but essentially on integrated variables, mostly cumulative outputs or quality criteria such as Nash-Sutcliffe model efficiency coefficient (Saltelli et al., 2019; Rouzies et al., 2023).

80

85

The outputs considered are the surface moisture outputs simulated over a two-month period on a catchment of assorted mesh elements (plots). In this study, the plots are observed one at a time. The sensitivity analysis is performed on their temporal outputs, which are rewritten through their principal components. Then, Sobol' indices are estimated through polynomial chaos expansion on each principal component.
The article is organized as follows. Section 2 describes the methodology, a brief description of Sobol' indices and their estimation via polynomial chaos expansion, then the functional principal component decomposition and two generalizations of sensitivity indices when the outputs are multidimensional. Section 3 presents the case study. The model PESHMELBA is introduced and the meaning of its input parameters explained. The results are presented in Section 4 in the following order. First, the resulting functional principal components are interpreted for each plot, then each PCE is validated, the Sobol' indices obtained are reported and commented by making a connection with the physics of the model. Finally, the Sobol' indices are aggregated temporally and spatially to the catchment scale.

## 2 Methodology

Let $\mathcal{M}: \mathbb{R}^K \to \mathbb{R}$ be a model with $K$ mutually independent scalar inputs $X_1, ..., X_K$ and a scalar output $Y$.

$$Y = \mathcal{M}(X_1, ..., X_K) = \mathcal{M}(\mathbf{X})$$

### 2.1 Sobol' Indices

The Sobol' indices are a widely accepted way of quantifying the sensitivity of the model output on the input parameters. The Sobol' approach is variational; it studies the impacts of input parameters on the output variance.

Let $D$ be the variance of the model output $Y$:

$$D = \mathrm{Var}(Y).$$

Let $D_i$ be the variance of the expectation of the output conditioned on input parameter $X_i$:

$$D_i = \mathrm{Var}(\mathbb{E}[Y|X_i]).$$

Then, the first order Sobol' indices of the input parameter $X_i$ are defined as the ratio between the conditional variance and the total variance of the output (Sobol', 2001):

$$S_i = D_i/D.$$

Let $D_{i_1, ..., i_s}$ be the variance of the expectation of the output conditioned on the set of input parameters $\{X_{i_1}, ..., X_{i_s}\}$ minus the variance of the expectation of the output conditioned on all other subsets of the given set:

$$D_{i_1, ..., i_s} = \mathrm{Var}(\mathbb{E}[Y|X_{i_1}, ..., X_{i_s}]) - \sum_{\mathcal{I}_s} D_{i_1, ..., i_s}, \\ \mathcal{I}_s = \{\{i_1, ..., i_k\} : \{i_1, ..., i_k\} \subseteq \{i_1, ..., i_s\}\}.$$

Similarly, the Sobol index of the set of input parameters $X_{i_1}, ..., X_{i_s}$ is defined as the ratio between the conditional variance of the output to those parameters minus the variances of the output conditioned to any other subset of the concerned input parameters, and the total variance of the output:

$$S_{i_1, ..., i_s} = D_{i_1, ..., i_s}/D.$$

The total Sobol' index of input parameter $X_i$ is defined as the sum of all Sobol' indices of those sets which contain it.

$$S_{T_i} = \sum_{\mathcal{I}_i} D_{i_1, ..., i_s}/D, \\ \mathcal{I}_i = \{\{i_1, ..., i_s\} : \exists k, 1 \le k \le s, i_k = i\}.$$

As all Sobol indices of first order and Sobol indices of higher orders sum up to 1, an alternative definition for the total Sobol index is:

$$S_{T_i} = 1 - S_{\sim i},$$
where $S_{\sim i}$ is the sum of all Sobol indices $S_{i_1, \dots, i_s}$ that do not include index $i$.

Sobol' indices have some nice properties such as explaining explicitly the influence of a parameter on its own (first order Sobol' indices) or in interaction with other parameters (total Sobol' indices). However, estimating them is not always an easy task, especially when it comes to large models with many input parameters and computationally expensive simulations.

## 2.2 Polynomial chaos expansion and the link with Sobol' indices

Except for simple analytical models, Sobol' indices are typically calculated via stochastic methods (Saltelli et al., 2008). Obtaining precise Sobol' indices via these methods requires an efficient exploration of the input parameters' space. One of the main challenges is dealing with large input spaces. Many points are needed to explore a highly dimensional space (curse of dimensionality). As each point is evaluated via the original model, this task becomes difficult when the original model is computationally expensive. In this work, polynomial chaos expansions (PCE) are used to reduce the cardinality of the experimental design (Sudret, 2008).

The model output can be rewritten via its polynomial chaos expansion decomposition as:

$$Y = \sum_{\alpha \in \mathbb{N}^K} y_{\alpha} \Psi_{\alpha}(\mathbf{X}) = \mathcal{M}(\mathbf{X}) \quad (1)$$

where $K$ is the number of input parameters, $\{\Psi_{\alpha}\}_{\alpha \in \mathbb{N}^K}$ is a basis of multivariate orthonormal polynomials, constructed according to the marginal probability density functions $f_{X_i}$ of the input parameters, $\alpha \in \mathbb{N}^K$ are multi-indices, i.e. $\alpha = (\alpha_1, \alpha_2, \dots, \alpha_K)$ where each $\alpha_i$ corresponds to the partial degree with which the $i$th entry parameter $X_i$ is represented in the basis component $\Psi_{\alpha}$ and $y_{\alpha}$ coefficients ($y_{\alpha}$ are the coordinates of the output $Y$ in the new vectorial space defined by the basis $\{\Psi_{\alpha}\}_{\alpha \in \mathbb{N}^K}$).

Once the model $\mathcal{M}$ is rewritten via its polynomial chaos expansion, the Sobol' indices are obtained analytically from the coefficients $y_{\alpha}$, (Sudret, 2008).

Indeed, by adding the squares of all coefficients $y_{\alpha}$ whose multi-index's $\alpha$ only non-null index is the $i$th index $\alpha_i \neq 0$, and dividing it by the total variance of the PCE metamodel $D$, we obtain the Sobol' index of first order of the input parameter $X_i$:

$$S_i = \sum_{\alpha \in \mathcal{I}_i} y_{\alpha}^2 / D, \\ \mathcal{I}_i = \{\alpha \in \mathbb{N}^K : \alpha_i > 0, \alpha_{j \neq i} = 0\}, \\ D = \mathbb{V}ar[\sum_{\alpha \in \mathcal{I}} y_{\alpha} \Psi_{\alpha}(\mathbf{X})] = \sum_{\substack{\alpha \in \mathcal{I} \\ \alpha \neq 0}} y_{\alpha}^2.$$

Similarly, by adding the squares of coefficients $y_{\alpha}$ whose multi-index $\alpha$ satisfies certain properties, the Sobol' indices of a set of input parameters and the total Sobol' indices are also obtained analytically:

$$S_{i_1, \dots, i_s} = \sum_{\alpha \in \mathcal{I}_{i_1, \dots, i_s}} y_{\alpha}^2 / D, \\ \mathcal{I}_{i_1, \dots, i_s} = \{\alpha \in \mathbb{N}^K : k \in \{i_1, \dots, i_s\} \Leftrightarrow \alpha_j \neq 0\}, \\ S_i^T = \sum_{\alpha \in \mathcal{I}_i^T} y_{\alpha}^2 / D, \\ \mathcal{I}_i^T = \{\alpha \in \mathbb{N}^K : \alpha_i > 0\}.$$

Thus, the difficulty and the computational cost of calculating the Sobol' indices for the original model $\mathcal{M}$ can be replaced with the difficulty and the computational cost of finding a polynomial chaos expansion $\mathcal{M}_{PCE}$ which approximates the behavior of the original model "well enough".

Section 2.3 explains how to estimate $\mathcal{M}_{PCE}$, how to limit the cardinality of the experimental design and lastly, how to evaluate if $\mathcal{M}_{PCE}$ is true to the original model.
## 2.3 Estimation of a polynomial chaos expansion

Let the model $\mathcal{M}: \mathbb{R}^K \to \mathbb{R}$ from the beginning of Section 2 be the original model for which we want to build a PCE. Let $f_{X_i}$ be the marginal probability density function of input parameter $X_i$. The goal is to estimate a $\mathcal{M}_{PCE}$ polynomial chaos expansion of the original model by using a limited number $N_{train}$ of original model evaluations:

$$ \mathcal{M}(\mathbf{X}) = Y \approx \sum_{\alpha \in \mathcal{A}_{p,q}^K} y_{\alpha} \Psi_{\alpha}(\mathbf{X}) = \mathcal{M}_{PCE}(\mathbf{X}). $$

The first step in PCE approximation is to choose the truncation scheme for $\alpha$, denoted by $\mathcal{A}_{p,q}^K$. Popular truncation schemes limit the maximum degree with which the input parameters $X_i$ are represented in the basis. This is done by setting an upper limit $p$ such that $\alpha_i < p$ for every $i$. The heuristics behind this truncation being that most of the influence of the input parameters is captured in low orders of the multivariate polynomials (Le Gratiet et al., 2015). For a PCE metamodel with $K$ input parameters and $p$ the maximum degree of the basis polynomial, the cardinality of the basis $card(\mathcal{A}_p^K) = \binom{K+p}{p}$ increases exponentially with respect to the number of input parameters $K$.

To further reduce the basis cardinality, the degree of interaction between the input parameters can be limited. This is done by introducing $q \in (0, 1)$ and imposing the $q$-norm of $\alpha$ to be limited by $p$, i.e. $||\alpha||_q = (\sum_{i=1}^K \alpha_i^q)^{1/q} < p$. Limiting the interactions between input parameters is justified with the sparsity of effects (Le Gratiet et al., 2015). This type of truncation shows a drastic decrease in cardinality of $\mathcal{A}_{p,q}^K$ for a fixed $p$.

Once the truncation scheme is chosen, the creation of the basis $\{\Psi_\alpha\}_{\alpha \in \mathcal{A}_{p,q}^K}$ is straightforward, as it is characterized by the marginal probability density functions $f_{X_i}$. The estimation of the PCE reduces to the estimation of the coefficients $y_\alpha$. To estimate the PCE coefficients $y_\alpha$ comes down to a linear regression on a training sample of pairs of inputs and model outputs $\{(\mathbf{x}_1, \mathcal{M}(\mathbf{x}_1)), (\mathbf{x}_2, \mathcal{M}(\mathbf{x}_2)), \dots, (\mathbf{x}_{N_{train}}, \mathcal{M}(\mathbf{x}_{N_{train}}))\}$.

The number of coefficients $y_\alpha$ to estimate is equal to $card(\mathcal{A}_{p,q}^K)$. In order for the linear regression to be successful, $N_{train}$ should be at least as big as the number of coefficients to estimate. However, the lack of the hypothesis $N_{train} \gg card(\mathcal{A}_{p,q}^K)$ can further be tackled via least angle regression (LARS, Blatman & Sudret (2011)).

Whether the approximated metamodel $\mathcal{M}_{PCE}$ is a good replacement of the original model $\mathcal{M}$ is evaluated via the coefficient of determination $R^2$ on a test set $\{(\mathbf{x}_1, \mathcal{M}(\mathbf{x}_1)), (\mathbf{x}_2, \mathcal{M}(\mathbf{x}_2)), \dots, (\mathbf{x}_{N_{test}}, \mathcal{M}(\mathbf{x}_{N_{test}}))\}$. These estimations are not used for the calculation of the PCE metamodel coefficients $y_\alpha$. The coefficient of determination $R^2$ measures the quality of the prediction of a linear regression:

$$ R^2 = 1 - \frac{\sum_{n=1}^{N_{test}} (\mathcal{M}(\mathbf{x}_n) - \mathcal{M}_{PCE}(\mathbf{x}_n))^2}{\sum_{n=1}^{N_{test}} (\mathcal{M}(\mathbf{x}_n) - \bar{\mathcal{M}})^2}, \quad (2) $$

where $\mathcal{M}(\mathbf{x}_n)$ is the evaluation of the original model for input parameters set to values $\mathbf{x}_n$, $\mathcal{M}_{PCE}(\mathbf{x}_n)$ is the evaluation of the PCE metamodel for input parameters set to values $\mathbf{x}_n$, and $\bar{\mathcal{M}}$ is the mean of all original model evaluations, $\bar{\mathcal{M}} = \sum_{n=1}^{N_{test}} \mathcal{M}(\mathbf{x}_n)$. The closer the $R^2$ is to 1, the better is $\mathcal{M}_{PCE}$ in approximating $\mathcal{M}$.

The estimated Sobol indices are calculated analytically from the coefficients of $\mathcal{M}_{PCE}$:

$$ \hat{S}_i = \sum_{\alpha \in \mathcal{I}_i^*} y_{\alpha}^2 / D, \quad (3) \\ \mathcal{I}_i^* = \{\alpha \in \mathcal{A}_{p,q}^K : \alpha_i > 0, \alpha_{j \neq i} = 0\}. $$

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.4 Sensitivity analysis when model outputs are multidimensional

In previous sections we have seen the definitions of Sobol' indices and their estimation for a model with $K$ independent scalar inputs and a scalar output. However, it often happens that the model whose sensitivity we wish to study produces outputs over a given time period.
Let $\mathcal{M}$ be a model with $K$ independent scalar inputs and a dynamic output over a time interval $\mathcal{T}$:

$$Y(t) = \mathcal{M}(\mathbf{X}, t), \quad t \in \mathcal{T}. \qquad (4)$$

Set $\mu(t) = \mathbb{E}[Y(t)]$ and assume the output random function to be square integrable: $\int_{\mathcal{T}} \mathbb{E}(Y^2(t)) dt < \infty$. We can rewrite the output random function via its principal component (PC) expansion:

$$Y(t) = \mu(t) + \sum_{j} H_{j} \mathbf{v}_{j}(t), \qquad (5)$$

where $\mathbf{v}_j(t)$ is the $j$th functional principal component (also called harmonic or mode), and $H_j \in \mathbb{R}$ is the $j$th functional principal component score (Ramsay & Silverman, 2005).

By decomposing the random variable $Y(t)$ into its principal components (Eq. 5), the variability of $Y(t)$ is rewritten in terms of the variabilities of the principal component scores $H_j$. The variability of the $j$th PC score $H_j$ is denoted with $\lambda_j$. The $\lambda_j$ are decreasing with respect to $j$ by construction.

By truncating the series to the first $J$ PCs we obtain an approximation of the output:

$$Y(t) \approx \mu(t) + \sum_{j=1}^{J} H_{j} \mathbf{v}_{j}(t), \qquad (6)$$

which accounts for $\sum_{j=1}^{J} \lambda_j / \sum_j \lambda_j$ of the total variance in the output $Y(t)$. Oftentimes, $J$ is chosen by first setting the percentage of total inertia observed in the output that we wish to account for. Then, the number of PCs $J$ that should be kept is deduced. Another choice is to directly set the number of PCs that will be studied. As the scores $H_j$ are scalar quantities, all theory developed in previous sections for models with scalar outputs can be applied. Thus, Sobol' indices can be obtained for one PC at a time. Furthermore, in order to have a synthetic sensitivity index for dynamic outputs, rather than sensitivity indices on each PC, the following definition of a generalized sensitivity index is proposed in Lamboni et al. (2011):

$$GSI_{w} = \frac{\sum_{j} \lambda_{j} SI_{w,j}}{\sum_{j} \lambda_{j}}, \qquad (7)$$

where $w$ is the group of the input parameters, $SI_{w,j}$ is the sensitivity index of inputs $w$ on the $j$th PC.

Gamboa et al. (2014) propose another way of aggregating sensitivity indices calculated on individual outputs of a model with multidimensional outputs:

$$ASI_{w} = \frac{\sum_{m=1}^{M} \mathrm{Var}(Y^{(m)}) SI_{w,m}}{\sum_{m=1}^{M} \mathrm{Var}(Y^{(m)})}, \qquad (8)$$

where $Y^{(m)}$ is the $m$th component of the $M$-dimensional output and $SI_{w,m}$ is the sensitivity index of the set of input parameters $w$ on the $m$th component of the output.

The two definitions coincide when the total variability of the output can be described with a finite number of principal components (De Lozzo & Marrel, 2017).

## 3 Case study

The PESHMELBA model simulates water and pesticide transfers at the catchment scale of the Morcille river in the Beaujolais wine region. The heterogeneity of the catchment's soil and vegetation types is represented by dividing the catchment in different plots. In the studied case, the virtual catchment is divided in $M = 14$ plots, each belonging to one of three soil units (SU): SU1 (sandy soil), SU2 (sandy soil on clay) or SU3 (heterogeneous sandy soils), illustrated in Figure 1. The virtual catchment is presented with more details in Rouzies et al. (2023).

PESHMELBA is a three-dimensional model, where the heterogeneity of the soil layers along the vertical axis is also taken into account. Indeed, the soil column under each plot is decomposed in three layers (soil horizons),
<figure><img src="image_2.png" /></figure>

**Figure 1:** The decomposition of the virtual catchment in 14 plots numbered with $m \in \{1, \dots, 14\}$. Each plot belongs to one of the three possible soil units, depicted by different fillings and has one of two possible vegetation types, depicted by the plot's borders (vegetative filter strips have a bold green border, whereas vineyard plots have fine black borders). The separation of the soil column under each plot into three soil horizons (surface, intermediary and deep) is illustrated in the upper right corner. The depths of the soil horizons varies from one SU to another. Image credit: Rouzies et al. (2023).

as in Figure 1. The depth of the soil horizons and their soil properties vary from one SU to another. The soil horizons are considered internally homogeneous when it comes to soil hydrodynamical properties.

PESHMELBA has a modular structure, simulating the physical processes of each soil compartment, then coupling them via the OpenPALM coupler (Buis et al., 2006). The coupling allows for water and pesticide transfers to occur via three types of transfer pathways: infiltration (Richards' equation solved by Ross, 2003), surface runoff (kinematic wave equation) and lateral exchanges (Darcy law). Surface runoff occurs on the surface of the catchment, whereas lateral exchanges refer to subsurface transfers occurring in saturated zones. One may refer to Rouzies et al. (2019) for more details.

Representing the physical properties of each soil horizon of the catchment results in a large number of input parameters. A typical step prior to the calculation of sensitivity indices, is to reduce the number of inputs of interest through a screening step on a small experimental design, before performing the variance-based analysis on a larger one, limited to the selected inputs, (Saltelli et al., 2008). In our case, the number of inputs kept for the calculation of Sobol' sensitivity indices is reduced from 145 to 52 with the elementary effects method (Morris, 1991) (not shown here). The marginal probability distributions of the $K = 52$ input parameters represent the variability of the parameters across the Morcille catchment and are determined from previous studies and available data on the soil, vegetation and pesticide properties of the Morcille catchment, (Frésard, 2010; Peyrard et al., 2016). The distributions are listed in Table 4 and justified in details in Rouzies et al. (2023). Furthermore, the hydrological input parameters are assumed independent (Song et al., 2015).

A soil unit is always made out of the same soil horizons; plots $\{1, 2, 3, 11\} \in SU1$ are all made of the same soil horizons, and the same goes for $\{4, 10, 14\} \in SU2$ and $\{7, 9, 12\} \in SU3$. The previously listed plots correspond to vineyard plots, while the remaining plots $\{5, 6, 8, 13\}$ correspond to vegetative filter strips. A vineyard plot and a vegetative filter strip belonging to the same soil type differ in their surface soil horizon. Indeed, the surface horizon of the vegetative filter strip has a higher infiltration capacity and adapted hydrodynamical properties, while the intermediary and deep horizon maintain the same properties as the vineyard plot.

As an example, the plot 8 and the plot 9 are both of soil type SU2, but 8 is a vegetative filter strip, whereas 9 is a vineyard plot. The physical properties of the deep and intermediary horizons of both plots are described by the input parameters with the suffix *interm_soil2* and *deep_soil2*, respectively. The physical properties of the surface soil horizon of plot 9 are described with the input parameters having the suffix *surf_soil2*; whereas the input parameters describing the properties of surface soil horizon of plot 8 are named with the suffix *surf_soil2(VFZ)*.
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Top row: rainfall time series. The same rainfall time series was used for all simulations and plots. Bottom row: surface moisture outputs of 20 PESHMELBA model simulations with different input parameter values. Different simulations are represented with various shades between orange and purple following no particular order. Outputs are reported for two plots, *m* = 13 on the left and *m* = 4 on the right.</figcaption></figure>

The output studied in this work is the surface moisture series simulated at 2712 consecutive time steps, corresponding to roughly 2500 hours or three months. The output consists in *M* = 14 time series, one for each plot. An illustration of these outputs is provided in Figure 2 for two plots, *m* = 13 and *m* = 4. It can be seen that the dynamics of the outputs is closely driven by rainfall events. Indeed, the peaks in surface moisture are reached right after the rain events.

Furthermore, we observe that most outputs of plot *m* = 13 result in the saturation of the surface moisture starting from around time *t* = 1500h. A study of the saturations of the fourteen plots is reported in Table 1. Two groups are identified. Plots {2, 3, 7, 9, 13} = *G*<sub>sat</sub> show a saturation in water retention during at least 100 hours for a significant part of the simulated outputs. On the other hand, the surface moisture of plots {1, 4, 5, 6, 8, 10, 11, 12, 14} ∉ *G*<sub>sat</sub> do not present long periods of saturation. The definition of the two groups does not affect the methodology, as the same approach is deployed on all plots. Rather, it is introduced for a further explanation of the results, which vary between the two groups.

<table><caption>**Table 1:** Number of simulations (out of *N*<sub>train</sub> = 3000) presenting a saturated surface moisture for more than 100 hours of the simulated period. The surface moisture is considered saturated when equal to the water content at saturation, *thetas_surf*. Figure 9 illustrates 100 simulated outputs of the fourteen plots.</caption><thead><tr><th>m</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th></tr></thead><tbody><tr><td>saturated simul.</td><td>0</td><td>1246</td><td>877</td><td>0</td><td>0</td><td>0</td><td>2943</td><td>0</td><td>3000</td><td>0</td><td>0</td><td>20</td><td>2483</td><td>0</td></tr></tbody></table>

Note that, the rainfall time series used for the simulation is the same across all simulations and plots; its value is fixed to a time series measured on the catchment in a typical winter scenario. Furthermore, as the time steps
are close, we consider the time to be continuous. The model can be written as:

$$ \begin{bmatrix} Y^{(1)}(t) \\ Y^{(2)}(t) \\ \vdots \\ Y^{(M)}(t) \end{bmatrix} = \mathcal{M}(\mathbf{X}, t), \quad t \in \mathcal{T}, \qquad (9) $$

with $M = 14$ the number of plots, $\mathbf{X} \in \mathbb{R}^K$ the $K = 52$ input parameters and $\mathcal{T}$ the considered time period.

# 4 Results and discussion

A space-filling experimental design of $N_{train} = 3000$ points is created over the $K = 52$ dimensional input space. The space-filling design used is a latin hypercube, optimized with the maximin criteria. The PESHMELBA model is evaluated in these points and this sample is used to calculate the first two principal components for each plot. Section 4.1 shows the principal component analysis and interprets the principal components. The same sample is used to build a metamodel for the first two principal components of each plot. A second latin hypercube sample of $N_{test} = 1000$ points is generated on the $K = 52$ dimensional input space for the evaluation of the metamodels. Section 4.2 studies the quality of the metamodels built. The Sobol' indices for the principal components of the dynamic outputs are found analytically from the metamodel coefficients following Section 2.2 and Section 2.3. In Section 4.3 the Sobol' indices are reported and their physical meaning is commented.

## 4.1 Functional principal components

One plot $m \in \{1, \dots, 14\}$ is fixed and its $N_{train} = 3000$ temporal outputs $Y^{(m)}$ are observed. The outputs represent the surface moisture time series of the plot $m$ over the period $\mathcal{T}$. For a fixed plot $m$, the model can thus be written as in Eq. 4. Then, a functional principal component basis is found for the plot $m$. The basis is truncated to the first two components:

$$ Y^{(m)}(t) \approx \mu^{(m)}(t) + H_{1}^{(m)}\mathbf{v}_{1}^{(m)}(t) + H_{2}^{(m)}\mathbf{v}_{2}^{(m)}(t). $$

The part of the total variance explained with the first two principal components PC1 and PC2 is listed in Table 2 per plot. Two groups of plots with similar behavior are identified. One group, with more than 90% of the total output variability explained with PC1 and less than 3% with PC2, and a second group, where the variability explained by PC1 is a bit lower, but still superior to 80%, whereas PC2 accounts for more than 7%. Interestingly, the two groups coincide with $\{2, 3, 7, 9, 13\} \in \mathcal{G}_{sat}$ and $\{1, 4, 5, 6, 8, 10, 11, 12, 14\} \notin \mathcal{G}_{sat}$, the groups identified in Section 3 by a study of their outputs. Furthermore, plots in $\mathcal{G}_{sat}$ have a larger total variability.

<table><caption>**Table 2:** Total variability per plot $\sum_j \lambda_j^{(m)}$ and its percentage captured with the jth PC, i.e. $\lambda_j^{(m)} / \sum_j \lambda_j^{(m)}$, for $j \in \{1, 2\}$.</caption><thead><tr><th>m</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th></tr></thead><tbody><tr><td>Σ<sub>j</sub> λ<sub>j</sub><sup>(m)</sup></td><td>1.85</td><td>2.36</td><td>2.10</td><td>1.94</td><td>1.87</td><td>2.04</td><td>2.21</td><td>1.89</td><td>2.66</td><td>1.99</td><td>1.87</td><td>1.85</td><td>2.80</td><td>1.83</td></tr><tr><td>PC1</td><td>0.99</td><td>0.82</td><td>0.89</td><td>0.96</td><td>0.99</td><td>0.95</td><td>0.85</td><td>0.99</td><td>0.81</td><td>0.96</td><td>0.99</td><td>0.94</td><td>0.83</td><td>0.97</td></tr><tr><td>PC2</td><td>&lt;0.01</td><td>0.16</td><td>0.08</td><td>0.02</td><td>&lt;0.01</td><td>0.02</td><td>0.07</td><td>&lt;0.01</td><td>0.10</td><td>0.02</td><td>&lt;0.01</td><td>0.03</td><td>0.13</td><td>0.02</td></tr><tr><td>PC1 + PC2</td><td>0.99</td><td>0.98</td><td>0.97</td><td>0.98</td><td>0.99</td><td>0.97</td><td>0.92</td><td>0.99</td><td>0.91</td><td>0.98</td><td>0.99</td><td>0.97</td><td>0.96</td><td>0.99</td></tr></tbody></table>

As for the interpretations of principal components, PC1 scores ($H_1$) show a strong positive correlation with the output at all times $t$, and relatively little variation across time (Figure 4). Such behavior is observed across all plots, with correlations not descending below 0.75, except for some time steps coinciding with the saturated periods of plots. PC1 can thus be interpreted as the global amount of surface moisture during the whole simulated period. Indeed, the perturbation on the mean in the direction of the eigenfunction $v_1(t)$ (upper row of Figure 3) manifests as a vertical shift of the whole output time series.
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Perturbations on the first two PCs of plots $m = 4 \notin G_{sat}$ on the left column and $m = 13 \in G_{sat}$ on the right column. The top row shows the perturbations in the direction of PC1, the bottom row shows the perturbations on PC2. The mean of the sample outputs is presented with a solid black line. The red dotted line and the blue dashed line represent the output that is three standard deviations away from the mean in the direction of the corresponding PC.</figcaption></figure>

The interpretation of PC2 is different in the two groups. For the group $m \in G_{sat}$, PC2 scores ($H_2$) present a low but positive correlation which stays constant up until around $t = 1500$, where it shows a change of sign. The correlation then stays negative during the next couple of hundreds of hours (or until the end of the simulation for $m = 13$). The period with negative correlation corresponds to the periods where a part of simulations of the indicated plots start to saturate (Figure 9). PC2 is thus interpreted as the contrast between the soil moisture observed during the saturated period and outside the period of saturation. The second group, $m \notin G_{sat}$, shows a succession of positive and negative correlations, where the changes of sign correspond to rainfall events. This represents a contrast between the surface moisture right after rainfall, and the surface moisture after a longer period of time has passed. It can further be interpreted as the surface soil layer to infiltration rate. The interpretations of PC2 are also illustrated as perturbation on the mean in the direction of the eigenfunction $v_2(t)$ (Figure 3). It can be seen that such perturbation results in a contrast during the saturation period for $m = 13$ and a change of amplitude for $m = 4$.

## 4.2 Polynomial chaos expansion

For each plot $m \in \{1, \dots, 14\}$, a PCE is built independently for the first principal component (PC1). A PCE is also built for the second principal components (PC2) of plots $m \in \{1, \dots, 14\} \setminus \{1, 8, 9, 11\}$. Indeed, the plots $\{1, 8, 9, 11\}$ are excluded as their PC1 alone already accounts for more than 99% of the total variability in the training sample (Table 2).

To sum up, the following input-output relations are estimated through PCE:

$$ 
\begin{aligned}
H_1^{(m)} &= M_1^{(m)}(\mathbf{X}) \quad \forall m \in \{1, \dots, 14\} \\
H_2^{(m)} &= M_2^{(m)}(\mathbf{X}) \quad \forall m \in \{1, \dots, 14\} \setminus \{1, 8, 9, 11\}
\end{aligned}
$$

Where $H_1^{(m)}$ (resp. $H_2^{(m)}$) are the scores obtained on PC1 (resp. PC2) of plot $m$ and $\mathbf{X}$ is the vector of the $K = 52$ independent input parameters with marginal distributions reported in Table 4.
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Correlations between the principal component scores $H_j^{(m)}$ and the output soil moisture at time t, for each PC (rows), and for each plot $m$ (curve color). The first column reports the plots $m \in \mathcal{G}_{sat}$ and the second column $m \notin \mathcal{G}_{sat}$.</figcaption></figure>

For the fitting of the metamodel, the truncation scheme $\mathcal{A}_{p,q}^K$ is set with $q = 0.5$ and $p = 10$ LARS is used. All metamodels are fitted and evaluated via UQLab (Marelli & Sudret, 2014; Marelli et al., 2022). Once all metamodels $\mathcal{M}_{j,PCE}^{(m)}$ are obtained, their coefficients of determination $R^2$ are calculated on a test set of $N_{test} = 1000$ samples.

<table><caption>**Table 3:** Coefficients of determination $R^2$ of the metamodels $\mathcal{M}_j^{(m)}$ built on the $j$th principal component of plot $m$ for $j \in \{1, 2\}$.</caption><thead><tr><th>m</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th></tr></thead><tbody><tr><td>R<sup>2</sup> on M<sub>1,PCE</sub><sup>(m)</sup></td><td>1.00</td><td>0.99</td><td>0.99</td><td>0.99</td><td>0.99</td><td>0.99</td><td>0.99</td><td>1.00</td><td>0.99</td><td>0.99</td><td>1.00</td><td>0.99</td><td>0.98</td><td>0.99</td></tr><tr><td>R<sup>2</sup> on M<sub>2,PCE</sub><sup>(m)</sup></td><td>NA</td><td>0.80</td><td>0.77</td><td>0.99</td><td>NA</td><td>0.97</td><td>0.82</td><td>NA</td><td>0.80</td><td>0.98</td><td>NA</td><td>0.95</td><td>0.78</td><td>0.99</td></tr></tbody></table>

The coefficients $R^2$ calculated on the test set are reported in Table 3. The values of $R^2$ are superior to 0.99 for all plots on the PC1. When it comes to the coefficients $R^2$ of the metamodels built on the PC2, $\mathcal{M}_{2,PCE}^{(m)}$, they are slightly lower for some plots. More in particular, the parcels $\{4, 6, 10, 12, 14\} \notin \mathcal{G}_{sat}$ have an $R^2 > 0.95$, whereas the parcels from $\mathcal{G}_{sat} = \{2, 3, 7, 9, 13\}$ have lower $R^2$ values, ranging from 0.77 for plot $m = 3$ up to 0.82 for plot $m = 7$.

These lower values on $\mathcal{M}_{2,PCE}^{(m)}$ can be explained by the choice of our methodology and the nature of the PC2. Indeed, plots belonging to group $\mathcal{G}_{sat}$ present mostly outputs that reach a saturation point, resulting in a plateau, at some point of the simulation period. However, not all simulations of the sample have this behavior, some simulations do not reach a saturation point, as it can be seen on Figure 2 for plot $m = 13$. Calculating the principal components of such outputs, while using the $L^2$ norm as the distance measure, results in a separation on the PC2 for the two types of behavior. Two distinct clusters are thus observed, one containing the simulations reaching a plateau (lower scores on the PC2), while the other containing simulations with no plateaus (scoring higher on the PC2). However, the cluster of simulations with lower scores on the PC2 has a larger cardinality (2800) than the cluster with the higher scores (200). Thus, when the PCE is evaluated on
these scores, it replicates essentially the behavior of the simulations reaching saturation. Figure 5 shows the comparison between the metamodel and the PESHMELBA model outputs. The $M_{PCE}$ estimates well the behavior of lower scores, but it fails in replicating the behavior producing the simulations that do not reach saturation.

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Comparison of the outputs generated by the original model $M_2^{(13)}$ versus the outputs generated by the polynomial chaos metamodel $M_{2,PCE}^{(13)}$ for the PC2 of plot $m = 13$. The test set is made out of $N_{test} = 1000$ points</figcaption></figure>

A slight increase of $R^2$ is achieved when the values of $p$ and $q$ are augmented. In spite of the improvement in the $R^2$, the Sobol' indices calculated in both cases remain qualitatively the same. In addition, no improvement is observed in the precision of their estimates. The results for the Sobol' indices presented in Section 4.3 are listed for the metamodels built for $p = 10, q = 0.5$ which is a good compromise for the numerical cost.

## 4.3 Sobol' indices

Once the metamodels are built as presented in Section 4.2, the Sobol' indices are calculated analytically via Eq. (3), for the PC1 of plots $m \in \{1, \dots, 14\}$ in Section 4.3.1 and for the PC2 of plots $m \in \{1, \dots, 14\} \setminus \{1, 5, 8, 11\}$ in Section 4.3.2.

### 4.3.1 Sobol' indices on the first principal component PC1

The Sobol' indices of total order for the PC1 of each plot of the catchment are reported in Figure 6. First order Sobol' indices are nearly identical to the total order indices, meaning that the direct influence of the input parameters on the PC1 surpasses the influence through interaction with other parameters.

* Soil type SU1 (sandy), PC1. The top row of Figure 6 represents the total Sobol' indices of the plots of soil type SU1 ($m \in \{1, 2, 3, 5, 11\}$). The largest Sobol indices correspond to parameters describing the surface soil horizon properties of the soil type SU1. More in particular, the input parameter *thetasemanual* (resp. *thetasemanual*($VFZ$)) accounts for more than 90% of the total output variance of the surface moisture of vineyard plots (resp. vegetative filter strips). This parameter describes the water content at saturation in a given soil, in this case the surface soil horizon of vineyard plots of soil type SU1 (resp. vegetative filter strips of soil type SU1). The second most influential input parameter, *mn_sens* (resp. *mn_sens*($VFZ$)) is a parameter governing the relation in Van Genuchten soil hydraulic property description of the surface soil horizon of the vineyard plots belonging to soil type SU1 (resp. vegetative filter strips of soil type SU1), (Van Genuchten, 1980).

* Soil types SU2 (sandy soil on clay) and SU3 (heterogeneous sandy soil), PC1. Total Sobol' indices of plots of
<figure><img src="image_7.png" /><figcaption>**Figure 6:** Sobol' indices of total order calculated on PC1 of each plot of the catchment $m \in \{1, ..., 14\}$, with corresponding bootstrap confidence limits. Rows correspond to soil types SU1, SU2 and SU3. The top seven parameters with the highest total Sobol' index are reported for each soil type. The color and the filling of the bars represent the nature of the input parameters. The light, medium and dark brown correspond respectively to those parameters that describe the physical properties of the surface, intermediate and deep soil horizon. The different fillings of the bars correspond to those parameters that describe the soil types SU1, SU2 and SU3. The colors and fillings are in line with those from Figure 1.</figcaption></figure>

soil type SU2 (resp SU3) are reported in the middle row (resp. bottom row) of Figure 6. The physical interpretation of the results is the same as for the soil type SU1. Once again, the output variability is mostly due to the variability of parameters describing the surface soil horizon properties. In particular, the input parameter $\theta_{\text{tas\_surf\_soil2}}$ (resp. $\theta_{\text{tas\_surf\_soil3}}$, $\theta_{\text{tas\_surf\_soil2}}(VFZ)$, $\theta_{\text{tas\_surf\_soil3}}(VFZ)$), describing the water content at saturation in the surface soil horizon of the vineyard plots belonging to SU2 (resp. surface of the vineyard plots belonging to SU3, surface of the vegetative filter stripes belonging to SU2 and the surface of the vegetative filter stripes belonging to SU3) has the largest total Sobol' index ($\approx 0.9$). The second most influential input parameter, $m_{\text{n\_surf\_soil2}}$ (resp. $m_{\text{n\_surf\_soil3}}$, $m_{\text{n\_surf\_soil2}}(VFZ)$, $m_{\text{n\_surf\_soil2}}$) is a parameter governing the relation in Van Genuchten soil hydraulic property description of the surface soil horizon belonging to the corresponding soil type (Van Genuchten, 1980).

We can conclude that, for a fixed plot, the variability of the PC1, interpreted in Section 4.1 as the average surface moisture over the observed time period, can be traced back to the variability of the input parameters describing the hydrodynamical properties (predominantly water content at saturation) of its surface soil horizon.

However, we find some cases which show influence of deeper soil parameters, such as $m \in \{2, 13\}$. Both plots $m = 2$ and $m = 13$ belong to group $G_{\text{sat}}$ and present plateaus in their output time series. The influence of deeper soil parameters is more prominent in PC2 and is commented in Section 4.3.2.

The Sobol' indices presented in Figure 6 were estimated with a high precision, as can be seen from the small standard deviations calculated with a bootstrap of 100 runs. This is in line with the high values of the determination coefficient $R^2$ previously seen in Table 3.

### 4.3.2 Sobol' indices on the second principal component PC2

For the second principal component PC2, the Sobol' indices of total order for all plots except $\{1, 5, 8, 11\}$ are reported in Figure 7. These plots were excluded from this study as their PC2 accounts for less than 1% of the total output variance (Table 2). Two types of behavior can be observed in Figure 7. We can identify plots $m \in \{4, 6, 10, 12, 14\} \notin G_{\text{sat}}$ where surface parameters are most influential, and plots $m \in \{2, 3, 7, 9, 13\} \in G_{\text{sat}}$ where deep parameters govern the output variability.
<figure><img src="image_8.png" /><figcaption>**Figure 7:** Sobol' indices of total order calculated on PC2 of plots $m \in \{2, 3, 7, 9, 12, 6, 4, 10, 14, 13\}$, with corresponding bootstrap confidence limits. Rows correspond to soil types SU1, SU2 and SU3. The top seven parameters with the highest total Sobol' index are reported for each soil type. The color and the filling of the bars represent the nature of the input parameters. The light, medium and dark brown correspond respectively to those parameters that describe the physical properties of the surface, intermediate and deep soil horizon. The different fillings of the bars correspond to those parameters that describe the soil types SU1, SU2 and SU3. The colors and fillings are in line with those from Figure 1.</figcaption></figure>

The dependence on deep soil properties of group $G_{sat}$ communicates that the saturation of the surface moisture is not limited to the surface layers. Rather, the whole soil column under the plot is saturated, thus the surface runoff observed is due to saturation excess overland flow.

We can conclude that the variability of the contrast in surface moisture during the plateau and outside of the plateau for plots $m \in G_{sat}$ can be traced back to the variability of the input parameters describing the hydrodynamical properties of its deep soil horizon. The amplitude of the surface moisture time series for plots $m \notin G_{sat}$) is governed by the input parameters describing the hydrodynamical properties of its surface soil horizon.

The standard deviations calculated with a bootstrap of 100 runs are depicted on Figure 7 are slightly larger than in Figure 6 for plots belonging to $G_{sat}$. This is in line with what was observed previously in Table 3; the determination coefficient $R^2$ was lower on the metamodels built on PC2, especially for plots belonging to $G_{sat}$. Nonetheless, the precision of Sobol' indices is still satisfactory.

### 4.3.3 Aggregation to the catchment scale

To obtain synthetic sensitivity indices of the input parameters over the whole spatio-temporal domain, the Sobol' indices on PC1 and PC2 from Sections 4.3.1 and 4.3.2 are first aggregated temporally using Eq. 7 where $j \in \{1, 2\}$ and for plot $m$: $\lambda_1 = \lambda_1^{(m)}$ and $\lambda_2 = \lambda_2^{(m)}$. The temporally aggregated sensitivity indices (denoted $GSI$) bear close resemblance to the Sobol' indices on PC1 (Figure 6) and are not reported. Indeed, given that on each plot the PC1 accounts for over 80% of the total inertia (Table 2), the Sobol' indices calculated on PC1 are given greater weights in Eq. 7 than PC2. Hence, the influence of deep soil layers on PC2 (Figure 7) is dwarfed.

Finally, the $GSI_w^{(m)}$ are aggregated spatially using Eq. 8 with $M = 14$, $Var(Y^{(m)}) = \lambda_1^{(m)} + \lambda_2^{(m)}$ and $SI_{w,m} = GSI_w^{(m)}$. The obtained aggregated sensitivity indices (ASI) over the whole spatio-temporal domain are reported on Figure 8. They are essentially a summary of the results presented in Sections 4.3.1 and 4.3.2. As PC1 explains over 80% of the total variability of each plot (Table 2), and as $theta\_surf$ has a Sobol' index superior to 0.85 (Figure 6), it comes as no surprise that by doing a weighted average (Eq. 7) and (Eq. 8), the most influential parameters are the $thetas\_surf$ of all soil types and vegetative filter strips. Note that, input
<figure><img src="image_9.png" /><figcaption>**Figure 8:** Aggregated sensitivity indices (ASI). Ten input parameters with the highest ASI are reported. The color and filling of the bars represents the nature of the input parameters as in Figures 1, 6 and 7.</figcaption></figure>

parameters influential on plots with high total inertia are given more weight in the ASI, Eq. 8. In the same fashion, the ranking of the ASI is influenced by the number of plots belonging to each SU.

# 5 Conclusion

In this work, variance-based sensitivity indices are calculated for the surface moisture outputs of the PESH-MELBA model. The sensitivity indices consider both the temporal and the spatial aspect of the output.

Firstly, the most influential input parameters on the output dynamics are obtained for the first two functional principal components of every plot on the catchment. The common aspects of the PCs are studied across plots. PC1 is interpreted as the average surface moisture, while PC2 is interpreted either as the infiltration rate or, in plots that present a long period of saturation of the surface layer, as the contrast between saturated and unsaturated periods. The Sobol' indices show that the average surface moisture and infiltration rate are governed by the surface soil hydrodynamic properties. Furthermore, the contrast between the surface moisture during saturated periods and unsaturated periods is influenced by deep soil hydrodynamic properties. The dependence on deep soil properties reflects the saturation of the whole soil column. Hence, the surface runoff observed on these plots is due to saturation excess overland flow.

Then, these sensitivity indices are aggregated temporally and spatially, enabling a synthetic view of the most influential parameters on the whole spatio-temporal domain. Notably, the surface hydrodynamic parameters have a predominant influence on the catchment's surface moisture. In particular, the input parameters describing the water content at saturation of all soil types constituting the catchment account for more than 80% of the total output variability. The aggregated indices provide a synthetic view of the most influential parameters at the catchment scale, however in order to understand local contributions to the dynamics, the Sobol' indices calculated on the principal components are more informative.

The use of PCE metamodels for the calculation of Sobol' indices lowers significantly the number of PESHMELBA model simulations needed. Indeed, only 3000 model simulations are used to calculate all sensitivity indices with high precision.

Further research could focus on the application of cluster based approaches, such as in Roux et al. (2021), thus distinguishing the inter and intra variability of the simulations that reach surface moisture saturation and those who don't. Additionally, other sensitivity indices, such as the HSIC could be adapted to spatio-temporal outputs and compared to the indices obtained in this work, bringing other information on dependence between inputs and outputs in the model. Also, the hypothesis of independence could be relaxed for some input parameters, and the Sobol' indices obtained here compared to a GSA considering correlated inputs, (Iooss & Prieur, 2019; Il Idrissi et al., 2021).

Lastly, it is important to note that in the performed study, the interpretation of the principal components and the sensitivity indices obtained are dependent on the rain time series, but also on other external conditions (catchment properties, evapotranspiration). Indeed, a significantly lighter rain, would not have necessarily led
to the soil saturation and thus would not have activated the same physical processes and model parameters.
The variabilities of the external conditions could be included in the GSA by considering them as stochastic
inputs, using stochastic GSA (Zhu & Sudret, 2021) or PCE based stochastic metamodels (Lüthen et al., 2023).

440

Declaration of competing interest

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

Acknowledgements

We thank Xujia Zhu, Matieyendou Lamboni and Céline Helbert for their insights which have enriched this work.

445

References

Antoniadis, A., Lambert-Lacroix, S., & Poggi, J.-M. (2021). Random forests for global sensitivity analysis: A selective review. *Reliability Engineering & System Safety*, 206:107312. 10.1016/j.ress.2020.107312.

Blatman, G. & Sudret, B. (2011). Adaptive sparse polynomial chaos expansion based on least angle regression. *Journal of Computational Physics*, 230(6):2345-2367. 10.1016/j.jcp.2010.12.021.

450

Buis, S., Piacentini, A., Déclat, D., & the PALM Group (2006). PALM: a computational framework for assembling high-performance computing applications. *Concurrency and Computation: Practice and Experience*, 18(2):231-245. 10.1002/cpe.914.

Campbell, K., McKay, M. D., & Williams, B. J. (2006). Sensitivity analysis when model outputs are functions. *Reliability Engineering & System Safety*, 91(10-11):1468-1472. 10.1016/j.ress.2005.11.049.

455

Da Veiga, S. (2015). Global sensitivity analysis with dependence measures. *Journal of Statistical Computation and Simulation*, 85(7):1283-1305. 10.1080/00949655.2014.945932.

De Lozzo, M. & Marrel, A. (2017). Sensitivity analysis with dependence and variance-based measures for spatio-temporal numerical simulators. *Stochastic Environmental Research and Risk Assessment*, 31(6):1437-1453. 10.1007/s00477-016-1245-3.

460

Dubreuil, S., Berveiller, M., Petitjean, F., & Salaün, M. (2014). Construction of bootstrap confidence intervals on sensitivity indices computed by polynomial chaos expansion. *Reliability Engineering & System Safety*, 121:263-275. 10.1016/j.ress.2013.09.011.

465

Fajraoui, N., Ramasomanana, F., Younes, A., Mara, T. A., Ackerer, P., & Guadagnini, A. (2011). Use of global sensitivity analysis and polynomial chaos expansion for interpretation of nonreactive transport experiments in laboratory-scale porous media. *Water Resources Research*, 47(2). 10.1029/2010WR009639.

Ferrer Savall, J., Franqueville, D., Barbillon, P., Benhamou, C., Durand, P., Taupin, M.-L., Monod, H., & Drouet, J.-L. (2019). Sensitivity analysis of spatio-temporal models describing nitrogen transfers, transformations and losses at the landscape scale. *Environmental Modelling & Software*, 111:356-367. 10.1016/j.envsoft.2018.09.010.

470

Frésard, F. (2010). Cartographie des sols d'un petit bassin versant en Beaujolais viticole, en appui à l'évaluation du risque de contamination des eaux par les pesticides. Master's thesis, Université de Franche Comté.

Gamboa, F., Janon, A., Klein, T., & Lagnoux, A. (2014). Sensitivity analysis for multidimensional and functional outputs. *Electronic Journal of Statistics*, 8(1). 10.1214/14-EJS895.

Garcia-Cabrejo, O. & Valocchi, A. (2014). Global Sensitivity Analysis for multivariate output using Polynomial Chaos Expansion. *Reliability Engineering & System Safety*, 126:25-36. 10.1016/j.ress.2014.01.005.

475

Il Idrissi, M., Chabridon, V., & Iooss, B. (2021). Developments and applications of Shapley effects to reliability-oriented sensitivity analysis with correlated inputs. *Environmental Modelling & Software*, 143:105115. 10.1016/j.envsoft.2021.105115.
looss, B. & Prieur, C. (2019). Shapley effects for sensitivity analysis with correlated inputs : comparison with sobol' indices, numerical estimation and applications. *International Journal for Uncertainty Quantification*, 9(5):493-514. 10.1615/Int.J.UncertaintyQuantification.2019028372.

480

Lamboni, M., Monod, H., & Makowski, D. (2011). Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models. *Reliability Engineering & System Safety*, 96(4):450-459. 10.1016/j.ress.2010.12.002.

Le Gratiet, L., Marelli, S., & Sudret, B. (2015). Metamodel-based sensitivity analysis: Polynomial chaos expansions and Gaussian processes. pages 1-37. arXiv:1606.04273.

485

Lüthen, N., Marelli, S., & Sudret, B. (2023). A spectral surrogate model for stochastic simulators computed from trajectory samples. *Computer Methods in Applied Mechanics and Engineering*, 406:115875.

Marelli, S., Lüthen, N., & Sudret, B. (2022). UQLab user manual - Polynomial chaos expansions. Technical report, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland. Report UQLab-V2.0-104.

Marelli, S. & Sudret, B. (2014). UQLab: A Framework for Uncertainty Quantification in Matlab. In *Vulnerability, Uncertainty, and Risk*, pages 2554-2563, Liverpool, UK. American Society of Civil Engineers. 10.1061/9780784413609.257.

490

Marrel, A., Perot, N., & Mottet, C. (2015). Development of a surrogate model and sensitivity analysis for spatio-temporal numerical simulators. *Stochastic Environmental Research and Risk Assessment*, 29(3):959-974. 10.1007/s00477-014-0927-y.

Morris, M. D. (1991). Factorial Sampling Plans for Preliminary Computational Experiments. *Technometrics*, 33(2):161-174. 10.1080/00401706.1991.10484804.

495

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. *Reliability Engineering & System Safety*, 195:106737. 10.1016/j.ress.2019.106737.

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. *Reliability Engineering & System Safety*, 211:107522. 10.1016/j.ress.2021.107522.

500

Peyrard, X., Liger, L., Guillemain, C., & Gouy, V. (2016). A trench study to assess transfer of pesticides in subsurface lateral flow for a soil with contrasting texture on a sloping vineyard in Beaujolais. *Environmental Science and Pollution Research*, 23(1):14-22. 10.1007/s11356-015-4917-5.

Ramsay, J. O. & Silverman, B. W. (2005). Functional data analysis. Springer series in statistics. Springer, New York, 2nd edition.

505

Ross, P. J. (2003). Modeling Soil Water and Solute Transport—Fast, Simplified Numerical Solutions. *Agronomy Journal*, 95(6):1352-1361. 10.2134/agronj2003.1352.

Roux, S., Buis, S., Lafolie, F., & Lamboni, M. (2021). Cluster-based GSA: Global sensitivity analysis of models with temporal or spatial outputs using clustering. *Environmental Modelling & Software*, 140:105046. 10.1016/j.envsoft.2021.105046.

510

Rouzies, E., Lauvernet, 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. *Science of The Total Environment*, 671:1144-1160. 10.1016/j.scitotenv.2019.03.060.

Rouzies, E., Lauvernet, 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. *Geoscientific Model Development*, 16(11):3137-3163. 10.5194/gmd-16-3137-2023.

515

Saint-Geours, N., Bailly, J.-S., Grelot, F., & Lavergne, C. (2014). Multi-scale spatial sensitivity analysis of a model for economic appraisal of flood risk management policies. *Environmental Modelling & Software*, 60:153-166. 10.1016/j.envsoft.2014.06.012.

Saltelli, A., Aleksankina, K., Becker, W., Fennell, P., Ferretti, F., Holst, N., Li, S., & Wu, Q. (2019). Why so many published sensitivity analyses are false: A systematic review of sensitivity analysis practices. *Environmental Modelling & Software*, 114:29-39. 10.1016/j.envsoft.2019.01.012.

520
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. 10.1002/9780470725184.

Sobol', I. (2001). Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. *Mathematics and Computers in Simulation*, 55(1-3):271-280. 10.1016/S0378-4754(00)00270-6.

525

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. *Journal of Hydrology*, 523:739-757. 10.1016/j.jhydrol.2015.02.013.

Sudret, B. (2008). Global sensitivity analysis using polynomial chaos expansions. *Reliability Engineering & System Safety*, 93(7):964-979. 10.1016/j.ress.2007.04.002.

530

Van Genuchten, M. T. (1980). A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. *Soil Science Society of America Journal*, 44(5):892-898. 10.2136/sssaj1980.03615995004400050002x.

Xiao, H. & Li, L. (2016). Discussion of paper by Matieyendou Lamboni, Hervé Monod, David Makowski "Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models", Reliab. Eng. Syst. Saf. 99 (2011) 450-459. *Reliability Engineering & System Safety*, 147:194-195. 10.1016/j.ress.2015.10.015.

535

Xu, L., Lu, Z., & Xiao, S. (2019). Generalized sensitivity indices based on vector projection for multivariate output. *Applied Mathematical Modelling*, 66:592-610. 10.1016/j.apm.2018.10.009.

Zhu, X. & Sudret, B. (2021). Global sensitivity analysis for stochastic simulators based on generalized lambda surrogate models. *Reliability Engineering & System Safety*, 214:107815. 10.1016/j.ress.2021.107815.
<figure><img src="image_10.png" /><figcaption>**Figure 9:** The first 100 output time series for all plots *m* = 1..14 (out of 3000 simulations). The number of simulations presenting a saturated soil surface for longer than 100 hours are denoted next to the plot number (previously seen in Table 1). The plots considered belonging to the group *G*<sub>sat</sub> have a bold black frame.</figcaption></figure>
<table><caption>**Table 4:** PESHMELBA parameters selected for the sensitivity analysis. The suffixes surf, interm and deep specify the soil horizon; soil1, soil2 and soil3 specify the soil type (SU1, SU2, SU3) and (VFZ) denotes vegetative buffer strips. Gaussian, lognormal and uniform distributions are N(μ, σ), LN(μ, σ) and U(a, b), thetar is bounded to [0, 1] and the PCE basis accounts for the truncated gaussian, following Marelli & Sudret (2014).</caption><thead><tr><th>Name</th><th>Definition</th><th>Unit</th><th>Distribution</th></tr></thead><tbody><tr><td>thetasSurfSoil1</td><td rowspan="10">water content at saturation</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.3375, 0.0338)</td></tr><tr><td>thetasIntermSoil1</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.3362, 0.0336)</td></tr><tr><td>thetasDeepSoil1</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.2844, 0.0284)</td></tr><tr><td>thetasSurfSoil2</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.3375, 0.0338)</td></tr><tr><td>thetasIntermSoil2</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.3537, 0.0354)</td></tr><tr><td>thetasDeepSoil2</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.4162, 0.0416)</td></tr><tr><td>thetasSurfSoil3</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.3375, 0.0338)</td></tr><tr><td>thetasIntermSoil3</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.3322, 0.0332)</td></tr><tr><td>thetasDeepSoil3</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.316, 0.0316)</td></tr><tr><td>thetasSurfSoil1(VFZ)</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.3375, 0.0338)</td></tr><tr><td>thetasSurfSoil2(VFZ)</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.3375, 0.0338)</td></tr><tr><td>thetasSurfSoil3(VFZ)</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.3375, 0.0338)</td></tr><tr><td>thetarSurfSoil1</td><td rowspan="10">residual water content</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.0372, 0.0093)</td></tr><tr><td>thetarDeepSoil1</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.0661, 0.0165)</td></tr><tr><td>thetarSurfSoil2</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.0372, 0.0093)</td></tr><tr><td>thetarIntermSoil2</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0, 0.0093)</td></tr><tr><td>thetarDeepSoil2</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0, 0.0093)</td></tr><tr><td>thetarSurfSoil3</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.0372, 0.0093)</td></tr><tr><td>thetarDeepSoil3</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.0612, 0.0153)</td></tr><tr><td>thetarSurfSoil1(VFZ)</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.0372, 0.0093)</td></tr><tr><td>thetarSurfSoil2(VFZ)</td><td>[L<sup>3</sup>L<sup>-3</sup>]</td><td>N (0.0372, 0.0093)</td></tr><tr><td>hgSurfSoil1</td><td rowspan="8">Van Genuchten water retention curve parameter</td><td>[-]</td><td>N (-9.69, 0.969)</td></tr><tr><td>hgSurfSoil2</td><td>[-]</td><td>N (-9.69, 0.969)</td></tr><tr><td>hgDeepSoil2</td><td>[-]</td><td>N (-30.18, 3.018)</td></tr><tr><td>hgSurfSoil3</td><td>[-]</td><td>N (-9.69, 0.969)</td></tr><tr><td>hgSurfSoil1(VFZ)</td><td>[-]</td><td>N (-9.69, 0.969)</td></tr><tr><td>hgSurfSoil2(VFZ)</td><td>[-]</td><td>N (-9.69, 0.969)</td></tr><tr><td>hgSurfSoil3(VFZ)</td><td>[-]</td><td>N (-9.69, 0.969)</td></tr><tr><td>mnSurfSoil1</td><td rowspan="10">Van Genuchten water retention curve parameter</td><td>[-]</td><td>N (0.2685, 0.0268)</td></tr><tr><td>mnDeepSoil1</td><td>[-]</td><td>N (0.2274, 0.0227)</td></tr><tr><td>mnSurfSoil2</td><td>[-]</td><td>N (0.2685, 0.0268)</td></tr><tr><td>mnIntermSoil2</td><td>[-]</td><td>N (0.1289, 0.0129)</td></tr><tr><td>mnDeepSoil2</td><td>[-]</td><td>N (0, 0.01)</td></tr><tr><td>mnSurfSoil3</td><td>[-]</td><td>N (0.2685, 0.0268)</td></tr><tr><td>mnDeepSoil3</td><td>[-]</td><td>N (0.1791, 0.0179)</td></tr><tr><td>mnSurfSoil1(VFZ)</td><td>[-]</td><td>N (0.2685, 0.0268)</td></tr><tr><td>mnSurfSoil2(VFZ)</td><td>[-]</td><td>N (0.2685, 0.0268)</td></tr><tr><td>mnSurfSoil3(VFZ)</td><td>[-]</td><td>N (0.2685, 0.0268)</td></tr><tr><td>KsSurfSoil2</td><td rowspan="8">hydraulic conductivity at total saturation</td><td>[cm/h]</td><td>LN (2.6291, 0.198)</td></tr><tr><td>KsIntermSoil2</td><td>[cm/h]</td><td>LN (2.0292, 0.198)</td></tr><tr><td>KsIntermSoil2</td><td>[cm/h]</td><td>LN (1.2206, 0.198)</td></tr><tr><td>KsDeepSoil2</td><td>[cm/h]</td><td>LN (0.3391, 0.198)</td></tr><tr><td>KsIntermSoil3</td><td>[cm/h]</td><td>LN (2.3762, 0.198)</td></tr><tr><td>KsSurfSoil1(VFZ)</td><td>[cm/h]</td><td>LN (2.6884, 0.198)</td></tr><tr><td>KsSurfSoil2(VFZ)</td><td>[cm/h]</td><td>LN (2.6884, 0.198)</td></tr><tr><td>KsSurfSoil3(VFZ)</td><td>[cm/h]</td><td>LN (2.6884, 0.198)</td></tr><tr><td>KxSurfSoil2</td><td rowspan="2">hydraulic conductivity at soil matrix saturation</td><td>[cm/h]</td><td>LN (-2.2926, 0.198)</td></tr><tr><td>KxIntermSoil2</td><td>[cm/h]</td><td>LN (-2.833, 0.198)</td></tr><tr><td>veget_LAImax_1</td><td>maximal Leaf Area Index</td><td>[-]</td><td>U (2, 3)</td></tr><tr><td>plot_hpond</td><td>maximal ponding height of the plots</td><td>[cm]</td><td>U (0.8, 1.2)</td></tr><tr><td>river_di</td><td>distance between the riverbed and the aquifer</td><td>[cm]</td><td>U (120, 180)</td></tr><tr><td>river_ks</td><td>hydraulic conductivity of the riverbed at saturation</td><td>[cm/h]</td><td>LN (2.1268, 0.198)</td></tr></tbody></table>