Application of multi-model ensemble averaging techniques for groundwater simulation: synthetic and real-world case studies

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.

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.

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.

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.

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 m²/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.

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.

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.

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.

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.

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

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.

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 R² 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).

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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.

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.

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.

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.

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

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.

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

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.

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.

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