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

: The migration and fate of pesticides in natural environments is highly complex. At the hillslope scale, the quantiﬁcation of contaminant ﬂuxes 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 ﬂow 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, inﬁltration, surface and subsurface solute transfers, and shallow groundwater ﬂow. 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.


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][2][3][4][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 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][12][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][28][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 sensivivity 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][33][34][35][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.

Model Description
2.1.1.Variably Saturated Flow and Surface Runoff in CATHY CATHY (CATchment HYdrology) 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: where S w (−) is the water saturation (S w = θ θ s ), θ (−) is the volumetric moisture content, θ s is the saturated moisture content or porosity, S s (L −1 ) is the aquifer specific storage, ψ (L) is the pressure head, t (T) is time, (L −1 ) is the gradient operator, K s (LT −1 ) is the saturated hydraulic conductivity, K r (−) is the relative conductivity, η z = (0, 0, 1), z (L) is the vertical coordinate directed upward, q ss (L 3 L −3 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 ψ), 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].

Advective-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: where C (ML −3 ) is the solute concentration, 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 2 T −1 ) is the surface solute diffusivity, and q ts (ML −1 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: where K d (L 3 M −1 ) is the equilibrium constant, C S (MM −1 ) and C (ML −3 ) are the solute concentrations in, respectively, the soil and water, K OC (L 3 M −1 ) is the soil organic carbon coefficient, F OC (MM −1 ) is the weight fraction of organic carbon in soil, and λ (T −1 ) 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].

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.

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.

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: 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 )): The first order sensitivity index of P i on Y is then defined as: 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)): with P ∼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 = j and j in {1, ..., k}) but are not taken into account by the direct effect of both of the parameters: 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 ≥ 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 ≥ 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 .

•
∑ 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 − ∑ i S i reflects the presence of interactions in the model.

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 , ..., 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 × M simulations are necessary, with M = q 2 and q a prime number such that q ≥ 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 boostrapping [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.

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.

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.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 −3 concentration is set (around 0.1 kg•ha −1 ).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.Most of the ouput 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 ouputs 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.

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][59][60].The pdfs of K s , θ s , n, α, 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.The saturated hydraulic conductivity is commonly represented by a lognormal distribution [61][62][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 θ r pdf is set according to field data and analyses from two site-specific studies [56,64].The coefficients of variation of α 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.

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.

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.

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 , α 1 , and α 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.The second-order indices confirm the important effect of n 2 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 2 , 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 2 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 α 1 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 s2 with other parameters for the cumulative volume exited, and of n 2 for breakthrough timing and cumulative surface runoff between V and VFS.

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 , θ s2 , θ 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.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 θ s1 , θ s2 , and α 1 (see Figure 8).

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 2n 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 of 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.

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 θ 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 × 10 −4 ; 1.5 × 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.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.

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 additionnal 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 s2 index from the method of [51] with N = 500 differs form 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. ) and mass balance error for tebuconazole (%) in the three methods: [51] with N = 1000, [51] with N = 500, and [52].

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 boostrapping.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 s 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.

Figure 1 .
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.

Figure 3 .
Figure 3.The 4 August 2004 rainfall event.The dotted horizontal line represents the mean saturated hydraulic conductivity of the vineyard soil.

Figure 4 .
Figure 4. Probability density of mass errors for water and the two solutes diuron and tebuconazole.

Figure 5 .
Figure 5. First-order indices S i (pink), total indices S T (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).

Figure 6 .
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.

Figure 7 .
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.

Figure 8 .
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.

Figure 9 .
Figure 9.Comparison of first-order and total indices for mass balance error for tebuconazole according to the chosen extraction threshold for extreme values.

Table 1 .
Input probability density functions for the model inputs.LN is lognormal, N is normal, T is triangular, and U is uniform.

Table 2 .
Sensitivity analysis methods with corresponding sampling size.