Coupling of WetSpass-M and MODFLOW Models for Groundwater Flow Assessment

Recharge is considered a key parameter in groundwater models for sustainable management of aquifers, which is in�uenced by factors such as land use, soil, weather, etc. The present study was conducted to couple WetSpass-M and MODFLOW models for evaluating Neyshabour aquifer condition in steady and transient states. To this aim, the simulated recharge by the WetSpass-M model was applied as an input of MODFLOW to assess the groundwater balance. The hydrodynamic coe�cients were determined by calibrating the model, evaluating and the model sensitivity to the hydraulic conductivity coe�cient, speci�c yield (S y ), and recharge. The results indicated that the annual average of surface runoff, actual evapotranspiration, interception, and recharge during 1991–2017 equaled 18, 36, 7.6, and 42.6% of the average annual precipitation in the basin, respectively, with the simulated water balance error 4.2%. The average annual recharge of the basin varies between 0-257.41 mm with an average of 105.25 mm/y. Accordingly, the maximum and minimum average monthly recharge occurs during March and July, respectively. The appropriate matching of the simulated and observed water levels and obtaining the suitable values of RMSE, R 2 , ME, and MAE evaluation criteria in steady and transient states indicate the adequate accuracy of the WetSpass-M model in estimating the recharge and success of the couple two models. Based on the simulated groundwater balance, the aquifer faces a de�cit of 421.3 MCM per year and 97.41 cm in the annual groundwater level. The model displayed more sensitivity to the hydraulic conductivity coe�cient compared to other parameters.


Introduction
Groundwater has emerged as a strategic resource not only in most arid and semi-arid countries but also in humid climates during the last 50 years (Shah, 2014).The high quality and relatively more stable state of such waters compared to surface ones are among the reasons for increasing the use of groundwaters (Siebert et al., 2010).Approximately 750-800 billion cubic meters of groundwater are consumed in agriculture on a global scale (Shah, 2014).Groundwater is regarded as a basic source of fresh water for humans, which constitutes more than 70% of water supply sources in arid and semi-arid regions.These resources have been reduced due to over-harvesting for irrigated agriculture at local and regional scales during the recent decades.In addition, the signi cance of groundwater resources is predicted to increase in dry areas including Iran due to future climate changes, periods of drought, and reduced access to surface resources (Hashemi et al., 2015).Groundwater is considered a strategic and necessary source over the world due to all of the aforementioned reasons.Thus, the effective management of such resources is regarded as a critical and urgent challenge.In addition, sustainable management of groundwater requires acquiring su cient knowledge regarding recharge (Tilahun and Merkel, 2009).Therefore, determining the amount of groundwater recharge is considered as a prerequisite for e cient and sustainable management of such resources since the amount of withdrawal from the groundwater in a region in the long term should not be higher than the long-term average of its recharge.Thus, estimating the recharge determines the appropriate level of groundwater extraction (de Vries and Simmers, 2002).
Groundwater models are often applied to simulate water ow in aquifers.Such models can be used to simulate the long-term behavior of an aquifer under different management plans when they are calibrated.
Recharge is regarded as a critical parameter in groundwater models.Lack of properly estimate recharge and its spatio-temporal distribution make these models unreliable (Sphocleous, 2005).
The direct recharge of groundwater due to precipitation is often offered to the model as a percentage of the aquifer rainfall in arid and semi-arid regions (Hashemi et al., 2015), while groundwater recharge exhibits a random behavior (Sen, 2008), which is in uenced by factors such as weather conditions, variety of land use, soil, and vegetated cover temporally and spatially (Kim et al., 2008).In this study, the monthly water balance of the WetSpass-M (Water and Energy Transfer between Soil, Plants and Atmosphere under quasi-Steady State-Monthly) distributed-hydrological model was utilized to estimate the groundwater recharge since the mentioned factors cannot be considered in groundwater models, and the WetSpass-M water balance distribution model can regard the above-mentioned factors and has shown acceptable performance as a successful instrument in estimating the spatio-temporal distribution of recharge in previous studies (Shrestha et   .Then, the groundwater ow was simulated in steady and transient states with the entry of WetSpass-M recharge to the MODFLOW.Such a method is applied in the Neyshabur watershed as one of the potential agricultural basins in Razavi Khorasan province, Iran.The Ministry of Energy declared this plain as a critical forbidden one since 1987 following the excessive exploitation and sharp decline of the groundwater resources, as well as the occurrence of a negative balance.
Few studies have been conducted on spatio-temporal estimation of groundwater recharge using the WetSpass-M model and its evaluation in the groundwater balance.Determining the water balance and knowing the characteristics of the Neyshabur watershed and aquifer can be a fundamental step to providing appropriate solutions in order to eliminate the obstacles of such critical agricultural basin in Razavi Khorasan province and compensate for its aquifer de cit.

Study area
The Neyshabur watershed is regarded as a part of Kavir-e Markazi basin in Iran and a part of the Kal Shur watershed located in the range of 58°13 -59°30 east longitude and 35°40 -36°39 north latitude and the center of Razavi Khorasan province.Totally, 2883.4 km 2 is occupied by the aquifer of the basin out of the 9449 km 2 of the area.The average slope of the watershed equals 6.86%, the altitude of the study area is between 1056-3309 m, and its lowest point is located at the outlet of the aquifer, namely the Hosseinabad Jangal region.Since the average annual precipitation and evaporation of this basin are 246.86 and 2357.23 mm, respectively, and its annual average temperature of 13.3°C, in terms of climatic conditions classi ed as semi-arid to the arid region.The network of surface streams ows towards the center of the plain from the northern, northwestern, eastern, and southern elevations in the form of permanent rivers and various intermittent streams, after joining the Kal Shur as the main drain rerouted to the west and exiting from the outlet of the study area.The study area is considered as a part of the southern edge of the geological structure of Binalood and the northeastern limit of the triangular zone of central Iran in terms of geological divisions (Izady, 2014).The Neyshabur watershed is among the agriculturally prone ones in Razavi Khorasan province, which faced a severe drop in the groundwater level due to excessive exploitation of groundwater resources in the agricultural sector, so the Ministry of Energy designated such plain as a critically forbidden one since 1987.Figure 1 shows the con ne of Razavi Khorasan province, basin, and aquifer, as well as the rain gauge, evapotranspiration, and hydrometric stations used in the study.
Figure 1 Location of Razavi Khorasan province in Iran, basin in the province, and aquifer in the basin, along with rain gauge, evapotranspiration, and hydrometric stations As indicated in Table 1, the number of wells exhibits a completely upward trend during 1968-2008.So that the number of wells increased almost 13 times and the total discharge from the water resources of the basin rose 2.5 times during 40 years.Implementing the groundwater restoration and balancing plan in the basin stopped the rising trend in the number of wells and total discharge, as well as control the exploitation of the aquifer to some extent.In order to use the nite difference method in groundwater problems, several practical models have been presented by the United States Geological Survey (USGS) and made available to everyone in recent decades.
The MODFLOW model has wide applications and is highly accepted among hydrogeologists.The MODFLOW is regarded as a three-dimensional (3D) nite difference model of groundwater ow developed by the USGS (1) where K xx , K yy and K zz indicate hydraulic conductivity along the axes X, Y, and Z parallel to the main axes of hydraulic conductivity (LT − 1 ).In addition, h, w, S s , and t represent hydraulic head (L), ow ux per unit volume which indicates recharge or discharge (T − 1 ), special storage of porous medium (L − 1 ), and time (T), respectively.Eq. 1 describes the ow of groundwater under the transient state in the non-homogeneous environment provided that the main axes of hydraulic conductivity are aligned with the coordinate axes.The aforementioned equation takes different forms in steady, transient states, uncon ned, or con ned aquifers.
Eq. 1 can rarely be solved analytically, except in simple cases.Thus, different numerical methods are applied to obtain approximate solutions in solving real problems.

Groundwater recharge
Some researchers have proposed recharge as the main parameter of groundwater (Aslam et al., 2022).Bene tting from accurate information regarding the values of recharge is among the most necessary steps to preparing the appropriate groundwater model (Dowlatabadi and Zomorodian, 2015).Groundwater recharge exhibits a stochastic behavior (Sen, 2008), which is in uenced by a large number of factors such as weather conditions, land use, soil, and vegetated cover temporally and spatially (Kim et al., 2008).Based on the literature review, modelers consider a uniform amount of precipitation for recharge because the above-mentioned factors cannot be applied in groundwater models.
The WetSpass model was developed to simulate the long-term average spatial patterns of groundwater recharge (Batelaan and De Smedt, 2007).The WetSpass-M model can simulate interception of vegetated cover, surface runoff, evapotranspiration, soil water balance, and recharge on a monthly scale.The monthly recharge is calculated as the last component of the water balance through the following equation (Abdollahi et al., 2012): 2 where P, I, S V , ET V , and R V represent amount of precipitation, interception, surface runoff, evapotranspiration, and groundwater recharge (mm/mounth), respectively.Therefore, spatial-distributed recharge is estimated based on the type of vegetated cover, soil type, slope, depth of groundwater, and climate variables such as precipitation, evapotranspiration, temperature, and wind speed (Abdollahi et al., 2012).Thus, the spatiotemporal distribution of recharge is estimated by the WetSpass-M and used in the MODFLOW model considering the presentation of the successful for estimating recharge, conducting few studies in the spatio-temporal estimation of groundwater recharge utilizing a distributed-hydrological model and its evaluation for estimating the groundwater balance in dry regions with the area of more than 100 km 2 (Hashemi et al., 2015).The cell size of both models is considered 300×300 m to increase the accuracy of the input of distributed recharge values extracted from the WetSpass-M model to MODFLOW.

WetSpass-M
The WetSpass-M model needs the Digital Elevation Method (DEM) maps, slope, land use, soil texture, and the monthly distribution maps of groundwater depth, Leaf Area Index (LAI), as well as climatic data including monthly precipitation, pan evaporation, wind speed, and mean temperature during the studied period in order to simulate the hydrological components of the basin.The DEM map of the SRTM satellite was prepared by the USGS with a resolution of 30 m, and the slope map of the basin was produced using the DEM map in the ArcGIS software.Based on the land use and soil texture maps taken from the Regional Water Company of Khorasan Razavi province, the studied area has 14 land use and 6 soil texture classes.Poor rangeland and loam texture are regarded as the dominant land use and soil texture in the area, respectively.The information in the lookup tables of the model was revised based on land use and soil in Iran, and a code was assigned to each class of land use and soil in the region.Daily climate data including precipitation, mean temperature, pan evaporation, wind speed of synoptic stations, climatology, rain gauge, and evapotranspiration during 1991-2017 were received from the Iran Meteorological Organization (IRIMO) and the Ministry of Energy.Finally, 48, 19, 17, and 1 stations were selected to produce monthly distribution maps of precipitation, mean temperature, evaporation, and wind speed in the study period after reviewing the received data.The monthly information of groundwater level measurement in all of the observation wells of the basin during the studied period was achieved from the Regional Water Company of Khorasan Razavi province, resulting in considering 29 ones for producing monthly distribution maps of the groundwater level.The monthly LAI map for the studied period was received from the Moderate Resolution Imaging Spectroradiometer (MODIS).All of the model inputs, except the number of rainy days were prepared monthly with a cell size of 300×300 m and 420 columns and 440 rows in raster map with ASCII format and introduced to the model.Finally, the model was prepared to start the hydrological simulation of the basin on a monthly scale by producing and introducing 1948 maps during 1991-2017.

MODFLOW
In GMS, there are two methods grid and conceptual model to simulate groundwater.The conceptual model method was used due to its acceptable accuracy in solving real and complex groundwater problems, despite more di culty in working with such a method.A conceptual model is considered a simpli ed image of the real world, the preparation of which is regarded as a fundamental step to developing a numerical model in modeling studies.This model is considered as the basis of a mathematical one, which is prepared based on the primary information of eld data and hydrogeological interpretations.Acquiring correct knowledge regarding the hydrogeological situation of the region, explaining the problem of groundwater to prepare a numerical model, helping in selecting an appropriate numerical model, and simplifying the problem rationally by utilizing proper hypotheses are among the most critical objectives of preparing a conceptual model.A suitable conceptual model is not necessarily prepared during the rst stage of studies.
Modeling is regarded as a continuous and dynamic process in which the conceptual model is evolved and its constituent hypotheses are re-examined, added, reduced, or modi ed in accordance with the continuation of studies, an increase of knowledge about the system, and the rise in the amount of available data (Iran Ministry of Energy, 2005).The aquifer boundary in the studied area was determined by reviewing the geological maps and its formations.The network of observation wells in the Neyshabur aquifer was established and leveled in 1977 in order to measure changes in the groundwater level (Iran Ministry of Energy, 2016).Currently, there are 59 observation wells in the entire Neyshabur aquifer, whose groundwater level is measured monthly.Here, 28 observation wells that were considered active and bene tted from information during 1991-2017 were selected and their water level data was prepared during the study period.In addition, October 1999 was regarded as the steady time after analyzing the uctuations of the groundwater level in the studied time period.The groundwater of the aquifer ows from its north, northeast, east, and southeast to the west and southwest, and follows the topography and network of streams.The boundaries related to the model in all of the parts, except the small ones in the south of the aquifer, were de ned as General Head Boundary by preparing and checking a map for the ow direction of the aquifer in the selected steady time (Fig. 2).The DEM of SRTM with a resolution of 30 m received from the USGS was applied for the topography of the aquifer.It is worth noting that the height of the selected observation wells with the DEM map was controlled and corrected when necessary.Figure 3 illustrates the three-dimensional map related to the topography of the aquifer.The bedrock map was prepared using the geophysical operation soundings, drilling logs of exploration wells, and its integration with geological and hydrogeological information.
In the steady state, there were generally 1270 water sources including wells, springs, and qanats in the investigated aquifer.The 25 and 75% of discharge values of these resources were introduced to the model as return water for agricultural and non-agricultural uses, respectively (Iran Ministry of Energy, 2010).Due to the high depth of the water table, evaporation from the groundwater was not performed and considered in the model during the studied period.The river package was utilized to introduce rivers and intermittent streams of the aquifer to the model.The initial values of hydrodynamic coe cients such as hydraulic conductivity coe cient and speci c yield were estimated by the pumping tests performed in the region and then determine during the modeling process and calibration in steady and transient states.Information indicated in the GMS 10.5.8 software was recalled, and creating the conceptual model.Neyshabour aquifer as a one-layer, the uncon ned aquifer was divided into a grid containing 247 rows and 336 columns (32230 active and 50762 inactive cells) with the size of 300×300 m and introduced into the model.Finally, the created nite difference grid and conceptual model are converted into the MODFLOW-2005 model to start simulating (Fig. 4).

Calibrating the WetSpass-M model
The inputs of the WetSpass-M model were prepared during 1991-2017, and the model ran on a monthly scale.The model outputs included values of direct runoff, groundwater recharge, evapotranspiration, and interception with their distribution maps in ASCII format were presented in each month of the simulation period.The model is calibrated by comparing the observed and simulated groundwater recharge (base ow) and direct runoff.Therefore, the hydrometric stations of the basin were assessed in terms of the length of the available statistical period, missing data, and outlier data test, and ve ones were selected to calibrate the model.In order to separate the base ow from the observed discharge, the different methods of hydrograph analysis including the local minimum, one parameter digital lter, recursive digital lter, Base Flow Index (BFI), and PART model were evaluated and the most appropriate method was selected in each station.The model was calibrated manually using the trial and error method by changing the parameters such as interception parameter (a), evapotranspiration coe cient (α), surface runoff (LP), runoff delay factor (x), base ow (β), and by considering coe cients of Nash-Sutcliffe (NS), coe cient of determination (R 2 ), the model bias for water balance (MB), and root mean square error (RMSE) as objective functions and model performance evaluation criteria in selected hydrometric stations.Table 2 shows the optimal values of calibration parameters with their allowed ranges.3 represents the values of the evaluation criteria of the model accuracy in simulating the monthly total runoff of the selected hydrometric stations.Eriyeh Chaharbagh hydrometric station was selected as an appropriate representative of the basin to evaluate the model validation due to its ability to cover the largest sub-basin of the studied area.Values of evaluation criteria such as NS, R 2 , MB, and RMSE during the validation period (2014-2017) obtained 0.64, 0.5, 0.17, and 0.35 m 3 /s, respectively.Moriasi et al. (2007) indicated that the values related to Nash-Sutcliffe should be greater than 0.5 in hydrological studies, as well as pollutant transfer processes on a monthly scale in order to consider the results of the model as acceptable.The same criterion is usually utilized for the coe cient of determination, as well.According to Table 3, the values of NS and R 2 in all of the stations are higher than 0.   Figure 6 demonstrates the portion of the monthly distribution of different water balance components in the basin from the total monthly precipitation of the area during the studied period.As perceived, the lowest amount of monthly recharge occurs during the summer season (July-September) probably due to the reduction of vegetated cover and precipitation in the area during this period.The maximum average monthly recharge is observed during March and April with 19.2 and 18.61 mm/month, while the minimum recharge occurred during August and July with 0.51 and 1.03 mm/month, respectively.
As displayed in Fig. 7a, the average annual recharge in the basin due to precipitation varies between 0-257.41mm with an average of 105.25 mm/year.Then, the spatial distribution map of groundwater recharge extracted from the hydrological model results on October 1999 (steady state) was prepared as input for the MODFLOW model (Fig. 7b).As observed, the northwest areas of the aquifer with irrigated farming land use have the maximum recharge in steady state.

Calibrating the MODFLOW model
The  shows sensitivity to a parameter in which small changes in its values affect the model signi cantly (Ghodrati and Barzegari, 2016).In the sensitivity process, the sensitivity of the model is evaluated by decreasing and increasing the hydraulic conductivity, speci c yield, and recharge at the levels of 10, 20, and 30% and transient state (Fig. 10).Based on the results, the model displays the highest and lowest sensitivity to the hydraulic conductivity coe cient and speci c yield, respectively.It is worth noting that the model did not converge while applying a 10% reduction in hydraulic conductivity and a 30% increase in recharge.The 3D topographic map of the Neyshabur aquifer The conceptual model of the study and nite difference grid introduced to the MODFLOW model Hydrographs of observed and simulated heads in piezometers 21 and 27 in the transient state Model sensitivity analysis to changes of hydraulic conductivity, speci c yield, and recharge in aquifer Aslam et al. (2022) examined WetSpass-M and MODFLOW models to investigate the effects of recharge and discharge on the Chaj Doab aquifer in Punjab region of Pakistan.Values of evaluation criteria (RMSE = 0.0022, MAE = 0.5158, and ME = 0.01836) during the calibration period of the MODFLOW model in the transient state indicate the success of utilizing the estimated recharge by WetSpass-M in the groundwater model.In the next step, three management scenarios were considered to predict the future state of the aquifer.Continuation of pumping during 2003-2019 until 2035 leads to 18.1 m decrease in the groundwater level.The above-mentioned scenario predicts the most critical situation in the future of the aquifer among the three management ones.Ghouili et al. (2017) employed the WetSpass model to determine the groundwater recharge values of Takelsa multilayered aquifer (Tunisia) on seasonal and annual time scales, and then the simulating the groundwater ow in steady state using the MODFLOW model.Based on the results, surface runoff, recharge, and evapotranspiration accounted for 17.8, 4.3, and 89% of the precipitation in the region, respectively.Considering the signi cance of recharge and its fundamental role in groundwater aquifers, the WetSpass-M and MODFLOW models synchronized to achieve the following general scienti c objectives: (1) estimating the spatial distribution of water balance components in the Neyshabur watershed in monthly time steps during 1991-2017, (2) evaluating the performance of the estimated recharge using WetSpass-M model to simulate the aquifer condition in steady and transient states by the groundwater ow model, (3) determining the hydrodynamic coe cients of the aquifer.
WetSpass-M is regarded as a spatially quasi-steady distributed water balance model, which separates precipitation into interception, surface runoff, evapotranspiration, and recharge for each cell.To consider the heterogeneity of land use in each cell, four subcells including vegetated cover, bare soil, open water, and impervious surface are de ned for each land use class(Abdollahi et al., 2017).The water balance of a unique raster is achieved by summing the independent water balances of the vegetated cover, bare soil, open water, and impervious surface of a raster cell, while the water balance of the studied area is acquired by summing the water balance calculated for each raster cell(Abdollahi et al., 2012).Precipitation is considered the starting point for calculating the water balance in a raster cell and the rest of the processes including interception, runoff, evapotranspiration, and recharge are followed regularly.Batelaan  and De-Smedt (2001) presented the WetSpass model for the rst time in the Department of Hydrology and Hydraulic Engineering of the University of Brussels in order to predict hydrological processes in seasonal and annual time steps.Then, the Hydrology and Hydraulic Programming Library (H2PL) provided the possibility of simulating water balance components including interception, surface runoff, evapotranspiration, and recharge on a monthly scale by presenting a new version of the model namely WetSpass-M under the Python program independent of the Arc view software (Abdollahi et al., 2012).

(
McDonald and Harbaugh, 1988) in 1984 utilizing FORTRAN66(Todd and Mays, 2005).Various graphical interfaces were prepared for MODFLOW code.Groundwater Modeling System (GMS) is considered a comprehensive and user-friendly graphical interface for conducting groundwater simulations.GMS includes several models such as MODFLOW, MODPATH, FEMWATER, MT3DMS, and the like.The 3D ow of groundwater in a porous medium with constant density is expressed by the following differential equation(Todd and Mays, 2005):

5 ,
indicating the satisfactory performance of the model in the simulation.Appropriate values of RMSE and MB as the criteria of model error in simulation in all the hydrometric stations indicate the high e ciency, low error, and acceptable accuracy of the model during the simulation.The appropriate values of the statistical indices during the calibration of the selected hydrometric stations, as well as their values during the validation of Eriyeh Chaharbagh station, indicate the high accuracy of the optimized model during calibration and validation.

Figure 5
Figure 5 compares the hydrograph of simulated runoff during the calibration period (1991-2010) and validation (2014-2017) with the observed runoff at Eriyeh Chaharbagh station, along with the average precipitation of the entire basin.As observed, the model operates successfully in simulating the process of changes in runoff values, time of occurrence, and peak and minimum discharge values, which play a key role in basin management to confront future oods and droughts in the basin, during the calibration and validation period, so that these changes follow the precipitation pattern of the basin.Based on the results, the average annual surface runoff, actual evapotranspiration, interception, and groundwater recharge during 1991-2017 are 18, 36, 7.6, and 42.6% of the average annual precipitation in the basin and the simulated water balance error is 4.2%.The results of this research are completely consistent with the results reported by Soleimani Motlagh (2017) and Bayati (2017).While, Zarei (2016), Ashaolu et al. (2020) and Yenehun et al. (2022) introduced actual evapotranspiration, and Zhang (2018) presented surface runoff as the component that the largest contribution of annual precipitation in their studied basin.
the outlet and center aquifer, respectively, are illustrated in Fig. 9.The results of the groundwater balance simulated by MODFLOW indicate the aquifer confronts 421.3 million cubic meters (MCM) per year of the reservoir de cit, which equals 97.41 cm of annual level drop in the entire aquifer.The total results obtained from the simulation of groundwater ow show a high and suitable ability of the prepared model in simulating of Neyshabur aquifer.Further, achieving the appropriate values of the evaluation criteria for the groundwater model in steady and transient states indicates the successful coupling of the WetSpass-M and MODFLOW models through the recharge parameter.Aslam et al. (2022) also reported the desirable performance for the couple of WetSpass-M and MODFLOW models in the Chaj Doab area of Pakistan.Most of the parameters used in groundwater modeling such as hydraulic conductivity, speci c yield coe cients, and recharge have uncertainty.Recognizing such uncertainties and quantifying their change intervals affect modeling and calibrating the mentioned parameters signi cantly.The calibrated model

Table 2
Allowed range and optimal values related to calibration parameters of WetSpass-M model

Table 3
The values of statistical indicators for evaluating the performance of the model in simulation the monthly total runoff Nazarieh et al. (2018), 1998) converted into a numerical model and entered the recharge values extracted from the WetSpass-M model in order to calibrate the hydraulic conductivity coe cient of the aquifer, and nally, the MODFLOW model was run in a steady state.The calibration continued manually (trial and error method) by a logical change and consistent with the characteristics of the region in the input parameters until getting an acceptable match between the observed and calculated water levels in 28 observation wells in the aquifer.It is noteworthy that the recharge values do not change during the calibration of the model.The values of evaluation criteria RMSE = 0.55, MAE = 0.48 (Mean Absolute Error), ME = 0.0018 (Mean Error) in meters, error variance = 0.3, and R 2 = 0.99 was obtained after the model calibration in steady state.Criteria such as error variance, RMSE, ME, and MAE indicate the amount of model error in the simulation, and the closeness of their values to zero means the desirable accuracy of modeling.Figure8shows the simulated water table map, with a chart of comparison between the simulated and observed water tables in a steady state after calibration.The values of the evaluation criteria, the greening of the levels in all of the observation wells (Fig.8a), and high agreement of observed and simulated groundwater tables (Fig.8b) indicate the appropriate accuracy of the prepared conceptual model, acceptable performance of the model in simulating steady state, and accurate estimation of the hydraulic conductivity coe cient of the aquifer.For calibration of speci c yield coe cient, the model was prepared from November 1999-September 2000 at monthly time steps as the transient state.Since the MOFLOW model needs the initial head values for beginning ow simulation(Chiang and Kinzelbach, 1998).Therefore, the simulated water table in a steady state is considered the initial head in the rst time step of the transient state, and then the calculated water table of the previous month is used as the initial head of the future month.The optimum values of the hydraulic conductivity coe cient and speci c yield in different zones, with the observation wells placed in each zone, are indicated in Table4.The values of the evaluation criteria in the transient state (Table5) illustrate the appropriate performance of the model prepared for the Neyshabur aquifer in simulating the real conditions of the aquifer.These ndings are regarded as more appropriate compared to those reported byIzady et al. (2015)in the Neyshabur watershed in order to apply the simulated recharge by SWAT in the MODFLOW model and alsoNazarieh et al. (2018), indicating the higher accuracy of the WetSpass-M model in estimating recharge.The results demonstrate that the model has appropriately simulated the variation process of observation wells water table.For instance, the observed and simulated hydrograph in piezometers 21 and 27 located in