Method for Identifying and Estimating Karst Groundwater Runoff Components Based on the Frequency Distributions of Conductivity and Discharge

Abstract

Identification and estimation of groundwater runoff components in karst groundwater systems to improve understanding of karst water circulation and water-rock interactions is essential for water resources assessment and development. A Gaussian mixture model is presented for identifying and estimating karst groundwater runoff components based on the frequency distributions of conductivity and discharge. Successful application of this method in the Heilongquan karst spring in South China showed that groundwater runoff components can be divided into 6–8 grades, corresponding to the grades of groundwater in karst fissures. The conductivity and discharge thresholds dividing fast and slow flow were determined to be approximately 300 μS cm⁻¹ and 0.3 m³ s⁻¹, respectively, with fast flow exhibiting lower conductivity and larger discharge. On an annual basis, fast flow occurred 9% of the time and accounted for 35% of total water volume. The results of the method compared favorably to that of hydrograph recession analysis. Estimation of groundwater runoff components based on frequency distributions of conductivity and discharge provides a novel alternative method for the quantitative evaluation of karst water resources.

1. Introduction

In all around of the world, karst aquifers are important freshwater resources and their development and protection requires improved quantitative understanding of groundwater runoff components in karst groundwater systems [1]. Identification and estimation of groundwater runoff components in karst groundwater systems to improve understanding of karst water circulation and water-rock interactions is essential for water resources development and protection. However, the heterogeneity of karst aquifers poses a barrier to this knowledge [2]. Traditional methods such as field surveys, fracture measurement, cave detection, hydrogeological drilling, and geophysical methods [3,4] are often costly, whereas separating karst groundwater runoff components through tracer tests, isotope analysis and hydrograph recession analysis [5,6,7,8] have some spatial and temporal scale limitations. High-frequency online monitoring of conductivity and discharge has been widely used in karst hydrogeological surveys [7,8,9,10,11]. A method for estimating groundwater runoff components in karst groundwater systems by the frequency distribution of conductivity has previously been proposed [12], thereby providing a novel auxiliary method for identifying the structure of karst groundwater systems.

The conductivity is a measure of the ability of water to pass an electrical current; it reflects the amount of total dissolved solids in the water, which in karst water systems is controlled mainly by calcium carbonate equilibria [12]. Rocks of different ages and lithology in the groundwater systems have different dissolution resistances and can provide different ion values to water, so the water ion balance varies from different system [13,14,15]. Hydrogeochemical modeling was often applied to recreate genesis and determine groundwater flow routes [16,17,18], but it was costly for the data. Although the conductivity is not directly equivalent to the concentrations of the major ionic species, it always was used to study water components as a synthesized parameter [12,19,20,21,22,23].

Conductivity can act as a comprehensive indicator of the hydrogeochemical characteristics of groundwater in karst water systems, and its variation reflects the contribution of different runoff components [19,24]. White et al. [20] used conductivity analysis to determine the proportions of water lost and infiltrating a river, and the water stored in the fracture-matrix. Neilson et al. [21] employed hydrogeochemical and hydrological analysis to quantify the interaction between karst aquifers and rivers, and estimated the contribution of conduit water and fissure water to a karst system. Previous studies [12,19] showed that the characteristics of groundwater runoff components in karst water systems can be analyzed by the frequency distribution of conductivity. Multiple peaks in conductivity have been identified in most karst groundwater systems through plotting the frequency distribution of conductivity [22]. Each peak represents a different water source and frequency, which can also reflect the degree of karstification [23]. This method has been used in several case studies of karst groundwater systems [13,25,26]. Although this method can determine the frequency of different groundwater runoff components in the karst water system, quantification of the contribution of each component to total water volume remains problematic.

The present study applied a novel Gaussian mixture model based on the frequency distribution of conductivity as proposed by Massei et al. [12] and the frequency distribution of discharge to Heilongquan Spring, a typical karst groundwater system in South China, for identifying and estimating karst groundwater runoff components.

2. Field Site Descriptions

The Heilongquan Spring is located in Xingshan County, western Hubei Province, Southern China (Figure 1). The study area falls within the northwestern wing of the Zigui anticline, and is characterized by a subtropical monsoon climate with abundant rainfall of 1100 mm a⁻¹, 68% of which occurs during the April–September rainy season [3]. Heilongquan Spring discharges from an altitude of 410 m–550 m a.s.l. whereas its recharge is from the northern area at 900 m–1100 m a.s.l., with a 17.8 km² catchment area faring 4–5 km away from it. While the geological formations lie relatively flat, the topography undulates sharply. Geological strata from the Sinian (late Precambrian) to the Silurian are exposed, with thick units of limestone and dolostone from the Sinian, Cambrian, and Ordovician ages widespread in the study area (Figure 1c). Non-carbonate rocks of mostly mudstone and shale of the Lower Cambrian age are also present. The Heilongquan karst groundwater system is the main aquifer of this area, and is composed of interconnected Ordovician limestone, Lower Cambrian limestone, and Middle-to-Upper Cambrian dolostone. This aquifer lies above a regional aquiclude constituted of Lower Cambrian mudstone and shale (Figure 2). The boundary of the Heilongquan karst groundwater system has been confirmed by field surveys, geophysical exploration, tracer tests, hydrogeochemical analysis, and groundwater isotope analysis [27].

3. Materials and Methods

3.1. Gaussian Mixture Model

A frequency distribution sorts a time series of values into groups of ranges, and shows the frequency of occurrence of values in the time series into each group [28]. A frequency distribution can be represented as a histogram, X-Y plot or table. Depending on the probable type of frequency distribution of the data, the shape of the frequency distribution graph can be approximated by an algorithm, which allows the mathematical characteristics of each group of the data to be obtained. In groundwater systems, conductivity and discharge are cumulative reflections of different runoff components in the system. Across all karst systems, an individual groundwater runoff component type would have a similar pathway, and the frequency distribution of conductivity of this component would conform to a single Gaussian model [12,19]. Then the frequency distribution of conductivity of a spring conforms to the Gaussian mixture model, the probability density equation can be expressed as:

$$p\left(x\right)={{\displaystyle \sum }}_{k=1}^K{\epsilon }_k·N\left(x|{\mu }_k,{\sigma }_k^2\right)$$

(1)

where $N\left(x|{\mu }_k,{\sigma }_k^2\right)$ is the normal distribution density equation of any runoff component; ${\mu }_k$ is the mathematical expectation of the runoff component; ${\sigma }_k$ is the standard deviation of the runoff component; and ${\epsilon }_k$ is the weight of the runoff component, which satisfies ${{\displaystyle \sum }}_{k=1}^K{\epsilon }_k=1$.

After completion of the basal statistical analysis for conductivity data, the normal distribution peaks of different groundwater runoff components are obtained using the peak-fitting software program PeakFit 4.12 (SPSS Inc., Chicago, IL, USA) to analyze the conductivity data in this study. Through combined analysis with meteorological, hydrological, and geological characteristics of the study area, the different groundwater runoff components can be identified. It should be noted that conductivity is a reflection of ionic charge rather than ion concentration; therefore, the proportion of conductivity represented by each peak does not represent the proportion of the volume of the runoff components [12]. However, the relative proportion of each peak represents the proportion of time that the groundwater runoff component appears over the study period, hereinafter referred to as “time proportion”, represented by weight “${\epsilon }_k$” in Equation (1).

3.2. Calculation of Groundwater Runoff Components

Since the sources of recharge in small and medium karst groundwater systems are often concentrated in a certain area, the difference in runoff pathways is the main reason for differences in conductivity and discharge [8,13,29]. Groundwater runoff components that rapidly pass through karst conduits and large-scale fractures are characterized by short stagnation periods, weak interaction with rocks, low conductivity, and large discharge [14,29,30]. In contrast, groundwater runoff components that pass more slowly through small fractures and matrix are characterized by long periods of stagnation, strong interaction with rocks, high conductivity, and small discharge [14,29,30]. Similar to conductivity, it can be assumed that discharge corresponding to a single runoff component is normally distributed, and therefore that the frequency distribution of discharge conforms to a Gaussian mixture model. In addition, the discharge and conductivity frequency distributions show similar characteristics. Under the assumption that the mathematical expectation in the discharge Gaussian mixture model represents the mean discharge of the runoff component and by calculating the duration of the runoff component in time proportion of the peak, the volumes of each groundwater runoff component can be calculated.

Figure 3 shows the frequency distribution of conductivity and discharge for a karst groundwater system over time T. For a certain groundwater runoff component P1, according to the Equation (1), the time proportion of the peak P1 in the conductivity frequency distribution graph is recorded as ${\epsilon }_c$, and the mathematical expectation of the runoff component in the discharge frequency distribution graph is recorded as ${\mu }_q$. The volume of the P1 groundwater runoff component in time T can then be expressed as:

$$V={\mu }_q·{\epsilon }_c·T$$

(2)

3.3. Hydrograph Recession Analysis

Previous studies [3,7,31,32,33] of relatively independent karst aquifer systems only recharged by rainfall concluded that discharge originates only from the water storage in the system under periods of no rainfall. During the initial stage of recession, karst water is drained from fractures of all sizes, although the recession time and recession coefficient are different for the different fracture sizes. The karst water drainage process is regarded as the superposition of the drainage process of different grades of karst fractures. The entire recession process can be categorized into several periods according to the differences in recession coefficient. An exponential equation can be used to describe the recession process:

$$Q_t=\{\begin{array}{c}{Q}_1·e^{-{\alpha }_{1}t} \left[0,t_1\right]\hfill \\ Q_2·e^{-{\alpha }_{2}t} \left[t_1,t_2\right]\hfill \\ Q_3·e^{-{\alpha }_{3}t} \left[t_2,t_3\right]\hfill \\ Q_4·e^{-{\alpha }_{4}t} \left[t_3,t_4\right]\hfill \\ \dots \dots \hfill \end{array}$$

(3)

where Q_t is the discharge at any moment; Q_i is the discharge at the initial moment of each recession period; α_i is the recession coefficient (coefficient of discharge) with dimension (T⁻¹), and $t_i$ is the end time of each recession period.

According to the duration and the coefficient of each recession period, the volume of each groundwater runoff component can be calculated by an integral equation (Figure 4), expressed as:

$$\{\begin{array}{c}V1={{\displaystyle \int }}_0^{t_1}\left(Q_1{e}^{-{\alpha }_{1}t}-Q_2{e}^{-{\alpha }_{2}t}\right)dt\\ V2={{\displaystyle \int }}_0^{t_2}\left(Q_2{e}^{-{\alpha }_{2}t}-Q_3{e}^{-{\alpha }_{3}t}\right)dt\\ V3={{\displaystyle \int }}_0^{t_3}\left(Q_3{e}^{-{\alpha }_{3}t}-Q_4{e}^{-{\alpha }_{4}t}\right)dt\\ V4={{\displaystyle \int }}_0^{\infty }\left(Q_4{e}^{-{\alpha }_{4}t}\right)dt\end{array}$$

(4)

In addition, the first and second recession periods are often considered to be fast flow, whereas the subsequent recession periods are often considered to be slow flow [32]. This method can be used to analyze all rainfall events in the study period and to calculate the volume of each groundwater component, and can therefore be used to verify the results of the frequency distribution.

3.4. Data and Materials

Water level (accuracy of 0.001 m), water temperature, and conductivity (accuracy of 0.1 μS cm⁻¹) of the Heilongquan Spring were monitoring every 15 min by an auto level logger (Model 3001 LTC Level logger, Solinst, Georgetown, ON, Canada) from 2014 to now. Due to the technological problem, some data in 2015 were missed, and hydrological years (2016–2018) with complete time series were analyzed. After the correction of the measurement results in Xiangxi river basin, the rating curve of Heilongquan Spring was established through measurements of the discharge at different water levels, and the monitored water level was converted to discharge, expressed in Equations (5) and (6).

$$Q=\frac{A}{n}·R^{2/3}·{i_h}^{1/2}$$

(5)

$$i_h=\{\begin{array}{c}0.001 \left[0,0.234\right]\\ 0.001\left(h-0.234\right)+0.001 \left[0.234,\infty \right]\end{array}$$

(6)

where Q is the discharge at any moment; A is the cross section of flow; n is the roughness coefficient; R is the hydraulic radius; i is the hydraulic gradient, and h is the water level at any moment.

4. Results and Discussion

4.1. Conductivity and Discharge of Heilongquan Spring

According to the conductivity and discharge for hydrological years 2016–2018 (Figure 5), the spring discharge in 2016 was predominantly below the median, and conductivity responded well to rainfall events. Rainfall in July–August resulted in higher than median spring discharge and the lowest conductivity in the year. The 2017 hydrological year was relatively wet, with rainfall in April–October often resulting in higher than median spring discharge and conductivity reached the lowest level in three years after the rainfalls in July. The 2018 hydrological year was relatively dry, with consistently high conductivity and spring discharge reaching its lowest level in September.

4.2. Identification of Groundwater Runoff Components

After analysis of conductivity and discharge data for hydrological years 2016–2018; seven, eight and six data peaks were processed by the Peakfit 4.12 software program in 2016, 2017, and 2018, with the goodness of fit (r²) larger than 0.75 and standard deviations (δ) less than 0.005, respectively (Figure 6). The mathematical characteristics of the Gaussian mixture models of conductivity and discharge are shown in

Table 1. Gaussian mixture models and the mean annual discharge of Heilongquan Spring.

Water Year		P1	P2	P3	P4	Pslow	P5	P6	P7	P8	Pfast	Qm/(m3·s−1)	r2	δ
Gaussian Mixture Model of the Conductivity
Time proportion εc/%	2016	53.24	18.05	12.79	7.24	91.31	5.12	3.29	0.27	/	8.69	0.1618	0.89	0.0038
	2017	35.23	26.16	14.56	10.13	86.08	5.91	3.80	2.95	1.27	13.92	0.1897	0.82	0.0045
	2018	34.02	26.07	19.24	15.42	94.75	3.34	1.91	/	/	5.25	0.1413	0.84	0.0043
Mathematical expectation μq/(μS·cm−1)	2016	352.18	341.62	328.68	318.59	344.14	293.82	272.35	245.58	/	284.19	0.1618	0.89	0.0038
	2017	351.18	339.59	325.55	307.04	338.13	277.54	247.03	220.38	195.57	249.64	0.1897	0.82	0.0045
	2018	359.02	348.54	334.50	321.07	344.98	296.26	276.73	/	/	289.16	0.1413	0.84	0.0043
Gaussian Mixture Model of the Discharge
Time proportion εc/%	2016	30.37	24.59	20.44	16.30	91.70	4.44	2.37	1.48	/	8.30	0.1618	0.87	0.0019
	2017	36.75	28.03	14.48	7.61	86.87	5.38	3.71	2.78	1.25	13.13	0.1897	0.76	0.0049
	2018	30.16	25.40	21.90	17.14	94.60	3.49	1.90	/	/	5.40	0.1413	0.79	0.0024
Mathematical expectation μq/(m3·s−1)	2016	0.07	0.13	0.20	0.25	0.13	0.55	0.82	1.06	/	0.72	0.1618	0.87	0.0019
	2017	0.07	0.14	0.20	0.27	0.13	0.57	0.83	1.10	1.27	0.82	0.1897	0.76	0.0049
	2018	0.04	0.10	0.15	0.20	0.11	0.50	0.76	/	/	0.59	0.1413	0.79	0.0024

The number of groundwater runoff components identified are different in each water year, and “/” indicates that no runoff components of this grade are present. The goodness of fit (r²) and standard deviations (δ) are about the frequency of conductivity and discharge measured and estimated.

.

The runoff components identified were consistent between the frequency distributions of conductivity and discharge. Although the shape of the peak of each runoff component identified by the conductivity Gaussian mixture model was consistent with that of the discharge Gaussian mixture model, the number of peaks varied from hydrological year 2016 to 2018. This can be explained by the recharge area of Heilongquan Spring only being concentrated in the northern part of the study area [27], which results in the different peaks not being representative of the different sources of recharge. The differences in conductivity and discharge was caused by the differences in runoff pathways, corresponding to the different grades of karst fissures in the Heilongquan karst groundwater system [32]. Therefore, it can be considered that different peaks are representative of groundwater components released from different grades of karst fissures in the system.

In the 2016 hydrological year, seven groups of groundwater runoff components were identified. Moderate total rainfall occurred in the year, with coexistence of groundwater in different grades of karst fissures. Therefore, the peaks of the runoff components identified by the frequency distributions of conductivity and discharge in 2016 were independent, and the number and shape of peaks identified by the two parameters were similar (Figure 6). There are several grades of karst fissures in the development of Heilongquan karst groundwater system, such as small fissures, medium fissures, wide fissures, and karst conduits [32]. It can be considered that the P1 to P4 peaks with higher conductivity and lower discharge represented the runoff components released from small and medium karst fissures, with an overall time proportion of over 91%. The P5, P6, and P7 peaks with lower conductivity and higher discharge represented the runoff components released from karst conduit and medium-wide to wide fissures, and their time proportion was less than 9% (Table 1).

The 2017 wet hydrological year was characterized by a large total amount of rainfall, high number of rainfall events, and multiple rainfall grades. Thus, eight groups of runoff components were identified, and the extra P8 peak represented the karst cave water with larger discharge, which only appeared during heavy rainfall. Due to the impact of higher rainfall, the overall time proportion of P5–P8 runoff components in the 2017 hydrological year reached 14%, whereas the proportion of the P1–P4 components fell to 86% (Table 1).

The 2018 dry hydrological year was characterized by lower total rainfall, fewer overall rainfall events and less rainfall grades. Thus, only six runoff components were identified, indicating that the runoff components released from karst cave and conduit did not manifest in this year, and other peaks still corresponded to the groundwater in different grades of karst fissures. In addition, the time proportion of P5–P6 runoff components in the 2018 hydrological year sharply dropped to 5%, whereas the proportion of the P1–P4 components rose to 95% (Table 1).

The shapes of the peaks in the conductivity and discharge frequency distributions were consistently flat in 2017 hydrological year, whereas they were consistently sharp in the 2018 hydrological year (Figure 6), which is mainly the result of rainfall. During the 2017 hydrological year, since there was a large total amount of rainfall and continuous rainfall grades, the groundwater runoff components in the system were continuously released during groundwater recession. The duration of the multiple groundwater runoff components released together were longer than those in other years, whereas the durations for single groundwater runoff component released independently were shorter, which resulted in a wider peak shape of groundwater runoff components and cross section area between the peaks in 2017. Since the total rainfall amount was lower in the 2018 hydrological year and the rainfall grades were differential, the durations of single groundwater runoff components released independently were longer than those in other years, which resulted in narrower peak shapes of groundwater runoff components and the cross-section areas between the peaks.

The mean conductivity and discharge of the P1–P4 groundwater runoff components released from small and medium karst fissures were consistently larger than 300 μS cm⁻¹ and less than 0.3 m³ s⁻¹, respectively (Figure 6 and Table 1). In contrast, those of the P5–P8 components released from medium-wide fissure and karst conduit were consistently less than 300 μS cm⁻¹ and larger than 0.3 m³ s⁻¹. Moreover, the groundwater runoff components with conductivity above 300 μS cm⁻¹ occurred at an annual frequency of a factor of 19–27 times that of the groundwater runoff components with conductivity below 300 μS cm⁻¹. The groundwater runoff component with discharge below 0.3 m³ s⁻¹ occurred at an annual frequency of a factor of 21–28 times that of the groundwater runoff components with discharge above 0.3 m³ s⁻¹ (Figure 6 and Figure 7). Therefore, 300 μS cm⁻¹ could represent a threshold of conductivity for the division of groundwater runoff components into fast flow below the threshold released from medium-wide and wide fissures, karst conduit, and karst cave, and slow flow above the threshold released from small and medium karst fissures. In addition, discharge of 0.3 m³ s⁻¹ of discharge could represent a threshold, with fast flow above the threshold and slow flow below the threshold.

The peaks of slow flow (P1–P4) in the Heilongquan karst groundwater system overlapped and tended to form a peak-cluster (Figure 6), particularly in 2018. In contrast, peaks of fast flow (P5–P8) in the system were independent. This was related to rainfall and grades of fissures. The slow flow components were continuously released from the small-medium karst fissures. The duration of multiple groundwater runoff components released together was longer than that of fast flow components, resulting in difficulty distinguishing between the components by conductivity and discharge. The fast flow groundwater runoff components were released from wide fissures, karst conduit, and karst cave, and they were significantly influenced by heavy rainfall. Thus, their peaks were relatively independent, and therefore could be easily distinguished by conductivity and discharge. Furthermore, during the wet 2017 year, abundant rainfall resulted in a shift of all peaks to the left with lower conductivity and larger discharge. During the dry 2018 year, low aquifer levels resulted in longer residence times and a shift of all peaks to the right with higher conductivity and smaller discharge. In addition, the average conductivity and discharge of fast flow were greatly affected by rainfall intensity, whereas the average conductivity and discharge of slow flow remained relatively stable (Figure 7). This can be explained by the fast flow having a shorter circulation path and residence time, and therefore being more susceptible to rainfall events, whereas slow flow has a longer flow path, a larger circulation depth, and a longer residence time, and is therefore less influenced by rainfall.

Additionally, the time proportion of fast flow and slow flow during 2017 were approximately 5% larger and 5% lower than that of 2016, respectively. The time proportion of fast flow and slow flow in 2018 was approximately 3% lower and 3% higher than that in 2016, respectively. However, regardless of hydrological conditions, the Heilongquan karst groundwater system consistently maintained slow flow as the main flow component, with its average time proportion over 2016–2018 being 91%.

4.3. Calculation of Groundwater Runoff Components

Over the three hydrological years, the peaks of groundwater runoff components in the conductivity frequency distribution showed greater independence than those in the discharge frequency distribution, indicating that the effects of the frequency distribution of conductivity were more easily identifiable. This indicates that the time proportion of each runoff component should be determined using the frequency distribution of conductivity, whereas each peak in the frequency distribution of discharge indicates the average discharge of each groundwater runoff component. The annual runoff of each groundwater runoff component in Heilongquan Spring for water years 2016 to 2018 was calculated according to the method described above, with results shown in

Table 2. The annual runoff of groundwater runoff components in Heilongquan Spring.

Water Year	2016				2017				2018
Slow flow	Time proportion/%	Average discharge/(m3·s−1)	Runoff volume/(104 m3)	Runoff proportion/%	Time proportion/%	Average discharge/(m3·s−1)	Runoff volume/(104 m3)	Runoff proportion/%	Time proportion/%	Average discharge/(m3·s−1)	Runoff volume/(104 m3)	Runoff proportion/%
P1	53.24	0.07	114.20	22.41	35.23	0.07	82.55	11.27	34.02	0.04	47.19	11.37
P2	18.05	0.13	74.37	14.59	26.16	0.14	112.39	15.34	26.07	0.10	78.95	19.01
P3	12.79	0.20	81.67	16.02	14.56	0.20	92.70	12.65	19.24	0.15	91.81	22.11
P4	7.24	0.25	57.69	11.32	10.13	0.27	85.98	11.74	15.42	0.20	98.93	23.83
Total up	91.31	0.10	327.93	64.34	86.08	0.12	373.61	51.00	94.75	0.10	316.88	76.32
Fast flow	Time proportion/%	Average discharge/(m3·s−1)	Runoff volume/(104 m3)	Runoff proportion/%	Time proportion/%	Average discharge/(m3·s−1)	Runoff volume/(104 m3)	Runoff proportion/%	Time proportion/%	Average discharge/(m3·s−1)	Runoff volume/(104 m3)	Runoff proportion/%
P5	5.12	0.55	88.09	17.28	5.91	0.57	106.58	14.55	3.34	0.50	52.60	12.67
P6	3.29	0.82	84.69	16.62	3.80	0.83	99.22	13.54	1.91	0.76	45.72	11.01
P7	0.27	1.06	8.97	1.76	2.95	1.10	102.55	14.00	/	/	/	/
P8	/	/	/	/	1.27	1.27	50.67	6.92	/	/	/	/
Total up	8.69	0.06	181.76	35.66	13.92	0.11	359.02	49.00	5.25	0.03	98.32	23.68
Total of all	100.00	0.16	509.69	100.00	100.00	0.23	732.63	100.00	100.00	0.13	415.20	100.00

The number of groundwater runoff components identified are different in each water year, and “/” indicates that no runoff components of this grade are present.

and Figure 8.

According to results, annual runoff of Heilongquan was 4.5 Mm³–7.5 Mm³. The annual runoff of Heilongquan and annual proportion of fast flow increased with increasing annual rainfall, whereas the annual proportion of slow flow decreased. The annual runoff of slow flow remained stable at approximately 3.5 Mm³, whereas the annual runoff of fast flow between wet and dry years showed higher variation, and was determined by the characteristics of the aquifer medium. The volumes of small and medium fissures contributing to slow flow were constant in the system. Since the groundwater in these small and medium fissures occurred at greater depths, the fissures remained in the phreatic water zone, even in dry years, resulting in the runoff of slow flow being less affected by rainfall. Therefore, the discharge from small and medium fissures represented the base flow of the groundwater to a certain extent. In contrast, the geological medium contributing to fast flow, such as wide fissures, karst conduits, and karst caves, occur in in the vertical infiltration zone or the seasonal variation zone in most cases, which is significantly influenced by rainfall. Moreover, although the average annual time proportion of the fast flow runoff components of Heilongquan Spring was 9%, the average proportion of total volume contributed by fast flow runoff reached 35%, which indicated that although the karst cave flow and karst conduit water occurred less frequently, it accounted for a significant proportion of total water volume.

4.4. Verification of Calculation Results

All rainfall events were analyzed hydrograph recession, and the annual runoff volumes of fast flow and slow flow components were calculated by the exponential equation. By comparing the results of the frequency distribution with the recession analysis, the applicability of the frequency distribution method for identifying and estimating karst groundwater runoff components can be assessed. The calculation results are shown in

Table 3. Groundwater runoff component resources and error statistics of Heilongquan Spring.

Water Year	2016			2017			2018
Runoff Volume/(104 m3)	QF	QS	QT	QF	QS	QT	QF	QS	QT
Frequency distribution	181.76	327.93	509.69	359.02	373.61	732.63	98.32	316.88	415.20
Recession analysis	169.70	311.41	481.11	328.25	350.60	678.85	85.13	302.57	387.70
Relative error/%	6.63	5.04	5.61	8.57	6.16	7.34	13.42	4.52	6.62

Q_F is the fast flow volume, Q_S is the slow flow volume, Q_T is the annual total runoff volume.

and Figure 9.

When comparing results between the two methods, the relative errors of annual runoff, fast flow runoff and slow flow runoff were less than 8%, less than 14%, and less than 7%, respectively. The methods were comparable to each other, only 10% error in hydrogeology is pretty minimal, indicating that the frequency distribution method performed well.

The runoff volumes of fast flow, slow flow and total flow calculated by the frequency distribution method were larger than that of the hydrograph recession analysis, with the relative error of fast flow runoff being larger than that of the slow flow analysis (Table 3, Figure 9). This was mainly due to limitations of the hydrograph recession analysis. The exponential equation of the recession analysis can only describe the retreating phase of the complete hydrologic processes, whereas the runoff of the increasing phase and the peak phase was calculated by the exponential equation of the last stage of the latest hydrologic processes, which resulted in runoff of fast flow being ignored during the two stages. This inevitably led to the hydrograph recession analysis calculating a smaller runoff, and was the main reason for the larger relative error of fast flow runoff in the validation. Moreover, the proportion of slow flow in Heilongquan was larger, which made it easier to identify and calculate karst groundwater runoff components by either method. However, since the proportion of fast flow in the system was small, some deviations in the frequency distribution fitting would occur, and the division of each recession stage in the hydrograph recession analysis can also be affected by human error, which can explain the above results.

The annual variability in rainfall also contributed to the change of relative error of runoff (Figure 9). The fluctuation in error of slow flow runoff was smaller than that of fast flow, and the relative error of slow flow runoff decreased with decreasing rainfall. The fluctuation of error in fast flow runoff was large, with the relative errors in the 2016 and 2018 hydrological years being the smallest and greatest across the study period, respectively. This was because slow flow is released from small and medium fissures at greater depths which remain in the phreatic zone regardless of hydrological conditions. Therefore, slow flow was less affected by rainfall, and the fitting effect was better compared to that for fast flow, resulting in a smaller fluctuation of error for slow flow runoff. On the other hand, the increase in rainfall inevitably increased runoff components with large discharge. However, too much rainfall resulted in peaks which were too steep to form a normal distribution and consequently resulted in a worse fitting effect and increased error.

In summary, the present study has validated the method for identifying and estimating karst groundwater runoff components based on the frequency distributions of conductivity and discharge, indicating that the method has good applicability to small and medium karst groundwater systems. Thereby the method can be used for the identification and calculation of groundwater runoff components in other similar karst groundwater systems developing several grades of karst fissures.

5. Conclusions

Groundwater runoff components in the Heilongquan karst groundwater system can be divided into 6–8 grades, corresponding to the grades of water in karst fissures. An increasing number of runoff component was identified with increasing rainfall. A conductivity threshold of 300 μS cm⁻¹ was identified separating lower conductivity fast flow released from karst cave, karst conduit, or wide fissures, and higher conductivity slow flow released from small and medium-sized karst fissures. A flow threshold of 0.3 m³ s⁻¹ of discharge was identified separating higher discharge fast flow and lower discharge slow flow.

The Heilongquan karst groundwater system consistently maintained slow flow as the predominant flow component, with its average annual time proportion for 2016–2018 being 91%. While the fast flow runoff component of Heilongquan Spring occurred at an annual frequency of 9%, the average proportion of total yearly volume contributed by fast flow runoff reached 35%, indicating that fast flow contributes significantly to water resources. The annual runoff of slow flow in the Heilongquan karst groundwater system remained stable at approximately 3.5 Mm³. This can be explained by the constant contribution of small and medium fissures to slow flow regardless of hydrological conditions due to their location at depth in the saturation zone. Thus, slow flow represents the base flow of groundwater to a certain extent.

A method based on the Gaussian mixture model for identifying and estimating karst groundwater runoff components based on the frequency distributions of conductivity and discharge is reasonable for small and medium karst groundwater systems. The method can be used for the identification and calculation of groundwater runoff components in other similar karst groundwater systems developing several grades of karst fissures.