Prerpints.org logo
Preprint
Article

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

Altmetrics

Downloads

168

Views

62

Comments

0

A peer-reviewed article of this preprint also exists.

supplementary.zip (709.45KB )

Submitted:

26 August 2024

Posted:

26 August 2024

You are already at the latest version

Alerts
Abstract
Groundwater flow and transport models are essential tools for assessing and quantifying the mi-gration of organic contaminants at polluted sites. Uncertainties in the hydrodynamic and transport parameters of the aquifer have a significant effect on model predictions. Uncertainties can be quantified with advanced sensitivity methods such as Sobol’s High Dimensional Model Reduction (HDMR) and Variogram Analysis of Response Surfaces (VARS). Here we present the application of VARS and HDMR to assess the global sensitivities of the outputs of a transient groundwater flow model of the Gállego alluvial aquifer which is located downstream the Sardas landfill in Huesca (Spain). The aquifer is subject to the tidal effects caused by the daily oscillations of the water level in the Sabiñánigo reservoir. Global sensitivities are analysed for hydraulic heads, aquifer/reservoir fluxes, groundwater Darcy velocity and hydraulic head calibration metrics. Input parameters in-clude aquifer hydraulic conductivities and specific storage, aquitard vertical hydraulic conductiv-ities, and boundary inflows and conductances. VARS, HDMR and graphical methods agree to identify the most influential parameters which for most of the outputs are the hydraulic conduc-tivities of the zones closest to the landfill, the vertical hydraulic conductivity of the most permeable zones of the aquitard, and the boundary inflow coming from the landfill. The sensitivity of heads and aquifer/reservoir fluxes with respect to specific storage change with time. The aquifer/reservoir flux when the reservoir level is high shows interactions between specific storage and aquitard conductivity. VARS and HDMR parameter rankings are similar for the most influential parameters. However, there are discrepancies for the less relevant parameters. The efficiency of VARS was demonstrated by achieving stable results with a relatively small number of simulations.
Keywords: 
Subject: Environmental and Earth Sciences  -   Water Science and Technology

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). Ns simulations of groundwater flow were performed with CORE2D for the selected combinations of input parameters. The outputs of the flow model were postprocessed to generate the tables containing the Ni input parameters and No 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·104 and 8·104 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 m3 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 i = 1 n · S y y · ∑ j = 1 i   ( y ~ π j − y ¯ )
where y ¯ and S y y are the sample mean and standard deviation of the outputs, respectively, n is sample size or the number of simulations, and   ( y ~ π j − 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, n , which guarantees that the cumulative sum is scaled relative to the sample size n; (2) A divisor term S y y (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 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, Si, measures the relative individual contribution of the i-th input variable on the variance of the output. The second order sensitivity indexes, Sij, measure the interactions between the i-th and the j-th inputs on the output variable [22]. Both terms, Si and Sij, are usually provided by HDMR software packages. The total sensitivity index, or total-effect index or total-order index, STi, 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 T i = S i + ∑ j S i j + ∑ j , k S i j k + …   ∑ j , … , n S i j … . . n  

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 h2 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 (IVARS50). IVARS50 provides a sensitivity index corresponding to intermediate scales and is calculated from the directional variograms according to:
I V A R S 50 = ∫ 0 0.5 γ h   d h
where γ h 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 (SS); (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 (KVs1), the silting sediments in the rest of the reservoir (KVs2) and the alluvial silts (KVs3); (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 (rc) and unconfined (ru) parts of the aquifer; and (6) Boundary inflows (Q6, Q7, Q9, Q2 and Q1). 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 ST1Ct1, ST1Ct2, ST1Ct3, PS19Bt1, PS19Bt2, PS19Bt3, SPN1t1, SPN1t2 and SPN1t3. Other outputs include the aquifer/reservoir fluxes at times t1, t2 and t3, which are denoted as Qt1, Qt2 and Qt3, and the average modulus of the Darcy velocity near well PS16C where a tracer test was performed (qav).
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:
M A E g = ∑ j M A E j N w  
N R M S E g = 1 N T ∑ j N j   R M S E j 2   σ T  
N S E g = 1 + ∑ j N j σ j 2 N S E j − 1 N T · σ T 2  
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.
  • σ j is the standard deviation of the measured piezometric heads in the j-th well.
  • M A E j , R M S E j and N S E j are the mean absolute error, the root mean squared error and the Nash–Sutcliffe index for the j-the well.
  • σ 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 CORE2D 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. CORE2D V5 has undergone extensive verification against analytical solutions and has been benchmarked against other reactive transport codes [57,58]. CORE2DV5 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 CORE2D 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, Ah and the time lag, tR, are very sensitive to the specific storage coefficient of the layer of sands and gravels (SS), and the vertical hydraulic conductivities of the aquitard (KVs1, KVs2 and KVs3) [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 ST1Ct2, PS19Bt2 and SPN1t2) and the aquifer/reservoir flow (Qt2) at time t2 corresponding the peak reservoir water level. The outputs are plotted versus the vertical conductivity of the aquitard (KVs1). 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 (SS): 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 SS. Generally, the computed piezometric heads in wells ST1C and SPN1 at t2 are low for large SS (p > 70%). However, simulations corresponding to extreme values of other input parameters result in high piezometric heads even for high values of SS. On the other hand, the computed heads are higher when SS 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 KVs1 and SS.
The scatterplot of Qt2 shows three branches. The most negative flows, which correspond to aquifer discharge into the reservoir, at time t2 are associated with low values of SS (p < 30%) while the largest positive water flows (flow from the reservoir into the aquifer) occur for high SS values (p > 70%). The water flow oscillates from positive (aquifer recharge) to negative (aquifer discharge) for intermediate values of SS (30% < p < 70%). The amplitude of the oscillations of the piezometric heads in the aquifer, Ah, is inversely proportional to the value of the specific storage of the aquifer. When SS is low, Ah 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 SS is large, Ah 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, Qt1, Qt2, Qt3 and qav.
From the CUSUNORO curves above, parameters K3 (green), Q2 (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), KVs1 (brown) and Ss (purple) are more influential. Q6 (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, Q2, K2 and αr. Q6 and Ss exhibit lower values. Parameters αr, KVs1 and Q2 have the largest influence on the aquifer/reservoir fluxes except at time t2, Qt2. At time t2, Ss is the most influential parameter. Again, Q6 and K1 consistently show the smallest influences.
Q2, K3 and K2 have the highest influence on the average groundwater Darcy velocity modulus near well PS16C, while Q6 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 Q6 and K1 (not shown here).

3.2.2. VARS Results

Figure 9 and Figure 10 show the IVARS50 indexes for the global mean absolute error (MAEg), and the average groundwater Darcy velocity modulus (qav) 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 IVARS50 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 Q2 which corresponds to groundwater flow coming from the Sardas landfill, the vertical conductivity of the silts (KVs1) and the aquifer/river conductance (αr). IVARS50 of the Darcy velocity, qav, is largest for Q2, K3 and K2. The rest of the parameters have much smaller indexes.
The robustness of ranking of input parameters K2, K3, SS, ru, and Q7 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 KVs1 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 SS, aquitard vertical conductivity KVs1 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, IVARS50 and VARS-ABE metrics are computed from the variograms. Figure 13 and Figure 14 show the VARS-ABE, IVARS50 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 Q2, followed by KVs1, αr and K2. The ranking in well PS19B is similar to the ranking in well ST1C. However, the sensitivity indexes of K3, K2 and Q2 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. KVs1 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 Q2 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 KVs1.
Despite not being an influential input, the relevance of SS 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 SS decrease when the reservoir water level descends at time t3. The time change of the sensitivity of SS 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 Q2, followed by KVs1 and αr. Most of the monitoring wells are located in material zones 2 and 3 and near the boundary zone corresponding to Q2. Aquifer/reservoir and aquifer/river interactions affect the computed head gradient, tR and Ah. SS 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 SS for the calibration metrics is very low compared to its relevance on the computed heads at specific times. SS is especially relevant during events of sudden rise of the reservoir water level.
To facilitate the interpretation of the sensitivities of reservoir/aquifer fluxes (Qt1, Qt2, Qt3), one should recall that the silting sediments (KVs1 and KVs2) are more permeable than alluvial silts (KVs3). Groundwater discharges mainly through the former course of the Gállego river, where only silting sediments confine the aquifer (KVs1). 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 (Q2), 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. KVs1 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 Q2 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 SS followed by KVs1. 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 (KVs2) and alluvial silts (KVs3). The relevance of KVs3 increases slightly when the reservoir floods more parts of the alluvial, but its relevance is still much smaller than those of KVs1 and SS. When the water level of the reservoir descends at time t3, the parameter sensitivities tend to be similar to those of Qt1. For Qt3, the sensitivity of KVs1 is much larger than those of αr and Q2. The sensitivity of SS remains, but it is much less relevant than for Qt2,
The most influential parameters for the average modulus of the Darcy velocity (qav) near well PS16C are Q2, 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 Q2 because this well is located in material zone 3 and near the boundary downstream the Sardas landfill.
The sensitivity indexes of SS 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 (Qt1, Qt2 and Qt3). The sensitivities of the computed head in well SPN1, located near the Gállego river, are largest for αr and KVs1 and SS 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 SS is smaller than in well SPN1. Finally, the sensitivity of the head in well PS19B to SS is irrelevant because the well is far away from the Gállego river and the reservoir flood area.
The IVARS50 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, Q2, αr, K2 and KVs1, 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, Q2, αr, K2 and KVs1. 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), Si, and 2nd order effects, Sij, 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 (qav); and (6) The computed aquifer/reservoir fluxes at times t1, t2 and t3 (Qt1, Q t2 and Qt3) by using 16 384 runs.
HDMR lower bound estimates for the total effects, (Si + Sij), 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 1st 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 Si 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 1st 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 Q2 and K3, between KVs1 and αr, and between Q2 and K2. Sobol total effects STi are typically larger and fall outside the ranges of the VARS total effects.
The largest main effect Si 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 1st 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 Q2 and K3. There are also interactions of smaller relevance between K3 and K2, between KVs1 and αr, and between Q2 and K2. Sobol total effects STi and those of VARS generally agree, although they show some discrepancies especially for the K1 and ru.
The largest main effects Si for the piezometric head in well ST1C at time t3 correspond to K3 and Q2 (see Table S3 in SM). The sum of the Sobol’s 1st 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 Q2, and KVs1 with αr. The Sobol total effects STi generally fall within the intervals of VARS total effects.
The largest main effects Si for the piezometric head in well PS19B at time t1 correspond to Q2, K3 and K2 (see Table S4 in SM). The sum of the Sobol’s 1st order indexes is equal to 0.71. The total effects STi 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 Q2, K3 and K2. The total effects STi fall within the intervals of VARS total effects.
The largest main effects Si for the piezometric head in well PS19B at time t2 correspond to Q2, K2, and K3 (see Table S5 in SM). The sum of the Sobol’s 1st order indexes is equal to 0.699. The total effects STi 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 Q2 and K2. There are also interactions of smaller relevance between Q2 and K3, between K3 and K2, and between αr and KVs1. The Sobol total effects STi are typically larger and fall outside the ranges of the VARS total effects.
The largest main effect Si for the piezometric head in well PS19B at time t3 corresponds to Q2. K2 and K3 are in second and third position (see Table S6 in SM). The sum of the Sobol’s 1st order indexes is equal to 0.703. Second order effects for this output are important especially due to the interaction between Q2 and K2. There are also interactions of smaller relevance between Q2 and K3, between K3 and K2, and between αr and KVs1. The Sobol total effects STi are generally larger and out of the intervals of the VARS estimates.
αr shows the largest main effect Si for the piezometric head in well SPN1 at time t1 which is equal to 0.394. The second largest effect corresponds to KVs1 (see Table S7 in SM). The sum of the Sobol’s 1st 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 Q2 and K3. The Sobol total effects STi are overall higher and beyond the limits of the VARS estimates.
The largest main effects Si for the piezometric head in well SPN1 at time t2 correspond to αr and SS (see Table S8 in SM). The sum of the Sobol’s 1st 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 KVs1. There are also smaller interactions of Q2 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 Si 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 1st 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 KVs1, and to a lower extent to the interactions of Q2 with K3 and αr. The Sobol total effects are generally larger and out of the intervals of the VARS estimates.
The largest main effects Si for the global mean absolute error correspond to K3, Q2 and αr (see Tables S10 and S11 in SM). The sums of Sobol’s 1st 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 Q2 and K3 and between K3 and K2. There are also interactions of smaller relevance between αr and KVs1, and between Q2 and K2. The SALib total effects STi fall within the intervals of VARS estimates or are slightly larger than VARS total effects.
K3 shows the largest main effect Si for the global normalized root mean squared error which is equal to 0.245 (SALib). The second largest effect corresponds to Q2 (see Tables S12 and S13 in SM). The sums of the Sobol’s 1st 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 Q2 and K3. The total effects STi from SALib are either within the ranges of VARS estimates or slightly exceed the VARS total effects.
The largest main effect Si for the global Nash–Sutcliffe index corresponds to K3. Q2 and K2 are in second and third position (see Tables S14 and S15 in SM). The sums of the Sobol’s 1st 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 Q2 and K3. There are also interactions of smaller relevance between K3 and K2, between Q2 and K2, and between αr and KVs1. The SALib total effects STi fall within the intervals of VARS estimates or are slightly larger than VARS total effects.
The largest main effects Si for the computed aquifer/reservoir flux at time t1 correspond to KVs1 and αr (see Tables S16 and S17 in SM). The sums of the Sobol’s 1st 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 KVs1 with αr, Ss, Q2 and K3. Additionally, there are less significant interactions between αr and K3, and between αd and K4. The Sobol total effects STi calculated by SALib are overall higher than the VARS estimates.
The largest main effect Si for the computed aquifer/reservoir flux at time t2 is equal to 0.362 (SALib) and corresponds to SS (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 SS and KVs1. Sobol total effects STi and those of VARS estimates show clear discrepancies, especially for K1, K2 and Q7.
The largest main effect Si of the computed aquifer/reservoir flux at time t3 is equal to 0.626 (SALib) and corresponds to KVs1 (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 KVs1 with αr, Ss, Q2, KVs3, K4 and K3. The Sobol total effects STi are generally larger and out of the intervals of the VARS estimates.
The largest main effects Si for the average modulus of the Darcy velocity near well PS16C correspond to Q2 and K3 (see Tables S22 and S23 in SM). The sums of the Sobol’s 1st 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 Q2 and K3. There are also smaller interactions between Q2 and K2, and between K3 and K2. The Sobol total effects STi are generally slightly larger and out of the intervals of the VARS total effects.
The 1st and 2nd order interactions among input parameters can also be presented visually through heatmaps. Figure S1 in SM shows the heatmaps of the 1st and 2nd 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 2nd 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 Q2 and K3 and between Ss and KVs1.

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, Si and 2nd order effects, Sij, calculated with SALib show that most outputs have Sobol total effects STi 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 (ST1Ct1), global mean absolute error (MAEg), aquifer/reservoir flux at time t1 (Qt1), and average groundwater Darcy velocity modulus (qav) across the following methods: IVARS50, 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 ST1Ct1, although they show different positions for some input parameters. The most influential parameters for ST1Ct1 are K3 and Q2. Q2 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) KVs1 (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 Qt1 are KVs1, αr and Q2. There are discrepancies among the methods on the ranking of the 4th and 5th positions (K3 and SS). SALib reveals the highest discrepancy for Qt1 compared to the other methods for the parameters ranging from inputs relevance from 4th to 10th rank.
All methods provide a similar ranking of the most influential parameters for qav. Q2 is the most influential parameter followed by K3. The third most influential is K2 and K1 is the fourth. αr is at the 5th position and KVs1 is the 6th 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 5th to 13th 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 Qt3 ranges from 7th to 15th position; KVs2 for SPN1t3 ranges from 7th to 13th 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, KVs1, SS, Q2 and αr.
  • While some parameter inputs such as K3 and Q2 are relevant for all the outputs, other parameter inputs such as K1 and SS are influential only for some outputs.
  • The sensitivity indexes of the computed heads in monitoring wells and aquifer/reservoir fluxes with respect to SS 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 Q2.
  • 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 SS and KVs1 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, Q2, from Sardas landfill. Leachate estimations based on groundwater flow models range from 17 to 32 m3/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].

Supplementary Materials

The following supporting information can be downloaded at the website of this paper posted on Preprints.org. Figure S1: Heatmaps of the HDMR Sobol 2nd order effects for: (1) Hydraulic heads in wells PS19B, SPN1 and ST1C at time intervals t1, t2 and t3; (2) MAEg; (3) NRMSEg; (4) NSEg; (5) qav; (6) Qt1; (7) Qt2 and (8) Qt3. The map visualizes both main and second order effects; Tables S1 to S23: 1st order (main effects), Si, and 2nd order effects, Sij, calculated with SALib and gui-HDMR of each output for the Sobol sequence. The largest 2nd order effects Sij, are listed in the off-diagonal boxes. Parameter rankings are indicated within brackets. Total effects (Si + Sij) are compared to the VARS interval total effects for each parameter); Tables S24 to S29: Ranking of the influence of parameters for each output.

Author Contributions

Conceptualization, J.S., B.S., and C.L.-V.; Data curation, B.S., and B.P; Funding acquisition, J.S.; Investigation, J.S., B.S., and A.M.; Methodology, J.S., B.S., B.P., A.M., C.L.-V. and J.S.-P.; Project administration, J.S.; Resources, J.S., B.S., and A.M.; Software, B.S., B.P., A.M., C.L.-V. and J.S.-P.; Validation, B.S., A.M., and C.L.-V.; Visualization, B.S., A.M., C.L.-V. and J.S.-P.; Writing—original draft, B.S., and C.L.-V.; Writing—review and editing, J.S., B.P., A.M., and J.S.-P. All authors have read and agreed to the published version of the manuscript.

Funding

Funding for this study was provided by EMGRISA, which was awarded a contract by the Aragon Regional Government for the Hydrogeology of the Sardas landfill site. Funding was also provided by the Ebro Water District. The work of the second author was funded by a research contract from the Spanish Government (Ref. FPU20/04135). Partial funding was also awarded from the Galician Regional Government (Grant ED431C2021/54).

Data Availability Statement

Data from the Sadas site and other lindane-affected areas in Sabiñánigo are available upon request from the Aragón Regional Government at suelos@aragon.es.

Acknowledgments

This study was performed within the framework of research contracts signed by the University of A Coruña Foundation with EMGRISA and the Ebro River Authority. We acknowledge the support received from the Aragon Government, EMGRISA, the Ebro Water District, the Spanish Ministry of Science and Innovation (PID2023-153202OB-I00) and the Galician Regional Government (Grant ED431C2021/54). The work of the second author was funded by research contracts from Xunta de Galicia (Ref. ED481A) and the Spanish Government (Ref. FPU20/04135). The work of the fourth author was funded by the Spanish Ministry of Science and Innovation (TED2021-130315B-I00). Computer runs were performed in the Galician Supercomputing Center (CESGA). We thank CESGA's personnel for their continuous support to our work. We are very thankful to Elmar Plischke from Institute for Resource Ecology, Helmholtz-Zentrum Dresden-Rossendorf (HZDR) for his support on the HDMR analysis of some outputs for which GUI-HDMR failed to provide estimates of Sobol indexes.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Cvetkovic, V.; Soltani, S.; Vigouroux, G. Global Sensitivity Analysis of Groundwater Transport. Journal of Hydrology 2015, 531, 142–148. [Google Scholar] [CrossRef]
  2. Bordbar, M.; Neshat, A.; Javadi, S. A New Hybrid Framework for Optimization and Modification of Groundwater Vulnerability in Coastal Aquifer. Environ Sci Pollut Res 2019, 26, 21808–21827. [Google Scholar] [CrossRef] [PubMed]
  3. Balakrishnan, S.; Roy, A.; Ierapetritou, M.G.; Flach, G.P.; Georgopoulos, P.G. A Comparative Assessment of Efficient Uncertainty Analysis Techniques for Environmental Fate and Transport Models: Application to the FACT Model. J. Hydrol. 2005, 307, 204–218. [Google Scholar] [CrossRef]
  4. Dai, H.; Chen, X.; Ye, M.; Song, X.; Zachara, J.M. A Geostatistics-Informed Hierarchical Sensitivity Analysis Method for Complex Groundwater Flow and Transport Modeling. Water Resources Research 2017, 53, 4327–4343. [Google Scholar] [CrossRef]
  5. Bianchi Janetti, E.; Guadagnini, L.; Riva, M.; Guadagnini, A. Global Sensitivity Analyses of Multiple Conceptual Models with Uncertain Parameters Driving Groundwater Flow in a Regional-Scale Sedimentary Aquifer. Journal of Hydrology 2019, 574, 544–556. [Google Scholar] [CrossRef]
  6. Maples, S.R.; Foglia, L.; Fogg, G.E.; Maxwell, R.M. Sensitivity of Hydrologic and Geologic Parameters on Recharge Processes in a Highly Heterogeneous, Semi-Confined Aquifer System. Hydrol. Earth Syst. Sci. 2020, 24, 2437–2456. [Google Scholar] [CrossRef]
  7. Zhang, X.; Ma, F.; Yin, S.; Wallace, C.D.; Soltanian, M.R.; Dai, Z.; Ritzi, R.W.; Ma, Z.; Zhan, C.; Lü, X. Application of Upscaling Methods for Fluid Flow and Mass Transport in Multi-Scale Heterogeneous Media: A Critical Review. Applied Energy 2021, 303, 117603. [Google Scholar] [CrossRef]
  8. Carrera, J.; Saaltink, M.W.; Soler-Sagarra, J.; Wang, J.; Valhondo, C. Reactive Transport: A Review of Basic Concepts with Emphasis on Biochemical Processes. Energies 2022, 15, 925. [Google Scholar] [CrossRef]
  9. He, X.; Jiang, L.; Moulton, J.D. A Stochastic Dimension Reduction Multiscale Finite Element Method for Groundwater Flow Problems in Heterogeneous Random Porous Media. J. Hydrol. 2013, 478, 77–88. [Google Scholar] [CrossRef]
  10. Tansar, H.; Duan, H.-F.; Mark, O. Global Sensitivity Analysis of Bioretention Cell Design for Stormwater System: A Comparison of VARS Framework and Sobol Method. Journal of Hydrology 2023, 617, 128895. [Google Scholar] [CrossRef]
  11. Greskowiak, J.; Prommer, H.; Liu, C.; Post, V.; Ma, R.; Zheng, C.; Zachara, J. Comparison of Parameter Sensitivities between a Laboratory and Field-Scale Model of Uranium Transport in a Dual Domain, Distributed Rate Reactive System. Water Resources Research - WATER RESOUR RES 2010, 46. [Google Scholar] [CrossRef]
  12. Abdelaziz, R.; Merkel, B.J. Sensitivity Analysis of Transport Modeling in a Fractured Gneiss Aquifer. Journal of African Earth Sciences 2015, 103, 121–127. [Google Scholar] [CrossRef]
  13. Samper, J.; Naves, A.; Montenegro, L.; Mon, A. Reactive Transport Modelling of the Long-Term Interactions of Corrosion Products and Compacted Bentonite in a HLW Repository in Granite: Uncertainties and Relevance for Performance Assessment. Applied Geochemistry 2016, 67, 42–51. [Google Scholar] [CrossRef]
  14. Montenegro, L.; Samper, J.; Mon, A.; De Windt, L.; Samper, A.-C.; García, E. A Non-Isothermal Reactive Transport Model of the Long-Term Geochemical Evolution at the Disposal Cell Scale in a HLW Repository in Granite. Applied Clay Science 2023, 242, 107018. [Google Scholar] [CrossRef]
  15. Saltelli, A.; Sobol’, I.M. About the Use of Rank Transformation in Sensitivity Analysis of Model Output. Reliability Engineering & System Safety 1995, 50, 225–239. [Google Scholar] [CrossRef]
  16. Saltelli, A.; Annoni, P.; Azzini, I.; Campolongo, F.; Ratto, M.; Tarantola, S. Variance Based Sensitivity Analysis of Model Output. Design and Estimator for the Total Sensitivity Index. Computer Physics Communications 2010, 181, 259–270. [Google Scholar] [CrossRef]
  17. Pianosi, F.; Beven, K.; Freer, J.; Hall, J.W.; Rougier, J.; Stephenson, D.B.; Wagener, T. Sensitivity Analysis of Environmental Models: A Systematic Review with Practical Workflow. Environmental Modelling & Software 2016, 79, 214–232. [Google Scholar] [CrossRef]
  18. Plischke, E. An Adaptive Correlation Ratio Method Using the Cumulative Sum of the Reordered Output. Reliability Engineering & System Safety 2012, 107, 149–156. [Google Scholar] [CrossRef]
  19. Morris, M.D. Factorial Sampling Plans for Preliminary Computational Experiments. Technometrics 1991, 33, 161–174. [Google Scholar] [CrossRef]
  20. Campolongo, F.; Cariboni, J.; Saltelli, A. An Effective Screening Design for Sensitivity Analysis of Large Models. Environmental Modelling & Software 2007, 22, 1509–1518. [Google Scholar] [CrossRef]
  21. Uddameri, V.; Hernandez, E.A.; Singaraju, S. A Successive Steady-State Model for Simulating Freshwater Discharges and Saltwater Wedge Profiles at Baffin Bay, Texas. Environ. Earth Sci. 2014, 71, 2535–2546. [Google Scholar] [CrossRef]
  22. Rabitz, H.; Aliş, Ö.F. General Foundations of High-Dimensional Model Representations. Journal of Mathematical Chemistry 1999, 25, 197–233. [Google Scholar] [CrossRef]
  23. Dai, H.; Ye, M. Variance-Based Global Sensitivity Analysis for Multiple Scenarios and Models with Implementation Using Sparse Grid Collocation. Journal of Hydrology 2015, 528, 286–300. [Google Scholar] [CrossRef]
  24. Gatel, L.; Lauvernet, C.; Carluer, N.; Weill, S.; Paniconi, C. Sobol Global Sensitivity Analysis of a Coupled Surface/Subsurface Water Flow and Reactive Solute Transfer Model on a Real Hillslope. Water 2020, 12, 121. [Google Scholar] [CrossRef]
  25. Dell’Oca, A.; Riva, M.; Guadagnini, A. Moment-Based Metrics for Global Sensitivity Analysis of Hydrological Systems. Hydrology and Earth System Sciences 2017, 21, 6219–6234. [Google Scholar] [CrossRef]
  26. Maina, F.Z.; Guadagnini, A. Uncertainty Quantification and Global Sensitivity Analysis of Subsurface Flow Parameters to Gravimetric Variations During Pumping Tests in Unconfined Aquifers. Water Resources Research 2018, 54, 501–518. [Google Scholar] [CrossRef]
  27. Razavi, S.; Gupta, H.V. A New Framework for Comprehensive, Robust, and Efficient Global Sensitivity Analysis: 1. Theory. Water Resources Research 2016, 52, 423–439. [Google Scholar] [CrossRef]
  28. Razavi, S.; Gupta, H.V. A New Framework for Comprehensive, Robust, and Efficient Global Sensitivity Analysis: 2. Application. Water Resources Research 2016, 52, 440–455. [Google Scholar] [CrossRef]
  29. Razavi, S.; Gupta, H.V. A Multi-Method Generalized Global Sensitivity Matrix Approach to Accounting for the Dynamical Nature of Earth and Environmental Systems Models. Environmental Modelling & Software 2019, 114, 1–11. [Google Scholar] [CrossRef]
  30. Puy, A.; Becker, W.; Lo Piano, S.; Saltelli, A. The Battle of Total-Order Sensitivity Estimators; 2020;
  31. Malaguerra, F.; Albrechtsen, H.-J.; Binning, P.J. Assessment of the Contamination of Drinking Water Supply Wells by Pesticides from Surface Water Resources Using a Finite Element Reactive Transport Model and Global Sensitivity Analysis Techniques. Journal of Hydrology 2013, 476, 321–331. [Google Scholar] [CrossRef]
  32. Zou, Y.; Yousaf, M.S.; Yang, F.; Deng, H.; He, Y. Surrogate-Based Uncertainty Analysis for Groundwater Contaminant Transport in a Chromium Residue Site Located in Southern China. Water 2024, 16, 638. [Google Scholar] [CrossRef]
  33. Wang, Y.; Bian, J.; Sun, X.; Ruan, D.; Gu, Z. Sensitivity-Dependent Dynamic Searching Approach Coupling Multi-Intelligent Surrogates in Homotopy Mechanism for Groundwater DNAPL-Source Inversion. J. Contam. Hydrol. 2023, 255, 104151. [Google Scholar] [CrossRef]
  34. Wang, Y.; Huang, J.; Tang, H. Global Sensitivity Analysis of the Hydraulic Parameters of the Reservoir Colluvial Landslides in the Three Gorges Reservoir Area, China. Landslides 2020, 17, 483–494. [Google Scholar] [CrossRef]
  35. Mishra, S.; Deeds, N.; Ruskauff, G. Global Sensitivity Analysis Techniques for Probabilistic Ground Water Modeling. Ground Water 2009, 47, 727–744. [Google Scholar] [CrossRef]
  36. Samper, J.; López-Vázquez, C.; Pisani, B.; Mon, A.; Samper-Pilar, A.-C.; Samper-Pilar, F.J. VARS Analysis of Reactive Transport Models for Radioactive Waste Disposal, Applied Geochemistry, (under Review). 2024.
  37. Santos, A.; Fernández, J.; Guadaño, J.; Lorenzo, D.; Romero, A. Chlorinated Organic Compounds in Liquid Wastes (DNAPL) from Lindane Production Dumped in Landfills in Sabiñanigo (Spain). Environmental Pollution 2018, 242, 1616–1624. [Google Scholar] [CrossRef] [PubMed]
  38. Sobral, B.; Samper, J.; Montenegro, L.; Mon, A.; Guadaño, J.; Gómez, J.; San Román, J.; Delgado, F.; Fernández, J. 2D Model of Groundwater Flow and Total Dissolved HCH Transport through the Gállego Alluvial Aquifer Downstream the Sardas Landfill (Huesca, Spain). Journal of Contaminant Hydrology 2024, 265, 104370. [Google Scholar] [CrossRef] [PubMed]
  39. Fernández, J.; Arjol, M.A.; Cacho, C. POP-Contaminated Sites from HCH Production in Sabiñánigo, Spain. Environ Sci Pollut Res 2013, 20, 1937–1950. [Google Scholar] [CrossRef]
  40. Biosca, B.; Arévalo-Lomas, L.; Izquierdo-Díaz, M.; Díaz-Curiel, J. Detection of Chlorinated Contaminants Coming from the Manufacture of Lindane in a Surface Detritic Aquifer by Electrical Resistivity Tomography. Journal of Applied Geophysics 2021, 191, 104358. [Google Scholar] [CrossRef]
  41. Samper, J.; Sobral, B.; Pisani, B.; Naves, A.; Guadaño, J.; Gómez, J.; Fernández, J. Groundwater Flow Model along a Vertical Profile of the Sardas Landfill in Sabiñánigo, Huesca, Spain. Water 2023, 15, 3457. [Google Scholar] [CrossRef]
  42. Guadaño, J.; Gómez, J.; Fernández, J.; Lorenzo, D.; Domínguez, C.M.; Cotillas, S.; García-Cervilla, R.; Santos, A. Remediation of the Alluvial Aquifer of the Sardas Landfill (Sabiñánigo, Huesca) by Surfactant Application. Sustainability 2022, 14, 16576. [Google Scholar] [CrossRef]
  43. Pankow, J.F.; Cherry, J.A. Dense Chlorinated Solvents and Other DNAPLs in Groundwater: History, Behavior, and Remediation; Waterloo Press: Portland, Or, 1996; ISBN 978-0-9648014-1-7. [Google Scholar]
  44. Casado, I.; Mahjoub, H.; Lovera, R.; Fernández, J.; Casas, A. Use of Electrical Tomography Methods to Determinate the Extension and Main Migration Routes of Uncontrolled Landfill Leachates in Fractured Areas. Science of The Total Environment 2015, 506–507, 546–553. [Google Scholar] [CrossRef]
  45. Julià, X.; González, G. ; Alonso, M Estudio Batimétrico y de Caracterización de Sedimentos Del Embalse de Sabiñánigo. URS. Technical Report for the Ebre River Water District. Zaragoza 2009.
  46. Shuai, P.; Chen, X.; Song, X.; Hammond, G.E.; Zachara, J.; Royer, P.; Ren, H.; Perkins, W.A.; Richmond, M.C.; Huang, M. Dam Operations and Subsurface Hydrogeology Control Dynamics of Hydrologic Exchange Flows in a Regulated River Reach. Water Resources Research 2019, 55, 2593–2612. [Google Scholar] [CrossRef]
  47. Espinha Marques, J.; Samper, J.; Pisani, B.; Alvares, D.; Carvalho, J.M.; Chaminé, H.I.; Marques, J.M.; Vieira, G.T.; Mora, C.; Sodré Borges, F. Evaluation of Water Resources in a High-Mountain Basin in Serra Da Estrela, Central Portugal, Using a Semi-Distributed Hydrological Model. Environ Earth Sci 2011, 62, 1219–1234. [Google Scholar] [CrossRef]
  48. Plischke, E.; Röhlig, K.-J. Methodological Approaches to Uncertainty and Sensitivity Analysis. Final Version as of 07.05.2024 of Deliverable D10.4 of the HORIZON 2020 Project EURAD. EC Grant Agreement No: 847593. Available online: https://www.ejp-eurad.eu/publications/eurad-d104-methodological-approaches-uncertainty-and-sensitivity-analysis (accessed on 2 August 2024).
  49. Sobol’, I.M. On the Distribution of Points in a Cube and the Approximate Evaluation of Integrals. USSR Computational Mathematics and Mathematical Physics 1967, 7, 86–112. [Google Scholar] [CrossRef]
  50. Sobol’, I.M. Sensitivity Estimates for Nonlinear Mathematical Models. Matematicheskoe Modelirovanie 1990, 2, 112–118. [Google Scholar]
  51. Sobol’, I.M. Sensitivity Analysis for Non-Linear Mathematical Models. Mathematical Modelling and Computational Experiment 1993, 1, 407–414, [English translation of the original paper of Sobol (1990) in Russian]. [Google Scholar]
  52. Homma, T.; Saltelli, A. Importance Measures in Global Sensitivity Analysis of Nonlinear Models. Reliability Engineering & System Safety 1996, 52, 1–17. [Google Scholar] [CrossRef]
  53. Razavi, S.; Sheikholeslami, R.; Gupta, H.V.; Haghnegahdar, A. VARS-TOOL: A Toolbox for Comprehensive, Efficient, and Robust Sensitivity and Uncertainty Analysis. Environmental Modelling & Software 2019, 112, 95–107. [Google Scholar] [CrossRef]
  54. Halton, J.H. On the Efficiency of Certain Quasi-Random Sequences of Points in Evaluating Multi-Dimensional Integrals. Numer. Math. 1960, 2, 84–90. [Google Scholar] [CrossRef]
  55. Krause, P.; Boyle, D.P.; Bäse, F. Comparison of Different Efficiency Criteria for Hydrological Model Assessment. Adv. Geosci. 2005, 5, 89–97. [Google Scholar] [CrossRef]
  56. Samper, J.; Xu, T.; Yang, C. A Sequential Partly Iterative Approach for Multicomponent Reactive Transport with CORE2D. Comput Geosci 2009, 13, 301–316. [Google Scholar] [CrossRef]
  57. Poonoosamy, J.; Wanner, C.; Alt Epping, P.; Águila, J.F.; Samper, J.; Montenegro, L.; Xie, M.; Su, D.; Mayer, K.U.; Mäder, U.; et al. Benchmarking of Reactive Transport Codes for 2D Simulations with Mineral Dissolution–Precipitation Reactions and Feedback on Transport Parameters. Comput Geosci 2021, 25, 1337–1358. [Google Scholar] [CrossRef]
  58. Águila, J.F.; Montoya, V.; Samper, J.; Montenegro, L.; Kosakowski, G.; Krejci, P.; Pfingsten, W. Modeling Cesium Migration through Opalinus Clay: A Benchmark for Single- and Multi-Species Sorption-Diffusion Models. Comput Geosci 2021, 25, 1405–1436. [Google Scholar] [CrossRef]
  59. Naves, A.; Samper, J.; Pisani, B.; Mon, A.; Dafonte, J.; Montenegro, L.; García-Tomillo, A. Hydrogeology and Groundwater Management in a Coastal Granitic Area with Steep Slopes in Galicia (Spain). Hydrogeol J 2021, 29, 2655–2669. [Google Scholar] [CrossRef]
  60. Mon, A.; Samper, J.; Montenegro, L.; Turrero, M.J.; Torres, E.; Cuevas, J.; Fernández, R.; De Windt, L. Reactive Transport Models of the Geochemical Interactions at the Iron/Bentonite Interface in Laboratory Corrosion Tests. Applied Clay Science 2023, 240, 106981. [Google Scholar] [CrossRef]
  61. González, P.; Prado-Rodriguez, R.; Gábor, A.; Saez-Rodriguez, J.; Banga, J.R.; Doallo, R. Parallel Ant Colony Optimization for the Training of Cell Signaling Networks. Expert Systems with Applications 2022, 208, 118199. [Google Scholar] [CrossRef]
  62. Vourlioti, P.; Kotsopoulos, S.; Mamouka, T.; Agrafiotis, A.; Nieto, F.J.; Sánchez, C.F.; Llerena, C.G.; García González, S. Maximizing the Potential of Numerical Weather Prediction Models: Lessons Learned from Combining High-Performance Computing and Cloud Computing. In Proceedings of the Advances in Science and Research; Copernicus GmbH, March 20 2023; Vol. 20, pp. 1–8.
  63. Herman, J.; Usher, W. SALib: An Open-Source Python Library for Sensitivity Analysis. Journal of Open Source Software 2017, 2, 97. [Google Scholar] [CrossRef]
  64. Ziehn, T.; Tomlin, A.S. GUI–HDMR – A Software Tool for Global Sensitivity Analysis of Complex Models. Environmental Modelling & Software 2009, 24, 775–785. [Google Scholar] [CrossRef]
  65. Kocis, L.; Whiten, W.J. Computational Investigations of Low-Discrepancy Sequences. ACM Trans. Math. Softw. 1997, 23, 266–294. [Google Scholar] [CrossRef]
Preprints 116249 g001
Figure 1. Flowchart of the methodology used in this study.
Figure 1. Flowchart of the methodology used in this study.
Preprints 116249 g001
Preprints 116249 g002
Figure 2. (a) Location of the study area; (b) enlargement showing the model domain, the Sabiñánigo reservoir, the Sardas landfill, the Gállego river course and the INQUINOSA former production site.
Figure 2. (a) Location of the study area; (b) enlargement showing the model domain, the Sabiñánigo reservoir, the Sardas landfill, the Gállego river course and the INQUINOSA former production site.
Preprints 116249 g002
Preprints 116249 g003
Figure 3. 2D finite element mesh, monitoring wells, material zones, boundary conditions, and GSA input parameters (top plot) and enlargement showing the area downstream the Sardas landfill (bottom plot). The confined storage coefficient (SS) is the same in the four material zones. The sands and gravels are assumed to be confined in the alluvial (rc), except in the wooded areas (ru). Unconfined areas are shown with a back hashed polygon.
Figure 3. 2D finite element mesh, monitoring wells, material zones, boundary conditions, and GSA input parameters (top plot) and enlargement showing the area downstream the Sardas landfill (bottom plot). The confined storage coefficient (SS) is the same in the four material zones. The sands and gravels are assumed to be confined in the alluvial (rc), except in the wooded areas (ru). Unconfined areas are shown with a back hashed polygon.
Preprints 116249 g003
Preprints 116249 g004
Figure 4. Map showing the reservoir tail area (hashed blue polygon) where aquifer/reservoir fluxes were calculated at times t1, t2 and t3, the monitoring wells whose piezometric data were used to calculate the calibration metrics, monitoring wells ST1C, PS19B, SPN1 and PS16C (where the average Darcy velocity is computed).
Figure 4. Map showing the reservoir tail area (hashed blue polygon) where aquifer/reservoir fluxes were calculated at times t1, t2 and t3, the monitoring wells whose piezometric data were used to calculate the calibration metrics, monitoring wells ST1C, PS19B, SPN1 and PS16C (where the average Darcy velocity is computed).
Preprints 116249 g004
Preprints 116249 g005
Figure 5. Measured reservoir hydrograph and piezometric heads in well ST1C from September 18, 2020, to September 20, 2020. The computed piezometric heads in monitoring wells ST1C, PS19B and SPN1 and the aquifer/reservoir fluxes are analyzed at the following times: 1) t1, September 18, 2020, 20:00 (low reservoir water level), 2) t2, September 18, 2020, 22:30 (peak reservoir water level) and 3) t3, September 19, 2020, 04:30 (descending reservoir water level).
Figure 5. Measured reservoir hydrograph and piezometric heads in well ST1C from September 18, 2020, to September 20, 2020. The computed piezometric heads in monitoring wells ST1C, PS19B and SPN1 and the aquifer/reservoir fluxes are analyzed at the following times: 1) t1, September 18, 2020, 20:00 (low reservoir water level), 2) t2, September 18, 2020, 22:30 (peak reservoir water level) and 3) t3, September 19, 2020, 04:30 (descending reservoir water level).
Preprints 116249 g005
Preprints 116249 g006
Figure 6. Scatterplots of the computed piezometric heads in wells ST1Ct2 (upper left plot), PS19Bt2 (upper right plot), SPN1t2 (lower left plot), and Qt2 (lower right plot) versus the vertical hydraulic conductivity of the silting sediments in the former river course (Kvs1). The sample of 16384 points was generated with a Sobol sequence. The clouds of plots are shown for the following three ranges of percentiles, p, of the specific storage coefficient (SS): 1) p < 30%; 2) 30% < p < 70% and 3) p > 70%.
Figure 6. Scatterplots of the computed piezometric heads in wells ST1Ct2 (upper left plot), PS19Bt2 (upper right plot), SPN1t2 (lower left plot), and Qt2 (lower right plot) versus the vertical hydraulic conductivity of the silting sediments in the former river course (Kvs1). The sample of 16384 points was generated with a Sobol sequence. The clouds of plots are shown for the following three ranges of percentiles, p, of the specific storage coefficient (SS): 1) p < 30%; 2) 30% < p < 70% and 3) p > 70%.
Preprints 116249 g006
Preprints 116249 g007
Figure 7. CUSUNORO curves of computed head in wells ST1C and PS19B at times t1, t2 and t3; and well SPN1 at times t1 and t2.
Figure 7. CUSUNORO curves of computed head in wells ST1C and PS19B at times t1, t2 and t3; and well SPN1 at times t1 and t2.
Preprints 116249 g007
Preprints 116249 g008
Figure 8. CUSUNORO curves of the computed head in well SPN1 at time t3, MAEg, NRMSEg, NSEg, Qt1, Qt2, Qt3 and qav.
Figure 8. CUSUNORO curves of the computed head in well SPN1 at time t3, MAEg, NRMSEg, NSEg, Qt1, Qt2, Qt3 and qav.
Preprints 116249 g008
Preprints 116249 g009
Figure 9. IVARS50 indexes of input parameters as a function of the number of star centers for MAEg (upper plot), and robustness of ranking as a function of the number of star centers (bottom plot).
Figure 9. IVARS50 indexes of input parameters as a function of the number of star centers for MAEg (upper plot), and robustness of ranking as a function of the number of star centers (bottom plot).
Preprints 116249 g009
Preprints 116249 g010
Figure 10. IVARS50 indexes of input parameters as a function of the number of star centers for the average Darcy velocity (qav) (upper plot), and robustness of ranking as a function of the number of star centers (bottom plot).
Figure 10. IVARS50 indexes of input parameters as a function of the number of star centers for the average Darcy velocity (qav) (upper plot), and robustness of ranking as a function of the number of star centers (bottom plot).
Preprints 116249 g010
Preprints 116249 g011
Figure 11. Sample variograms of the computed heads in monitoring wells ST1C and PS19B at times t1, t2 and t3 and monitoring well SPN1 at times t1 and t2. Only the variograms of the five most influential parameters are shown in the plots.
Figure 11. Sample variograms of the computed heads in monitoring wells ST1C and PS19B at times t1, t2 and t3 and monitoring well SPN1 at times t1 and t2. Only the variograms of the five most influential parameters are shown in the plots.
Preprints 116249 g011
Preprints 116249 g012
Figure 12. Sample variograms of the computed head in well SPN1 at time t3, MAEg, NRMSEg, NSEg, Qt1, Qt2, Qt3 and qav. Only the variograms of the five most influential parameters are shown in the plots.
Figure 12. Sample variograms of the computed head in well SPN1 at time t3, MAEg, NRMSEg, NSEg, Qt1, Qt2, Qt3 and qav. Only the variograms of the five most influential parameters are shown in the plots.
Preprints 116249 g012
Preprints 116249 g013
Figure 13. VARS-TO, IVARS50 and VARS-ABE indexes for the computed heads in wells ST1C and PS19B at times t1, t2 and t3 and well SPN1 at times t1 and t2.
Figure 13. VARS-TO, IVARS50 and VARS-ABE indexes for the computed heads in wells ST1C and PS19B at times t1, t2 and t3 and well SPN1 at times t1 and t2.
Preprints 116249 g013
Preprints 116249 g014
Figure 14. VARS-TO, IVARS50 and VARS-ABE indexes for the computed head in well SPN1 at time t3, MAEg, NRMSEg, NSEg, Qt1, Qt2, Qt3 and qav.
Figure 14. VARS-TO, IVARS50 and VARS-ABE indexes for the computed head in well SPN1 at time t3, MAEg, NRMSEg, NSEg, Qt1, Qt2, Qt3 and qav.
Preprints 116249 g014
Preprints 116249 g015
Figure 15. IVARS50 sensitivity indexes for computed heads in wells ST1C (top left plot), PS19B (top right plot), and SPN1 (bottom left plot) and aquifer/reservoir flow (bottom right plot) at times t1, t,2 and t3.
Figure 15. IVARS50 sensitivity indexes for computed heads in wells ST1C (top left plot), PS19B (top right plot), and SPN1 (bottom left plot) and aquifer/reservoir flow (bottom right plot) at times t1, t,2 and t3.
Preprints 116249 g015
Preprints 116249 g016
Figure 16. IVARS50 sensitivity indexes for calibration metrics MAEg, NRMSEg and NSEg.
Figure 16. IVARS50 sensitivity indexes for calibration metrics MAEg, NRMSEg and NSEg.
Preprints 116249 g016
Table 1. Ranges and statistical distributions of the input parameters.
Table 1. Ranges and statistical distributions of the input parameters.
Parameter Minimum Maximum Unit Distribution
Aquifer conductivity K1 10 103 m/d Log-uniform
Aquifer conductivity K2 10 103 m/d Log-uniform
Aquifer conductivity K3 10 103 m/d Log-uniform
Aquifer conductivity K4 10 103 m/d Log-uniform
Storage coefficient Ss 10-5 10-3 1/m Log-uniform
Aquitard conductivity KVs1 10-3 1 m/d Log-uniform
Aquitard conductivity KVs2 10-3 1 m/d Log-uniform
Aquitard conductivity KVs3 10-4 10-1 m/d Log-uniform
Leakage coefficient αr 10 103 m2/d Log-uniform
Conductance αd 1 100 m2/d Log-uniform
Boundary inflow Q6 3·10-3 0.05 m3/d/m Uniform
Boundary inflow Q7 2·10-3 0.20 m3/d/m Uniform
Boundary inflow Q9 0.25 1.00 m3/d/m Uniform
Boundary inflow Q2 1.70·10-2 1.70 m3/d/m Uniform
Boundary inflow Q1 2·10-3 10-1 m3/d/m Uniform
Recharge rc 5 200 mm/year Uniform
Recharge ru 20 401.5 mm/year Uniform
Table 2. Ranking comparison of sensitivity results across various methods for the computed piezometric head in well ST1C at time t1 (ST1Ct1), the global mean absolute error (MAEg), the aquifer/reservoir flux at time t1 (Qt1), and the average groundwater Darcy velocity modulus (qav).
Table 2. Ranking comparison of sensitivity results across various methods for the computed piezometric head in well ST1C at time t1 (ST1Ct1), the global mean absolute error (MAEg), the aquifer/reservoir flux at time t1 (Qt1), and the average groundwater Darcy velocity modulus (qav).
Output Methods K1 K2 K3 K4 Ss KVs1 KVs2 KVs3 αr αd Q6 Q7 Q9 Q2 Q1 rc ru
ST1Ct1 VARS-TO 17 4 1 6 13 3 8 7 5 9 16 11 14 2 10 12 15
IVARS50 16 4 1 6 12 3 8 7 5 9 17 11 14 2 10 13 15
VARS-ABE 17 5 1 6 9 4 11 10 3 8 16 12 14 2 7 13 15
SALib 17 5 1 10 6 3 13 7 4 8 15 12 14 2 9 11 16
GUI-HDMR - - - - - - - - - - - - - - - - -
CUSUNORO 17 5 1 6 13 3 11 9 4 7 16 10 14 2 8 12 15
MAEg VARS-TO 14 3 1 7 16 5 8 6 4 9 17 12 13 2 11 10 15
IVARS50 14 4 1 7 16 5 8 6 3 9 17 12 13 2 11 10 15
VARS-ABE 14 4 1 6 16 5 11 8 3 10 17 12 13 2 9 7 15
SALib 14 5 1 6 17 4 12 7 3 10 16 11 13 2 9 8 15
GUI-HDMR 14 4 1 7 17 5 11 6 3 8 15 12 13 2 10 9 16
CUSUNORO 14 5 1 6 17 4 12 7 3 10 16 11 13 2 8 9 15
Qt1 VARS-TO 17 9 4 6 5 1 14 12 2 8 16 7 10 3 11 13 15
IVARS50 17 9 4 6 5 1 14 11 2 7 16 8 10 3 12 13 15
VARS-ABE 17 11 5 7 4 1 14 13 2 8 16 6 12 3 9 10 15
SALib 17 9 13 10 6 1 15 12 2 5 16 4 11 3 7 8 14
GUI-HDMR 17 9 5 8 4 1 14 13 2 7 16 6 11 3 10 12 15
CUSUNORO 17 8 13 10 6 1 15 12 2 5 16 4 11 3 7 9 14
qav VARS-TO 4 3 2 7 16 5 9 10 6 8 17 14 13 1 11 12 15
IVARS50 4 3 2 8 16 6 9 10 5 7 17 14 13 1 11 12 15
VARS-ABE 4 3 2 7 16 6 12 11 5 9 17 14 13 1 8 10 15
SALib 4 3 2 9 16 6 13 10 5 8 15 17 12 1 7 11 14
GUI-HDMR 4 3 2 7 14 6 12 10 5 8 15 17 13 1 9 11 16
CUSUNORO 4 3 2 8 17 6 13 11 5 9 16 15 12 1 7 10 14
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated