Use of major ion chemistry and trace and rare earth elements to characterize hydraulic relations, mixing processes and water–rock interaction in the Dong Van karst aquifer system, Northern Vietnam

The effectiveness of using a groundwater geochemistry approach in karst hydrogeologic research is highlighted. In particular, this approach is useful for preliminary investigations, such as for the study described here on the Dong Van karst aquifer system in Northern Vietnam. Analyses of different groundwater chemistry parameters complement each other, to clarify hydrochemical processes that are occurring in the karst system. The results of this study show that major ion composition can be used to clarify water chemistry signatures, as well as to identify the mixing processes and water–rock interactions in aquifers. Meanwhile, trace element concentrations and rare earth element patterns can be used as potential natural tracers when some processes are not revealed through conventional hydrochemical methods. These natural tracers can also be used to identify contaminant sources and/or contaminant transport pathways in karst aquifers. Viewed holistically, the groundwater geochemistry approach provides scientific information to establish a basic hydrogeological conceptual model and to estimate the water balance, which has implications for water resources protection and management in karstic systems.


Introduction
Groundwater from karst aquifers is the major freshwater source for drinking-water supply and agricultural irrigation in many countries and regions throughout the world.The number of karst water consumers is estimated at 678 million, or 9.2% of the world's population (Stevanovic 2019), while the population living in karst areas is estimated at 1.3 billion people (Goldscheider et al. 2020).As pointed out by Do 2001, in Vietnam, limestone accounts for nearly 18% (about 60,000 km 2 ) of the total area, mostly in the northern provinces, while an updated result by Goldscheider et al. 2020 shows that carbonate rocks cover an area of 90,000 km 2 (occupying 27.2% of the total area).The Dong Van karst plateau (DVKP), located in Ha Giang province, is one of the most prominent karst regions in northern Vietnam, and belongs to the extended part of the South China karst belt (Yunnan karst plateau; Fig. 1a).The Na De spring system (the source of the Nhiem River) is the largest spring in the area, which could be the local base level for karst erosion in the DVKP (Tran et al. 2022) and is probably among the largest karst springs in Southeast Asia (Fig. 1c).There is no regular discharge monitoring at this spring.During a one-time measurement during the rainy season, a discharge of 44 m 3 /s was observed (Ender, Institute of Applied Geosciences (AGW), Karlsruhe Institute of Technology (KIT), unpublished data from a field trip report for Dong Van Karst Plateau, 2014), but much higher discharge values, probably several hundred m 3 /s, are expected to occur during extreme high-flow events.There is no specific research on the karst hydrogeological system associated with this major spring.Therefore, research on the hydrogeological characteristics of the Na De spring system and its surrounding area is critical for sustainable groundwater resources management, protection, and exploitation on the Dong Van karst plateau.Groundwater geochemistry has emerged as an effective tool that can elucidate water sources, hydrochemical processes, and groundwater dynamics, and to identify contamination input sources in karst aquifers (Lambrakis et al. 2004;Holland and Witthüser 2009;Goldscheider and Drew 2014).The use of major ions to identify the processes controlling groundwater chemistry has become a common method aided by conventional geochemical displays such as Gibbs diagrams and ion-ratio graphs (Gibbs 1970;White 1988;Goldscheider 2015 andGao et al. 2019).Previous studies by Güler et al. 2002 andHeikinien et al. 2002 applied multivariate statistical methods-e.g.principal component analysis (PCA) and hierarchical cluster analysis (HCA)-to compare and partition water samples into homogeneous groups based on water physicochemical parameters, providing information on the physicochemical processes controlling water chemistry, as well as providing insight into underground flow connectivity.In recent years, many studies have shown that trace elements can be used as potential natural tracers of groundwater origin when some processes are not revealed through conventional hydrochemical methods such as signatures of major ions and/or stable isotopes.As indicated already (Lee et al. 2003;Möller et al. 2004;Han and Liu 2007;Zhou et al. 2012), the rare-earth elements (REE) abundance and patterns can be exploited effectively as fingerprints of host aquifers or bedrock as well as indicating paleo-hydrological conditions.Gill et al. 2018 andCholet et al. 2019 used temporal variability of REE concentrations and patterns to evaluate possible links between different water sources and to quantify mixing contributions.In the same aspect, a study by Benischke 2021 provided an up-to-date and comprehensive overview of tracers for karst system analysis, whereas Sivelle et al. 2022 pointed out the use of excess air as representing trace gases for karst connectivity studies.In addition, Schilling et al. 2023 demonstrated the enormous potential of combining vanadium, eDNA, and helium as tracers of fault-facilitated fast groundwater flow paths.
The main aim of this work is to evaluate the effectiveness of the groundwater geochemistry approach for karst hydrogeological research, with a particular focus on using trace and rare earth elements as natural tracers.The work uses these methods to contribute to a better understanding of the Dong Van karst aquifer system.The specific objectives of this study are to: (1) use major ion chemistry to identify and characterize different water types, mixing processes, and water-rock interactions in the aquifers; (2) test the use of trace and rare earth elements as natural tracers for characterizing underground flow connectivity and identifying contaminant transport processes; and (3) establish a conceptual model of the karst aquifer system and produce a water balance of the study region to provide scientific information for water resources management.

Geological and hydrogeological settings
The study area is situated in the Dong Van Karst Plateau UNESCO Global Geopark, covering an area of about 620 km 2 , which is composed mainly of carbonate rocks (more than 60%).According to previous and recently updated geological maps, the strata range mainly from Paleozoic to Mesozoic, with local outcrops of Triassic gabbro and diabase along deep faults (Fig. 1b).Carbonate rock formations, including Bac Son (C-P bs), Dong Dang (P 2 đđ) and Hong Ngai (T 1 hn), are widely distributed in the central part of karst plateau.The Devonian formations including Mia Le (D 1 ml), Si Phai (D 1-2 sp) and Toc Tat (D 3 tt) occur in narrow bands, located along the north and northeast.In addition, Triassic formations, including Song Hien (T 1 sh) and Yen Binh (T 2 a yb), are distributed in the south and SW of this area.Tectonically, the investigated area is dissected by many fault systems (Fig. 1b), but the main system trends NW-SE.There are also sub-latitude and sub-meridian fault systems that trend in a NE-SW direction that were formed during a Neotectonic period.
The cave system developed a strong relationship to the tectonic and geomorphologic development processes in this area (Tran et al. 2004(Tran et al. , 2013(Tran et al. , 2022)).There are some stream caves such as Na Luong cave, Ma Le cave system, Tia Sang cave and Hang Rong cave, that are important water sources for domestic use and irrigation that have been explored by Belgian experts from University of Leuven's caving club (SPEKUL) and Vietnam Institute of Geosciences and Mineral Resources (VIGMR) speleologists.
The Bac Son Formation (Fm.) is composed of homogeneous, massive, partly dolomitized limestone, totaling 700-1,000 m in thickness, and is considered to be one of the most hydrologically important karst aquifers in this area, after the Dong Dang Fm. and Hong Ngai Fm. karst aquifers.The fractured Si Phai Fm. aquifer is also present (Ender et al. 2018;Nguyen et al. 2013).
According to published geologic data, Na De spring is located along the boundary between the limestone of Bac Son Fm. (C-Pbs) and terrigenous sediments of Song Hien Fm. (T 1 sh; Hoang 1976).A study based on cave survey data by Tran et al. 2022 indicates that the Na De spring could be the local base level for karst erosion in the Dong Van karst plateau, which probably drains large parts of this area.

Climatic features
Northern Vietnam's climate is generally of the subtropical monsoon type (Pham 1985); however, the Dong Van karst plateau belongs to the high-elevation region with strongly dissected terrain, and thus the climate is characterized as a temperate zone (Tran 2010) and is divided into two seasons: rainy season (May to October) and dry season (November to April).According to the data collected at the climate station Dong Van District in the period 2015-2020, operated by Vietnam National Centre for Hydro-Meteorological Forecasting (NCHMF), the annual average precipitation is 1,537 mm, with the highest rainfall in June (302 mm average) and lowest in February (12 mm).The annual average temperature is 18.4 °C, with the highest reaching 23.6 °C in June and the lowest 11.7 °C in December.

Sampling campaigns and sampling sites
Two sampling campaigns were organized between the dry (March) and rainy (June/July) seasons, 2020.A total of 40 water samples were collected from springs (Na De spring system and Southwest springs), in caves (Na Luong cave, Hang Rong cave), the recharge stream (surface flow) of Tia Sang cave (TSC), and hanging pools in the mining area.(Fig. 1b).
The Na De spring system (NDSS) is situated in Na De Valley in the south of DVKP.This system includes 16 springs, of which 4 springs appear continuously (BM1, BM2, BM6, BM7.3) and 12 springs only flow during the rainy season.(Figs.1c and 2a).In this spring system, BM1 is the largest spring and forms the upstream end of the Nhiem River (Fig. 2c).
Located in the Na De Valley, Na Luong Cave (NLC) is a river cave developed horizontally in the carbonate rocks of the Bac Son Fm. (C-Pbs) at an elevation of 238 m above sea level (asl) (Figs.1c and 2b).Elsewhere, Hang Rong cave (HRC) has developed in the Hong Ngai Fm. (T 1 hn), which consists of layers of argillaceous limestone, siliceous limestone, and medium-to thick-bedded limestone (Figs.1b and 2d).
Southwest springs (SW springs) are springs draining in the southwest area of DVKP.This group includes five springs namely DA7, DA9, DA11, DA12, and DA13 that appear along the NW-SE fault zone, which is also the boundary between the carbonate rocks of the Bac Son Fm. (C-P bs) and the terrigenous deposits of the Yen Binh Fm. (T 2 a yb), Ma Le Fm. (D 1 ml), as well as the shale and siliceous limestone of the Si Phai Fm. (D 1-2 sp) (Fig. 1b).
The mining area is located northwest of the Na De spring system.Based on survey data, the mined ores are mainly zincwolframite (Zn-W), which occur in the soil layers that are derived by in-situ weathering or weathering plus gravitational movement or accumulation overlying the carbonate rock erosion surface of the Bac Son Fm. (C-Pbs) and Dong Dang Fm. (P 2 đđ), or in the intermontane karst valleys and karst sinkholes.(Figs.1b and 2e,f).

Measurements and analysis
A portable multi-parameter device (Multi 3430 IDS, WTW, Weilheim, Germany) was used to measure the pH (IDS SenTix 980).The probe was calibrated on a regular basis with standard buffer solutions.Electrical conductivity (EC) and temperature (T) were measured with a TetraCon 325 sensor.
Water samples for major ion analysis were filtered with cellulose acetate membrane filters (0.45 μm, 25 mm, Sartorius AG, Göttingen, Germany).Cation samples were collected in conical centrifuge tubes (SuperClear™ 15 ml, Labcon North America, Petaluma, California, USA), while anion samples were taken using high density polyethylene bottles (30 ml, Rixius AG, Mannheim Germany).Samples for cation analyses were stabilized with nitric acid, and all samples were stored in a refrigerator before being transported to Germany for analysis.

Presentation of hydrochemical data
Gibbs diagrams were used to identify the major processes controlling natural groundwater chemistry (Gibbs 1970).A semi-logarithmic plot of the weight ratio of cations [Na + / (Na + +Ca 2+ )] or anions [Cl − /(Cl − +HCO 3 − )] versus total dissolved solids (TDS) provides information on the mechanisms controlling water chemistry, such as rock weathering, atmospheric precipitation, and evapotranspiration.All ionic concentrations are expressed in milliequivalents per liter (meq/L).In this study, Gibbs diagrams were established using Origin software (Origin 2020b).

Multivariate statistical analysis
Hierarchical cluster analysis (HCA) is a statistical technique that is proving to be a powerful grouping mechanism for determining whether samples can be grouped into distinct hydrochemical groups that may be significant in the geologic context (Güler et al. 2002;Heikinien et al. 2002;Lambrakis et al. 2004;Holland and Witthüser 2009).In this study, clustering is performed using normalized datasets (over two sampling campaigns), and the results are presented in dendrograms using Origin software.Ward's method (Ward 1963) was used to analyze the distances among linkages for the entire group of observations and the squared Euclidean distances were used to determine the distance between observations.The data were standardized by calculating their standard scores (z-scores) as follows in Eq. ( 1): where z i = standard score of the sample i, x i = value of sample, x = mean, s = standard deviation.Standardization scales the log-transformed or raw data to a range of approximately ±3 standard deviations (SD), centered about a mean of zero.Therefore, each variable has equal weight in the statistical analyses.

Normalization and calculation of REE patterns
Rare earth elements (REE) including the lanthanide group and yttrium can be separated into three subgroups according to their atomic weight: light REE (LREE, La-Sm), medium REE (MREE, Eu-Tb), and heavy REE (HREE, Dy-Lu).REE data are typically presented as distribution patterns with the individual REE, listed in the order of their atomic number, as categories on the x-axis and with the corresponding concentrations as normalized logarithmic values on the y-axis.For all water samples in this study, REE is normalized to the PAAS standard (Post Archean Australian Shales), which is a proxy for average upper continental crust (Taylor and McLennan 1985).
A comprehensive illustration of REE patterns can be identified with a quantitative expression of anomalies that are positive or negative peaks in the REE distribution patterns.Commonly occurring anomalies are cerium (Ce) anomalies and europium (Eu) anomalies, which are quantified by Ce/Ce* and Eu/Eu* ratios where Ce* and Eu* are interpolated values based on the measured concentrations of the neighbor elements of Ce and Eu, respectively, that include lanthanum (La), neodymium (Nd), samarium (Sm), and gadolinium (Gd).These ratios are calculated by the following formulas (McLennan 1989): where the subscript n stands for normalized values.Anomalies with Ce/Ce* and Eu/Eu* ratios >1 indicate "positive (1)  1989).By contrast, the Eu anomaly is strongly lithologydependent (Tricca et al. 1999;Möller et al. 2004;Han and Liu 2007); thus, water originating from any aquifer will have Eu anomalies similar to the host rock of that aquifer.In this study, the analysis data of limestones and shale samples in the Lung Cam Permian-Triassic boundary section located in the Dong Van karst plateau from a previous study by Trinh et al. 2007 was used to understand the relationship between the groundwater and exposed lithology.

Physicochemical characteristics
The physicochemical analyses of major parameters of all water samples show significant spatial and temporal differences (Table 1).
In the NDSS group and the caves group, the values of EC, TDS, Ca 2+ , and HCO 3 -are high and do not show significant variation between dry and wet seasons.The water in the SW springs group is characterized by low Ca 2+ , and HCO 3 -content.Table 1 also shows that the concentrations of major parameters in water samples collected at hanging pools in the mining area vary widely depending on the season and location of each sampling site.
The concentrations of SO 4 2-and NO 3 -ions show large differences among the sampling areas.The abnormally high value of NO 3 -found in BM2 (15.3 mg/L) is thought to be related to fertilizer pollution due to agricultural activities.Whereas the highest SO 4 2-value of 9.4 mg/L in the HRC during the dry season is probably related to the alkaline lake structure and evaporation; the oxidation of metal sulfides (e.g., FeS 2 , ZnS, Al 2 S 3 , PbS, etc.) in the mining area could be a possible source of SO 4 2-.These values are lower than the guideline value set by the World Health Organization for drinking water (WHO 2022).
The Gibbs diagrams (Fig. 3a,b) show that most of the water samples of the NDSS and the caves groups are located in the rock-weathering-dominated field.An exception is a sample taken at the recharge stream of Tia Sang Cave, which is in the dominant precipitation field in the rainy season.Although

Water-rock interaction
The interaction of water with its host environment, as well as atmospheric and human inputs, play an important role in forming the chemical compositions of groundwater.The ratios between the major-ion-equivalent concentrations further indicate the sources and the main hydrogeochemical processes that control groundwater chemistry (Gibbs 1970;Elango and Kannan, 2007;Goldscheider 2015;Gao et al. 2019).
As mentioned by Hunkeler and Mudry 2014, the dissolution of carbonate and gypsum can provide Ca 2+ , Mg 2+ , SO 4 2− and HCO 3 − in a groundwater system.If gypsum is the major source of hydrogeochemical reactions, the equivalent ratio of (Ca 2+ /SO 4 2− ) should be equal to 1. Conversely, if the weathering of carbonates (calcite and dolomite) is the major hydrogeochemical reaction, the ratios of HCO 3 − /(Ca 2+ +Mg 2+ ) should be around 1. /Ca 2+ ratios close to 0, while Fig. 4b shows that the HCO 3 − /(Ca 2+ +Mg 2+ ) ratios are distributed around a 1:1 dissolution line, except for samples from the mining area that have higher SO 4 2− contents.This means that the dissolution of calcite and dolomite is the major process influencing the water chemistry in all sampling areas, whereas the dissolution of gypsum or the oxidation of metal sulfides occurs in the mining area.
The (Mg 2+ /Ca 2+ ) index is a useful parameter that can reflect the relative proportion of calcite and dolomite in carbonate rock through which the water sample circulated (White 1988;Vu et al. 2016).A ratio of equivalent concentrations (Mg 2+ /Ca 2+ ) = 1 represents the dissolution of dolomite.In contrast, a ratio (Mg 2+ /Ca 2+ ) = 0 represents the dissolution of pure calcite.A ratio of (Mg 2+ /Ca 2+ ) = 0.5 indicates the equilibrium of water with both calcite and dolomite.As shown in Fig. 4c, most of the water samples collected at the Na De spring system, in caves, and in the mining area have Mg 2+ /Ca 2+ ratios below 0.5, suggesting dominant calcite in carbonate rocks of the areas.Meanwhile, most of the SW springs and TSC also have Mg 2+ /Ca 2+ ratios less than 0.5 and Ca 2+ content lower than 1.2 meq/L, which means that dissolution of calcite dominates in the hydrochemical reactions in these springs, but in low concentration.In addition, in this group, sample DA12 is located above the line (Mg 2+ /Ca 2+ ) = 1, suggesting that this spring flows through an aquifer that is affected by the dissolution of dolomite.
A study by Ender et al. (2018) indicates that the relation between (Ca 2+ /Mg 2+ ) and Ca 2+ can be used to illustrate the hydrochemical processes in different aquifers in the area of Dong Van.Accordingly, the Bac Son Fm. is characterized by high Ca 2+ concentration and the (Ca 2+ /Mg 2+ ) index ranges Fig. 3 Gibbs diagrams, used to explain the predominant processes of water hydrochemistry in the study area: a for cations; b for anions from ~5 to 25, meanwhile, the 'Others' group includes springs of the Song Hien, Yen Binh, and Mia Le Fm., which were defined based on their low concentration of Ca 2+ and their Ca 2+ /Mg 2+ ratios.Figure 4d shows that the samples at NDSS, NLC, and HRC are located in the Bac Son Fm. karst aquifer, while the SW springs and TSC fall into the 'Others' group.

Hierarchical Cluster Analysis
Figure 5 presents the dendrograms for 17 variables in the DVKP during both dry and rainy seasons.The nodes represent the grouped clusters, and the lengths of the branches indicate the distances between the groups.On this basis, water samples were classified into two clusters in both dry and wet seasons.Cluster 1 (CL1) is further subdivided into two subclusters SCL1a and SCL1b during the dry season (Fig. 5a), while in the rainy season, cluster 2 (CL2) was separated into two subgroups SCL2a and SCL2b (Fig. 5b).A summary of the physical-chemical data for each cluster is presented in Table 2.
The dendrogram for the dry season (Fig. 5a; Table 2) shows that cluster 1 combines the samples with signatures of calcium bicarbonate water type.In addition, this group is also characterized by elevated arsenic and aluminum concentrations indicating the influence of a combination of anthropogenic sources.The subdivision of cluster 1 is based on variations in ion sulfate and heavy metal concentrations, with SCL1b showing highly elevated sulfate, arsenic, iron and aluminum concentrations compared with SCL1a.Cluster 2 is a collection of samples characterized by low conductivity and little mineralization.These springs are typical of shallow groundwater movement, possibly in the weathered surface zone with a short residence time.
The dendrogram in Fig. 5b shows that in the rainy season cluster 1 is the largest group, with a total of 18 samples (Table 2) that display similar physicochemical parameters to SCL1a in the dry season due to the influence of the same geogenic factors.Cluster 2 is distinguished by low values of conductivity, TDS, Ca 2+ , and HCO 3 -, while water temperature and the concentrations of arsenic, iron and aluminum were higher than cluster 1. Cluster 2 was separated into two subclusters representing variations in values of pH, EC, TDS, Ca 2+ , and HCO 3 -, as well as heavy metal concentrations, with SCL2b having significantly higher arsenic and aluminum concentrations compared with SCL2a.

Presentation and interpretation of REE data
Normalized patterns of REE for the various water samples from the dry and wet season sampling campaigns are presented in Fig. 6, expressed both as raw normalized concentrations (parts a and b), and the average REE pattern for each sampling area as well as each cave (parts c and d).
In the dry season (March 2020, Fig. 6a), the karst groundwater samples at NDSS contrast significantly with the SW springs with lower ΣREE, and steeper, highly LREEdepleted REE patterns, as well as more pronounced negative Ce anomalies and stronger positive Eu anomalies.Normalized REE patterns at the cave group showed HREE enrichment at NLC, while slight enrichment at TSC and HRC.Ce anomalies are slightly negative at HRC and strongly negative at NLC and TSC, while Eu anomalies were moderately positive.In addition, samples taken from hanging pools at the mining area showed significant differences in ΣREE.Normalized REE patterns show HREE enrichment at DA5 and moderate HREE enrichment at DA6. Ce anomalies are strongly negative, while Eu anomalies are slightly negative at DA6 and moderately positive at DA5.
During the wet season (June 2020, Fig. 6b), total REE concentrations were significantly higher for the NDSS, caves, SW spring samples, and DA5 in the mining area, except for DA6, which has a lower ΣREE concentration, compared to the dry season.Normalized REE patterns from the NDSS group indicate high HREE enrichment, slightly positive Eu anomalies, and moderately negative Ce anomalies.The cave samples show high HREE enrichment, slightly negative (at NLC) to positive Eu anomalies (at HRC), and slightly negative Ce anomalies.Compared to the dry season, the SW springs exhibit slightly to highly enriched HREE, except that DA11 has flatter REE patterns.The Eu anomalies vary from mildly negative to slightly positive.The Ce anomalies are mildly positive at DA9, slightly negative at DA7, DA12 and DA13, and sharply negative at DA11.The samples at the mining site during the rainy season exhibited a fairly similar REE pattern, characterized by high HREE enrichment, slight negative Eu anomalies and moderately negative Ce anomalies.
A comparison of average REE patterns (Fig. 6c vs. Fig.6d) shows that the mining area samples, HRC, and SW springs differ least between dry and rainy seasons, whereas NLC patterns differ most.The NDSS samples and TSC also show some small differences in the Eu and Ce anomalies.Meanwhile, the average ΣREE comparison indicates that the mean total REE concentrations in the rainy season at NDSS, SW springs, NLC, and TSC were significantly higher than in the dry season, while there was only a slight increase at HRC.However, these values in the wet season were significantly lower than in the dry season in the mining area, mainly due to a sharp decrease in DA6, while they increased strongly in DA5.These results indicate that rainwater replenishment, as well as a longer water-rock interaction time, resulted in higher REE concentrations (Gill et al. 2018;Berglund et al. 2019).Meanwhile, sample DA6 in the mining area suggests that the accumulation of wastewater during mining, combined with evaporation in the long term and high pH values, significantly increased the concentration of trace elements in the dry season.

Evaluation of clustering results
For all permutations in the analysis of major ions, and timing for when the data set was expanded to include physical parameters and common trace elements, there was a clear (statistically significant) separation between samples collected from the Na De spring system (NDSS), the caves, the mining area, and the SW springs of the Dong Van karst plateau.A reasonable explanation for the similarity of physicochemical parameters between NLC and NDSS is that they belong to the same zones of recharge and discharge.Meanwhile, the SW springs and TSC were identified as typical of shallow groundwater flowing through the weathered surface zone with a short residence time in nonkarst formations.
Although Ca 2+ , Mg 2+ , and HCO 3 -concentrations do not vary much, HRC is distinct from the NDSS based on its higher SO 4 2-content.Meanwhile, the samples taken at the mining site were distinguished from the NDSS by significantly higher concentrations of SO 4 2-and heavy metals.However, the abnormally high arsenic concentrations in the mining area and NLC, combined with their location and geologic data, indicate potential connectivity between these sampling sites when arsenic is thought to be transported by groundwater from the mining area to NLC.
There is a large change in clustering patterns during the rainy season when the mining area switches to the same branch with SW springs and TSC (Fig. 5b).This is explained by the fact that in the wet season, most of the samples in Cluster 2 are located in the precipitation-dominant field (Fig. 3), with TDS, EC, and Ca 2+ concentrations lower than cluster 1 (Table 2).The mining area samples are still distinguished by their high content of heavy metals.It can be seen that these clustering results once again confirm the shallow groundwater pattern of the SW springs and TSC, and also show that rainwater is the main water source for the ore extraction process in the mining area.In summary, the clustering results are meaningful in grouping and linking water sources; however, the limited number of samples (40) and sampling campaigns (seasonal over a single hydrological cycle) pose mathematical issues which can affect cluster structure.Thus, even if hierarchical cluster analysis yields statistically significant results, these results are not significant unless they can be reasonably explained by the conceptual groundwater system model and hydrochemical dynamics, as well as geologic data.

Using the concentration of trace and rare earth elements as a natural tracer
As mentioned by Gill et al. (2018), the consistent contrast in REE features establishes a basis for quantifying the  (c-d) show the mean concentration of each group as well as each cave to emphasize the features in terms of their pattern, slope, concentrations, and variations in primary water anomalies (Ce, Eu).Water samples were collected from the Na De spring system (NDSS; orange-filled squares); from caves (caves; blue-filled circles); in the southwest of DVKP (SW springs; green-filled triangles); and at the mining area (red-filled inverted triangles) contributions of autogenic and allogenic inputs to groundwater.Therefore, REE patterns can be used as a natural tracer to provide additional information about the characteristics of aquifer and recharge sources.
The simple comparison of REE systematics exposes differences between SW springs and the karst-dominated waters (Na De spring system, NLC).Specifically, the normalized REE patterns (Fig. 6c,d) show that the karst groundwater samples have inherited their REE features from the carbonate rocks, including low ΣREE concentrations, strongly negative Ce anomalies, and high enrichment in HREEs.By contrast, the SW springs have much higher REE concentrations, a much smaller negative Ce anomaly, and a flatter overall REE pattern.
Figure 6c shows great similarities in the REE patterns of NLC and NDSS in the dry season, while in the rainy season (Fig. 6d), the negative Ce anomaly of NLC was significantly lower, and the HREE enrichment was higher than NDSS (more details in the Appendix Table 4).These results, combined with the ratios of major ion concentrations and spatial relationships strongly suggest that NDSS and NLC belong to the discharge zone of the same karst catchment, but do not preclude the possibility that there is also a connection between NLC and nearby flow systems.In addition, the Ce anomalies of NDSS and NLC become less negative, coupled with increasing total REE concentrations during the rainy season, basically reflecting the mixing with sources that were lacking negative Ce anomalies, or that exhibited Ce anomalies such as rainwaters and/or water from another catchment.
Excluding the total REE, the REE patterns between the mining area and the NLC in both seasons show great similarities, and there is only a difference in Ce anomaly in the rainy season (more details in the Appendix Table 4. Whereas, the abnormally high arsenic concentrations in the mining area and NCL combined with their location and geologic data indicate a connection between these sampling sites in the karst system when arsenic is thought to be transported by water flow from the mining area through fractures or cave systems in the unsaturated zone to NLC (Fig. 7a,c).
From a speleological and geological point of view, the HRC sample was taken at the lake (endpoint) of Hang Rong Cave, which developed in the Triassic siliceous limestone of the Hong Ngai Fm.Thus, the groundwater has the REE distribution patterns of the host rock (Cholet et al. 2019); however, the REE patterns of HRC are fundamentally different from those of the Permian-Triassic limestone and shale samples at the Lung Cam boundary (located in the Dong Van karst plateau) analyzed by Trinh et al. 2007.In addition, the major ion composition data in the water samples from HRC have similar characteristics to those from the NDSS, which belongs to the Bac Son karst aquifer (section 'Water-rock interaction'; Fig. 4d).Moreover, there are also similarities between HRC and NDSS in terms of REE patterns, except for lower negative Ce anomalies, which are thought to be the result of oxidation of Ce 3+ to Ce 4+ and subsequent stabilization of carbonate Ce 4+ complexes in alkaline lake water, the same as in the case of Lake Van (Möller and Bau 1993) or the Abert and Mono lakes (Johannesson et al. 1994).These results would indicate that the HRC has the characteristics of an alkaline lake developed in Triassic siliceous limestones (T 1 hn), but mainly associated with the Bac Son (C-Pbs) karst aquifer (Fig. 7a,b).
The TSC sample was collected at the drainage into Tia Sang Cave, which is located at the boundary between the terrigenous sediments of the Song Hien Fm. (T 1 sh) and the limestone of the Bac Son Fm. (C-Pbs).Notwithstanding the significant differences in major ion composition and water type (as presented in section 'Water-rock interaction'), there are several commonalities in the REE patterns between TSC and NDSS during the dry and wet seasons.In particular, the REE patterns in the dry season at both sampling sites show a strong negative Ce anomaly, a moderate positive Eu anomaly, and a significant HREE enrichment (Fig. 6c), while in the wet season, the REE patterns show weaker negative Ce anomalies, lower positive Eu anomalies, and flatter REE distribution patterns (Fig. 6d).These results, combined with the geologic data, lead one to speculate that TSC could be one of the recharge sources of the Na De spring karst catchment.(Fig. 7a).

Estimated water balance
Calculation of the water balance is never straightforward for a karst aquifer.It requires deep understanding of the hydrogeologic processes that took place in that karst system and the catchment area of the discharge zone.In this section, an estimated water balance is computed based on the established hydrogeologic conceptual model to provide preliminary ideas about the potential catchment area of the Na De discharge zone.The approximate water balance is calculated using the following general equation: where Q = annual average discharge of the spring (m 3 /s); R = estimated recharge (mm/year), i.e., the difference between precipitation (P in mm/year) and evapotranspiration (ET in mm/year); A = recharge area (km 2 ), i.e., the potential catchment area of the Na De discharge zone. (4) According to the available data (Table 3), this calculation includes precipitation during 2015-2020 at the climate station Dong Van District (NCHMF); potential evapotranspiration was calculated based on the hourly potential evapotranspiration (hPET) dataset for the global land surface at 0.1° spatial resolution (Singer et al. 2021) over the period 2015-2020 due to the lack of actual evaporation data; the potential catchment area is assumed to be the entire area of the Dong Van karst plateau; and the Na De spring system discharge is a one-time measurement in the wet season-Ender, Institute of Applied Geosciences (AGW), Karlsruhe Institute of Technology (KIT), unpublished data from the field trip report for Dong Van Karst Plateau, 2014.
The estimated average annual discharge of the spring is 25.6 m 3 /s, which is lower than the discharge of the Na De spring system, which was recorded at 44 m 3 /s.This means that the Na De discharge zone catchment could be much larger than the potential catchment boundary delineated in the Dong Van karst plateau.
It can be seen that all elements of this water balance are associated with a high degree of uncertainty, particularly the potential catchment area (A), and evapotranspiration (ET).However, this rough estimation could provide an overview of the magnitude and importance of the Na De spring system in the Dong Van karst plateau.Hence, detailed studies need to be implemented in the future in order to improve knowledge of this karst system.

Conclusions
This article highlights the effectiveness of the groundwater geochemistry approach in karst hydrogeologic research.In particular, this approach is very useful for preliminary studies in areas like the Dong Van karst aquifer system.This study shows that analytical results of different groundwater parameters (major ions, trace and rare earth elements) can complement each other to clarify hydrochemical processes that are occurring in a karst system.This can form the basis of a hydrogeologic conceptual model, which can inform strategies for the protection and management of karst water resources.In addition to these general conclusions, this work has led to several specific conclusions about the study area.
Multivariate statistical analysis (hierarchical cluster analysis) of major ion composition has clarified the water chemistry signatures and provided a scientific basis for comparing and partitioning water samples into distinct hydrochemical groups that may be significant in a geologic context.Specifically, the Na De spring system (NDSS) and Na Luong Cave (NLC) belong to the same discharge zone of a karst catchment of the Bac Son karst aquifer (C-Pbs), whereas the SW springs and Tia Sang cave (TSC) were identified as typical of shallow groundwater in nonkarst formations.
Trace and rare earth element concentrations and patterns can be used as potential natural tracers of groundwater origin when some processes are not revealed through conventional hydrochemical methods, such as signatures of major ions in the case of Dong Van karst aquifer system.Specifically, it is possible to show a connection between the mining area and Na Luong Cave by the abnormally high arsenic concentrations, and similarity in the REE pattern, combined with location and geologic data.Meanwhile, Hang Rong cave has the features of an alkaline lake that developed in Triassic siliceous limestones of the Hong Ngai Fm., but is mainly associated with the Bac Son karst aquifer due to similarity in the REE pattern and major ion composition.The commonalities in REE patterns between TSC and the NDSS, combined with geologic data, lead to speculation that it could be one of the recharge sources of the Na De discharge zone.
An initial water balance, while subject to uncertainties, indicates that the Na De discharge zone catchment may extend beyond the boundary of the Dong Van karst plateau.
These findings are subject to some limitations due to the number of samples and sampling campaigns (seasonal over a single hydrological cycle).However, these results are very encouraging and deliver a valuable base for a better understanding of the Dong Van karst aquifer system.Based on that, detailed studies could be considered and implemented in the future, such as an artificial tracer test and long-term monitoring of groundwater.

Fig. 1
Fig. 1 Map of the investigation area; a Location of Dong Van karst plateau in Vietnam and South China karst belt (based on the World Karst Aquifer Map, Goldscheider et al. 2020); b Map of sampling sites on top of a geo-

Fig. 2
Fig. 2 Photos showing different sampling areas in the study area: a Na De Valley, located on the boundary between the Bac Son Fm. (C-Pbs) and Song Hien Fm. (T 1 sh) (Photo: David Lagrou); b Discharge area of the underground stream at the entrance of Na Luong cave (NLC) during the rainy season (Photo: Diep Anh Tran); c BM1 spring, the main spring of Na De spring system of the SW springs group, two springs (DA7 and DA13) fall in the rock-weathering-dominated zone.The DA13 water sample is in the dominant precipitation field in the rainy season.Also in this group are three springs (DA9, DA11, and DA12) that are in the precipitation-dominated field.In addition, Fig.3also shows that the DA5 sample, which is part of the mining area, falls in the rock-weathering-dominated field, and sample DA6 is located on the boundary between the precipitation and the rock-weathering-dominated zone.It can be seen that the water source for ore extraction is mainly rainwater in the rainy season, while in the dry season, water may be obtained from other water caves near this area.Water is stored in the hanging pools and reused for the ore-extraction process.

Figure
Figure 4a indicates that most samples have SO 4 2−

Fig. 5
Fig. 5 Dendrograms classify the groundwater samples according to physical-chemical composition.a dry season; b rainy season

Fig. 6
Fig. 6 REE abundances in waters normalized to the PAAS (Post Archean Australian Shales) standard; Taylor and McLennan 1985) and divided into the samples collected during the dry season in March 2020 (a and c) and wet season in June 2020 (b and d).Panels (a-b) show all measured water data without adjustment, while panels (c-d) show the mean concentration of each group as well as

Fig. 7
Fig. 7 Conceptual model of the Dong Van karst aquifer system: a Cross-section displaying the distribution and association ability between sampling areas based on physicochemical parameters Gd n anomalies".Anomalies with Ce/Ce* and Eu/Eu* ratios <1 are called "negative anomalies".In addition, the Pr n /Yb n ratio is used as a general descriptor of the slope of the REE pattern, with the ratio Pr n /Yb n > 1 indicating LREE enrichment, and Pr n /Yb n < 1 indicating HREE enrichment (neither Pr nor Yb behaves anomalously, thus are suitable representative elements).The occurrence of positive or negative Ce anomalies is typical for oxidizing conditions where Ce is present as Ce 4+ and is, therefore, less soluble than the other REE (McLennan

Table 1
Descriptive statistics of the major physicochemical parameters of water samples in the study area

Table 2
Mean chemical values for seasonal groups distinguished by hierarchical cluster analysis.

Table 4
REE concentration (ng/L) and REE distribution pattern parameters of water samples in DVKP.Mean values are given in italic n Dry season