# Sobol Global Sensitivity Analysis of a Coupled Surface/Subsurface Water Flow and Reactive Solute Transfer Model on a Real Hillslope

Laure-An Gatel, Claire Lauvernet, Nadia Carluer, S. Weill, Claudio Paniconi

► To cite this version:

Laure-An Gatel, Claire Lauvernet, Nadia Carluer, S. Weill, Claudio Paniconi. Sobol Global Sensitivity Analysis of a Coupled Surface/Subsurface Water Flow and Reactive Solute Transfer Model on a Real Hillslope. Water, 2020, 12 (121), pp.21. 10.3390/w12010121. hal-02610007

HAL Id: hal-02610007

https://hal.inrae.fr/hal-02610007

Submitted on 16 May 2020

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
Article

# Sobol Global Sensitivity Analysis of a Coupled Surface/Subsurface Water Flow and Reactive Solute Transfer Model on a Real Hillslope

Laura Gatel <fnref n="1" />,<fnref n="2" />*,†, Claire Lauvernet <fnref n="1" />, Nadia Carluer <fnref n="1" />, Sylvain Weill <fnref n="3" /> and Claudio Paniconi <fnref n="2" />,*

Received: 22 November 2019; Accepted: 22 December 2019; Published: 30 December 2019

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

check for
updates

**Abstract:** The migration and fate of pesticides in natural environments is highly complex. At the hillslope scale, the quantification of contaminant fluxes and concentrations requires a physically based model. This class of model has recently been extended to include coupling between the surface and the subsurface domains for both the water flow and solute transport regimes. Due to their novelty, the relative importance of and interactions between the main model parameters has not yet been fully investigated. In this study, a global Sobol sensitivity analysis is performed on a vineyard hillslope for a one hour intensive rain event with the CATHY (CATchment HYdrology) integrated surface/subsurface model. The event-based simulation involves runoff generation, infiltration, surface and subsurface solute transfers, and shallow groundwater flow. The results highlight the importance of the saturated hydraulic conductivity $K_s$ and the retention curve shape parameter $n$ and they reveal a strong role for parameter interactions associated with the exchange processes represented in the model. The mass conservation errors generated by the model are lower than 1% in 99.7% of the simulations. Boostrapping analysis of sampling methods and errors associated with the Sobol indices highlights the relevance of choosing a large sampling size (at least $N = 1000$) and raises issues associated with rare but extreme output results.

**Keywords:** global sensitivity analysis; CATHY model; surface/subsurface hydrology; reactive solute transfer; pesticides; hillslope scale; Sobol indices; uncertainty analysis

## 1. Introduction

Pesticide use in agricultural catchments leads to widespread surface and subsurface water contamination. Given that a sustained decrease in pesticide use necessitates long-term changes in agricultural practices, it is of crucial importance to concurrently limit transfers from fields to rivers. To do this, it is necessary to quantify governing processes and their potential interactions (surface/subsurface, water/soil/solutes, etc.), with the help of physically based modeling. Existing solute transfer models can be roughly classified into two groups. Some models represent comprehensive transfer processes, including preferential flow, reactive processes, and root uptake, at the plot scale [1–5]. These models generally simulate transfers in one dimension (soil column) and do not represent processes that are significant at the catchment scale such as lateral fluxes and surface/subsurface interactions. They also lack connectivity and discontinuity representation, which plays an important role for pesticide transfer. A second group, the integrated surface/subsurface hydrological models (ISSHMs), is less detailed on reactive solutes at the fine scale but resolves surface

<footnote n="1">¹ Irstea, 5 rue de la Doua, 69100 Villeurbanne, France; claire.lauvernet@irstea.fr (C.L.); nadia.carluer@irstea.fr (N.C.)</footnote>

<footnote n="2">² INRS-ETE, Université du Québec, 490 rue de la Couronne, Quebec City, QC G1K 9A9, Canada</footnote>

<footnote n="3">³ LHyGeS, 1 rue Blessig, 67084 Strasbourg, France; sylvain.weill@engees.unistra.fr</footnote>

<footnote>* Correspondence: lau.gatel@gmail.com (L.G.); claudio.paniconi@ete.inrs.ca (C.P.)</footnote>

<footnote>† Current address: Université Laval, 2325 rue de l'Université, Quebec, QC G1V 0A6, Canada.</footnote>
and subsurface fluxes and exchanges at the catchment scale (e.g., CATchment HYdrology (CATHY) [6], HydroGeoSphere [7,8], and ParFlow [9]. Few models as yet allow simulation of detailed reactive transfer processes in a 3D coupled (surface/subsurface) context.

Pesticide transfer modeling at the hillslope scale raises several challenges. Solute flux coupling between the surface and subsurface can be managed with various methods, amongst them, the conceptual exchange layer [10] and the diffusion coefficient [11–13] are largely used. They both require a calibration to give efficient results [10,14]. Parameter setting is also a critical issue. At the hillslope scale, the representation of heterogeneity [15,16] or pore connectivity [17] can have non negligible effects on the output results. These variables often have nonlinear effects with strong interactions, sometimes resulting in unexpected outputs [18]. Whatever numerical coupling and parametrization strategies are adopted, model evaluation is crucial to assess the relevance of the results.

Model evaluation can target several key points: the ability to reproduce field data; computational efficiency; validity of numerical choices against the physics; effect of parameter uncertainty on output variables; and quantification and ranking of the output sensitivity to a specific parameter. Each type of evaluation can be associated with a particular methodology, and a combination of several aspects of model evaluation is recommended to ensure a coherent assessment of the model [19]. Comparison with observed data is the most direct evaluation approach and ensures that the model is able to represent variables measured in the field, based on site representative parameters. However this is not a sufficient validation test, and requires complementary analysis such as the assessment of the impact of numerical choices on output variables through model intercomparison [20,21] and the investigation of the effects of input factor variations on the same output variables via sensitivity analysis.

Sensitivity analysis can be classified as local, screening, or global. The local ("one-at-a-time") strategy (e.g., [22,23]) consists in varying only one input parameter at a time and allows only local variations in one dimension to be characterized, without consideration of parameter interactions [24]. Screening methods are based on elementary effects [25] and qualitatively detect parameter effects and interactions [26]. This strategy is an efficient compromise between comprehensive results and reasonable computing costs, and thus is widely used with complex hydrological models [27–29]. Global sensitivity analysis (GSA) is based on variance decomposition and aims at identifying the contribution of each input parameter to the variance of a model output variable, and its interactions with all other input parameters [19]. Global sensitivity analysis can produce a ranking of input parameters, which can be useful in checking for overparameterization and equifinality [30,31].

Among GSA methods, the Sobol method has been successfully used in hydrological studies based on conceptual or non-coupled models [32–36], but its computational cost is prohibitive for complex distributed models such as ISSHMs, which are thus generally evaluated through screening methods [37]. A rare example of GSA applied to ISSHMs is the study by [38], who used the Sobol method on HGS at the lysimeter scale to investigate flow processes such as evapotranspiration, water content, and seepage dynamics. GSA has not been applied yet to investigate reactive solute transport in addition to flow processes at the hillslope scale in ISSHMs.

The implementation of water quality related processes in ISSHMs allows an integrated and distributed representation of pesticide transfer. The complexity of the processes involved, as well as their coupling (surface/subsurface and water/soil/solute) suggests significant interactions between parameters. The main objective of the present study is to evaluate these interactions and to assess the importance of certain input parameters in surface/subsurface coupling processes, with a focus on processes involved during a rain event and directly after. A global sensitivity analysis is performed on the CATHY model with pesticide transfer [29,39]. This analysis is performed on a virtual hillslope based on a real vineyard in the Beaujolais region in France that includes a vegetative filter strip used to attenuate pesticide migration.

The paper is organized as follows: in Section 2 the CATHY model is described, focusing, in particular, on newly included reactive processes, the surface/subsurface solute transfer coupling procedure, and a new mixing module. The GSA methods and setup are also described. In Section 3,
the uncertainty and sensitivity analysis results are presented. The results and methodology are discussed in Section 4, with particular attention given to the use of two sampling methods and a systematic analysis of Sobol indices along with their associated errors.

## 2. Material and Methods

### 2.1. Model Description

#### 2.1.1. Variably Saturated Flow and Surface Runoff in CATHY

CATHY (CAtchment Hydralogy) is a 3D physically based model for the simulation of surface and subsurface water flow and solute transport [6,39,40]. Subsurface flow in CATHY is described by the Richards equation and surface runoff is described by the diffusive wave equation:

$$S_w S_s \frac{\partial \psi}{\partial t} + \theta \frac{\partial S_w}{\partial t} = \nabla [K_s K_r (\nabla \psi + \eta_z)] + q_{ss} \quad (1)$$

$$\frac{\partial Q}{\partial t} + c_k \frac{\partial Q}{\partial s} = D_h \frac{\partial^2 Q}{\partial s^2} + c_k q_s(h, \psi), \quad (2)$$

where $S_w$ (-) is the water saturation ($S_w = \theta_s$), $\theta$ (-) is the volumetric moisture content, $\theta_s$ is the saturated moisture content or porosity, $S_s$ (L$^{-1}$) is the aquifer specific storage, $\psi$ (L) is the pressure head, $t$ (T) is time, $\nabla$ (L$^{-1}$) is the gradient operator, $K_s$ (LT$^{-1}$) is the saturated hydraulic conductivity, $K_r$ (-) is the relative conductivity, $\eta_z$ = (0, 0, 1), $z$ (L) is the vertical coordinate directed upward, $q_{ss}$ (L$^3$L$^{-1}$T) is a source (positive) or sink (negative) term that includes the exchange fluxes from the surface to the subsurface, Q (L$^3$T$^{-1}$) is the surface discharge, s (L) is the longitudinal coordinate system used to describe the overland network, $c_k$ (LT$^{-1}$) is the kinematic celerity, $D_h$ (L$^2$T$^{-1}$) is the hydraulic diffusivity, $h$ (L) is the surface water level (a state variable that is continuous with the subsurface pressure head $\psi$), and $q_s$ (L$^3$L$^{-1}$T) is the inflow (positive) or outflow (negative) exchange rate from the subsurface to the surface.

The subsurface equation is solved on a regular 2D mesh and replicated vertically to form a 3D tetrahedral mesh. The vertical layers can be of varying thickness. Boundary conditions and atmospheric forcing can be dynamically prescribed. The subsurface 3D equation (1) is discretized in space with the linear Galerkin finite element method and in time with a weighted finite difference scheme. A velocity field reconstruction module [40] allows the correction of nonconservative velocity fields produced by the Galerkin method. The surface Equation (2) is solved with the explicit Muskingum-Cunge algorithm, yielding, for each rectangular surface cell, the discharge Q and the source term $q_s$. The drainage network, and consequently the surface flow direction s within each cell, is obtained via pre-processing of the digital elevation model (DEM) for the catchment [41]. Flow coupling across the land surface/subsurface interface is based on a boundary condition switching procedure [42] that is mathematically robust and ensures mass conservation [43].

#### 2.1.2. Advection-Reactive Transport

Solute transport in CATHY is represented by a 3D advection-dispersion equation in the subsurface and a 1D advection-diffusion equation for the surface:

$$\frac{\partial \theta C}{\partial t} = \nabla \cdot (-\vec{U}C + D\vec{\nabla}C) + q_{tss} \quad (3)$$

$$\frac{\partial Q_m}{\partial t} + c_t \frac{\partial Q_m}{\partial s} = D_c \frac{\partial^2 Q_m}{\partial s^2} + c_t q_{ts}, \quad (4)$$

where C (ML$^{-3}$) is the solute concentration, $\vec{U}$ (LT$^{-1}$) is the Darcy velocity vector, $D$ (L$^2$T$^{-1}$) is the tensor for both dispersion and diffusion, $q_{tss}$ (ML$^{-3}$T$^{-1}$) is a solute mass source (positive) or sink
(negative) term, $Q_m$ (MT⁻) is the solute mass discharge, $c_t$ (LT⁻) is the kinematic solute celerity, $D_c$ (L²T⁻¹) is the surface solute diffusivity, and $q_{ts}$ (ML⁻¹T) is the solute mass inflow (positive) or outflow (negative) exchange rate from the subsurface to the surface. Transport in the subsurface is solved with a time-splitting technique combining flux-limited finite volumes with an explicit time discretization for advection and a finite element-based implicit formulation that solves the dispersive part [39,44]. This implies that the concentration values are transferred from volumes to nodes and back at each time step, which creates a non negligible numerical dispersion. For this study, only the advective part of subsurface transport is considered.

Reactions of active molecules in the soil are considered and are simulated as reversible instantaneous equilibrium sorption and first-order decay:

$$K_d = \frac{C_S}{C} \quad \text{with} \quad K_d = K_{OC} F_{OC} \qquad (5)$$

$$\frac{\partial C}{\partial t} = -\lambda C, \qquad (6)$$

where $K_d$ (L³M⁻¹) is the equilibrium constant, $C_S$ (MM⁻¹) and C (ML⁻³) are the solute concentrations in, respectively, the soil and water, $K_{OC}$ (L³M⁻¹) is the soil organic carbon coefficient, $F_{OC}$ (MM⁻¹) is the weight fraction of organic carbon in soil, and $\lambda$ (T⁻¹) is the decay constant, which is inversely proportional to the half-life of the decaying solute. These equations are known to adequately represent the major pesticide reactions in soils [45,46]. In the model representation, the reaction calculations are performed after each advection time step [29].

### 2.1.3. Mixing Module

Surface/subsurface coupling for the transport module in CATHY is based on advective exchange [39], which does not account for all exchange processes, in particular solute remobilisation from soil to surface runoff. A conceptual mixing module is thus implemented to handle such processes, which arise either from diffusive transport between the surface and subsurface or from solute adsorption/desorption at the soil surface [47,48]. Indeed it is not possible to simulate surface adsorption in CATHY as the variable that carries solute information is a concentration, therefore solute presence is only acceptable when the soil surface is ponded. At the end of each time step, therefore, for each node under ponding, the solute concentration at the interface is set equal between the surface water and the first computational layer (see Figure 1). The mixing subroutine thereby produces a source or sink term for surface transport, and the subsurface concentration is updated, ensuring mass conservation.

<figure><img src="image_3.png" /><figcaption>Figure 1. Example of mixing subroutine activated by an incoming surface runoff: If a surface cell is ponded, solute concentration is set equal between the surface water and the first numerical layer.</figcaption></figure>
## 2.2. Global Sensitivity Analysis

Sensitivity analysis aims at studying how uncertainty in the output of a model can be apportioned to different sources of uncertainty in the model inputs [19]. After selecting the output variables to be studied and the sensitivity analysis method to be used, the input parameters and their uncertainty must be defined properly. For each parameter, a probability density function (pdf) is set, reflecting the potential variability in the simulation context due, for example, to measurement uncertainty or parameter heterogeneity in the domain. The definition of the pdf for input parameters is based, as far as possible, on data, literature, and expertise.

### 2.2.1. Variance Decomposition Method

The variance decomposition of an output, based on a large sampling of input parameters, can be used to determine the effect of each parameter, and their interactions, on the considered output. Sobol first-order indices and total indices are computed for each parameter and correspond, respectively, to the parameter’s direct effect and to its direct effect plus interactions with all other parameters.

Considering k input parameters $P(P_1, P_2, ..., P_{k-1}, P_k)$ of a model resulting in a response Y, the objective of variance-based GSA quantifies how much the variation of each input parameter contributes to the variation of the output. The variance of Y can be decomposed as the sum of direct effects of the k parameters on Y and the sum of their interactions with the other parameters:

$$V(Y) = V_1 + ... + V_k + V_{\{1,2\}} + V_{\{1,3\}} + ... + V_{\{k-1,k\}} + V_{\{1,2,3\}} + ... + ... + V_{\{1,...,k\}}, \quad (7)$$

with $V_i = V(E(Y|P_i))$ the variance part assigned to the direct effect of $P_i$ and $V_{\{i,j\}}$ the variance part assigned to the interaction between the two input parameters ($P_i, P_j$). Two parameters interact if their cumulative effect on the variance of Y is not equal to the sum of their individual direct effects on the variance of Y. The variance of Y can be written as the sum of the variance of Y due to the $P_i$ parameters and the residuals, namely the interactions $E(V(Y|P_i))$:

$$V(Y) = V(E(Y|P_i)) + E(V(Y|P_i)). \quad (8)$$

The first order sensitivity index of $P_i$ on Y is then defined as:

$$S_i = \frac{V_i}{V(Y)} = \frac{V(E(Y|P_i))}{V(Y)}. \quad (9)$$

Similarly, $V(Y)$ can be defined relative to all parameters except $P_i$. Then $S_{Ti}$, the total index of parameter $P_i$, takes into account direct effects (first order sensitivity) and the effect of all interactions between $P_i$ and the other parameters (Equations (10) and (11)):

$$S_{Ti} = S_i + \sum S_{i,j} + ... + \sum S_{1,...,k} \quad (10)$$

$$S_{Ti} = 1 - \frac{V(E(Y|P_{\sim i}))}{V(Y)}, \quad (11)$$

with $P_{\sim i} = (P_1, P_2, ..., P_{i-1}, P_{i+1}, ..., P_k)$. If $S_{Ti}$ is null, it means that the parameter $P_i$ could take any value in its defined domain without affecting the output Y [49].

The second order indices reflect an effect on Y due to $P_i$ and $P_j$ ($i \neq j$ and $j$ in $\{1, ..., k\}$) but are not taken into account by the direct effect of both of the parameters:

$$S_{ij} = \frac{V_{ij}}{V(Y)}. \quad (12)$$

The useful properties of the $S_i$ and $S_{Ti}$ indices are [19]:
* All $S_i$, $S_{Ti}$, and $S_{ij}$ are positive or null.

* $S_i \ge 0$ is the range of the direct effect, and whatever the interactions, $S_i$ reflects the variance decrease if the parameter $P_i$ was fixed.

* $S_{Ti} \ge S_i$, and $S_{Ti} - S_i$ represents the extent of interactions of $P_i$. The equality $S_{Ti} = S_i$ means that $P_i$ does not interact with any other tested parameter.

* If $S_{Ti} = 0$, then $P_i$ can be fixed and does not influence the variance of the output variable, i.e., this variable is not sensitive to the parameter $P_i$.

* $\sum_i S_i$ is equal to 1 if the model is purely additive (i.e., the total variance is a sum of first-order effects, and there is no interactive effect between input parameters). Thus the difference $1 - \sum_i S_i$ reflects the presence of interactions in the model.

### 2.2.2. Sampling Methods

Several GSA methods based on the Sobol method are combined in this study to set optimal information on the indices [24,50,51], for example, to allow the calculation of first order indices and total indices of the chosen input parameters. The sampling is based on two independent samples and combines them by subset. The two initial samples are generated by a random Monte Carlo draw. Among these three studies, the approach used by [51] is chosen as it is able to deal with NA (not assigned) values, in the case where a simulation has not succeeded, which can happen in practice with numerically complex models such as CATHY. Considering *k* input parameters ($P_1, P_2, ..., P_{i_r}, ..., P_{k-1}, P_k$), the size of the sampling is $N = M(k+2)$, with *M* a number between 500 and 2000 [19]. In theory, the higher *N* is, the lower the uncertainty of the computed indices is.

In order to better understand the Sobol total indices $S_T$, second order indices can be necessary. The [52] method computes first-order and second-order indices based on latin hypercube sampling (LHS) and replicated orthogonal arrays (ROA). The LHS and ROA methods ensure that the random sampling covers the entire defined space. In a two-dimensional space, for example, the LHS method takes into account the row and column of each sample point and ensures that no point overlaps the position of any other one. The ROA method takes the sampling another step forward in defining equally probable subspaces: The sampling points respect the LHS concept, and in addition all subspaces are sampled with the same probability density [53].

For the second-order indices, the method is generalised with orthogonal arrays, maximizing the range coverage and minimizing the number of simulations. $N = 2 \times M$ simulations are necessary, with $M = q^2$ and $q$ a prime number such that $q \ge k-1$. The total indices $S_T$ are not computed. For any sampling method, it is possible and recommended to evaluate the precision of the Sobol indices, for example by bootstrapping [54]. The GSA method applied in this study is based on these two complementary methods ([51,52]) and has been computed with R packages [55], including the bootstrapping evaluation.

### 2.3. Model Setup

The simulations are based on the 0.3 ha (20 m wide, 150 m long) St- Joseph experimental hillslope in the Beaujolais region in southeastern France that contains a 25 m long vegetative strip at its downslope end. The St-Joseph site is located in the Morcille catchment, characterized by a predominant vineyard land cover, a granitic bedrock 3 to 6 m deep, and sandy loam soil [56]. The site has been instrumented since 2003 with the aim of gaining a better understanding of contaminant transfer processes and vegetative filter strip efficiency [16,57]. A large number of piezometers are located on the vegetative filter strip (18) and on the vineyard field (1), monitoring groundwater over specific periods of time in correspondence to active projects and field campaigns. In addition, a single pluviometer records rain data.

Based on the hillslope data, the soil characteristics are simplified for the CATHY model into two zones of contrasting properties (see Figure 2): the first horizon of the vegetative filter strip (designated as Zone 1), and the rest of the domain (designated as Zone 2). The domain is discretized into a surface
mesh of 2 m × 2 m cells and a uniform depth of 3 m, subdivided into five layers of varying thickness (0.001 m, 0.249 m, 0.75 m, 1 m, and 1 m from top to bottom). The discretization yields a 3D mesh of 5016 nodes and 22,500 tetrahedral elements. A no-flux boundary condition is applied on all sides (including the bottom) except on the surface, which recieves rain. The outlet surface point is set downstream in the center of the hillslope width.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The St-Joseph hillslope domain as simulated in CATchment HYdrology (CATHY).</figcaption></figure>

## Initial and Boundary Conditions

The 4 August 2004 event (Figure 3) was chosen for the simulations as it is representative of summer rain events in Beaujolais: A short and intense rain producing significant surface runoff. This particular event produced 31 mm of cumulative rainfall in 45 min. The temporal scale allows some process simplifications: Evapotranspiration is neglected, and the bedrock is considered impermeable.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The 4 August 2004 rainfall event. The dotted horizontal line represents the mean saturated hydraulic conductivity of the vineyard soil.</figcaption></figure>

Water level data for the summer of 2004 are not available, therefore piezometric data for 2012, a year that was hydroclimatically similar to 2004, are used to establish the initial conditions. The initial water table is 3 m below the surface in the upslope half of the domain and gradually rises to 1 m below the surface downslope. The initial conditions for the transport model correspond to a pesticide post-treatment scenario: Solute concentration is null everywhere except in the first 10 cm of the vineyard soil, where a 1 g·m⁻³ concentration is set (around 0.1 kg·ha⁻¹). The transfers of two solutes
are analysed: diuron and tebuconazole. These correspond to pesticide substances actually used on the Saint-Joseph vineyard at different times of the year, and they have very different adsorption coefficients. We simulate them on the same event here in order to evaluate the influence of their contrasting properties.

## 2.4. Input Parameters and Output Variables for Sensitivity Analysis

### 2.4.1. Output Variables and Analysed Criteria

Most of the output variables studied are cumulative in order to integrate global processes and to report the global dynamics of the system. The limit between the vineyard and the vegetative filter strip (VFS) is particularly interesting in the global context of buffer strip efficiency, as well as for the model physics, as it constitutes a frontier between two hydrodynamically distinct zones. The analysed outputs are: cumulative surface and subsurface transfers between the vineyard and the VFS, cumulative outlet fluxes and breakthrough time (defined as the time at which 5% of the total exiting volume or mass has exited the system), subsurface stock evolution from initial to final (t = 60 min) state between the vineyard and the VFS, and global mass balance. All variables are examined for water and solutes.

### 2.4.2. Input Probability Density Functions for the Model

The input parameters selected correspond to the most influential subsurface parameters for pesticide transfer as previously determined from a Morris sensitivity analysis on the CATHY model [29], plus some additional parameters linked to surface/subsurface coupling and surface runoff. For the 60 min event simulated here, the effect of pesticide decay is negligible, therefore only adsorption is considered as a reactive process. All parameters vary independently, which is a common hypothesis used in Sobol sensitivity analysis [34,58–60]. The pdfs of $K_s$, $\theta_s$, $n$, $\alpha$, and $K_d$ are distinct for Zone 1 and Zone 2 (see Table 1), and the total number of parameters varying in the GSA is 15.

<table><caption>**Table 1.** Input probability density functions for the model inputs. LN is lognormal, N is normal, T is triangular, and U is uniform.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Parameter<br>Zone 1</th><th colspan="2">Probability Distribution</th></tr><tr><th>Zone 2</th><th></th></tr></thead><tbody><tr><td>Saturated conductivity</td><td>$K_s$ (m·s<sup>-1</sup>)</td><td>LN (-8.817; 0.69)</td><td>LN (-10.652; 0.69)</td></tr><tr><td rowspan="4">Van Genuchten<br>parameters</td><td>$\theta_s$ (-)</td><td>N (0.54; 0.054)</td><td>N (0.42; 0.042)</td></tr><tr><td>$\theta_r$ (-)</td><td colspan="2">N (0.15; 0.0375)</td></tr><tr><td>$n$ (-)</td><td>N (1.46; 0.146)</td><td>N (1.52; 0.152)</td></tr><tr><td>$\alpha$ (m<sup>-1</sup>)</td><td>N (3.2; 0.96)</td><td>N (10; 3)</td></tr><tr><td>Initial groundwater level</td><td>$GW_i$ (m)</td><td colspan="2">T (-0.24; 0.24; 0)</td></tr><tr><td rowspan="2">Equilibrium sorption<br>constant</td><td>$K_d$ diuron (m<sup>3</sup>·g<sup>-1</sup>)</td><td>T (2.44 × 10<sup>-5</sup>; 3.06 × 10<sup>-5</sup>; 2.75 × 10<sup>-5</sup>)</td><td>T (3.9 × 10<sup>-6</sup>; 16.3 × 10<sup>-6</sup>; 5.1 × 10<sup>-6</sup>)</td></tr><tr><td>$K_d$ tebuconazole (m<sup>3</sup>·g<sup>-1</sup>)</td><td>T (2.2 × 10<sup>-6</sup>; 6.8 × 10<sup>-5</sup>; 4.51 × 10<sup>-5</sup>)</td><td>T (4.5 × 10<sup>-7</sup>; 3.41 × 10<sup>-6</sup>; 1.98 × 10<sup>-6</sup>)</td></tr><tr><td>Roughness</td><td>$K_{strickler}$ (m<sup>1/3</sup>·s<sup>-1</sup>)</td><td colspan="2">T (20; 40; 30)</td></tr><tr><td>Minimum ponding height</td><td>$Pond_{min}$ (m)</td><td>U (0.0025; 0.0075)</td><td>-</td></tr><tr><td>Mixing layer thickness</td><td>$Mix$ (m)</td><td colspan="2">U (0.005; 0.015)</td></tr></tbody></table>

The saturated hydraulic conductivity is commonly represented by a lognormal distribution [61–63]. The mean and standard deviation of the $K_s$ probability density functions are based on measurements performed in Zone 1 (17 observations) and Zone 2 (eight observations) in the work of [64]. The porosity and van Genuchten soil hydraulic parameters follow a normal distribution [63,65,66]. The porosity pdf is also based on field data: two and five measurements for Zones 1 and 2, respectively, giving a
coefficient of variation close to those found in the literature [67,68]. The residual porosity $\theta_r$ pdf is set according to field data and analyses from two site-specific studies [56,64]. The coefficients of variation of $\alpha$ and $n$ are set according to the literature [63,65,69].

The piezometric observations reveal that the water table data follows a triangular distribution with an average variation for each piezometer of ±24 cm against the mean during one year. The piezometers are mainly located in the downstream part of the vineyard and in the buffer strip, however, without further information, the initial water table level is applied homogeneously on the entire domain. The adsorption coefficients $K_d$ for the two tested solutes diuron and tebuconazole were measured for various soils and depths on the Morcille catchment [70] and the pdf is assumed to follow a triangular distribution [28,71]. As there is little information in the literature about pdfs for the Strickler land surface roughness coefficient, a triangular distribution is chosen. The upper and lower limits are set according to [72] and correspond respectively to high grass pasture and a mature row crop.

A minimum ponding water height is used for Zone 1 only, representing the dense vegetation in this buffer strip: a uniform distribution from 2.5 mm to 7.5 mm, based on field observations, is used. For the vineyard land surface, which is part of Zone 2, the minimum ponding water height is kept fixed at zero, meaning that overland flow is triggered in the model as soon as the surface becomes saturated. According to [10,73], the mixing layer thickness is, on average, 1 cm. We selected a uniform distribution with a range from 0.5 cm to 1.5 cm. In the CATHY model, the mixing layer is the first numerical layer, therefore the variation of this parameter influences the mesh.

## 3. Results

Table 2 summarizes the different analyses carried out for the study. The methods from [51,52] are used to determine, respectively, the first order and total Sobol indices, and the first order and second order Sobol indices. For each sampling method, two analyses are performed, one for each of the solutes diuron and tebuconazole. In the following sections, most of the results come from the [51] analysis, and are completed in some cases with results from the [52] second order indices. First order indices and total indices results are presented after the convergence of Sobol indices, namely with the extraction of 57 simulations whose results are replaced with NA (not assigned). A discussion of this index stabilization is presented in Section 4.1, and of the various sampling methods and their relevance in this particular simulation context in Section 4.2.

<table><caption>**Table 2.** Sensitivity analysis methods with corresponding sampling size.</caption><thead><tr><td>Indices</td><td>Sampling Method</td><td>M Parameter</td><td>Sampling Size</td></tr></thead><tbody><tr><td>First order and total indices</td><td>[51]</td><td>1000</td><td>17,000</td></tr><tr><td>First order indices</td><td>[52]</td><td>31</td><td>1922</td></tr><tr><td>Second order indices</td><td></td><td>31</td><td>1922</td></tr></tbody></table>

### 3.1. *Uncertainty Analysis*

The effect of parameter uncertainty on output variables can be studied via the probability density functions of these variables, achieved by running CATHY on the Sobol sampling. The uncertainty analysis on the balance error shows that the probability density is not centered at zero for water and that the model tends to an artificial mass diminution (Figure 4). However, all results for water are very accurate and between -0.5% and 1.5% of error. For solutes, mass conservation is accurate as well and the error does not exceed 0.1% in most cases, centered at zero. There is no notable difference between diuron and tebuconazole on this output variable. Generally, the mass balance errors presented Figure 4 highlight the model robustness and the mass conservation in the coupling.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Probability density of mass errors for water and the two solutes diuron and tebuconazole.</figcaption></figure>

## 3.2. *Sensitivity Indices on Hydrodynamical Outputs*

Sensitivity indices of hydrodynamic parameters of the four hydrodynamic output variables show strong interactions between parameters (Figure 5, $S_T - S_i$). As expected, the vineyard saturated conductivity ($K_{s2}$) drives surface transfers between the two zones. This is a major parameter, influencing infiltration capacity and surface runoff generation. $K_{s2}$ is the saturated conductivity of the entire second horizon of the area, and thus influences greatly the subsurface flow, and the cumulative volume exiting the outlet. As the outlet point is at the surface, one process influencing the breakthrough timing is certainly the water table rising downslope. $n_2$, as a retention curve parameter, is particularly influential on this output. The influence of $n_1$ is however negligible, which may suggest that the first horizon (0–0.3 m) is too thin to have an effect on subsurface flow. The evolution in the vineyard subsurface volume between final and initial states integrates the response of both surface and subsurface flow, as illustrated by the parameters playing a role in its variation: mostly $n_2$, but $K_{s1}$, $K_{s2}$, $a_1$, and $a_2$ as well.

For each of the four output variables, the sum of first-order indices is between 0.55 and 0.65, showing that the model is not additive and confirming that interactions between parameters may explain a large part of the variance. However, this nonadditivity could also be due to other parameters that were not taken into account here. Globally, we note a reasonable uncertainty on the indices (less than 0.25) and a strong presence of interactions. The evolution in the vineyard subsurface volume between initial and final states in particular presents important and quite steady total indices for all parameters (except reaction parameters and minimum ponding height). For all these reasons, second-order indices are necessary to deepen the analysis of this output variable, and are presented next.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** First-order indices *S*<sub>*i*</sub> (pink), total indices *S*<sub>*T*</sub> (blue), and 95% confidence intervals for hydrodynamic output variables: cumulative volume exiting the outlet, breakthrough timing at the outlet, volume difference in the vineyard subsurface between final and initial states, and cumulative surface runoff fluxes between the vineyard (V) and the vegetative filter strip (VFS). The suffix “1” refers to zone 1 and the suffix “2” to zone 2 (see Section 2.4.2).</figcaption></figure>

The second-order indices confirm the important effect of *n*₂ on the evolution in the vineyard subsurface volume, interacting clearly with all other parameters (see Figure 6). The sum of first-order and second-order indices for *n*₂, namely 0.33 (Figure 5) and 0.4 (Figure 6), is near 0.79, which is the total index. This means that for this output variable, the parameter *n*₂ contributes mostly to direct effects and second-order interactions, and little to higher-order interactions. To a lesser extent, surface runoff parameters *Pondh*min and *K*str interact, and *Pondh*min and *α*₁ as well, potentially via surface/subsurface exchanges processes. Other parameters do not show significant interactions. Second-order indices on the other three hydrodynamical outputs confirm the clear trend of their first-order indices (Figure 5): major interactions of *K*<sub>s2</sub> with other parameters for the cumulative volume exited, and of *n*₂ for breakthrough timing and cumulative surface runoff between V and VFS.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Second order indices for the subsurface volume difference between each input parameters with all other ones. Grey lines represent 95% confidence intervals associated with each index.</figcaption></figure>

### 3.3. *Sensitivity Indices on Solute Transfer Outputs*

Sensitivity indices on the solute are illustrated only for tebuconazole as the diuron results are very similar (Figure 7). All tested parameters influence the mass balance error at the first order, essentially with a dominance for the hydrodynamical properties of the second zone (first horizon of the vineyard and second horizon of the whole hillslope): $K_{s2}$, $\theta_{s2}$, $\theta_r$, and $n_2$. As expected, mass breakthrough time results are close to volume breakthrough time results (Figure 5). Results for the evolution in the vineyard subsurface mass between initial and final states show that the mixing layer width has a nonnegligible influence. Part of the solutes move from very shallow subsurfaces into ponded water and are advected by surface runoff. The mixing process accelerates solute transfer, and influences the breakthrough timing in interaction with other parameters.

For the mass difference between the final and initial state, and the cumulative surface fluxes at the vineyard-VFS limit (and to a lesser extent for the breakthrough time), the sum of first-order indices is largely less than 1, confirming the importance of interactions in the model. For the mass balance error output variable, the sum of first-order indices is superior to 1. Theoretically this is not possible; however, in practice it can happen when the indices have a large uncertainty, which is the case as shown by the 95% confidence interval computed by bootstrapping. The trend shows that interactions
are weaker for mass balance error output than for the other outputs, although with larger uncertainty on first-order indices for this variable.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** First-order indices $S_i$ (pink), total indices $S_T$ (blue), and 95% confidence intervals for tebuconazole transfer output variables: mass balance error, mass breakthrough time at the outlet, mass difference in VFS subsurface between initial and final state, and cumulative surface mass fluxes between the vineyard (V) and VFS. The suffix "1" refers to Zone 1 and the suffix "2" to Zone 2.</figcaption></figure>

The analysis of second-order indices performed with the [52] method highlights that $K_{s2}$ has most of the second-order effect. It interacts pairwise with all other parameters, but there are also nonnegligible interactions between $\theta_{s1}$, $\theta_{s2}$, and $\alpha_1$ (see Figure 8).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Second order indices for the cumulative surface mass fluxes between the vineyard and VFS between each input parameters with all other ones. Grey lines represent 95% confidence intervals associated to each index.</figcaption></figure>

### 3.4. *Outcomes of the Global Sensitivity Analysis*

Analysed output variables concern mainly the surface water flow and solute transport, and the subsurface up to the outlet. For water flow, as well as for solute transport, the Zone 2 saturated conductivity and Zone 2*n* van Genuchten parameter are the most influential parameters. The *n* parameter is related to the distribution of pores sizes in soils and intervenes in the retention curve shapes along with the saturated conductivity. The processes occurring in the unsaturated zone seem to have a major influence on the analysed outputs of the particular setup. Except for the mass balance, the sum of first order indices for each output variable are far from 1, and the total indices are large and concern all water flow linked parameters. These findings highlight the complexity of transfers at the hillslope scale. The processes occurring at the hillslope scale are interacting, thus it is reasonable to infer that their numerical representations will do the same. For solute transfer, the interpretation of some results is non conclusive due to the influence of extreme outputs generated by non-extreme parameter samplings. In the following discussion this question will be further elaborated, as will the various sensitivity analyses methods used in this study.
# 4. Discussion

## 4.1. On Complex Interactions in ISSHMs

In the analysis of the Sobol GSA, we observed that index calculation and accuracy are very sensitive to the presence or not of extreme values. In order to get stable sensitivity indices, a few simulation results from the [51] analysis have been excluded, so that no index has an uncertainty superior to 0.5. The final list of extracted simulations is the union of extracted simulations for each considered output variable and is composed of 57 simulations, that is 0.33% of the 17,000 simulations. Given the fact that for each output variable, nearly 30 simulations are extracted, we note that extracted lists superimpose well and most of the extracted simulations are problematic for several output variables. Figure 9 illustrates the problem of extreme values.

In this figure, first-order and total indices are computed based on the same sample, but with various extractions of extreme values: on the left, only the most extreme values are set to NA; in the middle, two extreme values are set to NA; and on the right, 57 simulations are set to NA. In the first two graphs, indices are very similar, except for $n_2$ whose first-order and total indices, respectively, change from -0.023 to 0.22 and from 0.95 to 0.35. The uncertainty of $n_2$ indices strongly decreases as well (see error bars in Figure 9). The same type of transformation arises between the case "2 extractions" and "57 extractions" for parameters $K_{s1}$ and $\theta_r$. In the graph where 57 simulations are set to NA, no index uncertainty is superior to 0.1. The extraction of extreme values allows the stabilization of results and a strong reduction of the uncertainty. The parameter setting of the 57 extracted simulations is shown in Figure 10. It shows miscellaneous parameter combinations with no clearly visible pattern.

For example, the mixing layer width (Mix) is distributed over the whole sampling range (uniform distribution ($5.0 \times 10^{-4}$; $1.5 \times 10^{-3}$)). We note that none of the extracted simulations are characterized by a large saturated conductivity in Zone 2 (first horizon of the vineyard and entire second horizon). If we can hazard that a large conductivity in Zone 2 does not generate extreme values, there is too little information to affirm that the reciprocal assertion regarding low saturated conductivity would be true.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison of first-order and total indices for mass balance error for tebuconazole according to the chosen extraction threshold for extreme values.</figcaption></figure>

The reality of solute transfers at the hillslope scale is complex, and even if some simplifications are always needed in models (in this study: omission of diffusion processes, use of first-order reaction equations, simplification of rock base geometry underlying the hillslope, and heterogeneity representation limited to two zones), ISSHMs such as CATHY reflect this complexity with the representation of interacting physical processes. The particular shape of output variables (very narrow distribution with few extreme values distant from the mean) as well as the unclear pattern shown by the parameter settings resulting in extreme values, highlights the difficulty of managing sensitivity analysis with ISSHMs. Yet the same major influential parameters ($K_s$, $n$) stand out for the majority of studied outputs and highlight a logical picture of model behaviour regarding input parameter uncertainty. On the other hand, even while the influence of classical parameters such as $K_s$
and *n* on flow is clearly established, the influence of other parameters involved in long responses may be underestimated due to the short simulation duration.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Parameter details for extracted simulations (57 over 17,000, i.e., 0.33%). Saturated conductivities are represented on a logarithmic scale.</figcaption></figure>

## 4.2. On the Choice of Relevant GSA Methods

First-order indices have been computed with the two sampling methods from [51] (17,000 simulations per solute) and [52] (1922 simulations per solute). An additional sensitivity analysis is performed with the method from [51], but with a medium sample size (*N* = 500, namely 8500 simulations per solute). Figure 11 shows a comparison of obtained results for the cumulative volume exiting the outlet and the mass balance error for tebuconazole. For nonnegligible indices (i.e., superior to 0.05), the three analyses are consistent for the mass balance error; however the hierarchies of the most influential parameters are not identical. For the cumulative volume exiting the system, the *K*<sub>s2</sub> index from the method of [51] with *N* = 500 differs from the two others, and more generally on the two outputs, the index uncertainties are strong for this method, in particular compared to the same sampling method with *N* = 1000. Concerning the low index value, the method of [51] with *N* = 500 is not consistent with the other methods for cumulative volume at the outlet. For the evaluation of complex models, the determination of total indices and of their uncertainty allows quantitative highlighting of parameter interactions. However, the sampling methods associated with total indices entail high computational costs. As Figure 11 shows, the user cannot circumvent these costs by using a smaller *N* value without risking inaccurate results. Concerning first-order indices of the most influential parameters, the sampling method of [52] presents an undeniable advantage in terms of computation cost. However, its accuracy for less influential parameters vis-à-vis the sampling method from [51] is quite low.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison of first-order indices for cumulative volume exiting the outlet (m³) and mass balance error for tebuconazole (%) in the three methods: [51] with *N* = 1000, [51] with *N* = 500, and [52].</figcaption></figure>

**5. Conclusions**

The present study has explored the interactions and the importance of input parameters on ISSHMs integrating pesticide transfer processes. A global sensitivity analysis was performed on the surface/subsurface reactive transfer CATHY model at the hillslope scale. The modelling of pesticide transfers with CATHY in this study was simplified to some extent. However, the simulation context and the studied outputs reflect complex couplings (solute-water-soil, surface/subsurface) not previously examined in the context of ISSHMs incorporating reactive solute transport. Several sampling methods were used, allowing a detailed analysis of the results, and for each sampling method the uncertainties on Sobol indices were computed by bootstrapping. The large number of simulations performed in order to determine Sobol total indices with the [51] method allowed an uncertainty analysis on the model. The performance of CATHY in terms of mass conservation is excellent, with 99.97% of the simulations showing a mass balance error inferior to 1%. The Sobol sensitivity indices show that *K*<sub>s</sub> and *n* most strongly influence hydrological outputs. Concerning solute transfer output variables, the conclusions are more delicate to draw given the influence of extreme results. The CATHY model is nonadditive and strongly subject to parameter interactions, except for the reactive solute parameters and the minimum ponding height.

The three sampling methods ([51] method with *N* = 500 and *N* = 1000 and [52] method) for first-order index results are consistent when the indices are strong. Therefore we conclude that the evaluation of complex models such as CATHY with various sources of interactions is relevant with a robust method such as a Sobol sensitivity analysis. However, this study highlights some challenges concerning large Sobol index uncertainties linked to the sampling size and the presence of some extreme results. Even if the detailed results are valid in a particular simulation context, some guidelines can be provided for more general studies:

* The computation of Sobol index uncertainties, for example by bootstrapping, is very helpful. This allows verification of the relevance of the index value, and it enables a more robust discussion of the results obtained.

* The [52] method provides a manifest computing cost advantage in the determination of first-order indices compared to the [51] sampling method, with a comparable result accuracy.

* The total Sobol indices are powerful tools in a model evaluation, providing quantitative parameter interaction results; however, they have a high computational cost. In the context of the present study, the sample size *N* should be no less than 1000.

On the sensitivity analysis of the CATHY model, further work could involve the integration of
more spatial heterogeneity, parameter interdependence (for example with group sampling), and a long
term simulation. This last suggestion, in particular, may bring into light a wider spectrum of influential parameters, as longer processes would gain in significance, for example, subsurface solute reactions. Depending on future developments of the CATHY model, other processes such as dispersion, diffusion, and flow in connected pores could be investigated as well.

**Author Contributions:** Conceptualization and methodology, L.G. and C.L.; software, C.P. and S.W.; formal analysis and writing—original draft preparation, L.G.; writing—review and editing, C.L., N.C. and C.P.; supervision, N.C. and C.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Larsbo, M.; Jarvis, N. MACRO 5.0: A Model of Water Flow and Solute Transport in Macroporous Soil: Technical Description; Department of Soil Sciences, Swedish University of Agricultural Sciences: Uppsala, Sweden, 2003.

2. Larsbo, M.; Jarvis, N. Simulating solute transport in a structured field soil. *J. Environ. Qual.* **2005**, *34*, 621–634.
[CrossRef]

3. Ahuja, L.; Rojas, K.; Hanson, J. Root Zone Water Quality Model: Modelling Management Effects on Water Quality and Crop Production; Water Resources Publications, Highlands Ranch: Littleton, CO, USA, 2000.

4. Malone, R.W.; Ahuja, L.R.; Liwang, M.; Don Wauchope, R.; Qingli, M.; Rojas, K.W. Application of the Root Zone Water Quality Model (RZWQM) to pesticide fate and transport: An overview. *Pest Manag. Sci.* **2004**, *60*, 205–221. [CrossRef]

5. Simunek, J.; van Genuchten, M.T.; Sejna, M. Development and applications of the HYDRUS and STANMOD software packages and related codes. *Vadose Zone J.* **2008**, *7*, 587–600. [CrossRef]

6. Camporese, M.; Paniconi, C.; Putti, M.; Orlandini, S. Surface-subsurface flow modeling with path-based runoff routing, boundary condition-based coupling, and assimilation of multisource observation data. *Water Resour. Res.* **2010**, *46*, W02512. [CrossRef]

7. Aquanty Inc. HydroGeoSphere User Manual; Aquanty Inc.: Waterloo, ON, Canada, 2013.

8. Brunner, P.; Simmons, C.T. HydroGeoSphere: A fully integrated, physically based hydrological model. *Groundwater* **2012**, *50*, 170–176. [CrossRef]

9. Kollet, S.J.; Maxwell, R.M. Integrated surface–groundwater flow modeling: A free-surface overland flow boundary condition in a parallel groundwater flow model. *Adv. Water Resour.* **2006**, *29*, 945–958. [CrossRef]

10. Gao, B.; Walter, M.T.; Steenhuis, T.S.; Hogarth, W.L.; Parlange, J.Y. Rainfall induced chemical transport from soil to runoff: Theory and experiments. *J. Hydrol.* **2004**, *295*, 291–304. [CrossRef]

11. Ahuja, L. Modeling soluble chemical transfer to runoff with rainfall impact as a diffusion process. *Soil Sci. Soc. Am. J.* **1990**, *54*, 312–321. [CrossRef]

12. Wallach, R.; Shabtai, R. Surface runoff contamination by chemicals initially incorporated below the soil surface. *Water Resour. Res.* **1993**, *29*, 697–704. [CrossRef]

13. Wallach, R.; Grigorin, G.; Rivlin, J. A comprehensive mathematical model for transport of soil-dissolved chemicals by overland flow. *J. Hydrol.* **2001**, *247*, 85–99. [CrossRef]

14. Zhang, X.; Norton, L.; Lei, T.; Nearing, M. Coupling mixing zone concept with convection-diffusion equation to predict chemical transfer to surface runoff. *Trans. ASAE* **1999**, *42*, 987. [CrossRef]

15. Fiori, A.; Russo, D. Travel time distribution in a hillslope: Insight from numerical simulations. *Water Resour. Res.* **2008**, *44*. [CrossRef]

16. Gatel, L.; Lauvernet, C.; Carluer, N.; Paniconi, C. Effect of surface and subsurface heterogeneity on the hydrological response of a grassed buffer zone. *J. Hydrol.* **2016**, *542*, 637–647. [CrossRef]

17. Russo, D. On the effect of connectivity on solute transport in spatially heterogeneous combined unsaturated-saturated flow systems. *Water Resour. Res.* **2015**, *51*, 3525–3542. [CrossRef]

18. Hopp, L.; McDonnell, J. Connectivity at the hillslope scale: Identifying interactions between storm size, bedrock permeability, slope angle and soil depth. *J. Hydrol.* **2009**, *376*, 378–391. [CrossRef]

19. Saltelli, A.; Ratto, M.; Andres, T.; Campolongo, F.; Cariboni, J.; Gatelli, D.; Saisana, M.; Tarantola, S. Global Sensitivity Analysis: The Primer; John Wiley & Sons: New York, NY, USA, 2008.
20. Maxwell, R.M.; Putti, M.; Meyerhoff, S.; Delfs, J.O.; Ferguson, I.M.; Ivanov, V.; Kim, J.; Kolditz, O.; Kollet, S.J.; Kumar, M.; et al. Surface-subsurface model intercomparison: A first set of benchmark results to diagnose integrated hydrology and feedbacks. *Water Resour. Res.* **2014**, *50*, 1531–1549. [CrossRef]

21. Kollet, S.; Sulis, M.; Maxwell, R.M.; Paniconi, C.; Putti, M.; Bertoldi, G.; Coon, E.T.; Cordano, E.; Endrizzi, S.; Kikinzon, E.; et al. The integrated hydrologic model intercomparison project, IH-MIP2: A second set of benchmark results to diagnose integrated hydrology and feedbacks. *Water Resour. Res.* **2017**, *53*, 867–890. [CrossRef]

22. Persicani, D. Pesticide leaching into field soils: Sensitivity analysis of four mathematical models. *Ecol. Model.* **1996**, *84*, 265–280. [CrossRef]

23. Lenhart, T.; Eckhardt, K.; Fohrer, N.; Frede, H.G. Comparison of two different approaches of sensitivity analysis. *Phys. Chem. Earth Parts A/B/C* **2002**, *27*, 645–654. [CrossRef]

24. Saltelli, A.; Annoni, P.; Azzini, I.; Campolongo, F.; Ratto, M.; Tarantola, S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. *Comput. Phys. Commun.* **2010**, *181*, 259–270. [CrossRef]

25. Morris, M.D. Factorial Sampling Plans for Preliminary Computational Experiments. *Technometrics* **1991**, *33*, 161–174. [CrossRef]

26. Campolongo, F.; Cariboni, J.; Saltelli, A. An effective screening design for sensitivity analysis of large models. *Environ. Model. Softw.* **2007**, *22*, 1509–1518. [CrossRef]

27. Zhou, J.; Cheng, G.; Li, X.; Hu, B.X.; Wang, G. Numerical modeling of wheat irrigation using coupled HYDRUS and WOFOST models. *Soil Sci. Soc. Am. J.* **2012**, *76*, 648–662. [CrossRef]

28. Lauvernet, C.; Muñoz-Carpena, R. 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.* **2018**, *22*, 71–87. [CrossRef]

29. Gatel, L.; Lauvernet, C.; Carluer, N.; Weill, S.; Tournebize, J.; Paniconi, C. Global evaluation and sensitivity analysis of a physically based flow and reactive transport model on a laboratory experiment. *Environ. Model. Softw.* **2019**, *113*, 73–83. [CrossRef]

30. Liu, Y.; Gupta, H.V.; Sorooshian, S.; Bastidas, L.A.; Shuttleworth, W.J. Exploring parameter sensitivities of the land surface using a locally coupled land-atmosphere model. *J. Geophys. Res. Atmos.* **2004**, *109*. [CrossRef]

31. Beven, K. A manifesto for the equifinality thesis. *J. Hydrol.* **2006**, *320*, 18–36. [CrossRef]

32. Van Werkhoven, K.; Wagener, T.; Reed, P.; Tang, Y. Characterization of watershed model behavior across a hydroclimatic gradient. *Water Resour. Res.* **2008**, *44*. [CrossRef]

33. Gan, Y.; Duan, Q.; Gong, W.; Tong, C.; Sun, Y.; Chu, W.; Ye, A.; Miao, C.; Di, Z. A comprehensive evaluation of various sensitivity analysis methods: A case study with a hydrological model. *Environ. Model. Softw.* **2014**, *51*, 269–285. [CrossRef]

34. Song, X.; Zhang, J.; Zhan, C.; Xuan, Y.; Ye, M.; Xu, C. Global sensitivity analysis in hydrological modeling: Review of concepts, methods, theoretical framework, and applications. *J. Hydrol.* **2015**, *523*, 739–757. [CrossRef]

35. Pianosi, F.; Beven, K.; Freer, J.; Hall, J.W.; Rougier, J.; Stephenson, D.B.; Wagener, T. Sensitivity analysis of environmental models: A systematic review with practical workflow. *Environ. Model. Softw.* **2016**, *79*, 214–232. [CrossRef]

36. Dai, H.; Chen, X.; Ye, M.; Song, X.; Zachara, J.M. A geostatistics-informed hierarchical sensitivity analysis method for complex groundwater flow and transport modeling. *Water Resour. Res.* **2017**, *53*, 4327–4343. [CrossRef]

37. Srivastava, V.; Graham, W.; Muñoz-Carpena, R.; Maxwell, R.M. Insights on geologic and vegetative controls over hydrologic behavior of a large complex basin—Global Sensitivity Analysis of an integrated parallel hydrologic model. *J. Hydrol.* **2014**, *519*, 2238–2257. [CrossRef]

38. Ghasemizade, M.; Baroni, G.; Abbaspour, K.; Schirmer, M. Combined analysis of time-varying sensitivity and identifiability indices to diagnose the response of a complex environmental model. *Environ. Model. Softw.* **2017**, *88*, 22–34. [CrossRef]

39. Weill, S.; Mazzia, A.; Putti, M.; Paniconi, C. Coupling water flow and solute transport into a physically-based surface–subsurface hydrological model. *Adv. Water Resour.* **2011**, *34*, 128–136. [CrossRef]

40. Scudeler, C.; Putti, M.; Paniconi, C. Mass-conservative reconstruction of Galerkin velocity fields for transport simulations. *Adv. Water Resour.* **2016**, *94*, 470–485. [CrossRef]
41. Orlandini, S.; Moretti, G. Determination of surface flow paths from gridded elevation data. *Water Resour. Res.* **2009**, *45*, [CrossRef]

42. Putti, M.; Paniconi, C. Time step and stability control for a coupled model of surface and subsurface flow. In *Computational Methods in Water Resources: Volume 2; Developments in Water Science; Miller, C.T., Pinder, G.F., Eds.; Elsevier: Amsterdam, The Netherlands, 2004; Volume 55, pp. 1391–1402.* [CrossRef]

43. Sochala, P.; Ern, A.; Piperno, S. Mass conservative BDF-discontinuous Galerkin/explicit finite volume schemes for coupling subsurface and overland flows. *Comput. Methods Appl. Mech. Eng.* **2009**, *198*, 2122–2136. [CrossRef]

44. Mazzia, A.; Putti, M. High order Godunov mixed methods on tetrahedral meshes for density driven flow simulations in porous media. *J. Comput. Phys.* **2005**, *208*, 154–174. [CrossRef]

45. Jury, W.A.; Horton, R. *Soil Physics*, 6th ed.; John Wiley & Sons: New York, NY, USA, 2004.

46. Kohne, J.M.; Kohne, S.; Simunek, J. A review of model applications for structured soils: A) Water flow and tracer transport. *J. Contam. Hydrol.* **2009**, *104*, 4–35. [CrossRef]

47. Kay, P.; Blackwell, P.A.; Boxall, A.B. Transport of veterinary antibiotics in overland flow following the application of slurry to arable land. *Chemosphere* **2005**, *59*, 951–959. [CrossRef] [PubMed]

48. Tang, X.; Zhu, B.; Katou, H. A review of rapid transport of pesticides from sloping farmland to surface waters: Processes and mitigation strategies. *J. Environ. Sci.* **2012**, *24*, 351–361. [CrossRef]

49. Sobol, I.M.; Tarantola, S.; Gatelli, D.; Kucherenko, S.; Mauntz, W. Estimating the approximation error when fixing unessential factors in global sensitivity analysis. *Reliab. Eng. Syst. Saf.* **2007**, *92*, 957–960. [CrossRef]

50. Jansen, M.J. Analysis of variance designs for model output. *Comput. Phys. Commun.* **1999**, *117*, 35–43. [CrossRef]

51. Martinez, J. *Analyse de Sensibilité Globale Par Décomposition de la Variance*; Ondes & Mascot Num. Institut Henri Poincare: Paris, France, 2011.

52. Tissot, J.Y.; Prieur, C. A randomized orthogonal array-based procedure for the estimation of first- and second-order Sobol indices. *J. Stat. Comput. Simul.* **2015**, *85*, 1358–1381. [CrossRef]

53. Iman, R.; Davenport, J.; Zeigler, D. *Latin Hypercube Sampling (Program User's Guide)*; Technical Report; Sandia Labs.: Albuquerque, NM, USA, 1980.

54. Lemaître, P.; Sergienko, E.; Arnaud, A.; Bousquet, N.; Gamboa, F.; Iooss, B. Density modification-based reliability sensitivity analysis. *J. Stat. Comput. Simul.* **2015**, *85*, 1200–1223. [CrossRef]

55. Iooss, B.; Janon, A.; Pujol, G.; Boumhaout, K.; Veiga, S.D.; Delage, T.; Fruth, J.; Gilquin, L.; Guillaume, J.; Le Gratiet, L.; et al. *Sensitivity: Global Sensitivity Analysis of Model Outputs*; R package version 1.15.1; R Package: Madison, WI, USA, 2018. Available online: https://cran.r-project.org/package=sensitivity (accessed on 29 December 2019).

56. Van den Bogaert, R. Typologie des Sols du Bassin Versant de la Morcille, Caractérisation de Leurs Propriétés Hydrauliques et Test de Fonctions de Pédotransfert. Ph.D. Thesis, Université Pierre et Marie Curie (Paris 6), Institut des Sciences et Industries du Vivant et de l'Environnement, Paris, France, 2011.

57. Lacas, J.G.; Carluer, N.; Voltz, M. Efficiency of a grass buffer strip for limiting diuron losses from an uphill vineyard towards surface and subsurface waters. *Pedosphere* **2012**, *22*, 580–592. [CrossRef]

58. Boateng, S. Probabilistic Unsaturated Flow along the Textural Interface in Three Capillary Barrier Models. *J. Environ. Eng.* **2007**, *133*, 1024–1031. [CrossRef]

59. Zhu, J.; Pohlmann, K.; Chapman, J.; Russell, C.; Carroll, R.; Shafer, D. Sensitivity of Solute Advective Travel Time to Porosities of Hydrogeologic Units. *Groundwater* **2010**, *48*, 442–447. [CrossRef]

60. Zhan, C.; Song, X.; Xia, J.; Tong, C. An efficient integrated approach for global sensitivity analysis of hydrological model parameters. *Environ. Model. Softw.* **2013**, *41*, 39–52. [CrossRef]

61. Coutadeur, C.; Coquet, Y.; Roger-Estrade, J. Variation of hydraulic conductivity in a tilled soil. *Eur. J. Soil Sci.* **2002**, *53*, 619–628. [CrossRef]

62. Fox, G.A.; Muñoz-Carpena, R.; Sabbagh, G.J. Influence of flow concentration on parameter importance and prediction uncertainty of pesticide trapping by vegetative filter strips. *J. Hydrol.* **2010**, *384*, 164–173. [CrossRef]

63. Schwen, A.; Bodner, G.; Scholl, P.; Buchan, G.D.; Loiskandl, W. Temporal dynamics of soil hydraulic properties and the water-conducting porosity under different tillage. *Soil Tillage Res.* **2011**, *113*, 89–98. [CrossRef]
64. Lacas, J.G. Processus de Dissipation des Produits Phytosanitaires Dans les Zones Tampons Enherbées: Étude Expérimentale et Modélisation en vue de Limiter la Contamination des Eaux de Surface. Ph.D. Thesis, Montpellier 2, Montpellier, France, 2005.

65. Alletto, L.; Pot, V.; Giuliano, S.; Costes, M.; Perdrieux, F.; Justes, E. Temporal variation in soil physical properties improves the water dynamics modeling in a conventionally-tilled soil. *Geoderma* **2015**, *243–244*, 18–28.
[CrossRef]

66. Dairon, R. Identification des Processus Dominants de Transfert des Produits Phytosanitaires Dans le sol et évaluation de Modèles Numériques Pour des Contextes Agro-Pédo-Climatiques Variés. Ph.D. Thesis, Université Claude Bernard (Lyon), Villeurbanne, France, 2015.

67. Meyer, P.; Rockhold, M.; Gee, G. *Uncertainty Analyses of Infiltration and Subsurface Flow And Transport for SDMP Sites*; Technical Report; Nuclear Regulatory Commission: Washington, DC, USA, 1997.

68. Baroni, G.; Tarantola, S. A General Probabilistic Framework for uncertainty and global sensitivity analysis of deterministic models: A hydrological case study. *Environ. Model. Softw.* **2014**, *51*, 26–34. [CrossRef]

69. Carsel, R.F.; Parrish, R.S. Developing joint probability distributions of soil water retention characteristics. *Water Resour. Res.* **1988**, *24*, 755–769. [CrossRef]

70. Randriambololohasinirina, P. Pesticide Dissipation Properties in Soils of a Wine-Growing Watershed. Ph.D. Thesis, Université Pierre et Marie Curie (Paris 6), Institut des Sciences et Industries du Vivant et de l'Environnement, Paris, France, 2012.

71. Muñoz-Carpena, R.; Fox, G.; Sabbagh, G. Parameter Importance and Uncertainty in Predicting Runoff Pesticide Reduction with Filter Strips. *J. Environ. Qual.* **2010**, *39*, 630–641. [CrossRef]

72. Arcement, G.J.; Schneider, V.R. *Guide for Selecting Manning's Roughness Coefficients for Natural Channels and Flood Plains*; Technical Report; USGS: Washington, DC, USA, 1989.

73. Walter, M.T.; Gao, B.; Parlange, J.Y. Modeling soil solute release into runoff with infiltration. *J. Hydrol.* **2007**, *347*, 430–437. [CrossRef]

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

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).