Abstract

Growing demands in arid regions have increased groundwater vulnerabilities necessitating appropriate modeling and management strategies to understand and sustain aquifer system behaviors. Sustainable management of aquifer systems, however, requires a proper understanding of groundwater dynamics and accurate estimates of recharge rates which often cause error and uncertainty in simulation. This study aims to quantify the uncertainty and error associated with groundwater simulation using various multi-model ensemble averaging (MEA) techniques such as simple model averaging, weighted averaging model, multi-model super ensemble, and modified MMSE. Two numerical solutions, i.e., finite difference and finite element (FE), were first coupled under three schemes such as explicit scheme (ES), implicit scheme, and Crank-Nicolson Scheme to numerically solve groundwater simulation problems across two case studies, synthetic and real-world (the Birjand aquifer in Iran) case studies. The MEA approach was considerably successful in calibrating a complex arid aquifer in a way that honors complex geological heterogeneity and stress configurations. Specifically, the MEA techniques skillfully reduced the error and uncertainties in simulation, particularly those errors associated with water table variability and fluctuation. Furthermore, a coupled FE-ES method outperformed other approaches and generated the best groundwater-level simulation for the synthetic case study, while stand-alone FE was particularly successful for the Birjand aquifer simulation as a real-world case study.

HIGHLIGHTS

  • Multi-model ensemble averaging (MEA) techniques are employed for improving groundwater numerical simulation.

  • Two finite differences and finite element models developed to simulate groundwater.

  • Numerical models' proficiency was derived in the different complex conditions obtained.

  • MEA approaches show promise in enhancing groundwater flow.

  • Some MEA techniques having weakness in estimating numerical models’ contributions.

INTRODUCTION

Accurate simulation of the groundwater flow processes across various spatial and temporal scales is critical to reach efficient, timely, and sustainable groundwater management, particularly in the regions with water scarcity problems such as arid areas (e.g., Krakauer et al. 2014). In these regions, groundwater resources are under a constant threat of being overdrawn due to irrigation, anthropogenic, and climate pressure among others (Mustafa et al. 2020). With these stressors, the lack of perennial sources of groundwater has been a concern for domestic water supply and economic prosperity in these areas. To provide useful guidelines for the stakeholders, modeling-based solutions and recommendations have been proposed to manage arid aquifer systems and reduce their vulnerabilities.

So far, several groundwater models have been developed to simulate aquifer flow including lumped, semi-distributed, and distributed models. The lumped models are simplified mathematical domains of an aquifer system that do not require incorporating physical parameters into the modeling process. These models indeed rely on the theoretical understanding of groundwater flow to simulate the fundamental physical concepts. While less computationally burden and fewer numbers of parameters for optimization make lumped approaches attractive for groundwater simulation (e.g., Mackay et al. 2014), these models often commit a certain amount of errors that challenge groundwater simulation.

As technologies advanced, more physics-based models were developed and applied for various case studies (e.g., Ajami et al. 2004). In these models, the physical principles of fluid dynamics are simplified numerically based on two- or three-dimensional spatial grids of groundwater flow (Mackay et al. 2014), although these models require a proper understanding of aquifer complex processes, including anisotropy and heterogeneity of the system (e.g., Shepley et al. 2012). Among several physics-based/distributed approaches, the finite difference (FD) and finite element (FE) methods are among the most widely used numerical solutions for groundwater studies. FD is the simplest numerical method that is widely used in groundwater studies. For example, Pollock (1988) applied FD to simulate particle tracing, while Narasimhan & Witherspoon (1976) employed this model to handle the multidimensional systems for an aquifer composed of heterogeneous material. Soon after, a MODular groundwater FLOW (MODFLOW) was developed in Fortran programing language by McDonald & Harbaugh (1984) based on the FD that was widely implemented on some interfaces, including Groundwater Modeling System, Processing Modflow for Windows (PMWIN), and Groundwater Vistas.

The successful application and widespread use of MODFLOW have made this model a benchmark approach to simulate groundwater systems (e.g., Guiguer & Franz 1996). So far, MODFLOW has been applied in various case studies such as stream depletion simulation of the aquifer (Ou et al. 2016; Minderhoud et al. 2017; Fletcher et al. 2019), calibration of hydrodynamic components (Xu et al. 2017), evaluation of the interaction effects of aquifer recharge and evapotranspiration (Maquin et al. 2017), simulation of drawdown in the inclined wells (Batu 2015), assessment of aquifer vulnerability (Pacheco et al. 2018), flood propagation and infiltration (Li et al. 2019), among others. However, the MODFLOW model can perform best in regular geometry with proper grid flexibility, whereas for real-world applications with irregular geometry, this approach may poorly perform due to its governing flow approximation which leads to structural uncertainty in simulation (e.g., Strack 2017). An excellent application of MODFLOW structural uncertainty was provided by Faust & Mercer (1980) and Alvarez. & Illman (2005) which led to designing a mesh consisting of arranged nodes with an equal distance, although Istok (1989) found designing a mesh with arranged nodes inappropriate for irregular aquifer simulation and suggested including some dummy columns/rows into the mesh and estimate additional nodes manually.

To address this shortcoming, other numerical methods such as the FE method were applied in groundwater simulation. The first application of the FE method is back to the early 1970s when Zienkiewicz et al. (1966) used this method to estimate seepage in an anisotropic aquifer. Subsequently, Javandel & Witherspoon (1968) use FE to simulate groundwater flow under transient conditions. Since then, many case studies have been conducted using the FE such as flow simulation under transient and steady states (Mahdavi 2019), simulation of groundwater level (Kim & Lee 2019), calibration of hydrodynamic components (Berglund et al. 2020), simulation of saltwater intrusion (Sowe et al. 2020), climate change impacts on groundwater recharge (Lamontagne-Hallé et al. 2018), modeling of solute transport (Hosseini et al. 2020), among others. Unlike the FD, the FE can simulate irregular geometry, without defining dummy nodes, although it might suffer from time-consuming and complex preprocess practices (computational burden) including the arrangement of the elements in the problem domain. Additionally, the distributed elements on the domain of irregular problems do not completely cover the entire area of the aquifer, so some parts of the domain may be dismissed in computation.

One of the most promising solutions to enhance simulation is through the use of a multi-model ensemble averaging (MEA) model. In this algorithm, the outputs of several models can be potentially combined to improve the accuracy of the simulation. The MEA methods comprise two major categories, namely statistical and probabilistic approaches. In the statistical domain (e.g., MMSE method or linear and nonlinear transfer functions), the MEA is computed as a weighted average of several individual models and their deviation (i.e., a statistical formulation), while a probability concept has been used as the probabilistic MEA.

MEA has been around since the 1970s when Bates & Granger (1969), Dickinson (1973), and Thompson (1977) used this algorithm in economic and weather forecasting assessment. Afterward, Shamseldin et al. (1997) applied this algorithm for rainfall-runoff simulation. They compared the MEA simulation with other methods such as the simple model average (SMA), weighted average method (WAM), and artificial neural network and found that the MEA could better describe runoff processes and simulation. Following these results, Krishnamurti et al. (1999) advanced MEA techniques by proposing a new approach called multi-model super ensemble (MMSE). They showed that MMSE has much higher performance compared with the conventional MEA. The MEA advancement continued by Ajami et al. (2006) who introduced a new method called Modified multi-model super ensemble (M3SE) and compared its performance against three other combination methods (SMA, WAM, and MMSE) across five watersheds in the USA.

Following successful implementations of the MEA algorithm, a technical question has been arisen concerning how to skillfully decrease different sources of uncertainties in groundwater modeling through MEA application?

Owing to the above studies, this research aims to address this question and to quantify the structural uncertainty arising in the groundwater simulation model using synthetic and real-world case studies. Our proposed approaches not only complement prior studies but also sheds light on how variability in flow properties can impact aquifer hydrodynamic characteristics. Building on our prior contributions to groundwater simulation research (see Hamraz et al. 2015; Sadeghi-Tabas et al. 2017), we programmed both FE and FD algorithms in the MATrix LABoratory (MATLAB) environment and coupled them with the MEA method. To the best of our knowledge, very few studies have been conducted on groundwater simulation using MEA. The exceptions are Barzegar et al. (2018) and Mustafa et al. (2020) who coupled the DRASTIC (Depth to water, Recharge, Aquifer media, Soil media, Topography, Impact of the vadose zone, Conductivity) and the PMWIN models with MEA to address aquifer vulnerability. While these studies shed light on the MEA implementation, they largely neglected structural uncertainty during groundwater simulation. To address this shortcoming, this study developed both FD and FE as two numerical solutions to simulate the groundwater level in a synthetic test case and a real-world application. Various MEA methods such as SMA, WAM, MMSE, and M3SE were coupled with FD and FE to improve the performance of numerical results and reduce structural uncertainty. Furthermore, the spatial variability of multiple models was considered and their weights on the final simulation results were achieved. The weighted multi-model ensembles were then validated across multiple observed wells in the study area.

This paper is organized as follows. In section 2, the numerical solution for groundwater modeling along with the MEA and case studies are provided. Section 3 presents the results and provides the numerical model application and MEA performance assessment. Section 4 summarizes the results and discusses future directions for groundwater modeling research.

METHODS

Groundwater flow-governing equations

Simultaneous application of mass balance and Darcy's law creates the governing equation of groundwater flow. The equation of two-dimensional, isotropic, and the homogenous aquifer is given by the following equation (see Arnold et al. 1993):
formula
(1)
Subject to the following auxiliary conditions:
formula
(2)
where h, Sy, T, and t are, respectively, potential head, specific yield, transmissivity in horizontal and vertical directions (m2/day), and time (day). Q, q, and are, respectively, source or sink function (m3/day), distributed rate of recharge or evapotranspiration over aquifer domain (m/day), and Dirac Delta function. , , and are known, respectively, as inflow rate (m3/day), constant head boundary (m), and initial head (m). Furthermore, , , and reflect the global, essential (dirichlet), and natural (Neuman) boundary conditions, respectively, while denotes the aquifer domain.
This study used forward difference approximation to estimate time derivative of potential head as well as the method of solution based on Equation (3) (Owais et al. 2008):
formula
(3)
represents the explicit scheme (ES), while the implicit scheme (IS) is produced by . In addition, the Crank-Nicolson (CN) is obtained by setting in Equation (3).

Finite differences

The simplest numerical approach for solving groundwater flow equations is using the FD approach, which was proposed by Richtmyer & Morton (1994). FD uses the Taylor series approximation to define the discretization of the derivatives of the flow variables. By considering forward difference approximation for time derivatives of the potential head, FD approximation can be expressed as follows:
formula
(4)

If takes 1, 0.5, and zero, Equation (4) is then converted, respectively, to IS-FD, CN-FD, and ES-FD. Readers are referred to Wang & Anderson (1995) for more details about the FD.

Finite element

FE is based on the methods of weighted residuals (MWR), in which the approximated solutions are obtained by a finite series of weight functions such that the residuals are forced to zero. The MWR makes the residual as small as possible using some average techniques such as setting weighted integrals of residuals to zero. The choice of weight function is performed based on some approaches such as Galerkin, least square, Petrov Galerkin, collection method, and subdomain. In the Galerkin method employed in the FE, the weight function is assumed to be equal to the shape function. This study used a linear triangular element due to its good performance and simplicity of implementation. The first step in the FE formulation is a trial solution definition:
formula
(5)
where is a trial solution, n is the node size in the aquifer mesh, is the node's head, and is weight function (which is equal to shape function). Formulating this equation requires determining values within the aquifer domain. By replacing the trial solution into Equation (1) and using , the integration of the weighted residuals may be expressed as the following equation:
formula
(6)
Multiplying by all internal components, a second-order approximation of the space derivatives in the trial solution can be decreased to a first order. The shape function can be then used to interpolate the potential head for each element. Equation (7) shows the final matrix form of FE:
formula
(7)
where [G] is conductance matrix, [P] represents mass matrix, {F} is a boundary flux vector, and {B} denotes load vector. Finally, the potential head of groundwater level for various schemes such as implicit, explicit, and CN can be expressed using the following equation:
formula
(8)

If is equal to one, then the above equation becomes IS-FE. Also, is equal to 0.5 and zero, respectively, for the ES-FE and the CS-FE cases. More details about the FE formulation may be found at Pinder & Gray (1977), Wang & Anderson (1995), and Desai et al. (2011). The schematic layout of the FE method is displayed in Figure 1.

Figure 1

Flowchart for the FE-based groundwater flow model.

Figure 1

Flowchart for the FE-based groundwater flow model.

Case studies

We used two case studies to evaluate the performance of groundwater numerical models. The synthetic test case contains required materials like the conceptual model, parameters, and input data with true values to test the applicability and validity of developed numerical models. Also, a real-world case was applied to address the ability of proposed numerical models in a complied problem. Low rainfall amount, poor recharge condition, and intense overexploitation caused severe deficits in groundwater storage. These challenges motivated the authors to test various numerical solutions for real-world conditions.

Synthetic test case

We first used a simple case to test FE and FD simulation performances. We employed a rectangular aquifer as a synthetic case study introduced by Illangasekare & Döll (1989) and then compared the aquifer's analytical and simulated values (see Figure 2). This synthetic aquifer has two no-flow boundaries on the left and right sides and two constant boundaries with a value of 100 m in the top and bottom portions. The hydrodynamic parameters such as specific yield and transmissivity were also considered 0.15 and 885.71 m2/day, respectively. Groundwater table drawdown was simulated for 210 consecutive days. The well located at the 85th node was considered for the calibration section (W85), while W129, W174, W156, W66, and W68 were considered for validation wells.

Figure 2

Schematic plan view of the hypothetic aquifer (synthetic case study).

Figure 2

Schematic plan view of the hypothetic aquifer (synthetic case study).

Field test case

Here, different characteristics of the Birjand aquifer are introduced as a real-world case study to demonstrate the performance of proposed numerical models. The Birjand plain is in the East part of Iran – an arid region with low annual rainfall and high-temperature amounts (see Figure 3). In this region, the groundwater system is the main source of water supply including drinking water, irrigation, etc.

Figure 3

Location of Birjand plain and the aquifer in Iran.

Figure 3

Location of Birjand plain and the aquifer in Iran.

To simulate groundwater flow in the Birjand aquifer, a grid model comprising 1,175 nodes with a vertical and horizontal cell size of 500 m, and 34 rows, as well as 96 columns, were imposed into the numerical model. The aquifer consists of nine input sections in the south, northeast, and northwest and one output section in the southwest (Figure 4; see Sadeghi-Tabas et al. 2017). The boundary conditions in these sections are considered as constant head values and the rest was considered as the no-flow boundary. Around 191 extraction wells have been drilled into a deep, confined aquifer to supply different kinds of demands. There are overall 11 piezometric wells in the Birjand aquifer that are programmed based on mesh point's numbers, i.e., piezometric wells 13, 53, 85, 212, 393, 568, 631, 644, 760, 811, and 995 (hereafter P13, P53, P85 P212, P393, P568, P631, P644, P760, P811, and P995; see Figure 3). The groundwater levels are monthly measured in these wells to monitor hydraulic head fluctuations and flow dynamics. We constructed the Birjand geological structure in the model based on our prior studies (see Sadeghi-Tabas et al. 2017). Specifically, we followed Sadeghi-Tabas et al. (2017) and used their proposed approaches for geological heterogeneity and anisotropy structure in the model.

Figure 4

Definition of hydraulic conductivity (a) and specific yield (b) zones across the Birjand aquifer (adopted from Sadeghi-Tabas et al. 2017).

Figure 4

Definition of hydraulic conductivity (a) and specific yield (b) zones across the Birjand aquifer (adopted from Sadeghi-Tabas et al. 2017).

Ensemble methods

Since averaging methods is the most fundamental combination method for numerical individual models (Zhou 2012; Shamshirband et al. 2019), the current study aimed to examine the effectiveness of some averaging techniques described as the following.

Simple model average

SMA is the simplest ensemble method proposed by Georgakakos et al. (2004). SMA formulation for any time step t is expressed as the following equation:
formula
(9)
where is the produced drawdown using SMA, m is the number of individual competing models (e.g., IS-FE), is the simulated drawdown of ith numerical method at time t, is the mean of the simulated drawdown in ith numerical method, and is the averaged analytical drawdown.

Weighted averaging model

The WAM presents a weighted average combining all numerical models based on the multiple linear regression (MLR) fit, while all weights can be positive, and their summation must be one. The WAM formulation can be expressed as:
formula
(10)
where is the generated drawdown by the WAM and is the weight of numerical method. Also, denotes the analytical solution for synthetic case or observed values for real-world aquifer.

Multi-model super ensemble

In this method, weights are calculated based on the same procedure as the WAM with the exception that weights are not constrained in the MMSE. This means weights can capture all possible real numbers, and there is not a specific limitation for their summation. The logic of MMSE can be defined using the following equation:
formula
(11)
where is the MMSE's drawdown at time t.

Modified multi-model super ensemble

M3SE was introduced by Ajami et al. (2006), and it is based on an interesting modification of MMSE. This method was equipped with frequency mapping for bias correction to improve the MMSE procedure. Hence, a frequency mapping is imposed on simulations of all competing models, and then a procedure like MMSE is used. Readers are referred to Ajami et al. (2006) for more information.

Testing multi-collinearity problem in simulation

The MLR is an approach to define weights in the model which is the basis for ensemble method construction. When there is a strong linear correlation between independent variables, the MLR likely estimates the unstable weight for predictors. This shortcoming is addressed in this study by including a new time-series data to the contributing models. This new time series is a multiplier of the analytical solution, and it was added to the computation as a new ensemble member to capture average model weight more precisely. All factors were computed based on a normal distribution with a mean of 0.95 and a variance of 0.01 (see Vrugt et al. 2005).

Groundwater modeling framework

This section describes our developed modeling framework for groundwater modeling assessment. We first discretized the governing equation of groundwater flow (Equation (1)) based on the FD and FE formulations. The two proposed numerical models have been developed to simulate groundwater levels for two case studies: the synthetic and real-world case studies. The formulation of the numerical modeling for the synthetic case study has been first described. This follows real-world case explanations. In the synthetic case, groundwater modeling was performed in the non-steady state under three different schemes. We captured the spatio-temporal variability of the aquifer properties by implementing MEA for well 85 (calibration well; hereafter W85), and then the weights associated with the MEA were utilized to simulate groundwater drawdown for the validation wells including well 129, well 174, well 156, well 66, and well 68 (hereafter W129, W174, W156, W66, and W68). In the real-world case, we focused on the steady-state modeling and numerical outputs were obtained based on the FD and FE without any scheme. Finally, the skill rating of the numerical solutions using various MEA methods was performed through the RMSE. Figure 5 illustrates the conceptual framework of our developed aquifer model. All computational steps were carried out using a laptop system with a Quad-Core A8 AMD 1.09 GHz CPU and 8 GB RAM. The groundwater numerical models were programmed in the MATLAB programming environment along with the MEA that can be accessed from the corresponding author upon request.

Figure 5

Proposed workflow for the groundwater modeling framework applied in this research.

Figure 5

Proposed workflow for the groundwater modeling framework applied in this research.

In addition, we used RMSE (root mean square error) to assess the performance of various simulation models. Since, in groundwater modeling, water level fluctuations have a lower dynamic than surface water simulation (i.e., rainfall-runoff simulation), the variance-based criterion such as Nash–Sutcliff efficiency may not assess the ability of simulation models well. Also, the correlation-based criterion such as the coefficient of determination R2 was dismissed because these factors merely investigate the correlation and the fluctuations of simulated groundwater flow (underestimation/overestimation) could largely disregard in evaluation (Krause et al. 2005).

RESULTS AND DISCUSSION

This study developed FD and FE for groundwater simulation in synthetic and real-world case studies. We coupled four MEA approaches with FD and FE to perform the groundwater modeling in the non-steady state for the synthetic case study, while both FE and FE were tested for the steady-state modeling of groundwater-level fluctuation across the Birjand aquifer. The performance of each modeling result is explained below.

Modeling performance

This section addresses the performance of various numerical solutions. At the first, we developed three different schemes for each FD and FE model to simulate the groundwater level over a theoretical aquifer domain. Specifically, we developed six groundwater models including IS-FD, CN-FD, ES-FD, IS-FE, CN-FE, and ES-FE for a synthetic aquifer system. Groundwater drawdown due to pumping extraction wells was then calculated for each numerical model and compared with an analytical drawdown in W85. Table 1 presents modeling performance for each simulation. As noted in Table 1, aquifer simulations were almost the same across six approaches. However, ES-FE and ES-FD were slightly superior to other techniques.

Table 1

RMSE values for the synthetic simulation

Numerical schemesRMSE (cm)
ES-FE 2.08 
CN-FE 2.11 
IS-FE 2.15 
ES-FD 2.39 
CN -FD 2.41 
IS -FD 2.44 
Numerical schemesRMSE (cm)
ES-FE 2.08 
CN-FE 2.11 
IS-FE 2.15 
ES-FD 2.39 
CN -FD 2.41 
IS -FD 2.44 

Good performances are bolded.

The time series of analytical and simulated drawdown for W85 is illustrated in Figure 6. Since there was no noticeable difference among different models, only one explicit result is displayed in Figure 6. The validity of numerical models has a close agreement with other related studies that employed a similar synthetic case study (see Illangasekare & Döll 1989; Kulkarni 2015). Also, the result of this study revealed that the FE slightly outperformed the FD, although there is an underestimation in both methods, especially at the end of the simulation period. This result is well compatible with the findings of some related studies (Simpson & Clement 2003; Akbarpour et al. 2020). Simpson & Clement (2003) reported that the FE model is more robust for groundwater modeling than FD, especially for coarsely discretized problems. They also stated that the FE model avoids the common erroneous that arises in FD.

Figure 6

Comparison of analytical and simulated drawdowns.

Figure 6

Comparison of analytical and simulated drawdowns.

As expected, the groundwater table showed more stability and less drawdown in the beginning days, and thereby the numerical simulation is in good agreement with the analytical solution. After testing the proposed numerical models in the synthetic case study, we implemented the modeling setup for the real-world case study to simulate groundwater levels across the Birjand arid aquifer system.

Since most of the modeling setups for groundwater modeling such as defining the boundary conditions, assuming the initial values, subsurface and surface recharge, and applying the extraction and observation wells are all implemented in a steady state, and in non-steady state merely the temporal variation of these employed factors is applied and simulation will be continued by the same conceptual model constructed in the steady state, the groundwater modeling in the Birjand aquifer was carried out only under steady state. For this setup, no time derivative of the potential head concerning time and hence no different modeling schemes were considered (see Wang & Anderson (1995) for more information).

Table 2 represents the observed groundwater level versus the simulated values including the RMSE. As shown, the developed numerical methods provided good results compared with observed groundwater levels across various locations. In terms of RMSE, the FE prediction seems to be more accurately reflecting the suitability of this model to simulate high-dimensional groundwater dynamics. Figure 7 displays the variability of potential heads, as well as the error associated with FE simulation (the last two columns in Table 2). As shown, the maximum error was obtained while simulating the middle portion of the aquifer system. Also, based on these simulation results, the dominant groundwater direction appears to be from east to west and to southwest regions. The average groundwater drawdown is expected to be about 0.7 m annually considering aquifer depth (distance between the bedrock and the water level) and the storage capacity (840 MCM approximately). Comparing this result with long-term groundwater levels, observations discovered a severe depletion of storage capacity (−18 MCM) during the 2000–2020 period that caused the aquifer to pass the stable situation and converted to a critical state.

Table 2

Error prediction for the FD and FE models in the observed wells

PiezometerObserved potential head (m)Error prediction
FD (m)FE (m)
P13 1,264.07 − 0.63 − 0.13 
P53 1,299.10 0.43 0.12 
P85 1,291.00 − 1.92 0.06 
P212 1,296.60 − 0.76 0.23 
P393 1,392.91 0.11 0.1 
P568 1,322.76 − 5.12 0.27 
P644 1,310.08 − 1.73 0.06 
P631 1,307.29 − 0.19 − 0.18 
P760 1,363.28 0.27 0.08 
P811 1,358.05 − 0.57 − 0.51 
P995 1,342.68 0.32 0.89 
RMSE (m) 1.77 0.34 
PiezometerObserved potential head (m)Error prediction
FD (m)FE (m)
P13 1,264.07 − 0.63 − 0.13 
P53 1,299.10 0.43 0.12 
P85 1,291.00 − 1.92 0.06 
P212 1,296.60 − 0.76 0.23 
P393 1,392.91 0.11 0.1 
P568 1,322.76 − 5.12 0.27 
P644 1,310.08 − 1.73 0.06 
P631 1,307.29 − 0.19 − 0.18 
P760 1,363.28 0.27 0.08 
P811 1,358.05 − 0.57 − 0.51 
P995 1,342.68 0.32 0.89 
RMSE (m) 1.77 0.34 
Figure 7

Simulated potential heads driven by the FE model for real-world case study across the Birjand aquifer.

Figure 7

Simulated potential heads driven by the FE model for real-world case study across the Birjand aquifer.

Despite incorporating a similar conceptual model into six groundwater numerical models, different performances were achieved. This reflects the fact that addressing structural uncertainty is important and seems to be significant in this arid groundwater simulation setup. Moreover, the FE performed best for P644 and P760, while the worst prediction was obtained for P811 and P995 where boundary conditions have lower impacts and extraction activities led to high fluctuations of the groundwater level.

Ensemble method application

This section focuses on coupling various groundwater models with the MAE methods. Table 3 shows the performance results of ensemble methods for the W85 of the synthetic case test. The values of the RMSE revealed the ability of ensemble approaches to correct biases and errors in simulation. Based on this result, the WAM result was unsatisfactory, while the SMA, MMSE, and particularly the M3SE proved to be superior models. However, the M3SE generated a more accurate simulation compared with others.

Table 3

The RMSE values for different ensemble methods in the calibration well for the synthetic case study

Ensemble methodsRMSE (cm)
SMA 1.199 
WAM 2.076 
MMSE 0.11 
M3SE 0.007 
Ensemble methodsRMSE (cm)
SMA 1.199 
WAM 2.076 
MMSE 0.11 
M3SE 0.007 

Figure 8 illustrates the outputs of ensemble methods for the ES-FE and ES-FD prediction models and compares them with the analytical solution. Results indicated that the MMSE and M3SE are skillful simulators, while the WAM was inclined to the ES-FE model and was unable to mimic groundwater flow fluctuations. Overall, the SMA slightly overestimated the groundwater level compared with the ES-FE and ES-FD. These results are in good agreement with other related studies in particular Shamshirband et al. (2019) where two ensemble models were developed to combine the forecasting of chlorophyll concentration. Also, our finding is compatible with the conclusion of Ajami et al. (2006) that combined the rainfall-runoff predictions through the same MEA methods.

Figure 8

Groundwater drawdown simulation using multiple ensemble methods. The best numerical scheme versus analytical drawdown is also illustrated for a synthetic aquifer.

Figure 8

Groundwater drawdown simulation using multiple ensemble methods. The best numerical scheme versus analytical drawdown is also illustrated for a synthetic aquifer.

Table 4 describes the performance of ensemble methods for the Birjand aquifer simulation using a steady-flow simulation. According to Tables 2 and 4, the ensemble methods provided better simulation results compare to the stand-alone models. For instance, the ability of SMA (as the simplest ensemble method) appears to be more accurate than the FD. Furthermore, the MMSE (a more complex technique) provided more reliable groundwater estimates compares to the FD and FE simulations. Since the performance of the ensemble methods in the Birjand aquifer was merely examined based on the steady-flow condition, the M3SE simulation was not possible due to its non-steady-state process and was left out from the computation. Accordingly, different ensemble techniques showed that reducing the total uncertainty of groundwater modeling is useful particularly in the arid setting where the boundary condition is highly variable due to different stressors.

Table 4

Observed and ensemble groundwater-level simulations across different piezometers for the Birjand aquifer

PiezometerObserved potential head (m)SMAWAMMMSE
P13 1,264.07 1,264.09 1,263.91 1,263.87 
P53 1,299.10 1,299.78 1,299.24 1,299.19 
P85 1,291.00 1,290.47 1,290.94 1,290.92 
P212 1,296.60 1,296.74 1,296.77 1,296.74 
P393 1,392.91 1,393.42 1,393.01 1,392.96 
P568 1,322.76 1,320.74 1,322.71 1,322.73 
P644 1,310.08 1,309.65 1,310.03 1,310.01 
P631 1,307.29 1,307.51 1,307.11 1,307.07 
P760 1,363.28 1,363.86 1,363.37 1,363.32 
P811 1,358.05 1,357.91 1,357.54 1,357.49 
P995 1,342.68 1,343.69 1,343.54 1,343.50 
RMSE (cm) 78.08 32.05 31.82 
PiezometerObserved potential head (m)SMAWAMMMSE
P13 1,264.07 1,264.09 1,263.91 1,263.87 
P53 1,299.10 1,299.78 1,299.24 1,299.19 
P85 1,291.00 1,290.47 1,290.94 1,290.92 
P212 1,296.60 1,296.74 1,296.77 1,296.74 
P393 1,392.91 1,393.42 1,393.01 1,392.96 
P568 1,322.76 1,320.74 1,322.71 1,322.73 
P644 1,310.08 1,309.65 1,310.03 1,310.01 
P631 1,307.29 1,307.51 1,307.11 1,307.07 
P760 1,363.28 1,363.86 1,363.37 1,363.32 
P811 1,358.05 1,357.91 1,357.54 1,357.49 
P995 1,342.68 1,343.69 1,343.54 1,343.50 
RMSE (cm) 78.08 32.05 31.82 

To validate the efficiency of ensemble methods, an optimized setting of these models for the W85 was applied to other wells (i.e., W129, W174, W156, W66, and W68) through the synthetic case study. Table 5 shows the RMSE values of MEA and the best numerical solution in the validation period.

Table 5

The RMSE values for different ensemble and numerical methods in the validation wells for the synthetic cast study

Ensemble and best numerical methodsWell number
1291741566866
ES-FE 4.04 2.21 2.52 1.45 1.51 
ES-FD 4.30 2.55 2.82 1.89 1.91 
SMA 1.54 1.36 1.38 1.23 0.98 
WAM 4.04 2.21 2.52 1.45 1.51 
MMSE 2.33 0.66 0.87 0.46 0.66 
M3SE 0.06 0.02 0.03 0.06 0.06 
Ensemble and best numerical methodsWell number
1291741566866
ES-FE 4.04 2.21 2.52 1.45 1.51 
ES-FD 4.30 2.55 2.82 1.89 1.91 
SMA 1.54 1.36 1.38 1.23 0.98 
WAM 4.04 2.21 2.52 1.45 1.51 
MMSE 2.33 0.66 0.87 0.46 0.66 
M3SE 0.06 0.02 0.03 0.06 0.06 

Compared with ES-FE and ES-FD numerical simulation results (first two rows of Table 5), the ensemble methods generated better simulations, indicating that the estimated weights for the W85 are a well-optimized setting to validate the groundwater level in other wells. Based on Table 5, the M3SE outperformed other methods. This revealed that M3SE can reduce model uncertainty significantly. As the numerical models have a spatial dependency, the success of the ensemble method also depends on capturing the spatial variability of aquifer parameters. For example, the MEA methods generated a more precise simulation for the W68, while the outputs for the W129 appeared to be unsatisfactory.

Estimated weights results

This section discusses the contribution of each ensemble method in estimating the average weights of groundwater-level simulation. Table 6 represents the calculated weight of each numerical scheme in different ensemble settings for the synthetic case study. These results revealed that the ensemble methods employed different strategies for weight calculation. For example, the WAM maximized the proportion of the ES-FE (as the best numerical method), while the ES-FE never obtained an appropriate weight in the MMSE or the M3SE setting. A key point inferring here is that, although the accuracy of different numerical schemes was approximately similar, their weights and contribution are very significant and different among the competing models. Additionally, the M3SE's estimated weights (shown in Table 6) were more consistent and showed more uniform distribution across the simulation period compared with those estimated by the MMSE and WAM. This result is in good agreement with Raftery et al. (2003), Shamshirband et al. (2019), and Ajami et al. (2006).

Table 6

Estimated weights for each numerical scheme in the synthetic case study

Numerical methodsWeights
WAMMMSEM3SE
ES-FE 0.9999 5.27 0.6 
CN-FE 7.28 × 10−5 −17.56 −0.87 
IS-FE 5.25 × 10−6 13.61 0.82 
ES-FD 2.75 × 10−6 1.14 × 105 5.45 
CN-FD 2.70 × 10−6 − 2.28 × 105 7.52 
IS-FD 2.67 × 10−6 1.14 × 105 −12.52 
Numerical methodsWeights
WAMMMSEM3SE
ES-FE 0.9999 5.27 0.6 
CN-FE 7.28 × 10−5 −17.56 −0.87 
IS-FE 5.25 × 10−6 13.61 0.82 
ES-FD 2.75 × 10−6 1.14 × 105 5.45 
CN-FD 2.70 × 10−6 − 2.28 × 105 7.52 
IS-FD 2.67 × 10−6 1.14 × 105 −12.52 

Since in the real-world case study we only focused on SMA, WAM, and MMSE simulations based on a steady state, the weights of each model were only obtained for those models displayed in Table 7. As presented, ensemble weights for the FD model seem to be insignificant, while both WAM and MMSE contributed meaningful weights to the FE simulation. Specifically, the estimated weights of numerical methods in the real-world case study have more consistency and reliability. That is because the numerical models' simulation in the Birjand aquifer is more independence than the synthetic case which leds to more appealing MEA results.

Table 7

Estimated weights for each numerical scheme in the Birjand aquifer

Numerical methodsWAMMMSE
FD 0.0598 0.0472 
FE 0.9402 0.9528 
Numerical methodsWAMMMSE
FD 0.0598 0.0472 
FE 0.9402 0.9528 

To address the multi-collinearity problem and its impacts on the ensemble method, a new time-series data were generated and added to the simulation as a multiplier for the analytical solution. Table 8 shows the computed weights and performances of groundwater simulations using the new dataset.

Table 8

Estimated weights and RMSE values for each competing method using a new multiplier approach

Numerical methodsWeights
SMAWAMMMSEM3SE
ES-FE 0.253 4.255 0.590 
CN-FE 5.88 × 10−15 −14.410 −0.855 
IS-FE 4.80 × 10−13 11.264 0.809 
ES-FD 2.99 × 10−13 9.821 × 104 6.022 
CN-FD 2.87 × 10−13 − 1.968 × 105 6.943 
IS-FD 2.76 × 10−13 9.861 × 104 − 12.511 
Multiplier 0.746 0.151 2.607 × 10 4 
RMSE 1.09 1.81 0.11 0.007 
Numerical methodsWeights
SMAWAMMMSEM3SE
ES-FE 0.253 4.255 0.590 
CN-FE 5.88 × 10−15 −14.410 −0.855 
IS-FE 4.80 × 10−13 11.264 0.809 
ES-FD 2.99 × 10−13 9.821 × 104 6.022 
CN-FD 2.87 × 10−13 − 1.968 × 105 6.943 
IS-FD 2.76 × 10−13 9.861 × 104 − 12.511 
Multiplier 0.746 0.151 2.607 × 10 4 
RMSE 1.09 1.81 0.11 0.007 

Despite including a new multiplier with high dependency on analytical data, the MMSE and M3SE weights (0.151 and 0.00026, respectively) were not significant. This reflects the fact that multi-collinearity exists in the MEA simulation results. Besides, SMA showed unsatisfactory results and provided a constant weight equal to 1 for all numerical schemes including a new multiplier method. Furthermore, the WAM presented a more uniformity of weights compared with prior results (see Table 6). Overall, including a new dataset in the ensemble simulations proved to be beneficial since this approach improved the convergence of the models and provided better results compared with simulating each numerical model individually. These outcomes further proved that multi-collinearity problems exist in the ensemble simulation when two and more modeling parameters are highly correlated. This may make ensemble simulation more vulnerable or may result in unstable and unreasonable estimates of the weights.

CONCLUSIONS AND FUTURE WORK

This study presented an ensemble simulation of groundwater modeling across a vulnerable arid aquifer system in the east portion of Iran. To enhance groundwater numerical simulation, various ensemble techniques such as SMA, WAM, MMSE, and M3SE were tested. The FE- and FD-based schemes such as the implicit, explicit, and CN were then coupled with the ensemble methods to simulate groundwater level and behavior using synthetic as well as real-world (the Birjand aquifer system) case studies. The simulation results of the synthetic and the Birjand aquifer system were compared in terms of how the accuracy of multi-model simulations varies across different methods.

The ES in both FE and FD models (ES-FE and ES-FD) outperformed other models in a synthetic aquifer, while FE prediction proved to be superior for the real-world application. The analytical groundwater drawdown was equal to 42.5 cm, while simulated drawdown through the ES-FE and ES-FE models were calculated 39 and 38.5 cm, respectively. Overall, these simulation results revealed that the M3SE approach is more skillful to predict groundwater levels for the synthetic case study. For the real-world case study, all ensemble methods enhanced groundwater simulation, although the MMSE outperformed other methods. The good fit of simulated groundwater level to observed water levels and reasonableness of the aquifer estimated parameter values indicate that MMSE is a good approach for combining various groundwater simulation models. The model is intended to be useful for combining the simulation results in an arid groundwater system with respect to multiple demands and supplies.

Focusing on various ensemble techniques presented here, the authors recommend using M3SE and MMSE when dealing with an aquifer system with complex hydraulic and geological structures. Both MMSE and M3SE contribute the regression coefficients (weights) computed over the training period to the simulation through the unconstrained least squares technique where they can be assigned to any real numbers. This seems to be appropriate and reliable when dealing with a highly nonlinear arid groundwater system.

Here we listed our recommendations for implementing the MEA approaches in groundwater studies. First, the application of FD requires regular mesh in which some dummy nodes must be included and estimated, although it may propagate the uncertainty in simulation. Second, the FE model considerably outperformed the FD method because writing the equations in a matrix form to include related stresses, forces, or strains in the simulation seems to be more effective than replacing the derivatives with simple differences. However, arid aquifer simulation using the FE model in an irregular geometry has some disadvantages. For example, the FE simulation requires extensive computational time on preprocessing practices to create the elements. Furthermore, the modeling setup may not fully cover the area of irregular domains, and this may increase structural uncertainty and bias in simulation. A probable alternative is to use the Meshfree (Mfree) method that relaxes the geometry domain proposed by Thomas et al. (2018). Furthermore, it would be beneficial to perform MEA capabilities in more unsteady and complex groundwater simulation problems.

Third, including a multiplier to the simulation models revealed the presence of multi-collinearity in simulation that may result in the unstable estimation of the weights. To reduce the vulnerability of the models to the multi-collinearity problem, the authors recommend using the penalized regression methods such as the Least Absolute Shrinkage and Selection Operator (LASSO) that seem to provide more reliable estimated weights (Tibshirani 1996; Hammami et al. 2012). Coupling groundwater numerical solution with the Bayesian model averaging techniques might be a good practice because the ensemble weights are proportional to the individual model's accuracy and can be computed recursively (Hoeting et al. 1999). Our future work will focus on implementing these suggested approaches that will be reported in our future research works.

Groundwater simulation can be complicated with the complexities of the geological formation and demand patterns that make aquifer simulation especially difficult. In an arid environment, complex aquifer settings, and overshooting stresses, the system becomes extremely fragile and sensitive making groundwater simulation even more challenging. We stressed that our result is site-specific and our findings are necessarily subject to data limitations. Although the development of all models has been driven by the needs of the Birjand aquifer, but it is an appropriate tool that can be used in any arid aquifer. All simulations were performed by programming in a MATLAB environment with a high computational run time. Our outcomes showed promising results for multi-model application as a superior alternative to single-model simulation and may improve the knowledge on the functioning of arid aquifer systems with multiple stressors and demands.

ACKNOWLEDGEMENTS

Th authors acknowledge the constructive comments of two reviewers. The data and related code are available upon a request to the corresponding author.

CONFLICTS OF INTEREST

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

FUNDING

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

RESEARCH DATA

The data are deployed at GitHub and are available at https://github.com/AhmadJafarzadeh. Also, the codes are available from the corresponding author upon request.

AUTHORS’ CONTRIBUTIONS

A.J.: conceptualization, methodology, software, visualization, investigation, writing-original draft preparation, data curation, writing-reviewing and editing. M.P.-B.: conceptualization, methodology, software, visualization, investigation, writing-original draft preparation. A.A.: supervision, investigation, validation, writing-reviewing and editing, data curation. A.K.-S.: conceptualization, methodology, software, visualization, investigation, writing-original draft preparation. S.Z.S.: supervision, validation, methodology, writing-reviewing and editing.

DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

REFERENCES

Ajami
N. K.
,
Gupta
H.
,
Wagener
T.
&
Sorooshian
S.
2004
Calibration of a semi-distributed hydrologic model for streamflow estimation along a river system
.
Journal of Hydrology
298
(
1–4
),
112
135
.
Akbarpour
A.
,
Zeynali
M. J.
&
Tahroudi
M. N.
2020
Locating optimal position of pumping wells in aquifer using meta-heuristic algorithms and finite element method
.
Water Resources Management
34
(
1
),
21
34
.
Alvarez
P. J.
&
Illman
W. A.
2005
Bioremediation and Natural Attenuation: Process Fundamentals and Mathematical Models, vol. 27. John Wiley & Sons. Hoboken, NJ, USA
.
Arnold
J. G.
,
Allen
P. M.
&
Bernhardt
G.
1993
A comprehensive surface-groundwater flow model
.
Journal of Hydrology
142
(
1–4
),
47
69
.
Barzegar
R.
,
Moghaddam
A. A.
,
Deo
R.
,
Fijani
E.
&
Tziritis
E.
2018
Mapping groundwater contamination risk of multiple aquifers using multi-model ensemble of machine learning algorithms
.
Science of the Total Environment
621
,
697
712
.
Bates
J. M.
&
Granger
C. W.
1969
The combination of forecasts
.
Journal of the Operational Research Society
20
(
4
),
451
468
.
Berglund
J. L.
,
Toran
L.
&
Herman
E. K.
2020
Can Karst conduit models be calibrated? A dual approach using dye tracing and temperature
.
Groundwater
58
(
5
),
924
937
.
Desai
Y. M.
,
Eldho
T. I.
&
Shah
A. H.
2011
Finite Element Method with Applications in Engineering
.
Pearson Education Limited
,
Delhi
.
Dickinson
J. P.
1973
Some statistical results in the combination of forecasts
.
Journal of the Operational Research Society
24
(
2
),
253
260
.
Faust
C. R.
&
Mercer
J. W.
1980
Ground-water modeling: numerical models
.
Groundwater
18
(
4
),
395
409
.
Fletcher
S.
,
Strzepek
K.
,
Alsaati
A.
&
de Weck
O.
2019
Learning and flexibility for water supply infrastructure planning under groundwater resource uncertainty
.
Environmental Research Letters
14
(
11
),
114022
.
Georgakakos
K. P.
,
Seo
D. J.
,
Gupta
H.
,
Schaake
J.
&
Butts
M. B.
2004
Towards the characterization of streamflow simulation uncertainty through multimodel ensembles
.
Journal of Hydrology
298
(
1–4
),
222
241
.
Guiguer
N.
&
Franz
T.
1996
User's Manual for Visual MODFLOW
.
Waterloo Hydrogeologic Inc.
,
Waterloo, Ontario
,
Canadá
.
Hammami
D.
,
Lee
T. S.
,
Ouarda
T. B.
&
Lee
J.
2012
Predictor selection for downscaling GCM data with LASSO
.
Journal of Geophysical Research: Atmospheres
117
(
D17
),
1
11
.
Hamraz
B. S.
,
Akbarpour
A.
,
Bilondi
M. P.
&
Tabas
S. S.
2015
On the assessment of ground water parameter uncertainty over an arid aquifer
.
Arabian Journal of Geosciences
8
(
12
),
10759
10773
.
Hoeting
J. A.
,
Madigan
D.
,
Raftery
A. E.
&
Volinsky
C. T.
1999
Bayesian model averaging: a tutorial
.
Statistical Science
14
,
382
417
.
Illangasekare
T. H.
&
Döll
P.
1989
A discrete kernel method of characteristics model of solute transport in water table aquifers
.
Water Resources Research
25
(
5
),
857
867
.
Istok
J.
1989
Groundwater Modeling by the Finite Element Method
. American Geophysical Union, Washington, DC.
Javandel
I.
&
Witherspoon
P. A.
1968
Application of the finite element method to transient flow in porous media
.
Society of Petroleum Engineers Journal
8
(
03
),
241
252
.
Krakauer
N. Y.
,
Li
H.
&
Fan
Y.
2014
Groundwater flow across spatial scales: importance for climate modeling
.
Environmental Research Letters
9
(
3
),
034003
.
Krause
P.
,
Boyle
D. P.
&
Bäse
F.
2005
Comparison of different efficiency criteria for hydrological model assessment
.
Advances in Geosciences
5
,
89
97
.
Krishnamurti
T. N.
,
Kishtawal
C. M.
,
LaRow
T. E.
,
Bachiochi
D. R.
,
Zhang
Z.
,
Williford
C. E.
,
Gadgil
S.
&
Surendran
S.
1999
Improved weather and seasonal climate forecasts from multimodel superensemble
.
Science
285
(
5433
),
1548
1550
.
Kulkarni
N. H.
2015
Numerical simulation of groundwater recharge from an injection well
.
International Journal of Water Resources and Environmental Engineering
7
(
5
),
75
83
.
Lamontagne-Hallé
P.
,
McKenzie
J. M.
,
Kurylyk
B. L.
&
Zipper
S. C.
2018
Changing groundwater discharge dynamics in permafrost regions
.
Environmental Research Letters
13
(
8
),
084017
.
Li
J.
,
Li
F.
,
Li
H.
,
Guo
C.
&
Dong
W.
2019
Analysis of rainfall infiltration and its influence on groundwater in rain gardens
.
Environmental Science and Pollution Research
26
(
22
),
22641
22655
.
Mackay
J. D.
,
Jackson
C. R.
&
Lei
W.
2014
A lumped conceptual model to simulate groundwater level time-series
.
Environmental Modelling & Software
61
,
229
245
.
Maquin
M.
,
Mouche
E.
,
Mügler
C.
,
Pierret
M. C.
&
Viville
D.
2017
A soil column model for predicting the interaction between water table and evapotranspiration
.
Water Resources Research
53
(
7
),
5877
5898
.
McDonald
M. G.
&
Harbaugh
A. W.
1984
A Modular Three-Dimensional Finite-Difference Groundwater Flow Model (Open-File Report 83-875)
.
U.S. Department of the Interior. U.S. Geological Survey
,
Reston, VA
.
Minderhoud
P. S. J.
,
Erkens
G.
,
Pham
V. H.
,
Bui
V. T.
,
Erban
L.
,
Kooi
H.
&
Stouthamer
E.
2017
Impacts of 25 years of groundwater extraction on subsidence in the Mekong delta, Vietnam
.
Environmental Research Letters
12
(
6
),
064006
.
Mustafa
S. M. T.
,
Nossent
J.
,
Ghysels
G.
&
Huysmans
M.
2020
Integrated Bayesian multi-model approach to quantify input, parameter and conceptual model structure uncertainty in groundwater modeling
.
Environmental Modelling & Software
126
,
104654
.
Narasimhan
T. N.
&
Witherspoon
P. A.
1976
An integrated finite difference method for analyzing fluid flow in porous media
.
Water Resources Research
12
(
1
),
57
64
.
Ou
G.
,
Li
R.
,
Pun
M.
,
Osborn
C.
,
Bradley
J.
,
Schneider
J.
&
Chen
X. H.
2016
A MODFLOW package to linearize stream depletion analysis
.
Journal of Hydrology
532
,
9
15
.
Owais
S.
,
Atal
S.
,
Sreedevi
P. D.
2008
Governing Equations of Groundwater Flow and Aquifer Modelling Using Finite Difference Method
. In:
Groundwater Dynamics in Hard Rock Aquifers
(
Ahmed
S.
,
Jayakumar
R.
&
Salih
A.
, eds).
Springer
,
Dordrecht
.
https://doi.org/10.1007/978-1-4020-6540-8_16
.
Pacheco
F. A. L.
,
Martins
L. M. O.
,
Quininha
M.
,
Oliveira
A. S.
&
Fernandes
L. S.
2018
Modification to the DRASTIC framework to assess groundwater contaminant risk in rural mountainous catchments
.
Journal of Hydrology
566
,
175
191
.
Pinder
G. f.
&
Gray
W. G.
1977
Finite Element Simulation in Surface and Subsurface Hydrology
.
Academic Press
,
New York, NY
.
Raftery
A. E.
,
Balabdaoui
F.
,
Gneiting
T.
&
Polakowski
M.
2003
Using Bayesian Model Averaging to Calibrate Forecast Ensembles. Technical Report No. 440
.
Department of Statistics, University of Washington
.
Richtmyer
R. D.
&
Morton
K. W.
1994
Difference Methods for Initial-Value Problems
. Krieger Publisher Company. Malabar, FL.
Sadeghi-Tabas
S.
,
Samadi
S. Z.
,
Akbarpour
A.
&
Pourreza-Bilondi
M.
2017
Sustainable groundwater modeling using single-and multi-objective optimization algorithms
.
Journal of Hydroinformatics
19
(
1
),
97
114
.
Shamseldin
A. Y.
,
O'Connor
K. M.
&
Liang
G. C.
1997
Methods for combining the outputs of different rainfall–runoff models
.
Journal of Hydrology
197
(
1–4
),
203
229
.
Shamshirband
S.
,
Jafari Nodoushan
E.
,
Adolf
J. E.
,
Abdul Manaf
A.
,
Mosavi
A.
&
Chau
K. W.
2019
Ensemble models with uncertainty analysis for multi-day ahead forecasting of chlorophyll a concentration in coastal waters
.
Engineering Applications of Computational Fluid Mechanics
13
(
1
),
91
101
.
Shepley
M. G.
,
Whiteman
M. I.
,
Hulme
P. J.
&
Grout
M. W.
2012
Introduction: groundwater resources modelling: a case study from the UK
.
Geological Society, London, Special Publications
364
(
1
),
1
6
.
Simpson
M. J.
&
Clement
T. P.
2003
Comparison of finite difference and finite element solutions to the variably saturated flow equation
.
Journal of Hydrology
270
(
1–2
),
49
64
.
Strack
O. D. L.
2017
Finite differences and finite elements
. In:
Analytical Groundwater Mechanics
.
Cambridge University Press
,
Cambridge
, pp.
403
422
.
http://doi.org/10.1017/9781316563144.011
.
Thomas
A.
,
Majumdar
P.
,
Eldho
T. I.
&
Rastogi
A. K.
2018
Simulation optimization model for aquifer parameter estimation using coupled meshfree point collocation method and cat swarm optimization
.
Engineering Analysis with Boundary Elements
91
,
60
72
.
Thompson
P. D.
1977
How to improve accuracy by combining independent forecasts
.
Monthly Weather Review
105
(
2
),
228
229
.
Tibshirani
R.
1996
Regression shrinkage and selection via the lasso
.
Journal of the Royal Statistical Society: Series B (Methodological)
58
(
1
),
267
288
.
Vrugt
J. A.
,
Diks
C. G.
,
Gupta
H. V.
,
Bouten
W.
&
Verstraten
J. M.
2005
Improved treatment of uncertainty in hydrologic modeling: combining the strengths of global optimization and data assimilation
.
Water Resources Research
41
(
1
),
1
17
.
Wang
F.
&
Anderson
M.
1995
Introduction to Groundwater Modeling Finite Difference and Finite Element
.
Academic Press
,
New York
.
Xu
T.
,
Valocchi
A. J.
,
Ye
M.
,
Liang
F.
&
Lin
Y. F.
2017
Bayesian calibration of groundwater models with input data uncertainty
.
Water Resources Research
53
(
4
),
3224
3245
.
Zhou
Z.-H.
2012
Ensemble Methods: Foundations and Algorithms
.
Chapman and Hall/CRC
,
London
.
Zienkiewicz
O.
,
Mayer
P.
&
Cheung
Y. K.
1966
Solution of anisotropic seepage by finite elements
.
Journal of the Engineering Mechanics Division
92
(
1
),
111
120
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).