Application of Geophysical and Hydrogeochemical Methods to the Protection of Drinking Groundwater in Karst Regions

To provide theoretical support for the protection of centralized drinking groundwater sources in karst areas, it is necessary to accurately identify the development of karst conduits and analyze the differences in hydrogeochemical characteristics of different karst systems. This provides a scientific basis for the accurate designation of risk zones that may cause drinking groundwater pollution. In this study, a geophysical survey, hydrogeological chemical process analysis and optimized fuzzy cluster analysis were used to gradually improve the understanding of karst water systems. AMT and HDR methods were used to calibrate the resistivity around the water-filling karst conduits, which ranged from 39 to 100 Ω·m. A total of seven karst systems were identified, including four karst systems in the north of the study area, one karst system in the west and two karst systems in the south. Analysis of the hydrochemical data showed that HCO3-Ca and HCO3-Mg-Ca types accounted for 90% of all samples. The δD and δ18O values of their main conduits were −51.70‰ to −38.30‰ and −7.99‰ to −5.96‰, respectively. The optimized fuzzy clustering analysis method based on the weight of variables assigned by AHP more accurately verified karst water systems. Based on these findings, the drinking groundwater source risk zone was designated with an area of 33.90 km2, accounting for 34.5% of the study area. This study effectively improved the rationality and accuracy of the designation of drinking groundwater source risk zones in karst areas, and provided a scientific basis for the identification of karst water systems and decision-making of drinking groundwater source protection in karst areas.


Introduction
Low temperature groundwater in karst development areas, which is mainly discharged in the form of karst springs or underground rivers, has become an important water source because of its stable quality and abundance [1][2][3]. However, the high heterogeneity of carbonate rocks in karst areas greatly increases the difficulty of pollution prevention and control [4]. Additionally, the strong hydraulic alternation conditions make this kind of water more susceptible to pollution and rapid migration, thereby posing health risks to the supply objects [5,6]. In view of the concealment of groundwater contamination and difficulty in remediation [7], prevention measures should be taken as the main means for groundwater protection when conducting industrial development, especially in karst areas [8][9][10]. Studies of heterogeneity, development degree and hydrogeochemical characteristics in typical karst aquifers can provide important information for sustainable development of groundwater and safety management of drinking groundwater sources [11,12], which is particularly important for developing the complex karst conduit systems and being the only drinking groundwater source for more than 15,000 people in the study area of this paper [13].
The geophysical survey is a qualitative or semi-quantitative exploration method based on obtaining differential information pertaining to electrical conductivity, electromagnetism or energy transmission, which is widely used studies pertaining to development and hydrogeological characteristics of a karst system [14,15]. McCormack et al. [16] overcame various difficulties in the application of traditional karst investigation techniques (such as the tracer method and the hydrogeochemistry analysis method) caused by a complex mixture of upland, lowland and coastal karst by using a high density electrical method, which is a common method for the construction of hydrogeological conceptual models, to identify the hydrogeological characteristics of a typical karst area. When compared with the high density resistivity method (HDR), which is suitable for media with a depth of less than 200 m [17,18], the audio frequency magnetotellurics method (AMT) can obtain resistivity information regarding deep media [19]. Wang et al. [20] used this method to successfully identify and forecast concealed karst conduits along a typical tunnel.
Hydrogeochemical methods, which investigate the distribution, migration and dispersion of chemical elements and isotopes in groundwater to interpret geochemical processes, are mainly applied to solve practical problems such as identification of groundwater sources and division of hydrogeological units [21,22]. Previous studies have shown that solutes reflect the hydrogeochemical background conditions of a study area, and are influenced by both natural and human activities [23]. Zhu et al., Sharma et al. and Chetelat et al. used multiple ion proportionality and chemical equilibrium analytical methods to analyze the relevant information pertaining to groundwater affected by precipitation input, the water-rock interaction, anthropogenic activity and other factors [24][25][26]. Based on analysis of the changes in characteristics of hydrogeochemical characteristics combined with isotopic information obtained along the groundwater runoff route and relative statistical methods, Leybourne et al. [27] and Yuan et al. [28] successfully interpreted the origin of groundwater and the dominant hydrogeochemical processes controlling groundwater water quality in a karst area.
Thus, this study was conducted to accurately identify karst conduits and effectively protect karst water systems that serve as drinking groundwater sources in a karst development area. AMT and HDR were used to identify the development of karst conduits, and groundwater samples were subjected to hydrogeochemical and isotopic analysis to identify the ionic component characteristics and hydrogeochemical process differences of groundwater in different karst systems. The optimized fuzzy cluster analysis method with weights assigned by variable indicators was then used to classify karst systems. The suitability of different karst systems in the area to serve as drinking groundwater sources was comprehensively studied and appropriate risk zones were designated according to the suitability of industrial activities.

Location and Hydrogeology of the Study Area
The study area was located in the hilly areas of the southeastern margin of the Sichuan Basin. It was defined based on hydrogeological conditions and karst development, corresponded to the recharge area of the Gu Song River upper reaches and had an area of 101.68 km 2 . Based on historical rainfall data, the rainfall is obviously controlled by seasonal factors, with 72% of the annual rainfall occurring from May-September, and only 8.8% occurring during winter.
According to the hydrogeology survey, karst depressions and sinkholes are densely distributed in the study area, and surface water systems are rarely developed. The only drainage datum plane, the Gu Song River, runs off from west to east across the study area. The geomorphology is obviously controlled by water erosion and accumulation as follows: the surrounding areas belong to the landform of a medium incision mountain in the shape of a spire. The corresponding emergence stratum is clastic rock strata of the Feixianguan and Xujiahe formation of the Triassic, with an exposed area of 30.84 km 2 , accounting for 30% of the study area. The emergence stratum in the center of the study area is carbonate of the Leikoupo and Jialinjiang formations of the Triassic, which have formed low Geological investigation of the karst development showed that the horizontal solution fissures in carbonate rocks were mainly developed along rock stratification, while the vertical solution fissures were mainly controlled by tectonic activities and vertical erosion of precipitation, and the maximum width of the fissure exceeded 1 m. Part of the solution fissures gradually eroded and developed into karst conduits or caves. The identified caves were 3-5 m high and 5-10 m wide, with a maximum length of several kilometers. According to data from 71 groups of karst springs observed in multiple field investigations in this study, the minimum mean flow and maximum mean flow of all springs were 0.1 and 1215.0 L/s respectively. The runoff process of groundwater in the study area was characterized by repeated alternation of surface runoff and underground runoff. It was not Geological investigation of the karst development showed that the horizontal solution fissures in carbonate rocks were mainly developed along rock stratification, while the vertical solution fissures were mainly controlled by tectonic activities and vertical erosion of precipitation, and the maximum width of the fissure exceeded 1 m. Part of the solution fissures gradually eroded and developed into karst conduits or caves. The identified caves were 3-5 m high and 5-10 m wide, with a maximum length of several kilometers. According to data from 71 groups of karst springs observed in multiple field investigations in this study, the minimum mean flow and maximum mean flow of all springs were 0.1 and 1215.0 L/s respectively. The runoff process of groundwater in the study area was characterized by repeated alternation of surface runoff and underground runoff. It was not feasible or necessary to divide the boundary of each spring drainage area. Therefore, this paper attempted to identify the main karst conduits by geophysical survey, and then conducted further hydrogeochemical study based on the identification results.

Geophysical Method
Based on the surface elevation and buried depth of the saturated zone of the karst aquifer, HDR and AMT were used to investigate the development of buried karst. In combination with the surface karst investigation, the developments of karst conduits in the study area were preliminarily identified. HDR adopted the Wenner four-electrode method of the electrical measurement system, with a current resolution of 0.01 µA and a current accuracy of 0.1%. Swedish software RES2DINV was applied to conduct the inversion of two-dimensional resistivity of the detection results. After topographic correction and removal of the discontinuity point, repeated parameter adjustments and iterative calculations were conducted until the root mean square error (RMSE) between the calculated resistivity and the measured resistivity was in the range of 3.6-8.79%. The reliable 2-D inversion resistivity profiles of HDR were obtained and used to interpret the buried depth of karst conduits. AMT adopted the MTU-5A GPS satellite synchronization magnetotelluric instrument produced by the Phoenix Company (Canada) and used its own AMT data interpretation software package for the inversion to obtain the electrical structure of underground media.
According to the hydrogeological investigation of surface karst and the application conditions of geophysical survey, HDR was mainly adopted in the middle of the study area with low altitude and shallow karst development depth. There were 20 survey lines (HDR-1-HDR-20), the distance between survey points was 5 or 7.5 m and the total length of the survey lines was 53.14 km, including 9435 survey points. AMT was mainly adopted south and northwest of the study area, where there was high altitude and deep karst development depth. There were six survey lines (AMT-1-AMT-6), the distance between survey points was 40 m and the total length of the survey lines was 25.84 km, including 684 survey points ( Figure 1).

Sampling and Laboratory Analysis
According to the development of karsts obtained by geophysical survey and the numerous groundwater outcrops distribution identified by surface hydrogeological survey, this study conducted optimal sampling based on the principle that the samples analysis results need to basically reflect the chemical characteristics of groundwater in the processes of recharge, runoff and discharge. As shown in Figure 1, 26 points, including springs, underground river outlets and water filled caves were selected as groundwater sampling points, and 4 surface water sections were selected as surface water sampling points. The selected sampling points of groundwater and surface water were sampled once in dry season (November 2018) and wet season (August 2019) respectively. In addition, a rain sample were collected in wet season. The collection, storage and monitoring of samples were strictly in accordance with the relevant standards of the American Public Health Association [29] and China [30,31]. Physical, chemical, metal, organic and inorganic factors were monitored. Table 1 shows the detailed information of each sampling point, and Table 2 shows the monitoring method and limit for each factor.   To further analyze the hydrogeochemical process and characteristics of the karst system in the study area, stable isotopes δ 18 O and δD were added as monitoring factors during the wet season. The δ 18 O and δD were measured using a stable isotope mass-spectrometer at the Key Laboratory of Karst Dynamics, Ministry of Natural Resources. The isotopic composition of oxygen was determined through the water-CO 2 equilibration technique and that of hydrogen through the water H 2 equilibration technique using a platinum catalyst. The reproducibility calculated from standards systematically interspersed in the analytical batches was ±0.2% for δ 18 O and ±1.5% for δD. The laboratory standards were regularly calibrated according to international standards (Vienna Standard Mean Ocean Water, VSMOW). The isotope ratios are reported in the denotation: where R is the molar concentration ratio of D to H or 18 O to 16 O. δ 18 O and δD were reported relative to VSMOW [32].

Calibration of Medium Resistivity
Controlled by the differences in mineral composition, porosity and groundwater occurrence, there are obvious differences in the resistivity between different rocks [33,34]. Therefore, to ensure accuracy of the identification, it was necessary to calibrate the resistivity of different types of rocks before conducting the geophysical survey in the study area to provide the basis for interpretation of karst conduits. HDR and AMT identification results showed that the resistivity of clastic rock with a buried depth of less than 40 m was between 10 2 -10 3 Ω·m because of the existence of the shallow weathered fracture aquifer, but the resistivity increased to more than 10 3 Ω·m as the buried depth increased.
In the carbonatite rocks area, the resistivity was between 251 and 1000 Ω·m from the vadose zone to the epiphreatic zone. Under the epiphreatic zone, the weak karst development area had a resistivity of more than 10 3 Ω·m, which was similar to the characteristics of clastic rock. Geophysical survey were carried out for several water-filled karst conduits identified by field hydrogeological investigation. The result showed that the resistivity of water-filled karst conduits was between 39 and 100 Ω·m. The unique conductivity of karst conduits, which was less than one order of magnitude than other identification areas, indicated the high heterogeneity of the carbonate area and ensured the accuracy of karst conduits identification.

Identification Results
Based on the results of resistivity calibration and considering the continuity of geophysical survey profiles and surface karst development, karst conduit 1-15 were identified in the study area ( Figure 1). Taking the Gu Song River as a reference point, karst conduit 1-6 were identified in the north, encompassing four underground river systems. Each underground river system in the study area is named after its underground river outlet and the names of the four underground river systems are Xin Zhai Spring (XZS), Long Xian Spring (LXS), Giant Salamander Spring (GSS) and Black Cave Spring (BCS). Karst conduit 7-10 were identified in the west of the Gu Song River, where the groundwater ran and drained through the Fish Well Spring (FWS). Karst conduit 11-15 were identified in the south, where there were two underground river outlets, Big Fish Spring (BFS) and Small Fish Spring (SFS). The detailed hydrogeological parameters of each karst conduit are shown in Table 3. The identification data in Table 3 shows that the development of karst conduit 1-6, which were included by XZS, LXS, GSS and BCS in the north, were similar, with an initial elevation of 851-967 m, outlet elevation of 454-514m, a possible main conduit length of 1.69-3.40 km and an possible average hydraulic gradient of 0.12-0.19. Karst conduit 14 and 15 included by SFS had initial elevations of 636-840 m, an outlet elevation of 459 m, the possible main conduit length of 5.73 km and an possible average hydraulic gradient of 0.07-0.04. It is worth noting that FWS and BFS are the centralized urban drinking water sources in the study area. The FWS, which is the drinking groundwater source for 3000 people, has a mean flow of 325.3 L/s during dry season and 813.2 L/s in the wet season, and karst conduit 7-10 within the FWS developing from the northwest to the southeast of the Gu Song River. The longest karst conduit in the area, 9, is possibly 7.21 km. The interpretation information of karst conduit 9 shows that the initial elevation is 1221 m, the buried depth of runoff area is 63-139 m and the elevation of the outlet is 528 m, as shown in Figure 2a. The BFS is the drinking groundwater source for 12,000 people, with a mean flow of 529.3 L/s in the dry season and 1215.0 L/s in the wet season. This area included karst conduit 11-13 developing from south to north. The longest conduit in the area 12 is possibly 7.21 km. The interpretation information for karst conduit (12) showed that the initial elevation was 887 m, the buried depth of the runoff area was 50-181 m and the elevation of the outlet was 475 m, as shown in Figure 2b. According to the function, FWS and BFS have become the focus of groundwater environmental protection in the study area. source for 12,000 people, with a mean flow of 529.3 L/s in the dry season and 1215.0 L/s in the wet season. This area included karst conduit 11-13 developing from south to north. The longest conduit in the area 12 is possibly 7.21 km. The interpretation information for karst conduit (12) showed that the initial elevation was 887 m, the buried depth of the runoff area was 50-181 m and the elevation of the outlet was 475 m, as shown in Figure 2b. According to the function, FWS and BFS have become the focus of groundwater environmental protection in the study area.

Basic Characteristics of Groundwater Quality and Chemistry
The collected water samples included rain, groundwater and surface water. The pH of rain was 7.56, the total hardness was 8.51 mg/L and the total dissolved solids content was 34.9 mg/L. The chemical types of groundwater and surface water were similar. According to data of groundwater (Table 4) and surface water collected during the dry and wet season, the pH value was between 7.02-8.46 and 7.15-8.25, respectively; the total hardness was 55-240 and 30-330 mg/L; and the total dissolved solids were between 71-276 and 71-504 mg/L. Na + , K + , Ca 2+ , Mg 2+ , Cl − , SO 4 2− , HCO 3 − and other factors accounted for more than 95% of the total dissolved solids in water samples. At the same time, all hydrochemical data had passed the ionic charge balance test of anions and cations, indicating that the measured data are reasonable and credible. According to the classification of Shoka Lev [35,36], the hydrochemical type of rainwater is HCO 3 -SO 4 -Na-Ca, while the chemical types of surface water and groundwater were similar, and mainly HCO 3 -Ca and HCO 3 -Mg-Ca, which accounted for 90% of the collected water samples. The low salinity and single hydrochemical type indicate that the circulation and alternation of groundwater are intensive and the rock type involved in the water rock interaction are not complex. In addition to the above factors, only Fe, Cd and Al were detected in the 10 metal factors, and the detection rate and concentration value did not fluctuate significantly during the dry and wet season. According to the data from the dry season, the detection rate and value of Fe were 56.7% and 3.1 × 10 −2 -6.8 × 10 −2 mg/L, respectively, while the detection rate and value of Cd were 46.7% and 1.9 × 10 −5 -1.0 × 10 −4 mg/L. Additionally, Al was detected in all water samples with a detection value of 8.0 × 10 −3 -9.8 × 10 −2 mg/L. No volatile phenols or petroleum were detected, and COD Mn value was 0.75-1.67. Three kinds of nitrogen compounds were detected in water samples, the detection rate of NO 3 − was 100% and the detection value was 3.11-14.37 mg/L during dry season and wet season. NO 2 − was only detected during dry season, with a detection rate of 6.7% and an average value of 0.27 mg/L. The detection values of NH 4 + were 46.7% and 60% in dry season and wet season, respectively, and the concentration was 0.017-0.19 mg/L. The concentration of nitrogen compounds in groundwater was higher during the wet season than the dry season, which may have been because of differences in the groundwater recharge intensity after rainfall leaching through surface materials between the dry season and wet season. Additionally, the strong hydraulic alternation conditions also increased the oxygen content in groundwater, resulting in the main nitrous oxide in groundwater being NO 3 − and almost no NO 2 − being detected. In general, 26 factors monitored in this study meet the drinking groundwater source standard of groundwater [37]. The concentrations of main components in groundwater were relatively low, especially for Cl − , SO 4 2− and NO 3 − , indicating that the input of exogenous substances caused by human activities is not the controlling factor of groundwater quality, while natural factors mainly control the solute composition and evolution of groundwater. Therefore, the primary measure of drinking groundwater source protection in the study area is pollution prevention and control of the groundwater environment.

Stable Oxygen and Hydrogen Isotopes
This section analyzes the differences between groundwater recharge sources based on hydrogen and oxygen stable isotopes. The δD and δ 18 O values of groundwater in the study area ranged from −51.70% to −27.60% and −7.99% to −5.15% , with average values of −40.60% and −6.74% , respectively. The values of δD and δ 18 O of the surface water ranged from −41.60% to −37.90% and −7.04% to −6.82% , with average values of −40.10% and −6.92% , respectively. The values of δD and δ 18 O of rain were −23.40% and −4.47% .
The differences in δ D and δ 18 O values were closely related to the water circulation process of the recharge source [38,39]. All spots of δ 18 O vs. δD of the karst water and surface water in the study area were located near the global metric water line (GMWL) [40], indicating that precipitation is the main recharge source of groundwater and surface water in the area ( Figure 3). It is not surprising that the distributions of δD and δ 18 O in the surface water and groundwater samples were similar, because they represent the strong hydraulic conductivity among precipitation, surface water and karst aquifer. Using the least square regression equation method to fit the δ 18 O vs. δD relationship of all water samples yielded δD = 7.33δ 18 O + 9.14. The slope of 7.33 was close to that of the GMWL equation δD = 8δ 18 O + 10, and the intercept of 9.14 was also close to that of the GMWL equation, indicating that δ D and δ 18 O are not affected by the secondary evaporation process in the study area. the distributions of δD and δ 18 O in the surface water and groundwater samples were similar, because they represent the strong hydraulic conductivity among precipitation, surface water and karst aquifer. Using the least square regression equation method to fit the δ 18 O vs. δD relationship of all water samples yielded δD = 7.33δ 18 O + 9.14. The slope of 7.33 was close to that of the GMWL equation δD = 8δ 18 O + 10, and the intercept of 9.14 was also close to that of the GMWL equation, indicating that δ D and δ 18 O are not affected by the secondary evaporation process in the study area.  The distribution of the spots of δ 18 O vs. δD in Figure 3 reveals the obvious difference between karst systems in the north and the FWS, BFS and SFS. The δD and δ 18 O values of the eight groups of water samples in the north karst system were between −43.50% and −36.20% and −7.23% and −6.12% , respectively. As shown in Figure 3, the distribution of the spots was relatively more concentrated, indicating that recharge and runoff areas are relatively small. The δD and δ 18 O values of the water samples of FWS, BFS and SF were between −51.70% and −38.30% and −7.99% and −5.96% , respectively, and the distribution range of δ D and δ 18 O values was obviously larger. In addition, there is a type of karst spring with small flow and typical characteristics of in-situ recharge and discharge, which has δD and δ 18 O values close to those of rain. Because of the rapid circulation of groundwater in karst voids, the water rock interaction between the groundwater and surrounding rock is weak. Sample GW5 is this type and the corresponding calcite saturation index (SIC) and dolomite saturation index (SID)are the minima of all samples, which were −1.551% and −4.075% , respectively.

Analysis of Hydrochemical Processes Affecting Groundwater Solute Components
Based on the hydrogeological conditions and the development of karsts obtained by geophysical prospecting, this study conducted optimal sampling and analysis of groundwater to further analyze the development of karst conduits by applying the chemical equilibrium analytical method and the multiple ion proportionality method. 5.1.1. Source Analysis of Cl − and Na + Cl − and Na + in groundwater are mainly derived from the input of precipitation and human agricultural activities. The molar ratio of Cl − to Na − is usually used to analyze the sources of Cl − and Na + in the study area [41]. The molar ratio of Cl − /Na + in the atmospheric precipitation input from the ocean is 1.16 [42]; therefore, the ratio of Cl − /Na + can be used to determine whether Cl − and Na + in the study area is from the atmospheric precipitation input. As shown in Figure 4, Cl − /Na + ranged from 0.03 to 1.36, with an average value of 0.58. The Cl − /Na + ratio of some water samples was close to 1:1.16, indicating that these water samples are heavily influenced by sea salt deposition. However, 90% of the samples had a molar ratio of Cl − /Na + less than 1.16, indicating that there are additional Na + ions involved in the chemical reactions in the study area. ocean is 1.16 [42]; therefore, the ratio of Cl − /Na + can be used to determine whether Cl − and Na + in the study area is from the atmospheric precipitation input. As shown in Figure 4, Cl − /Na + ranged from 0.03 to 1.36, with an average value of 0.58. The Cl − /Na + ratio of some water samples was close to 1:1.16, indicating that these water samples are heavily influenced by sea salt deposition. However, 90% of the samples had a molar ratio of Cl − /Na + less than 1.16, indicating that there are additional Na + ions involved in the chemical reactions in the study area. Carbonate rock strata were mainly distributed in the study area, while Triassic clastic rock strata were distributed in the hilly area on the outer edge of the study area, and the mineral compositions of the rocks distributed in the study area were mainly carbonate and silicate. No unique stratum containing Na + minerals was found in the study area. Therefore, the weathering and dissolution of carbonate and silicate strata were the main sources of Ca 2+ , Mg 2+ , Na + and K + . If it is assumed that exogenous acids such as sulfide oxidation are only used to balance Ca 2+   Carbonate rock strata were mainly distributed in the study area, while Triassic clastic rock strata were distributed in the hilly area on the outer edge of the study area, and the mineral compositions of the rocks distributed in the study area were mainly carbonate and silicate. No unique stratum containing Na + minerals was found in the study area. Therefore, the weathering and dissolution of carbonate and silicate strata were the main sources of Ca 2+ , Mg 2+ , Na + and K + . If it is assumed that exogenous acids such as sulfide oxidation are only used to balance Ca 2+

Carbonate Dissolution Analysis
The carbonate rock strata in the study area included the lower Triassic strata dominated by limestone and the middle Triassic strata with limestone and dolomite interbedding. To further explore the relative contributions of limestone and dolomite dissolution to the chemical ions of water in the basin, diagrams of Mg 2+ /Ca 2+ vs. HCO 3 − and Mg 2+ /Ca 2+ vs. SO 4 2− were drawn. The molar ratio of Mg 2+ /Ca 2+ can reflect the lithology of the carbonate aquifer through which the groundwater flows. When flowing through the limestone aquifer with calcite as the main mineral, the molar ratio of Mg 2+ /Ca 2+ was between 0.01 and 0.26, while when flowing through dolomite aquifer, the molar ratio will be greater than 0.85 [44]. As shown in Figure 6, the ratio of Mg 2+ /Ca 2+ of the samples ranged from 0.058 to 0.82, and the spots of Mg 2+ /Ca 2+ vs. HCO 3 − were basically distributed in the calcite and calcite-dolomite dissolution area. As shown in Figure 1, the karst system in the north developed in the lower Triassic aquifer dominated by limestone, and the spots were mainly distributed in the calcite dissolution area. The hydrogeological units of FWS, BFS and SFS, which were found to have a longer runoff path and wider recharge area, all contained the middle Triassic aquifer dominated by dolomite, of which the molar ratio of Mg 2+ /Ca 2+ in groundwater was positively correlated with HCO 3 − , which can indicate the water rock interaction intensity. In addition to the process of water rock interaction, the dissolution of dolomite will gradually increase after the dissolution of more soluble calcite, which is the main factor controlling the hydrochemical characteristics of FWS, BFS and SFS. The molar ratios of Mg 2+ /Ca 2+ showed no obvious correlation with SO 4 2− (Figure 7), indicating that the mineral source in groundwater is not closely correlated with sulfate dissolution.

Carbonate Dissolution Analysis
The carbonate rock strata in the study area included the lower Triassic strata dominated by limestone and the middle Triassic strata with limestone and dolomite interbedding. To further explore the relative contributions of limestone and dolomite dissolution to the chemical ions of water in the basin, diagrams of Mg 2+ /Ca 2+ vs. HCO3 − and Mg 2+ /Ca 2+ vs. SO4 2− were drawn. The molar ratio of Mg 2+ /Ca 2+ can reflect the lithology of the carbonate aquifer through which the groundwater flows. When flowing through the limestone aquifer with calcite as the main mineral, the molar ratio of Mg 2+ /Ca 2+ was between 0.01 and 0.26, while when flowing through dolomite aquifer, the molar ratio will be greater than 0.85 [44]. As shown in Figure 6, the ratio of Mg 2+ /Ca 2+ of the samples ranged from 0.058 to 0.82, and the spots of Mg 2+ /Ca 2+ vs. HCO3 − were basically distributed in the calcite and calcitedolomite dissolution area. As shown in Figure 1, the karst system in the north developed in the lower Triassic aquifer dominated by limestone, and the spots were mainly distributed in the calcite dissolution area. The hydrogeological units of FWS, BFS and SFS, which were found to have a longer runoff path and wider recharge area, all contained the middle Triassic aquifer dominated by dolomite, of which the molar ratio of Mg 2+ /Ca 2+ in groundwater was positively correlated with HCO3 − , which can indicate the water rock interaction intensity. In addition to the process of water rock interaction, the dissolution of dolomite will gradually increase after the dissolution of more soluble calcite, which is the main factor controlling the hydrochemical characteristics of FWS, BFS and SFS. The molar ratios of Mg 2+ /Ca 2+ showed no obvious correlation with SO4 2− (Figure 7), indicating that the mineral source in groundwater is not closely correlated with sulfate dissolution.

Analysis of Hydrochemical Characteristics along the Underground River Systems
Analysis of the composition and proportion of the main ions revealed that the karst conduit systems in the north of the Gu Song River were obviously different from those in the FWS, BFS and SFS. Considering the natural separation of the Gu Song River from the karst systems, this section focused on analysis of the hydrochemistry characteristics along the runoff path of the adjacent FWS and north karst conduit systems to further distinguish the differences between them.
The karst conduits in the north are densely developed, with short runoff paths and limited recharge areas. These characteristics result in small-scale development of branches in the north area, and the branches have no obvious impact on the water quality of the main conduits. As shown in Figure 8a, the concentrations of total dissolved solids (TDS), HCO3 − and other main ions in the north karst conduits are positively related to the runoff distance and negatively related to the elevation. The TDS was found to increase from 190 to 399 mg/L and the SIC from −0.0891 to 0.3024 from the recharge areas to the outlets. The carbonate rock strata along the karst conduits are dominated by limestone; therefore, even though the SID fluctuates, it was not found to be correlated with other elements. There are four karst conduits with a length of more than 2 km in the FSW, and the impact of branches on the water quality of the main conduits cannot be ignored. As shown in Figure 8b, from the recharge area of FWS to sampling point GW13 (elevation of 920 m), the concentrations of TDS, HCO3 − and other main ions were basically positively correlated with the runoff distance, and the peaks of TDS and HCO3 − were 358 and 220 mg/L, respectively. These findings were influenced by the water inflows from conduit branches, while TDS and HCO3 − were reduced to 237 and 149 mg/L, respectively, at the outlet. The elevation of the branch sampling points GW15 and GW16 was 710-

Analysis of Hydrochemical Characteristics along the Underground River Systems
Analysis of the composition and proportion of the main ions revealed that the karst conduit systems in the north of the Gu Song River were obviously different from those in the FWS, BFS and SFS. Considering the natural separation of the Gu Song River from the karst systems, this section focused on analysis of the hydrochemistry characteristics along the runoff path of the adjacent FWS and north karst conduit systems to further distinguish the differences between them.
The karst conduits in the north are densely developed, with short runoff paths and limited recharge areas. These characteristics result in small-scale development of branches in the north area, and the branches have no obvious impact on the water quality of the main conduits. As shown in Figure 8a, the concentrations of total dissolved solids (TDS), HCO 3 − and other main ions in the north karst conduits are positively related to the runoff distance and negatively related to the elevation. The TDS was found to increase from 190 to 399 mg/L and the SIC from −0.0891 to 0.3024 from the recharge areas to the outlets. The carbonate rock strata along the karst conduits are dominated by limestone; therefore, even though the SID fluctuates, it was not found to be correlated with other elements. There are four karst conduits with a length of more than 2 km in the FSW, and the impact of branches on the water quality of the main conduits cannot be ignored. As shown in Figure  respectively. The basin of FWS contains limestone and dolomite, and the variations in SIC and SID were found to be similar to those of TDS and HCO3 − : increasing and then decreasing. The gypsum saturation index (SIG) values and fluctuations of the two karst systems were relatively stable, ranging from −1.914 to −2.572 and −1.952 to −2.712, respectively. The features of SIG values and fluctuations, and the lack of an obvious correlation with the molar ratio of Mg 2+ /Ca 2+ and SO4 2− supported each other.
(a) (b) Figure 8. (a) North karst system hydrochemistry variation along the main conduit; (b) FWS karst system hydrochemistry variation along the main conduit and the influence of tributaries.

Optimized Fuzzy Cluster Analysis
Based on the above discussion regarding differences of hydrochemical and stable isotope data in different karst systems, the fuzzy clustering analysis method was used to classify karst systems quantitatively. Fuzzy cluster analysis is a mathematical method used to investigate and deal with classifications based on the variables of observation samples. This method can classify items that lack reliable historical data, but have similar properties, into one group [45,46]. Fuzzy cluster analysis is widely used in geological exploration, water pollution, pattern recognition and other fields [47]. In traditional fuzzy cluster analysis, it is assumed that each variable has the same control over the results, while in the karst area the weights of different variables such as hydrochemical factors and stable isotopes are obviously different. Based on the characteristics of hydrogeochemical classification in karst area, this study attempted to assign the weight to variables and obtain the optimized fuzzy clustering method to quantitatively analyze the differences between each karst system.

Optimized Fuzzy Cluster Analysis
Based on the above discussion regarding differences of hydrochemical and stable isotope data in different karst systems, the fuzzy clustering analysis method was used to classify karst systems quantitatively. Fuzzy cluster analysis is a mathematical method used to investigate and deal with classifications based on the variables of observation samples. This method can classify items that lack reliable historical data, but have similar properties, into one group [45,46]. Fuzzy cluster analysis is widely used in geological exploration, water pollution, pattern recognition and other fields [47]. In traditional fuzzy cluster analysis, it is assumed that each variable has the same control over the results, while in the karst area the weights of different variables such as hydrochemical factors and stable isotopes are obviously different. Based on the characteristics of hydrogeochemical classification in karst area, this study attempted to assign the weight to variables and obtain the optimized fuzzy clustering method to quantitatively analyze the differences between each karst system.

(2) Standardization
Different variables may have different dimensions or orders of magnitude. To eliminate the above differences and meet the requirements of the fuzzy matrix analysis, the data were standardized and compressed to the interval of 0 to 1. It is usually necessary to change the translation standard deviation first: x ik = x ik − x k s k (i = 1, 2, · · · , n; k = 1, 2, · · · , m;) After the change, the mean value of variables in the sample was 0, the standard deviation was 1, and the influence of dimension and the order of magnitude between variables was eliminated. The changed x ik may not be in the interval of 0 to 1, and the translation-range transformation 0 ≤ x ik ≤ 1 is made, which is suitable for the construction of a fuzzy similar matrix. The translation-range transformation formula is shown in Equation (2), where x ik replaces x ik for clustering analysis.
(2) Construction of similarity matrix The traditional clustering method considers the variables in the sample to have the same importance to the sample, and constructs the similarity matrix by calculating the similarity between x i and x j . The similarity r x i , x j is generally characterized by the Euclidean distance method: The importance of different variable factors in groundwater in karst areas is obviously different. The weight (ω k ) refers to the quantitative allocation of the importance of different indicators of the analysis samples, and the different treatment of the role of each variable in the clustering process. In the process of multi-index evaluation and analysis, ω k has the function of highlighting key indexes, making the multi-index structure reasonable and realizing the overall optimization. Under the condition of considering the influence of different variables on the similarity degree, the Euclidean distance formula is changed into [48]: (3) Assignment of weights At present, the determination of ω k has entered the stage of combining qualitative and quantitative analysis. By introducing mathematical methods, the theoretical model of determining weights has been developed and improved. There are two kinds of weight determination: subjective and objective [49]. The subjective weighting method is mostly based on the knowledge or experience of experts or individuals, and adopts the qualitative method of comprehensive consulting scoring to determine the weight, and then synthesizes the standardized data. The commonly used methods include the comprehensive index method, the Delphi method, and an analytic hierarchy process [50,51]. The objective weighting method is used to determine the weight according to the correlation between each index or the variation degree of each index value. The objective weight is determined from the data obtained from the survey; therefore, there is no need to consult experts. Common methods of determining this weight include principal component analysis, factor analysis, the variation coefficient method and determination of complex correlation coefficients [52,53]. In this study, the Delphi method and AHP method were used to determine the weight of each index. A group of experts in environmental engineering, hydrogeology and other fields was also employed to score the weights of variables in the sample, and the improved AHP method was used for calculation and analysis. The weight assignment steps were as follows: (1) Experts rank the importance of evaluation indicators, x 1 > x 2 · · · x m ; (2) Experts determine the rational assignment r k of the ratio of the importance degree of adjacent indexes x k−1 and x k .
where for r k = 1, x k and x k−1 have the same importance; for r k = 1.2, x k is slightly more important than x k−1 ; for r k = 1.4, x k is significantly more important than x k−1 ; for r k = 1.6, x k is strongly more important than x k−1 ; where r k = 1.8, x k is extremely more important than x k−1 . In this analysis, 1. According to the rational assignment r determined by the experts, the weight of the m evaluation index ω m can be calculated as: According to the weight ω m , the weight calculation formula of the m−1, m−2, . . . 3, 2, 1 indicators is as follows: (4) Output of cluster analysis chart The fuzzy matrix constructed by the distance obtained from the optimized calculation formula is a fuzzy similarity matrix R. Because it is not necessarily transitive, the intra group connection method is used to achieve its transitivity [54]. In the fuzzy matrix, all samples in the sample set U are classified and gradually merged to form a dynamic clustering diagram.

Sample Analysis
Based on the results of physical survey and hydrogeochemical data analysis, the outlets of the karst system in the north area, FWS, BFS and SFS, and some karst springs distributed in their runoff area were selected as the sample set. The results showed that the solute composition of groundwater in the study area is mainly controlled by water rock interactions, and hardly affected by exogenous materials. The main dissolved minerals were carbonate, followed by silicate. K + , Na + , Ca 2+ Figure 9 shows the results of traditional and weighted optimization cluster analysis. Overall, it can be divided into the north karst system in the study area represented by the outlet samples GW03, GW02 and GW06,;the karst systems with long runoff paths and large scales represented by the outlet samples GW25, GW20 and GW14; and rain, which shows the maximum distance from the outlets of each karst system. The samples located in the recharge-runoff area were between rain and outlet points. The samples collected from near the watershed in areas with high elevation had properties similar to those of rain, while the samples collected from near the outlets had properties similar to the outlets. Comparison of the analyses conducted using the two methods revealed that the optimized method produced a classification that is more hierarchical and recognition that is more accurate. Additionally, in the optimization calculation, samples GW21 and GW18 were identified and belonged to the runoff area of karst area represented by GW25, GW20 and GW14, and they are obviously different from other samples belonging to the runoff area of karst area represented by GW03, GW02 and GW06. These findings objectively reflect the characteristics of the groundwater cycle process. points. The samples collected from near the watershed in areas with high elevation had properties similar to those of rain, while the samples collected from near the outlets had properties similar to the outlets. Comparison of the analyses conducted using the two methods revealed that the optimized method produced a classification that is more hierarchical and recognition that is more accurate. Additionally, in the optimization calculation, samples GW21 and GW18 were identified and belonged to the runoff area of karst area represented by GW25, GW20 and GW14, and they are obviously different from other samples belonging to the runoff area of karst area represented by GW03, GW02 and GW06. These findings objectively reflect the characteristics of the groundwater cycle process.  Table 5 shows the Euclidean distances of water samples. The optimized Euclidean distance was 0.039-0.267, and the maximum distance was observed in GW03 and the rain sample. In addition to rain, the maximum distance of groundwater samples was 0.255, which indicates the difference of the recharge and runoff area between BFS represented by GW17 and north karst systems represented by GW05. Additionally, GW03, GW02 and GW06 represent the north karst systems and have a Euclidean distance of 0.042-0.068, while GW25, GW20 and GW14 represent SFS, BFS and FWS, respectively, and have a Euclidean distance of 0.043-0.056. The distance among the north karst systems, FWS, BFS and SFS, was 0.086-0.118. These findings are supported by the results of geophysical survey and hydrogeochemical data analysis.    Table 5 shows the Euclidean distances of water samples. The optimized Euclidean distance was 0.039-0.267, and the maximum distance was observed in GW03 and the rain sample. In addition to rain, the maximum distance of groundwater samples was 0.255, which indicates the difference of the recharge and runoff area between BFS represented by GW17 and north karst systems represented by GW05. Additionally, GW03, GW02 and GW06 represent the north karst systems and have a Euclidean distance of 0.042-0.068, while GW25, GW20 and GW14 represent SFS, BFS and FWS, respectively, and have a Euclidean distance of 0.043-0.056. The distance among the north karst systems, FWS, BFS and SFS, was 0.086-0.118. These findings are supported by the results of geophysical survey and hydrogeochemical data analysis.

Groundwater Environmental Risk Zone Identification
This study defined the karst conduits' development area and recharge-runoff area of the drinking groundwater source together as a risk zone. In the risk zone, human activities that may produce pollution, such as chemical industry, mining, centralized livestock and poultry farming will be restricted.
The results of geophysical survey, hydrochemical data and cluster analysis clearly show that FWS and BFS differ from the north karst systems. However, no clear boundary was found between conduit 6 and 7 or conduit 13 and 7. Therefore, in view of the safety of the drinking groundwater source, the risk zone was designated with the recharge-runoff area of FWS and BFS; conduit 6 and conduit 14 are the boundaries, giving an area of 33.90 km 2 and accounting for 34.5% of the study area ( Figure 10). In the risk zone, the industrial activities that may have an impact on the groundwater environment should be prohibited or restricted.

Groundwater Environmental Risk Zone Identification
This study defined the karst conduits' development area and recharge-runoff area of the drinking groundwater source together as a risk zone. In the risk zone, human activities that may produce pollution, such as chemical industry, mining, centralized livestock and poultry farming will be restricted.
The results of geophysical survey, hydrochemical data and cluster analysis clearly show that FWS and BFS differ from the north karst systems. However, no clear boundary was found between conduit 6 and 7 or conduit 13 and 7. Therefore, in view of the safety of the drinking groundwater source, the risk zone was designated with the recharge-runoff area of FWS and BFS; conduit 6 and conduit 14 are the boundaries, giving an area of 33.90 km 2 and accounting for 34.5% of the study area ( Figure 10). In the risk zone, the industrial activities that may have an impact on the groundwater environment should be prohibited or restricted.

Conclusions
This study carried out a detailed investigation based on the objective of drinking groundwater source protection in typical karst areas in the southwest of China. Geophysical survey, hydrogeological chemical process analysis and optimized fuzzy cluster analysis were applied to conduct a systematic study on the development of karst conduits. Each step employed in this study was mutually dependent and corroborated. Results of each step gradually enhanced the understanding of karst system, and effectively improved the rationality and accuracy of drinking groundwater source risk zone designation of karst aquifer. Finally, 34.5% of the study area was designated as risk zone. In the risk zone, controlled by the high gradient between 0.04-0.23 of karst

Conclusions
This study carried out a detailed investigation based on the objective of drinking groundwater source protection in typical karst areas in the southwest of China. Geophysical survey, hydrogeological chemical process analysis and optimized fuzzy cluster analysis were applied to conduct a systematic study on the development of karst conduits. Each step employed in this study was mutually dependent and corroborated. Results of each step gradually enhanced the understanding of karst system, and effectively improved the rationality and accuracy of drinking groundwater source risk zone designation of karst aquifer. Finally, 34.5% of the study area was designated as risk zone. In the risk zone, controlled by the high gradient between 0.04-0.23 of karst conduits system, once the pollutants enter the groundwater system, they will quickly migrate to the drinking groundwater source. There is not enough response time to carry out remediation to prevent the drinking groundwater source from being polluted. Therefore, in order to avoid the human health risk caused by the drinking groundwater sources pollution, results of this study can be submitted to the local government. The result will be an important scientific basis for the construction of local regulations on the protection of drinking groundwater sources. The regulations need to specify that human activities, which may produce pollutants or have a negative impact on the groundwater environment, shall be prohibited or restricted in the risk zone.