Hydrochemical and Isotopic Difference of Spring Water Depending on Flow Type in a Stratigraphically Complex Karst Area of South Korea

Characterizing the subsurface flow in karstic areas is challenging due to distinct flow paths coexisting, and lithologic heterogeneity makes it more difficult. A combined use of hydrochemical, environmental isotopic, and hydrograph separation study was performed to understand the subsurface flow in a karst terrain where Ordovician carbonate rocks overlie Jurassic sandstone and shale along thrusts. Spring water collected was divided into Type Ⅰ (n = 11) and Ⅱ (n = 30) based on flow patterns (i.e., low and high discharge, respectively). In addition, groundwater (n = 20) was examined for comparison. Three Type Ⅱ springs were additionally collected during a storm event to construct hydrographs using δ18O and δD. As a result, Type Ⅱ had higher electrical conductivity, Mg2+, HCO3 −, and Ca2+/(Na+ + K+) than Type Ⅰ and was mostly saturated with calcite, similar to deep groundwater. The hydrochemical difference between Types Ⅰ and Ⅱ was opposite to the expectation that Type Ⅱ would be undersaturated given fast flow and small storage, which could be explained by the distinct geology and water sources. Most Type Ⅱ springs and deep groundwater occurred in carbonate rocks, whereas Type Ⅰ and shallow groundwater occurred in various geological settings. The carbonate rocks seemed to provide conduit flow paths for Type Ⅱ given high solubility and faults, resulting in 1) relatively high tritium and NO3 − and Cl− via short-circuiting flow paths and 2) the similar hydrochemistry and δ18O and δD to deep groundwater via upwelling from deep flow paths. The deep groundwater contributed to 83–87% of the discharge at three Type Ⅱ springs in the dry season. In contrast, Type Ⅰ showed low Ca2+ + Mg2+ and Ca2+/(Na+ + K+) discharging diffuse sources passing through shallow depths in a matrix in mountain areas. Delayed responses to rainfall and the increased concentrations of contaminants (e.g., NO3 −) during a typhoon at Type Ⅱ implied storage in the vadose zone. This study shows that hydrochemical and isotopic investigations are effective to characterize flow paths, when combined with hydrograph separation because the heterogenous geology affects both flow paths and the hydrochemistry of spring water passing through each pathway.

In the karst area, groundwater recharge occurs by either diffuse infiltration through the fissured matrix, called diffuse-dominated, or direct point infiltration into the conduit network, called conduitdominant. The diffuse-dominated flow is associated with slow and delayed responses of water seepage from the aquifer matrix, less connected fractures, or overlying soils, whereas the conduitdominated flow is fast through the interconnected voids and passages of the aquifer (Schilling and Helmers, 2008). Theoretically, the diffuse laminar flow is distinguished from the turbulent conduit flow based on flow velocities (e.g., a few centimeters or meters per day for diffuse flow in Atkinson, 1977; a few m/h to several hundred m/h for conduits in Rusjan et al., 2019) or Reynolds numbers (e.g., less than 2,000 for laminar flow in Grasso et al., 2003;Ghasemizadeh et al., 2012; over 100 for conduit flow in Lee et al., 2006). The diffuse flow in general is assumed to undergo water-rock interactions, increasing the contents of dissolved ions (e.g., Adji et al., 2016), whereas the conduit-dominant flow is habitually undersaturated with respect to calcite (e.g., De Rooij and Graham, 2017). Of them, conduit flow makes the karst aquifer vulnerable to contamination through direct and fast infiltration (Stueber and Criss, 2005;Ghasemizadeh et al., 2012;Parise et al., 2015) or to rock collapse due to open caves (Waltham and Lu, 2007;Pogačnik et al., 2017). Thus, an understanding of conduit flow paths is essential when managing water resources or building facilities (e.g., dams) in the karst area (Milanović, 2021).
Subsurface flow in karst aquifers is mainly characterized by monitoring the hydrological (e.g., discharge), hydrochemical (e.g., cations and anions), and isotopic (mainly, δD, δ 18 O, and δ 13 C) variability of spring water that is groundwater naturally emerging to the surface (e.g., Lee and Krothe, 2001;Aquilina et al., 2005;Barbieri et al., 2005;Aquilina et al., 2006;Doctor et al., 2006;Moore et al., 2009;Jiang et al., 2013;Lorette et al., 2018;Gil-Márquez et al., 2019;Rusjan et al., 2019;Palcsu et al., 2021). For instance, Moore et al. (2009) identified water sources (e.g., deepwater upwelling) and their mixing and interactions with aquifer materials based on hydrochemical data to understand the relationship between spring characteristics and the complexities of karst aquifers. Doctor et al. (2006) used the δD and δ 18 O of water and δ 13 C of dissolved inorganic carbon (DIC) as well as major ion chemistry to estimate mixing proportions among multiple sources, and found that water released from storage within the epikarst may comprise as much as two-thirds of conduit flow in a karst aquifer following rainfall. Palcsu et al. (2021) found the 10% of the discharge water with a residence time of half a year and the contribution of an older component based on the seasonality of δ 18 O in 2015-2018 at a karst spring with an average discharge of 20 L/min in Hungary. Lorette et al. (2018) found high Mg 2+ and low dissolved oxygen (DO) during low-water periods and explained the uncommon hydrochemistry using the mixture with a Jurassic confined aquifer in western France.
In those studies, geochemical modelling and saturation index (SI) computation were often applied to interpret the monitoring data (e.g., Barbieri et al., 2005;Peyraube et al., 2012) because the source of acidity for carbonate dissolution and the geochemical conditions for geochemical facies as well as the impact of human activities can be evaluated by examining the relationships between solution concentrations, partial pressure of CO 2 (P CO2 ), and SI of calcite (SI calcite ). Also, principal component analysis (PCA) was extensively used to infer hydrogeochemical processes or water sources (e.g., Doctor et al., 2006;Moore et al., 2009;Mudarra et al., 2012;Lorette et al., 2018;Gil-Márquez et al., 2019). Mudarra et al. (2012) included P CO2 and SI calcite in PCA.
Another method to understand the groundwater flow in a karst area is an examination of spring responses to rainfall to quantify the amount of rainwater (or baseflow) contributing to the discharge. Hydrograph separation is performed using various data including hydrochemical parameters (e.g., NO 3 − , Cl − , and specific conductance by Schilling and Helmers (2008)) and stable isotopes (e.g., δD and δ 18 O by Aquilina et al. (2005), Klaus and MacDonnell (2013)) to apportion the discharge among multiple water sources. Vesper and White (2004) noted that karst spring hydrographs reflect not only the changing mix of baseflow and storm flow, but also a shift in recharge sources during a storm event. Grasso et al. (2003) coupled the analysis of karst spring hydrographs and chemographs and evaluated the geometric dimensions of submerged karstic networks.
We noted that the combined use of hydrochemical and isotopic compositions and hydrograph separation with the help of PCA would be effective to describe the subsurface flow in karstic areas as in Gil-Márquez et al. (2019), who used hydrodynamics, hydrochemistry, and environmental isotopes to devise a hydrogeological conceptual model of evaporitekarst springs. In particular, the multi-approach is expected to help delineate flow paths when the karstic flow system is intricate with lithologic heterogeneity. Besides, spatial studies were rarely conducted for the hydrochemical and isotopic differences of spring water depending on flow types (e.g., low and high discharge), whereas the temporal variation of discharge and hydrochemistry as a function of the hydrological regime was performed by many researchers, including Bicalho et al. (2012) Adji et al. (2016, Lorette et al. (2018), Gil-Márquez et al. (2019), Frank et al. (2019), and de Filippi et al. (2021. Therefore, we conducted hydrochemical, environmental isotopic, and hydrograph separation studies to assess the subsurface flow in a stratigraphically complex karst area with thrusts where both low discharge and high discharge springs occurred ( Figure 1). It was assumed that the low discharge and high discharge come through diffuse and conduit pathways, respectively, because the flow in karst aquifers is turbulent in conduits such as caves and sinkholes, while diffusive through the fissured matrix or intergranular unfractured bedrock (Demiroglu, 2016;Chang et al., 2019). Specifically, springs were divided into Type Ⅰ and Ⅱ when sampled based on the manner by which the discharge occurred (i.e., low and high, respectively) as in Mocior et al. (2015) who surveyed the spatial distribution of springs in mountains and categorized them by discharge. Then the hydrochemical and isotopic properties of spring water were characterized depending on the flow type and compared with those of groundwater to understand the flow path. Then baseflow separation was conducted for three Type Ⅱ springs by monitoring δ 18 O and δD data during a storm event to estimate the proportion of baseflow in the conduit flow-dominated spring discharge. Lastly, the subsurface flow and the storage of water and solutes in the vadose zone were discussed.
The contribution of this work is the characterization of groundwater flow systems in a stratigraphically complex karst aquifer with thrusts using the spatial distribution of hydrochemical and isotopic data with statistical and hydrographic separation analysis in South Korea which has few studies on karstic aquifers. It should be noted that a multipurpose dam was planned to be built on the Dong River on the Bansong Formation (Jbs and Jbc; Figure 1) in the early 1990s for water supply, flood management, and hydropower generation, but the plan was canceled in 2000 because of environmental issues. The field work for this study was conducted during the feasibility study of the dam since karst structures in heterogeneous lithological units may cause geohazards for civil infrastructures (Pogačnik et al., 2017). The data has not been properly interpreted until this study because of social sensitivity. Although the data are >20 years, the data interpretation is expected to provide information on the subsurface flow in this stratigraphically complex karst aquifer given the geologic time scale and no significant land use change over 20 years in the area (Supplementary Figure S1).

Geography and Climate
The study area is located in the Yeongwol-Danyang basin at the middle eastern part of South Korea (Figure 1). The topographic elevation ranges from 200 to 800 m above sea level (asl). In terms of geomorphology, high and steep mountains are located within the study area, and flat terranes are developed on a mountain ridge at high altitude (Supplementary Figure S2), indicating repeated cycles of uplifting and erosion in this site over geological time. The Dong River meanders along narrow and deep valleys formed between high and steep mountains ( Figure 1A). The land along the river has been used for agriculture since the late 1980s (Supplementary Figure S1), which implies that the study area was exposed to agrochemicals when water samples were collected as it is today.
The study area has four distinct seasons. The summer in June to August is hot and humid, whereas the winter in December to February is cold and dry. In summer, typhoons bring strong winds and heavy rain. Specifically, in 1999, the highest monthly average temperature was 23.6°C in July and August, while the lowest was −2.7°C in January. Similarly, the highest and lowest in 2000 were 24.4°C in July and −3.3°C in February, respectively. The total annual precipitation was 1,360.6 mm in 1999 and 1,060.8 mm in 2000, of which 70 and 78% occurred in June to September in 1999 and 2000, respectively. The precipitation during the study period is shown in Figure 2. The total  Figure S2. Springs were divided into Type Ⅰ and Ⅱ when sampled based on the manner by which the discharge occurred (i.e., low and high, respectively).
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712865 amount of rainfall was 115 mm over approximately 24 h and the maximum rainfall intensity was <10 mm/h during a typhoon called Saomai in 2000 ( Figure 2B).

Geology
The geology of the study area consists of the Cambro-Ordovician Joseon Supergroup (mainly the Maggol (Omg) and Yeongheung (Oy) Formations), Carboniferous to Permian Pyeongan Supergroup (e.g., Hongjeom series (Ch)), and Late Triassic to Early Jurassic Daedong Supergroup (Bansong Formation (Jbs and Jbc)) as in Figure 1 (Han et al., 2006). The age of Kosong limestone (Ogl) on the east (the white area in Figure 1A) is unknown, and had been suggested to be Ordovician or later (Seo, 1997). The Maggol Formation (Omg) comprises a thick (250-400 m) sequence of carbonate rocks with a variety of lithology such as calcareous mudstone, dolostone, lime conglomerate, lime breccia, bioclastic grainstone, and oolitic grainstone, while the Yeongheung Formation (Oy; ca. 400 m thick) consists of massive to thick-bedded, light to dark gray dolostone in the lower part and bluish gray limestone in the upper part (Chough et al., 2000). These carbonate rocks occupy most of the study area ( Figure 1A), and are characterized by fractures ( Figure 1B) and sinkholes, caverns, and dolines (Supplementary Figure S2). The carbonate rock aquifers yield a usable quantity of water in the study area (http://www.gims.go.kr). Park et al. (2011) reported the hydraulic conductivity of 0.004-1.1 m/day (mean of 0.078 m/ day) in carbonate aquifers near the study area. Recently a tracer test was performed using uranine at the west of the study area (out of the map in Figure 1A) in January 2021 by a private company which plans to build a waste landfill in a historical limestone mine site. The uranine leaked to a small stream nearby 3 days after being released at the abandoned mine site (https:// news.kbs.co.kr/news/view.do?ncd 5158987), implying hydraulic connectivity.
Meanwhile, the Bansong Formation (Jbs and Jbc) consists of a conglomerate succession at the base and metamorphosed shale and sandstone with thin coal seams at the upper part, and unconformably overlies the carbonate sequence of the Joseon Supergroup or the siliciclastic sequence of the Pyeongan Supergroup (Han et al., 2006). The Bansong Group is highly metamorphosed and deformed in places such as the northeasttrending regional folds and thrusts developed during the Daebo tectonic event (Han et al., 2006;Ree et al., 2009), and represents syntectonic sedimentation during the late Early Jurassic in a compressional regime. Specifically, the sequence is structurally overlain by the Maggol Formation (Omg) along the Gongsuweon thrust and by the Yeongheung Formation (Oy) along the Deokpori thrust ( Figure 1B). According to Ree et al. (2009), the deposition of the Bansong Group and the thrusting were coeval during the middle Mesozoic Daebo tectonic event.
Given that springs in general originate along geologic contacts, faults, and ground depression (Fiorillo et al., 2018) and the layering of rocks relative to the slope of mountains affect both the number and discharge of springs in mountains (Mocior et al., 2015), the overturned strata (e.g., Cambro-Ordovician Yeongheung Formation (Oy) overlying Jurassic Bansong Formation (Jbs) along the thrust in Figure 1B) in this mountainous area can provide a pathway for deep groundwater upwelling as in Lorette et al. (2018) who showed deep water rising at springs through the faulted anticline structure of Périgueux, France. Similarly, Zhu et al. (2020) observed ascending springs along secondary tectonic fissures in the contact zone of limestone and magmatic rocks invaded as well as contact springs where the permeability of the local contact zone was enhanced. According to Zhu et al. (2020), the two contact springs with relatively large and stable flow rates reflected the high proportion of karst water replenishment from the deep circulation.

Sampling
Water samples were collected at two steps. First, various kinds of water were collected for major chemical compositions in the dry season (May/June; n 44) and wet season (July; n 27) in 1999 (Figure 1), including rainwater (n 2 from RW1 in Figure 1A by installing a rainfall collector), spring water (n 41 from 30 sites), surface water (n 8 from two sites at Dong River and four sites at streams), shallow groundwater (n 11 from seven preexisting wells), and deep groundwater (n 9 from six preexisting wells) to assess the hydrochemistry of each water group. Sampling dates were summarized in Supplementary Table S1 where sampling sites visited in both dry and wet seasons were in bold. In addition, water samples were collected for δD and δ 18 O of water (n 18) in the dry season and for δ 13 C of DIC (δ 13 C DIC ; n 10) and tritium (n 11) in the wet season to identify water sources. The lithology of each sampling site was assessed based on a geologic map (Figure 1), although the surface geology might not represent the subsurface geology through which the water passed, in particular in this stratigraphically complex study area with thrusts.
Spring water samples were divided into Type Ⅰ (n 11; much lower than the rates in Figure 2A) and Type Ⅱ (n 30; discharge rates similar to those in Figure 2A (> 480 L/min)) based on the flow type observed during sampling to assess the hydrochemical and isotopic differences of spring water depending on flow patterns and to determine the factors contributing to the differences in the study area with stratigraphic overturning along thrusts. Type Ⅱ springs were expected to exhibit large variations in discharge and chemical composition through time (Moore et al., 2009). In fact, Type Ⅱ quickly responded to rain events as in Figure 2B, and the discharge rate was over 1.7 m 3 /s at S46 (Type Ⅱ) during the storm event in 2000. However, Reynolds numbers were not calculated to confirm the laminar or turbulent flow due to little information on conduit sizes, and discharge rates were measured only at three Type Ⅱ springs in Figure 2 due to the limited accessibility to the springs in the wild areas, which is the limitation of this study.
In the second stage, a field investigation was conducted from September 15 to 19 in 2000 to characterize the short-term variability of isotopic compositions over the course of the typhoon event named Saomai. Three Type Ⅱ springs were chosen (S41, S45, and S46 in Figure 1A) from the Yeongheung (Oy) Formation to assess water sources through conduits. They were selected because of relatively easy accessibility and discharge measurement. Spring water samples for stable isotopes (δD and δ 18 O) were collected two or three times a day from the three springs (a total of 10 samples for each spring), and discharge rates were obtained through a constant cross section using a flowmeter (General Oceanics). In addition, a rainwater sample was collected at a relatively high elevation point (RW1 in Figure 1A).
Water samples were collected in 60 ml Nalgene ® bottles for the analysis of major compositions (Na + , K + , Ca 2+ , Mg 2+ , SO 4 2− , NO 3 − , Cl − , F − , and SiO 2 ) and δD and δ 18 O, while 1 L of water was collected for tritium analysis. Water samples were filtered through 0.45 μm pore size cellulose nitrate membrane filters (Whatman ® ) to eliminate suspended materials and then stored in the polyethylene bottles. To the samples for cation analysis, concentrated nitric acid was added to keep pH <2. All water samples were stored at a constant low temperature of 4°C until analysis in the laboratory. In addition, DIC was obtained for δ 13 C DIC analysis by precipitating barium carbonate (BaCO 3 ) with adding sodium hydroxide (NaOH) powder and barium chloride (BaCl 2 ) liquid to a 1 L water sample.

Analysis
In-situ measurements, including temperature, pH, redox potential (Eh), DO, and electrical conductivity (EC), were determined using a pH/Eh meter (Model Orion 290A), DO meter (Model Orion 835), and conductivity meter (Model Orion 130) in the field, while alkalinity was determined by titration in the field and re-examined with the Gran method in the laboratory.
Major cations and SiO 2 and anions were analyzed by inductively coupled plasma atomic emission spectrometry (ICP-AES; Model Perkin Elmer Optima 3000XL) and ion chromatography (IC; Model Dionex 120), respectively, at the Center for Mineral Resource Research (CMR) in Korea University. The δ 18 O and δD of water and δ 13 C DIC were determined with an isotope ratio mass spectrometer (Finnigan MAT 252) at CMR. The δ 18 O and δD values were reported related to the Vienna Standard Mean Ocean Water (V-SMOW) with a precision of ± 0.1 and ± 1‰, respectively, while δ 13 C DIC was expressed related to the Pee Dee Belemnite (PDB) and its standard deviation was ± 0.05‰. The tritium ( 3 H) contents in water samples were measured at the Korea Atomic Energy Research Institute with a liquid scintillation counter (Model Packard 3255) after a ca. 600 g water sample was electrolytically condensed to be 20 g. The analytical error was within approximately 0.2 TU.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712865 5 Thermodynamic, Statistical, and Hydrographic Separation Analysis Partial pressure of CO 2 (P CO2 ) and the saturation index of calcite (SI calcite ) and dolomite (SI dolomite ) in water samples were calculated using the physicochemical data of water samples and the PHREEQC program (Parkhurst and Appelo, 2013), by assuming an equilibrium state. Then, the binary relations between SI calcite and SI dolomite (Langmuir, 1971), between SI calcite and P CO2 (Doctor et al., 2006;Lacelle et al., 2008;Frondini et al., 2019), and between pH and HCO 3 − (Langmuir, 1971;Lacelle et al., 2008) were drawn to assess the hydrochemical evolution with carbonate minerals. Groundwater flows at relatively high velocities in a karst area, and it is difficult to assume the equilibrium state. However, equilibrium thermodynamics provide quantitative information about geochemical conditions, including the open or closed system dissolution defined by Langmuir (1971). Once the equilibrium is interpreted, the hydrochemical and hydrodynamic evolution can be deduced.
Principal component analysis (PCA) was applied to the physicochemical parameters and SI calcite and SI dolomite of water samples (n 69), excluding rainwater to assess major hydrogeochemical processes. PC loading and scores are considered together to assess major hydrochemical processes in the studied area, while water samples were plotted in the plot of PC scores to differentiate a major hydrochemical process for each water group (Moore et al., 2009;Bicalho et al., 2012;Lorette et al., 2018;Kim et al., 2021). PCA was conducted using the latest MATLAB software (https://www.mathworks.com) after log-transformation except pH.
Hydrograph separation curves were made for the three Type Ⅱ springs using δ 18 O and δD to better understand the discharge at Type Ⅱ springs. A two-component mixing model was applied with rainwater (R) and pre-storm water (PS), and then the contribution of R in a spring water sample (M) was evaluated as: where Q is discharge (L/s) and δ is the isotopic composition. Changes in Q M and Q PS were plotted with time to see the different proportion of Q PS in the total discharge during a storm event. Note that soil water and epikarst water were excluded due to little information despite their contributions to discharge (Lee and Krothe, 2001;Aquilina et al., 2005;Aquilina et al., 2006). Instead, the water and solutes stored in the vadose zone were discussed with the concentrations of anthropogenic contaminants (e.g., NO 3 − ) given the contaminant storage occurring in the vadose zone including soil and epikarst (Ghasemizadeh et al., 2012) as well as water storage (Perrin et al., 2003;Aquilina et al., 2006;Doctor et al., 2006;Jacob et al., 2008).

Hydrochemical and Isotopic Characteristics
Physicochemical data for each water group are summarized in Table 1. The water samples from both dry and wet seasons were combined to characterize the hydrochemistry because of no temporal variation between two seasons on the Piper's diagram with NO 3 − on the Cl − axis (Figure 3), which showed that most of the water samples had the Ca-(Mg)-HCO 3 type. Given that Types Ⅰ and Ⅱ showed different hydrochemical and isotopic compositions (Table 1 and Figures 3-6), the division seemed acceptable, although it was not confirmed using Reynolds numbers.

Groundwater and Surface Water
The deep groundwater, whose average depth to water level and average water table were respectively 102.0 and 258.0 m asl, contained higher Mg 2+ and HCO 3 − than the shallow groundwater with an average depth to water level of 11.3 m and an average water table of 265.6 m asl. Whereas the shallow groundwater contained higher Cl − and NO 3 − than the deep groundwater (Table 1; Figure 4). The hydrochemical difference indicated substantial carbonate dissolution in the deep groundwater and the shallow groundwater vulnerable to contaminants from the ground surface because Mg 2+ increases as the water-rock interaction proceeds in dolomitic geologic settings (e.g., Lorette et al., 2018;De Filippi et al., 2021), while Cl − and NO 3 − can be indicators of agricultural contaminants (Schilling and Helmers, 2008;Kim et al., 2021). Consistently, the deep groundwater had lower tritium contents (3.3-6.7 TU) than the shallow groundwater (8.7-9.4 TU) ( Figure 5A), indicating relatively old ages (Palcsu et al., 2021), and high δ 13 C DIC up to −10.0‰ ( Figure 5B), implying a larger carbonate dissolution effect. According to Marfia et al. (2004), the DIC in groundwater evolves initially under open system conditions (δ 13 C DIC of −25 to −17‰), and later under closed system conditions as the δ 13 C DIC increases (δ 13 C DIC of -17 to -13‰) as a result of carbonate dissolution.
In addition, all deep groundwater samples were in equilibrium or supersaturated with calcite and dolomite, while some of the shallow groundwater was undersaturated with both ( Figure 6A). Most of the Ca 2+ seemed to originate from calcite dissolution in the shallow groundwater, while from dolomite dissolution in the deep groundwater ( Figure 6B), according to Lacelle et al. (2008) who assumed that all Mg 2+ in water is derived from dolomite dissolution and the relative contribution of Ca 2+ from calcite dissolution can be estimated by subtracting the Mg 2+ from Ca 2+ concentrations in water.
Furthermore, the shallow groundwater had higher Ca 2+ /Mg 2+ but lower Ca 2+ /(Na + + K + ) ( Figures 6C,D) than the deep groundwater, probably due to the geology, as seen in Figures  6E,F, since most of the deep groundwater samples were obtained from limestone and dolomitic limestone and thus had high Ca 2+ + Mg 2+ and Ca 2+ /(Na + + K + ), whereas the shallow groundwater was obtained from various geological sources including the Bansong Formation (Figure 1; Supplementary Table S1). As the water-silicate interaction proceeded in the Bansong Formation, Ca 2+ /(Na + + K + ) decreased despite high Ca 2+ + Mg 2+ as Y37, MW3, and MW7 in Figures 6D,F. Meanwhile, the surface water occurred geologically in the Yeongheung Formation (Oy in Figure 1) and had relatively high Mg 2+ concentrations ( Figure 4A), implying the hydraulic Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712865 6  Rainwater had a charge balance out of ± 15%, possibly because of a failure to analyze for all dissolved species and/or analytical errors, but was included for comparison.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712865 connectivity between surface water and groundwater as shown in the tracer test with uranine at the west side of the study area (out of the map in Figure 1A) in January 2021, and was saturated with calcite and dolomite ( Figure 6A) except for a sample (ST4) from the Bansong Formation (Jbc in Figure 1), probably due to CO 2 degassing given low P CO2 (average 10 −2.93 atm in Table 1). A surface water sample (ST2) had a slightly low δ 13 C DIC (−14.2‰) compared to the groundwater in Figure 5B, but still high compared to the biogenic CO 2 (−26 to −21‰ in Jiang et al., 2013;Kang et al., 2020), indicating the effect of other carbon sources (e.g., calcite).

Spring Water
Geographically, Type I springs mainly occurred in the Jurassic Bansong Formation occupying topographically high areas  Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712865 8 (220-610 m asl) and had an average water level of 333.5 m asl (n 11), while Type II mainly discharged in the Yeongheung Formation and occasionally occurred in the Bansong Formation (Figures 1, 7) and had an average water level of 301.2 m asl (n 30). Type Ⅱ was also hydrochemically distinct from Type Ⅰ by high EC, Mg 2+ , HCO 3 − , and to a lesser degree, Ca 2+ ( Figure 4; Table 1). In addition, most of the Type Ⅱ springs were close to saturation with both calcite and dolomite ( Figure 6A) and had high Ca 2+ + Mg 2+ and Ca 2+ /(Na + + K + ) ( Figure 6D), while Type Ⅰ was undersaturated with both carbonate minerals ( Figure 6A) and had most of the Ca 2+ from calcite dissolution ( Figure 6B) and low Ca 2+ + Mg 2+ and Ca 2+ /(Na + + K + ) ( Figure 6D).
The hydrochemical differences between Types Ⅰ and Ⅱ were conflicting with the expectation that conduit flow (Type Ⅱ in this study) is characterized by low contents of dissolved elements and is undersaturated with respect to calcite (e.g., De Rooij and Graham, 2017), which was probably due to the different geology ( Figures 6F, 7) and water sources (Figure 8) as in Lorette et al. (2018) and supported by the PCA result ( Figure 9A). Four principal components (PCs) were extracted, explaining 71% of the total variance (Table 2). PC1, which explained the highest fraction of explained variance, was highly correlated with EC, Ca 2+ , Mg 2+ , HCO 3 − , SI calcite , and SI dolomite ( Table 2), reflecting the dissolution of carbonate rocks. Most of the Type Ⅱ springs were plotted on the positive side of PC1 with deep groundwater (Figure 9A), implying the discharge of deep groundwater at Type Ⅱ springs in carbonate rocks (Figure 8), while Type Ⅰ was on the left side ( Figure 9A), implying relatively little water-carbonate interaction (Figure 8). In addition, Type Ⅱ had relatively higher δ 18 O and δD values than Type Ⅰ in the dry season ( Figure 5C) except for one sample (Y25-1 from the Kosong limestone in Figure 1A), and relatively higher temperature and lower DO (Table 1), which also supported the discharge of deep groundwater at Type Ⅱ (Figure 8).
Similarly, Lorette et al. (2018) found high Mg 2+ up to 12.8 mg/ L and equilibrium with calcite during low-water periods at a Vauclusian-type spring located on a major faulted anticlinal structure, and explained the unusual hydrochemistry using the mixture with a Jurassic confined aquifer since high Mg 2+ concentrations were only observed in Jurassic confined aquifers (> 20 mg/L as a result of the presence of magnesiumcalcite) as in the study area (Table 1; Figure 3). Moore et al. (2009) suggested the deep-water source (Ca-Mg-SO 4 -type water from Well 2 with Mg 2+ up to 49.6 mg/L) and local diffuse recharge (Ca-HCO 3 -type water from Well 4 with low Mg 2+ between 1.2 and 2.2 mg/L resulting from rain water equilibrating with the Ocala limestone) in a system dominated by conduit flow to explain the changes in chemical compositions (e.g., Na + , Mg 2+ Figure 11. Sample names are given for some Type Ⅱ springs. Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712865 to deep groundwater rising and to superficial infiltration water at the Lez spring in the Mediterranean basin, and showed a decrease of the deep compartment participation to the Lez spring outflow due to intense exploitation that decreases the total hydraulic head in the aquifer. Fiorillo et al. (2018) deducted upwelling groundwater (with high 222 Rn activity) flux feeding the springs by the increasing of the hydraulic head in depth. The high flow condition is sustained by high ascendant hydraulic gradient, under which higher water velocity and higher radon activity occur. Meanwhile, spring water contained high NO 3 − compared to groundwater (Figures 3, 4C). The contamination of spring water with NO 3 − was supported by PC3 which was positively related with Cl − and NO 3 − ( Table 2), indicating anthropogenic contamination given that Cl − and NO 3 − are major contaminants in agricultural areas (e.g., Schilling and Helmers, 2008;Kim et al., 2021). Most of the groundwater was plotted on the negative side of PC3, and a shallow groundwater well (Y37-1) occurring on the Bansong Formation with 0.2-0.8 mg/L for NO 3 − and 1.6-2.1 mg/L for Cl − had the lowest PC3 scores ( Figure 9B). In contrast, about half of the spring water samples were on the positive side of PC3. In particular, Type Ⅱ had high NO 3 − and Cl − concentrations ( Figure 4C), implying the fast infiltration of contaminants through conduits (Figure 8). The fast infiltration from surface was supported by the tritium concentrations ( Figure 5A) similar to the 3 H levels around 10 TU in precipitation during the study period in South Korea (Chae and Kim, 2019). Type Ⅰ also showed relatively high tritium concentrations but less than Type Ⅱ, which suggests the recent or short-circulating discharge of meteoric water occurring at both springs (Palcsu et al., 2021).  Figure 1A, Ogl Kosong limestone (the white area at the east in Figure 1A). The letters on the diagram indicate the names of distinct samples. In summary, the Type Ⅱ springs seemed to be fed by both deep groundwater and recent conduit flow (Figure 8) as in Lorette et al. (2018) who suggested that Toulon Springs are fed by two aquifers: (i) the Jurassic confined karst aquifer responsible for low-water period and related to dolomitic hydrochemical characteristics, and (ii) the upper Cretaceous unconfined aquifer responsible for quick hydrodynamic and hydrochemical variations. In contrast, Type Ⅰ seemed to discharge diffuse sources with low residence time ( Figure 8) as in the Tufa Springs system provided by Tobin and Schwartz (2012) who calculated specific conductance in unsampled diffuse sources in mountain hydrologic systems that ranged from 34 to 257 mS/cm. δ 13 C DIC was in a wide range (−15.3-−11.2‰ in Figure 5B), implying multiple carbon sources from calcite, atmospheric CO 2 , and soil CO 2 given the δ 13 C of carbonate rocks (∼0.4‰), atmospheric CO 2 (−9-−7‰), and soil CO 2 (−26-−21‰) (Jiang et al., 2013;Kang et al., 2020), and land use (Supplementary Figures S1, S2).

Major Hydrochemical Processes
Based on the PCA results (Table 2 and Figure 9) and hydrochemical and isotopic characteristics, water-rock interactions with carbonate (PC1) and silicate rocks (PC2) and anthropogenic contamination (PC3) can be suggested as major hydrochemical processes in the study area. PC4 was positive with pH but negative with SiO 2 , and surface water and some deep groundwater samples had high PC4 scores, whereas PC4 could not be clearly explained given the information and thus was excluded from further discussion.

Effects of Geology
Type Ⅱ springs and deep groundwater had high EC, Mg 2+ , HCO 3 − , SI calcite , and SI dolomite , which were positively correlated with PC1 (Table 2), and also had high Ca 2+ + Mg 2+ and Ca 2+ / (Na + + K + ) ( Figure 6D). These properties can be explained by the carbonate dissolution given that most of Type Ⅱ springs and deep groundwater occurred in the limestone and dolomitic limestone (Figures 1, 6F, 7), whereas a surface water sample (ST4) and Type Ⅰ spring sample (S68-3) with negative PC1 scores occurred in the Bansong Formation ( Figure 1) and were undersaturated with calcite and dolomite ( Figure 6A).
However, the deep groundwater had Ca 2+ from dolomite dissolution ( Figure 6B) and low Ca 2+ /Mg 2+ close to one ( Figure 6C), whereas the Type Ⅱ springs contained higher Ca 2+ than the deep groundwater (Table 1; Figure 4B), probably because of the extensive water-rock interaction in deep groundwater, causing calcite precipitation ( Figure 6A) and thus decreasing Ca 2+ /Mg 2+ . According to Langmuir (1971), the molar Ca 2+ /Mg 2+ ratio as low as 0.6 results from incongruent dissolution of dolomite. Besides, calcite is precipitated at low temperature when the molar Ca 2+ /Mg 2+   ratio is greater than one given a large uncertainty (Drever, 1997).
The deep groundwater showed the molar Ca 2+ /Mg 2+ ratios between 0.8 and 1.4 (average of 1.1; Figure 6C), which indicates the possible condition for the incongruent dissolution of dolomite and calcite precipitation. The lower average P CO2 of deep groundwater (10 −2.46 atm) than that of Type Ⅱ springs (10 −2.28 atm) also supported the prolonged dissolution of calcite in deep groundwater in a closed system ( Figure 10) since CO 2 is consumed for calcite dissolution (Langmuir, 1971;Lacelle et al., 2008;Frondini et al., 2019). Meanwhile, the inconsistent Mg 2+ contents in carbonate rocks might affect the Ca 2+ /Mg 2+ ratios. For instance, the water sample from the Kosong limestone had high Ca 2+ /Mg 2+ regardless of water groups ( Figure 6E), probably because the Kosong limestone consists mainly of calcite. PC2 was positively correlated with Na + and K + but negatively with Eh ( Table 2). Some groundwater samples had high positive PC2 scores, while most springs were plotted on the negative side of PC2 ( Figure 9A), indicating that PC2 addressed silicate dissolution. In particular, a deep groundwater sample (MW7) from the Bansong Formation had the highest PC2 score. This deep groundwater sample had the highest Na + and F − (Figure 7). In addition, some groundwater had relatively low Ca 2+ /(Na + + K + ) values due to high Na + and K + in the Bansong Formation (e.g., Y37, MW3, and MW7; Figures 6F, 7), indicating silicate dissolution and thus geological heterogeneity in this carbonate aquifer. Especially, two deep groundwater samples (MW7, MW3) imply extensive water-silicate interaction but little anthropogenic contamination in the deep aquifer based on low Ca 2+ /(Na + + K + ) ( Figure 6D), high Na + and F − concentrations (Figure 7), and low NO 3 − and Cl − concentrations ( Figure 4C). Note that Type Ⅰ springs also had low Ca 2+ /(Na + + K + ) values and occurred in the Bansong Formation ( Figures 6D,F). However, the effect of silicate dissolution seemed slight in Type Ⅰ given the similar concentrations of Na + , K + , and F − in both types of springs (Table 1) despite the different geology (Figures 1, 7). The low contents of Na + and K + as well as EC and Ca 2+ would be probably because of the discharge of diffuse sources with low residence time at Type Ⅰ springs (Figure 8) as in the mountain systems provided by Tobin and Schwartz (2012), which also can explain Type Ⅰ containing lower concentrations than the shallow groundwater for all major ions except NO 3 (Table 1). Similarly, Moore et al. (2009) showed low and relatively constant concentrations (e.g., Mg 2+ and SO 4 2− ) of diffuse recharge at a well in a karstic aquifer due to low inputs of deep water.

Anthropogenic Contamination
Spring water contained higher NO 3 concentrations ( Figure 4C) and PC3 scores ( Figure 9B) than groundwater, indicating that the contaminants stored in the vadose zone discharged through springs. In particular, the point infiltration through conduits seemed to deliver more contaminants, causing higher Cl − and NO 3 − concentrations at Type Ⅱ ( Figure 4C) given that contaminant storage occurs in the rock matrix and epikarst, while contaminants are mostly transported along preferential pathways in karst systems (Ghasemizadeh et al., 2012). It should be noted that the Type Ⅰ springs occasionally (e.g., S37) had higher ratios of Cl − + NO 3 − than Type Ⅱ on the Piper's diagram in Figure 3, which was because the Type Ⅰ springs contained relatively low HCO 3 − concentrations similar to rainwater ( Figure 4A) implying the discharge of diffuse sources with little water-rock interaction ( Figure 8).
Shallow groundwater was also affected by the influx of anthropogenic contaminants from the surface, given the high NO 3 − or Cl − concentrations ( Figure 4C). In particular, the Y28 had the highest PC3 score ( Figure 9B) due to a high level of both NO 3 − and Cl − ( Figure 4C) and showed a high level of Mg 2+ , Na + , and SO 4 2− (Figure 7), which indicates the combined effect of carbonate and silicate dissolution and anthropogenic contamination since the subsurface can be the Bansong Formation at the geologic boundary, although the surface geology of Y28 was the Yeongheung Formation ( Figure 1). Also, the sampling site (Y28) was close to the Dong River and agricultural land (Supplementary Figure S2). In addition, Y37 had high K + and Cl − (Figure 7) and thus low Ca 2+ /(Na + + K + ) ( Figure 6D). MW2 also had high K + (Figure 7). High K + , Na + , and F − were occasionally found in groundwater (Figure 7) and the average Na + and K + concentrations of groundwater were larger than the corresponding concentration of spring water (Table 1), implying the effect of silicate dissolution based on their (e.g., MW7, Y37) occurrence in the Bansong Formation. Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712865 However, the effect of agrochemicals cannot be ignored for high K + as in Stueber and Criss (2005) and Moore et al. (2009) Figure S2) and the fact that K + increased with Cl − during the wet season (July) in shallow groundwater (e.g., Y37 and MW2 in Figure 7 and Supplementary Figure S3). Furthermore, HCO 3 − concentrations were higher than those predicted from calcite dissolution at a given P CO2 (Supplementary Figure S4a), which indicates additional proton (H + ) sources. There was no evidence of pyrite to produce H + in the study area, while nitrate contamination might enhance carbonate dissolution . In fact, the water samples were plotted between the enhanced dissolution of carbonate rocks by nitrification (i.e., the 1:1 of Ca 2+ + Mg 2+ (mmol) and HCO 3 − (mmol)) and natural dissolution (i.e., the 1:2 of Ca 2+ + Mg 2+ (mmol) and HCO 3 − (mmol)) (Supplementary Figure S4b). PC3 was negatively related with pH but positively with NO 3 − ( Table 2), implying the acidification and nitrate contamination of water.

Subsurface Flow in the Study Area
Based on the hydrochemical difference of spring water depending on the flow type and PC1 (carbonate dissolution) in Figure 9A and their spatial distribution on the geologic map in Figure 1A, it can be deduced that the heterogeneous geology causes different flow patterns in the study area. The carbonate rocks seemed to provide the conduits for Type Ⅱ due to their relatively high solubility and fractures ( Figure 1B), whereas the sandstone and shale provided a matrix for the diffusive flow of Type Ⅰ (Figure 8). The different hydrological conditions may further differentiate the hydrochemistry, by causing different contributions of endmembers (e.g., groundwater, rainfall) in spring discharge, as in Bicalho et al. (2012) who found that hydrological conditions induce the different proportion of end-members (deep groundwater rising and superficial infiltration water) participating in groundwater mixing of a spring.

Type Ⅰ Spring
Type Ⅰ springs seemed to slowly discharge the rainwater infiltrated from soil cover in the Bansong Formation, carrying contaminants (e.g., NO 3 and Cl) and dissolving calcite via shallow (near-surface) circulating pathways (Figure 8) based on the following results: low concentrations of major compositions including HCO 3 − and EC ( Figure 4; Table 1), undersaturation with respect to calcite ( Figure 6A) and Ca 2+ from calcite dissolution ( Figure 6B), and little water-rock interaction (i.e., negative PC1 and PC2 scores in Figure 9A). Low P CO2 similar to that in rainwater despite slow flow also suggests a short residence time or bare soil given soil P CO2 around 10 −1.5 -10 −2.5 atm (Appelo and Postma, 2005) or higher due to an accumulation of CO 2 in the epikarst (Peyraube et al., 2012). Type Ⅰ seemed to represent springs in mountains following a relatively shallow flow path with little contribution of groundwater because of steep slopes and shallow soil development (Somers and McKenzie, 2020).
Most of Type Ⅰ springs were undersaturated with respect to calcite and showed P CO2 near typical biologically respired values (10 −2.5 atm in Lacelle et al., 2008; Figure 10A). In addition, the plotting of the Type Ⅰ springs on pH against HCO 3 − suggests an evolutionary path by the dissolution of carbonates under open system conditions with P CO2 near 10 −2.5 atm ( Figure 10B), which increases both the pH and HCO 3 − of spring water with subsurface residence times. The significant amounts of NO 3 − and high DO also suggested an open system (Table 1).

Type Ⅱ Spring
The relatively high levels of NO 3 − , Cl − ( Figure 4C), tritium, and low δ 13 C DIC down to -15.3‰ in Type Ⅱ springs ( Figure 5) supported the fast flow of recent or short-circulating recharge through conduits and the admixture of biologically sourced carbon (Figure 8). In addition, Type Ⅱ had high Mg 2+ , HCO 3 − ( Figure 4A), δ 18 O, δD ( Figure 5C), SI calcite , and SI dolomite ( Figure 6A) similar to the deep groundwater, which can be explained by the discharge of deep groundwater at Type Ⅱ springs (Figure 8). Other possible causes such as CO 2 exsolution, increases in temperature, or dissolution of soluble minerals were excluded given the higher P CO2 (average 10 −2.28 atm) and lower temperature (14.7°C) in Type Ⅱ springs than those in groundwater ( Table 1) and the fast flow rates (Figure 2). A negative relation between P CO2 and SI calcite ( Figure 10A) and the relationship between pH and HCO 3 − along the saturation curve ( Figure 10B) also support the hydrochemical changes of Type Ⅱ springs due to changes in water chemistry that affect P CO2 (e.g., carbonate dissolution or precipitation), similar to the deep groundwater.
Type Ⅱ seemed to discharge both the flow component through the epikarst zone and the baseflow from the lower reservoir, as defined by Tritz et al. (2011) who used two reservoirs (i.e., the upper epikarst/soil zone and the lower vadose and saturated zone) to simulate the behavior of a karst system catchment and assumed a fast flow through the epikarst and slow baseflow from the lower reservoir. The deep groundwater rising at springs was also observed by other researchers including Moore et al. (2009), Bicalho et al. (2012, Demiroglu (2016), Gil-Márquez et al. (2019), Lorette et al. (2018), andZhu et al. (2020). Some implied fast transfer through conduits unlike Tritz et al. (2011). For instance, Lorette et al. (2018) showed Toulon Springs with a mean annual daily discharge of 450 L/s fed by deep confined karst aquifers during low-water periods. Gil-Márquez et al. (2019) mentioned the piston-flow effect at the beginning of floods, which causes the drainage of deep ascending flows through the saturated zone and subsequently the rise in mineralization at a spring. Baudement et al. (2017) mentioned chimney-shafts, by which deep phreatic karst systems are connected to vauclusian springs. Based on these previous study results and the results of this study, it can be inferred that karst conduit networks including palaeo-karstic systems or abandoned cavities seemed to still be communicating with the main drainage axis in the study area as in Aquilina et al. (2005) and Pogačnik et al. (2017).

Proportion of Groundwater at Type Ⅱ Springs
The proportion of baseflow at Type Ⅱ springs was evaluated using the hydrograph separation. δ 18 O and δD decreased as it began to rain heavily at three Type Ⅱ springs and increased back ( Figure 11A), reflecting the changing of water sources (e.g., rainwater (R) and pre-storm water (PS)) in the spring discharge during the typhoon event. The PS was represented using the average δD and δ 18 O of deep groundwater collected during the dry season (δ 18 O −9.6 and δD −70.2‰; n 3; Figure 5C), given the hydrochemistry of Type Ⅱ similar to that of deep groundwater. It was assumed that the deep groundwater in the dry season records the isotopic compositions of baseflow at Type Ⅱ. In fact, the spring water collected during the same period (n 11; Figure 5C) had a similar isotopic composition to the groundwater on average, indicating baseflow discharge in the dry season. All the δD and δ 18 O were plotted slightly below the meteoric water line ( Figure 5C), indicating evaporation or water mass mixing (Serno et al., 2017). Also, all three springs had initial values similar to the average of deep groundwater in Figures 5D, 11A, reflecting baseflow discharge before the storm.
In contrast, a rainwater sample (R) collected during the storm event had −12.0‰ for δ 18 O and −85.0‰ for δD ( Figure 5D). These values were quite different from those of PS. This large isotopic difference in the two end-members (i.e., 2.4‰ for δ 18 O and 14.3‰ for δD) and the plotting of spring water between R and PS supported the applicability of the two-component hydrograph separation in the study area. Given the low δ 18 O and δD of R ( Figure 5D), the lower δ 18 O and δD value of spring water indicates a higher proportion of R in discharge; thus, high proportions of R were expected in discharge at S46 given the great decreases in δ 18 O and δD ( Figure 11A).
We chose the bigger proportion of R between those using δ 18 O and δD (Table 3; Figure 11B). As a result, the proportion of R was similarly 17, 13, and 13% of the total discharge at S41, S45, and S46, respectively, before the storm. At peak flow, however, R accounted for 53, 39, and 87% at S41, S45, and S46, respectively, while the maximum contribution of R was 60, 53, and 87%.
This result indicates that the proportion of groundwater was up to 83-87% in the three Type Ⅱ springs in the dry season, similar to Moore et al. (2009) who identified that upwelling from deep flow paths provides significant contributions of water to spring discharge and emphasized that restricted monitoring of springs limits the interpretation of karst systems by masking critical components of the aquifer. The proportion of baseflow decreased during the storm, while the contribution of rainwater increased at Type Ⅱ (Table 3; Figure 11) as in Moore et al. (2009) who showed that the contribution from the deep source depends inversely on flow conditions. Note the high tritium concentrations at Type Ⅱ ( Figure 5A) despite the high proportion of deep groundwater in the dry season was probably because the water samples for tritium were obtained in wet season including S46.

Storage Effect
S46 had the maximum contribution of rainfall (87%) at peak flow without a delay, while S41 and S45 had a delay between the peak flow and maximum contribution of R (Table 3). At the end of monitoring, the proportion of R decreased to 21 and 17% at S41 and S45, respectively, indicating the recover to the baseflow condition, whereas it was still 60% at S46 unlike before the storm (i.e., 13%).
In addition, S46 did not show an increase in NO 3 − during the typhoon in 2000. Rather, the NO 3 − was diluted as other major compositions during the wet season in 1999 ( Figure 12; Supplementary Figure S5). In contrast, NO 3 − increased at S41 and S45 during the typhoon event in 2000 as well as wet season in 1999 compared to during the dry season in 1999, indicating the flushing out of NO 3 stored in the vadose zone, as observed by other researchers including Doctor et al. (2006) and Lorette et al. (2018). The influx of NO 3 − from the vadose zone seemed to be more influential than the dilution effect of rainwater at S41 and S45, as in Lorette et al. (2018) who observed that during recharge events the first flood event leads to an increase of NO 3 − concentrations, resulting in the mobilization of NO 3 − storage in soil and the unsaturated zone of the aquifer.
Given the similar NO 3 levels ( Figure 12) as well as geographical locations (Figure 1), geology (Yeongheung Formation), and flow type (Type Ⅱ), the distinguished behavior of NO 3 − ( Figure 12) and water isotopes ( Figure 11) at S46 was probably due to the low storage capacity of the vadose zone for water and contaminants. Water and solutes stored in the vadose zones can be discharged through springs following rainfall events in karst areas (Perrin et al., 2003;Aquilina et al., 2005), which may cause buffered rain isotope responses despite highly peaking hydraulic response (Perrin et al., 2003), increase P CO2 and anthropogenic components (e.g., Cl − ) and decrease δ 13 C DIC by the contribution of water held in storage in the vadose zone to the conduit system in storm events (Doctor et al., 2006), or increase NO 3 − at the first flood event during recharge events, with the decrease in NO 3 − in next flood events due to the stock of NO 3 − storage in the soil and the unsaturated zone moved out of the karst system during the previous flood event (Lorette et al., 2018). Jacob et al. (2008) showed spatial and temporal variations of vadose zone water storage in a karst system using groundbased absolute gravimetry measurements.
The low storage seemed to cause the discharge quickly responding to rainfall not only during the typhoon event in 2000, but also in the small precipitation events in 1999 (Figure 2), and the rainfall accounting for 87% of the discharge at peak flow and continuing to flow out and contributing more than half of the discharge after the heavy rainfall stopped at S46 (Figure 11; Table 3). Similarly, Demiroglu (2016) reported that developed karst sinkholes allow fast percolation of heavy rainfall into an aquifer up to 80%, and explained that the very low storage combined with the high transmissivity means that most of the recharge will not be retained by the karst system, but will rapidly flow out to springs, rivers, lakes, or seas. In contrast, the delayed response to the storm and quick recovery to the baseflow condition after it stopped raining at S41 and S45 ( Figure 11; Table 3) were probably due to the storage effect of the vadose zone, given that all three springs had flow characteristics through permeable conduits. Doctor et al. (2006) considered the anthropogenic component stored within the epikarst as one of the recharge sources, and found that an anthropogenic component derived from epikarstic storage affects the well under conditions of elevated hydraulic head, accounting for the chemical response in the well during wet conditions. A similar increase in NO 3 − and Cl − during the wet season was observed at the wells L6-1 (deep) and Y37-1 (shallow) as well as the increases in K + and Cl − at MW2 and Y37 (shallow) in the study area (Supplementary Figure S3).
We acknowledged that the storage effect of the subsurface in the study area needs to be further investigated since the storage degree was evaluated only with the response process at this current study, although the subsurface storage is involved with many processes, including response, recession, and discharge, which are also affected by aquifer structures such as the development of conduits, degree of weathering, and recharge 3 | Contribution (%) of rainfall (R) to discharge at three Type Ⅱ springs during the storm event in Figure 2B. The remaining was assumed to come from pre-storm water (PS; deep groundwater). The bigger proportion of R between those using δ 18 O and δD was in bold. Time for peak flow was shaded in gray.

S41
S45 S46 areas. In addition, there is an opposite interpretation from Demiroglu (2016) who mentioned that the variability of physicochemical characteristics is higher in conduit permeability with low storage. Stueber and Criss (2005) suggested that less soluble ions increase with flow as they are mobilized from fields to karst conduits under storm conditions, while highly soluble ions supplied by diffuse groundwater are diluted by high flows.

CONCLUSION
Spring water quality was mainly controlled by three hydrogeological processes in the studied karst area: 1) carbonate dissolution (PC1 in Figure 9) anthropogenic contamination (PC3 in Figure 9) rainfall infiltration ( Figure 11). Precisely, the Type Ⅱ springs, which was assumed to be conduit flow-dominated, were greatly affected by carbonate dissolution and had high EC, Ca 2+ , Mg 2+ , HCO 3 − , Ca 2+ + Mg 2+ , Ca 2+ /(Na + + K + ), and SI calcite , similar to the deep groundwater, probably because of their spatial occurrence in the carbonate rocks and Type Ⅱ springs supplied by the deep groundwater ( Figure 8). Hydrograph separation using δ 18 O and δD data showed that the deep groundwater contributed to 83-87% of the discharge in three Type Ⅱ springs in the dry season. The contribution of baseflow decreased as the contribution of rainfall increased during storm events at Type Ⅱ springs. A Type II spring (S46) rapidly responded to rainfall and showed the contribution of rainwater exceeding 50% 2.5 days after the rain stopped, with little NO 3 released from the vadose zone during the storm, implying low storage. In contrast, the other two Type II springs (S41 and S45) had delayed responses to precipitation and had NO 3 flushed out from the vadose zone during the wet season, indicating that the vadose zone may store rainwater and release contaminants stored depending on the type of precipitation.
In contrast, the effect of silicate dissolution (PC2 in Figure 9) appeared in groundwater but not in spring water, probably due to the short circulation of Type I springs despite their occurrence in the metamorphosed sandstone and shale. Type Ⅰ springs seemed to slowly discharge infiltrated rainwater, dissolving carbonate minerals in the relatively high-altitude area. Both types of springs were contaminated with NO 3 − and Cl − compared to the groundwater, and in particular, Type Ⅱ had high average concentrations of Cl − and NO 3 − , indicating severe contamination through conduits.
Based on the study results, the heterogeneous geology affected both the hydrochemistry and flow types of spring water in this stratigraphically complex karst area with thrusts, causing unusual hydrochemical characteristics: elevated concentrations of Mg 2+ at Type Ⅱ and low concentrations at Type Ⅰ. The carbonate rocks provided the conduits for Type Ⅱ, whereas the silicate rocks provided a matrix for Type Ⅰ, which further caused hydrochemical differences depending on the flow type.
The study results suggest that hydrochemical and isotopic data are applicable to understand the flow paths of springs (e.g., conduit flow through carbonate rocks) and the storage characteristics in the vadose zone in a karst area since the hydrological condition (e.g., upwelling of deep groundwater, rainfall recharge through a matrix) affects the participation of water sources and consequently the hydrochemical and isotopic composition of springs. Mg 2+ , Ca 2+ /Mg 2+ , Ca 2+ /(Na + + K + ), and SI calcite were good indicators to distinguish the geology that the water passed through (e.g., limestone, dolomitic limestone, sandstone, and shale) and the hydrochemical evolution of water samples (e.g., incongruent dissolution of dolomite) in the tectonically complex karst area. PCA was useful to characterize the hydrochemical processes in each type of spring. Also, the present study showed the successful applicability of water isotopic tracers for evaluating the FIGURE 12 | The temporal variation of major compositions at three Type Ⅱ springs in the precipitation events in Figure 2. The data in 2000 were the average of 10 samples collected for 5 days (Kim et al., 2001).
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 712865 hydrodynamic properties of subsurface groundwater flow, although this study failed to quantify the effect of epikarst or soil water compared to rainfall and groundwater during storm events due to little information, which remains future work. Besides, upwelling of deep groundwater at Type Ⅱ springs estimated by the hydrochemical and isotopic analysis should be confirmed by a hydrogeological investigation. In particular, the main drainage network for deep groundwater upwelling needs to be identified in the near future given a new plan to build a waste landfill in a historical limestone mine site around the study area. Thrust faults should be first assessed as a potential pathway. In addition, the classification of springs is needed to be confirmed based on Reynolds numbers.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
SY interpreted the data and wrote the manuscript, GC collected the data and drew figures, JO conducted the principal component analysis and drew figures, S-HK and D-IK collected and interpreted the data, and S-TY initiated the project and wrote the manuscript.