Natural vs. artificial groundwater recharge, quantification through inverse modeling

. Estimating the change in groundwater recharge from an introduced artiﬁcial recharge system is important in order to evaluate future water availability. This paper presents an inverse modeling approach to quantify the recharge contribution from both an ephemeral river channel and an introduced artiﬁcial recharge system based on ﬂood-water spreading in arid Iran. The study used the MODFLOW-2000 to estimate recharge for both steady-and unsteady-state conditions. The model was calibrated and veriﬁed based on the observed hydraulic head in observation wells and model precision, uncertainty, and model sensitivity were analyzed in all modeling steps. The results showed that in a normal year without extreme events, the ﬂoodwater spreading system is the main contributor to recharge with 80 % and the ephemeral river channel with 20 % of total recharge in the studied area. Uncertainty analysis revealed that the river channel recharge estimation represents relatively more uncertainty in comparison to the artiﬁcial recharge zones. The model is also less sensitive to the river channel. The re-sults show that by expanding the artiﬁcial recharge system, the recharge volume can be increased even for small ﬂood events, while the recharge through the river channel increases only for major ﬂood events.


Introduction
Per capita, water resource availability has dwindled rapidly during the last four decades in the Middle East.Especially groundwater has undergone dramatic changes in arid areas because of higher demand.Here, groundwater is often the main source of both drinking and irrigation water, and thus rapidly decreasing groundwater levels calls for new methods to restore water availability.
In arid regions, recharge often occurs by intermittent flow through the ephemeral river course.In Iran, for example, the magnitude of flood volume resulting from ephemeral streams is in the order of 65 billion cubic meters out of 127 billion cubic meters of the total surface water flow, most of which ends up in swamps, deserts, and the sea (Ghayoumian et al., 2007).In general, floodwater that is generated in the upper catchment flows into the ephemeral river channel, crosses the downstream plain, and flows out to the down-adjacent catchment.Lower magnitude events generally do not result in complete channel submersion on reaches where the channel is very broad (e.g.Bull and Kirkby, 2002).Hence the infiltration surface and consequently the recharge are limited.In ephemeral channels, water infiltrates into the permeable bed and bank of the channel, and is usually named transmission losses.However, most of the water fills up the unsaturated zone before reaching the groundwater (Gheith and Sultan, 2002;S ¸en, 2008).Transmission losses and recharge depend on several factors, such as underlying geology, temporal and spatial variability of flood events, and soil moisture characteristics.
Besides natural recharge to the groundwater, water can be stored (1) in a reservoir upstream of a dam, or (2) artificially infiltrated into an underground aquifer.Due to high evaporation rates from the surface in arid regions and the costs to construct storage dams, artificial recharge of water to an underground aquifer is often an efficient solution for water scarcity problems.Thus, underground storage Published by Copernicus Publications on behalf of the European Geosciences Union.
is becoming a major alternative for overcoming seasonal groundwater deficiency.
Artificial recharge is a method to balance and recover groundwater resources through floodwater spreading systems and injection wells (e.g.Bouwer, 1996).Recharge is a fundamental component of groundwater systems (Sanford, 2002).However, often adequate estimation of recharged water is difficult due to complex geophysical features and the large temporal and spatial variability of runoff (Rushton and Ward, 1979;Sophocleous, 1991;Arnold et al., 2000;Jyrkama et al., 2002;Flint et al., 2002).These issues are important to consider in order to assess the efficiency of artificial recharge systems and the reliability of estimated recharged water.Determining which of a wide variety of techniques that is likely to provide recharge estimates is often also difficult, and many factors need to be considered when choosing a method of quantifying recharge (Scanlon et al., 2002).In addition, an artificial recharge area is often constituted by both natural and artificial recharge that must be taken into account to adequately estimate the influence on the groundwater resources.
Many different methods have been proposed to adequately estimate groundwater recharge.Examples include water table fluctuation methods (Butterworth et al., 1999;Healy and Cook, 2002), Darcy's law (Butterworth et al., 1999;Healy and Cook, 2002), tracer techniques (Healy and Cook, 2002), mathematical models (Taylor and Howard, 1996) and/or a combination of several methods (Lin and Anderson, 2003;Sutanudjaja et al., 2011), e.g.mathematical model and tracer technique.However, each method has its own limitations considering available data, space and timescales, and range and reliability of recharge estimates (Zhu, 2000).
The scarcity of data in many arid regions, especially in the Middle East, has necessitated the use of combined mathematical models and field observations to estimate recharge.Mathematical groundwater models are used to simulate aquifer conditions, to estimate aquifer parameters, and to predict groundwater condition.In addition, as groundwater is essentially a hidden resource, studies on groundwater under both natural and artificial boundary conditions require modeling techniques (Scanlon et al., 2002).
The ability to use groundwater models to estimate recharge has been made easier by the development of inverse modeling techniques (Sanford, 2002).Many studies have used groundwater inverse modeling to investigate the effect of total recharge on groundwater storage (Poeter and Hill, 1997;Samper-Calvete and García-Vera, 1998;Flint et al., 2002;Sonnenborg et al., 2003;Dahan et al., 2004;Hendricks Franssen et al., 2009;Karlsen et al., 2012), but very few studies have been conducted on quantifying the contribution of different sources of recharge, including both natural and artificial recharge systems (Vázquez-Suñé et al., 2010).To address this issue, Vázquez-Suñé et al. (2010) concluded that the inverse modeling approach may prove a good tool for total recharge evaluation, but does not help in identifying the  contribution of each particular source to the total recharge.To the authors' knowledge, very few studies have aimed at quantifying the recharge contribution from different sources in arid and semiarid areas.This motivated us to use an inverse groundwater modeling approach to quantify the contribution of both natural and artificial recharge in an arid area located in southern Iran.
In view of the above, the main objective of the present study was, first, to apply inverse groundwater modeling to estimate the recharge rate of an artificial floodwater spreading (FWS) system and a natural river channel for both steady and unsteady conditions based on observed groundwater levels.Second, we wanted to compare the estimated recharged water through the artificial recharge system with the estimated recharged water through the natural ephemeral river channel.In this paper, we describe the conceptual groundwater recharge model and present a preliminary verification of the steady-state model.Partial objectives were also to analyze model precision, uncertainty/reliability, and model sensitivity.

Study area and observations
The Gareh-Bygone Plain (GBP) is located between 53 • 53 and 53 • 57 longitude and 28 • 35 and 28 • 41 latitude at an altitude ranging from 1125 m to 1185 m above mean sea level, and 190 km southeast of Shiraz city (Fig. 1).The study area is composed of an alluvial fan with an average thickness of 25 to 30 m on red clay bedrock and is mainly covered by sand deposit.This unconsolidated medium has created an unconfined aquifer with an area of 6000 ha.
The north to northeastern part of the GBP is delimited by impermeable layers of marl and sand stone that constitute the Agha-Jari formation of the Gar Mountain, created during the Mio-Pliocene period (Fig. 2).The eastern to southeastern part of the studied area is delimited by a series of discontinuous hills that approximately coincide with the Tchah-Qootch ephemeral river from the upper adjacent Tchah-Qootch subbasin (Fig. 2).This river and the ephemeral river coming from the upper Bisheh-Zard Basin comprise the main sources of surface inflow water into the GBP (Fig. 2).The 28-km long river is the main source of water input to the FWS system.However, the Tchah-Qootch River also contributes to the artificial aquifer recharge.These two rivers join in the lower southeastern part of the GBP and discharge into the Persian Gulf.According to aerial photographs of the studied area, there has been no shifting of the confluence or branching points of the existing ephemeral rivers during the study period.
Rainfall is the main provider of regional water input to the studied area, which is influenced by the Mediterranean synoptic system.Rainfall often starts in October and can continue to late April.The rest of the year is more or less dry.Annual average rainfall obtained from the newly established meteorological station (1995) in the GBP is about 212 mm and annual potential class A-pan evaporation for the same period is about 2555 mm.
In the GBP, groundwater levels have been monitored monthly by the Fasa District Water Organization since 1993.There are only four observation wells (P1, P2, P3, and P4) in the GBP, and these are not evenly distributed over the area (Fig. 2).Two newer wells (P5 and P6) were established in 2007 (Fig. 2).These wells were used to verify the steady-state model.

Floodwater spreading (FWS) system
FWS is the collection of runoff water from the upper basin(s) for irrigation and groundwater recharge.FWS is a technique where runoff from upland areas is collected and redistributed on a smaller area to artificially recharge the groundwater (Barksdale and Debuchananne, 1946).Artificial recharge by water spreading has been practiced at 36 multipurpose FWS sites in Iran since 1983.The technique is an inexpensive method for flood mitigation and artificial recharge of aquifers that results in a large economic return for relatively small investment (e.g.Ghayoumian et al., 2005).However, floodwa-ter in arid environments is usually characterized by a high load of suspended material.Hence, reduction in the stream water energy in the recharge zone results in accumulation of clay particles on the reservoir bottom.Therefore, floodwater spreading over an extensive area rather than in a small pond, injection well, or storage dam can be considered a preferred method.Along with this, an FWS system consists of several different basins, the first two basins act mainly as sedimentation pools and the rest act as infiltration pools.In principal, the investigated FWS system serves as sedimentation basins and infiltration ponds for the artificial recharge of groundwater and also as experimental plots for improvement, moving sand stabilization, flood mitigation, and reforestation (Kowsar, 1992;Arzani, 2010).The main objective of the system in Iran is to improve groundwater quantity and quality.
The FWS system in the studied area was first established in 1983 with an area of about 1400 ha and developed to 2000 ha in 1996.FWS is an improved water harvesting system consisting of six major components: (1) diversion dam, (2) conveyance channel, (3) earth embankment situated along contour lines, which consists of water gateway(s), (4) spreader channel, (5) basins which represent either a sediment or infiltration pool, and (6) outflow channel.
FWS system obtains floodwater directly from the ephemeral river by means of a diversion dam and a conveyance channel.Floodwater enters the first spreader channel and overflow spreads uniformly on the ground until it reaches a height of 25 cm, after which the water flows to downstream basins by one or several gateways.This process is repeated for all basins until excess water may be returned back into the river.The temporary 25-cm ponding in the basins allows water to infiltrate and groundwater to recharge.

Conceptual model and data
Groundwater modeling requires five components: data, conceptualization, simulation, calibration (Yang et al., 2010), and verification.The data required for the recharge model, to simulate natural and artificially recharged water in an unconfined aquifer, include: (1) geology and soil properties of the aquifer; (2) ground surface and bedrock elevation; (3) observation points; (4) borehole records; (5) boundary conditions; (6) location of recharge, e.g.river network, and discharge, e.g.pumping well, cells within the model; and (7) starting head which can be simulated within the first model run based on groundwater level at the observation points (Fig. 3).Moreover, the data (measured or estimated through inverse modeling) include: (1) specific storage; (2) specific yield; (3) horizontal and vertical hydraulic conductivity (based on the homogeneity concept both are assumed to be equal).Todd and Mays (2005)    of the key flow processes of the system.A conceptual model is a simplified representation of the groundwater system and the accuracy of a numerical model depends on how well the conceptual model, upon which it is based, captures the conditions within the aquifer (Shaki and Adeloye, 2007).Accordingly, a three-dimensional finite-difference approach of the groundwater model was employed using the MODFLOW-2000 software (Harbaugh et al., 2000).Moreover, PEST (Doherty et al., 2004), which is the automated parameter estimation module for MODFLOW, was used to optimize parameters for best agreement between simulated and observed groundwater elevations.The modeling period spanned 14 yr, from January 1993 to January 2007.
In terms of modeling discretization, if the grid size is large, hydraulic details contained in the region might not be considered, and if the grid size is too small, details are well exposed but with high computational load.Based on the finitedifference schematization, the studied area was discretized into a 250 m grid size, taking into account the distance between groundwater discharge sources, e.g.pumping wells, which enabled the model to directly calculate the impact of each well for the groundwater storage (Fig. 4).
The anisotropy factor in the MODFLOW-2000 is defined for each layer as the ratio of hydraulic conductivity along a column to hydraulic conductivity along a row (McDonald and Harbaugh, 1988;Restrepo et al., 1998).In reality, aquifer properties are generally heterogeneous with different degrees of variation (e.g.Yang et al., 2010).However, based on the alluvial fan properties, existing borehole records, and in order to simplify the modeling process, the model domain was discretized vertically as a single-layer, unconfined aquifer (Fig. 4), and the anisotropy factor was assumed to be one means by which the horizontal and vertical hydraulic conductivity are equal.For this, a single aquifer property value was associated with the aquifer layer (Yang et al., 2010).In addition, due to monthly records of the hydraulic head in the observation wells, the time-discretization within the model was assumed as one month.One crucial issue in groundwater model design is to determine the exact elevation of the ground surface and bedrock.Therefore, an existing 1/5000 scale topographic map was used to enter multiple elevation points into the model domain.Then, ground surface elevation between these points was obtained through inverse distance interpolation.The lowest and highest elevation points of the GBP are 1130 m (at the lower part of the plain) and 1184 m (at the foothill) above mean sea level, respectively.
As a consequence of the FWS system, the number of legal and illegal pumping wells in the area increased to over 120 wells in 1996 (Ghahari and Pakparvar, 2007), about 10 times the number before establishing the system.However, this number decreased to 85 in 2006.An inventory of existing pumping wells in the studied area revealed that most wells (with 10-15 m depth in 1995 or earlier) either were dried out or had a low water level due to excessive exploitation, especially during the dry period between 1997 and 2003.Extending well depths by drilling exposed the bedrock at 30-40 m depth.Thus, the depths of 85 existing pumping wells were used to determine the elevation of the bedrock to produce a bedrock contour map.The lowest and highest elevation points of the GBP's bedrock are 1090 m and 1155 m above mean sea level, respectively.
Besides gathering well pumping discharge data from the Fasa district water organization, a lot of time was spent measuring the discharge rate in each pumping well using a water discharge gauge (Jet ruler) and determining the exact location of the pumping wells using a GPS.Collected data (average discharge rate for each well) was used for the unsteady model simulation.Since there were no available data for agricultural return flow, 10 % of total discharge from the pumping wells were subtracted as a discharge input to the unsteady model (e.g.Jafari et al., 2012).However, it was modified during the calibration process of each unsteadystate model in order to achieve the best agreement between observations and simulations.

Boundary conditions
Based on geological maps, aerial photographs, and numerous field surveys of the studied area, three different types of boundary conditions were defined in the model design; (1) no-flow boundary was assigned to the northern part of the model area coinciding with the Gar Mountain (Fig. 2).Observed hydraulic head outside but close to the northwestern border of the model domain was used to assign the (2) general-head boundary condition at this border.In this case the groundwater was exchanged between the model domain and adjacent aquifer.Based on the topography and bedrock map of the GBP, the general trend for groundwater flow is from north to southwest.Therefore, a (3) timevariant specified head boundary was defined along the south and southwestern border of the GBP as a discharge/drain zone.Also, due to a few observed water levels in some of the pumping wells located in the southwest part of the GBP and outside the model domain, a very short time-variant specified head boundary condition was defined for this area as a recharge zone.
Due to the large hydraulic head at observation well number one (P1) as compared to other observation wells (Fig. 2), it was assumed that there is a direct connection between P1 and an external source.After analyzing satellite images and aerial photographs of the study area, the existing geological map was modified with a fault affecting the hydraulic head at P1 (Hashemi et al., 2012).

Steady-and unsteady-state groundwater model
In order to simplify the modeling approach, the conceptual model can be divided into steady-state and unsteadystate modeling.First, aquifer parameters such as hydraulic conductivity are estimated in the steady-state modeling then estimated parameters can be transferred to the transient or unsteady-state model (Middlemis et al., 2000) in order to estimate other aquifer parameters, such as storativity and recharge fluxes through inverse modeling.
In principal, if the model is run for steady-state conditions, the only input data observed is water level at each observation well during the steady period.However, steadystate conditions rarely occur during intense recharge or discharge periods.For the unsteady-state or transient condition, input-observed data are the water levels in observation wells at different time intervals and recharge and discharge of the groundwater system.
Steady flow implies that no change (no recharge or discharge) occurs with time (Graham and Tankersley, 1994).Consequently, when the hydraulic head variation for observed wells during successive months is at minimum, steady-state conditions can be assumed.Based on a modeling point of view, in the steady-state condition, the net input and output from each cell is set to zero, which means no water is added or taken from the internally stored water in a cell.Thus, the aquifer storage coefficient is not involved (Eqs. 1 and 2).In this case, the only unknown parameter which influences the groundwater flow is hydraulic conductivity that can be estimated in separate zones.Additionally, for steady-state simulations, there is no direct requirement for initial head because the time derivative (∂h/∂t) is removed from the flow equation (Eq.1).Therefore, a steady-state simulation is represented by a single stress period having a single time step with the storage term set to zero (Harbaugh, 2005).It should be mentioned that in the adequately chosen steady-state condition the system response to the boundary conditions in the aquifer allows for elaboration.
In principal, the three-dimensional saturated groundwater flow of constant density through porous earth material is described by the partial differential equation (Kresic, 2006): where K xx , K yy , and K zz are hydraulic conductivity along the x, y, and z coordinate axes assumed to be parallel to the major axes of hydraulic conductivity (L/T); h is the potentiometric head (L); W is a volumetric flux per unit volume representing sources and/or sinks of water, with W < 0.0 for flow out of the groundwater system, and W > 0.0 for flow into the system (T −1 ); S is the storativity of the porous material (L −1 ); and t is time (T).Equation (1) was developed assuming transient condition.However, the transient flow equation becomes the steadystate flow when the storage term and therefore, the right-hand side of the equation is set to zero (Eq.2) (Harbaugh, 2005;Kresic, 2006).
In view of the above, and according to the hydrogeological characteristics of the studied aquifer, a steady-state groundwater model was developed to simulate the groundwater flow, boundary condition, and estimation of hydraulic conductivity.To estimate the hydraulic conductivity, time steps corresponding to the steady-state flow were used for model calibration.For this purpose, ten different steady-state conditions during the entire model period were selected when the absolute difference in water level in all observation wells and between 2 successive months was at minimum or less  1).Also, observations showed that during the selected steady-state periods, the discharge through the pumping wells was at minimum and neither recharge through the river channel nor the FWS system occurred.Accordingly, the model was calibrated for each steady flow to estimate hydraulic conductivity and determine the appropriate boundary condition.Calibration for different independent time periods was done to achieve an acceptable confidence of accuracy of the estimated hydraulic conductivity.Then the average estimated hydraulic conductivity in each zone was transferred to the unsteady model to estimate specific yield and recharge rate.
In mathematical groundwater modeling unsteady-state conditions occur when the net flow into or out of a cell are not equal.In this case, besides the estimated hydraulic conductivity, it is possible to estimate the storage coefficient and recharge rate.To simplify the calculations, the time period with no recharge into the aquifer should be chosen in order to minimize the number of unknown parameters.Therefore, firstly, three different time periods of 6, 8, and 12 months' duration with no recharge were selected.Then the model was run and calibrated for these three transient periods to estimate specific yield in each zone.Finally, the values were transferred to the next step of transient model to estimate the recharge rate.
In transient flow, simulations can be mixed transient or steady-state.Thus, a simulation can start with a steadystate stress period and continue with one or more transient stress periods (Harbaugh, 2005).Accordingly, ten different transient models were assigned and calibrated starting from each steady-state and ending with the next one.Moreover, in MODFLOW-2000 the input parameters are constant during the transient modeling.Consequently, by dividing the entire model period into different sub-transient models, one can deal with time variable parameters within the entire model period but constant parameters in each sub-transient model.Also, to achieve best agreement between observed and simulated values, the average estimated hydraulic conductivities from the steady models were modified slightly in each subtransient model.

Calibration and verification
Monthly observations from four wells (P1, P2, P3, and P4) located within the GBP during 14 yr between 1993 and 2007 were used for calibration in ten different steady-state and transient periods (Fig. 2).Monthly observations from two newer wells (P5 and P6) during 2007 and 2009 were used for model verification in the case of steady-state condition (Fig. 2).It should be mentioned that in the verification process normally the data from the same wells are used and this is an advantage if the model can be verified by the data from more wells.
Due to the limited number of observation wells but assuming the homogeneity of the geologic formation in the GBP, the Thiessen method was used to find representative zones around each observation well.Therefore, the entire area was divided into four zones, and the model was calibrated with observed hydraulic head at the observation well in each zone in both steady and transient condition.In the case of verification, the model was verified using both old and new observation wells (six wells) but still four hydraulic zones in one more steady-state period (during the period 2008).

Model sensitivity and uncertainty assessment method
Estimated parameters are subject to a degree of uncertainty that can be associated with the conceptual model or with data and parameters for the various components of the model (Zhu, 2000;Walker et al., 2003;Arnold et al., 2009;Hildebrandt et al., 2010;Rojas et al., 2010).In conceptual models, uncertainty estimation is dealt with calibration and estimation procedures (Ratto et al., 2007;Dams et al., 2008;Rojas et al., 2010).
Calibration may be affected by uncertainty due to weakness in determining the exact parameters.This weakness could be due to either inadequate or to inappropriate spatial and temporal distribution of the observed data and boundary conditions.Uncertainty also can be due to errors in observations.
In general, to determine the degree of uncertainty of the estimated parameters, sensitivity analysis and confidence interval are used.Uncertainty analysis includes a number of techniques for determining the reliability of model estimations/predictions (Hill and Tiedeman, 2007).For this reason, confidence intervals are used to indicate the reliability of estimated parameters.Sensitivity analysis can be done for both parameters/unknown and observed/known data within a model.Determining the sensitivity of a specific parameter to the observed data means the change in the estimated parameter corresponding to a unit change in the observed data.If the desired parameter value is zero, this parameter is not sensitive to the observed data.Parameters that generally are insensitive to observed data can be omitted from the model or considered constant.
Based on the contents of the Jacobian matrix (Hill and Tiedeman, 2007), PEST calculates a figure related to the sensitivity of each parameter with respect to all observations.The composite sensitivity of parameter i is defined as (Doherty et al., 2004) where J is the Jacobian matrix (also called the sensitivity matrix); Q is the cofactor matrix, in most instances the latter will be a diagonal matrix whose elements are comprised of the squared observation weights; m is the number of observations with non-zero weights.Thus, the composite sensitivity of the ith parameter is the normalized (with respect to the number of observations) St. 9 St. 10 magnitude of the column of the Jacobian matrix pertaining to that parameter, with each element of that column multiplied by the weight pertaining to respective observation (Doherty et al., 2004).
In view of the above, sensitivity analysis was done for the parameters at all modeling levels in order to find the best estimation.Furthermore, the hydraulic head estimation was followed by a sensitivity analysis aiming to identify the parameters, which have a large influence on the modeling results.
PEST also calculates 95 % confidence limits for the adjustable parameters.The presentation of 95 % confidence limits provides a useful means of comparing the certainty of different parameter values as estimated by PEST (Doherty et al., 2004).It should be mentioned that confidence limits can only be obtained if the covariance matrix has been calculated.Hence, in case of any over-or underparameterization, the confidence interval is not produced (see PEST user manual).

Results and discussion
The main objective of the present study was to estimate the recharged water from the surface runoff via both the natural ephemeral river channel and artificial recharge system through an inverse groundwater modeling approach.Therefore, the groundwater model was run for both steady-state and unsteady-state conditions to simulate the groundwater flow and estimate aquifer hydraulic parameters for the period 1993-2007.The estimated parameters included horizontal hydraulic conductivity (K h ), specific yield (S y ), and recharge rate (RCH).

Steady-state model results
In the inverse modeling technique, aquifer parameters such as hydraulic conductivity and recharge rate are estimated simultaneously.Hence, reliable estimation of recharge rate depends on adequate estimation of the hydraulic conductivity.With regard to this, ten different steady-state simulation periods were modeled for which the hydraulic head difference between successive months was negligible.These ten steadystate simulations were calibrated, using the PEST module, regarding their K h values.Then the model was verified using a new steady-state period in 2008 (Hashemi et al., 2012).The summary results of the model calibrations in ten steady-states and its verification are presented in Table 2.The results show that the calibrated and verified K h value are quite close, however, the residual and standard deviation are somewhat larger for the verified K h .
As a result of steady-state simulation, the general trend of groundwater flow is from north to south/southwest.The steady-state groundwater flow with no recharge from floodwater is greatly influenced by the water coming from the upper catchment via the fault.The average estimated K h in the GBP was nearly 0.1 m day −1 , which is in the range of values of K h (0.001 and 1 m day −1 ) for the alluvial fan (e.g.Freeze and Cherry, 1979

Unsteady-state model results
As previously described, simulation and estimation of the hydrodynamic aquifer parameters were done in three different model steps.To estimate the S y in each zone as a second step, the model was calibrated in three transient period intervals with no recharge but active pumping wells.Besides model calibration based on monthly observed groundwater levels, groundwater contour maps and direction of groundwater flow were also simulated.Groundwater contour maps (Fig. 5) were developed at each model step, taking into account the initial groundwater level in the beginning of the simulation period and water abstraction through pumping wells during each simulation period.The estimated S y in each zone ranged from 0.008 to 0.10 with an average of 0.045.The average value of S y was then transferred to the next transient interval to estimate the recharge rate.During these periods there was no recharge from the surface water and water was exploited through the pumping wells.As obtained from simulations, the groundwater flow direction in some parts of the aquifer is moving toward pumping wells  to the extent that in some areas a sharp cone of depression around the wells is formed (Fig. 5). Figure 5 shows that in the transient periods with no surface water recharge and with high groundwater exploitation rate, saline water intrusion to the groundwater from Shur River of Jahrom (saline water) (Figs. 2 and 5) may be a problem.The base flow of this river is quite saline with electrical conductivity ranging from 0.6 to 4.5 S m −1 (Kowsar and Pakparvar, 2003).Therefore, it appears that the groundwater flow direction was changed during these dry periods and was from the Shur River of Jahrom toward the central and southern parts of the aquifer (Fig. 5).
In order to estimate the recharge rate through both the river channel and FWS system, the area was divided into three different recharge zones.Artificial recharge zones 1 and 2 (AR 1 and AR 2) correspond to the FWS system, and zone 3 represents the natural ephemeral river (ER) channel (Fig. 6).Simulations and recharge estimations were conducted for the entire model period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) taking into account the establishment dates of different recharge zones.AR 1, including zones 1a and 1b (Fig. 6), only receives diverted floodwater from the Bisheh-Zard ephemeral river and was developed from 420 ha in 1983 to 1020 ha in 1996.AR 2 receives diverted floodwater from both the Bisheh-Zard and Tchah-Qootch Ephemeral Rivers covering an area of about 970 ha.It should be noted that the Tchah-Qootch River is located at the southern border of the GBP and due to the hydraulic gradient of the plain, the natural recharge of this river does not contribute to aquifer recharge.In addition, it was assumed all recharge zones received surface water at the same time and flooding duration was equal in each flood event considering flood date but not magnitude of flood.To check model fit, simulated hydraulic heads were compared to observed heads for the entire model period (Fig. 7).The recharge rate was estimated when the best agreement was obtained between observed and simulated hydraulic head using parameter estimation software (PEST).The average observed and simulated hydraulic head for all observation wells is shown in Fig. 7.In addition, Table 3 shows the statistics for observed and simulated values for all observation wells in all unsteady-state model periods between 1993 and 2007.
It should be noted that in order to balance inflow and outflow to the groundwater storage, the model adds or withdraws some water into the storage term but does not increase the recharge rate to the considered recharge zone.This makes the model quite precise and the residuals are almost the same in the case of modeling one, two, and all recharge zones.In addition, due to the monthly observation data, the duration of recharge for each event was assumed to be within one month.

Uncertainty and sensitivity analysis
Figure 8 shows the estimated recharge and 95 % confidence limits of all estimated values for the ten different simultaneous transient periods during 1993 and 2007.In these simulations, the estimated values tend to increase as the size of the intervals increase and vice versa.For instance, for AR 2 all intervals are quite small and, subsequently, the estimated recharge is relatively small in comparison to other recharge zones.The result shows that the uncertainty and confidence intervals of ER (river channel zone) is the highest.Also, the estimated recharge in this zone is relatively high for some of the model periods.These high values with high uncertainties match well with time periods representing extreme rainfall events.Thus, the river channel recharge estimation represents relatively more uncertainty, which can be due either to less model sensitivity for parameters of ER or to major influence of extreme flood events on the estimated recharge in the channel.Also possible explanations may be the higher infiltration rate and percolation velocities under the river channel as compared to the flooding zones.This is a localized phe-   nomenon that is expected to act during periods of time and therefore will be reflected in a larger uncertainty.In general, the recharge rate for the plain area (artificial recharge area) was much lower than that of the river bed which was close to zero as shown in Fig. 8.However, considering the large infiltration surface of the FWS system, it results in a large volume of recharged water.
Sensitivity analysis was conducted to evaluate the effect of parameters on groundwater modeling and to identify which parameters have the most effect on recharge.Figure 9 shows the model sensitivity for all recharge areas (AR 1, AR 2, and ER) for the ten model periods.In principal, the model is sensitive to all areas for all model periods.The model, in general, is more sensitive to AR 2 and less sensitive to ER for all recharge periods.Since AR 2 receives floodwater from two different ephemeral rivers, the model is quite sensitive to this zone, and this sensitivity reaches a maximum in period P2, P3, and P9 due to extreme rainfall events (Fig. 10).In contrast, the model is much less sensitive to ER, which is the river zone, and that can be due to the large uncertainty for this zone, as discussed above.This shows the significant influence of AR 2 on the model results.

Recharge contribution to groundwater from both artificial and natural recharge area
S ¸en (2008) stated that groundwater recharge has a random behavior depending on the sporadic, irregular, chaotic, and complex features of storm rainfall occurrences, geological  composition, and geomorphologic features.To quantify the amount of recharged water for each recharge zone for the 14-year modeling period, the groundwater simulation and recharge estimate was done in three different steps: (a) modeling AR 1, AR 2 and the ER, (b) modeling only AR 1 and AR 2, and (c) modeling only AR 2 (Fig. 6). Figure 10 shows the estimated recharged water volume between 1993 and 2007.
The figure shows that before 1996, more water was recharged through the river channel rather than the FWS systems.So, the difference between total recharged water (including AR 1 (a, b), AR 2, and ER) and recharged water through the FWS systems (AR 1 (a, b) and AR 2) is significant for recharge cases between 1993 and 1996 (Fig. 10).In some periods with extreme precipitation events this difference reached a maximum, e.g.2003 and 2005.As is seen in Fig. 10, there is, however, no significant difference in the precipitation for either the 1993-1997 or the 2003-2006 periods, but the recharge in 2003-2006 was about double that in 1993-1997 with the exception of one event in early 1993.However, it is noted that from 1997, after the extension of FWS system area 1 (AR1) from 420 to 1020 ha, water inflow to the system increased and more water was recharged through the FWS system.In case of clogging of the system, the recharge difference between the AR (1 and 2) and ephemeral river (ER) should be at maximum, which is not the case, because more water recharged through the AR (1 and 2) e.g. during 2006.
During extreme events, since the FWS system is not capable of handling more water than its capacity and perhaps because of a complex interaction between different recharge zones, and duration of stream flow in the river channel, a greater volume of floodwater was transferred via the river channel and recharge took place more through the river bed and bank.However, in the year 2004, despite the extreme rainfall event, in terms of duration, the difference between total recharged water and recharged water through the FWS systems was small.This can be due to the small flood resulting from low precipitation intensity or soil-moisture characteristics of the river channel which was not taken into account in the model.
Less water recharged through AR 2 (area of 970 ha).This is probably due to poor vegetation cover (observed in-field,   Fig. 11) or clogging in this zone and consequently less recharge occurred as compared to AR 1 (a, b).In order to stabilize and increase the infiltration rate in the FWS systems, indigenous trees adapted to this environment have been tested, but as their growth is too slow, two species belonging to the exotic genera (Eucalyptus and Acacia) have been experimentally planted, mainly in AR 1, along with three indigenous trees and bushes (Khanmirzaei et al., 2011).The trees were planted on plots of 4 × 4 m spacing covering about 15 % of AR 1 and also behind each embankment in the entire area of AR 1 and part of the area of AR 2 in 1986 (Fig. 11).
The Agha-Jari formation covers most of the Tchah-Qootch Basin and represents early to middle Miocene and greater sediment volumes, as compared with the Bisheh-Zard Basin, being transported to the GBP during flood periods.Since AR 2 receives water from both Tchah-Qootch and Bisheh-Zard Ephemeral Rivers, sediments also are transported that can be a cause of clogging and infiltration reduction in this recharge zone.To address this issue, the research carried out by Mohammadnia and Kowsar (2003) shows increase in the concentration of three different types of clay at the depth of 7.5 m after 13 yr operation of the system, which proves the gradual impermeability and eventual clogging of the vadose zone by translocation of extremely small particles.Also, research reveals that the average infiltration rate of an afforested system was 9.3 cm h −1 as opposed to 3.8 cm h −1 for a treeless system, and 7.7 cm h −1 for the area outside of the FWS systems (Kowsar, 2005).It is well known that plant roots, together with organic material, increase infiltration rates through physical and biological processes.Consequently, reforestation of areas also means that the effects  of siltation on permeability to some extent are balancing the recharge rates.
As a result, after development in 1996, more water recharged through the artificial recharge system rather than through the river channel, and AR 1 (a, b) is the main recharge zone among all FWS zones.

Recharge contribution to groundwater from the river channel
Figure 12 shows the cumulative recharged volume separated into two parts, namely with and without the river channel.In 1995 and the beginning of 1996, due to the development of FWS system, the inflow channel to some of the FWS systems had been blocked and more floodwater passed the area via the river channel.Therefore, more water recharged through the river channel.Also, due to an extreme event in 2003 the modeling result shows a small rise in groundwater recharge from the river channel in comparison with the FWS systems, and this rise reached a maximum in 2005 (Fig. 12) due to the occurrence of major floods in that year.
In order to determine the contribution of each recharge zone in the GBP, the percentage of recharge in each recharge zone, considering all years of extreme events and all years without extreme events, is listed in Table 4. Based on the second row of Table 4, the river channel is the main recharge zone in the studied area during three different years of extreme rainfall events mentioned above.According to Table 4, the recharge contribution from AR 1 (a, b), AR 2, and ER were 31 %, 10 %, and 59 %, respectively.There is a large difference between the cumulative natural and artificial recharge by the end of 2006, which resulted from a few extreme events in 1995, 1996, and 2005 (Figs. 10 and 12) (Figs. 10 and 12).By subtracting these extreme events and related recharged water through the river channel, the contribution of the river channel decreased from 59 % to 20 %, and the contribution of the FWS systems increased from 41 % to 80 % (Table 4).
It appears that in a normal year without extreme events, the recharge contribution to the groundwater storage from the river channel is about 20 %, while the recharge contribution from the FWS systems is about 80 %; this means that the FWS system is the main recharge zone in the studied area.In contrast, the river channel area is just a minor portion of the plain (Alencoão and Pacheco, 2006) that only covers about 10 to 15 % of the total recharge area.These results clearly show the value of artificial recharge systems in arid and semiarid areas.Therefore, extension and development of FWS systems in the studied area would be a main and permanent solution to increase the floodwater recharge in case of any extreme or small flash flood event.Bull and Kirkby (2002) stated that ephemeral river channels exhibit an inverse relationship between the magnitude and frequency of events.By expanding the FWS systems there is a possibility to divert more floodwater into the artificial recharge system and increase recharge even in a small flood event, while the recharge through the river channel increases only for large flood events.Therefore, larger infiltration surfaces result in less water lost via river channels even for small flash flood events.

Summary
The main objective of the FWS system is to increase the recharge by increasing the area of flooding.Consequently, the increased area governs the recharge volume.Assuming similar recharge per area for both channel and artificial area means that the volume recharged by far is from the artificial recharge areas.Similarly, if the areas are maintained then recharge must be greater in the recharge areas.However, even if there is no doubt that many of the stream channels may have a high capacity of recharge, and as mentioned before, then lower magnitude flood events generally do not result in complete channel submersion.Therefore, the infiltration and consequently the recharge are limited.Thus, recharge through the river channel is a natural phenomenon and most probably with a high infiltration rate but with limited and often insufficient volumes for domestic use due to the limited infiltration surface area.
There are two common ways for providing water for the rural community (1) storage dams about which several studies have indicated their inefficiency due to the high potential for evaporation, silting, filling up of the reservoir with sediments, cost of construction, environmental impact, and (2) artificial recharge.However, arid zone flash floods carry a high load of suspended materials that may cause clogging in artificial recharge systems.Studies have demonstrated that many floodwater harvesting systems in arid areas were abandoned after several years of operation due to clogging, but, for most of them the author(s) did not give any reason for site selection, inspection, and maintenance of the system.The FWS system constitutes a low-cost passive technology that may give high return in terms of increased groundwater supply.Owing to that, the FWS system, as well as all water/hydraulic infrastructure projects and techniques need regular inspection and maintenance.
One reason among many for abandoning FWS systems is improper site selection.The water harvesting method applied strongly depends on local conditions and includes widely differing practices, such as bunding, pitting, micro-catchment water harvesting, floodwater and groundwater harvesting (e.g.Critchley et al., 1991;Prinz, 1996).In order to find the most appropriate technique, the area must be studied in terms of soil properties, geology of the upstream catchment, topography, floodwater quantity and quality, climatology, hydrology, etc., before establishing such systems.In addition, since the FWS system requires an extensive area, the sediment may spread over a large area with relatively limited thickness.However, this property creates limitations for site selection, which tends to be site-specific.This is an area for future research.
Another discussion point is the unspecified budget at the government level for such projects in many arid/semiarid countries.Unfortunately, once the system has been established, no or insufficient budget may be allocated to maintain the system.Our studies show that after the first year of operation about 20 % of the initial cost of the project should be allocated for maintenance and from then, 10 % is required for an annual maintenance.However, this may increase after each significant flood.Hence, continuous maintenance is a basic requirement to have a fully functional system.Sediments, just like other problems, can be managed, for example silt can be removed every few years.
In summary, clogging and siltation in all artificial recharge systems as well as all water structures are inevitable and result in less efficiency.According to this, the studied system in the beginning was assigned a 12 to 15 yr lifetime.But our studies and particularly the one presented here show that the system is still working well after about 30 yr.Hence, with regular inspection and maintenance of such systems, their lifetime and efficiency will be much longer, and they will yield a high return in terms of small yearly investments.

Conclusions
Artificial recharge plays a significant role in improving the groundwater resources in the arid and semiarid Middle East.In order to understand the function and efficiency of artificial recharge systems, the contribution of such systems in total recharge needs to be quantified using a proper technique.
There are a great many advantages in modeling long-term natural and artificial recharge for complex natural systems to overcome some of the problems associated with the lack of data in Middle Eastern countries.However, there is still a lack of understanding of the recharge processes for both natural and artificial recharge zones in arid areas.Monitoring the natural and artificial recharge systems combined with modeling techniques is an appropriate way to better estimate arid and semiarid region recharge.The methodology presented in this study was based on groundwater inverse modeling to estimate the contribution of natural and artificial recharge in a complex system.The fundamental benefit of inverse groundwater modeling is its ability to automatically calculate parameter values that produce the best fit between observed and simulated hydraulic heads and flows (e.g.Poeter and Hill, 1997).
The only available groundwater data for this study were observed water levels for 6 monitoring wells (the oldest installed around 1993).For this reason, an inverse modeling approach was developed in which a confidence interval of the estimated parameters and sensitivity to the observations and parameters carefully were checked in all model steps in order to decrease the level of uncertainty.Our aim was to provide an approach capable of estimating the recharge parameters in such complex systems even with the scarcity of available data.
The modeling results show that, assuming the observation wells are representative of the behavior of the studied area's aquifer, in a normal year without an extreme event the contribution of natural recharge to the groundwater storage from the river channel is about 20 % and the contribution of artificial recharge from the FWS systems is about 80 %.Therefore, the FWS system is the main source of recharge in the studied area.
There are very few studies and inadequate information about the recharge estimates through the ephemeral channels where there is an extra effect on the groundwater system by an artificial recharge system.But the methodology presented here can be used as representative for such complex systems in estimating recharge and quantifying the contribution of different recharge sources in an area.In the arid Middle East there are numerous newly developed artificial (as well as traditional) recharge techniques.The model approach developed in this paper could be used to assess the improvement in recharge by any kind of modification and a way to optimize recharge systems.
Finally, in order to reduce the uncertainty of the estimated recharge, water movement through the unsaturated zone under both river bed and flood spreading area needs to be studied.For this reason, the research by Dahan et al. (2003Dahan et al. ( , 2009) ) can be a preferred method to monitor the water content in vadose zone.

1Fig. 1 .
Fig. 1.Map of Iran and location of study area.

Fig. 1 .
Fig. 1.Map of Iran and location of study area.

Fig. 2 .
Fig. 2. Floodwater spreading system, river network, and observation well distribution in the Gareh-Bygone Plain.Fig. 2. Floodwater spreading system, river network, and observation well distribution in the Gareh-Bygone Plain.

4Fig. 4 .
Fig. 4. Spatial discretization of the modeled area and the distribution of groundwater system (meters above mean sea level).

Fig. 4 .
Fig. 4. Spatial discretization of the modeled area and the distribution of groundwater system (meters above mean sea level).

Fig. 5 .Fig. 5 .
Fig. 5. Average groundwater contour map (meter above mean sea le April 2005 and March 2006, with active pumping wells but with no s Fig. 5. Average groundwater contour map (meter above mean sea level) of the GBP during 12 months between April 2005 and March 2006, with active pumping wells but with no surface recharge.

Fig. 6 .
Fig. 6.Location of recharge zones in the GBP, artificial recharge (AR) 1a with 420 ha area established in 1983, AR 1b with 600 ha area established in 1996, AR 2 with 970 ha area established in 1983, and ephemeral river channel (ER).

Fig. 6 .
Fig. 6.Location of recharge zones in the GBP, artificial recharge (AR) 1a with 420 ha area established in 1983; AR 1b with 600 ha area established in 1996; AR 2 with 970 ha area established in 1983; and ephemeral river channel (ER).

7Fig. 7 .
Fig. 7.The average observed and simulated groundwater level (GWL) for all observation wells (due to missing data in two observation wells during 2000 to 2001, no simulation was done for this period).

Fig. 7 .
Fig. 7.The average observed and simulated groundwater level (GWL) for all observation wells (due to missing data in two observation wells during 2000 to 2001, no simulation was done for this period).

Fig. 8 .
Fig. 8. Simultaneous confidence limits (upper and lower limits) on recharge rate estimated by 10 different unsteady-state groundwater flow model of three different recharge zones (AR 1, AR 2, and ER).

Fig. 8 .
Fig. 8. Simultaneous confidence limits (upper and lower limits) on recharge rate estimated by 10 different unsteady-state groundwater flow model of three different recharge zones (AR 1, AR 2, and ER).

Fig. 10 .
Fig. 10.Rainfall (right axis) and recharged water volume estimated by groundwater model (left axis) into the groundwater storage through both FWS systems (artificial recharge, AR 1 (a, b) and AR 2) and ephemeral river channel (natural recharge, ER) for the years between 1993 and 2007. 11

Fig. 12 .
Fig. 12. Cumulative recharged water estimated by groundwater model from the FWS system zone 2 (AR 2), all FWS systems (AR 1 (a, b) and AR 2) and all FWS zones plus ephemeral river channel (AR 1 (a, b), AR 2, and ER) for the entire model period.

Fig. 12 .
Fig. 12. Cumulative recharged water estimated by groundwater model from the FWS system zone 2 (AR 2), all FWS systems (AR 1 (a, b) and AR 2) and all FWS zones plus ephemeral river channel (AR 1 (a, b), AR 2, and ER) for the entire model period.

Table 1 .
Steady-state conditions for all observation wells.Horizontal lines indicate the assumed steady-state (St. 1 to St. 10) periods for two successive months between 1993 and 2007 in model calibration.

Table 1 .
Steady-state conditions for all observation wells.Horizontal lines indicate the assumed steady-state 927 (St.1 to St. 10) periods for two successive months between 1993 and 2007 in model calibration.

Table 2 .
(Hashemi et al., 2012)r ten steady-state modeled periods between 1993 and 2007 and verification result for 2008.In model verification, two new observation wells were used(Hashemi et al., 2012).

Table 3 .
Statistics for observed and simulated values for all observation wells (P1 to P4) in all unsteady-state model periods between 1993 and 2007.

Table 4 .
Contribution of artificial recharge (AR 1 (a, b) and AR 2) and natural recharge (ER) zones in groundwater recharge in the period from 1993 to 2007 with and without consideration of extreme events.