Analysis of the Erosion Law of Karst Groundwater Using Hydrogeochemical Theory in Liulin Spring Area, North China

Abstract

The comprehensive geological, hydrogeological and hydrogeochemical model of the Liulin karstic spring area in the eastern limb of the Ordos syncline was established by a combination of chemical thermodynamics, chemical kinetics and hydrogeology. The study area was divided into four zones based on the saturation indices of calcite, dolomite and gypsum, which were computed by the groundwater-chemical simulation software PHREEQC (a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations), with consideration of the geological and hydrogeological conditions and hydro-geochemical reactions. The weight and volume modulus of carbonate rocks and sulphate rocks in each zone were calculated by the method of correlation analysis to evaluate the dissolution law of karst groundwater. The results showed that in the zone I (saturation index of calcite βc ≤ 1) the dissolution of calcite was the major geochemical reaction, the weight modulus of calcite was higher than that of dolomite and gypsum, and the pore space generated by the dissolution of calcite was one order of magnitude larger than that of dolomite and gypsum. In zone II (saturation index of calcite βc > 1 to saturation index of dolomite βd ≤ 1) the corrosion moduli were all smaller than that in zone I, the solubility of dolomite and gypsum increased, and calcite reached saturation. The space occupied by the calcite sediment was less than that dissolved by dolomite and gypsum. In zone III (saturation index of dolomite βd > 1 to saturation index of gypsum βg ≤ 1), calcite and dolomite had reached saturation, accompanied by dedolomitization, and the amount of gypsum dissolution increased obviously. The conclusions indicate that the hydrogeochemical environment plays an important role in mineral dissolution.

1. Introduction

Carbonate rocks are widely distributed in China, and its acreage accounts for about a third of the land area. In northern China, the location of karst spring is not only an important resource for local water supply but is also a famous tourist site [1]. Hence, the government sector and local residents have been attaching importance to water quality problems for many years. Water quality is closely associated with hydrochemical compositions [2], and carbonate rock dissolution plays an important role in the variation of such hydrochemical compositions. Hence, the study of the law of carbonate rock dissolution is of significant importance for obtaining a more comprehensive understanding of the groundwater environment. Scholars at home and abroad usually relate to qualitative and quantitative analysis of karstification and karst dissolution when they analyze the hydrogeology, engineering geology, and mineral geology problems in karst spring area [3]. It is necessary for analysis of karst hydrogeology characteristics to make a quantitative analysis and calculation of karst dissolution. For a long time, problems about magnitude and distribution of carbonate rocks karst dissolution—especially in the field of karst hydrogeochemistry and groundwater mixing—have been discussed by many scholars [4,5,6,7,8].

Currently, the methods for determining carbonate rock dissolution quantities mainly include the following. Firstly, to measure the karstic dissolution quantities repeatedly on the surface of the karst with micro erosion measurement instruments in the research field and then to calculate the actual amount of dissolution quantities [9,10]. Mottershead [9] had measured the eroded quantities of weathered bedrock for seven years and computed the average decrease rate (0.625 mm/a). Secondly, to put the standard test pieces of carbonate rocks into the air, ground, and soil layer in the study filed or the indoor specific experimental facility [11], then to measure the weight of the test pieces after a period of time and calculate the dissolved quantities during this process based on the quantity difference of test pieces before and after the carbonate rocks dissolution. Plan [12] investigated the influences of altitude, subsoil vs. sub-aerial exposure, vegetation, karst morphology, soil humidity, sample lithology, and sample surface morphology to carbonate dissolution by measuring the mass difference of carbonate tablets exposed to natural dissolution for one year. Thirdly, to calculate the total dissolved quantities (X, viz. X = T.Q) of the catchment basin in the karst area by measuring the outlet flow (Q) of the spring or underground river mouth and the concentrations (T) of mineral that carried by water [13,14,15]. Utilizing the formula of Corbel—which was one of the methods—Li et al. [16] computed the dissolved quantity of karst and analyzed the relationship of several factors such as temperature, annual rainfall, lithological characteristics and geomorphologic conditions based on a set of data of hydrochemistry points and the corresponding runoff depth of these points in Hongshui river basin.

In the 1990s, Cao et al. [17] proposed a whole set of theories and methods to determine the dissolved quantity of karst. However, so far these methods had not been applied widely in practical problems. The dissolution law of underground karst water in the Liulin spring area are calculated and analyzed by Cao’s theory and method in this paper. It will provide some basis for the formation and evolution of karst groundwater characteristics.

The remainder of the paper is organized as follows. Section 2 describes the information about Liulin spring area, four zones which divided by saturation indices are discussed in Section 3. Section 4 presents the main results and discussions of dissolution law in three zones. The summary and conclusions are given in Section 5.

2. Liulin Spring Area

2.1. Physical Geography

Liulin spring is one of the famous karst springs in Shanxi Province, China, and expose dispersedly in Sanchuan river valley in the village of Xuejiawan—Zhaidong, 3 km from Liulin county (showed in Figure 1). The study area is west of the middle of Lvliang Mountain, there are high and middle mountains where bedrock is exposed, and low and middle mountains with loess hilly regions in the east and west of the research area respectively. The climate in the area is continental semi-arid with the characteristics of a windy spring, hot summer, intensive precipitation in fall and drought and cold winter. The average value of temperature and precipitation are 9.2 °C and 506 mm respectively. The rivers in the area belong to the Yellow River system. Beichuan, Dongchuan and Nanchuan Rivers collect as the Sanchuan River at the Jiaokou County, then flows into the Yellow River. The basin area of the Sanchuan River is 4161 km², and its average annual runoff is 194 million m³/a.

2.2. Geology and Hydrogeology

Regular Precambrian metamorphic rock outcroppings can be found in the northern, eastern, and southern parts of the spring area, while Cambrian-Ordovician carbonate strata outcroppings are found mostly in the south-eastern and central parts of the spring area. Carboniferous-Triassic clastic rocks are mainly outcropped in the western and central regions of the spring area. Tertiary and Quaternary rickle overlies various ages of the bedrock. From a geotectonic point of view, the spring area is situated in a uniclinal structure which tilts to the west with a 2- to 8-degree dip angle, within which secondary folds and fractures develop. The geologic structure dominates the directions of the groundwater flow.

The northern and eastern boundaries of spring area are composed of drainage divides in metamorphic rock mountains, whereas the southern and south-eastern boundaries consist of subsurface divides in exposed areas of carbonate rocks. In the western part of the study area, carbonate rocks incline toward the syncline core of the Ordos Basin and are covered by thick-bedded clastic rocks for hundreds of meters. There is no discharge path for groundwater. Hence the groundwater is considered stagnant in this region. This boundary is defined as the lines where the roof depth of carbonate rocks reaches 1000 m. Cambrian and Ordovician marine carbonate rocks constitute the main aquifer for the water supply, especially the middle series of the Ordovician, which consists of limestone, dolomite, and some of gypsum. The upper and the lower confining beds are composed of Carboniferous-Triassic clastic rocks and Precambrian metamorphic rocks, respectively.

The karst groundwater system is mainly recharged in exposed areas of carbonate rocks, as defined by the recharge area. Leakage of river water from the Sanchuan River and its tributaries in exposed areas of carbonate rocks serve as additional recharge sources. In the center of the spring area, carbonate rocks are covered by thick-bedded clastic rocks which are composed of an interbed of sandstone and mudstone, preventing the precipitation recharging the aquifer. This serves as the flow-through area, with the waters flowing from the north, east, and south parts of study area towards the Liulin spring. There are two major discharge ways in this groundwater system: spring discharge and well exploitation. Places where the springs emerge in the Sanchuan River valley constitute the dominant discharge regions of the groundwater system.

2.3. Sampling and Analysis

Sample collection and component analysis of groundwater is the basis of the following analysis, the precision of which may directly influence reliability and rationality of the results. The groundwater samples were collected along predominant groundwater flow paths. The testing of water chemical components including temperature, pH value, Eh value, HCO₃⁻, SO₄²⁻, Cl, Na⁺ + K⁺, Ca²⁺, Mg²⁺ and TDS (Total Dissolved Solids) was accomplished in the water environmental monitoring center of Shanxi province.

Temperature, pH value and ORP (Oxidation-Reduction Potential) value were measured in situ using potable REX PHB-4 pH meters and a 501 ORP pole. Water samples for cation and anion analysis were filtered through 0.45 μm membranes in the field and collected in 500 mL polyethylene bottles. HCO₃⁻ was measured using the Gran titration method in the sampling day. Samples for cation analysis were acidified using 1:1 nitric acid to pH < 2 in field. The cations (K⁺ + Na⁺, Ca²⁺, Mg²⁺) and anions (Cl⁻, SO₄²⁻) were determined by atomic absorption spectrometry (PERSEE TAS990F) and ion chromatography (CIC300), respectively. The uncertainties of measurements were all within 5%.

3. Model Building and Feature Analysis of Hydrogeochemistry

3.1. Hydrogeochemical Partition and the Construction of Its Model

The minerals which formed the karst spring aquifer system mainly are calcite, dolomite, gypsum, etc. By studying the nature and characteristics of these minerals’ dissolution and precipitation, the corrosion properties of carbonate rocks and sulphate rocks can be understood very well [18,19]. The hydrogeochemical reaction proceeded in the system of solid and liquid phases that consisted of different mineral compositions, and the chemical reactions in the hydrogeochemical reaction would change along the direction of groundwater flow [20].

The saturation index (β) can not only determine the direction of dissolution and precipitation of different minerals but is also an important index for dividing the spring area into different hydrogeological zones and different hydrogeochemical environments. In the study area, the features of geology, hydrogeological conditions, hydrogeochemical reactions in the spring area, and the saturation indices of water points in the research area—which are calculated by the water chemical equilibrium simulation software PHREEQC—are comprehensively considered for dividing the four zones shown in Figure 1 (viz. βc ≤ 1, βc > 1 ~ βd ≤ 1, βd > 1 ~ βg ≤ 1 and βg > 1, where, βc, βd and βg represent the saturation index of calcite, dolomite and gypsum respectively). The hydrogeological-hydrogeochemical model for Liulin spring, as shown in Figure 2, is set up based on these horizontal zones.

3.2. Geochemical Characteristics Analysis

3.2.1. βc ≤ 1 (Zone I)

This zone, where the limestone in the layer of Cambrian and Ordovician is exposed or half-exposed and the surface karst developed well, is mainly distributed in the east, north and south of the study area. The groundwater form in zone I is phreatic water with the characteristics of vertical seepage, atmospheric precipitation recharge for groundwater, and the total dissolved solids (TDS) equals 260 mg/L. Calcite, dolomite and gypsum are in the dissolved state. The phreatic water surface and the line where βc equals 1 are the upper bound and lower bound of zone I respectively. Since the upper bound (phreatic water surface) is in contact with the air in the soil directly, this zone is called the CO₂ open system. The karst water chemistry reactions, including the reaction of CO₂ dissolved in phreatic water and the reaction between H⁺ generated by CO₂ dissolution and the calcite and dolomite in carbonate rocks, nevertheless, have little effect on the process of gypsum dissolution. The main equations of chemical reaction models in zone I are as follows:

CO₂ + H₂O = H₂CO₃

H₂CO₃ = H⁺ + HCO₃⁻

HCO₃⁻ = H⁺ + CO₃²⁻

CaCO₃ + H⁺ = Ca²⁺ + HCO₃⁻

CaMg(CO₃)₂ + 2H⁺ = Ca²⁺ + Mg²⁺ + 2HCO₃⁻

CaSO₄·2H₂O = Ca²⁺ + SO₄²⁻ + 2H₂O

3.2.2. βc > 1 ~ βd ≤ 1 (Zone II)

This zone is distributed in the Sanchuan River and the valley of its main tributaries. The main aquifer is covered by quaternary system strata. The velocity of karst water circulation is relatively quick because of the shallow buried layer. In this zone, the direction of groundwater flow transformed from vertical to horizontal, the value of TDS is 334 mg/L and calcite is in the saturated state, while dolomite and gypsum are in the dissolved state. The line where βc equals 1 and the line where βd equals 1 are respectively the upper and lower bounds of zone II. The form of groundwater transforms from phreatic water to confined water which CO₂ cannot directly pass into, so the hydrogeochemical environment of this zone is a relatively closed system. The karst water chemistry reaction mainly consumes CO₂ that dissolved in zone I, including the precipitation of calcite and dissolution of dolomite and gypsum. The main equations of chemical reaction models of zone II are as follows:

H₂CO₃ = H⁺ + HCO₃⁻

HCO₃⁻ = H⁺ + CO₃²⁻

CaMg(CO₃)₂ + 2H⁺ = Ca²⁺ + Mg²⁺ + 2HCO₃⁻

CaSO₄·2H₂O = Ca²⁺ + SO₄²⁻ + 2H₂O

Ca²⁺ + CO₃²⁻ = CaCO₃

3.2.3. βd > 1 ~ βg ≤ 1 (Zone III)

The zone III, whose aquifer is covered by the Carboniferous, Permian and Triassic strata, is mainly distributed in the midwest and northwest of the Liulin spring area. It should be noted that the aquifer in zone III has the following characteristics: TDS equals 846 mg/L, there is little variation in the range of water level, as well as slower groundwater circulation, flow velocity and permeability. Mineral ions in carbonate rocks have all reached the supersaturated state then precipitated from the water, leaving sulphate minerals to continue dissolving. The upper and lower bound of zone III are respectively the line where βd equals 1 and the line where βg equals 1. In addition to the precipitation of carbonate rocks and dissolution of sulphate rocks, the karst chemistry reaction is also accompanied by the reactions of de-sulphidation and de-dolomization. The main equations of chemical reaction models in zone III are as follows:

CaSO₄·2H₂O = Ca²⁺ + SO₄²⁻ + 2H₂O

Ca²⁺ + CO₃²⁻ = CaCO₃

Ca²⁺ + Mg²⁺ + 2CO₃²⁻ = CaMg(CO₃)₂

3.2.4. βg > 1 (Zone IV)

This region, located in the west of the Liulin spring area, not only has a relative good sealing performance but also is a reduction zone. The permeability of the aquifer in this area becomes worse and the groundwater circulation becomes very slow; zone IV is a stagnant area, in which the mineral ions of carbonate rocks and sulphate rocks are in the supersaturated state. The main equations of the chemical reaction models in zone IV are as follows:

Ca²⁺ + CO₃²⁻ = CaCO₃

Ca²⁺ + Mg²⁺ + 2CO₃²⁻ = CaMg(CO₃)₂

Ca²⁺ + SO₄²⁻ = CaSO₄

4. Calculation and Discussion of Corrosion Modulus

The corrosion modulus which consists of the weight modulus (M) and volume modulus (N) in each karst zones mentioned above are calculated by the method of correlation analysis with combining the simplified model of geological-hydrological geochemistry and the theory of the chemical thermodynamics and chemical kinetics. The weight modulus is defined as the weight of minerals dissolved by a unit volume of water, while the volume modulus is defined as the void volume dissolved by a unit volume of groundwater flowing through a certain karst area, or is filled by certain ions’ precipitation. In the process of calculating for M and N, the dissolution content of some ions in certain zones is computed by deducting the average ionic content in the rainfall from the total ionic content in groundwater. The calculated results are shown in

Table 1. Mineral ion concentration and the dissolved quantity of each water sample point.

Zone	Water Sample	Mineral ion Concentration (mmol/L)				Mineral ion Dissolved Quantity (mmol/L)
		HCO3−	SO42−	Ca2+	Mg2+	HCO3−	SO42−	Ca2+	Mg2+
βc ≤ 1	Jiuliwan	3.20	0.20	1.15	0.61	2.79	0.09	0.99	0.52
	Liujiagou	4.35	0.18	1.58	0.61	3.94	0.07	1.42	0.52
	Sanjiaozhuang	4.20	0.18	1.75	0.63	3.79	0.07	1.59	0.54
	Shizhuang	3.90	0.38	1.58	0.63	3.49	0.27	1.42	0.54
	Chemingyu	4.65	0.35	1.83	0.73	4.24	0.24	1.67	0.64
βc > 1 ~ βd ≤ 1	Youpingfang	4.45	0.45	1.55	0.96	4.04	0.34	1.39	0.87
	Chejiawan	3.95	0.12	1.40	0.61	3.54	0.01	1.24	0.52
	Tianjiahui	3.85	0.30	1.63	0.73	3.44	0.19	1.47	0.64
	Nancun	4.15	0.15	1.58	0.58	3.74	0.04	1.42	0.49
	Liujiagou	4.00	0.25	1.25	0.71	3.59	0.14	1.09	0.62
	Shang’an	4.49	0.86	1.61	1.07	4.08	0.75	1.45	0.98
	Tuanshuitou	4.90	2.25	2.46	1.42	4.49	2.14	2.30	1.33
βd > 1 ~ βg ≤ 1	Mazewan	5.35	1.30	1.80	1.17	4.94	1.19	1.64	1.08
	Loumenhui	4.05	0.68	1.63	0.96	3.64	0.57	1.47	0.87
	Shangqinglong	4.10	1.00	1.83	1.09	3.69	0.89	1.67	1.00
	Zhaidong	4.35	1.40	1.91	1.12	3.94	1.29	1.75	1.03
	Yangjiagang	4.40	0.95	1.75	0.99	3.99	0.84	1.59	0.90
	Liujiagata	3.75	2.18	2.51	1.95	3.34	2.07	2.35	1.86
	Yangjiagang	3.60	6.48	5.36	2.41	3.19	6.37	5.20	2.32
	Zhangjiazhuang	4.10	1.55	2.23	1.60	3.69	1.44	2.07	1.51
	Luotuoju	4.35	4.33	4.81	2.05	3.94	4.22	4.65	1.96
	Jucaita	4.85	2.88	3.61	2.05	4.44	2.77	3.45	1.96
ion quantity in rainfall		0.41	0.11	0.16	0.09

. It is necessary to mention that three zones (I, II and III) are analyzed in this section, since the mineral dissolution of carbonate rocks and sulphate rocks mainly occurs in these three zones, and all of the ions in zone IV are in the supersaturated state without dissolution, which means that corrosion modulus cannot be calculated.

4.1. βc ≤ 1 (Zone I)

According to the antecedent model, three minerals (calcite, dolomite, gypsum) can be totally dissolved in zone I. Based on the material equivalent principle of chemical combination and decomposition, the value of Ca²⁺ computed by water quality analysis is supposed to be the total quantities of Ca²⁺ (T_Ca2+) dissolved from calcite, dolomite and gypsum. The amount of Ca²⁺ and SO₄²⁻ dissolved from gypsum is equivalent, so that this part of Ca²⁺ can be expressed as (T_Ca2+/SO42−). Based on the previous analysis, the content of Ca²⁺ dissolved from carbonate rocks, which has correlativity with the content of HCO₃⁻, can be computed by deducting T_Ca2+/SO42− from T_Ca2+. Assuming that the contents of gypsum, calcite and dolomite dissolved into groundwater are m₁, m₂ and m₃ respectively, the dissolution equations of three kinds of minerals are presented as follows:

m₁CaSO₄·2H₂O = m₁Ca²⁺ + m₁SO₄²⁻ + 2m₁H₂O

(1)

m₂CaCO₃ = m₂Ca²⁺ + m₂CO₃²⁻

(2)

m₃CaMg(CO₃)₂ = m₃Ca²⁺ + m₃Mg²⁺ + 2m₃CO₃²⁻

(3)

Combining Equations (1)–(3) and balancing the coefficient, the following equation is obtained:

m₁CaSO₄·2H₂O + m₂CaCO₃ + m₃CaMg(CO₃)₂ + (2m₃ + m₂)CO₂ + (2m₃ + m₂)H₂O
=(m₁ + m₂ + m₃)Ca²⁺ + m₁SO₄²⁻ + m₃Mg²⁺ + (2m₂ + 4m₃)HCO₃⁻ + 2m₁H₂O

(4)

Equation (4) reveals that the product dissolved by sulphate has no effect on the dissolution of carbonate rocks, so the content of ion dissolved from sulphate and carbonate rocks could be calculated independently. In this zone, little ion is accumulated in the runoff flow and the human and other natural factors have a smaller influence on the aquifer. Therefore, the hypothesis that the dissolved quantities of Ca²⁺ and SO₄²⁻ are identical is reasonable, namely:

y₁ = A₁X₁ + B₁

(5)

where, y₁ denotes the average content of SO₄²⁻ (mmol/L), X₁ represents the average content of T_Ca2+/SO42- (mmol/L), and A₁ and B₁ are a constant that usually equal 1 and 0, respectively.

Based on material equivalent principle of decomposition, a linear relationship exists between the quantities of HCO₃⁻, Ca²⁺ and Mg²⁺, so the correlation equations are established as follows:

y₂ = A₂X₂ + B₂

(6)

y₃ = A₃X₂ + B₃

(7)

where, y₂ and y₃ are the contents of Ca²⁺ dissolved from calcite and Mg²⁺ dissolved from dolomite (mmol/L), X₂ denotes the average content of HCO₃⁻ (mmol/L), A₂ and A₃ are the regression coefficients, B₂ and B₃ are constant.

The values of dissolved minerals are utilized to draw the relationship plots of Ca²⁺ with SO₄²⁻, Ca²⁺ (Ca²⁺ − SO₄²⁻) with HCO₃⁻ and Mg with HCO₃⁻ (shown in Figure 3). The following correlation equations can be observed in Figure 3.

y₂ = 0.4012X₂ − 0.1921

y₃ = 0.0603X₂ + 0.3328

To substitute A₁X₁ for m₁, A₂X₂ for m₂, A₃X₂ for m₃, Equation (4) changes to:

A₁X₁CaSO₄·2H₂O + A₂X₂CaCO₃ + A₃X₂CaMg(CO₃)₂ + (2A₃ + A₂)X₂CO₂ + (2A₃ + A₂)X₂H₂O + ∑C = A₁X₁Ca²⁺ + A₁X₁SO₄²⁻ + 2A₁X₁H₂O + (A₂ + A₃)X₂Ca²⁺ + A₃X₂Mg²⁺ + (2A₂ + 4A₃)X₂HCO₃⁻ + B₂Ca²⁺ + B₃Mg²⁺

(8)

where, ∑C ≤ B₂Ca²⁺ + B₃Mg²⁺, ∑C can be seen as the constant in Equations (6) and (7), and is also the product of carbonate rock dissolution and has nothing to do with the dissolution of sulphate rocks. The following comprehensive chemical equation is established in the form of Equation (4):

(B₂ − B₃)CaCO₃ + B₃CaMg(CO₃)₂ + (B₂ + 2B₃)CO₂ + (B₂ + 2B₃)H₂O = B₂Ca²⁺ + B₃Mg²⁺ + (2B₂ + 4B₃)HCO₃⁻

(9)

where, B₂Ca²⁺ represents the total content of Ca²⁺ dissolved from calcite and dolomite, B₃ represents the content of Ca²⁺ dissolved from the dolomite, and (B₂ − B₃) represents the content of Ca²⁺ dissolved from calcite.

The Equations (8) and (9) are combined to obtain the following equations:

A₁X₁CaSO₄·2H₂O = A₁X₁Ca²⁺ + A₁X₁SO₄²⁻ + 2A₁X₁H₂O

(10)

A₂X₂CaCO₃ + (B₂ − B₃)CaCO₃ = (A₂X₂ + B₂ − B₃)CaCO₃

(11)

A₃X₂CaMg(CO₃)₂ + B₃CaMg(CO₃)₂ = (A₃X₂ + B₃)CaMg(CO₃)₂

(12)

The Equations (10)–(12) are employed to compute the weight modulus of gypsum, dolomite and calcite (M_g, M_d, M_c):

M_g = 0.5A₁X₁CaSO₄·2H₂O = 86.09X₁ = 86.09 × 0.15 = 12.914(g/m³)

M_d = 0.25(A₃X₂ + B₃)CaMg(CO₃)₂ = 46.1(A₃X₂ + B₃) = 46.1 × (0.0603 × 3.65 + 0.3328) = 25.489(g/m³)

M_c = 0.5(A₂X₂ + B₂ − B₃)CaCO₃ = 50.045(A₂X₂ + B₂ − B₃) = 50.045 × (0.4012 × 3.65 − 0.1921 − 0.3328) = 47.016(g/m³)

The total weight modulus (M₁) is the sum of each mineral modulus. Namely:

M₁ = M_c + M_d + M_g = 12.914 + 47.016 + 25.489 = 85.418(g/m³)

The volume modulus (Ng, N_d, N_c) of gypsum, dolomite and calcite are computed as follows:

N_g = M_g/γ_g = 86.09X₁/2.3 × 10⁻⁶ = 3.743X₁ × 10⁻⁵ = 3.743 × 10⁻⁵ × 0.15 = 5.615 × 10⁻⁶

N_d = M_d/γ_d = 46.1/2.8 × (A₃X₂ + B₃) × 10⁻⁶ = 1.646(A₃X₂ + B₃) × 10⁻⁵

= 1.646 × (0.0603 × 3.65 + 0.3328) × 10⁻⁵ = 9.103 × 10⁻⁶

N_c = M_c/γ_c = 50.045/2.7 × (A₂X₂ + B₂ − B₃) × 10⁻⁶ = 1.845(A₂X₂ + B₂ − B₃) × 10⁻⁵

= 1.854 × (0.4012 × 3.65 − 0.1921 − 0.3328) × 10⁻⁵ = 1.741 × 10⁻⁵

Then the total volume modulus (N₁) is computed as:

N₁ = N_g + N_d + N_c = 5.615 × 10⁻⁶ + 9.103 × 10⁻⁶ + 1.741 × 10⁻⁵ = 3.213 × 10⁻⁵

In the exposed area of carbonate rocks where βc ≤ 1, since the atmospheric precipitation has directly contact with aquifer the velocity of water exchange becomes rapid, which intensify the extent of water dissolution of the karst . Calcite, dolomite and gypsum are in the dissolved state, where the dissolved quantities of gypsum in the unit volume of water is minimum; M_g equals 12.914 g/m³ and N_g equals 5.615 × 10⁻⁶, and the dissolved quantities of dolomite are such that M_d equals 25.489 g/m³ and N_d equals 9.103 × 10⁻⁶. The dissolved quantities of calcite in the unit volume of water is maximal with M_c equal to 47.016 g/m³ and N_c equals 1.741 × 10⁻⁵. The total weight modulus and volume modulus of the three minerals are 85.418 g/m³ and 3.213 × 10⁻⁵, respectively. Evidently, calcite dissolution, whose dissolved quantity is larger than that of dolomite and gypsum, is more likely given priority in this zone. The mineral dissolution lead to the new space increasing, and the analysis which is mentioned above indicates that the volume modulus of new space created by calcite dissolution is one order of magnitude larger than that of dolomite and gypsum dissolution.

4.2. βc >1 ~ βd ≤ 1 (Zone II)

When it comes to zone II, the dissolved quantities of Ca²⁺ and CO₃²⁻ are in the supersaturated state, so that the dissolution of calcite does not occur in this zone but the dissolution of dolomite and gypsum still exist. Based on the same principle as zone I, the comprehensive chemical equations (in zone II) about the dissolution of dolomite and gypsum are as follows:

m₄CaMg(CO₃)₂ = m₄Ca²⁺ + m₄Mg²⁺ + 2m₄CO₃²⁻

(13)

m₅CaSO₄·2H₂O = m₅Ca²⁺ + m₅SO₄²⁻ + 2m₅H₂O

(14)

m₄CaMg(CO₃)₂ + m₅CaSO₄·2H₂O = m₄CaCO₃ + (m₅Ca²⁺ + m₄CO₃²⁻) + m₄Mg²⁺ + m₅SO₄²⁻ + 2m₅H₂O

(15)

According to Equation (15), 2m₄CaCO₃ (i.e., m₄CaCO₃ + m₄Ca²⁺ + m₄CO₃²⁻) represents the precipitation of calcite when m₅Ca²⁺ ≥ m₄CO₃²⁻ and the extra quantity of (m₅ − m₄) Ca²⁺ saves in the karst water, while (m₄ + m₅) CaCO₃ (i.e., m₄CaCO₃ + m₅Ca²⁺ + m₅CO₃²⁻) represents the quantity of calcite precipitation when m₅Ca²⁺ > m₄CO₃²⁻ and the extra quantity of (m_{4 −} m₅) CO₃²⁻ saves in the karst water.

After deducting the content of Ca²⁺ dissolved from gypsum from total content in zone II, the rest is the sum content of Ca²⁺ which consists of the part inputted from zone I and the other part dissolved from dolomite in zone II. It is difficult to determine the quantity of CaCO₃ which formed by Ca²⁺ and CO₃²⁻, so the values of A₄ and B₄ are computed by the relationship between Mg²⁺ and HCO₃⁻ based on the equation as follows:

y₄ = A₄X₄ + B₄

(16)

where, y₄ and X₄ are the content of Mg²⁺ and HCO₃⁻ (mmol/L). A₄ and B₄ are the coefficient of regression and constant, respectively.

Both the content of Ca²⁺ and SO₄²⁻ dissolved from gypsum and accumulated from zone I are still few. Assuming that the dissolved quantities of Ca²⁺ and SO₄²⁻ are identical, the correlation equation of Ca²⁺ and SO₄²⁻ is presented as follows:

y₅ = A₅X₅ + B₅

(17)

where, y₅ and X₅ denote the content of SO₄²⁻ and Ca²⁺, which are equivalent (mmol/L). A₅ and B₅ are constants of 1 and 0, respectively.

According to the values of dissolved minerals, the relationship plots of Ca²⁺ with SO₄²⁻ and Mg with HCO₃⁻ are drawn (Figure 4). The following correlation equations can be observed from Figure 4.

y₄ = 0.7429X₄ − 2.0786

The following equation can be obtained by putting the Equations (16) and (17) into Equation (15):

A₄X₄CaMg(CO₃)₂ + A5X5CaSO₄·2H₂O + ∑C = A₄X₄CaCO₃↓ + A₄X₄CO₃²⁻ + A₅X₅Ca²⁺ + A₄X₄Mg²⁺ + A₅X₅SO₄²⁻ + 2ª₅X₅H₂O + B₄Ca²⁺ + B₄Mg²⁺

(18)

where ∑C ≤ B₄Ca²⁺ + B₄Mg²⁺ is the product of carbonate rock dissolution and has nothing to do with the dissolution of sulphate:

B₄CaMg(CO₃)₂ = B₄Ca²⁺ + B₄Mg²⁺ + B₄CO₃²⁻

(19)

Equations (18) and (19) are combined to obtain the following equations:

A₄X₄CaMg(CO₃)₂ + B₄CaMg(CO₃)₂ = (A₄X₄ + B₄)CaMg(CO₃)₂

(20)

A₅X₅CaSO₄·2H₂O = A₅X₅Ca²⁺ + A₅X₅SO₄²⁻ + 2A₅X₅H₂O

(21)

Equations (20) and (21) are used to compute M_g and M_d:

M_g = 0.5A₅X₅CaSO₄·2H₂O = 86.09X₅ = 86.09 × 0.25 = 21.523(g/m³)

M_d = 0.25(A₄X₄ + B₄)CaMg(CO₃)₂ = 46.1(A₄X₄ + B₄) = 46.1 × (0.7429 × 3.85 − 2.0786) = 36.030(g/m³)

According to the calculated results of mineral dissolution, the content of SO₄²⁻ is two times less than that of Mg²⁺, so M_c can be computed as:

M_c = −50.045 × (0.7429 × 3.85 + 0.25 − 2.0786) = −51.625(g/m³)

The total weight modulus (M₂) can be calculated as:

M₂ = M_g + M_d + M_c = 21.523 + 36.030 − 51.625 = 5.928(g/m³)

N_g, N_d, N_c are computed as follows:

N_g = 3.743 × 10⁻⁵X₅ = 3.743 × 10⁻⁵ × 0.25 = 9.358 × 10⁻⁶

N_d = 1.646(A4X4 + B4) × 10⁻⁵ = 1.646 × (0.7429 × 3.85 − 2.0786) × 10⁻⁵ = 1.287 × 10⁻⁵

N_c = −1.854(A₄X₄ + A₅X₅ + B₄) × 10⁻⁵ = −1.854 × (0.7429 × 3.85 + 0.25 − 2.0786) × 10⁻⁵ = −1.912 × 10⁻⁵

The total volume modulus (N₂) can be calculated as:

N₂ = N_g + N_d + N_c = 9.358 × 10⁻⁶ + 1.287 × 10⁻⁵ − 1.912 × 10⁻⁵ = 3.105 × 10⁻⁶

In zone II βc >1 > βd ≤ 1, the dissolved quantity of gypsum in the unit volume of water is still minimal, M_g equals 21.523 g/m³ and N_g equals 9.358 × 10⁻⁶. The dissolved quantity of dolomite is increased, which can be seen from M_d (36.030 g/m³) and N_d (1.278 × 10⁻⁵). The dissolution of calcite begins to be subside since it has already reached the supersaturated state and the values of M_c and N_c equal −51.625 g/m³ and −1.912 × 10⁻⁵ respectively. The total weight modulus and volume modulus of the three minerals are 5.928 g/m³ and 3.105 × 10⁻⁶. Evidently, dolomite dissolution, whose dissolved quantity is larger than that of gypsum, is more likely given priority in this zone while calcite starts to precipitate. Generally speaking, the new space created by the dissolution of dolomite and gypsum is larger than the space filled by the precipitation of calcite, which can be seen from the value of N₂ (3.105 × 10⁻⁶ > 0).

4.3. βd >1 ~ βg ≤ 1 (Zone III)

In this zone, the carbonate rocks cannot be dissolved by groundwater any more, and the chemical reactions mainly consist of gypsum dissolution, calcite precipitation and the reactions of desulphurization and dedolomitization, which can be presented as the following equations:

m₆CaSO₄·2H₂O = m₆Ca²⁺ + m₆SO₄²⁻ + 2m₆H₂O

(22)

m₇CaMg(CO₃)₂ + m₇Ca²⁺ = m₇Mg²⁺ + 2m₇CaCO₃↓

(23)

The solubility of gypsum is high, and dedolomitization needs to consider some conditions. In the Equation (23), the exchange quantity of Ca²⁺ (in Equation (23)) is far lower than the quantity of Ca²⁺ (in Equation (22)), that is to say, m₆ > m₇. Combining the Equation (22) and (23) then obtains the following equation:

m₆CaSO₄·2H₂O + m₇CaMg(CO₃)₂ = (m₆·− m₇)Ca²⁺ + m₆SO₄²⁻ + 2m₆H₂O + m₇Mg²⁺ + 2m₇CaCO₃

(24)

According to Equation (24), the relationship plot is drawn (Figure 5), and the following correlation equations can be observed from Figure 5.

y₆ = 0.7253X₆ + 1.0131

y₇ = 0.2590X₇ + 0.8869

The following contents of different materials can be obtained according to Equation (24):

(A₆X₆ + B₆)CaSO₄·2H₂O

(25)

(A₇X₇ + B₇)CaMg(CO₃)₂

(26)

(A₆X₆ + B₇)CaCO₃

(27)

In this zone, the dedolomitization, a physical chemistry reaction, is defined as the exchange between dolomite and Ca²⁺ in groundwater. The two kinds of chemistry reaction of the dissolution of gypsum and precipitation of calcite are main considerations for the rock reconstruction in zone III. Equations (25) and (27) are used to compute M_g and M_c:

M_g = 86.09(A₆X₆ + B₆) = 86.09 × (0.7253 × 2.16 + 1.0131) = 222.091(g/m³)

M_c = −100.089(A₆X₆ + B₇) = −100.089 × (0.7253 × 2.16 + 0.8869) = −245.573(g/m³)

The total weight modulus (M₃) can be calculated as:

M₃ = M_g + M_c = 222.091 − 245.573 = −23.482(g/m³)

N_g and N_c are computed based on the weight modulus:

N_g = 3.743(A₆X₆ + B₆) × 10⁻⁵ = 3.743 × (0.7253 × 2.16 + 1.0131) × 10⁻⁵ = 9.656 × 10⁻⁵

N_c = −3.707(A₆X₆ + B₇) × 10⁻⁵ = −3.707 × (0.7253 × 2.16 + 0.8869) × 10⁻⁵ = −9.095 × 10⁻⁵

The total volume modulus (N₃) can be calculated as:

N₃ = N_g + N_c = 9.656 × 10⁻⁵ − 9.095 × 10⁻⁵ = 5.61 × 10⁻⁶

In the zone that βc > 1 ~ βd ≤ 1, the dissolved quantity of gypsum in the unit volume water increases dramatically, with M_g equals to 222.091 g/m³ and N_g equal to 9.656 × 10⁻⁵. Calcite (M_c and N_c are −245.573 g/m³ and −9.095 × 10⁻⁵, respectively) is still in the supersaturated state and is accompanied by dedolomitization. The total weight modulus in this zone is −23.482 g/m³, so the quantity of calcite precipitation is greater than that of gypsum dissolution. The value of the total volume modulus (5.61 × 10⁻⁶) indicates that the groundwater influence on the aquifer will still form little amounts of new pore space. It is of significant importance to get a more comprehensive understanding of interaction between water and rocks by combining these two kinds of moduli (the weight modulus and volume modulus) for the evaluation of the dissolution quantities in zone III.

In conclusion, there are some difference from zone I to zone II, then to zone III. Firstly, the value of M_c reduced from 47.016 g/m³ to −51.625 g/m³, then to −245.573 g/m³, and the value of N_c reduced from 1.7 × 10⁻⁵ to −1.912 × 10⁻⁵, then to −9.095 × 10⁻⁵, which indicates that the corrosion modulus of calcite decreased and the chemical reaction transformed dissolution into precipitation, and the ability of calcite precipitation continued strengthening from zone II to III. Secondly, the value of M_d increased from 25.489 g/m³ to 36.030 g/m³, and N_d increased from 9.103 × 10⁻⁶ to 1.287 × 10⁻⁵, which reveals that the corrosion modulus of dolomite was increasing gradually, and the ability of dolomite dissolution continued strengthening from zone I to zone II, then dolomite reached the state of supersaturation with weak dedolomitization. Thirdly, the value of M_g increased from 12.914 g/m³ to 21.523 g/m³, then to 222.091 g/m³ and N_g increased from 5.615 × 10⁻⁶ to 9.358 × 10⁻⁶, then to 9.656 × 10⁻⁵, which shows that the corrosion modulus of gypsum increased gradually, however, its dissolution ability is weaker than that of carbonate rocks. Fourthly, the total weight modulus reduced from 85.418 g/m³ to 5.928 g/m³, and finally to −23.482 g/m³, and the total volume modulus reduced from 3.213 × 10⁻⁵ to 3.105 × 10⁻⁶, then increased to 5.61 × 10⁻⁶, which manifests that the quantities of precipitation of carbonate rocks minerals was greater than that of sulphate rocks, but a tiny amount of new pore space would still be produced by the interaction between the groundwater and karst aquifer.

5. Conclusions

This paper considers calcite, dolomite and gypsum as representative minerals of karst rocks in Liulin spring and combines the geology, hydrogeological conditions and hydrogeochemical reactions within the different zones which were divided by saturation indices of calcite, dolomite and gypsum to develop the hydrogeological-hydrogeochemical model of Liulin spring. The two prime issues can be summarized:

(1)

In zone I where βc ≤ 1, three kinds of minerals (viz. calcite, dolomite and gypsum) are in the dissolved state. The corrosion modulus of calcite is maximum and the second is that of dolomite, and the value of corrosion modulus of gypsum is minimum; in zone II where βc ~ 1 > βd ≤ 1, the value of corrosion modulus of calcite is negative, which indicates that the reaction of calcite transforms dissolution to precipitation. The corrosion modulus of dolomite and gypsum increased, and the corrosion reaction was dominated by carbonate rocks dissolution; in the third zone where βd > 1 ~ βg ≤ 1, calcite and dolomite achieved a supersaturated state and the corrosion modulus of gypsum increased remarkably.

(2)

From zone I to zone III, the corrosion modulus of calcite decreased gradually along the direction of groundwater flow, however, the corrosion modulus of dolomite increased at first then decreased, and the corrosion modulus of gypsum was augmented with the increase in aquifer depth. The dissolution law of karst groundwater in the Liulin spring area shows that the hydrogeochemical environment plays an important role in mineral corrosion. However, there is less work in the hydrodynamic analysis of the spring system, and more efforts in this regard should be made in the future.