Method for calculating surface karst water storage of karst underground river basin
By collecting and analyzing the physical, chemical, and isotopic indicators of karst water samples, and using EMMA and MIX iterative calculations to separate karst water flow rates, the high cost and low accuracy problems of surface karst water storage calculation in karst underground river basins have been solved, achieving low-cost, high-accuracy storage calculation and long-term monitoring.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2023-02-10
- Publication Date
- 2026-03-24
AI Technical Summary
Existing technologies are difficult to use efficiently and cost-effectively to calculate the surface karst water reserves in karst underground river basins, especially in basins with high karstification levels, where the accuracy is poor and the workload is large.
By collecting samples of surface karst water, saturated zone fissure water, and surface water, and testing their physical, chemical, and isotopic indices, the relationships between representative components were established. Using EMMA analysis and MIX iterative calculations, the flow process curve of surface karst water was separated, and the reserves were calculated by combining the monitoring data of underground river outlet flow.
It enables low-cost and high-precision calculation of surface karst water reserves, reducing workload and improving accuracy, and is suitable for long-term monitoring of seasonal and interannual variations of karst water.
Smart Images

Figure CN116227178B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of water resource exploration technology, and in particular to a method for calculating the surface karst water reserves in karst underground river basins. Background Technology
[0002] Karst water is an important water resource, playing a vital role in global economic and social development and the ecological environment. Approximately 25% of the world's population relies on karst water for drinking water. China has the widest distribution area and the most diverse types of karst formations in the world, and it also has the most widespread exploitation of karst water resources and the most related environmental problems. In the karst regions of southwestern my country, under the influence of subtropical monsoon climate conditions, with simultaneous rainfall and heat and favorable hydrothermal conditions and strong karst dynamics, there exists a fractured zone in the shallow part of the surface. This zone exhibits superior karst dynamic conditions, intense karstification, and a high degree of karstification, with a relatively intact carbonate rock base. This karst development zone is called the surface karst zone. In areas with well-developed surface karst zones, intense karstification creates a dual karst water cycle, coupling shallow karst water circulation with deep underground conduit circulation. This forms a dual hydrogeological structure of surface and underground. Groundwater in the surface karst zone drains through two pathways: shallow surface karst springs and deep underground river channels via shafts and sinkholes. Surface karst water is easy and economical to develop, readily available for use by residents in karst mountain areas, effectively addressing engineering-related water shortages. Understanding surface karst water reserves is crucial for the coordinated development and rational protection of groundwater resources.
[0003] Currently, there are few methods for calculating surface karst water reserves. Existing technical methods mainly fall into two categories:
[0004] Geophysical exploration methods, including ground-penetrating radar, spontaneous potential, magnetotellurics, gravity, and acoustic detection, work by utilizing the differences in physical properties of different media to determine the thickness of the surface karst zone, and then using this thickness to estimate the surface karst water reserves. However, the channels and individual spaces where groundwater is stored in the surface karst zone are small and their distribution patterns are complex, requiring high-resolution geophysical exploration techniques for effective identification of surface karst water, which is costly.
[0005] Surface karst spring survey: This method estimates surface karst water reserves by investigating and measuring the flow rates of all surface karst springs. However, this method involves a large workload and is tedious. Furthermore, because karst areas have a dual structure with shallow karst water circulation and deep underground pipeline water circulation coupled together, surface karst springs are less common in highly developed karst underground river basins. Surface karst water is mostly discharged into underground river pipelines through sinkholes and shafts. Therefore, this survey method is not suitable for watersheds with high karst development and has poor accuracy. Summary of the Invention
[0006] The purpose of this invention is to address the aforementioned shortcomings of existing technologies by proposing a low-cost, convenient, and highly accurate method for calculating the surface karst water storage in karst underground river basins. Guided by karst water system theory, this method uses parameters such as hydrochemistry and isotopes to separate the surface karst water flow process curve from the flow process curve after a rainstorm in an underground river, thereby calculating the surface karst water storage.
[0007] The present invention provides a method for calculating the surface karst water storage in a karst underground river basin, comprising the following steps:
[0008] S1: Collect surface karst water (EW), saturated zone fissure water (FW), and surface water (PW) samples, test their physical, chemical, and isotopic indicators, analyze and determine their representative components, and establish the relationship between the representative components of the three samples.
[0009] S2: Select periods of heavy rainfall to continuously monitor the flow rate, physical, chemical, and isotopic indicators of the underground river outlet;
[0010] S3: Perform EMMA analysis on the monitoring data of the underground river outlet in step 2 to determine the final physical, chemical and isotopic components used in the mixing ratio calculation;
[0011] S4: Treat EW, FW and PW as three endmembers, recalculate the final component values of the components determined in step S3 using MIX iteration, and use the representative component relationship of the three types of water established in step S1 as a constraint.
[0012] S5: Calculate the mixing ratio of the three endmembers in each sample at the underground river outlet using the component values calculated in step S4 and the monitoring data of the final components determined in step S3.
[0013] S6: Reconstruct the mixing ratio calculated in step S5 into the flow data of the three terminal elements to obtain the flow process curve of each terminal element;
[0014] S7: By integrating the EW flow process curve obtained in step S5, the total discharge of EW during the rainstorm event is obtained, which is equal to the surface karst water storage.
[0015] Furthermore, the specific operation of step S3 is as follows: perform EMMA analysis on the monitoring data of the underground river outlet in step 2 to analyze whether two principal components can be extracted. If two principal components can be extracted, the monitored components are used as the final physical, chemical, and isotopic components used to determine the mixing ratio calculation. If more than two principal components are extracted, all analyzed components are screened, and the representative components of EW, FW, and PW determined in step 1 are retained first. Then, EMMA analysis is performed until two principal components are extracted. The finally screened components are used as the final physical, chemical, and isotopic components used to determine the mixing ratio calculation.
[0016] Furthermore, in step S6, the proportion of each end-member in each sample of the underground river water is multiplied by the underground river flow rate at the corresponding time to obtain the flow process curve of each end-member.
[0017] Furthermore, in step S1, the typical physical, chemical, and isotopic characteristics that distinguish EW from PW and FW are manifested in the representative components.
[0018] Furthermore, in step S1, the magnitude relationship of the test values of the representative components of the three samples is established.
[0019] Changes in flow rate and hydrochemistry have been widely used in the study of karst water systems. In particular, the study of changes in flow rate and hydrochemical composition during rainstorm events is considered an important way to understand the characteristics of karst water systems and an important method for inferring the source and composition of groundwater. Taking into account burial conditions and aquifer type, karst water can be divided into surface karst water (EW), saturated zone fissure water (FW), and surface water (PW). Karst underground river water is a mixture of EW, FW, and PW. Due to differences in burial conditions and water cycle conditions, EW, FW, and PW exhibit differences in their physical, chemical, and isotopic properties. This invention considers EW, FW, and PW as three endmembers and utilizes their physical, chemical, and isotopic differences. By employing endmember mixing analysis (EMMA) and MIX mixing ratio calculations, the flow process curve of EW is separated from the underground river flow process curve. Using the separated EW flow process curve, the total discharge of EW during rainfall events can be calculated. If the rainfall amount and intensity are sufficiently large, it can be assumed that the rainfall completely replaces EW from the surface karst zone, and the total discharge of EW is equal to the storage capacity of the surface karst water.
[0020] This invention utilizes the hydrodynamic characteristics of rapid discharge from surface karst zones into underground river channels after heavy rainfall. Based on the physicochemical differences between surface karst water and other water sources, it decouples the surface karst water flow curve from the underground river flow curve after heavy rainfall, indirectly obtaining the surface karst water storage capacity. Compared to existing methods for estimating surface karst water storage capacity, this method significantly reduces workload, saves substantial costs, and offers high accuracy. Furthermore, this method can be used for long-term monitoring of surface karst water storage capacity to understand its seasonal and interannual variations. Attached Figure Description
[0021] Figure 1 Monitoring curves for underground river outlet flow, physical and chemical indicators;
[0022] Figure 2 The proportions of EW, PW, and FW in each sample at the underground river outlet;
[0023] Figure 3 Flow curve of surface karst water after heavy rain;
[0024] Figure 4 Total discharge of surface karst water after heavy rain. Detailed Implementation
[0025] The following are specific embodiments of the present invention, which are described in conjunction with the accompanying drawings to further illustrate the technical solutions of the present invention. However, the present invention is not limited to these embodiments.
[0026] Step 1: Test and analyze the physical and chemical properties of EW, FW, and PW.
[0027] EW, PW, and FW were sampled and analyzed in a karst underground river basin. EW was characterized by surface karst spring water, PW by post-rainfall surface runoff, and FW by hydrogeological borehole water or large karst spring water that penetrates the karst saturation zone. The typical physical, chemical, and isotopic characteristics that distinguish EW water from PW and FW are mainly reflected in electrical conductivity (EC) and HCO3-. - Concentration, Sr 2+ Concentration, Ca 2+ Concentration and Mg 2+ The specific relationship between the concentrations of these five components is as follows: c[HCO3] - ] EW >c[HCO3] - ] FW >c[HCO3] - ] PW ,c[Sr 2+ ] FW >c[Sr 2+ ] EW >c[Sr 2+ ] PW ,c[Ca2+ ] PW <c[Ca 2+ ] EW and c[Ca 2+ ] FW ,c[Mg 2+ ] PW <c[Mg 2+ ] EW and c[Mg 2+ ] FW EC PW <EC EW and EC FW (c[HCO3 - ]、c[Sr 2+ ]、c[Ca 2+ ]、c[Mg 2+ [ ] represents the corresponding ion concentration, EC represents conductivity, and the subscripts correspond to the three end members.
[0028] Step 2: Continuous monitoring of flow, physical, chemical, and isotopic parameters at the underground river outlet during heavy rainfall periods.
[0029] Continuous monitoring of flow, physical, chemical, and isotopic parameters at the outlet of an underground river was conducted during a selected rainstorm event. The monitoring interval was 1–2 hours, with a total monitoring period of 15 days. Continuous monitoring curves for flow, physical, and chemical components were obtained, such as… Figure 1 As shown.
[0030] Step 3: Determine the final physical, chemical, and isotopic components used in the mixing ratio calculation.
[0031] Five hydrochemical indicators from the collected underground river outlet samples were standardized to form a dataset matrix. The selected indicators were HCO3-, HCO3-, and HCO3-. - Ca 2+ Mg 2+ and Sr 2+ Concentration and electrical conductivity (EC). EMMA analysis of the standardized dataset matrix showed that samples from both storm events retained two principal components with eigenvalues greater than 1, explaining over 90% of the differences. These two principal components could be used for the mixed calculation of the three endmembers. Therefore, HCO3- was finally determined. - Ca 2+ Mg 2+ and Sr 2+ Five components, including concentration and conductivity (EC), are used to calculate the mixing ratio.
[0032] Table 1. Principal component eigenvalues and explained variances identified by EMMA
[0033] principal component Eigenvalues Explained variance (%) Grand total(%) 1 3.196 63.9 63.9 2 1.366 27.3 91.2 3 0.346 6.9 98.2 4 0.076 1.5 99.7 5 0.016 0.3 100.0
[0034] Step 4: Recalculate the values of each component of the endmembers.
[0035] Most methods for calculating mixing ratios rely on the assumption that endmember concentrations are fully known, a situation that is rare. Endmember samples are often difficult to collect, and their concentrations vary temporally and spatially. Therefore, we use the MIX code to calculate mixing ratios and endmember concentrations. The MIX code, proposed by Carrera et al., is a maximum likelihood method for estimating mixing ratios. A key feature of this code is that it acknowledges that the uncertainty in endmember components may be greater than the uncertainty in the actual sample. Essentially, it assumes that the measured concentration C of the species in sample j is... mij It is composed of n e The conservatism of endmembers, coupled with a measurement error, contributes to this, and the measured concentration C of the endmembers is... mei It also includes errors. That is to say:
[0036]
[0037] C mei =C ei +ε ei i = 1, n s e = 1, n m (2)
[0038] Where εmij and εei are the measurement errors of sample concentration and endmember concentration, respectively, and the mixing ratio λej and endmember concentration Cei are obtained by minimizing the objective function:
[0039]
[0040] In the formula, σij and σei are the standard deviations of εmij and εei, respectively. By adjusting the variance to expand or limit the concentration calculation, eigenvector projection allows us to estimate the range of variance.
[0041] For the underground river outlet sample and the initial endmember, we assigned a standard deviation that is n times the standard deviation of component i across all sample measurements, i.e., σ ij =σ si , σ ei =nσ si (here σ) si Let n be the standard deviation of component i across all samples, and n be the number of iterations. The selection of recalculated results is based on the relationships between components determined by the three endmembers established in the first step, i.e., c[HCO3] - ] EW >c[HCO3] - ] FW >c[HCO3] - ] PW ,c[Sr 2+ ] FW>c[Sr 2+ ] EW >c[Sr 2+ ] PW ,c[Ca 2+ ] PW <c[Ca 2+ ] EW and c[Ca 2+ ] FW ,c[Mg 2+ ] PW <c[Mg 2+ ] EW and c[Mg 2+ ] FW EC PW <EC EW and EC FW (c[HCO3-], c[Sr) 2+ ]、c[Ca 2+ ]、c[Mg 2+ [These represent the corresponding ion concentrations, EC represents conductivity, and the subscripts correspond to the three endmembers]. Based on this, the average of the linear regression determination coefficients between the observed and calculated concentrations is calculated. Pick The maximum calculated value. Finally, the recalculated values of each component of the three endmembers are obtained.
[0042] Step 5: Calculate the mixing ratio of the three endmembers in each sample from the underground river outlet.
[0043] The recalculated values of each component of each endmember and the underground river outlet monitoring data of each component of each endmember are input into the MIX code to calculate the mixing ratio of the three endmembers in each sample at the underground river outlet, such as... Figure 2 As shown.
[0044] Step 6: Obtain the surface karst water flow process curve
[0045] Multiplying the proportion of surface karst water in each sample of underground river water by the corresponding underground river flow rate at that time yields the flow rate curve of the surface karst water, such as... Figure 3 As shown.
[0046] Step 7: Calculate the surface karst water reserves
[0047] By integrating the EW flow process curve, the total discharge of EW during the rainstorm event is obtained. This total discharge is approximately equal to the surface karst water storage. Figure 4 As shown in the medium gray area.
[0048] For any points not covered above, existing technologies shall apply.
[0049] Although specific embodiments of the present invention have been described in detail by way of examples, those skilled in the art should understand that the above examples are for illustrative purposes only and are not intended to limit the scope of the invention. Those skilled in the art can make various modifications or additions to the described specific embodiments or use similar methods to replace them, without departing from the direction of the invention or exceeding the scope defined by the appended claims. Those skilled in the art should understand that any modifications, equivalent substitutions, improvements, etc., made to the above embodiments based on the technical essence of the present invention should be included within the protection scope of the present invention.
Claims
1. A method for calculating the surface karst water storage in a karst underground river basin, characterized in that: Includes the following steps: S1: Collect surface karst water (EW), saturated zone fissure water (FW), and surface water (PW) samples, test their physical, chemical, and isotopic indicators, analyze and determine their representative components, and establish the relationship between the representative components of the three samples. S2: Select periods of heavy rainfall to continuously monitor the flow rate, physical, chemical, and isotopic indicators of the underground river outlet; S3: Perform EMMA analysis on the monitoring data of the underground river outlet in step 2 to determine the final physical, chemical and isotopic components used in the mixing ratio calculation; S4: Treat EW, FW and PW as three endmembers, recalculate the final component values of the components determined in step S3 using MIX iteration, and use the representative component relationship of the three types of water established in step S1 as a constraint. S5: Calculate the mixing ratio of the three endmembers in each sample at the underground river outlet using the component values calculated in step S4 and the monitoring data of the final components determined in step S3. S6: Reconstruct the mixing ratio calculated in step S5 into the flow data of the three terminal elements to obtain the flow process curve of each terminal element; S7: By integrating the EW flow process curve obtained in step S6, the total discharge of EW during the rainstorm event is obtained, which is equal to the surface karst water storage.
2. The method for calculating the surface karst water reserves in a karst underground river basin as described in claim 1, characterized in that: The specific operation of step S3 is as follows: perform EMMA analysis on the monitoring data of the underground river outlet in step 2 to analyze whether two principal components can be extracted. If two principal components can be extracted, the monitored components are used as the final physical, chemical and isotopic components used to determine the mixing ratio calculation. If more than two principal components are extracted, all the analyzed components are screened, and the representative components of EW, FW and PW determined in step 1 are retained first. Then, EMMA analysis is performed until two principal components are extracted. The finally screened components are used as the final physical, chemical and isotopic components used to determine the mixing ratio calculation.
3. The method for calculating the surface karst water storage in a karst underground river basin as described in claim 1, characterized in that: In step S6, the proportion of each end-member in each sample of the underground river water is multiplied by the underground river flow rate at the corresponding time to obtain the flow process curve of each end-member.
4. The method for calculating the surface karst water storage in a karst underground river basin as described in claim 1, characterized in that: In step S1, the typical physical, chemical, and isotopic characteristics that distinguish EW from PW and FW are manifested in the representative components.
5. The method for calculating the surface karst water storage in a karst underground river basin as described in claim 1, characterized in that: In step S1, the magnitude relationship of the test values of representative components of the three samples is established.