Automatic calibration of the groundwater simulation model with high parameter dimensionality using sequential uncertainty fitting approach

Groundwater plays an imperative role in supplying water for agriculture, industry, and drinking needs. Meanwhile, modeling groundwater has substantial effects on the economy, notably on water resource management (Dagès et al. 2012). Recent variations in the hydrologic cycle due to changes in global climate, population growth, and agricultural area expansion have led to the overuse of groundwater resources in most parts of the world. Accordingly, the management of these resources is challenging due to problems caused by overuse, including reduced groundwater levels, diminished water quality, and increased harvest costs (Wada et al. 2010).

Undoubtedly, numerical models of groundwater are essential tools for groundwater resource management. In recent years, different researchers have applied numerical models to manage groundwater resources (e.g., Yihdego et al. 2016; Kamali & Niksokhan 2017). The main challenge researchers face in utilizing numerical models is the calibration of these models. In general, these models can be calibrated by trial and error and automatically. In recent decades, studies have proposed various methods based on inverse modeling for the automatic estimation of hydrological parameters (e.g., Dai & Samper 2006; Shang et al. 2016). Automated calibration methods are far faster than trial-and-error calibration and offer better outcomes due to their capability of searching more parameter space (Droogers et al. 2001). Further, different software applications have been developed for the automatic calibration of groundwater models, including PEST (Doherty et al. 1994), UCODE (Poeter & Hill 1998), and HydroPSO (Zambrano-Bigiarini & Rojas 2013). Although these applications accelerate the calibration process, they may obtain irrational values for parameters since they mainly aim to match observational and computational values regardless of physical realities (Delottier et al. 2017).

The findings of the applied models for groundwater modeling heavily rely on the input parameters of the model. On the other hand, the scarcity of field data, poor knowledge of the conceptual model, various modeling scenarios, and hydrogeological complexities cause uncertainties in the model (Wu & Zeng 2013). Thus, automated calibration methods based on uncertainty analysis have been of great interest to researchers in recent years. Interval estimates of parameter values are provided in these methods, unlike other optimization methods that generate a point and definite parameter estimates. Hydrologists have utilized uncertainty-based methods for several reasons such as errors in measuring observations (e.g., temperature and precipitation), the uncertainty of the output variables (e.g., the discharge and sedimentation rate), along with the applied model assumptions (Abbaspour et al. 2007) and the existence of different responses with identical values of the objective function. Uncertainty analysis, based on input variable parameters, allows the description of all possible outputs, along with their probability of occurrence. In the uncertainty-based automatic calibration method, a range of values is assumed for each parameter, and then parameter intervals are detailed during the calculation process. Currently, the most commonly used methods for parameter uncertainty testing are the general dissimilarity method of Generalized Likelihood Uncertainty Estimation (GLUE) by Hassan et al. (2008), and Markov Chain Monte Carlo (MCMC) by Hassan et al. (2009). In addition, other methods include similar algorithms such as Blocking Monte Carlo Markov Chain (Fu & Gómez-Hernández 2009), Null Space Monte Carlo Markov Chain (Sepúlveda & Doherty 2015), Bisen structure (Li & Tsai 2009), parameter solution (ParaSol) by van Griensven & Meixner (2006), and sequential uncertainty fitting (SUFI-2) by Abbaspour et al. (2004). Masoumi et al. (2019) developed the SUFI-2 algorithm for groundwater model calibration. They used hydraulic conductivity and specific yield in the calibration process. They showed that SUFI-2 has excellent performance in calibrating the MODFLOW model compared with PEST methods. Based on the groundwater flow equation, although hydraulic conductivity and specific yield are essential parameters, many other input parameters can be considered as uncertain parameters in calibration.

In one comprehensive study, Yang et al. (2008) compared five different methods (including GLUE, ParaSol, SUFI-2, and Importance Sampling) for the uncertainty analysis of the Soil and Water Assessment Tool (SWAT) model in China. The advantages and disadvantages of these models are summarized in Table 1. Based on this study, uncertainty methods are divided into three categories. The first one includes methods based on informal Bayesian approaches such as GLUE and SUFI-2. These methods are simple in terms of implementation and statistical approaches. On the other hand, all model uncertainties can be presented by the uncertainty of the parameters. The second category is the formal Bayesian-based approach, such as MCMC. The effect of inputs and the model structure is considered by adding an error value to the output. It is noteworthy that the time correlation typically gives the remaining values to the output. The last category includes approaches that consider all uncertainty sources (the input and model structure) by the improved likelihood function (Yang et al. 2008). Although the models of the last two categories can account for other sources of uncertainty, their implementation can be more complicated, and thus other methods may be considered, including GLUE or SUFI-2.

SUFI-2 is one of the methods that perform calibration and model uncertainty analysis simultaneously. This method has been widely used in the automatic calibration and uncertainty determination of hydrological models, although it has not been applied for calibration and analysis of uncertainty in groundwater resources. For example, Mousavi et al. (2012) utilized this algorithm for surface water resource systems. In that study, the SUFI-2 algorithm coupled with the HEC-HMS model was used for the automatic calibration and uncertainty analysis of the Tamer River Basin in northern Iran. The findings demonstrated the excellent performance of the SUFI-2 algorithm in automatic calibration based on uncertainty. Further, Cao et al. (2018) used the SUFI-2 algorithm to analyze the sensitivity and uncertainty of a daily-scale model for a river basin and reported the acceptable performance of the calibrated SWAT model in daily flow forecasting. Similarly, Elçi (2017) applied the SUFI-2 algorithm in the SWAT-CUP software to evaluate the nutrient retention in vegetated filter strips and used it for sensitivity analysis and model calibration. Based on the findings, the model adequately simulated hydrologic flow components and nutrient transport processes in the studied stream basin. Moreover, He et al. (2019) utilized SWAT software to analyze nitrogen load in a catchment with paddy field in China and the uncertainty of simulation using the SUFI-2 algorithm. The SUFI-2 algorithm was further used to calibrate product development models. Kamali et al. (2018) developed the SUFI-2 algorithm with product modeling software (EPIC) to calibrate corn yields. It eventually resulted in user-friendly software (EPIC⁺) for calibrating product development models.

In the current research, an uncertainty-based structure was applied to calibrate the groundwater model automatically. In addition, the MODFLOW code was employed for groundwater modeling and the SUFI-2 algorithm for automatic calibration, along with the uncertainty analysis of hydraulic conductivity and specific yield of the aquifer, hydraulic conductance of river–aquifer interaction, boundary conductance, and the recharge rate. Further, the Groundwater Modeling System (GMS) user-friendly software was utilized for the pre-processing and post-processing of the MODFLOW model. The MODFLOW-SUFI-2 linked model was developed in the MATLAB software environment as well. Furthermore, the groundwater model was calibrated using particle swarm optimization (PSO), and then the findings were compared with the SUFI-2 algorithm. Although various optimization algorithms have been developed for calibrating groundwater models, the SUFI-2 algorithm has not been so far used for the calibration of groundwater models. One of the advantages of this algorithm over other conventional optimization algorithms is the time to consider parameters at the end of the calibration process and the need for fewer simulations. Accordingly, the present study sought to evaluate the effect of the proposed model, and the Ardabil plain aquifer was selected as a case study.

According to Figure 1, the Ghara-Su river basin is located in the northwest of Iran in Ardabil province (37°46′–38° 36′ N and 47°46′–48°42′ E) with an area of about 4,193 km². The Ardabil plain aquifer covers an area of 1,074 km² with a precious groundwater resource. Moreover, it has permanently been an excellent resource for supplying drinking and agricultural water over the last half-century. The topography of the surrounding area is mountainous with Sabalan mountain located in the western part of the plain at 4,850 m above mean sea level (AMSL), the Bozqush mountains in the southern region at 3,306 m AMSL, and the Talesh mountains in the eastern and southeastern parts at about 3,197 m AMSL. Additionally, the lowest height is located in the northwest of the plain at 1,170 m AMSL. Similarly, this aquifer is connected to the basin of the Ghara-Su river, which arises from surrounding mountains, generally flows from the east to the west in the Ardabil aquifer, and finally, in the northern part of Ardabil province, joins the Araz river. Balikhli-Chay, Ghori-Chay, Noran-Chay, and Saghezchi-Chay are the other rivers of this area that flow into the Ghara-Su river. The mean annual precipitation and temperature in the Ardabil plain are 273.7 mm and 9.7 °C, respectively (Nadiri et al. 2017; Taie Semiromi & Koch 2019). All data used in this study are provided by Ardabil Regional Water Authority (2013), Kord (2014), and Najjar-Ghabel et al. (2019).

Considering the diversity of layers and changes in their types in different regions, the Ardabil plain aquifer is considered an unconfined aquifer since there is a complete connection of its various parts across the plain, where only the upper boundary is enclosed in some areas. The thickness of the aquifer deposits in the Ardabil plain aquifer varies from 10 to 200 m. In the study area, there are 1,500 operating wells, two springs, and 13 qanats. The groundwater level has decreased by about 12 m due to the high withdrawal rate of groundwater in the past 25 years.

As shown in Figure 2, the groundwater level of September 2007, which was obtained from observation wells distributed in the model domain, is considered as the initial groundwater level. The input, output, and non-flow boundaries were determined according to the equipotential map of this month. Based on this figure, the output boundary is only located in the northwestern part of the aquifer, and the remaining boundaries are considered approximately as input boundaries. In addition, observation wells were used for calibrating the model. Further, the transmissivity range for the study domain, which was obtained from the pumping test, was 50–2,200 ⁠. Based on the field data, the specific yield range was 0.02–0.2. To simulate the real conditions of the aquifer, the boundary of the aquifer was divided into 22 segments with different boundary conductances.

In this study, the MODFLOW code was used to simulate the groundwater flow. This code was developed by the United States Geological Survey. The original version of this model was first prepared by McDonald & Harbaugh (1988). This code uses the finite difference method for modeling the groundwater flow. The GMS software was also applied as a graphical interface for groundwater pre- and post-processing. The GMS software provides the inputs of the model in an appropriate format by using different modules. Then, the main core of the model (which is MODFLOW in this study) is called for executing the program. Finally, the output of the model is arranged in the appropriate format for the users.

The GMS stores files in Hierarchical Data Format (HDF) when executing the MODFLOW code. Furthermore, this format stores data in a binary format and is normally used to manage and classify large data. This extension was developed by the National Center for Supercomputing Applications. In this study, MATLAB software was employed to access HDF files and link the SUFI-2 algorithm to the MODFLOW code.

MODFLOW uses different packages to solve the groundwater-flow equation. The present study utilized Well (WEL), River (RIV), Recharge (RCH), and General Head Boundary (GHB) packages as the stress package for the Ardabil aquifer. The MODFLOW–Newton Solution Method (NWT) is applied for solving the groundwater-flow equation. In addition, MODFLOW-NWT, as an updated version of MODFLOW-2005, is a program that uses the NWT solution method and Upstream-Weighting package for calculating intercell conductance. The MODFLOW-NWT version is used to solve the nonlinearity problem of the groundwater equation in the unconfined aquifer, which is caused by drying–wetting of cells (Niswonger et al. 2011).

In Equations (7)–(9), is the computed head for observation well i obtained from the MODFLOW model. Furthermore, and n denote the observation head for observation well i and the number of observation heads, respectively.

PSO is a technique in which this optimization method has been enlarged in applications in recent years since its inception by Eberhart & Kennedy (1995). This method considers the regulation of solutions, which is identified by its position and velocity vectors. The best solution of each particle is banked as experience in memory, and the best solution of all particles is attained as the best global particles.

The general applied framework of this study is shown in Figure 4. According to this figure, a simulation–calibration model was used for the groundwater calibration model. For this purpose, the MODFLOW code was utilized to simulate the groundwater flow and determine the values of the objective function for different values of parameters. In addition, the MODFLOW model was developed with writing and reading access in MATLAB software. Further, the SUFI-2 algorithm was applied for the automatic calibration and uncertainty analysis of the simulation model. In this method, by changing the values of parameters (used as the input to the groundwater simulation model), the SUFI-2 algorithm repeatedly calls the groundwater simulation model and sampling depending on the objective function of the evaluation. Then, the uncertainty criteria were calculated, and finally, optimal values were determined for each parameter.

Moreover, GMS software was used for the pre-processing and post-processing of the MODFLOW code. This software saves input files in the HDF5 format, which can be edited (write and read) in MATLAB. Furthermore, MODFLOW can be executed in MATLAB by the MODFLOW application file. Finally, the hydraulic head of each cell is achievable by the binary head file that is produced by MODFLOW. In this research, this approach was employed to present the uncertainty-based automatic calibration of the groundwater model by linking MODFLOW and the SUFI-2 algorithm.

The Ardabil MODFLOW model was developed for 12 monthly time steps (from October 2007 to September 2008) under transient conditions. The study domain was discretized by 109 columns and 95 rows with 350 × 350 (⁠⁠) cell size in one layer. As mentioned earlier, Well (WEL), River (RIV), Recharge (RCH), and General Head Boundary (GHB) packages were used in the MODFLOW model. Additionally, field data and interpolation were utilized for the entire study area. Then, the Inverse Distance Weighting method was employed for the interpolation of field data, which was previously provided by Najjar-Ghabel et al. (2019). Based on interpolation, hydraulic conductivity and specific yield were divided into four and five zones, respectively. Figure 5(a) and 5(b) display the zonation of the hydraulic conductivity and specific yield, respectively.

The river–aquifer interaction is considered based on Equations (3) and (4), and is calculated by the river hydraulic conductance and river water elevation coefficients. In this study, an inverse solution was applied to determine the initial river hydraulic conductance. In addition, the river–aquifer water exchange was estimated, and hydrometric stations were used to determine the river water elevation. Table 2 presents the initial range of these parameters. As illustrated in Figure 1, there are five river flows in the aquifer domain. The hydraulic conductance of Ghara-su (⁠⁠), Balikhli-Chay (⁠⁠)), Ghori-Chay (⁠⁠), Noran-Chay (⁠⁠), and Saghezchi-Chay (⁠⁠) were considered as uncertain parameters adjusted during the calibration of the model.

The recharge rate in Equation (2) was calculated based on effective rainfall and returned water from agriculture, industry, and municipal sewage. The amount of agricultural, industrial, and municipal return water to the aquifer relies on the water allocated to these sectors. Considering the available allocated data of the Ardabil area, 70% of the allocated water for domestic and industrial use and 25% of the allocated water to agriculture were assumed as the returned water (Kord 2014). Further, a considerable part of the rainfall in the area infiltrates and feeds on the alluvial aquifer due to the low ground surface slope of the plains and alluvial grading features. Accordingly, 15% of the rainfall was also considered for aquifer infiltration. It should be noted that the spatial variation of the recharge rate was constant due to the similarity of the soil material on the top layer and the small variation of rainfall in different locations of the study area. Thus, the amount of rainfall varied at different times in the study area.

Based on the conceptual model of the Ardabil aquifer (Figure 2), the inflow boundary of the aquifer was divided into 22 segments. The initial value of the boundary conductance for each segment was obtained from the field data interpolation of hydraulic conductivity. The initial range for segments is provided in Table 2. These parameters were adjusted during the calibration.

In Equation (12), and represent the average d-factor and p-factor of each iteration. Moreover, and are considered 0 and 1, respectively. Eventually, and are the maximum d-factor in each iteration and 0, respectively.

In Figure 7, the unsteady conditions of the 95% confidence interval for different observation wells (i.e., Obs1, Obs2, Ob19, and Obs25), along with observational data from the Ardabil aquifer are illustrated as a time series for the duration of the modeling period. As these figures show, it seems that the SUFI-2 algorithm could correctly predict the range of calibration values. In most cases, the observed values lie within the range of values predicted by the MODFLOW-SUFI-2 model.

The d-factor is calculated to have a relatively high value due to the existence of uncertainty in the structure of the model. The values for the d-factor and p-factor are also reported in Table 3. If we do not have a good quality of data, the p-factor may be low, and the result may be satisfied when 50% of data lie within 95PPU (Abbaspour et al. 2007). According to this table, most observation points were bracketed by 95PPU (on average, 60%). In some observation wells (Obs7 and Obs22), more geological information is required to reduce uncertainty. Additionally, many factors, such as measurement error, uncertainty in the structure of the conceptual model, and the unknown number of illegal wells, can affect these values. Thus, the calculated values for the uncertainty index can be acceptable due to the lack of observational data.

Appendix A displays the groundwater level of the Ardabil plain aquifer simulated and calibrated by the MODFLOW-SUFI-2 linked model. Based on this figure, the overall flow direction in the Ardabil plain aquifer is to the northwest of the aquifer, which is known as the outflow boundary in this study. The obtained flow pattern from this figure is highly similar to the pattern expected from the conceptual model. In addition, the error bar in this figure is shown in green and yellow. Errors less than one metre are represented by green color, while those between one metre and two metres are displayed in yellow. Most of the observation points show less than 1 m error. There are 35 observation wells in the study area, although only 11 of them have complete time-series data. The location of all observation wells is depicted in Figure 2.

Table 4 provides the best values for parameters after the calibration of the Ardabil groundwater model by the SUFI-2 algorithm. In the present study, the SUFI-2 algorithm was run for 500 rounds in each iteration and calibrated 48 parameters. Given the field data used for determining the initial range of parameters, the optimum value for parameters was obtained in iteration 3.

In this paper, the PSO algorithm was used for calibrating the groundwater simulation model. In this algorithm, RMSE was used as the objective function. To use this method, an initial range of parameters is required with these values (Table 2). Table 4 summarizes the final values for parameters after the calibration of the groundwater model with the above-mentioned method. As shown, the best value of parameters seems irrational since this method only aims at adjusting the computed head value with the measuring head of the observation well. The RMSE of 0.8 m was obtained based on the results obtained from model implementation. Thus, the convergence rate of the answer improved due to applying the field data for the initial value of parameters. It should be noted that the best value for the objective function was obtained in the Number of Function Evaluation (NFE) 220 (Appendix B) because of a negligible specific termination condition after 200 NFE.

The comparison of the PSO algorithm and the SUFI-2 method revealed that although the PSO method led to good coverage and a low value of RMSE, the optimum values obtained for some parameters were irrational and the maximum or minimum of the initial range was selected in most cases. On the other hand, the SUFI-2 method considers a reliable value for parameters with acceptable RMSE. For example, the optimum values for HH4 and Sy3 are obtained as 19.54 m/day and 0.14 by the SUFI-2 algorithm, respectively, which is closer to the observed data (10.5 m/day and 0.13) in comparison with the PSO algorithm. In other words, this method can better predict parameter values compared with the PSO algorithm because of using different criteria such as d-factor and p-factor in the SUFI-2 algorithm.

Figure 8 shows the simulated hydraulic head in Obs19 by the PSO and SUFI-2 approach with observational data. As depicted, these approaches have similar behavior, although the SUFI-2 method can better simulate the fluctuation of the groundwater head compared with the PSO algorithm. Table 5 presents the performance of these methods evaluated by other statistical criteria (i.e., NSE and PBIAS). According to this table, both RMSE and SUFI-2 methods demonstrate similar performance in NSE and PBIAS criteria. The NSE coefficient was calculated as 0.75 and 0.67 for PSO and SUFI-2, respectively, indicating the good capability of these methods in the groundwater simulation model. Moreover, PBIAS coefficients for the methods are satisfactory (less than 1% in each method). Although the statistical criteria are identical for the two algorithms, the SUFI-2 algorithm represents a high capability in terms of the coverage rate because of the fundamental theoretical and sampling method embedded in this approach.

In general, an uncertainty-based automatic calibration of the groundwater model was developed in the present study. The MODFLOW model was used for groundwater modeling, and then the sequential uncertainty fitting (SUFI-2) algorithm was employed for the calibration and uncertainty analysis of the model. Based on the findings, the MODFLOW-SUFI-2 linked model provided an efficient tool for calibrating the Ardabil groundwater model. In addition, the performance of this method was evaluated by some statistical criteria (i.e., RMSE, NSE, and PBIS), which showed good accuracy in simulating the groundwater hydraulic head (0.86, 0.67, and −0.95, respectively). Furthermore, the result of the proposed method was compared with those of the developed MODFLOW–particle swarm optimization (PSO) model for the auto-calibration of the groundwater model. The findings showed that the SUFI-2 method could more reliably predict parameters compared with the PSO algorithm due to the use of different criteria (d-factor and p-factor). Moreover, although the statistical criteria were similar in these methods, SUFI-2 indicated a small number of iterations in comparison with PSO. The findings further revealed that the d-factor was calculated to be somewhat higher. However, the findings seem acceptable by considering a large magnitude of uncertainty and error in the structure of the model due to the size, complexity of the aquifer, and the lack of field data. Finally, this MODFLOW-SUFI-2 was generally efficient in calibrating the Ardabil groundwater model compared with the PSO algorithm.

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