## Abstract

Calibration of the groundwater simulation model is one of the main challenges in the modeling process. In addition, hydrogeological complexities and the lack of field data in terms of time and space lead to uncertainty in the model. Therefore, the present study linked the groundwater simulation model (MODFLOW) and sequential uncertainty fitting approach (SUFI-2) to the uncertainty-based automatic calibration of the Ardabil groundwater model located in northwestern Iran. Hydraulic conductivity, specific yield, recharge rate, the hydraulic conductivity of the riverbed material, and the boundary conductance of the aquifer were considered as the uncertain parameters. Furthermore, the Newton solution method for the unconfined aquifer was used for solving the groundwater flow equation. A Normalized Total Uncertainty Index was defined to evaluate the performance of the SUFI-2 algorithm. According to the MODFLOW-SUFI-2 calibration findings, 60% of observational data was bracketed by a 95% confidence interval, on average. The Ardabil groundwater model was also calibrated with the PSO algorithm. In comparison with SUFI-2, although this method resulted in good coverage of the solution, it obtained irrational values for most parameters since they only aimed to match observational and computational values. Eventually, SUFI-2 showed a small number of simulation runs compared with the PSO algorithm.

## HIGHLIGHTS

Uncertainty-based auto-calibration of the groundwater simulation model.

MODFLOW-SUFI-2 linked model.

High parameter dimensionality in the calibration process.

A small number of simulation runs compared with other conventional algorithms.

## INTRODUCTION

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.

Method . | Theoretical basis . | Sampling method . | Uncertainty described by parameter uncertainty . | Difficulty of implementation . |
---|---|---|---|---|

GLUE | Normalization of generalized likelihood measure | Primitive random | All sources of uncertainty | Very easy |

ParaSol | Least squares (probability theory) | SCE-UA-based | Parameter uncertainty only | Easy |

SUFI-2 | Generalized objective function | Latin hypercube | All sources of uncertainty | Easy |

MCMC | Likelihood function (Probability theory) | Optimal parameter set based on SCE-UA | Parameter uncertainty only | More complicated |

IS | Likelihood function (Probability theory) | Primitive random | Parameter uncertainty only | More complicated |

Method . | Theoretical basis . | Sampling method . | Uncertainty described by parameter uncertainty . | Difficulty of implementation . |
---|---|---|---|---|

GLUE | Normalization of generalized likelihood measure | Primitive random | All sources of uncertainty | Very easy |

ParaSol | Least squares (probability theory) | SCE-UA-based | Parameter uncertainty only | Easy |

SUFI-2 | Generalized objective function | Latin hypercube | All sources of uncertainty | Easy |

MCMC | Likelihood function (Probability theory) | Optimal parameter set based on SCE-UA | Parameter uncertainty only | More complicated |

IS | Likelihood function (Probability theory) | Primitive random | Parameter uncertainty only | More complicated |

*Source:*Yang *et al*. (2008).

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.

## STUDY AREA

### Location and physical details

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^{2}. The Ardabil plain aquifer covers an area of 1,074 km^{2} 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).

### Hydrogeology

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.

## METHOD

### Groundwater simulation model

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.

*x*,

*y*, and

*z*directions, respectively. Moreover,

*h*and

*W*are the hydraulic head, and the recharge (

*W*is negative in this case) or discharge (

*W*is positive in this case) terms, respectively. Additionally, stands for the specific storage coefficient with , which is the specific yield coefficient;

*b*and

*t*denote the aquifer thicknesses and the time, respectively. In general, , , and denote functions of space, while

*W*can be a function of space and time (Todd & Mays 2005).

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

*h*and denote the groundwater head and the water level in the river (), respectively. Moreover, is the bottom elevation of the river () and represents the river hydraulic conductance of the river–aquifer interaction (), which is defined by Equation (5):where , , , and represent the hydraulic conductivity of the riverbed material, the length of the river, the river width, and the thickness of the riverbed layer, respectively.

*n*) per unit time ():where , , and are the boundary conductance (), the head assigned to the boundary condition, and the head in the cell

*n*, respectively.

### SUFI-2 algorithm

*et al.*(1997). Based on the Bayesian framework, this algorithm finds the best range of parameters and analyzes parameter uncertainty simultaneously (Abbaspour

*et al.*2004). As shown in Figure 3, this approach uses different steps in each round in order to fit observational data with the simulated value in the minimum number of iterations. First, the objective function and initial range of parameters are defined in this regard. In this study, the Root-Mean-Square Error (RMSE) function was used as the objective function. In the second step, sampling, which is provided by Latin Hypercube Sampling (LHS) (McKay

*et al.*1979) with uniform distributions, is evaluated by the objective function. Finally, this algorithm analyzes the uncertainty by p-factor and d-factor indices. The number of observational data, which are bracketed by the 95% prediction uncertainty (95PPU) band to total observations, is known as the p-factor (Abbaspour

*et al.*2004). Additionally, the average bandwidth of 95PPU to the standard deviation of the observational data is defined as the d-factor index. Further detail about the basic formulations of the SUFI-2 algorithm is introduced in the Appendix.

### Particle swarm optimization (PSO) algorithm

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.

*W*, , and denote inertia weight, the personal learning coefficient, and the social learning coefficient, respectively. Moreover, and indicate uniform random numbers in and reveals the best position of the swarm of the dimension. Finally, is the position of the global best of all swarms.

### MODFLOW-SUFI-2 model

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.

## RESULTS AND DISCUSSION

### The MODFLOW model setup

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.

Number . | Parameter . | First range . | Number . | Parameter . | First range . | ||
---|---|---|---|---|---|---|---|

1 | Hydrodynamic parameters | HK1 | 2–25 | 25 | 300–2,000 | ||

2 | HK2 | 2–25 | 26 | 50–500 | |||

3 | HK3 | 2–25 | 27 | The parameter corresponding to the inlet boundaries of groundwater flow | 100–200 | ||

4 | HK4 | 2–25 | 28 | 200–400 | |||

5 | Sy1 | 0.01–0.22 | 29 | 300–500 | |||

6 | Sy2 | 0.01–0.22 | 30 | 300–550 | |||

7 | Sy3 | 0.01–0.22 | 31 | 300–600 | |||

8 | Sy4 | 0.01–0.22 | 32 | 300–500 | |||

9 | Sy5 | 0.01–0.22 | 33 | 350–650 | |||

10 | Recharge parameter (due to precipitation and return of water from agriculture, drinking and industry) | 0.0001–0.0009 | 34 | 300–550 | |||

11 | 0.0001–0.0009 | 35 | 100–200 | ||||

12 | 0.0001–0.0009 | 36 | 500–1,000 | ||||

13 | 0.0001–0.0009 | 37 | 300–500 | ||||

14 | 0.0001–0.0009 | 38 | 300–550 | ||||

15 | 0.0001–0.0009 | 39 | 400–700 | ||||

16 | 0.0001–0.0009 | 40 | 400–700 | ||||

17 | 0.0001–0.0009 | 41 | 100–200 | ||||

18 | 0.0001–0.0009 | 42 | 150–300 | ||||

19 | 0.0001–0.0009 | 43 | 400–800 | ||||

20 | 0.0001–0.0009 | 44 | 400–800 | ||||

21 | 0.0001–0.0009 | 45 | 200–300 | ||||

22 | River parameter | 50–500 | 46 | 100–200 | |||

23 | 300–2,000 | 47 | 80–200 | ||||

24 | 300–2,000 | 48 | 50–100 |

Number . | Parameter . | First range . | Number . | Parameter . | First range . | ||
---|---|---|---|---|---|---|---|

1 | Hydrodynamic parameters | HK1 | 2–25 | 25 | 300–2,000 | ||

2 | HK2 | 2–25 | 26 | 50–500 | |||

3 | HK3 | 2–25 | 27 | The parameter corresponding to the inlet boundaries of groundwater flow | 100–200 | ||

4 | HK4 | 2–25 | 28 | 200–400 | |||

5 | Sy1 | 0.01–0.22 | 29 | 300–500 | |||

6 | Sy2 | 0.01–0.22 | 30 | 300–550 | |||

7 | Sy3 | 0.01–0.22 | 31 | 300–600 | |||

8 | Sy4 | 0.01–0.22 | 32 | 300–500 | |||

9 | Sy5 | 0.01–0.22 | 33 | 350–650 | |||

10 | Recharge parameter (due to precipitation and return of water from agriculture, drinking and industry) | 0.0001–0.0009 | 34 | 300–550 | |||

11 | 0.0001–0.0009 | 35 | 100–200 | ||||

12 | 0.0001–0.0009 | 36 | 500–1,000 | ||||

13 | 0.0001–0.0009 | 37 | 300–500 | ||||

14 | 0.0001–0.0009 | 38 | 300–550 | ||||

15 | 0.0001–0.0009 | 39 | 400–700 | ||||

16 | 0.0001–0.0009 | 40 | 400–700 | ||||

17 | 0.0001–0.0009 | 41 | 100–200 | ||||

18 | 0.0001–0.0009 | 42 | 150–300 | ||||

19 | 0.0001–0.0009 | 43 | 400–800 | ||||

20 | 0.0001–0.0009 | 44 | 400–800 | ||||

21 | 0.0001–0.0009 | 45 | 200–300 | ||||

22 | River parameter | 50–500 | 46 | 100–200 | |||

23 | 300–2,000 | 47 | 80–200 | ||||

24 | 300–2,000 | 48 | 50–100 |

*Note*. HK: hydraulic conductivity; Sy: specific yield; *I*: recharge rate; *C*: hydraulic conductivity of the riverbed material; *T*: boundary conductance of the aquifer.

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.

### Calibrating the MODFLOW model

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.

Name of observation well . | d-factor . | p-factor . |
---|---|---|

Obs1 | 2.80 | 66% |

Obs2 | 2.18 | 83% |

Obs3 | 2.21 | 83% |

Obs4 | 3.20 | 66% |

Obs7 | 1.71 | 25% |

Obs19 | 1.57 | 50% |

Obs22 | 1.42 | 33% |

Obs25 | 1.26 | 50% |

Obs27 | 1.20 | 50% |

Obs29 | 1.10 | 66% |

Obs33 | 2.89 | 83% |

Name of observation well . | d-factor . | p-factor . |
---|---|---|

Obs1 | 2.80 | 66% |

Obs2 | 2.18 | 83% |

Obs3 | 2.21 | 83% |

Obs4 | 3.20 | 66% |

Obs7 | 1.71 | 25% |

Obs19 | 1.57 | 50% |

Obs22 | 1.42 | 33% |

Obs25 | 1.26 | 50% |

Obs27 | 1.20 | 50% |

Obs29 | 1.10 | 66% |

Obs33 | 2.89 | 83% |

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.

Number . | Parameter . | Best value (SUFI-2) . | Best value (PSO) . | Number . | Parameter . | Best value (SUFI-2) . | Best value (PSO) . | ||
---|---|---|---|---|---|---|---|---|---|

1 | Hydrodynamic parameters | HK1 | 2.11 | 2 | 25 | 598.40 | 300 | ||

2 | HK2 | 2.45 | 2 | 26 | 432.78 | 500 | |||

3 | HK3 | 6.51 | 2 | 27 | The parameter corresponding to the inlet boundaries of groundwater flow | 144.60 | 200 | ||

4 | HK4 | 19.54 | 2 | 28 | 265.16 | 400 | |||

5 | Sy1 | 0.12 | 0.2 | 29 | 349.18 | 300 | |||

6 | Sy2 | 0.19 | 0.2 | 30 | 483.91 | 300 | |||

7 | Sy3 | 0.14 | 0.2 | 31 | 375.25 | 300 | |||

8 | Sy4 | 0.10 | 0.2 | 32 | 464.29 | 500 | |||

9 | Sy5 | 0.19 | 0.2 | 33 | 548.83 | 350 | |||

10 | Recharge parameter (due to precipitation and return of water from agriculture, drinking and industry) | 0.0005 | 0.0004 | 34 | . | 464.53 | 300 | ||

11 | I_{Re}(2) | 0.0008 | 0009 | 35 | 148.90 | 200 | |||

12 | 0.0007 | 0.0009 | 36 | 942.21 | 1,000 | ||||

13 | 0.0004 | 0.0001 | 37 | 366.96 | 300 | ||||

14 | 0.0007 | 0.0001 | 38 | 514.77 | 300 | ||||

15 | 0.0003 | 0.0005 | 39 | 621.65 | 700 | ||||

16 | 0.0002 | 0.0001 | 40 | 608.36 | 700 | ||||

17 | 0.0002 | 0.0009 | 41 | 156.70 | 131.50 | ||||

18 | 0.0005 | 0.0001 | 42 | 170.34 | 300 | ||||

19 | 0.0002 | 0.0001 | 43 | 501.83 | 800 | ||||

20 | 0.0001 | 0.0001 | 44 | 665.89 | 800 | ||||

21 | 0.00014 | 0.0001 | 45 | 279.72 | 300 | ||||

22 | River parameter | 414.68 | 500 | 46 | 179.35 | 109.21 | |||

23 | 1,132.70 | 300 | 47 | 121.78 | 200 | ||||

24 | 1,113.10 | 300 | 48 | 69.48 | 50 |

Number . | Parameter . | Best value (SUFI-2) . | Best value (PSO) . | Number . | Parameter . | Best value (SUFI-2) . | Best value (PSO) . | ||
---|---|---|---|---|---|---|---|---|---|

1 | Hydrodynamic parameters | HK1 | 2.11 | 2 | 25 | 598.40 | 300 | ||

2 | HK2 | 2.45 | 2 | 26 | 432.78 | 500 | |||

3 | HK3 | 6.51 | 2 | 27 | The parameter corresponding to the inlet boundaries of groundwater flow | 144.60 | 200 | ||

4 | HK4 | 19.54 | 2 | 28 | 265.16 | 400 | |||

5 | Sy1 | 0.12 | 0.2 | 29 | 349.18 | 300 | |||

6 | Sy2 | 0.19 | 0.2 | 30 | 483.91 | 300 | |||

7 | Sy3 | 0.14 | 0.2 | 31 | 375.25 | 300 | |||

8 | Sy4 | 0.10 | 0.2 | 32 | 464.29 | 500 | |||

9 | Sy5 | 0.19 | 0.2 | 33 | 548.83 | 350 | |||

10 | 0.0005 | 0.0004 | 34 | . | 464.53 | 300 | |||

11 | I_{Re}(2) | 0.0008 | 0009 | 35 | 148.90 | 200 | |||

12 | 0.0007 | 0.0009 | 36 | 942.21 | 1,000 | ||||

13 | 0.0004 | 0.0001 | 37 | 366.96 | 300 | ||||

14 | 0.0007 | 0.0001 | 38 | 514.77 | 300 | ||||

15 | 0.0003 | 0.0005 | 39 | 621.65 | 700 | ||||

16 | 0.0002 | 0.0001 | 40 | 608.36 | 700 | ||||

17 | 0.0002 | 0.0009 | 41 | 156.70 | 131.50 | ||||

18 | 0.0005 | 0.0001 | 42 | 170.34 | 300 | ||||

19 | 0.0002 | 0.0001 | 43 | 501.83 | 800 | ||||

20 | 0.0001 | 0.0001 | 44 | 665.89 | 800 | ||||

21 | 0.00014 | 0.0001 | 45 | 279.72 | 300 | ||||

22 | River parameter | 414.68 | 500 | 46 | 179.35 | 109.21 | |||

23 | 1,132.70 | 300 | 47 | 121.78 | 200 | ||||

24 | 1,113.10 | 300 | 48 | 69.48 | 50 |

*Note*. HK: hydraulic conductivity; Sy: specific yield; *I*: recharge rate; *C*: hydraulic conductivity of the riverbed material; *T*: boundary conductance of the aquifer.

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.

Criteria . | PSO . | SUFI-2 . |
---|---|---|

RMSE (m) | 0.8 | 0.86 |

NSE | 0.75 | 0.67 |

PBIAS (%) | −0.63 | −0.95 |

Criteria . | PSO . | SUFI-2 . |
---|---|---|

RMSE (m) | 0.8 | 0.86 |

NSE | 0.75 | 0.67 |

PBIAS (%) | −0.63 | −0.95 |

## CONCLUSION

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.

## DATA AVAILABILITY STATEMENT

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