Mapping Interflow Potential and the Validation of Index-Overlay Weightings by Using Coupled Surface Water and Groundwater Flow Model

: Interﬂow is an important water source contributing to river ﬂow. It directly inﬂuences the near-surface water cycles for water resource management. This study focuses on assessing the interﬂow potential and quantifying the interﬂow in the downstream area along the Kaoping River in southern Taiwan. The interﬂow potential is ﬁrst determined based on the modiﬁed index-overlay model, which employs the analytical hierarchy process (AHP) to calculate the ratings and weightings of the selected factors. The groundwater and surface water ﬂow (GSFLOW) numerical model is then used to link the index-overlay model to quantify the interﬂow potential for practical applications. This study uses the Monte Carlo simulations to assess the inﬂuence of rainfall-induced variations on the interﬂow uncertainty in the study area. Results show that the high potential interﬂow zones are located in the high to middle elevation regions along the Kaoping River. Numerical simulations of the GSFLOW model show an interﬂow variation pattern that is similar to the interﬂow potential results obtained from the index-overlay model. The average interﬂow rates are approximately 3.5 × 10 4 (m 3 /d) in the high elevation zones and 2.0 × 10 4 (m 3 /d) near the coastal zones. The rainfall uncertainty strongly inﬂuences interﬂow rates in the wet seasons, especially the peaks of the storms or heavy rainfall events. Interﬂow rates are relatively stable in the dry seasons, indicating that interﬂow is a reliable water resource in the study area.


Introduction
In hydrology, interflow is the lateral movement of water in the unsaturated zone. It can return to the surface or enter a stream before becoming groundwater. Interflow is also the water flow zone in a porous media right below river beds, and the interflow is the interface for the interactions between surface water and groundwater. Interflow occurs when water infiltrates into the subsurface. When hydraulic conductivity decreases with depth, lateral flow proceeds downslope. It may exfiltrate as return flows, becoming overland flow. Such interface is important to the bioenvironmental system near the rivers. Exchange of water across the interface can result from local and upstream to downstream water movement in and out of the river beds and banks. Water exchange across the surface water-groundwater interface has been explored in some detail in the past few decades, with most studies focused on rivers, and it is increasingly studied concerning the effects on the chemical composition of surface and subsurface water and the distribution of biota [1][2][3][4]. In many arid and semiarid areas, interflow also plays an essential role in water resource management in dry seasons when the water levels in rivers are relatively low [5][6][7]. With varied climate conditions, high river discharge variations in dry and wet seasons are generally observed. Therefore, the assessment of available water resources from interflow has become a key process to support land-use planning and regional development. This available water resource from interflow can be considered as the interflow potential for specific reaches of rivers.
The concept of index-overlay models to create different types of potential or vulnerability maps has been widely used for many implementations. Such potential or vulnerability maps are critical references to decision-making for land-use planning, groundwater management, monitoring, and remediation [8][9][10][11]. Previous investigations have shown that the advantages of the index-overlay models are their ease of use, minimal data demand, and precise definition of parameters for site-specific problems. With the recent advances in Geographical Information Systems (GISs), index-overlay models are especially helpful for creating various types of potential maps in the fields of groundwater resource potential [12][13][14]. A GIS associated with useful tools enables the overlay and integration of multiple maps, representing the selected parameters to account for the groundwater potential [15][16][17]. The development of interflow potential maps can follow the same concept as groundwater potential. In an interflow potential map, the value of the potential represents the relative availability of interflow in river segments.
The interflow potential maps are fundamental for assessing groundwater and surface water interaction in a local water cycle. For most practical problems, developing interflow potential maps is efficient because the intrinsic parameters such as the river slope, river stage, river bed material, and the groundwater level near the river are usually the available hydrogeological data from the specific sites [18,19]. Users determine the weightings and ratings of different selected parameters. The summation of the weightings and ratings then creates the relative potential for a site-specific problem [20][21][22]. It is essential to validate the developed potential maps for practical cases. However, the validation of the interflow potential might require interflow observations, and the direct observations are usually limited. In many studies, the observations of interflow are based on the calculation of water balance in different reaches of the rivers [23][24][25]. The river discharge values across different river sections must account for the discharge differences of interflow, specifically between two river sections. The net interflow rates in different river sections are then the basis to validate the interflow potential map. River bed infiltration tests could show direct measurements of water flow in and out of the river bed [26][27][28][29]. However, the relatively small scale of the infiltration tests might not represent well the interflow dynamics at the target site. The question of how to upscale and integrate different types of observations is challenging because it is necessary to validate the interflow potential obtained from index-overlay models. A physical-based model that could reproduce the surface water and groundwater dynamics might be useful to solve this challenge. Previous studies such as Grodzka-Łukaszewska et al. (2017) [30] or Diaz et al. (2020) [31] used different models to account for the interactions between surface water and groundwater. However, the interflow was not included in the conceptual model in their studies. The study aims to develop a workflow that uses the physical-based numerical model to validate the indexoverlay model and quantify the interflow in the downstream area along the Kaoping River. The study uses the analytic hierarchy process (AHP) to modify the weightings and ratings in the index-overlay model [10,11]. The quantification of the interflow potential relies on the integrated GSFLOW model to account for groundwater and surface water interactions. The study developed a basin-scale of the modular groundwater flow (MODFLOW) model for the entire Pingtung Plain groundwater basin (PPGB). The results of the basin-scale model then provide head conditions for a refined sub-model in the downstream area along the Kaoping River. The refined model then uses the GSFLOW model to quantify the interaction between surface water and groundwater in the local area. With the available river and groundwater observations, the basin-and local-scale transient models are calibrated and validated based on the three-year observations from January 2015 to December 2017.
The accuracy evaluations for the numerical models use different statistical parameters such as the coefficient of determination (R 2 ), coefficient of correlation (RE), Nash-Sutcliffe efficiency coefficient (NSE), normalized root-mean-square errors (NRMSE), and percent bias (Pbias) [32][33][34][35]. The interflow from the GSFLOW model is then used to evaluate the accuracy of the interflow potential map obtained from the index-overlay model. The study also uses Monte Carlo simulations (MCS) to quantify the rainfall-induced uncertainty for the refined modeling area along the Kaoping River.

Study Area
The PPGB covers approximately 1300 km 2, and land surface elevation varies from the north to the south, located at Pingtung County in Taiwan ( Figure 1). The annual precipitation in the PPGB varies significantly, from 1431 to 3597 mm, and the average precipitation is 2331 mm. Intense precipitation related to typhoons occurs three or four times per year on average [36]. Directions of major flow paths on the PPGB are similar to the terrain slope that is from the northeast to the southwest. In the PPGB, the Kaoping River is the largest river system, located in the south-west of the country, including part of the foothills in the west and flat plains in the south. Three tributaries contribute flow to the Kaoping River in its upper reaches, including the Laonung, Ailiao, and Qishan Creeks, and the discharge basins cover Kaohsiung City and Pingtung County. The upper and middle reaches of the Kaoping River originate from the high slopes of the mountainous area. The riverbed material is deposited in the middle and lower reaches, where the water flows through the mountainous zones into the delta plain. This sedimentation forms the main soil texture of the PPGB. These sediments are the primary materials distributed on the surface of the PPGB. In the PPGB, hydrogeological surveys have been conducted since 1995, and field hydraulic test data have been produced for the aquifer system [37]. The monitoring network was developed to monitor groundwater levels in different aquifer layers. To date, a total of 51 boreholes and 126 monitoring wells are available in the groundwater basin. The hydrogeological surveys showed that the sediment deposits in the downstream fan had formed four aquifer layers that are separated by low permeable marine sequences [37]. Field geological logs showed that there are no marine sequences observed in aquifers in the upstream fan area.
The study focuses on the downstream area of the Kaoping River, where considerable interflow potential has been recognized to support regional water resources (see the marked area in Figure 1) [38]. Figure 1 also shows the variations of the land surface in the PPGB. In the PPGB, there are foothills to the west and a relatively flat area near the Kaoping River. The high precipitation variations in the wet and dry seasons in the PPGB have made limited usage of surface water resources. Such difficulty leads to a challenging task in integrating available water resources for efficient and sustainable management. The downstream area along the Kaoping River (marked with a black polygon in Figure 1) has been recognized as the high interflow potential zone for regional water resources [39]. In this area, there are many interflow collection stations beside the river, and water collection pipes have been installed under the river to collect water for public use. Such public instruments support the alternative strategy to use the water resource from the river (i.e., the surface water), the groundwater, or both. The water collected from the porous media under the river is beneficial because of the relatively stable flow rate in dry seasons and high water quality in wet seasons.

Data Sources and Assumptions
The study aims to develop a workflow that uses a physically-based numerical model to validate an index-overlay model and quantify interflow in the downstream area along the Kaoping River. The basin-scale and local-scale numerical models comprise three years of simulations, particularly a calibration from January 2015 to December 2015 and a validation from January 2016 to December 2017. The steady-state model was calibrated based on the monthly averaged data in December 2014, which is a month earlier than the starting time of the transient simulations. In the refined sub-model along the downstream of the Kaoping River, there are three streamflow stations, namely Li-lin, Gaoping, and Wanda stations (see Figure 1). Because of the missing data at the Gaoping station, the calibration and validation processes mainly used the data from the Li-lin and Wanda stations. Based on the simulated GSFLOW results for the surface water and groundwater interactions, the volumetric interflow rates, such as one of the relative flows, are determined. In the GSFLOW model, the simulations rely on the predefined hydrologic response units (HRUs). In the study, the basin-scale MODFLOW model was developed with conditions similar to previous work in Tran et al. (2020) [40]. In this study, a sub-model was developed based on the specified heads obtained from the basin-scale MODFLOW model. In the sub-model, there are refined HRUs in the downstream area of the Kaoping River. Additional river discharge observations at the Wanda station were included in the sub-model. Based on the observations of groundwater levels and river discharge, calibration and validation of the GSFLOW model are required to represent the system response in the sub-model. In the study, the refined HRUs enabled the quantification of detailed interflow variation. Results could be used to assess the interflow potential maps created by the concept of index-overlay.
The study focuses on the high interflow potential zone. Interflow potential zones were assessed by the AHP technique that is used to modify the weightings and ratings of the selected factors in the index-overlay model. In the study, the selected factors in creating interflow potential maps are the land surface slope (S), transmissivity (T), the difference between river and groundwater levels (DH), depth to groundwater level (DG), and drainage density of the river system (DD). Specifically, the determination of the factors is based on the three-year averaged observations of the time-series observations in the study area. The slope map of the PPGB was derived from the digital elevation model (DEM) using the spatial analysis package in ArcGIS software. The DEM data were collected from the numerical terrain model of Taiwan, with a grid spacing of 20 m. Transmissivity describes the general ability of an aquifer to transmit water. The concept of transmissivity is similar to hydraulic conductivity, but the main difference is that transmissivity is a measurement that applies across the thickness of an aquifer. In the study, the transmissivity data were obtained from a report by the Taiwan Central Geological Survey [37].
The parameter DH represents the differences between river stages and groundwater levels along the river. DH indicates that the available water volume near the river channel controls the local interactions between groundwater and surface water. In addition, the water resource from interflow relies on the available water to be used for the water supply. Therefore, a small discrepancy between groundwater levels and river stages might yield high scores in calculating the interflow potential. The factor of DG was calculated based on the groundwater levels obtained from the monitoring network in the PPGB. In the study, data related to river and groundwater levels were provided by the Water Resources Agency in Taiwan and other studies [41,42]. A lower DG value signifies that more water volume is available near the land surface. Drainage density (DD) defines the total drainage length in a unit area, and the value indicates the closeness of the spacing of channels. In the study area, there is no artificial lining for the streams and channels. Therefore, we have assumed that the stream beds connect well with the aquifer system. Higher DD values can lead to higher interactions between groundwater and surface water. The factor represents the degree of fluvial dissection and is influenced by the resistance to erosion, infiltration capacity, vegetation cover, surface roughness and runoff indices, and climatic condition [43]. The DD is an important indicator that provides a numerical measurement of landscape dissection and surface runoff potential. The DD has been computed by using the line density analysis tool in the software ArcGIS.

Index-Overlay Model
The present study considers multiple parameters for delineating the interflow potential map by using the Geographic Information System (GIS). The weightings for the potential map are determined based on the analytic hierarchy process (AHP) technique. In the study, the selected thematic parameters that influence interflow properties include the land surface slope (S), transmissivity (T), the difference between river and groundwater levels (DH), depth to groundwater level (DG), and drainage density of the river system (DD). These index parameters are represented with spatial distribution maps based on the weighted analysis of the AHP technique [44]. Each thematic factor in pairwise comparison is assessed for the relative rate by its influences on the interflow. All selected thematic factors were represented as raster data with a fixed pixel size of 20 × 20 m for all the map layers. The concept of the AHP is to break down the complex multicriteria problem into a hierarchy. It is based on the pairwise comparison of different criteria [45]. The pairwise comparison was then developed to identify the relative influence of a specific parameter on the interflow potential. The rescaled weightings were then determined based on the results from the numerical model.
In the study, the major and minor effects were assigned the scores of 1 and 0.5 following the concept proposed by previous studies (e.g., [20,21]). Note that judgment of the major and minor effects is based on the direct or indirect influence of the factors on the interflow potential. The influence between selected factors was also determined based on a similar concept. Figure 2 shows the influences of the selected thematic factors for generating the interflow potential map [20,21]. Table 1 summarizes the relative ratings for the thematic factors. Among the proposed relative rating of the five selected thematic factors, slope and transmissivity factors have shown relatively high ratings to calculate the interflow potential. The difference between river and groundwater levels has a mild rating compared to the slope and the transmissivity. The depth to groundwater and the drainage density show less influence on the interflow potential ratings. The ratings associated with the weightings for the factors were then integrated to create the potential map.   Figure 3 shows the distribution of the selected factors for generating the interflow potential map. Figure 3a is the distribution of the slope for the downstream area of the Kaoping River. The slope is estimated as an important parameter that influences the interflow potential in the study area. In Figure 3a, the slope map clearly shows the relatively flat areas with the slope varies from 0 • to 2 • along the Kaoping River. In particular, the slope in the study area varies from 0 • to 62 • . The local high slope is mainly distributed in the hill zone on the western side along the Kaoping River. In the low slope areas, we also obtained low DH and DG values (see Figure 3c,d), which might show relatively large water volumes in the aquifer system. Such large water volumes near the river system could indicate the relatively high interaction between groundwater and surface water.  Figure 3b shows the distribution of the transmissivity interpolated based on the field hydraulic tests conducted by Taiwan Central Geological Survey (Taiwan CGS). The study employed the ordinary Kriging method to create the interpolated transmissivity map. In general, the transmissivity varies from upstream to downstream along the Kaoping River. Figure 3c,d shows the DH and DG maps based on a similar interpolation algorithm by considering three-year averaged groundwater levels and river stages. In the study, the defined cell size in the DEM has been used to calculate the drainage density. Figure 3e shows the drainage density for estimating the interflow potential map. We found that the slope of the river basin also influences the development of the drainage pattern. In the foothill zone, there are local high drainage density areas. However, the high and continuous drainage density areas are mainly near the main channel of the Kaoping River (see Figure 3e).
The data shown in Figure 3 are fundamental to create the interflow potential map. We further classified the data map into five different levels in the study, including very high, high, moderate, low, and very low. The determination of the thresholds depends on the specific characteristic of each selected factor. Table 2 lists the selected thematic factors and the determined weightings for the interflow potential. The weightings, which vary from 5 to 1, represent the five influenced levels. In the study, the assigned weights (AWs) refer to the fundamental scale of nine levels proposed by Satty (1980) [46]. Based on the normalized weightings, the interflow potential could be obtained by summing the products of the respective factors and weights. The equation for the calculation of the interflow potential (IP) has the following formula: where S, T, DH, DH, DG, and DD represent the five hydrogeological parameters included in assessing interflow potential. Subscriptions r and w denote the corresponding ratings and weightings for the parameters. Note that the minimum and maximum values of the different factors are usually varying in different ranges. Therefore, the study used the normalized weights (Table 2) for the presentation. The advantage of normalization is to break down the relations of anomalies among different parameter maps.

GSFLOW Numerical Model
The GSFLOW model proposed by Markstrom et al. (2008) [47] is an integrated hydrological model that couples the Precipitation Runoff Modeling System (PRMS) with the Modular Groundwater Flow Model (MODFLOW). Specifically, the modeling framework in GSFLOW was divided into different conceptual regions, including the region between the plant canopy and the soil zone, the region of streams and lakes, and the region beneath the soil zone (i.e., under the water level). GSFLOW relies on PRMS to simulate the responses of the hydrologic processes in the regions between the plant canopy and the soil zone. Other regions require the solutions obtained from the MODFLOW model and the supporting packages. The rivers and lakes regions arecrucial for representing the interactions between surface water and groundwater. In the MODFLOW, the groundwater flow with the constant density can be described by the following formula [48]: where K x , K y , K z [LT −1 ] are the hydraulic conductivity in the x-, y-, and z-directions, h = h(x, y, z, t) [L] is the hydraulic head. W = W(x, y, z, t)[T −1 ] is a volumetric flux per unit volume, representing the sources and sinks of an aquifer system. The parameter S is the storage coefficient of the aquifer, and t [T] is time. The solution to the above equation is based on a three-dimensional finite-difference scheme. Flow associated with external stresses (i.e., the W in Equation (2)), such as wells, areal recharge, evapotranspiration, drains, and rivers, can be simulated with steady and transient conditions.
In the present study, assessment of the interflow is focused on the downstream area of the Kaoping River. The challenging task is to quantify the interflow in a local area and consider the regional groundwater dynamics that might influence the local water budget. The concept of sub-model is therefore employed to account for the interflow in a local area. The study developed a basin-scale groundwater model (i.e., the MODFLOW model) and a sub-model that focused on the downstream area of the Kaoping River (see Figure 1 for details). The refined HRUs in the sub-model have been conducted by the GSFLOW model (see Figure 4). Figure 4a,b shows the conceptual models for the basin-scale model and the sub-model. The boundary conditions of the groundwater flow in the sub-model were assigned based on the head results obtained from the basin-scale MODFLOW model. Note that the basin-scale groundwater model is a transient state model with a fixed time step of one month. The sub-model was specified to use the MODFLOW transient head results and thus be consistent with the basin-scale model. In the sub-model, the GSFLOW model considered the time step of one day. The input data such as the precipitation and river discharge were used for observations and specified at a one-day time step. In the study, the groundwater flow based on the monthly time step was converted to the daily groundwater flow for the PRMS simulation. For practical use of the process, the groundwater flow is usually fixed in different months for the PRMS model. In the sub-model, the calibration and validation of the GSFLOW model relied on the observations of daily river discharge at the Li-lin and Wanda stations (please see Figure 1 for the locations). The simulation of GSFLOW requires the definition of the HRUs for the simulation domain. Figure 4b shows the specified HRUs in the sub-model. We are interested in quantifying the interflow obtained from the GSFLOW model. The results were used to compare the interflow potential map developed using the AHP technique. The calculation of the interflow in GSFLOW uses the following mass balance formula: where I [L 3 T −1 ] is the interflow and f [L 3 T −1 ] is the infiltration rate for a specific HRU. Variable ET [L 3 T −1 ] in Equation (3) is the evapotranspiration and S s [L 3 T −1 ] represents the soil moisture storage. In the GSFLOW model, the moisture storage comprises the capillary, gravity, and preferential-flow reservoirs. The variable R [L 3 T −1 ] is the groundwater recharge. Note that the groundwater flow processes, such as areal recharge, leakage from streamflow to aquifers, groundwater discharge to streams, and groundwater pumping, were simulated using the MODFLOW model. The unsaturated zone flow (UZF) package in MODFLOW is required to simulate the infiltration through the unsaturated zone. The package, associated with the PRMS module, can adequately account for the complex water interactions between the land surface and unsaturated zones above the groundwater levels.
Precipitation is the most important parameter that influences the interflow. We have found that the precipitation in the PPGB varies significantly in wet and dry seasons. Based on the calibrated GSFLOW model in the sub-model, the study further conducted stochastic modeling to assess the interflow uncertainty induced by the precipitation variations in the PPGB. The stochastic modeling used the Monte Carlo simulation (MCS) approach. In the MCS, realizations of precipitation time series for three-year were generated and applied to the calibrated GSFLOW model in the sub-model. There are 100 realizations generated for the assessment of the uncertainty induced by the precipitation. Figure 5 shows an example of the precipitation realizations at a station for three years. Note that the precipitation observations from 2015 to 2017 are the basis to generate the precipitation time series. The study assumed that the observations of the precipitation time series are the mean precipita-tion at each rainfall station. Random fluctuations based on a fixed standard deviation (STD) of 21.14 (mm/d) in a Gaussian model are used to generate the precipitation fluctuations. In this study, the STD value represents an averaged fluctuations of precipitation time series collected from the rainfall station in the PPGB. The Gaussian distribution model is used to generate the fluctuations of the rainfall dataset. The fluctuation time-series of three-year length were then added back to the mean precipitation at different stations for the modeling. The generated precipitation data at stations were used with the Thiessen-weighted method to obtain daily precipitation in the simulation domain. The spatial and temporal variations of the precipitation might not be stationary conditions. It could be interesting to assess the effect of similar precipitation fluctuations on the interflow variations in the wet and dry seasons in the PPGB.

Interflow Potential Based on the Index-Overlay Model
The interflow potential zones were determined by the weighted overlay method. The weighted overlay techniques in the ArcGIS software have been integrated into the package of spatial analysis. Based on the available raster data and the ratings and weightings obtained from the AHP technique, Figure 6 shows the results of the interflow potential for the downstream area of the Kaoping River. The map clearly shows that the high interflow potential zones are mainly distributed along the Kaoping River, and the potential varies from upstream to downstream. We found that the very low and low levels of interflow potential are in the hill area and some flat areas away from the main channel of the Kaoping River. In the hill area, the index very low level dominates the area. However, there are low-level zones in the hill area. The patterns of the low-level zones are similar to the patterns of the drainage density (i.e., the DD parameter), indicating that DD has improved the interflow potential in the hill zone. Note that the transmissivity (T) might control the groundwater flow. There are no hydraulic test data in the hill zone in the study, and the T map in the hill zone relies on the interpolated T near the Kaoping River. Therefore, the interflow potential might be overestimated in the hill zone because the T might be relatively small in the hill zone. In Figure 6, there is a clear interface between the hill zone and the river deposition zone in the north of the sub-model. This clear interface is consistent with the distribution of the parameters S, DH, and DG (see Figure 3). In this study, the values of inflow potential summarize the ratings and weightings obtained from the AHP technique. Slope has been recognized as the most influenced factor with the highest weighting to the interflow potential ( Table 2). Low DH and DG values indicate a better connection between surface water and groundwater. To use the water resource, the large water volume near the land surface could facilitate the development of water supply systems. In summary, the interflow potential in Figure 6 has shown the availability of the near-surface water resource along the Kaoping River. However, the quantification of the interflow potential might require a physically-based model that accounts for the interaction between surface water and groundwater, and the physically-based GSFLOW model assists efficiently to validate the developed interflow potential map and quantify the interflow in the sub-model of Kaoping River.

Quantification of Interflow Potential Based on the GSFLOW Model
In this study, the proposed numerical model was developed through the calibration and validation procedures. These procedures used the data obtained from the hydrogeological monitoring network in the PPGB. The calibration procedure used the data in the year 2015. However, the validation considered the data obtained from 2016 to 2017. Note that the study applied the MODFLOW model for the entire PPGB, i.e., the basin-scale model. However, the GSFLOW was applied to the sub-model in the downstream area of the Kaoping River. The transient solutions of the basin-scale groundwater flow were then used for the flow boundaries in the sub-model. The results and discussion are based on the model outputs covering the simulation time from January 2015 to December 2017.

Model Calibration and Validation
The calibration and validation of the GSFLOW relied on the groundwater and surface water observations in the sub-model. In the sub-model, the Li-lin station is the upstream controlling point for the river discharge observations. The downstream river discharge observations are from the Wanda station (see Figure 1). In addition, the length of the Kaoping River in the sub-model is approximately 35 km from the Li-lin station to the Taiwan Strait. Based on these classified characteristics, the study area was divided into 15 HRUs (see Figure 4b) to estimate the detailed interflow dynamics. Figures 7 and 8 show the results of the simulated river discharge at two different stations. The data in 2015 were used for the calibration, and the data from 2016 to 2017 were used for the model validation. The GSFLOW results and the observations are consistent in the three years, indicating that the developed GSFLOW can reasonably reproduce the system dynamics of the surface water and groundwater interaction. The variations of the river discharge in wet and dry seasons have been accurately captured. In Figures 7a and 8a, we found that the river discharge patterns at Li-lin and Wanda stations are similar. This is because of the relatively small drainage zones in the sub-model along the Kaoping River. Figures 7b and 8b further show the direct comparison of the observed and the simulated results. We obtained high correlations between the observed and the simulated river discharge.

Interflow Dynamics in the Sub-Model
In the GSFLOW model, interflow is one of the flow components that account for the net flow contributing to the stream. The time-varied calculation of the interflow in each HRU relies on Equation (3), which considers the water cycle components near the interface between the surface water and groundwater. In addition to surface runoff, interflow is the important flow that supplies the flow to the stream. Figure 9a shows the temporal variations of the interflow in the sub-model and the entire PPGB. The results of the interflow for the entire PPGB were proposed by Tran et al. (2020). The relatively consistent rainfall influenced similar temporal variations of the interflow in the PPGB. In the sub-model, the interflow rates are approximately 1/10 of the interflow rates in the PPGB in the dry season. In general, we obtain high interflow rates in the wet seasons and low flow rates in the dry seasons. The interflow rates in the sub-model can vary within two different magnitudes, respectively in wet and dry seasons. In addition, the interflow rates in wet seasons in the sub-model show a similar amount of interflow rates that are obtained from the entire PPGB (see Figure 9a). The similar interflow magnitudes in the basin and local models represent the important role of the interflow potential in the downstream area of the Kaoping River. In the wet seasons, the interflow in the downstream area of the Kaoping River could be an important water resource because there is always high turbidity of the river water. The water treatments for public use have become a challenging task.  In the sub-model, the high interflow locations are generally distributed in the upstream area of the sub-model. These areas have relatively high altitude, steepness, and permeability properties. In the GSFLOW results, we found that the highest interflow rate is approximately 4.5 × 10 4 (m 3 /d), corresponding to the volumetric flow per unit HRU area of 1.42 × 10 −3 (m/d). Interflow decreases in the middle and middle-lower zones along the Kaoping River in the sub-model. The interflow rates near the coastal zone are slightly higher than those in the middle and middle-lower zones. In the sub-model, the lowest value of the interflow rate is about 1.5 × 10 4 (m 3 /d), corresponding to 1.27 × 10 −3 (m/d) of the volumetric flow per unit HRU area. However, the variations of the averaged interflow rates were insignificant in the HRUs along the Kaoping River.
We have obtained the interflow potential based on the index-overlay method and the associated AHP techniques to modify the ratings and weightings. Figure 9b could be a reference to evaluate the accuracy of the developed interflow potential map in Figure 6. The interflow potential map in Figure 6 had shown that the highest interflow potential is near the Li-lin station and the interflow potential gradually decreases along the Kaoping River to the downstream area. The lowest interflow potential zone is near the Wanda station. Figure 6 also shows a slightly lower interflow potential near the coastal zone as compared with that near the Wanda station. The results shown in Figure 6 are consistent with the distributions of the averaged interflow rates obtained from the GSFLOW simulation (see Figure 9b). To quantify the interflow potential based on the GSFLOW results, we have estimated the averaged interflow for different HRUs. In Figure 9b, the highest average interflow rate is approximately 3.5 × 10 4 (m 3 /d) at high elevation zones in the sub-model. The value is 2.0 × 10 4 (m 3 /d) near the coastal zones in the downstream Kaoping River. In this study, the developed interflow potential maps based on the selected factors can generally represent the physically-based numerical model, indicating that the selected thematic factors and the associated AHP ratings and weightings are acceptable.
It is worth discussing here the contribution of the selected factors to create the interflow potential map. Note that the interflow potential map has integrated the hydrogeological parameters to evaluate the relative potential for the study area. In this study, intrinsic parameters such as the slope (S) could significantly influence the behavior of T, DH, DG, and DD near the main channel of the river (see Figure 3c,d). It is recognized that all the factors might be influenced by other factors based on different geological processes and different temporal and spatial scales. For the value of T, it is clear that the surficial deposit near the river also depends on the surface slope. A higher slope would produce depositional sediments with relatively larger grain sizes, which lead to higher T values for the aquifer system. In the study, the distribution of T in the hill zone might not be accurate because of the limited T observations in the hill zone (see Figure 3b). The DH, DG, and DD are the indicators that show the available water volume near the land surface. The high slope areas in the hill zone might obtain high ratings and weightings. However, the available water volume is limited because of the high DH and DG values in the hill zone. A similar concept could be applied to the condition of the parameter DD. The limited water volume (i.e., high DH and DG values) in the hill zone gives a relatively low interflow potential in the hill zone. We found that the factor DD can play an important role in modifying the interflow potential map based on the AHP technique. Figure 10 presents the results of the interflow potential when the factor DD is excluded. The result without the DD factor in Figure 10 clearly shows the relatively high interflow potential upstream when compared to the map in Figure 6. Specifically, the result in Figure 6 also shows more detailed variations of the interflow potential near the main channel of the Kaoping River than the results shown in Figure 10.

Interflow Uncertainty Induced by Precipitation
Based on the validated GSFLOW model, uncertainty analysis was conducted by using the concept of Monte Carlo simulations. One hundred realizations of precipitation time series were generated for the simulations (see Figure 5). In the study, we aim to quantify the influence of the input uncertainty (i.e., the precipitation) on the output uncertainty (i.e., the interflow). Figure 11 shows the variation of interflow STDs with the increased number of realizations. The result indicates a relatively small variation of the averaged STD for all the HRUs when the number of the realizations is greater than 50. Therefore, the 100 realizations in the Monte Carlo simulations could be representable to quantify the uncertainty of the interflow in the sub-model.  Figure 12 demonstrates the temporal variations of total interflow rates based on the precipitation realizations applied to the sub-model. Figure 12a shows the mean and realizations for the three-year simulations period. Figure 12b shows the variations of the mean and STD in a small period defined in Figure 12a (i.e., from December 2016 to April 2017) for demonstration purposes. In general, the variation of the precipitation in the PPGB has significantly influenced the variation of the interflow in the downstream area of the Kaoping River. With the high non-uniform precipitation in wet and dry seasons, the uncertainty of the interflow also varies significantly in wet and dry seasons. The uncertainty of the interflow is relatively higher in the wet seasons than that in the dry seasons. The result in Figure 12b has shown that a relatively stable variation of the interflow in the dry season from December 2016 to April 2017. The storm events such as those on December 6 to 8 and December 29 to 31 have clearly shown high variations based on the predefined mean precipitation in the model. Similar behavior is obtained in the wet seasons through the three-year simulations. However, the uncertainty is much higher in the storm events in the wet seasons. In this study, we previously calculated an STD of 21.14 (mm/d) for precipitation. Such mean variation of precipitation has produced the mean interflow STD of 1.26 × 10 4 (m 3 /d), corresponding to STD of 1.42 × 10 4 (m 3 /d) in the wet seasons and STD of 0.84 × 10 4 (m 3 /d) in the dry seasons. The results show that the overall fluctuation of the interflow rate in the wet seasons is slightly higher than that in the dry seasons, indicating a relatively stable interflow rate in the study area. Note that we have determined the mean precipitation based on the results of the validated GSFLOW model. The variation pattern of the precipitation might be constrained by the original observations of the precipitation in the PPGB. Further implementation could use a design rainfall pattern for the mean precipitation in the model.

Conclusions
This study developed an index-overlay and a numerical model to assess the interflow potential and quantify the interflow in the downstream area along the Kaoping River in southern Taiwan. The analytical hierarchy process was used to modify the ratings and weightings for the interflow potential map. The GSFLOW numerical model was then employed to quantify the interflow potential for practical applications. The rainfallinduced uncertainty was assessed to characterize the critical water cycle components in the study area.
The results of the interflow potential show that the high interflow potential zones are mainly distributed along the Kaoping River, and the potential varies from the upstream to the downstream areas. The very low and low levels of interflow potential are in the hill area and some flat areas away from the main channel of the Kaoping River. The study obtained a clear interface of interflow potential between the hill zone and the river deposition zone in the study area. This clear interface is consistent with the distribution of the parameter S, DH, and DG, indicating the better connection between surface water and groundwater.
In the study, the calibration and validation of GSFLOW were based on 15 specified HRUs and were used to estimate the detailed interflow dynamics. The results of the simulated river discharge at two different stations are fitted well, indicating that the developed GSFLOW can reproduce the system dynamics of the surface water and groundwater interaction. In the sub-model, the interflow rates are approximately 1/10 of the interflow rates in the PPGB. In general, the interflow rates are high in the wet seasons and low in the dry seasons. The interflow rates in the model can vary within two orders of magnitudes. The high interflow locations are generally distributed in the upstream area and the interflow rate decreases in the middle and middle-lower zones along the Kaoping River. However, the variations of the averaged interflow rates were insignificant in the HRUs along the Kaoping River.
We have obtained the interflow potential based on the index-overlay model and the associated AHP techniques to modify the ratings and weightings. The interflow potential map shows consistently the distributions of the averaged interflow rates obtained from the GSFLOW simulation. The highest average interflow rate is approximately 3.5 × 10 4 (m 3 /d) at high elevation zones, and the value is 2.0 × 10 4 (m 3 /d) near the coastal zones in the downstream Kaoping River. The interflow potential maps developed based on the selected factors can reproduce well the physical-based numerical model, indicating that the selected thematic factors and the associated AHP ratings and weightings are acceptable. The factor DD can play an important role in modifying the interflow potential map based on the AHP technique.
The uncertainty analysis was conducted by using the Monte Carlo simulations. With 100 realizations of precipitation time series, the study quantifies the influence of precipitation uncertainty on interflow uncertainty. In general, the variation of the precipitation in the PPGB significantly influences the variation of interflow in the downstream area of the Kaoping River. The uncertainty of interflow is relatively higher in the wet season than in the dry season. The stochastic simulation obtained the averaged interflow STD of 1.89 × 10 4 (m 3 /d) for the study area. From 2015 to 2017, the average STD was 2.35 × 10 4 (m 3 /d) for the wet seasons, and the average STD was 1.42 × 10 4 (m 3 /d) for the dry seasons. This is based on the calculated precipitation STD of 21.14 (mm/d). Such variation pattern of the precipitation might be constrained by the original observations of the precipitation in the PPGB. Further work could use a design rainfall pattern for the mean precipitation in the model.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.