Integrating In-Situ Data and RS-GIS Techniques to Identify Groundwater Potential Sites in Mountainous Regions of Taiwan

Due to rapid urbanization, the development of megacities and metropolises worldwide is creating water scarcity, social-environmental risk, and challenges to the regions where water supply from rivers and alluvial aquifers is insufficient and unstable. Groundwater exploration in fractured bedrock of mountainous regions is thus a crucial issue in the search for substitute water resources. To achieve cost effectiveness on groundwater exploration, the use of comprehensive remote sensing (RS)and geographic information system (GIS)-based models appears feasible. The required parameters selected and analyzed from the literature depend on the hydrogeological characteristics. This study intends to investigate and improve the proposed parameters and data sources upon those presented in the literature. A total of 17 hydrogeological units of concern was delineated from 105 complex geological formations of the geological sections and main rock types. The other parameters related to groundwater potential were derived from the digital elevation model and Landsat imagery. In addition, 118 drilling cores were inspected and in-situ well yield data from 72 wells were employed to assess the normalized groundwater potential index in the raster-based empirical GIS model with a higher spatial resolution. The results show that the accuracy of the interpretation of groundwater potential sites improved from 48.6% to 84.7%. The three-dimensional (3D) visualization of a thematic map integrated with satellite imagery is useful as a cost-effective approach for assessing groundwater potential.


Introduction
Water resources are indispensable to mankind, society, and countries and have a tremendous impact on people's livelihood as well as national agriculture, industry, and economy. When the water supply from rivers and alluvial aquifers is insufficient and unstable, it is crucial to explore substitute water supply from unexploited areas, such as mountainous regions [1][2][3][4][5][6]. However, most mountainous regions are situated in geologically complex terrains (GCTs) with heterogeneous hydrogeological features. For example, Taiwan is located in the subtropical monsoon climate zone and arc-continent collision (ACC) between the Luzon volcanic arc and the Eurasian continent [7][8][9]. Due to the rugged terrain and uneven spatial and temporal distribution of rainfall, Taiwan's annual precipitation is about 2.5 times higher than the world's average, but the water supply available to the people per capital is only 20% of the world's average. Thus, exploration of groundwater potential in such a geological environment is a challenge. The high cost of pointwise exploration of groundwater resources and inconvenience of conducting on-site tasks in remote areas are often noticed in engineering practice. Thus, an alternative approach to achieve cost-effective groundwater exploration, and therefore identify groundwater potential sites (GWPS), is of great importance. Integrating remote sensing (RS) and geographic information system (GIS) is a possible alternative [10][11][12][13][14][15][16][17].
This study lists case studies in bedrock, and three categories of parameters favorable to groundwater (i.e., geology, topography, and remote sensing index), as shown in Table 1 [1,5,12,[18][19][20][21][22][23][24]. The literature provides useful information on the characteristics of regional hydrogeology and geomorphology. However, certain case study areas were delineated by the boundaries of the administrative area or random rectangle. They were mostly located in terrains with simple geological units and slope degrees. At these local scales, the high-ratio units of each parameter may have more influence on the results, such as rock types and slope degrees; even the adjacent fracture network of bedrock outside the study area was neglected. It needs to consider the hydrogeological model at the scale of the entire drainage system and the regional groundwater flow system. Additionally, most data sources are in vector format with a scale of 1:50,000 or less, such as geological maps and landform maps; certain data can be analyzed as raster files from given spatial resolutions of the digital elevation model (DEM) and remote imagery. In order to achieve consistency of data sources, a raster-based empirical GIS model is recommended to calculate the thematic map and validate it with investigated data.
In this paper, the above issues are solved at the scale of catchment management in the GCTs of Taiwan, where the seven categories of terrain were defined by topographic position index (TPI), topographic wetness index (TWI) and slope degree (SD), and the 17 hydrogeological units (HGUs) were delineated from digital geological maps. In addition, the regional lineament and abnormal surface temperature favorable to groundwater were considered and derived from RS. Comparing insitu well yield data, the accuracy of the interpretation of GWPS varied with rock types and terrains. However, assessing the normalized groundwater potential index (NGPI) in the comprehensive RS-GIS model improved the results to achieve a cost-effective method to identify GWPS.

Regolith and Fractured Bedrock
According to the hydrogeological conceptual model (Figure 1), there are two major layers (i.e., regolith and fractured bedrock) in the mountainous regions. Regolith materials overlaid on bedrock are defined by the conditions of weathering, from slightly weathered to residual soil [25,26]. The specific rock types are named in the literature [27][28][29]. In this study, the average depth of regolith was 15.9 m. There were three kinds of materials in the regolith layer inspected in 118 sites (100 m in depth) and samples. The first kind consists of soil, backfill, alluvium, and colluvium as inspected by the degree of weathering, sphericity, angularity, and roundness in Holocene. The other two kinds were saprolite and saprock, weathered from bedrock with the weathering between slight and high degrees, respectively. Faintly weathered and fresh bedrock belongs to fractured bedrock ranging from Eocene to Pleistocene.
The hydrogeological settings are mainly dominated by the lithological properties of the selected study area where porosity or secondary porosity (e.g., fracture) can provide space to store and trap water [30][31][32][33][34][35]. This information is included in the map of the hydrogeological unit to provide a conceptual hydrogeological model and characteristics [36][37][38][39].

Hydrogeological Properties
In terms of experimental hydraulic parameters (Figure 2), the storage, permeability, and porosity of regolith are generally higher than those of bedrock [27,[40][41][42]. However, groundwater in bedrock is distributed in fracture networks and tends to exist and flow at secondary weak planes or weathered zones with higher permeability [33,34,[43][44][45].
The effects of surface topography on the groundwater flow system are described and simulated by theoretical and numerical methods [46,47]. There are three types of groundwater flow systems (i.e., local, intermediate, and regional flows) affected by topographic relief in a basin. Surface topography also generally reflects the spatial distribution of soil moisture and groundwater levels in the process of water infiltration, recharge, and discharge. The properties of SD, river gradient, flow accumulation (FA), drainage density (DD), and TWI can be analyzed from a DEM. The regional geomorphology, land cover or land use, groundwater recharge, and water bodies can be determined by ground investigation or derived from RS [48][49][50][51][52]. Furthermore, geological lineament structures provide information on regional connectivity of a fracture network, where the groundwater potential could be higher [53][54][55]; in addition, spatiotemporal land surface temperature (LST), normalized difference vegetation index (NDVI), and soil moisture index (SMI) provide dynamic evidence on groundwater recharge and discharge [56][57][58].
Therefore, hydrogeological characteristics change with the evolution of regional geology and topography. The target of groundwater exploration in this study was discharge areas that are composed of regolith and fractured bedrock aquifers. It becomes feasible to determine remotely sensed surface characteristics linked to groundwater by using the proposed RS and GIS model.

The Study Area
The island of Taiwan ( Figure 3) covers an area of 35,873 km 2 , and its central geographic coordinates are 23°58′ N and 120°58′ E. In Taiwan, 25 central government-governed rivers originate in the steep topography, which has an elevation ranging from sea level to 3952 m. The island is 394 km long and 144 km wide. The geological sections [7,59] from east to west are delineated by major tectonic structures and faults, including the Coastal Range (CR), Longitudinal Valley (LV), Taroko and Yuli Belt (TY) of the East Central Range, the Hsueshan Range (HS) and the Backbone Range (BR) of the West Central Range, the Western Foothills (WF), and the Coastal Plain (CP). In the evolution of plate tectonics in Taiwan, the Tananao Schist complex is composed of a metamorphic pre-Tertiary basement of the Eurasian passive margin, which has evidence of past orogenic events in recent metamorphic events. For instance, the TY belt was strongly deformed and metamorphosed in the Mesozoic period. Due to the Cenozoic arc-continent collision and orogeny located on the boundary of Eurasia and the Philippines Sea plates, the slate belt in the HS and BR is connected with the Chinese passive continental margin, the CR corresponds to the accreted Luzon arc, and the LV is considered a suture zone between the arc and continental margin. In the western foreland basin, the WF at lower altitudes has been accreted and deformed with syn-orogenic sediments, and the CP is part of the present foreland basin. As a result of this orogenic movement, the terrain of Taiwan consists of approximately 1/3 plains and 2/3 mountain regions ( Figure 3). The plain areas (e.g., CP and LV), where the formations of unconsolidated rocks are divided into aquifers and aquitards according to the hydrogeological characteristics, are delineated into 9 major groundwater areas and are the locations of most of the major cities. In the mountain regions (e.g., WF, HS, BR, TY and CR), the formation of consolidated rock mainly includes sedimentary, metamorphic, and volcanic rocks ranging from Eocene to Pleistocene with complex fracture networks of folds and faults.
Due to the limited number of investigations to date, the central and southern mountainous regions of Taiwan were chosen as the study area. The area is about 19,098 km 2 in size, and it is 280 km long and 110 km wide. The eastern and western watersheds originate from the Central Range, below the elevation of 3952 m, and include 11 rivers: the Dajia, Wu, Jhuoshuei, Pachang, Zengwun, Gaoping, Sizhong, Hoping, Hualien, Xiuguluan, and Beinan rivers. The 118 drilling sites and 72 groundwater wells were distributed in different topographic terrains and in each geological section of the WF, HR, BR, TY, and CR, as shown in Figure 3.

Materials
The primary data used in this study consist of: (1) A digital geological map with a scale of 1:25,000 established by the Central Geological Survey (CGS), including distributions of geological formations, faults, and folds.

Methodology
To achieve the goal of cost-effective groundwater exploration and recognize the influence of the likely parameters on the ACC, the unique features related to groundwater potential were configured and described by integrating RS and GIS techniques with favorable groundwater parameters in the same spatial resolution of 30 m × 30 m, including lineaments and SMI derived from Landsat imagery, and SD, DD, and TWI analyzed from DEM. HGU was considered as the parameter of lithology to simultaneously reflect the complex spatial distribution of geological structures and rock types in the ACC. Therefore, six specific parameters, namely, HGU, LD, SD, DD, TWI, and a variation of the SMI (VSMI) were selected to compute the NGPI in the ArcGIS Spatial Analyst module using Equation (1). The values of ranking and weight may change in each case study. To avoid the subjective effect on each parameter, its equal weight and appropriate 5 classifications of ranking were used in this study. That is, the weighting of each parameter was designed to be the same value of 1. In addition, the ranking of each parameter was based on the low to high groundwater potential and computed by the scale value of 0-5. Except for the parameter of HGU, the data formats of the other 5 parameters were digital numbers to determine the best arrangement of values by Jenks Natural Breaks (JNB) classification method [60,61] into 5 designed rankings. Therefore, it is helpful to understand the groundwater potential of each selected parameter in the RS-GIS model. Finally, the NGPI with 6 proposed parameters in each pixel was computed with the values of 0-30 and normalized from 0 to 1 in the final synthesized map of groundwater potential. The results were evaluated according to insitu well yield data from 72 wells.
where Pi is the ranking of each parameter given by the scale value of 0-5, and Wi is the weight of each parameter given by the same value of 1 in this study.
As illustrated in the study framework in Figure 4, the selected parameters and products include: (1) generating a geological map with a scale of 1:25,000 established by the CGS and delineating HGUs to complete the post-produced hydrogeological map; (2) deriving LD and multiple temporal SMI from Landsat imagery; and (3) analyzing DEM to calculate SD, DD, and TWI. The processing of the parameters (Table 2 and Figure 5) and in-situ well yield data are introduced in Appendix A.  0.00~0. 23 Lower potential 1 1 The proposed six parameters were hydrogeological units (HGUs), lineament density (LD), slope degree (SD), drainage density (DD), topographic wetness index (TWI), and variation of the soil moisture index (VSMI), After assigning the ranking of HGUs in the vector format, they were reclassified as a raster format. The other data formats were digital numbers to determine the best arrangement of values by Jenks Natural Breaks classification method [60,61] into 5 designed rankings.

Results
According to the method described, six group variables, namely, HGU, LD, SD, DD, TWI, and VSMI ( Figure 5), were computed in the RS-GIS model. The final synthesized map with NGPI was produced. In addition, the correlations between groundwater potential, NGPI, and the findings are described below.

HGU Thematic Map and Distribution of Well Yields
Delineation of HGUs is a primary task of hydrogeological investigations. For each case of an HGU thematic map, the information depends on the scale and purpose of the study. In the ACC of Taiwan, delineation of HGUs by geological sections and rock types simplified the complex geological units into 17 HGUs and provided the general geological and hydraulic properties. In addition, limited in-situ investigation in each HGU outlined the groundwater potential, as shown in Figure 6. The 72 wells investigated were located in the WF, the Central Range, and CR. Results showed that most average well yields of the HGUs belonged to WFSS in the Western Foothills, and HRQS, HRSA, HRSL, BRSP, and TYSC in the Central Range. However, the HGUs of WFSS, BRSP, and TYSC had large variability in their well yield values as well as the lithologic units in the Lawrenceville of Georgia [62]. Other dominant factors also control groundwater potential, such as fracture and geomorphology [7,12]. In contrast, WFSH and HRAR had lower groundwater potential due to the composition being mostly mud and less-fractured; consistently the well yield of 30 L per minute was distributed in the shale and slate [20]. This finding strongly supports the correct delineation and ranking of HGUs in the ACC. The trend of groundwater potential provided the primary information and references for the candidates of water supply.

LD Thematic Map
Tectonic movement may cause numerous weak planes, fractures, and lineaments near the geological structures of faults and folds. In the LD thematic map (Figure 7), fault and fold lines were distributed around the areas with high LD, such as near the boundary of the WF and Central Range, the east part of the Central Range, and the middle of the CR. These results indicated that the given parameters of lineament extraction in the PCI-LINE module were appropriate. The values of RADI, GTHR, LTHR, FTHR, ATHR, and DTHR were tested and adjusted according to a case study in the literature (Table A1). Linking the LD thematic map to groundwater potential, some in-situ sites with higher well yields were located in the areas with higher LD; however, some of those near the boundary of the WF and Central Range, where the rock types are deposited with more mud matrix, did not match the expected results. In the fold-bend fault region of the WF, the mud matrix of Miocene strata might fill the fracture space and block the groundwater flow from reaching connected fractures. Thus, a negative correlation was found between LD and groundwater potential in the WF of the ACC, and the local mud matrix should be considered and examined in future explorations, especially in the WF area.

Geomorphic Thematic Map (TWI, SD, and DD)
Since groundwater may flow in accordance with variations of topography, the TWI is useful to describe the drainage system and runoff (Figure 8). Class 1 is distributed near the boundary of a watershed, class 2 belongs to slope flow, and class 3 is the downstream of slope flow. These flows may accumulate in the areas of classes 4 and 5, identified as streams and branches. The drainage calculated DD with the threshold of FA is distributed within classes 4 and 5 of the TWI. It helps to rank higher values in the sources of basins (i.e., lower DD) as the recharge areas on a regional scale. In addition, five classes of SD provide the general categories of terrain on a local scale. Therefore, the geomorphic thematic map was developed by overlaying the TWI, SD, and DD thematic maps to present preliminary information on groundwater occurrences.

Seasonal SMI Thematic Map
To understand the correlation between the variation of soil moisture and groundwater potential, the SMI thematic maps of two wet and dry seasons and their standard deviation statistics were compared within the regions with different well yield levels ( Figure 9). The results indicated that the seasonal trends in each year were similar, with the exception of specific regions, such as the B100-21 site. This site was located near a region with higher variation of SMI, and its well yield was comparatively higher than that of the B100-23 site. Therefore, SMI derived from Landsat was found to be a reliable source of information on groundwater potential.

Classification of Groundwater Potential
In the above sections, all the selected parameters were discussed separately from well yields. Before mapping of groundwater potential, the four categories of GWPS (Table 3) were defined by the values of well yields [63]. Groundwater sites with well yields of more than 600 L per minute are suitable as regional water supplies for humans and irrigation. Those with yields between 60 and 600 L per minute are considered local water supplies. Therefore, a GWPS is defined as having a well yield of greater than 60 L per minute; otherwise, it belongs to non-GWPS where groundwater is insufficient for water supply. According to the well yield data, the 72 groundwater sites were classified into two groups, GWPS and non-GWPS, as shown in Table 4, whose percentages were 48.6% and 51.4%, respectively. Table 3. Classifications of groundwater potential and groundwater potential sites (GWPS) (after [63]).

Degree of Groundwater Potential
Well Yield (L/min) Description GWPS High Potential >600 Regional water supply for humans and irrigation Yes Medium Potential 60-600 Local water supply for humans and irrigation  Average value 203.9 -15.9 --0.60 1 A GWPS is defined as having a well yield of greater than 60 L per minute; otherwise, it belongs to non-GWPS where groundwater is insufficient for water supply. Abbreviations: near roof (Rf), at ridge (Rg), steep slope (Ss), flat slope (Sf), valley or creek bottom (Vb), alluvial fan of downstream the valley (Vbf), main riverbed deposit (Vm).

Final Synthesized Map of Groundwater Potential in the ACC
The resulting synthesized map of groundwater potential is shown in Figure 10. NGPI values were classified into four levels by JNB as low (0.20~0.50 NGPI), medium (0.25~0.56 NGPI), high (0.56~0.63 NGPI), and very high (0.63~0.97 NGPI) groundwater potential. The percentages of the low, medium, high, and very high potentials were 26%, 32%, 27%, and 15%, respectively. NGPI values, which range from 0 to 1, were compared with the well yield data, as shown in Table 4. There was a positive correlation between the NGPI and well yields.

Application of the NGPI to Predict Groundwater Potential and Verification of In-Situ Data
Since the very high level of NGPI was above 0.63, the threshold of the NGPI for identifying GWPS and non-GWPS was determined to be 0.63. To apply and inspect the NGPI for groundwater exploration, the well yield data of 48 sites in central Taiwan were employed as training data. The value achieved an accuracy of 83.3%. According to previous experience, 24 sites in southern Taiwan were predicted. The estimated results achieved an accuracy of 87.5%, while that for all the tested sites was 84.7% (Table 5). On the other hand, using a threshold of 0.63 to predict more well yields (higher than 300 or 600 L per minute) as GWPS resulted in lower accuracy. Therefore, the definitions of GWPS were met by the degree of 0.63~0.97 NGPI (very high) and verified with in-situ data. The verified NGPI value of 0.63 is a significant threshold for identifying local and regional water supplies in the follow-up groundwater exploration in the ACC.

Proposed Parameters of Spatial Resolution
Although GIS spatial analyst tool provides the fast result of a thematic map of weighted sum with proposed parameters, the result may be less useful due to the inconsistency of data sources and spatial resolution. However, the use of a finer grid for pointwise exploration is expensive and may still not be able to determine in areas where traffic cannot be reached. Under the issue of achieving the goal of a cost-effective method in remote areas by integrating RS and GIS techniques, the global Landsat imagery and ASTER GDEM can be obtained from the National Aeronautics and Space Administration (NASA) website with a uniform spatial resolution of 30 m × 30 m. The proposed parameters can be produced from the raw data. Thus, a raster-based empirical GIS model is recommended to calculate the thematic map of groundwater potential and validate it with investigated data.

Geological Characteristic of Groundwater Potential
For inspection of the possibility of groundwater potential in each geological section and rock type, Table 6 shows the percentage of GWPS at investigated sites. The geological section of the Central Range (e.g., HR, BR, and TY) is mostly composed of metamorphic rock (e.g., slate, schist, phyllite, and quartz) with a higher groundwater potential than that of WF (e.g., alternations of sandstone and shale). For the rock types of mud, shale, and argillite, the ingredients of very fine sand and mud caused the frequency of groundwater potential to be lower, whether they belonged to metamorphic rock or not. In addition, the percentages of medium, high, and very high levels of groundwater potentials in each geological section were higher in the Central Range ( Figure 11). These findings correspond to the high possibility of GWPS in the metamorphic rock of the ACC ( Figure 6). Thus, appropriate rankings and selected parameters of HGU and LD are sufficient to support the information on fracture networks affecting the groundwater potential in the active ACC area of Taiwan.  Figure 11. The percentages of WF, the Central Range, and CR in low to very high levels of groundwater potential.

Topographic Characteristic of Groundwater Potential
Revised from a previous study [64], the seven categories of terrain are defined as near roof (Rf), at ridge (Rg), steep slope (Ss), flat slope (Sf), valley or creek bottom (Vb), alluvial fan of downstream the valley (Vbf), and main riverbed deposit (Vm), as shown in Table 7. To illustrate those terrains, the TPI, TWI, and SD were implemented. As for SD, classes 1 and 2 belong to Ss, classes 3 and 4 represent Sf, Vb, and Vbf, and class 5 is Vm. Furthermore, Vb and Vbf can be extracted by the class 4 of TWI. Finally, the topography tool of ArcGIS was used to derive the TPI and delineate the geomorphic terrain from ridge to valley (Figure 12) where the areas with TPI values above 28 were defined as Rf and Rg. Within the data from the 72 investigated wells, less data on the Rf, Rg, and Ss terrains were available due to the difficulty of drilling and testing. The lower well yields may be located in those terrains as well as Vb. However, other terrains with higher percentage of GWPS included Vm, Vbf, and Sf, the average regolith depths of which were thicker in sequence. In addition, there was a positive correlation between regolith depth and well yield. Thus, for groundwater exploration based on topographic characteristics, the possibility of GWPS was about 44.4% to 76% (Table 7). With integration of RS and GIS, the possibility increased to 84.7% (Table 5). Therefore, another key parameter of GWPS is the characteristics of geology and remote sensing index (e.g., the proposed parameters of HGUs, LD, and VSMI in this study).

Visualization of the NGPI to Identify Groundwater Potential
The visualization of three-dimensional (3D) results provides information about potential sites [65]. However, in agricultural eras preceding modern computer analysis, people relied on the characteristics of terrain and their personal experience of groundwater exploration. For instance, the Chinese proverbs related to the locations of GWPS describe the confluences of two rivers, intersections of two mountains, concave sides of riverbeds at the foot of mountains, and so on. It is not easy to understand these descriptions of groundwater potential from the ancient proverbs without comparing the 3D visualization of the NGPI map with a Google map (Figures 13 and 14). At site B103-02, the distribution of higher NGPI values shows a pattern consistent with the confluence of two upstream rivers. Site B102-08, located in a downstream deposited alluvial fan, is shaped into the mouth where two mountains intersect, and is also illustrated by a higher NGPI value. Site B106-02, located on meandering terrain and the concave side of a riverbed, is full of groundwater due to reduction of the flow velocity and increases in recharged water resources. Those predominant types of terrain are favorable to GWPS. Oppositely, an area composed of nearly horizontal bedding has less groundwater flow, accumulation, and potential. In order to solve the problem of water scarcity and assess the substitute water resources, the regions located in geological units of metamorphic rock, and topographic units of Vm, Vbf, and Sf have more groundwater potential for regional water supply. In addition, the RS data-and GIS-based results provide visualized and comprehensive information for groundwater exploration, which improved the accuracy of the interpretation of GWPS from 48.6% to 84.7% (Tables 5 and 7).    Figure 13). Without comparing the 3D visualization of the NGPI maps with a Google map, it is difficult to link the results of GWPS with these descriptions of groundwater potential from the ancient proverbs: (a) more groundwater at confluences of two rivers (at B103-02 GWPS); (b) more groundwater at intersection of two mountains (at B102-08 GWPS); (c) more groundwater at concave side of riverbed and foot of mountain (at B106-02 GWPS); (d) less groundwater at nearly horizontal bedding of rocks (at B099-15 non-GWPS).

Uncertainty of Groundwater Potential and Suggestion
In the terrain and geology of the ACC, groundwater flow of fractured bedrock is complex. To achieve the goal of regional groundwater exploration and sustainability, it is useful to integrate RS and GIS techniques to identify the groundwater occurrence in places where the geological profile of the tested sites may consist of regolith and fractured bedrock aquifers within 40 m depth. However, relying on those parameters of surface information entails some uncertainty. For example, regolith depth changes greatly with variation in terrain, but external collapse deposits may occur and accumulate more groundwater. In addition, an area with a high LD represents a well-developed fracture network under the movement of the ACC, which leads to brittle fractures and open fractures favorable to groundwater flow in metamorphic rock. In contrast, it causes unfavorable conditions of GWPS due to the high number of shear planes or shear zones filled with mud and less open fracture in sedimentary rock. Thus, it is possible to have no groundwater in areas of high groundwater potential, while it is impossible to expend further effort to conduct an entire in-situ investigation. To figure out the key parameters of groundwater potential, the possibility of GWPS in a zone of metamorphic rock is about 62.2%, and except for argillite, those in each rock type composed of metaphoric rock are about 50 to 80% (Table 6). In the geomorphic thematic map considering SD, DD, and TWI, the possibility of GWPS is about 44% to 76% (Table 7). Geological properties may have greater influence on GWPS. Thus, the parameters of HGU, LD, and SMI were implemented in the RS-and GIS-based models with equal rankings, values, and weighting, as well as the topographic parameters of SD, DD, and TWI. The results showed that the accuracy of identifying GWPS reached 84.7%. The parameters of HGU, LD, and SMI derived from hydrogeological mapping and multitemporal Landsat imagery increased the predicted result. They are important factors to affect the groundwater potential in the ACC. In future work, the proposed method may be applied to other test areas over a large scale. In addition, future work can increase the precision of HGU mapping, LD analysis, DEM models, and RS imagery, as well as the quantity of in-situ data, to illustrate groundwater flow clearly and reduce its internal uncertainty. Therefore, it is suggested that finding the key impact factors initially by integrating RS and GIS data for targeted area is needed on the issue of groundwater exploration.

Worldwide Potential Applicability
Global hydrogeological characteristics change with the development of geology and geomorphology. The proposed RS and GIS method was implemented in GCTs, as a case study for tectonic movement of Taiwan. On the basis of in-situ investigation, the key parameters conducive to groundwater were determined. Similar regions worldwide include Japan, the Philippines, the western and eastern United States, western South America, western Asia, southern Europe, and so on. The adjustment of parameters obviously needs to be considered. For example, the ranking of slate rock type is lower in the ancient platform due to lower groundwater storage and poor permeability, but it is higher in the developed fracture network of Taiwan. In this case, the weathering zone (i.e., soil and regolith) may be one of the selected parameters because it has a greater impact on the shallow groundwater potential than the rock type of bedrock. By applying the developed method worldwide and discussing the results, it will be useful to understand the impact parameters of groundwater potential.

Conclusions
Groundwater hydrogeological characteristics vary in different geological regions, leading to uncertainty in mapping groundwater potential with integrated RS and GIS. Literature showed the precise dominated parameters in each specific area on this issue. In this study, a new approach was used to identify the key parameters of groundwater potential in the ACC of Taiwan. For example, HGUs delineated by geological sections and main rock types of geological formations were considered as the parameters of lithology; lineaments were extracted by the PCI-LINE module from Landsat imagery with a high spatial resolution; and remote seasonal detection of SMI was used to identify the dynamic evidences on groundwater recharge and discharge. Integrating geomorphic thematic maps by overlaying TWI, SD, and DD thematic maps produced the final synthesized map of the NGPI. Results showed positive correlations between in-situ well yield data and the NGPI, and the possibility of groundwater potential mapping. Therefore, the comprehensive RS-and GIS-based model appears to be a useful tool for the follow-up groundwater investigation. The predicted accuracy also achieves the goal of developing a cost-effective method for exploring complex geological and topographic landforms.

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

Appendix A Appendix A.1. Lithology/Hydrogeological Unit (HGU)
Lithology is one of the important and recognized parameters for mapping of groundwater potential in bedrock. Each rock type has a specific hydraulic property that affects the characteristics of groundwater. Due to the active collision in Taiwan, the study area contains many geological events and formations. To simplify and reflect the geological and hydraulic properties, 17 HGUs (Figure 5a) were considered as the parameter of lithology, and they were delineated by the geological sections and main rock types of geological formations. The geological sections (Figure 3b) were originally identified by the boundaries of regional geological structures and formations in typical geological times and events. The distributions from east to west in Taiwan were the CR, LV, TY of the East Central Range, the HS and BR of the West Central Range, the WF, and the CP. The geomorphologic features, such as the plains area (e.g., CP and LV), piedmont (e.g., WF), and mountainous region (e.g., HR, BR, TY and CR), were also reflected. Additionally, each rock type involves a range of porosity, hydraulic conductivity, storage, and so on. Thus, HGUs could provide primary quantification of the groundwater potential in these complex geological and topographic areas. To rank the value of each HGU, the maximum values of geological section and main rock type were distributed to two, for a total of five. Since higher porosity, permeability, and hydraulic conductivity are favorable to groundwater flow, three categories of groundwater potential were ranked in the following order: porous medium of unconsolidated rock (e.g., CP and LV), fractured bedrock (e.g., HR, BR, and TY), and sedimentary rock and breccias (e.g., WF and CR). Rock types with higher porosity and welldeveloped fracture networks (i.e., primary and secondary weakness of the bedrock), such as gravel, alternations of sandstone and shale, quartz, slate, phyllite, schist, gneiss, and volcanic breccias, were assigned higher rankings and values; mud, shale, argillite, and marble were assigned lower rankings and values. Details of the rankings and values of HGUs are listed in Table 2.

Appendix A.2. Lineament Density (LD)
Due to the uncertain distribution of geological structures, the geological lineament structure directly derived from RS data was employed for the GIS model. Lineaments provide information on regional connectivity of a fracture network, where the groundwater potential could be higher, as well as lithology. In order to extract the lineament, the Landsat imagery was processed with enhanced image processing and edge filtering techniques. In addition, the parameters related to the line features in the software were used to achieve automatic lineament extraction. Ultimately, post processing of the extracted lineaments with ground-based data was required.
(1) Enhanced image processing For completion of the radiometric calibration and atmospheric correction on selected highquality images with less clouds and noise, PCA was used to enhance the extraction of geological information. Since the PC-1 images contain most of the information as they are based on maximum variance, they were implemented as the enhanced images and the line features had clear edges.
(2) Automatic lineament extraction and parameter test Based on the above-mentioned PC-1 images, edge filtering and line extraction techniques were performed in MATLAB or RS software. The most widely used software tool is the LINE module of PCI Geomatica. The parameters of concern were the radius of the filter in pixels (RADI), threshold for edge gradient (GTHR), threshold for curve length (LTHR), threshold for line fitting error (FTHR), threshold for angular difference (ATHR), and threshold for linking distance (DTHR). The parameters can determine the lineament forms and circular (e.g., curve) structures. The parameter tests were based on the given bedding (e.g., layer of stratum) and extracted lineaments. Thus, compared with the geological structures of a geological map (Figure 7), the results of this study and the literature are shown in Table A1. (3) Post-processing of verification and density map Verification by in-situ data was needed to exclude artificial linear features, such as roads, power grids, rivers, and image backgrounds. In addition, automatically extracted lines with lower RADI values had more details (e.g., short lines) and noise, and those with DTHR values made the extracted lines linked with each other. Thus, in the 30 m × 30 m spatial resolution Landsat imagery, extracted lines longer than 1000 m were considered as geological lineaments in post processing. Furthermore, the lineament density (LD) was analyzed with the spatial tools of ArcGIS software. The calculation of LD was defined by Equation (A1), and the results were ranked by JNB. A higher ranking and value of LD indicates the conditions of weathering or vulnerability of the bedrock, which present more potential for groundwater flow and storage.
where Li is the length of the lineament and A is the unit of area.

Appendix A.3. Topographic Wetness Index (TWI)
For the identification of assumed groundwater flow by topography, the TWI was analyzed from a DEM with 30 m × 30 m resolution in the ArcGIS Raster Calculator tool using Equation (A2), with the properties of SD and FA included. The significant distribution patterns of the TWI thematic map with five categories ranked by JNB provided basic information for groundwater exploration, and higher TWI values indicated the discharge areas, riverbeds, and flat slope areas with more potential for regolith depth and water content.

. Slope Degree (SD)
To sketch the regions of local and regional groundwater systems, SD and DD ranked by five categories of JNB were additionally overlaid in the TWI thematic map. SD was expressed as the angle of inclination from a horizontal plane, and the maximum rate of change between each cell and its neighbors was calculated with the GIS slope tool.

Appendix A.5. Drainage Density (DD)
For DD, suitable drainage was extracted with the GIS spatial analyst tool. Since drainage systems are produced as the terrain develops, the detail patterns were similar to the distribution of higher TWI values. The major drainage in each watershed was extracted by a threshold value of FA. The maximum, mean, and standard deviation values of FA were 1,461,926, 1017, and 23,682, respectively, in this study area. Testing results showed that a lower threshold value of FA led to the display of more creeks, and higher stream ordering was defined. In this discrimination, the DD thematic map was not able to show the favorable distribution of groundwater potential clearly. Ultimately, a suitable threshold FA value of 8046 was presented and inspected with the method of stream ordering proposed by the literature [65], resulting in five stream orderings, as shown in Figure 3b. The analyzed results of drainage were used to calculate the DD map by the spatial tools of ArcGIS software. In a drainage basin, the calculation of DD, defined by Equation (A3), is the total length of all river segments divided by the total area of the basin.
where Li is the length of the drainage and A is the unit of area.
Integrating the TWI, SD, and DD maps into a geomorphic thematic map, the areas simultaneously with lower SD and DD values indicated flatter terrain with higher permeability, less runoff, and recharge potential. Therefore, surface runoff would decrease, and downward infiltration would increase. Therefore, groundwater potential was significantly distinct in the topographic conceptual model of the geomorphic thematic map.

. Soil Moisture Index (SMI)/Variation of SMI (VSMI)
The uncertainties of mapping groundwater potential integrating RS and GIS included complex seasonal groundwater fluctuations and insufficient underground investigation. Thus, RS played an important role in detecting the regional surface appearances of groundwater flow systems, such as green biomass (i.e., NDVI), LST, and soil moisture. Soil moisture is a direct indicator of subsurface water that is found in the unsaturated zone above the water table. Rainfall and groundwater flow may contribute to soil moisture content. Seasonal detection of SMI is useful for identifying the variation of groundwater potential. Therefore, multiple thematic maps of SMI derived from Landsat imagery were implemented in two wet and dry seasons by using data from 9 June and 2 December 2015, and 27 June and 4 December 2016. SMI was defined as Equation (A4) based on the scatter plot of LST and NDVI forming a trapezoid in the LST-NDVI space [66][67][68][69][70]. Consequently, the standard deviations in the four seasons of SMI results were compiled as statistics for the parameter of VSMI.
where Tmax and Tmin are the maximum and minimum surface temperature for a given NDVI, respectively. These values were obtained by a linear regression of known remotely-sensed data for both dry and wet edges in the LST-NDVI space.