VARS and HDMR Sensitivity Analysis of Groundwater Flow Modeling through an Alluvial Aquifer Subject to Tidal Effects

1. Introduction

Groundwater flow and transport numerical models are used to simulate the migration of persistent pollutants in contaminated sites. The relationship between the input variables and the outputs of these models over the entire ranges of inputs is non-linear. Numerical models are subject to numerous uncertainties which hinder the groundwater and contaminant transport modeling. Quantifying the impact of these uncertainties is crucial to improve the accuracy of model predictions [1,2]. Uncertainty in groundwater models may arise from different sources, including: (1) The lack of hydrogeological and hydrodynamic data in the study area; (2) Experimental and data measurement errors; (3) Conceptual or mathematical model oversimplification; (4) Heterogeneity of the hydrodynamic parameters; (5) Boundary conditions; (6) Sparse estimations of aquifer properties derived from tests; and (7) Scale effects [3,4,5,6,7,8,9].

Sensitivity analysis provides a powerful tool to evaluate the impact of uncertainties in input parameters on model outputs and identify the most influential parameters. Measures of sensitivity can be local or global. Local sensitivity methods quantify the sensitivity of the model to one-at-a-time changes in parameters around a reference set of parameters [10,11,12,13,14]. On the other hand, Global Sensitivity Analysis (GSA) evaluate the model sensitivity for wide ranges of parameter and quantify also the interactions among input parameters [15,16,17].

Graphical methods are often used in the first steps of the sensitivity analysis. Scatterplots or tow-variable plots of an output versus an input parameter can be useful sometimes for identifying parameter interactions. Some graphical methods are based on the analysis of the cumulative sum of the normalized reordered model output (CUSUNORO) curve plots [18]. CUSUNORO curves provide a compact and fast way to rank the input parameters, identify the sign of the parameter sensitivity, assess the monotonicity of the dependence of the input parameter and the output, and identify nonlinear relationships.

The method of Morris “elementary effects” is a derivative-based method which consists on perturbing each parameter independently and averaging either the partial finite differences [19] or their absolute values [20]. The method of Morris and its variants are often used to identify and discard the least influential parameters. Morris method is less reliable in determining the most relevant parameters [21].

Variance-based methods such as the Sobol method also known as High Dimensional Model Representation, HDMR, are useful to rank the relevance of input parameters [22], quantify their importance and identify parameters having linear additive effects or nonlinear interactive effects [4,23]. The Sobol method has been found useful to study the uncertainties of groundwater flow and solute transport models [3,4]. However, several studies have highlighted some challenges related to variance-based GSA implementation for complex numerical models [4]. Morris and Campolongo methods [19,20] do not address the dependencies between various sources of uncertainty [4]. Furthermore, sampling size and extreme results can hinder the identification of uncertainty by using the Sobol indexes [24]. Variance-based methods require many model simulations with high computational cost to evaluate the variance of the outputs [4,24]. This is especially challenging when considering high-dimensional spatially-distributed inputs such as hydraulic conductivity, areal recharge and boundary fluxes [4]. In addition, variance-based methods assume that the uncertainty of the model outputs is fully characterized by its variance [25]. The ranking of the influential parameters based only on Sobol indexes may exclude important information [26].

Morris and HDMR methods require a sufficiently large number of model simulations to achieve reliable results. Running thousands of simulations of complex numerical models can be challenging. Alternative solutions have been proposed such as the Variogram Analysis of Response Surfaces (VARS). VARS is a GSA method based on the properties of variograms [27,28]. It combines local sensitivities such as those coming from derivative-based Morris methods and variance-based approaches such as HDMR [29]. Therefore, Morris and Sobol methods are particular cases of VARS. The spatial structure of model outputs and the sensitivity analysis across a wide range of scales are characterized from the directional variograms of the model outputs [27]. According to Razavi and Gupta (2016), VARS is between 1 and 2 orders of magnitude more efficient compared to the Sobol method, while still providing consistent results [28]. This efficiency derives from reducing the number of elementary effects used to compute the total-order index in a given number of simulations [30]. The reduction to explore the input space is compensated when the model is dominated by main effects (as most physical models are) by using a star-based sampling design [30]. In addition, VARS provides metrics that accurately estimate Sobol total-order effects [28] and, for high-dimensional complex models, VARS estimators achieve a high performance even for a low number of runs [30].

GSA methods have been used to quantify the output uncertainties caused by parameter uncertainties for hydrogeological models. Malaguerra et al. [31] applied the Morris method to rank the influence of parameters in a tracer test. Zou et al. [32] proposed a surrogate model to improve the computational efficiency of the analysis and to approximate the results of a numerical contaminant transport model. The most influential parameters identified with VARS were the hydraulic conductivity, the recharge and the porosity. The Sobol analysis of a synthetic heterogeneous phreatic aquifer contaminated by chlorobenzene [33] concluded that the contaminant distribution in the aquifer depends mainly on: (1) Transverse coordinate of the contamination source; (2) Porosity; (3) Hydraulic conductivity; (4) Hydraulic gradient; and (5) Longitudinal and transverse dispersivities. Wang et al. [34] presented a GSA of a groundwater flow model of a colluvial landslide. The results highlighted the importance of selecting an appropriate range of input parameters [34].

Using more than a single GSA method is advisable to increase the confidence in the ranking of input parameters [21,35]. Mishra et al. [35] presented the application of two GSA methods to a synthetic groundwater flow model and a groundwater flow and transport model of a nuclear testing site. The two GSA methods provided consistent results and supplemented one another for both models. VARS and HDMR methods were applied recently to a reactive transport model of a high-level radioactive waste repository in a granitic host rock [36]. The study concluded that parameter rankings of both methods are nearly identical for the 5 input parameters.

The Sardas site is near the Sabiñánigo (Huesca, Spain) reservoir and is heavily affected by lindane and its degradation products released from the Sardas landfill which contains solid lindane production wastes and chlorinated organic contaminants forming a dense non-aqueous phase liquid (DNAPL) [37]. The Sabiñánigo reservoir fluctuations produce a tidal effect on the piezometric heads of the alluvial aquifer [38]. Understanding the dynamics of the tidal effect and the aquifer/reservoir interactions is crucial for quantifying groundwater and contaminant fluxes and proposing remediation techniques. Sobral et al. [38] presented a two-dimensional horizontal groundwater flow model through the Gállego alluvial aquifer. The model reproduced the oscillations of the piezometric head in the aquifer caused by the Sabiñánigo reservoir. The groundwater flow model of the Gállego alluvial is subject to uncertainties in aquifer and aquitard parameters, boundary fluxes and recharge rates. The local sensitivity analysis presented by Sobral et al. [38] was useful to identify the most influential parameters for hydraulic heads in the aquifer. However, their sensitivity analysis was limited to the combination of parameter values corresponding to the calibrated conditions and did not account for the interactions among parameters. The limitations of the local sensitivity analysis are overcome here by performing global sensitivity analyses by using VARS and HDMR methods. Input uncertain parameters include (1) Aquifer parameters; (2) Aquitard vertical conductivities; (3) Boundary inflows; 4) Conductances or leakage coefficients for aquifer/river and aquifer/dam interactions and 5) Areal recharge. The outputs include the computed piezometric heads in 3 monitoring wells, calibration metrics, aquifer/reservoir fluxes, and the average groundwater Darcy velocity modulus. HDMR and VARS methods are used to: (1) Identify the most influential input parameters on the model outputs; and (2) Quantify parameter interactions. The paper starts by describing the study area and the groundwater flow model. Then, the global sensitivity methods (VARS and HDMR) are presented, and the input variables and the outputs are listed. GSA results and discussion are presented afterwards. The paper ends with the main conclusions.

2. Materials and Methods

2.1. Methodological Framework

The methodology used in this study is outlined in Figure 1. A conceptual model was developed based on the analysis of the study area, available hydrogeological and hydrodynamic data, the initial and boundary conditions, and material zones as defined in Sobral et al. [38]. A finite element mesh of triangular elements was used to solve the partial differential equations of groundwater flow [38]. Uncertainties in model input aquifer and boundary parameters lead to uncertainties in model outputs. Thus, global sensitivity analysis methods (VARS and HDMR) are employed here to rank the most important input parameters, quantify the contribution of each parameter to the variance of the results, and quantify uncertainties in model outputs. The input parameters were selected first together with their ranges and probability distribution functions (Table 1). N_s simulations of groundwater flow were performed with CORE^2D for the selected combinations of input parameters. The outputs of the flow model were postprocessed to generate the tables containing the N_i input parameters and N_o outputs. A Halton sequence was selected to generate the input parameters for the VARS analysis and a Sobol sequence was adopted for the HDMR analysis. VARS and HDMR results were compared in terms of parameter rankings.

2.2. Site Description

The study area is located in Sabiñánigo (Huesca, Spain) where a lindane-producing factory (INQUINOSA) operated from 1975 to 1992 on the right bank of the Sabiñánigo reservoir [39] (Figure 2). The Sabiñánigo dam was built in the Gállego River course in 1963 to provide enough hydroelectric power to chemical factories [38]. INQUINOSA deposited solid and liquid hexachlorocyclohexane (HCH) wastes in an uncontrolled manner in the Sardas landfill until approximately 1984 [39]. The landfill is located in the left bank of the Gállego river at 500 m from the Sabiñánigo reservoir.

The floodplain of the Gállego River downstream the Sardas landfill is heavily affected by HCH wastes. By the 1980s, the Sardas landfill was completely filled with urban, construction, and industrial solid wastes including lindane production wastes (between 3·10⁴ and 8·10⁴ tons of solid HCH wastes) [39]. Since the landfill lacks a bottom liner system [39], DNAPL and leachates flowed freely from the landfill into the alluvial of the Gállego river until 1995. In addition, during the construction of the road N-330 in the early 1990s, 50 000 m³ of landfill wastes were deposited on the ground surface of the alluvial plain [40]. In 1995, the landfill was sealed superficially with a PEAD sheet and laterally with a front slurry wall. However, the slurry wall does not prevent leachates from flowing into the alluvial aquifer [38,41].

The alluvial of the Gállego river consists of quaternary silts overlying a layer of sands and gravels. The Larrés marls underlie the quaternary sediments of the Gállego river alluvial [40]. A geological profile across the Sabiñánigo reservoir and the Gállego river floodplain is presented by Sobral et al. [38]. The sands and gravels are much more permeable than the silts and the marls [42].

Liquid HCH wastes were detected downstream the Sardas landfill in 2009 [37], prompting hydrogeological and chemical studies aimed at identifying groundwater contaminant sources and proposing treatment options [39]. DNAPL has migrated by gravity through the alluvial due to its high density [43], and is mainly located on top of the marl layer [42] and inside its fractures [44]. This poses a significant risk since the aquifer and the reservoir are partially connected [38].

The layer of alluvial sands and gravels is confined by quaternary alluvial silts. Since its construction, the reservoir has undergone a siltation process [45] which deposited silting sediments and reducing greatly the reservoir capacity [38]. The alluvial silts and silting sediments act as an aquitard and play the role of a barrier for pollutants by retaining and slowing the arrival of contaminants to the reservoir [38]. However, the presence of DNAPL and HCH sorbed in the soil [42] constitutes a persistent source of organic pollutants.

2.3. Conceptual Model

The groundwater flow model presented here focuses on the alluvial aquifer downstream the Sardas landfill with a large presence of solid and liquid HCH wastes. According to the monthly chemical analyses by the Ebro River Authority, the Sabiñánigo reservoir is the main receptor of chlorinated organic contaminants in the Sardas site.

The model assumes that the aquifer recharges through rainfall infiltration, from the surrounding fluvioglacial terraces on the right bank and from the Larrés marls on the left bank (Figure 3). The hydraulic conductivity of the sands and gravels layer is extremely large compared to that of the silts, silting sediments and marl formations. Therefore, groundwater flow takes place mostly through the sands and gravels of the alluvial and discharges into the Gállego river, the Sabiñánigo reservoir and underneath the Sabiñánigo dam in the downstream part of the study area. The tidal effect caused by the reservoir water level fluctuations has a significant effect on water transfer between the alluvial aquifer and the reservoir [38,46]. Aquifer groundwater flows from the aquifer into the reservoir for normal and low reservoir water levels. However, the flow reverses when the reservoir level rises above the piezometric head in the aquifer [38]. The sands and gravels located underneath the Sabiñánigo reservoir are assumed to be confined by the reservoir silting sediments and the alluvial silts.

2.4. Numerical Groundwater Flow Model

The numerical model presented here is an updated version of the groundwater flow model through the Gállego alluvial aquifer reported in Sobral et al. [38]. Contour maps of thickness of the geological formations were slightly improved by incorporating information from recently drilled wells. Leakage coefficients for the aquifer/reservoir interactions were also updated.

The limits of the model domain include the right bank of the Gállego river along the west, the left bank of the river where the Sardas landfill is located, and the Sabiñánigo dam. The northern limit is located about 500 m north of the tail of the Sabiñánigo reservoir (Figure 3). The model was performed with a non-uniform finite element mesh of triangular elements made of 4399 nodes and 8655 elements (Figure 3). Groundwater flow simulation spans 109 days, starting July 15, 2020 and ending on October 31, 2020. Time increments are all equal to 30 minutes. This time increment is equal to the frequency of the automatic reservoir level measurements.

Aquifer inflows include recharge from rainfall infiltration which was estimated with a hydrological water balance model [47], and inflows along the boundaries in both left and right banks, These inflows were simulated with a Neuman condition [38]. River/aquifer, reservoir/aquifer and aquifer/dam interactions were simulated with Cauchy conditions as described by Sobral et al. [38].

2.5. Global Sensitivity Methods

2.5.1. Graphical Methods

Graphical methods provide a simple, compact and informative tool to study global sensitivities. They are especially adequate for communicating results. They provide a good compromise between the complexity of other methods and compactness [48]. Two-variable interaction plots are helpful to identify visually parameter interactions.

The method of cumulative sum of the normalized reordered model output (CUSUNORO) curve plots condenses in a compact form the sensitivity of an output, y, to a set of input parameters [18]. The CUSUNORO curve function, z(i), is the scaled cumulative sum of the ordered residuals of the output, y, [18] and is given by:

$$z\left(i\right)={\displaystyle \frac{1}{\sqrt{n·S_{yy}}}}·\sum _{j=1}^i({\tilde{y}}_{\pi \left(j\right)}-\overline{y})$$

(1)

where $\overline{y}$ and $S_{yy}$ are the sample mean and standard deviation of the outputs, respectively, $n$ is sample size or the number of simulations, and$({\tilde{y}}_{\pi \left(j\right)}-\overline{y})$ is the j-th ordered residual. The function z(i) can be interpreted based on the following terms: (1) A divisor term or scaling factor,$\sqrt{n}$, which guarantees that the cumulative sum is scaled relative to the sample size n; (2) A divisor term $S_{yy}$ (standard deviation) which ensures that z(i) is dimensionless and standardized; and (3) A summation term representing the cumulative sum of the deviations of the ordered output values from the mean $\overline{y}$.

2.5.2. Sobol Method

The High Dimensional Model Representation (HDMR) method evaluates the input-output relationships of complex models with many input parameters [22]. The Sobol method for global sensitivity analysis requires to compute complex quadratures in a high dimensional space. Random and quasi-random rules are usually applied, which involve a set of sample points where model outputs are analyzed. The Sobol sequence was selected for this study [49]. The input parameters are normalized to range from 0 to 1.

The variance of the output is split as the sum of non-negative variance terms, each associated with combinations of input parameters [22]. Sobol indexes quantify the combined contributions of the input parameters to the variance of the output. The first-order sensitivity index or "main effect” index, S_i, measures the relative individual contribution of the i-th input variable on the variance of the output. The second order sensitivity indexes, S_ij, measure the interactions between the i-th and the j-th inputs on the output variable [22]. Both terms, S_i and S_ij, are usually provided by HDMR software packages. The total sensitivity index, or total-effect index or total-order index, S_Ti, is a measure of the total contribution of the i-th input​ to the output variance [50,51,52], including all the interactions with the rest of the inputs of order k ranging from 2 to n. Its expression is given by:

$$S_{Ti}=S_i+\sum _j{S}_{ij}+\sum _{j,k}{S}_{ijk}+\dots \sum _{j,\dots ,n}{S}_{ij\dots ..n}$$

(2)

2.5.3. Variogram Analysis of Response Surfaces

The variogram analysis of response surfaces (VARS) is a method of characterizing the structure and variability of model outputs within the input parameter space [27]. VARS computes several global sensitivity metrics, including Morris and total sensitivity Sobol indexes by using a tailored sampling strategy [27]. VARS star-based sampling strategy requires all input parameters to be normalized to range between 0 and 1. Then, the star centers are selected by using a random sampling strategy. A cross section of equally spaced points with a selected resolution Δh is performed for each star center and for each input variable. Each set of points constitutes a “star”. There are different methods to sample the star centers [53]. The Halton sequence was selected for this study [54]. The outputs are calculated at all points for all stars. Directional variograms and covariograms are calculated along the lines of the input parameters.

Like other GSA methods, VARS requires performing a sufficiently large number of model simulations to obtain estimates of the directional variograms for each input variable. The variograms along the directions of the parameters provide the sensitivity indexes of the GSA. The limit of the directional variogram divided by h² as h approaches 0, corresponds to Campolongo's version of the Morris effects [28]. On the other hand, the limit of the variogram for large h corresponds to the total variance linked to the Sobol total effect or total sensitivity index [30].

Morris and Sobol methods are limiting cases of VARS. The directional variograms of VARS provide metrics that accurately estimate Sobol and Morris indexes [27]. In addition, VARS calculates the Integrated Variogram Across a Range of Scales index (IVARS₅₀). IVARS₅₀ provides a sensitivity index corresponding to intermediate scales and is calculated from the directional variograms according to:

$${IVARS}_{50}={\int }_0^{0.5}\gamma \left(h\right)dh$$

(3)

where $\gamma \left(h\right)$ is the directional variogram and $h$ is the normalized input parameter. It should be recalled that the VARS method requires that all the input parameters are normalized to range from 0 to 1.

2.6. Inputs and Outputs

Input parameters were selected based on the conclusions of the local sensitivity analyses performed by Sobral et al. [38]. The final list of parameters includes also other potentially relevant parameters such as boundary fluxes. The final list of input parameters considered in the GSA includes (Figure 3): (1) The hydraulic conductivities of the aquifer in four material zones which were defined to account for the spatial heterogeneity of the alluvial aquifer (K1, K2, K2 and K4); (2) The aquifer specific storage (S_S); (3) The vertical hydraulic conductivities of the aquitard interposed between the aquifer and the reservoir which includes the silting sediments settled along the former course of the Gállego river (K_Vs1), the silting sediments in the rest of the reservoir (K_Vs2) and the alluvial silts (K_Vs3); (4) The leakage coefficients of the aquifer/river interactions (α_r) and the conductance for the flow underneath the dam (α_d); (5) Aquifer recharge rates from precipitation in confined (r_c) and unconfined (r_u) parts of the aquifer; and (6) Boundary inflows (Q₆, Q₇, Q₉, Q₂ and Q₁). The ranges and statistical distributions of the input parameters are listed in Table 1. Log-uniform probability density functions were used for the hydrodynamic parameters and leakage coefficients, and uniform distributions were adopted for recharge rates and boundary inflows (Table 1).

The outputs of the GSA were selected to gain insight on the sensitivity of model calibration metrics to changes in input parameters and evaluate the time change in sensitivities of hydraulic heads, h, aquifer/reservoir fluxes, Q, and Darcy velocity, q, caused by the tidal effects. The outputs related to calibration metrics include the global mean absolute error in hydraulic heads (MAEg), the global root mean square error in hydraulic heads normalized by the standard deviation of the measured head data (NRMSEg) and the global Nash–Sutcliffe index [55] (NSEg) of the head residuals (measured minus computed heads) at 8 monitoring wells (Figure 4). Other outputs include the computed piezometric heads in wells ST1C, PS19 and SPN1 at 3 selected times, t1, t2 and t3, which correspond to an event of sudden rise of the reservoir water level which took place from September 18, 2020 at 20:30 to September 19, 2020, at 04:30 (Figure 5). The heads in these 3 wells at the 3 selected times are denoted as ST1C_t1, ST1C_t2, ST1C_t3, PS19B_t1, PS19B_t2, PS19B_t3, SPN1_t1, SPN1_t2 and SPN1_t3. Other outputs include the aquifer/reservoir fluxes at times t1, t2 and t3, which are denoted as Q_t1, Q_t2 and Q_t3, and the average modulus of the Darcy velocity near well PS16C where a tracer test was performed (q_av).

The model head residuals were calculated as the differences between measured and computed heads using recorded piezometric head data from ST1B, ST1C, ST2, SPN1, PS16C, PS19B, PS25B, and PS26 wells. The mean absolute error (MAE), the root mean squared error normalized by the standard deviation of the measured data (NRMSE), and the Nash–Sutcliffe index (NSE) were computed for each well (j). The metrics were globalized to evaluate goodness of fit in the eight wells equipped with divers:

$${MAE}_g={\displaystyle \frac{\sum _j{MAE}_j}{N_w}}$$

(4)

$${NRMSE}_g={\displaystyle \frac{\sqrt{\frac{1}{N_T}\sum _j{N}_j{{RMSE}_j}^2}}{{\sigma }_T}}$$

(5)

$${NSE}_g=1+{\displaystyle \frac{\sum _j{N}_j{\sigma }_j^2\left({NSE}_j-1\right)}{N_T·{\sigma }_T^2}}$$

(6)

where:

• $N_w$ is the number of monitoring wells.
• $N_j$ is the number of measured piezometric heads in the j-th well.
• $N_T$ is the total number of measured piezometric heads in all the wells.
• ${\sigma }_j$ is the standard deviation of the measured piezometric heads in the j-th well.
• ${MAE}_j$, ${RMSE}_j$ and ${NSE}_j$ are the mean absolute error, the root mean squared error and the Nash–Sutcliffe index for the j-the well.
• ${\sigma }_T$ is the standard deviation of the measured piezometric heads in all wells.

2.7. Global Sensitivity Simulation Runs

Global sensitivity simulation runs were generated by using two different sequences. A first set of 30 800 runs were prepared for VARS by using the Halton sequence to generate the sequence of star centers [54]. The second sequence of 16 384 runs was obtained by using the Sobol method [49]. The runs were performed at the Galician Supercomputing Center (CESGA).

2.8. Software

Model simulations were performed with CORE^2D V5, a finite element code for transient saturated and unsaturated water flow and heat transport in heterogeneous and anisotropic porous and fractured media [56]. The code solves for groundwater flow and solute transport by using Galerkin triangular finite elements and a Euler time discretization scheme. CORE^2D V5 has undergone extensive verification against analytical solutions and has been benchmarked against other reactive transport codes [57,58]. CORE^2DV5 has been used extensively for modeling groundwater flow and solute transport in aquifers and lab and in situ experiments [59,60].

The single processor code CORE^2D V5 was compiled at the CESGA by using the FinisTerrae-III supercomputer [61,62]. This advanced computing infrastructure supports the concurrent execution of hundreds of simultaneous and parallel model runs. The simulation of the VARS-Halton and Sobol sequences runs took 6 days of computing wall time. The same task would have taken more than three years running on a personal computer.

VARS version 2.1 was run under MATLAB® version 2022a in an UBUNTU 20.04.6 LTS system. HDMR computations were carried out by using both SALib [63] and GUI-HDMR V1.1 [64].

3. Results and Discussion

3.1. Groundwater Flow Model Results

The computed hydraulic gradient at the Sardas site downstream the landfill is very small due to the large hydraulic conductivity of the layer of sands and gravels, and the small groundwater inflows [38]. The daily fluctuations of the reservoir water level induce a tidal effect on the piezometric heads of the aquifer. The fluctuations in the aquifer are damped and delayed because the aquitard acts as a dumping barrier between the reservoir and the aquifer [38]. The amplitude of the oscillations of the piezometric heads, A_h and the time lag, t_R, are very sensitive to the specific storage coefficient of the layer of sands and gravels (S_S), and the vertical hydraulic conductivities of the aquitard (K_Vs1, K_Vs2 and K_Vs3) [38].

The reservoir water level oscillates between 764.31 m and 765.43 m daily in the modelling period. The average water level in the reservoir (764.80 m) is slightly smaller than the piezometric head in the aquifer (764.90 m). Groundwater generally flows from the aquifer to the silting sediments and the alluvial silts of the reservoir. However, when the reservoir water level rises above the piezometric head in the aquifer, the Sabiñánigo reservoir recharges the aquifer, resulting in a rise in piezometric heads. The modulus of the Darcy velocity in the aquifer decreases during these high-reservoir water level events. If the duration of the high water level is long enough, the direction of groundwater flow may reverse near the reservoir [38,41].

3.2. GSA Results for the Groundwater Flow Model of the Gállego Alluvial Aquifer

3.2.1. Graphical Methods

Two-variable scatterplots are useful to identify and illustrate the interactions of two input parameters. Figure 6 shows two-variable scatterplots of the computed heads in 3 monitoring wells s (outputs ST1C_t2, PS19B_t2 and SPN1_t2) and the aquifer/reservoir flow (Q_t2) at time t2 corresponding the peak reservoir water level. The outputs are plotted versus the vertical conductivity of the aquitard (K_Vs1). These plots correspond to the 16 384 runs of the Sobol sequence. The clouds of plots are shown for the following three ranges of percentiles, p, of the specific storage coefficient (S_S): 1) p < 30%; 2) 30% < p < 70% and 3) p > 70%.

The scatterplots show that the computed heads and the aquifer/reservoir flow at time t2 depend on S_S. Generally, the computed piezometric heads in wells ST1C and SPN1 at t2 are low for large S_S (p > 70%). However, simulations corresponding to extreme values of other input parameters result in high piezometric heads even for high values of S_S. On the other hand, the computed heads are higher when S_S is in the lower percentile (p < 30%). It should be noticed that the well PS19B is located further away from the reservoir, and thus, its piezometric head is less affected by the interaction between K_Vs1 and S_S.

The scatterplot of Q_t2 shows three branches. The most negative flows, which correspond to aquifer discharge into the reservoir, at time t2 are associated with low values of S_S (p < 30%) while the largest positive water flows (flow from the reservoir into the aquifer) occur for high S_S values (p > 70%). The water flow oscillates from positive (aquifer recharge) to negative (aquifer discharge) for intermediate values of S_S (30% < p < 70%). The amplitude of the oscillations of the piezometric heads in the aquifer, A_h, is inversely proportional to the value of the specific storage of the aquifer. When S_S is low, A_h is large and the piezometric heads in the alluvial increase quickly and may rise above the reservoir water level, thus leading to groundwater discharge to the reservoir. On the other hand, when S_S is large, A_h is small and the piezometric heads increase slowly. Here, piezometric heads can hardly rise above the reservoir water level during the events of rise of the reservoir water level, Then, the reservoir recharges the aquifer.

Figure 7 and Figure 8 show CUSUNORO curves for computed heads in wells ST1C, PS19B and SPN1 at times t1, t2 and t3 and MAEg, NRMSEg, NSEg, Q_t1, Q_t2, Q_t3 and q_av.

From the CUSUNORO curves above, parameters K3 (green), Q₂ (red), and K2 (orange) consistently exhibit the highest maximum absolute values across all hydraulic head plots, except for the hydraulic head at well SPN1, where parameters α_r (yellow), K_Vs1 (brown) and Ss (purple) are more influential. Q₆ (blue) and K1 (blue) typically show the lowest maximum absolute values. The most influential parameters across all performance metrics (MAEg, NRMSEg and NSEg) are K3, Q₂, K2 and α_r. Q₆ and Ss exhibit lower values. Parameters α_r, K_Vs1 and Q₂ have the largest influence on the aquifer/reservoir fluxes except at time t2, Q_t2. At time t2, Ss is the most influential parameter. Again, Q6 and K1 consistently show the smallest influences.

Q₂, K3 and K2 have the highest influence on the average groundwater Darcy velocity modulus near well PS16C, while Q₆ and Ss are the least influential parameters for this output.

Most of the CUSUNOTO curves do no show crossings of the x-axis. This attests the expected monotonicity of the outputs versus the input parameters. Some crossings of the x-axis are found in the curves of the least influential parameters, especially for Q₆ and K1 (not shown here).

3.2.2. VARS Results

Figure 9 and Figure 10 show the IVARS₅₀ indexes for the global mean absolute error (MAEg), and the average groundwater Darcy velocity modulus (q_av) as a function of the number of star centers. They also show the ranking of the input parameters and the robustness of ranking. It should be noticed that IVARS₅₀ achieves stable values after just 50 star centers, which amounts to 7700 runs. The largest sensitivity indexes for MAEg correspond to the aquifer hydraulic conductivity in material zone 3 (K3), the boundary inflow Q₂ which corresponds to groundwater flow coming from the Sardas landfill, the vertical conductivity of the silts (K_Vs1) and the aquifer/river conductance (α_r). IVARS₅₀ of the Darcy velocity, q_av, is largest for Q₂, K3 and K2. The rest of the parameters have much smaller indexes.

The robustness of ranking of input parameters K2, K3, S_S, r_u, and Q₇ for both outputs is greater than 90% after 50 stars. However, for the rest of the input parameters, the robustness is smaller than 80 %, and, sometimes, not stable even after 200 stars. Input parameters α_r and K_Vs1 are very influential for MAEg, but they show a robustness measure that ranges from 50 to 70 % after 100 star centers. Robustness does not directly increase with the number of star centers for some input parameters. This could be caused by interactions among S_S, aquitard vertical conductivity K_Vs1 and aquifer/river conductance α_r at some water level rise events. Another reason for the lack of stability of the robustness of rankings is that the least relevant input variables interfere with the calculation of sensitivity indexes of the most relevant parameters. However, the rankings of the most and the least significant input parameters are stable with just 50 star centers (7700 runs). Input parameters with intermediate influence are also well identified, despite not being ranked reliably.

Sample variograms along the directions of the 17 input parameters were computed for the piezometric heads in wells ST1C, PS19B and SPN1 at times t1, t2 and t3, the calibration metrics, the aquifer/reservoir fluxes at times t1, t2 and t3 and the average groundwater Darcy velocity modulus near well PS16C. Only a limited number of parameters are relevant for each output variance. Figure 11 and Figure 12 show the sample variograms along the directions of the 5 most influential parameters. VARS-TO, IVARS₅₀ and VARS-ABE metrics are computed from the variograms. Figure 13 and Figure 14 show the VARS-ABE, IVARS₅₀ and VARS-TO indexes for each output considering all the input parameters.

The most influential parameters for the computed head in well ST1C at times t1, t2 and t3 are K3 and Q₂, followed by K_Vs1, α_r and K2. The ranking in well PS19B is similar to the ranking in well ST1C. However, the sensitivity indexes of K3, K2 and Q₂ in well ST1C are much larger than those of well PS19B. K3 is much more relevant than K2 in well ST1C while K2 and K3 have almost similar sensitivity indexes in well PS19B. The sensitivity indexes of the hydraulic conductivities depend on the location of the wells. Material zone 3 is the largest zone of the Sardas site. Material zone 2 is located between material zones 1 and 3. Material zone 1 is much smaller than the other two and is located just downstream the Sardas landfill. K_Vs1 is more influential for the head in well ST1C than for the well PS19B. ST1C is located right next to the reservoir maximum flood area and PS19B is located 135 m to the east of the reservoir. Tidal effects on the aquifer depend on the duration of the high reservoir level, its amplitude, and the distance to the reservoir. Likewise, the head in well PS19B is most sensitive to Q₂ because this well is near the boundary just downstream the Sardas landfill.

The sensitivity indexes for the computed heads in well SPN1 differ from those of other monitoring wells because well SPN1 is near the Gállego riverbed. α_r is the most influential parameter in well SPN1 followed by K_Vs1.

Despite not being an influential input, the relevance of S_S for the computed heads in the monitoring wells increases from time t1 (low level) to time t2 (peak reservoir water level). The sensitivity indexes of S_S decrease when the reservoir water level descends at time t3. The time change of the sensitivity of S_S in well SPN1 is more significant than in wells ST1C and PS19B.

The most influential parameters for the calibration metrics (MAEg, NRMSEg and NSEg) are K3, K2 and Q₂, followed by K_Vs1 and α_r. Most of the monitoring wells are located in material zones 2 and 3 and near the boundary zone corresponding to Q₂. Aquifer/reservoir and aquifer/river interactions affect the computed head gradient, t_R and A_h. S_S only affects the temporal variability of computed heads and not their average values. Since calibration metrics evaluate the average fit of the computed heads in the monitoring wells, the relevance of S_S for the calibration metrics is very low compared to its relevance on the computed heads at specific times. S_S is especially relevant during events of sudden rise of the reservoir water level.

To facilitate the interpretation of the sensitivities of reservoir/aquifer fluxes (Q_t1, Q_t2, Q_t3), one should recall that the silting sediments (K_Vs1 and K_Vs2) are more permeable than alluvial silts (K_Vs3). Groundwater discharges mainly through the former course of the Gállego river, where only silting sediments confine the aquifer (K_Vs1). If the aquifer is less connected to the Gállego river (lower values of α_r), groundwater discharges to the reservoir. On the other hand, if the aquifer and the river are more connected (higher values of α_r), groundwater discharges to the reservoir and the river. As expected, the greater the inflow of water from the landfill (Q₂), the more groundwater discharge to the reservoir.

The aquifer/reservoir groundwater flow changes with time due to the reservoir tidal effect. The most influential parameters for the aquifer/reservoir groundwater flow are also time dependent. K_Vs1 is the parameter with the largest sensitivity index at time t1 when the reservoir water level is low. The second influential parameter is the aquifer/river leakage coefficient, α_r. The boundary inflow Q₂ is the third most relevant parameter. However, when the reservoir water level rises suddenly rises at time t2, the parameter sensitivities change. The most influential parameter becomes S_S followed by K_Vs1. We recall that in the model reference run [38], there is flow from the reservoir into the aquifer when the reservoir water level rises above the piezometric head in the aquifer. If the specific storage of the aquifer is high, a larger part of the flow coming from the reservoir is stored in the aquifer, which results in a smaller rise of the piezometric head in the aquifer.

The reservoir flood area depends on its water level, so when the water level rises at time t2, the area where the alluvial silts confine the aquifer underneath the reservoir increases. The areas of the reservoir away from the former course of the Gállego river are assumed to be confined by both silting sediments (K_Vs2) and alluvial silts (K_Vs3). The relevance of K_Vs3 increases slightly when the reservoir floods more parts of the alluvial, but its relevance is still much smaller than those of K_Vs1 and S_S. When the water level of the reservoir descends at time t3, the parameter sensitivities tend to be similar to those of Q_t1. For Q_t3, the sensitivity of K_Vs1 is much larger than those of α_r and Q₂. The sensitivity of S_S remains, but it is much less relevant than for Q_t2,

The most influential parameters for the average modulus of the Darcy velocity (q_av) near well PS16C are Q₂, K3 and K2. K1 is slightly relevant, and the contribution of the rest of the parameters is negligible. Darcy velocity in the aquifer mainly depends on the boundary inflow Q₂ because this well is located in material zone 3 and near the boundary downstream the Sardas landfill.

The sensitivity indexes of S_S are time dependent for the computed heads in the wells, and the aquifer/reservoir fluxes (Figure 15). When the reservoir water level rises above the piezometric head in the aquifer at time t2, the reservoir starts recharging the aquifer, and the piezometric head in the alluvial starts rising. If the specific storage is small, the amplitude of the oscillation of the piezometric head increases. This affects especially the aquifer/reservoir fluxes (Q_t1, Q_t2 and Q_t3). The sensitivities of the computed head in well SPN1, located near the Gállego river, are largest for α_r and K_Vs1 and S_S at time t2.

On the other hand, well ST1C is closest to the reservoir and further away from the Gállego river. In this well the sensitivity index of S_S is smaller than in well SPN1. Finally, the sensitivity of the head in well PS19B to S_S is irrelevant because the well is far away from the Gállego river and the reservoir flood area.

The IVARS₅₀ sensitivity indexes for calibration metrics, MAEg, NRMSEg and NSEg are very similar (Figure 16). The most influential parameters are the same for all the metrics. The rankings show some slight differences. Input ranking for the MAEg is K3, Q₂, α_r, K2 and K_Vs1, while α_r and K2 switch places for NRMSEg and NSEg. NRMSEg and NSEg are more prone to the presence of outliers because their formulas include squared residuals.

Sensitivity indexes are large for K3, Q₂, α_r, K2 and K_Vs1. Some input variables have very small sensitivity indexes.

3.2.3. HDMR Results and Analysis of Interactions for the Sobol Sequence

Tables S1 to S23 in the SM show the values of the first order (main effects), S_i, and 2^nd order effects, S_ij, calculated with SALib [63] and GUI-HDMR [64] and the parameter ranking for: (1) The computed piezometric heads in wells PS19B, SPN1 and ST1C at times t1, t2 and t3; (2) The global mean absolute error (MAEg); (3) The global Nash–Sutcliffe index (NSEg); (4) The global normalized root mean squared error (NRMSEg); (5) The average groundwater Darcy velocity modulus near well PS16C (q_av); and (6) The computed aquifer/reservoir fluxes at times t1, t2 and t3 (Q_t1, Q _t2 and Q_t3) by using 16 384 runs.

HDMR lower bound estimates for the total effects, (S_i + S_ij), are compared to the VARS intervals of the total effects for each parameter. SALib and GUI-HDMR indexes generally agree, although they show some small differences. The analysis of the interactions among parameters is an important part of the global sensitivity analysis. Interactions represent the joint influence of parameters on the model outputs. Usually, interactions are revealed when the sum of the Sobol’s 1^st order indexes are significantly smaller than 1.

In the next paragraphs the Sobol based lower bound estimates for the total effects will be denoted simply as “total effects” to shorten the presentation of the HDMR results. The largest main effect S_i for the piezometric head in well ST1C at time t1 is equal to 0.242 and corresponds to K3 (see Table S1 in SM). The sum of the Sobol’s 1^st order indexes is equal to 0.644. The total effects of the piezometric head in well ST1C at time t1 are slightly larger than 1 (see Table S1 in SM). Second order effects for this output are important especially due to the interaction between K3 and K2. Interactions of smaller relevance occur between Q₂ and K3, between K_Vs1 and α_r, and between Q₂ and K2. Sobol total effects S_Ti are typically larger and fall outside the ranges of the VARS total effects.

The largest main effect S_i of the piezometric head in well ST1C at time t2 is equal to 0.194 and corresponds to K3 (see Table S2 in SM). The sum of the Sobol’s 1^st order indexes is equal to 0.628. The total effects of the piezometric head in well ST1C at time t2 are slightly larger than 1 (see Table S2 in SM). Second order effects here are relevant especially due to the interaction between Q₂ and K3. There are also interactions of smaller relevance between K3 and K2, between K_Vs1 and α_r, and between Q₂ and K2. Sobol total effects S_Ti and those of VARS generally agree, although they show some discrepancies especially for the K1 and r_u.

The largest main effects S_i for the piezometric head in well ST1C at time t3 correspond to K3 and Q₂ (see Table S3 in SM). The sum of the Sobol’s 1^st order indexes is equal to 0.644. The total effects of the piezometric head in well ST1C at time t3 are equal to 1.223 (see Table S3 in SM). Second order effects for this output are relevant due to interactions of K3 with K2 and Q₂, and K_Vs1 with α_r. The Sobol total effects S_Ti generally fall within the intervals of VARS total effects.

The largest main effects S_i for the piezometric head in well PS19B at time t1 correspond to Q₂, K3 and K2 (see Table S4 in SM). The sum of the Sobol’s 1^st order indexes is equal to 0.71. The total effects S_Ti for the piezometric head in well PS19B at time t1 are slightly greater than 1 (see Table S4 in SM). Second order effects here are relevant especially due to the interactions among Q₂, K3 and K2. The total effects S_Ti fall within the intervals of VARS total effects.

The largest main effects S_i for the piezometric head in well PS19B at time t2 correspond to Q₂, K2, and K3 (see Table S5 in SM). The sum of the Sobol’s 1^st order indexes is equal to 0.699. The total effects S_Ti for the piezometric head in well PS19B at time t2 are slightly greater than 1 (see Table S5 in SM). Second order effects for this output are important especially due to the interaction between Q₂ and K2. There are also interactions of smaller relevance between Q₂ and K3, between K3 and K2, and between α_r and K_Vs1. The Sobol total effects S_Ti are typically larger and fall outside the ranges of the VARS total effects.

The largest main effect S_i for the piezometric head in well PS19B at time t3 corresponds to Q₂. K2 and K3 are in second and third position (see Table S6 in SM). The sum of the Sobol’s 1^st order indexes is equal to 0.703. Second order effects for this output are important especially due to the interaction between Q₂ and K2. There are also interactions of smaller relevance between Q₂ and K3, between K3 and K2, and between α_r and K_Vs1. The Sobol total effects S_Ti are generally larger and out of the intervals of the VARS estimates.

α_r shows the largest main effect S_i for the piezometric head in well SPN1 at time t1 which is equal to 0.394. The second largest effect corresponds to K_Vs1 (see Table S7 in SM). The sum of the Sobol’s 1^st order indexes is equal to 0.69. The total effects of the piezometric head in well SPN1 at time t1 are close to 1 (see Table S7 in SM). Second order effects for this output are relevant especially due to the interaction between Q₂ and K3. The Sobol total effects S_Ti are overall higher and beyond the limits of the VARS estimates.

The largest main effects S_i for the piezometric head in well SPN1 at time t2 correspond to α_r and S_S (see Table S8 in SM). The sum of the Sobol’s 1^st order indexes is equal to 0.655. The total effects of the piezometric head in well SPN1 at time t2 are slightly larger than 1 (see Table S8 in SM). Second order effects for this output are relevant especially due to the interaction between α_r and K_Vs1. There are also smaller interactions of Q₂ with K3 and with α_r. The Sobol total effects are generally larger and out of the intervals of the VARS estimates.

The largest main effects S_i for the piezometric head in well SPN1 at time t3 correspond to α_r and K3 (see Table S9 in SM). The sum of the Sobol 1^st order indexes is equal to 0.679. The total effects of the piezometric head in well SPN1 t3 are close to 1 (see Table S9 in SM). Second order effects are relevant mainly due to the interaction between α_r and K_Vs1, and to a lower extent to the interactions of Q₂ with K3 and α_r. The Sobol total effects are generally larger and out of the intervals of the VARS estimates.

The largest main effects S_i for the global mean absolute error correspond to K3, Q₂ and α_r (see Tables S10 and S11 in SM). The sums of Sobol’s 1^st order indexes calculated with SALib and GUI-HDMR are equal to 0.682 and 0.734, respectively. The total effect of the global mean absolute error computed with SALib is equal to that of GUI-HDMR (see Tables S10 and S11 in SM). Second order effects for this output are important especially due to the interactions between Q₂ and K3 and between K3 and K2. There are also interactions of smaller relevance between α_r and K_Vs1, and between Q₂ and K2. The SALib total effects S_Ti fall within the intervals of VARS estimates or are slightly larger than VARS total effects.

K3 shows the largest main effect S_i for the global normalized root mean squared error which is equal to 0.245 (SALib). The second largest effect corresponds to Q₂ (see Tables S12 and S13 in SM). The sums of the Sobol’s 1^st order indexes calculated with SALib and GUI-HDMR are equal to 0.675 and 0.731, respectively. The total effects of the global normalized root mean squared error are close to 1 (see Tables S12 and S13 in SM). Second order effects for this output are relevant especially due to the interaction between K3 and K2, and between Q₂ and K3. The total effects S_Ti from SALib are either within the ranges of VARS estimates or slightly exceed the VARS total effects.

The largest main effect S_i for the global Nash–Sutcliffe index corresponds to K3. Q₂ and K2 are in second and third position (see Tables S14 and S15 in SM). The sums of the Sobol’s 1^st order indexes calculated with SALib and GUI-HDMR are equal to 0.411 and 0.477, respectively. The total effect of the global Nash–Sutcliffe index computed with SALib is smaller than that of GUI-HDMR (see Tables S14 and S15 in SM). Second order effects for this output are important especially due to the interaction between Q₂ and K3. There are also interactions of smaller relevance between K3 and K2, between Q₂ and K2, and between α_r and K_Vs1. The SALib total effects S_Ti fall within the intervals of VARS estimates or are slightly larger than VARS total effects.

The largest main effects S_i for the computed aquifer/reservoir flux at time t1 correspond to K_Vs1 and α_r (see Tables S16 and S17 in SM). The sums of the Sobol’s 1^st order indexes calculated with SALib and GUI-HDMR are equal to 0.716 and 0.746, respectively. The total effects of the computed aquifer/reservoir flux at time t1 are close to 1 (see Tables S16 and S17 in SM). The second order effects for this output are significant, particularly because of the interactions of K_Vs1 with α_r, Ss, Q₂ and K3. Additionally, there are less significant interactions between α_r and K3, and between α_d and K4. The Sobol total effects S_Ti calculated by SALib are overall higher than the VARS estimates.

The largest main effect S_i for the computed aquifer/reservoir flux at time t2 is equal to 0.362 (SALib) and corresponds to S_S (see Table S18 in SM). The sums of the Sobol’s 1st order indexes calculated with SALib and GUI-HDMR are equal to 0.514 and 0.469, respectively. The total effects of the computed aquifer/reservoir flux at time t2 are greater than 1 (see Tables S18 and S19 in SM). For this output, the second order effects are notable, mainly due to the interaction between S_S and K_Vs1. Sobol total effects S_Ti and those of VARS estimates show clear discrepancies, especially for K1, K2 and Q₇.

The largest main effect S_i of the computed aquifer/reservoir flux at time t3 is equal to 0.626 (SALib) and corresponds to K_Vs1 (see Table S20 in SM). The sums of the Sobol’s 1st order indexes calculated with SALib and GUI-HDMR are 0.812 and 0.819, respectively. The total effects of the computed aquifer/reservoir flux at time t3 are close to 1 (see Tables S20 and S21 in SM). Second order effects for this output are relevant to the interactions of K_Vs1 with α_r, Ss, Q₂, K_Vs3, K4 and K3. The Sobol total effects S_Ti are generally larger and out of the intervals of the VARS estimates.

The largest main effects S_i for the average modulus of the Darcy velocity near well PS16C correspond to Q₂ and K3 (see Tables S22 and S23 in SM). The sums of the Sobol’s 1^st order indexes calculated with SALib and GUI-HDMR are equal to 0.855 and 0.861, respectively. The total effects of the average modulus of Darcy velocity are slightly larger than 1 (see Tables S22 and S23 in SM). Second order effects for this output are relevant to the interaction between Q₂ and K3. There are also smaller interactions between Q₂ and K2, and between K3 and K2. The Sobol total effects S_Ti are generally slightly larger and out of the intervals of the VARS total effects.

The 1^st and 2^nd order interactions among input parameters can also be presented visually through heatmaps. Figure S1 in SM shows the heatmaps of the 1^st and 2^nd order HDMR Sobol effects for all the outputs. The analysis of the heatmaps confirms that: (1) The largest components correspond to the main effects (which have been conveniently located along the diagonal of each of the maps); and (2) Off-diagonal terms corresponding to the 2^nd order effects are important for all input variables, especially the global Nash–Sutcliffe index and groundwater flow between aquifer and reservoir at time t2 due to the interactions between Q₂ and K3 and between Ss and K_Vs1.

3.2.4. HDMR Results for the VARS Runs by Using the Halton Sequence

HDMR analyses can be performed reusing the simulations performed with the VARS-Halton sequence. However, preliminary analysis of the main effects, S_i and 2^nd order effects, S_ij, calculated with SALib show that most outputs have Sobol total effects S_Ti that are consistently greater than 2, with some of them even reaching values of 3. Since the parameters are normalized, the sum of main and interaction effects of any order should equal the total variance, which should theoretically be equal to 1. It should be noticed that the VARS procedure to locate the star centers is either the Halton or the Sobol sequences, both of which generate low-discrepancy sequences [65]. They are intended to optimize quasi-Monte Carlo approaches by maximizing the efficiency to fill the parameter hypercube with sparsely located points. As the number of parameters increases, the capacity of the Halton sequence [54] to distribute points uniformly rapidly decreases [65], while Sobol sequence is more stable. This results in higher maximum errors of integration for high dimensional data related to the Halton sequence. The Halton sequence gives good quality, near uniform distributions when the number of parameters is lower than 10 [65], as it was demonstrated recently in a reactive transport model for a high-level radioactive repository in granitic rock [36].

3.3. Input Parameter Rankings

Table 2 shows the rankings of the parameters for the computed piezometric heads in well ST1C at time t1 (ST1C_t1), global mean absolute error (MAEg), aquifer/reservoir flux at time t1 (Q_t1), and average groundwater Darcy velocity modulus (q_av) across the following methods: IVARS₅₀, VARS-TO (equivalent to Sobol), VARS-ABE (equivalent to Morris), HDMR (SALib), GUI-HDMR and CUSUNORO plots.

All the methods agree in identifying the most influential parameters for ST1C_t1, although they show different positions for some input parameters. The most influential parameters for ST1C_t1 are K3 and Q₂. Q₂ is the most influential parameter, followed by K3.

Output MAEg shows similar results with the following order of ranking: (1) K3; (2) Q2; (3) α_r (except for VARS-TO); (4) K2 (except for VARS-TO and SALib); and (5) K_Vs1 (except for SALib). GSA methods also agree on the ranking of the least relevant input parameters but switching in some positions.

The three most influential parameters for Q_t1 are K_Vs1, α_r and Q₂. There are discrepancies among the methods on the ranking of the 4^th and 5^th positions (K3 and S_S). SALib reveals the highest discrepancy for Q_t1 compared to the other methods for the parameters ranging from inputs relevance from 4^th to 10^th rank.

All methods provide a similar ranking of the most influential parameters for q_av. Q₂ is the most influential parameter followed by K3. The third most influential is K2 and K1 is the fourth. α_r is at the 5^th position and K_Vs1 is the 6^th most influential parameter, except for the VARS-TO method which interchanges the positions of these two parameters.

Table S24 to Table S29 present the rankings of the parameters for all the outputs. For the most part, all methods agree in identifying the most and the least relevant inputs for all outputs. Some methods switch the order within the most and within the least influential inputs. On the other hand, the inputs located at intermediate positions (from 5^th to 13^th place) show less consistency across the different GSA methods. Their ranking does not usually change more than three places, but some outputs show important differences (K4 for Q_t3 ranges from 7^th to 15^th position; K_Vs2 for SPN1_t3 ranges from 7^th to 13^th position).

The rankings of the input parameters derived from CUSUNORO curves for computed heads in wells ST1C, PS19B and SPN1 at times t1, t2 and t3 agree with the rankings of VARS and HDMR methods. Unlike the Morris method, CUSUNORO plots are well suited to identify the most influential input parameters.

4. Conclusions

VARS and HDMR GSA of the groundwater flow model of the Gállego alluvial aquifer has been presented. Computed piezometric heads in monitoring wells and aquifer/reservoir fluxes change with time due to the tidal effect caused by the daily oscillations of the water level in the Sabiñánigo reservoir. Therefore, the sensitivities of the heads and fluxes change wit time. The results of the GSA lead to following conclusions:

• The most influential parameters for the selected outputs are consistently detected by all methods. They include: K2, K3, K_Vs1, S_S, Q₂ and α_r.
• While some parameter inputs such as K3 and Q₂ are relevant for all the outputs, other parameter inputs such as K1 and S_S are influential only for some outputs.
• The sensitivity indexes of the computed heads in monitoring wells and aquifer/reservoir fluxes with respect to S_S change with time.
• Sensitivity indexes of the calibration metrics are similar. MAEg is less prone to model result outliers.
• The average groundwater Darcy velocity near well PS16C depends mainly on the boundary inflow Q₂.
• VARS achieves stable values for the most important and the least influential input parameters after 50 star centers, which amounts to 7700 runs. For other inputs, the robustness of the ranking does not increase monotonically with the number of star centers.
• VARS and HDMR methods provide similar results in terms of rankings and significance of the most influential parameters. However, they show slight differences in the ranking of parameters of intermediate and low influence. The ranking of the least relevant input variables with the different methods is less consistent.
• Graphical methods and HDMR results highlight that the most important input parameter interactions occur between S_S and K_Vs1 for groundwater flow between aquifer/reservoir groundwater flux when the water level of the reservoir is high at time t2.

Future work should be devoted to extending GSA to other model outputs such as aquifer/river fluxes, discharges underneath the dam foundation or aquifer/reservoir fluxes in other parts of the reservoir. Sensitivity indexes of the heads and aquifer/reservoir fluxes are time dependent due to the tidal effect of the Sabiñánigo reservoir. It might be advisable to include more time-dependent outputs to capture the effect of the oscillations of the water level of the reservoir.

The influence of extreme results on some calibration metrics such as NRMSEg and NSEg could be overcome by using calibration metrics immune to the outliers, such as the median absolute deviation. Furthermore, ranges in some inputs could be revised as more data becomes available. There are no reliable data on the boundary inflow, Q₂, from Sardas landfill. Leachate estimations based on groundwater flow models range from 17 to 32 m³/d. This range is much smaller than the range used in the GSA presented here. Moreover, recent monitoring wells closer to the Sabiñánigo reservoir reveal that the hydraulic conductivity of the aquitard (silting sediments and alluvial silts) is very heterogeneous.

On the other hand, only a limited number of parameters are relevant for each output variance, and some inputs’ contribution to the variance of the results is negligible. It might be advisable to reduce the number of parameters and analyze further the results of the HDMR analyses for the Halton sequence.

Parameter ranking is useful to identify the most and the least influential input parameters. However, parameter ranking only provides information on the ordering, not on the values of the sensitivity indexes. The results presented in the preceding sections suggest that most outputs are mostly sensitive to 5 input parameters. The impact on model outputs of the uncertainty of the least relevant parameters is almost irrelevant.

GSA methods presented here provide a quantitative tool to assess the impact of the uncertainty of parameters on the groundwater flow model outputs in alluvial aquifers. The findings of this study can guide future management and data acquisition in polluted sites to reduce uncertainties related to the most relevant parameters. Moreover, these methods can guide future uncertainty analysis of the total dissolved hexachlorocyclohexane transport model through the Gállego alluvial aquifer presented in Sobral et al. [38].