Assessment of Spatiotemporal Groundwater Recharge Distribution Using SWAT-MODFLOW Model and Transient Water Table Fluctuation Method

: Recharge is a crucial section of water balance for both surface and subsurface models in water resource assessment. However, quantifying its spatiotemporal distribution at a regional scale poses a signiﬁcant challenge. Empirical and numerical modeling are the most commonly used methods at the watershed scales. However, integrated models inherently contain a vast number of unknowns and uncertainties, which can limit their accuracy and reliability. In this work, we have proposed integrated SWAT-MODFLOW and Transient Water Table Fluctuation Method (TWTFM) to evaluate the spatiotemporal distribution of groundwater recharge in Anyang watershed, South Korea. Since TWTFM also uses SWAT model percolation output data, calibration was performed for individual models and a coupled model. The coupled model was calibrated using daily streamﬂow and hydraulic head. The SWAT-MODFLOW model performed well during the simulation of streamﬂow compared to the SWAT model. The study output showed that the study watershed had signiﬁcant groundwater recharge variations during the simulated period. A signiﬁcant amount of recharge happens in the wet season. It contributes a signiﬁcant amount of the average annual precipitation of the region. The direct ﬂow components (surface and lateral) showed signiﬁcant contributions when the water balance components were evaluated in the region. TWTFM showed a glimpse to estimate recharge, which requires representative monitoring wells in the study region. Comprehensively, the SWAT-MODFLOW model estimated groundwater recharge with reasonable accuracy in the region.


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]. 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's emphasis is 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]. However, each method has its limitation while assessing recharge [11]. The numerical groundwater contributes a high role, as stated by [18,51]. In recent hydrological studies, parametrization tools have been applied for the SWAT-MODFLOW model, which reduces the uncertainty of parameters applied for the calibration of the model [30,52]. 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; for example, in the case of streamflow evaluation, it showed high precision [9,30,40,44]. However, numerical models like SWAT-MODFLOW always have uncertainty [53]. 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 [42]. For example, Guevara-Ochoa et al. [43] mentioned that the 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, the 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,54]. The method requires observed groundwater level data series and aquifer parameters, which can be determined using various methods [11,12,55]. To reduce subjectivity and increase efficiency, several modifications have been made to this method [56][57][58]. One of the approaches was to link the WTF method with a hydrological model, such as the SWAT model, which helps to overcome local representation of the method. For instance, Chung et al. [57] estimated recharge by linking the 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,54,59,60]. This study aims to apply the 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 a 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.

Description of Study Region
The Anyang watershed is located in northwestern South Korea and is displayed in Figure 1. The study watershed covers an approximately 137 km 2 area and drains into the Han River, one of the longest rivers in the Korean Peninsula. The altitude ranges from 11-591 m above mean sea level at 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 • C, respectively. For a while, the streams of the study region have been gauged; this started before 2002. During the summer period, the streamflow achieves its maximum value.
Anyang watershed land use/land cover consists of mostly forest and urban areas, with forest enclosing approximately 50% of the whole area and urban industrial areas covering 43%. The remaining 7% is composed 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. Anyang watershed land use/land cover consists of mostly forest and urban areas, with forest enclosing approximately 50% of the whole area and urban industrial areas covering 43%. The remaining 7% is composed 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 hydrological soil groups of the study are displayed in Figure 3. In the study area, most soil types have 4 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. It is 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 [61] and Metrological Agency Weather Data Service [62], respectively. As shown in Figure 1, the weather station is outside of the watershed boundary location and it may have some uncertainties if simulated weather data used for hydrological model. Therefore, measured weather data were used to model the study watershed. The hydrological soil groups of the study are displayed in Figure 3. In the study area, most soil types have 4 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. It is 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 [61] and Metrological Agency Weather Data Service [62], respectively. As shown in Figure 1, the weather station is outside of the watershed boundary location and it may have some uncertainties if simulated weather data used for hydrological model. Therefore, measured weather data were used to model the study watershed.   The SWAT model's ability to address groundwater flow is restricted by its aggregated nature, while MODFLOW struggles to estimate the surface water condition, which is a critical input for the groundwater model because it lacks a land surface hydrology model. However, the integrated SWAT-MODFLOW model overcomes the limitations by replacing the SWAT model module for groundwater with the MODFLOW [42,46].

Evaluation of Recharge Using Integrated Model and Transient Water
During simulating the groundwater recharge using the SWAT-MODFLOW model, data transfer from individual models is vital for accurate estimation. The data transfer between individual models uses a mapping scheme. The simulated deep percolation of the SWAT model (HRUs) mapped as recharge to the Disaggregated HRUs (DHRUs),  The SWAT model's ability to address groundwater flow is restricted by its aggregated nature, while MODFLOW struggles to estimate the surface water condition, which is a critical input for the groundwater model because it lacks a land surface hydrology model. However, the integrated SWAT-MODFLOW model overcomes the limitations by replacing the SWAT model module for groundwater with the MODFLOW [42,46].

Evaluation of Recharge Using Integrated Model and Transient Water
During simulating the groundwater recharge using the SWAT-MODFLOW model, data transfer from individual models is vital for accurate estimation. The data transfer between individual models uses a mapping scheme. The simulated deep percolation of the SWAT model (HRUs) mapped as recharge to the Disaggregated HRUs (DHRUs), MODFLOW in each grid cell, while the SWAT model will obtain groundwater discharge in each subbasin from MODFLOW [40]. The SWAT model uses Equation (1) to simulate deep percolation (groundwater recharge).
where ω rchrg,i is the quantity of water inflow to the aquifers as recharge on day i (mm), δ gw (days) is the 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 of 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. The MODFLOW model discretizes the aquifer into grid cell level, which can be defined in terms of rows, columns, and layers. Simulations such as groundwater head are presented at grid cell level for the model. After the SWAT-MODFLOW setup is completed, the SWAT groundwater module will be deactivated and the MODFLOW groundwater flow equation used in the case of recharge.
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 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.

Transient Water Table Fluctuation Method (TWTFM)
The Water Table Fluctuation (WTF) method has been utilized in various studies due to its ease and straightforwardness to apply [64][65][66], and the method is described in detail by [12]. The basic concept of this method is that 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,67,68]. Some studies use regression approaches to estimate specific yield, by developing a relationship between groundwater fluctuation in response to series event precipitation [69]. To calculate recharge with time series, researchers acquire inverse modeling of the water table elevation with temporally changing specific yield [58,68]. However, determining the aquifer parameters, including specific yield from known recharge (especially from a numerical model approach), reduces the uncertainty of aquifer parameters values [57,70].
De Zeeuw and Hellinga (1958) suggested analytical Equation (3) to determine recharge and discharge during the transient state subsurface drainage process [71]. 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. 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.   The reaction factor (α) can be calculated based on Equation (4), which depends only on the groundwater level data and assumes no recharge in the aquifer [71]. Specific yield (µ) was determined using Equation (3) and iteration procedure. The detailed procedure for the TWTFM approach is reported by Chung et al. [57].
For this study, seven monitoring wells with reliable and long-period records were selected to apply the TWTFM approach to, 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 are stated in Figure 4.

SWAT Model
To build the study area, SWAT2012 (ver.636) was employed. The 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). One 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) [72]. To adjust the model size without missing essential data, setting the proper number of HRUs is a vital procedure [73]. The 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 options with threshold values of 5% were 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). The stream network was constructed to duplicate the water network of the study area. Then 1679 HRUs were created with multiple HRUs options. For streamflow calibration, a station gauge (SG2) (Figure 1) was used for this study. Weather data were imported and

SWAT Model
To build the study area, SWAT2012 (ver.636) was employed. The 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). One 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) [72]. To adjust the model size without missing essential data, setting the proper number of HRUs is a vital procedure [73]. The 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 options with threshold values of 5% were 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). The stream network was constructed to duplicate the water network of the study area. Then 1679 HRUs were created with multiple HRUs options. For streamflow calibration, a 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)(2007)(2008) were designated as a warm-up period, during which the model did not provide results for analysis or interpretation.

MODFLOW Model
The study region conceptual model was developed using the GMS software MOD-FLOW, version 10.5.8. MODFLOW is employed to solve the groundwater flow of the aquifer using a three-dimensional finite difference approach [74]. 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 [41]. In addition, the 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 subsurface water flow for the study watershed [41].
The Anyang watershed (with 137 km 2 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 were obtained from the SWAT model and 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 the 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. The 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, the interpolation approach was used. As anticipated, the initial simulation revealed a discrepancy between the observed and calculated head. The conceptual model was calibrated using the Parameters ESTimation (PEST) algorithm, which involved adding more pilot points in areas where data were 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.

SWAT-MODFLOW Model Integration
The linking code, which integrates both SWAT and MODFLOW models in daily time steps, was introduced by Bailey et al. [40]. 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 in conjunction with a MODFLOW grid shape file, provide spatially explicit information for the model. The procedure made available by Park et al. [75] explains the entire process for creating and transferring the current linking files and integrating the model.
In the present study, QSWATMOD, a QGIS-based graphical interface, was applied before processing and configuration, and after processing for the 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. [40].
To perform integrating using QSWATMOD, the following inputs are required: SWAT output file, shape files of subbasin, HRU, river network, and a MODFLOW model folder as native text as a pre-processing procedure. The stream network can be selected from either the 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.

SWAT Calibration
The SWAT model is suitable for simulating streamflow; it can be validated with observed streamflow data. Parameters that influence simulated streamflow should be calibrated until they match with observed streamflow values. To do that, the 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 [76]. Detailed descriptions of SUFI-2 and other algorithms were explained by [77]. 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 such as the Nash-Sutcliffe efficiency (NSE) and the regression coefficient (R 2 ) values [78]. An NSE value higher than 0.5 is acceptable for calibration performance. An NSE nearer to one indicates that the model outputs are significantly nearer to observed data. The regression coefficient (R 2 ) is applied to determine how closely the 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 that SG2 did not cover, the regionalization approach was used [79].

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 the coupled model [30,80].
SWAT-MODFLOW calibration using PEST framework was first introduced by [81]. In the present work, a parallel implementation of PEST called BEOPEST was used, which minimizes 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). Detailed descriptions about the procedure can be found in [30].
To improve streamflow simulation, a 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.

Calibrated Parameters
In this section, the final automated parameters used for integrated SWAT-MODFLOW model are 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 [82]. Zonal polygons were created preceding to coupling process based on hydraulic head distribution to calibrate aquifer parameters. Figure 6 shows the final automated aquifer parameter value ranges of the study region. Calibrated hydraulic conductivity values 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 3 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.01 m −1 , and specific yield range from 0.03 to 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.     Table 4 summarizes model performance statistics of NSE and R 2 between measured and calculated streamflow for the SWAT and SWAT-MOD-FLOW model. During the integrated model, better performance statistics are shown for streamflow than the SWAT model, as stated in Table 4. The 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 [78], both models accomplished well in simulating the streamflow during calibration and validation set time.   Table 4 summarizes model performance statistics of NSE and R 2 between measured and calculated streamflow for the SWAT and SWAT-MODFLOW model. During the integrated model, better performance statistics are shown for streamflow than the SWAT model, as stated in Table 4. The 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 [78], both models accomplished well in simulating the streamflow during calibration and validation set time.

Model Performance during Calibration and Validation
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 below the ground level. As shown in Figure 8a, good trends were visualized between the 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. 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 below the ground level. As shown in Figure 8a, good trends were visualized between the 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.     Figure 9 displays the distribution of average groundwater recharge during the whole simulation period (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(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 to 22.4 m 3 /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 to 58.6 m 3 /day. In the region, the dry season (January-March) shows the lowest values among a yearlong recharge. During this season, most areas obtain the lowest recharge and the highest value is 4.5 m 3 /day (Figure 9b). Topography with high elevation shows relatively high recharge values.  Table 5 shows the simulated recharge amount using the TWTFM approach. In addition, it lists 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 the 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 to 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 and 1020 mm for this method, which 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 the month of July 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, a numerically similar recharge pattern was shown graphically in this study region by Yifru et al. [9]. In contrast, 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 lists the estimated parameters (specific yield and reaction factor) in the selected subbasin of the study region.

Groundwater Recharge
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 the 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 to 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 and 1020 mm for this method, which is shown in Table 5. TWTFM provides almost similar results of recharge with the 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).  Figure 10 illustrates the major water balance segments from the 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 a 11.7 mm value where the precipitation was also the 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 is influenced by the aquifer parameter properties and values. SWAT-MODFLOW model calibrates the groundwater volume and hydraulic head based on the optimized parameters. It shows a decreasing pattern to adjust the initial hydraulic head of the region. is shown in Table 5. TWTFM provides almost similar results of recharge with the 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). Figure 10 illustrates the major water balance segments from the 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 a 11.7 mm value where the precipitation was also the 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 is influenced by the aquifer parameter properties and values. SWAT-MODFLOW model calibrates the groundwater volume and hydraulic head based on the optimized parameters. It shows a decreasing pattern to adjust the initial hydraulic head of the region.

Conclusions
A numerical model and empirical approach were utilized to evaluate the groundwater recharge of the Anyang watershed. A numerical model, SWAT-MODFLOW, was used to assess the spatiotemporal groundwater recharge distribution of the study region. The observed streamflow and groundwater level were calibrated for the coupled model.

Conclusions
A numerical model and empirical approach were utilized to evaluate the groundwater recharge of the Anyang watershed. A numerical model, SWAT-MODFLOW, was used to assess the spatiotemporal groundwater recharge distribution of the study region. The observed streamflow and groundwater level were calibrated for the coupled model. Streamflow was calibrated from the years 2013-2018 and validated from 2009 to 2010. The performance of the model was evaluated using the statistical parameter function of R 2 and NSE. The integrated model showed better performance in simulating the streamflow during the calibration and validation period. Recalibration of selected SWAT parameters was a proven reason for better streamflow simulation for the integrated model, other than aquifer parameters calibration. Monitoring wells were evaluated from May 2017 to September 2018 with calculated water level below the ground level, which showed a good trend with the measured value.
The 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 to 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 the 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.

Conflicts of Interest:
The authors declare no conflict of interest.