Assessment of Spatiotemporal Groundwater Recharge Distri-bution Using SWAT-MODFLOW Model and Transient Water Fluctuation Method

1. Introduction

Sustainable development and management of water resources are vital in watersheds that include urban and agricultural land use due to the rising demand for water supply, climate change, and improper management of water resources [1,2,3,4]. An Accurate understanding of the hydrological water cycle and its kinship system is mandatory to sustain and manage the water resources of the watershed [5,6]. It influences watershed recharge-discharge [6,7], water supply [8], and water quality. Further, a proper method for estimating water balance components is important [5,6,9]. Hence, this paper emphasis on assessing the groundwater recharge of the study watershed.

Preceding studies used several techniques to estimate groundwater recharge, which can be categorized into physical [10,11,12], tracer [13,14,15], and numerical methods [2,16,17]. Nonetheless, each method has its limitation while assessing recharge [11]. The numerical approach has become increasingly popular among researchers for evaluating groundwater recharge, which involves analyzing either surface or aquifer conditions, or both [18,19].

The Soil and Water Assessment Tool (SWAT) is usually applied to simulate surface water flow, while Modular Three-Dimensional Finite-Difference Ground-Water Flow Model (MODFLOW) is frequently employed for assessing unconfined groundwater flow. SWAT [20,21] has been extensively applied to determine surface and groundwater quality and quantity across various catchment scales. Further, it predicts how land use/cover and climate change influence designated catchment [9,18,22]. SWAT model have been applied in various hydrological watersheds for wide ranges of scale and multi-senior environmental studies [23,24,25]. For instance, SWAT model was applied in some parts of South Korea watersheds for water quality assessment [26], water balance analysis [27], and sustainability of water resources [28,29]. SWAT model places a strong emphasis on surface processes, but its ability to simulate groundwater flow is inadequate by its lumped nature [30]. Various strategies have been proposed to address the limitations of the original SWAT groundwater module, including altering the code of the module [31], implementing a more intricate groundwater storage module [32], or substituting the non-shallow aquifer with the grid cell aquifer [33]. Though, limited to simulating groundwater flow accurately due to semi-distribution approach. MODFLOW model is a suitable option for groundwater flow simulation since it is distributed and physically based, which was the limitation of SWAT model. This model considers various boundary conditions to simulate groundwater flow processes; such package includes Recharge, River, Reservoir, Well, and others [34]. It was also applied to investigate the influence of sub-surface water abstraction on surface water resources [35,36]. However, it is inadequate in hypothesizing surface water flow [9,18]. Hence, coupled SWAT-MODFLOW approach is critical to overcoming the drawbacks of models and for a better representation of the hydrological system [37].

Integrating SWAT and MODFLOW models was initially introduced by Sophocleous et al. [38] and has been upgraded while comparing the model performance with SWAT model [18,39,40]. In the integrated model framework by Bailey et al. [41], SWAT model focuses on simulating the surface hydrological processes, while MODFLOW-NWT (Newton-Raphson formulation for MODFLOW) [42], which is capable of solving unconfined aquifer sub-surface water flow problems, simulates groundwater flow processes. According to Wang et al. [43], preceding studies by the SWAT-MODFLOW model can be classified into non-situational and situational simulations. Non-situational state studies are primarily focused on streamflow simulation, groundwater discharge, surface water-groundwater interaction, and regional water balance [39,41,44,45,46,47], while situational scenario simulation studies using SWAT-MODFLOW consider the influences of groundwater abstraction on the surface water resources [9,48,49].

As specified by Ntona et al.[34], SWAT-MODFLOW is the most applied model for surface-groundwater analysis from the period 1992 to 2020 worldwide and it continues to be the keystone for recent studies. For instance, a study by Mosase et al. [50] estimates the spatio-temporal distribution of recharge, water table level, and interaction between surface and sub-surface water in the Limpopo River Basin, Africa. Sophocleous and Perkins [51] use SWAT-MODFLOW to three basins in Kansas, USA watershed to explore abstraction impacts on streamflow and water table levels. They also analyze the outcome of various scenarios, including climate change on water resources. A study by Gao et al. [39] analyzed spatiotemporal variations of the surface water resources applying integrated model in the River catchment, USA. Taie Semiromi and Koch [18] utilize the integrated model to inspect surface and subsurface water interaction in an agricultural watershed in Iran. They point out that climate change alone does not create future water scarcity but rather the overutilization of the groundwater source in the study of the agricultural watershed. Chunn et al. [52] apply the same model to evaluate consequence of climate variation and groundwater abstractions on the surface water-groundwater interaction in western Canada. They pinpoint that in the reduction of river flow and groundwater level, the overuse of groundwater contributes a high role, as stated by [45,53]. In recent hydrological studies, parametrization tools have been applied for SWAT-MODFLOW model, which reduces the uncertainty of parameters applied for the calibration of the model [30,54]. For instance, Liu et al. [30] applied PEST (parameter estimation by sequential testing) as calibration tools for integrated model and evaluate the streamflow response to groundwater abstractions for the Uggerby River Catchment in northern Denmark. SWAT-MODFLOW has a significant advantage in simulating interconnected surface-groundwater scenarios; such as in the case of streamflow evaluation, it showed high precision [9,30,41,46]. However, numerical models like SWAT-MODFLOW always have uncertainty [55]. Environmental models, especially those at the regional scale often suffer from significant uncertainty due to the large number of unknown factors. The source can be from model structure, input database, or parameterization options [43]. For example, Guevara-Ochoa et al. [44] mentioned SWAT model transfers a higher recharge amount to MODFLOW during the wet season since SWAT-MODFLOW does not have a module to associate the alteration among groundwater and soil saturation. While there are many other recharge estimation methods, water table fluctuation method is optioned as additional approach.

Among the many methods for estimating groundwater recharge, water table fluctuation (WTF) is a straightforward, accurate, and simple method [11,12,56]. The method required observed groundwater level data series and aquifer parameters, which can be determined using various methods [11,12,57]. To reduce subjectivity and increase efficiency several modifications have been taken to this method [58,59,60]. One of the approaches was to link the WTF method with a hydrological model, such as SWAT model, which helps to overcome local representation of the method. For instance, Chung et al. [59] estimated recharge by linking SWAT model with TWTFM in Hancheon watershed, Jeju Island, South Korea. Moreover, deep percolation (recharge from SWAT) was used as input for TWTFM to determine aquifer parameters (reaction factor and specific yield) for TWTFM. Applying an alternative approach to estimating groundwater recharge is recommended since it assures results [12,56,61,62]. This study aims to apply SWAT-MODFLOW model and TWTFM method to assess the groundwater recharge in our study area. As per our knowledge, no research has been reported to use coupled model and TWTFM to evaluate groundwater recharge and it will be presented for the first time in our study.

In this study, we proposed numerical and empirical methods to evaluate the spatiotemporal distribution of groundwater recharge. The spatial and temporal distribution of groundwater recharge using SWAT-MODFLOW model is evaluated for the study watershed. Further, we approach recharge estimation using TWTFM as an alternative option and evaluate the output with numerical model approach. Seasonal variation of groundwater recharge of the region is presented. Furthermore, the water balance segments are acquired and discussed.

2. Materials and Methods

2.1. Description of Study Region

Anyang watershed is located in northwestern South Korea which is displayed in Figure 1. The study watershed covers an approximately 137 km² area and drains into the Han River, one of the longest rivers in the Korean Peninsula. The altitude ranges from 11–59l m above mean sea level and an average of 116 m for the region. Between 2002 and 2018, the mean annual precipitation was 1266.8 mm, and the yearly averages for the lowest and peak daily temperatures were 8.5 and 17.5 ℃, respectively. For a while, the streams of the study region are gauged that started before 2002. During the summer period, the streamflow gets its maximum value.

Anyang watershed land use/land cover consist of most forest and urban, with forest enclosing approximately 50% of the whole area and urban-industrial areas covering 43%. The remaining 7% is comprised of pasture, water bodies, wetlands, and agricultural land use. The land use was classified into 14 classes by the SWAT model, as represented in Figure 2.

The soil types and hydrological soil groups of the study are displayed in Figure 3. In the study area, most soil types have four layers, the highest depth varies between 60 and 2000 mm. Hydrological soil group “A” dominates the region with 67.1%, which characterizes the soil texture of sand and sandy loam with low runoff potential. Then followed by soils that have a hydrological group of “B” with 27% coverage in the region, which encompasses loam and silt loam soils with moderate infiltration potential characteristics. Soil hydrological group “C” includes sandy clay loam soil, soil group “D” comprises clay, clay loam, and sandy clay soils. Hydrological soil group “C” comprises approximately 3.2% and “D” 2.7%. Both hydrological soil groups “C” and “D” are characterized by low seepage and a high potential for surface runoff. They mainly occupy the lowland zone, including wetlands and water bodies.

The information of Soil types and LULC of the study watershed were obtained from the National Institute of Agricultural Sciences and the Ministry of Environment of Korea, respectively. Daily streamflow data and weather data were provided by Water Resources Management Information System [63] and Metrological Agency Weather Data Service [64] respectively.

Figure 3. Anyang watershed soil types: (a) Soil classes applied in SWAT model (the whole information can be retrieved from [65] website); (b) hydrological soil groups.

Figure 3. Anyang watershed soil types: (a) Soil classes applied in SWAT model (the whole information can be retrieved from [65] website); (b) hydrological soil groups.

2.2. Evaluation of Recharge Using Integrated Model and Transient Water Table Fluctuation Method

2.2.1. SWAT-MODFLOW Model

SWAT model ability to address groundwater flow is restricted by its aggregated nature, while MODFLOW struggles to estimate surface water condition, which is a critical input for the groundwater model because it lacks a land surface hydrology model. However, integrated SWAT-MODFLOW model overcomes the limitation by replacing the SWAT model module for groundwater with the MODFLOW [43,48].

During simulating the groundwater recharge using SWAT-MODFLOW model, data transfer from individual model is vital for accurate estimation. The data transfer between individual models uses a mapping scheme. The simulated deep percolation of SWAT model (HRUs) mapped as recharge to the Disaggregated HRUs (DHRUs), MODFLOW in each grid cell, while SWAT model will get groundwater discharge in each sub-basin from MODFLOW [41]. SWAT model uses equation (1) to simulate deep percolation (groundwater recharge).

(1)

where ω_rchrg,i is the quantity of water inflow to the aquifers as recharge on day i (mm), δ_gw (days) is lapse period of the aquifer materials, ω_seep is the sum quantity of water exiting below soil profile on day i(mm), and ω_rchrg,i-1 is the quantity water inflow to the aquifers as recharge on day_i-1 (mm).

MODFLOW uses a finite difference approach as stated in equation (2) to describe the groundwater flow at grid level. MODFLOW model discretizes the aquifer into grid cell level, which can be defined in terms of rows, columns, and layers. Simulation such as groundwater head presented in grid cell level for the model. After the SWAT-MODFLOW setup is completed, SWAT groundwater module will be deactivated and MODFLOW groundwater flow equation used in the case of recharge.

(2)

where h is the hydraulic head (m). K_xx, K_yy, and K_zz signify the hydraulic conductivity (m.day^-1) along x, y, and z directions, respectively. The S_s is the specific storage of the aquifer that is a dimensionless quantity, t (day) represents time, and where W (day^-1) signifies the source or sink. If the W is positive, it represents the recharge into the aquifer, while if it is negative it represents abstraction.

2.2.2. Transient Water Table Fluctuation Method (TWTFM)

Water table fluctuation (WTF) method has been utilized in various studies due to its easiness and straightforwardness to apply [67,68,69] and the method described in detail by [12]. The basic concept of this method, the rise of groundwater in the aquifers occurs due to the percolation process in the region. Several monitoring groundwater level data are required to apply this method with specific yield values of the aquifer. For accurate estimation of recharge, specific yield value should be acceptable [12,70,71]. Some studies use regression approaches to estimate specific yield, by developing a relationship between groundwater fluctuation in response to series event precipitation [72]. To calculate recharge with time series, researchers acquire inverse modeling of the water table elevation with temporally changing specific yield [60,71]. However, determining the aquifer parameters, including specific yield from known recharge (especially from numerical model approach), reduces the uncertainty of aquifer parameters value [59,73].

De Zeeuw and Hellinga (1958) suggested analytical equation (3) to determine recharge and discharge during transient state subsurface drainage process [74]. In the present study, we implemented this equation to estimate recharge. The procedure to determine recharge using this method is described in Figure 4. First, prepare selected groundwater (GW) level series data (daily), which should be continuous and reliable. Then determine reaction factor (α) parameter using equation (4) which is solely dependent on observed water level data series. Sum up simulated percolation from the SWAT model and use it for alteration input. Iterate specific yield (μ) value until the sum of SWAT percolation and computed recharge using TWTFM values close to each other.

(3)

where α denotes the reaction factor (d^-1), μ is the specific yield of the aquifer, h_i is the GW level (L) on the i day, h_i-1 is the GW level (L) on the i-1 day, Δt is the time interval (T), and R∆t is average recharge for time interval t-1 to t, and was presumed constant.

(4)

Reaction factor (α) can be calculated based on equation (4), which depends only on the groundwater level data and assumes no recharge in the aquifer [74]. Specific yield (μ) was determined using equation (3) and iteration procedure. Detailed procedure the for TWTFM approach reported by Chung et al. [59].

For this study, seven monitoring wells with reliable and long-period records were selected to apply TWTFM approach as shown in Figure 5. We applied this method even if the selected monitoring wells were not well distributed in the study watershed. Therefore, the number of wells and their distribution were considered a limitation of the TWTFM approach for this study. The whole process and workflow of the method is stated in Figure 4.

2.3. Model Setup and Coupling

2.3.1. SWAT Model

To build the study area, SWAT2012 (ver.636) was employed. SWAT model, is a semi-distributed model, which functions on a daily time series. A digital elevation model (DEM) was applied to delineate into watershed and sub-watershed (subbasins). For this model, the landscape is classified into conceptual units named Hydrological Response Units (HRUs). A piece HRU signifies an area of the landscape that has similar soil, land use, slope, and management characteristics. The SWAT model simulates the quantity and quality of water resources based on the conceptual units (HRUs)[75]. To adjust the model size without missing essential data, setting the proper number of HRUs is a vital procedure [76]. SWAT model has multiple HRU definition options, which allow users to represent the variability of LULC, soil types, and topography within the study watershed. In the present study, multiple HRUs option with threshold values of 5% was optioned to create the HRUs of the watershed.

QSWAT is a graphical user interface (GUI, version 3.22) that was utilized to construct and run the SWAT model of the study watershed. The region was delineated and divided into 21 subbasins (Figure 1) using a 30 m resolution digital elevation model (DEM). Stream network was constructed to duplicate the water network of the study area. Then 1679 HRUs were created with multiple HRUs option. For streamflow calibration, station gauge (SG2) (Figure 1) was used for this study. Weather data were imported and written in the model. From 2006 to 2018, the SWAT model simulation was performed; the initial three years (2006–2008) were designated as a warm-up period, during which the model did not provide results for analysis or interpretation.

2.3.2. MODFLOW Model

The study region conceptual model was developed using the GMS software MODFLOW, version 10.5.8. MODFLOW is employed to solve the groundwater flow of the aquifer in three dimensions finite difference approach [77]. To account for changes in the wetness and confinement of the grid cells over time, this study utilized the MODFLOW-NWT (Newtonian Solver) engine. MODFLOW-NWT is preferred over other versions of MODFLOW due to its ability to simulate hydraulic heads for dry cells without changing them to inactive cells in an unconfined aquifer [42]. In addition, upstream weighting (UPW) package was used, which enables zero magnitude for dry cells flow out and determines the head from the inflow to the dry cells. Overall, utilization of MODFLOW-NWT enables a more accurate analysis of sub-surface water flow for the study watershed [42].

Anyang watershed (with 137 km² coverage) was discretized into grid cells with a lateral resolution of 100 m by 100 m, which resulted in 17,465 active grid cells. Conceptual top elevation data was obtained from SWAT model, DEM raster file. The bottom elevation was determined using data from electrical resistivity tests conducted in the study watershed; elevation ranges 50–200 m below the surface. The stream network created during delineation SWAT process was introduced to MODFLOW as a shape file. Using the river package 825 stream cells were created for the watershed.

During the simulation of steady-state groundwater flow, parameters include hydraulic conductivity, river conductance, initial head, and sources and sinks. Recharge rate was obtained from SWAT percolation output, which provides a single value for MODFLOW. The initial head values were acquired by interpolating from 35 monitoring wells data collected in the study region. To attain initial values for hydraulic conductivity interpolation approach was used. As anticipated, the initial simulation revealed a discrepancy between the observed and calculated head. The conceptual model was calibrated by the Parameters ESTimation (PEST) algorithm, which involved adding more pilot points in areas where data was limited. During calibration, the pilot points were comprised by fixing the minimum and maximum hydraulic conductivity to the same values.

After achieving reasonable results between the observed and computed head during the steady-state simulation, a simulation setup for the transient state was initiated. Before proceeding with the coupling procedure between SWAT and MODFLOW, specific storage and specific yield parameters were established for the transient simulation.

2.3.3. SWAT-MODFLOW Model Integration

The linking code, which integrates both SWAT and MODFLOW models in daily time steps, was introduced by Bailey et al. [41]. It has significant advantages for the models since the coupling approach enables the models to share their computed results without having to rewrite or import their input data. One of the significant steps during the coupling procedure is converting SWAT HRUs to disaggregated Hydrologic Response Units (DHRUs), which conjunction with a MODFLOW grid shape file, to provide spatially explicit information for the model. The procedure made available by Park et al. [78] explains the entire step for creating and transferring the current linking files and integrating the model.

In the present study, QSWATMOD, a QGIS-based graphical interface, was applied to do before processing, configuration, and after processing for coupled model. It is written in Python code, which is suited to create linkage files between SWAT and MODFLOW. The fundamental procedures are described by Bailey et al. [41].

To perform integrating using QSWATMOD, the following inputs are required, SWAT output file, shape files of subbasin, HRU, and river network, and a MODFLOW model folder as native text as a pre-processing procedure. Stream network can be selected from either SWAT or MODFLOW model. For this study, the MODFLOW stream was chosen due to modification of the riverbed conductance. To link the models, DHRUs were created to intersect with MODFLOW grid cells. A zonal polygon was created for aquifer parameters in the MODFLOW model before the coupling process, which can be used during the calibration process of SWAT-MODFLOW. All integrated model files were stored in a single working directory.

2.4. Model Calibration and Performance Evaluation

2.4.1. SWAT Calibration

SWAT model is suitable for simulating streamflow, it can be validated with observed streamflow data. Parameters that influence simulated streamflow should be calibrated until it sounds with observed streamflow values. To do that SWAT-CUP program was applied with the Sequential Uncertainty Fitting program (SUFI-2) algorithm, which considers all roots of uncertainty related to the model, parameters, and input data [79]. Details description of SUFI-2 and other algorithms was explained by [80]. In this work, the simulation of the SWAT model covered 13 years, with the calibration period from 2013 to 2018 and the validation period from 2009 to 2010. In our study, calibration and validation for the observed streamflow were in the daily time series.

From several SWAT hydrological parameters, sensitive parameters were chosen to calibrate the model by providing the acceptable range value. Selected parameters control simulated soil water, surface runoff, evaporation, and groundwater. Detailed descriptions of selected parameters were stated in Table 1.

During calibration, the execution of the model can be assessed using statistical parameters like the Nash-Sutcliffe efficiency (NSE) and the regression coefficient (R²) values [81]. NSE Value higher than 0.5 is acceptable for calibration performance. For NSE nearer to one indicates that the model outputs are significantly nearer to observed data. The regression coefficient (R²) is applied to determine how closely observed and simulated data align. In this study, NSE was selected as the objective function to evaluate the model simulation during the calibration and validation period in SWAT-CUP program. As stated in Figure 1, the study watershed has two stream gauges (SG1 and SG2). Calibration was implemented at SG2 since the other stream gauge has missing data, which makes it difficult to select. For subbasins, where SG2 did not cover regionalization approach was used [82].

2.4.2. SWAT-MODFLOW Model Calibration

The process of acquiring suitable parameters value is essential in modeling process since it helps in the accuracy of the model output. Previous studies have approached several calibration and validation frameworks for SWAT-MODFLOW. The first is to calibrate using automated parameter estimation tools for both SWAT and MODFLOW independently and then recalibrate the integrated model manually [9,22]. The other option is to use automated and semi-automated calibration approaches for coupled model [30,40,83].

SWAT-MODFLOW calibration using PEST framework was first introduced by [84]. In the present work, a parallel implementation of PEST called BEOPEST was used, which minimize the time needed for calibration. The PEST framework requires five types of files to run; control file, PEST-executable file (BEOPEST), model batch file, model input template file, and instruction file (for model output). Detail descriptions about the procedure can be found [30].

To improve streamflow simulation further five parameters from SWAT were selected for recalibration. Aquifer parameters (H_K, S_S, and S_y) were included in the calibration processes of SWAT-MODFLOW. Table 2 describes the selected parameters for calibration with their value ranges.

3. Results and Discussions

3.1. Calibrated Parameters

In this section, the final automated parameters used for integrated SWAT-MODFLOW model were presented. Parameters in Table 3 have lower values compared to the SWAT model calibration parameters value displayed in Table 1. One of the possible reasons for this is to overcome the simulated streamflow lower base flow and peak flow conditions and the statement supported by [85]. Zonal polygons were created preceding to coupling process based on hydraulic head distribution to calibrate aquifer parameters. Figure 6 showed the final automated aquifer parameter value ranges of the study region. Calibrated hydraulic conductivity value range from 0.014 m/day to 0.68 m/day. Yifru et al. [9] reported the hydraulic conductivity range between 0.001 m/day to 26 m/day for three watersheds that include the present work region. They also reported specific storage and specific yield values as 0.0001 m^-1 and 0.18, respectively. For this study, specific storage values range from 0.0017 m^-1 to 0.01m^-1, and specific yield range from 0.03–0.41. Referring to Figure 6, we can draw a relationship between individual aquifer parameters. For example, specific yield and specific storage have inverse relations during calibration of the hydraulic head of the study watershed.

3.2. Model Performance during Calibration and Validation

Figure 7 describes daily streamflow between simulated and measured for SWAT and integrated model. The streamflow was calibrated from 2013–2018 and validated from 2009–2010. Table 4 summarizes model performance statistics of NSE and R² between measured and calculated streamflow for SWAT and SWAT–MODFLOW model. During the integrated model, better performance statistics showed for streamflow than the SWAT model, as stated in Table 4. SWAT-MODFLOW model performs better simulation during low base flow and peak flow conditions of observed streamflow, as shown in Figure 7. According to the evaluation criteria mentioned by [86], both models accomplished well in simulating the streamflow during calibration and validation set time.

The hydraulic head simulation was performed based on the continuity of recorded monitoring wells. Several of the observed wells in our study do not have continuous daily time series data, which is considered a limitation of our study in the region. The selected monitoring wells have records from May 2017 to September 2018. Figure 8 describes the hydraulic head of the study watershed. The calculated depth to water values was compared with the observed, as displayed in Figure 8a. The negative sign indicates that the water level is level below the ground level. As shown in Figure 8a, good trends were visualized between observed and calculated groundwater levels, with the difference between the heads less than 1 m. Figure 8b displays the spatial distribution of simulated hydraulic head of our study watershed. The hydraulic head ranges from 9.7 m to 453.5 m (Figure 8b). The middle part of the study region shows low hydraulic heads while most of the north and east area have higher values.

3.3. Groundwater Recharge

Figure 9 displays the distribution of average groundwater recharge during the whole simulation period (2009–2018). The figure demonstrates a considerable spatial and temporal distribution of groundwater recharge in the watershed. The monthly average recharge distribution describes these patterns in the area that ranges from 0–22.4 m³/day (Figure 9a). In the present work, the groundwater recharge amount is complicatedly distributed and problematic to connect with the physical topographies of the region, even though the northeast has a maximum recharge. The study area precipitation was intense during the wet season of the summer period (June–August), and the groundwater recharge shows this tendency (Figure 9c). During the wet season, the average recharge ranges from 0–58.6 m³/day. In the region, the dry season (January–March) shows the lowest value among a yearlong recharge. During this season, most areas get the lowest recharge and the highest value is 4.5 m³/day (Figure 9b). Topography with high elevation shows relatively high recharge values.

Yifru et al. [9] stated that the recharge covers 18% of the annual average precipitation in the Han River watershed which includes our study region. Further, the peak and significant amount of recharge happen in the course of the summer period, and July month has the largest recharge value. In the present study, the groundwater recharge covers 34 % of the average annual precipitation of the region. Even though, it is not indicated numerically similar recharge pattern was shown graphically in this study region by Yifru et al. [9]. Whereas, the minimum recharge occurs during a month of the winter season with the highest value of 8.85 mm.

Table 5 shows the simulated recharge amount using the TWTFM approach. In addition, it listed the estimated parameters (specific yield and reaction factor) in the selected subbasin of the study region.

Groundwater level records and initial recharge input data are essential to apply the transient water table fluctuation method. From available monitoring wells in our study, seven were selected to apply TWTFM approach. For initial recharge input, SWAT percolation values at subbasin level where observed wells exist. For this method, recharge was estimated from May 2017 to September 2018. Table 5 describes the calculated aquifer parameters and recharge for selected monitoring wells in the study watershed. Specific yield ranges from 0.0078–0.022, while the reaction factor ranges from 0.04 d^-1 to 0.2 d^-1. The specific yield values were considered as initial input values in calibration of coupled model. Estimated recharge ranges between 553.6 mm to 1020 mm for this method, which is shown in Table 5. TWTFM provides almost similar results of recharge with SWAT model. For SWAT-MODFLOW, the average annual recharge was 430 mm for the whole watershed. For this method to draw a representative recharge value for the study watershed, well-distributed monitoring wells are required, which consider a limitation (Figure 5).

3.4. Water Balance Components of Study Watershed

Figure 10 illustrates the major water balance segments from integrated model for the period 2009–2018 (until September). Percolation is the major contributor to the water cycle of the study region. Based on Figure 10, the maximum monthly percolation occurs in June 2011 with 11.7 mm value where the precipitation was also maximum. Further, surface runoff and lateral flow show similar peak values during this time. The yearly average evapotranspiration is around 252 mm, almost 20% of the annual precipitation of the region. The direct flow components (surface and lateral) were significant, which is around 45% of the annual precipitation in the watershed. The groundwater volume shows decreasing pattern to adjust the initial hydraulic head of the region.

4. Conclusions

Numerical model and empirical approach were utilized to evaluate the groundwater recharge of Anyang watershed. A numerical model, SWAT-MODFLOW, was used to assess the spatiotemporal groundwater recharge distribution of the study region. Observed streamflow and groundwater level were calibrated for the coupled model. Streamflow was calibrated from the year 2013–2018 and validated from 2009–2010. The performance of the model was evaluated using statistical parameter function of R² and NSE. The integrated model showed better performance in simulating the streamflow during calibration and validation period. Re-calibration of selected SWAT parameters was a proven reason for better streamflow simulation for integrated model, other than aquifer parameters calibration. Monitoring wells were evaluated with calculated water level below the ground level, which showed a good trend with the measured value.

SWAT-MODFLOW model quantified the average monthly groundwater recharge in the region and showed the spatiotemporal distribution variations yearlong. Seasonal recharge in the region was well spotted in the model and a significant difference was shown between the dry and wet seasons. The groundwater recharge contributes 34% of the annual average precipitation even though the region has large urban area coverage. The contribution of the direct flow component (surface and lateral) was significant, which relates to the geological characteristics of the study region.

The other method applied to evaluate recharge was the transient water fluctuation method (TWTFM). Seven monitoring wells with long-term record daily series were selected for this method. It provides recharge that ranges from 553.6–1020 mm and aquifer parameters that were used as initial input during calibration of coupled.

In general, the present work showed the utilization and performance of the SWAT-MODFLOW model and TWTFM approach to analyze groundwater recharge. SWAT-MODFLOW shows the spatial and temporal distribution of the groundwater recharge of the study watershed. Hence, it provides a better representation of the surface water and groundwater of the region. TWTFM shows a glimpse as an alternative option for recharge estimation, but it requires well-distributed groundwater level data to represent recharge condition of the region.