Simulating Quar-Maharlu Aquifer in Fars Province, Iran, and Optimizing Artificial Recharge Using PMWIN Model

1. Introduction
Today, the improper operation, management, and restriction of water resources are serious concerns. Population growth, rising per capita water use, and the occurrence of droughts and destructive floods highlighted the need to improve the current water resources management [1, 2]. In this management, a better understanding of aquifer performance in natural conditions and anticipating the effects of water harvesting and feeding can help find a solution to save and balance the groundwater [3, 4]. The uncontrolled water exploitation may have consequences, such as permanent static water-level drop, reduction of the groundwater storage volume, and undermined groundwater quality [5]. To compensate for the aquifer depletion, one of the management procedures is artificial groundwater recharge using winter and spring floods, as well as surface water collection during the rainy seasons [2, 6].
Artificial recharge, known as aquifer reinjection, is the process by which groundwater storage in an aquifer is augmented at a rate much higher than those under the natural condition of percolation [7]. In addition to the water storage for future requirements, the artificial recharge of aquifers affects water quality by diluting and making it suitable for reuse. It can also improve land subsidence, the hydrologic regime of an aquifer, groundwater flow directions and rates, movement of saline groundwater, and leakage from other sources [5]. There are several artificial recharge methods, including surface or spreading method, flooding, basin or percolation tanks, stream augmentation, injection wells, ditch, furrow systems, and dug wells [6, 8, 9]. By simulating an aquifer in compliance with natural conditions, we can easily predict the future conditions on the substrate by considering different scenarios and examining the impact of different management methods on developing groundwater level [10].
Groundwater simulation models play a decisive role in developing and managing groundwater resources. They also have an essential role in predicting the outcomes of water management [11].
Various software applications have been used for this purpose. Emace et al. identified a hydrological system and predicted water levels and fluctuations relative to pumping and future drought potential [12]. They used the MODFLOW-96 3D numerical model to simulate the upper and middle trench aquifer in the Hill County area of southern Texas, USA.12 Ren X and Minsker B studied the effect of network size on hypothetical models of MT3D and MODFLOW. They used a genetic algorithm to model and optimized the system for the correction of pumping contamination and purification [13]. Another study was conducted to validate hydrodynamic parameters of the eastern aquifer of Bakhtegan Lake (Neyriz) in Iran using V.MODFLOW simulation [14]. Moreover, in several studies, the MODFLOW software was used [15, 16, 17, 18, 19, 20].
Paul combined the WETPASS and MODFLOW models to study the effects of land use on groundwater discharge and recharge in Eastern China in 2006. The results indicated that urban development and agricultural lands are generally the main causes of recharge decline in the studied area. However, forests can effectively maintain the aquifer and its dependent ecosystems [21]. Also, some studies have been conducted on artificial recharge and the application of PMWIN. For example, the effects of construction on artificial recharge on the groundwater level of the Marvast plain were monitored by analyzing the data from 1986 to 2009 in Yazd, Iran [22]. Choopani et al. used MODFLOW to examine artificial and natural recharge of the Sarkhahan plain underground water in Hormozgan Province, Iran [23].
In addition, the PMVIN software code has been used in several studies to model the Zarghan plain aquifer numerically in terms of hydraulic and propagation and spreading of contamination, [24]. to investigate the effects of artificial recharge on the hydraulic head in different regions of Bakan aquifer, [25] to model Ray City aquifer quantitatively, [26] to study the phreatic zone of Jiroft plain Groundwater, [27] and to evaluate the hydraulic behavior of the aquifer in Shiraz plain (all in Iran) [28]. 
Because of the recent droughts, uncontrolled abstraction of aquifers, and the need for artificial recharge in this plain, it is essential to identify best management practices and tools. The present study aims to simulate the Quar Maharlu aquifer (in Fars Province, Iran) quantitatively, optimize the location of artificial recharge, apply an artificial recharge plan at selected locations, and determine the recharge effect on the upstream and downstream heads using PMWIN software.
2. Materials and Methods
Study area
The area of Bakhtegan and Maharlou lakes is located in Fars Province. This area extends from north to Karabakh and east to Sarvestan (Figure 1).

The size of the study area is 323 km2, with a height of 175.9 ft and a plain of 147.1 km2. Also, the area of the alluvium aquifer is 131.5 km2. The mean temperature is 15.9°C for the whole area, 12.73°C in heights, and 17.1°C in plains. In the entire area, there was neither surface/underground flow nor artificial water transfer or construction from other neighboring areas. The maximum height of this area is 2202 m in the south, and the minimum is 1500 m in the plains of the center. Also, the average height of the plains and mountains are 1593.3 and 1948.91 m, respectively. The most important city in this area is Quar.
At the altitudes of Quar Maharlu, there are some formations such as Asmari-Jahrom, Sarvak, Tharbur, Darian, Saghun, Gurpi, and Kazhdumi. Also, the average annual rainfalls in the plain and altitude were 491.1 and 520.9 mm, respectively. The average annual evaporations from the plain and altitudes of this range were 2257.8 and 1587.2 mm, respectively. The wells are the primary utilization in this area. The general trend of groundwater hydrology was descending, which indicates the continued loss of groundwater resources [29].
Modeling
Conceptual Scheme: A mathematical model is a management tool for aquifer studies and predicting the effects of various hydrological and hydrogeological stresses [30]. The PMWIN software is used for the 2D modeling of a region. Therefore, it is impossible to use the digital elevation model (DEM) map for the topographic file of the region. According to the geological data of this area, drinking water wells log information, and existing reports, there is an unconfined aquifer type. Also, according to the steady state modeling, there is only one height for the floor (altitude=0) and one for the high level (height=1600 m). Therefore, the digits related to the altitude of the surface of the observed well are included from the high level. Besides Darcy’s law, the continuity equation or conservation of mass is the second most important law in ground water equations. By combining these two laws, a partial second-order equation is obtained for groundwater. If we know the dimensions of the aquifer, boundary conditions, initial conditions, and hydrological characteristics, it is possible to find the hydraulic head in different locations by solving these equations. Assuming a homogeneous aquifer, k is independent of X, Y, Z, and Equation 1 is given as follows:

Equation 1 is known as the Laplace Equation, [31], which governs groundwater flow in a homogeneous environment or steady-state condition that can be solved using analytical methods. Under no steady circumstances, it is impossible to use the Laplace method. Depending on the finite and free aquifer, the phrase would be different. In a free aquifer, two separate processes occur. First, the aquifer and water compaction will change the potential of groundwater. Also, the special storage coefficients are being applied to all the elements. In addition, the loss of free water levels leads to recharge from the aquifer. In a free aquifer, the special drainage is often equivalent to special porosity. Equation 2 illustrates the governing of the 3D flow of groundwater in an open aquifer:

Among the numerical methods for solving these equations, the finite-difference method is the most widely used one, in which each equation leads to several partial equations. For each node, a multiplicative equation is obtained in each period. For example, if a model has N nodes and M is the period, then the N*M equations are obtained. The solution matrix is used to solve the equations. To model the range, 245×245 m cells were used to construct the main grid. The number of rows and columns of the model were both 100, and the number of layers was considered to be one unit due to the existence of only one aquifer layer in the region. The total number of cells is 10000, and the number of model active cells is 5439. The area of the model is 326 km2. The base year for the scenario was 2014-2015.
Determining the initial conditions and selecting the time steps: To solve the equations governing the groundwater by a numerical method, one of the initial conditions for the computation model was being in a steady state. In this case, only the initial conjecture for the repetitive solution is introduced to the system. However, the selection of inappropriate initial conjecture when the initial level in the cell is less than the level of the floor in the same section will lead to cell drying in the calculations. In the PMWIN model, various options have been proposed to select the initial conditions, such as the provision of initial conditions (e.g., cell property inputs), using the introduced file into the model, or using a similar presented level for all cells. Inputs for the one-year time interval (due to availability of average annual discharge of wells) were entered as average, and the outputs were obtained based on inputs of the previous year. Therefore, we will have a model that gives us the status in the next year if we enter the average annual information.
Determining the values of recharge and discharge of cells: In this study, based on the reports from the upgraded Water Resources Atlas of the catchment area of Tashk-Bakhtegan and Maharlo Lakes, the water balance was as follows: 
The average annual discharge of this study area is 0.56 m3/s (0.37 and 0.19 m3/s in the highlands and plain, respectively), equivalent to 17.72 mm3/y. Based on the hydro-climatological balance (Table 1), the total volume of rainfall is reached by penetration, evapotranspiration, and surface flow (runoff).

The percentages of the three parameters are variable according to geological characteristics, slope, vegetation, and other factors. To estimate the penetration percentage, the sum of evapotranspiration and runoff percentage has been deducted from 100. Runoff and evaporation rates were determined by the equations of the Indian agriculture association (Equation 3) and the White Torrent (Equation 4), respectively:

In Equation 3, T, P, and A are the temperature, rainfall, and domain area, respectively. Also, in Equation 4, ETp is the monthly potential of the evapotranspiration (mm), and tn is the average monthly temperature (°C). Furthermore, I, in, and a are the annual and monthly heat index and dependent coefficient, respectively.
Rainfall was considered the main input factor in the aquifer at the heights, equivalent to 91.63 mm3. Also, the evaporation values of the free surface, spring, and the qanat were 0 at heights. The underground drainage flow to alluvium and the groundwater transportation by a well was 36.82 and 2.30 mm3/year, respectively, with the evapotranspiration of 91.63 mm3. The total amount of recharge and discharge of the alluvial aquifer is 91.16 and 100.04 mm3/y, respectively. The total annual discharge from the alluvial aquifer (98.89 mm3) and the hard formations (2.3 mm3) was estimated in this studied area to be 101.9 mm3. The general balance of water showed that the total amount of inlet and outlet water was estimated to be 163.87 and 172.75 mm3/y, respectively. Moreover, the variation in the mm3 reserve was 88.8% [29].
As PMWIN does not distinguish the wells between alluvial and mountainous regions, 882 wells (872 in alluvial sources and 10 sources in hard formations) were imported annually in the model (pumping rate was introduced as m3/d). The negative pumping rate indicates the discharge, and the positive pumping rate expresses the recharge to the well. Because of countless wells, the dimensions in the study area, and the size of the cells, more than one well has been located in some cells, in which the pumping rate of the target cell was cumulative. In the cells where the wells are located, the values of discharge are separately entered.
In this study, data from 5 rain gauge stations were analyzed by the Surfer program and mapping the isohyetal line. Since precipitation values are entered into the PMWIN software in the recharge pack, rainfall data should be calculated according to penetration and evapotranspiration coefficients. Afterward, the penetration rates in mountainous and alluvial regions were introduced into the software using Equations 3 and 4. A similar trend was repeated in the following years based on the rainfall of the same year.
Agricultural, drinking, and industrial water return rates were entered into the software based on the existing data in the water balance document. Accordingly, consuming values were applied to the cells as recharge. It is noteworthy that the software shows the well at a point when the cell recharge is applied. The effects of evaporation (discharge sources), runoff, transitional flow regime of water, and input flow from adjacent basins (feeding or discharge sources) were calculated and entered into the software. To apply the calculated recharge package, the area must be divided into 14 parts. Table 2 lists the applied recharging values.

The values in Table 2 are derived from the aquifer groundwater balance. So, according to the different input and output items from each area and their precipitation, the net recharge value for that area is estimated. Recharge values are negative in the mountainous and positive in the aquifer areas. Also, because only hydraulic conductivity (K value) is considered in the solution of the model in the steady-state condition, the storage coefficient does not play a role in the calculations. Several restricted areas (13 areas) were defined for the hydraulic conductivity coefficient at the basin.
Calibration and verification of the model
Calibration minimizes the difference between the calculated value and the observed value.1 To simulate the model each year, the piezoelectric data were used at the beginning. The end of the year was also introduced in the software to do the calibration process in 2012, 2013, and 2014 and verification in 2015. Although the water level at the beginning of the period in a steady state is not necessary for the simulation, the initial conditions must be entered into the model to start the calculation. In this step, we calculated Root Mean Square Deviation (RMSD) (Equation 5), a statistical quantity for the variable values:

where n is the number of observed wells and x1, x2, … xn are the differences between observed and estimated values in the first, second, …, and the nth wells. In the verification process, plain information in 2015, including average wells discharge, rainfall, water exchange at the borders, increase in the number of wells, drying wells, return water per year, drinking and industry return water, were defined in the software and calibrated by the results for each year. 
Determining the aquifer boundary conditions
Solving partially differential equations of groundwater through numerical methods requires hydraulic determination of boundary conditions and the use of information at these boundaries. The aquifer boundaries are often presented as boundary lines without flow in limited parts (in the vicinity of adjacent aquifers) as the general head boundary (GHB) boundary conditions. By evaluating the equilibrium maps of the surface water, groundwater topographic maps, basin geological maps, and general satellite images of the region, the exact location of boundaries and height of the identified basin and impermeable structures are determined as non-flow boundaries.
Investigating the artificial recharge effect on hydraulic condition in the aquifer 
In this study, the flood spreading method, the easiest and most cost-effective method for artificially recharging aquifers, was applied. Site selection of artificial recharge is the main challenge in this method. At first, 16 points at the aquifer were selected based on the problems in agricultural land ownership that could solve the ownership problem (Figure 2).

Afterward, these points were gradually eliminated concerning other environmental factors, such as distance between the selected site and city or village, proximity to the road and the industry at the plains, the existence of qualitative monitoring wells, distance to drinking water, as well as the consideration of the hydraulic conductivity coefficient (K).
Because of the proximity of urban and rural areas, three points were removed. In addition, a point was removed due to the proximity to the Pasargad Quar smelting plant. Due to inappropriate location and distance from the monitoring wells, two downstream points were also eliminated. Two points on the upstream were removed due to their lack of effect on drinking wells. Among the four remaining points in upstream, two points were in the east, between regions 9 and 10, where the division of the hydraulic conductivity existed. Finally, the point located in area 10 was selected because it had a higher K value than the others. Regarding the points located in zones 6 and 8, we selected area 6 due to its higher K value. Four locations were considered for the ditch construction (Figure 3).

A good area for flood spreading in each area was 9 cells equivalent to 0.54 km2.
Since the runoff coefficient was 12% at the heights, the total amount of runoff from mountainous areas was then calculated first to determine the total rate of recharge. Then, depending on the location of the 4 selected points and the altitude line of the heights, these areas were divided into 6 sections. The final runoff will reach the 4 selected points from one of these six areas. Total runoff was equal to 1.055 m3/y, calculated using average rainfall of 300 mm/y, an area of 175.9 km2, and dividing the mountainous areas into six parts. To optimize the recharge rate, the model was evaluated for each scenario by assuming 25%, 35%, and 50% of the total runoff volume for all selected points.
General head boundary 
GHB is determined in the model using the distance from the nearest observed head, the observed head surface at the closest hydraulic counter head, and the average hydraulic conductivity at the boundary.
K-zonal sensitivity analysis 
In this study, the effect of network-shifting and the hydraulic conductivity zoning parameter was considered in the calibration. As mentioned, this area was divided into 20 sections, including 2500 cells which contained 50 rows and 50 columns (490 m at 490 m) for each cell based on the geological maps and the hydraulic conductivity parameter. The recharge and discharge data of the area and the exchange of water with surrounding boundaries were input into the software for the year 2012. Also, groundwater level at the beginning and end of the period was used in observation wells to help calibration process and preliminary estimate hydraulic coefficients.
3. Results and Discussion 
Calibration
Because of the proper dispersion of available data in the basin and the acceptable number of observational data, it is possible to use the calibration and verification processes of the model desirably. Thus, the years 2011-2012, 2012-2013 were used for calibration and 2012-2014 for verification. The calibrated results by the reliable model have been shown in Figure 4 and Table 3 (Kx=Ky=10 Kz).

According to Table 3, the most relevant results in 2012, 2013, and 2014 were obtained in wells 13, 7, and 3, respectively.

Also, the most unsuccessful results in the same years were obtained in wells 3, 6, and 5, respectively. RMSD values were 1.55, 1.49, 1.50 m in 2011, 2012, 2013, respectively. The total RMSD was equal to 1.519 m which was acceptable according to the observed and estimated heads value (i.e., 1500 m). Finally, the average hydraulic conductivity coefficients estimated during three calibrated years were considered as the hydraulic conductivity coefficients.
The hydraulic conductivity coefficients (i.e., the coefficients obtained from different calibrated years) and the final values are presented in Table 4.

Verification
In the calibration phase, it was expected that the differences between the observed and estimated heads were low. The model’s accuracy is shown by the results of verification in the following year. Figure 5 shows the results.

The minimum and maximum differences between the observed and estimated heads in 2014-2015 were 0.847 and 2.568 m, respectively. The RMSD in this year was 10766 m, which was approximately equal to 1500 m, which seems to be acceptable according to the observed and estimated head. Table 5 lists the difference in the observed and estimated heads for piezometric wells in 2014-2015.

The results show that despite fluctuations in water levels in different months, the predicted water level in the following year gives us an appropriate approximation of the actual water level. The reason for this is the use of the average annual data, which according to the monthly fluctuations, the water level in the next year is predicted to be at an acceptable limit. The sensitivity analysis of the K zoning was defined as the effect of variation in the network and region of the hydraulic conductivity parameter in 20 parts on the results. For example, the minimum and maximum differences between the observed and estimated heads in 2011-2012 were 0.776 and 2.703 m, respectively. Also, the RMSD was 1.866 m, indicating an increase in the number of K zoning developed model’s error which is not suitable.
Figure 6 shows the hydraulic head values of different parts of the plains according to the color layout.

Light green indicates the largest head, and the dark green represents the smallest hydraulic head. Therefore, as we move from light green to dark green, the hydraulic head is reduced. As seen, the head decreases from the south and southwest towards the north and east of the plain. It can be concluded that the head is the highest around the city of Quar and its nearby heights. This finding is confirmed by the existence of most drinking wells in this range of plain. Also, the results of velocity vectors showed that the magnitude of vectors is higher in mountainous regions than in others. GBH is shown in Figure 7.

Evaluating artificial recharge effect on the hydraulic status of the Aquifer
The position of the quadruple points selected for artificial recharge using the flood spreading method was as follows:
• Area 1 is located upstream of the aquifer and close to the Thalweg of western heights of the plains. It is not in the path and has a reasonable distance from the drinking wells (K =337 m/d).
• Area 2 is located at the bottom of the Thalweg at the western heights of the plain and close to the Shiraz-Firoozabad route. There is a small distance from the rural areas and a reasonable distance from the drinking wells (K=273 m/d).
• Area 3 is upstream and located in the Thalweg of the eastern heights. It is close to the Shiraz route of Jahrom, a small distance from the rural and cultivated areas (K=413 m/d).
• Area 4 is located upstream in the Thalweg of the eastern heights. It is near the Shiraz-Jahrom route, a short distance from rural areas (K=503 m/d). Figure 8 illustrates the head contours caused by feeding in quadruple areas.

Assuming 25% runoff control as recharge rate at points 1, 2, 3, and 4, the average head rise in the observed wells was 3.5, 4.6, 1.5, and 7.3cm at the plain surface, respectively.
The results showed that if the rate of 35% is applied at points 1, 2, 3, and 4, the average hydraulic levels of the plain increase to 9, 8, 10.2, and 6 cm, respectively. Also, the average hydrological levels of the plain increase by 15.2, 11, 13.7, and 8.4cm if the 50% restraining rate is applied for the four points. These values are presented in Table 6.

Applying 25%, 35%, and 50% rates at point 1 resulted in an increase of 5.3, 9, and 15.2cm in the mean hydraulic level of the plain.
According to Table 6, the largest and smallest values result in an increase in the head at the hydraulic level of the observed wells:
• When recharge is carried out at the first point, the maximum and minimum head variations will be in wells 3 and 6, respectively.
• By reinjecting at the second point, the maximum head variation is seen in wells 4 and 7, and the minimum will be at point 2.
• In the case of recharging at the third point, the maximum head variation is in well 3, and the minimum will be at points 4 and 5.
• The recharging at the fourth point results in the maximum head variation in well 5. No change in the height head was observed in wells from 9 to 13.
Groundwater flow models have been used as assessment tools for evaluating recharge, discharge, quantifying sustainable yield, and aquifer storage processes. In addition, as a supporting tool for field data collection and practical solution designs, this model can predict future conditions or effects of human activities [32]. In this study, modeling is a screening tool for evaluating the groundwater development scenarios. To determine the areas for recharge, it is crucial to be in locations where groundwater levels are declining due to over-exploitation or where groundwater quality is poor, and there is no alternative source for water. In addition, water availability from wells and hand pumps is inadequate during the lean months in these locations. The recharge process causes a reduction in runoff and an increase in water availability, especially in the summer months. Also, ground availability would remain a limitation to the full development of this process. The quality of raw waters that are available for recharge is essential. Therefore, before using, some sources require treatments. A significant requirement is being silt free that can be used for the recharge process. As seen, with the rates of 25% and 50%, point 1 is the best option, while with the rate of 35%, point 3 is the best. These points are both upstream. The main reason for the results of these points can be due to the movement of groundwater and the proper position of these points relative to the hydraulic head gradient.
These points also have a positive effect on downstream drinking wells. As we move to the north, the artificial recharging effect decreases. Also, the closer to the recharging point, the more positive effect will be observed in the points. A study investigating the effects of Dauder dam construction on the aquifer of the Ladiz plain and aquifer management using the Visual MODFLOW premium 4.2 model found that an artificial recharging basin for water injection into the aquifer in March, April, May, and June are the best options [33]. A modeling study of the artificial recharging of the Dasht Naz in Sari City, Iran, showed that the nearby southwest areas of Neka and the south Tajan rivers were the most suitable options for recharging in the Dasht Naz [34].
Previous studies have shown that the annual drop in groundwater levels in the following years of artificial recharging has been significantly lower than in previous years.6 Performance evaluation of the artificial recharging plan of Ab-barik aquifer (Bam) using MODFLOW software revealed that the annual increase of 12.3 mm3 of water could not stop the decline in the phreatic zone in the aquifer of narrow Bam water [35]. According to the mentioned issues, the result of the present study about the ability and efficiency of the MODFLOW model to optimize the location of artificial recharge is similar to that of Mahdavi et al. [36].
4. Conclusion 
Improper exploitation and management of groundwater resources are the main factors causing the water crisis in Fars Province. The cooperation of responsible organizations can prevent this crisis. Artificial recharge as a solution to the issue is a design developed by humans in which the water is injected into the ground in a controlled way. In this regard, the models can determine the optimum recharge locations. The results of this study indicated the necessity of an artificial recharge process and simulation using the PMVIN software that can be well suited for determining the best alternatives among the defined scenarios.

Ethical Considerations
Compliance with ethical guidelines
There were no ethical considerations for this research. 

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

Authors' contributions
All authors equally contributed to preparing this article. 

Conflict of interest
The author declared no conflict of interest.

Acknowledgments
The authors would like to thank the Department of Civil Engineering, Shiraz University, for their assistance.

References

1. Ehteshami M, Sharifi A. [Finite difference modeling of rey province ground water (Persian)]. J Water Soil Sci. 2008; 22(2):257-70. [Link]
2. Mahdavi A, Nouri Emamzadei M, Mahdavi Najafabadi R, Tabatabaei H. Identification of artificial recharge sites using fuzzy logic in Shahrekord basin. Water Soil Sci. 2011; 15(56):63-78. [Link]
3. Ghobadiani R, Fatahi A, Majidi S, Mohamadzare M. Groundwater reservoir simulation of Kermanshah Miandarband Plain by GMS Model. The First National Conference on Challenges on Water Resources & Agriculture. 13 February 2014; Esfahan, Iran. [Link]
4. Yousefi Sahzabi H, Zahedi S, Niksokhan MH. Ten year prediction of groundwater level for the purpose of of determining reasonable policies for exploitation from aquifer. Ecohydrology. 2017; 3(3):405-14. [DOI:22059/IJE.2016.60028]
5. Ministry of Energy. Drinking water quality criteria. In: Power Mo, editor. Technical office; 1988.
6. Ghordemi Milan M, Karami GH. Evaluation of the effects of Khoy artificial recharge In: University MF, editor. 8th Conference of the Iranian Association of Engineering Geology and the Environment. 14 February 2013; Mashhad, Iran. [Link]
7. Asano T. Artificial recharge of groundwater. Amsterdam: Elsevier; 1985. [DOI:10.1016/B978-0-250-40549-7.50006-0]
8. Prickett TA. Modeling techniques for groundwater evaluation. Adv Hydrosci. 1975; 10:1-143. [DOI:10.1016/B978-0-12-021810-3.50006-0]
9. Ayers RS, Westcot Water quality for agriculture: Food and Agriculture Organization of the United Nations Rome; 1985. [Link]
10. Javadi S. Predict and determine the risk of contamination of groundwater aquifers by rating analytical and Monte Carlo technique. Tehran: Tarbiat Modarres University; 2010.
11. McDonald MG, Harbaugh AW. A modular three-dimensional finite difference ground-water flow model. Reston: U.S. Geological Survey; 1988. [Link]
12. Emace R, Chodhury A, Anaya R, Way SC. A numerical groundwater flow model of the upper and middle Trinity Aquifer, Hill Country area. Texas: Texas Water Development Board; 2000. [Link]
13. Ren X, Minsker B. Which groundwater remediation objective is better: A realistic one or a simple one? J Water Resour Plann Manag. 2005; 131(5):351-61. [DOI:10.1061/(ASCE)0733-9496(2005)131:5(351)]
14. Majedi H, Sorori M. Application of simulation process to certainty the hydrodynamic parameters of aquifer using a V.MODFLOW mathematical model in the aquifer of Bakhtegan Lake East (Neyriz). 05th Iranian Hydraulic Conference; Kerman2005. [Link]
15. Ayenew T, Demlie M, Wohnlich S. Application of numerical modeling for groundwater flow system analysis in the Akaki catchment, central Ethiopia. Math Geosci. 2008; 40(8):887-906. [DOI:10.1007/s11004-008-9144-x]
16. Joodi A. Investigation the possibility to balance groundwater Faizabad aquifer by the mathematical model. Shiraz: Shiraz university; 2009.
17. Pourjenabi A. Determining of hydrogeological regions and groundwater balance calculating at stable condition by MODFLOW model (Case Study: Sarzahrzvan Plain, Hormozgan). 5th National Conference on Watershed Management and Soil and Water Resources Management. 29 February 2012; Kerman, Iran. [Link]
18. Asghari k, Sourinejhad J, Zolanvar AK. [Optimal Operation of Borkhar Aquifer Using Simulation - Optimization Methods (Persian)]. J Water Soil S 2005; 9(3):13-25. [Link]
19. El Yaouti F, El Mandour A, Khattach D, Kaufmann O. Modelling groundwater flow and advective contaminant transport in the Bou-Areg unconfined aquifer (NE Morocco). J Hydro-Environ Res. 2008; 2(3):192-209. [DOI:10.1016/j.jher.2008.08.003]
20. Thorley M, Call P, Williams H, Lough H, Kininmonth M, Henderson B. Christchurch city groundwater model. 2005.
21. Paul MJ. Impact of land-use patterns on distributed groundwater recharge and discharge. Chin Geogr Sci. 2006; 16(3):229-35. [DOI:10.1007/s11769-006-0229-5]
22. Raesai R, Shabani M, Abaszadeh M. The effect of dams artificial recharge of ground water surface Case Study: Marvast Plain, Yazd. 1st Regional Water Resources Development Conference. 19 May 2011; Abarkuh, Iran. [Link]
23. Choopani S, Rastegar H, Hoseinipoor H. Evaluation of quantitative changes in ground water table affected by Sarchahan’s Flood spreading plan by mathematical model. Water Manag Arid Lands. 2014; 1(1):63-73.
24. Mojarrad M, Rakhshandehrou GR, Monadi R, Ghorbani Musellu M. Investigating wellhead protection area (whpa) assessment methods, case study: Beyza-Zarghan drinking wells. 2nd International Conference - Water resources and wetlands; Tulcea (Romania). 2014; Shiraz, Iran.
25. Peirovi F. The effects of artificial recharge on groundwater quality and quantity of Bakan aquifer using the PMWIN software Shiraz. Shiraz: Shiraz University; 2016.
26. Ehteshami M, Sharifi A. Evaluation of Ray city aquifer quality model. Journal of Environmental Science and Technology. 2007; 8(4):1-10.
27. Pourseidi A, Kashkouli HA. Studying of Groundwater Conditions in Jiroft Basin with Modflow. Irrig Sci Eng. 2012; 35(2):510-63. [DOR:20.1001.1.25885952.1391.35.2.6.2]
28. Karimipour AR, Rakhshandehroo GR. Sensitivity analysis for hydraulic behavior of shiraz plain aquifer using pmwin. J Water Wastewater. 2011; 22(2):102-11. [Link]
29. Ministry of Energy. Report of balance the study area Kovar Maharloo Ministry of Power, 1389.
30. Vafaiyan M. Embankment dams: Jihad Daneshgahi Esfahan University; 2003.
31. Nikbakht J, Najib Z. The effect of increasing the efficiency of irrigation on groundwater level fluctuations (Case study:Ajabshir plain-East Azerbaijan). J Water Irrig Manag. 2015; 5(1):115-27. [DOI:10.22059/JWIM.2015.55221]
32. Zhou Y, Li W. A review of regional groundwater flow modeling. Geosci Front. 2011; 2(2):205-14. [DOI:10.1016/j.gsf.2011.03.003]
33. Salari M. Evaluation of the effects of Duder dam construction on Ladiz aquifer and management of aquifer by mathematical models. Sistan and Baluchestan: University of Sistan and Baluchestan; 2008.
34. Ashtianimoghadam G, Ahmadi MZ, Kholghi M, Safaiian N. Artificial recharge model of Dasht Naz in Sari in order to prevent advancing of salt water. Journal of Agricultural Sciences and Natural Resources. 2003; 1(4).
35. Katibeh H, Hafezi S. Application of Modflow in Groundwater Management and Evaluation of Artificial Recharge Project of Ab-barik Aquifer (Bam). J Water Wastewater. 2005; 15(2):45-58. [Link]
36. Mahdavi M, Farokhzadeh B, Salajeghe A, Malakian A, Souri M. [Simulation of Hamedan-Bahar aquifer and investigation of management scenarios by using PMWIN (Persian)]. Watershed Manag Res (Pajouhesh & Sazandegi). 2012; 1392(98):108-16. [Link]