Hydrogeochemical evolution and processes based on multivariable statistical and inverse simulation modeling: A coal mine in the Northern Coalfield, China

It is a worrying fact that drinking freshwater accounts for only 3% of the world's water resources, and the water demand for the development of industry and agriculture, and residents' lives cannot be satisfied fully (Hoekstra 2014; Holland et al. 2015; Bonsal et al. 2020; Stewart-Harawira 2020). Thus, the reasonable utilization and assessment of freshwater resources have raised concern in the past few decades (Ferng 2007; Vollmer et al. 2018; Nurtazin et al. 2020). In terms of the water in the mining areas, underground aquifers and surface water systems will be affected after the extraction of mineral resources. Once wastewater generated by mining is poured into the ground surface, the quality of the surface water and below will be polluted (Zhang et al. 2019; Ju et al. 2022). Besides, the water table of aquifers influenced by water-flowing fissures resulting from the excavation of tunnels and working faces is decreasing year by year (Arkoc et al. 2016; Sahoo & Khaoash 2020; Sun et al. 2021). As such, the current situation of water resource management in mining areas is particularly prominent (Prathap & Chakraborty 2019; Xu et al. 2020). Considering the importance of groundwater to human health and industrial production, the hydro-chemical characteristics and groundwater assessment of mining areas have raised more attention in the scientific community.

Researchers have proposed some practical methods and techniques in qualitative analysis and a quantitative evaluation for the hydrogeochemical processes and evolution characteristics of groundwater, such as multivariate statistical analysis, stable isotopes and tracer elements, laboratory tests, and geochemical modeling (Liu et al. 2019; Roy et al. 2020; Qiu et al. 2021; Wijewardhana et al. 2022). Besides, the hydrogeochemical characteristics and the spatial distribution of groundwater in Bangladesh were investigated using a geographic information system (GIS) with geochemical modeling. Areas with the availability of the best-quality drinking water were identified (Das et al. 2021, 2022). Principle component analysis and hydrogeochemical modeling were utilized to investigate the hydrogeochemical evolution in groundwater of the target area (Hu et al. 2023; Jiang et al. 2023). The chemical quality of groundwater in Bushehr, southwest of Iran, was assessed for its suitability for drinking purposes through analysis of geochemical processes and their relation to water quality (Sarikhani et al. 2015). To investigate the comparative performance of two cover systems to reclaim acid mine drainage in Quebec, Canada, coupling with hydrogeological and geochemical numerical simulations was selected to evaluate the reclaimed performance of acid-generating tailling areas (Pabst et al. 2018). The hydrochemistry and water–rock interaction of ground surface water at the Mandi Baha-Ud-Din district in Pakistan was analyzed through GIS and bivariate Gibbs, and it was confirmed that most of the area was affected by water–rock interaction. The water quality in most of the area met drinking standards (Abbas et al. 2021). Water–rock interaction and chemical compositions of deep Taiyuan Formation limestone aquifer were discussed in coal mines, and extensive groundwater drainage has a great influence on water quality, accelerating the water–rock interaction (Chen et al. 2023). Based on previous relevant studies, the complicated rock–water interaction processes between different strata and aquifers have been analyzed and revealed (Morán-Ramírez et al. 2016; Yousif & El-Aassar 2018). Due to the intensive and large-scale mining, some toxic and harmful wastewater and suspended material (rock dust, coal cinder, and mechanical oil) as well as other contaminants that influence the groundwater environment are produced and may flow into groundwater through runoff and infiltration.

Scholars have made achievements in the study of hydrogeochemical processes and water–rock interaction. However, more attention is paid to the Mesozoic aquifers or surface water bodies. As a kind of important water source, the research on the whole Cenozoic loose sandy aquifer water of the extra thickness loose layer is relatively few. In this study, hydrogeochemical processes and evolution among the three aquifers distributed in the huge thick loose strata of the Gubei coal mine, Huainan mining area, are assessed by using the multivariate statistical analysis, ArcGIS 10.4, and inverse geochemical numerical simulation method. It is intended to provide insights into water–rock interaction controlling hydrogeochemical processes by evaluating the general groundwater chemistry characteristics and establishing geochemical evolution inverse modeling.

More than 20 coal-bearing layers are developed in the study area, most of which are distributed in Permian strata overlaid by the extra-thick Cenozoic strata. Main aquifers include the Cenozoic loose sandy aquifers, the Permian coal strata fractured sandstone aquifers, the Carboniferous Taiyuan Formation fractured limestone aquifers, and the Ordovician fractured limestone aquifers according to the sedimentary period of the stratum. The Cenozoic loose sandy aquifers consisting of the upper aquifer, middle aquifer, and lower aquifer are the main discussed targets in the present study.

Water samples were obtained from ground supply water wells, surface pumping wells, and auxiliary shaft and ventilation shaft including 8 samples from the upper aquifer, 20 samples from the middle aquifer, and 13 samples from the lower aquifer. The sampling location is shown in Figure 1(c). Before sampling, the collectors were preprocessed to prevent the original water quality from being polluted during sampling. Namely, collectors were treated with the diluted HCI solution with a concentration of about 2.5% to remove impurities and washed several times with pure water, and then they are dried and tightly sealed with clean caps. Each time while onsite sampling, we labeled the surface of the sample collectors in time. The information contents of the label mainly included collection location, depth, time, and sample number. Generally, some key physiochemical parameters such as temperature, pH value, electric conductivity, and dissolved inorganic carbon can be measured during the sampling (Liu et al. 2019). To prevent changes in water quality at higher outdoor temperatures, all the samples were stored in the refrigerator at the laboratory.

After that, the water samples were tested in the laboratory by using inductively coupled plasma-optical emission spectrometry for cations including calcium, magnesium, sodium, and potassium, and ion chromatography for anions including chloride and sulfate, whereas bicarbonate was analyzed through the titration method. Besides, the water quality test was required to be completed within 24 h after the water samples were collected. The data relating to the main ions have been acquired to further be analyzed. The statistical description of the test results of three groups water samples is shown in Table 1.

Table 1

Descriptive statistics for three aquifers water samples

Item .	Upper aquifer (n = 8) .				Middle aquifer (n = 20) .				Lower aquifer (n = 13) .
	Min. .	Max. .	Mean .	SD .	Min. .	Max. .	Mean .	SD .	Min. .	Max. .	Mean .	SD .
Ca2+	24.05	65.73	41.92	16.91	16.43	151.44	51.96	30.02	16.04	61.86	42.25	12.48
Mg2+	9.72	24.31	16.06	5.52	2.52	33.36	18.66	9.49	12.40	40.31	23.11	7.54
K+ + Na+	42.52	296.13	178.86	86.96	335.42	2,170.26	765.43	369.93	27.82	1,040.38	714.99	253.50
	0.00	0.00	0.00	0.00	0.00	73.81	12.42	19.26	0.00	18.45	5.95	6.88
	323.41	506.47	446.06	61.14	56.29	443.30	212.52	89.92	176.10	347.20	282.60	52.38
Cl−	35.00	177.25	72.29	47.77	3,35.71	2,987.73	912.15	547.48	13.83	1,047.55	797.39	298.32
	19.21	164.70	63.09	46.69	30.35	979.81	366.93	218.25	2.68	795.57	377.51	222.69
PH	7.90	8.45	8.23	0.22	8.18	10.60	8.70	0.61	7.62	8.70	8.29	0.32
TDS	354.00	831.00	597.38	155.69	924.00	6,017.00	2,288.95	1,074.73	290.00	3,112.00	2,122.92	698.84

Item .	Upper aquifer (n = 8) .				Middle aquifer (n = 20) .				Lower aquifer (n = 13) .
	Min. .	Max. .	Mean .	SD .	Min. .	Max. .	Mean .	SD .	Min. .	Max. .	Mean .	SD .
Ca2+	24.05	65.73	41.92	16.91	16.43	151.44	51.96	30.02	16.04	61.86	42.25	12.48
Mg2+	9.72	24.31	16.06	5.52	2.52	33.36	18.66	9.49	12.40	40.31	23.11	7.54
K+ + Na+	42.52	296.13	178.86	86.96	335.42	2,170.26	765.43	369.93	27.82	1,040.38	714.99	253.50
	0.00	0.00	0.00	0.00	0.00	73.81	12.42	19.26	0.00	18.45	5.95	6.88
	323.41	506.47	446.06	61.14	56.29	443.30	212.52	89.92	176.10	347.20	282.60	52.38
Cl−	35.00	177.25	72.29	47.77	3,35.71	2,987.73	912.15	547.48	13.83	1,047.55	797.39	298.32
	19.21	164.70	63.09	46.69	30.35	979.81	366.93	218.25	2.68	795.57	377.51	222.69
PH	7.90	8.45	8.23	0.22	8.18	10.60	8.70	0.61	7.62	8.70	8.29	0.32
TDS	354.00	831.00	597.38	155.69	924.00	6,017.00	2,288.95	1,074.73	290.00	3,112.00	2,122.92	698.84

Note: The units of ions concentration and TDS-mg/l; SD, standard deviation.

View Large

In this study, we conducted a preliminary analysis using the hydro-chemical type analysis method on the basic hydro-chemical data from the laboratory test. It is commonly accepted that a Piper trilinear diagram is applied to identify the hydrogeochemical characteristics (Piper 1944; Hill 1940). The ratio of major cations and anions can be visualized on the graph, where the diamond-shaped and triangle-shaped areas are classified into several subareas with individual hydro-chemical characteristics. Besides, we used a GIS with a powerful spatial interpolation function to investigate the distribution characteristics of the total dissolved solids (TDS) for three aquifers. Multivariate statistical approaches including correlation analysis and cluster analysis were taken to analyze the chemical characteristics of groundwater and its formation mechanism. The correlation between the main ions can reflect the groundwater chemical evolution processes and water–rock interaction. Another multivariate statistical method-cluster analysis is widely used in various fields, such as geochemistry and bioengineering. The process of categorizing the distance proximity of conventional ionic components allows for the assessment of water chemistry.

Finally, we adopted the reverse hydrogeochemical simulation to reproduce the chemical evolution characteristics of groundwater and water–rock interaction, and to further verify the accuracy of the conclusions obtained by multivariate statistical approaches. PHREEQC software used for geochemical calculation and inverse simulation is a powerful and widely applied program in modeling, and hydrogeochemical processes have been confirmed in the previous study (Triantafyllidis & Psarraki 2020; Islam & Mostafa 2022; Mosai et al. 2022). It is based on the investigated runoff condition, namely, the determination of recharge area, discharge area, and ‘possible mineral phases’, and the study of the differences in water chemical composition between the two areas, to speculate the possible mineral dissolution and precipitation along the specific path (Zhang et al. 2020, 2022; Jiang et al. 2022a; Lu et al. 2023).

The results provide confirmations to the conceptual sight that two flow compartments exist in the study area, and the hydro-chemical boundary probably lies between the upper aquifer (Quaternary aquifer) and middle aquifer (Neogene aquifer), indicating hydraulic connectivity is poor. Relative to the middle and lower aquifers, the shallow upper aquifer is influenced by natural climate and human activities, which result in a lower extent of leaching and concentrations of major ions in groundwater.

In the upper aquifer, Ca²⁺ and Mg²⁺, K⁺ + Na⁺, and /Cl⁻ are characterized by a positive correlation. Ca²⁺, K⁺ + Na⁺, ⁠, and Ca²⁺/Mg²⁺ have a negative correlation. Ca²⁺ and ⁠, K⁺ + Na⁺, and are characterized by the positive correlation in the middle aquifer. Mg²⁺ and ⁠, K⁺ + Na⁺, and Cl⁻ are characterized by a positive correlation in the lower aquifer. The correlation of ions implied complex chemical processes involving cation exchange and adsorption, dissolution and filtration of carbonate minerals and halite minerals, and desulfurization have occurred in groundwater.

Cluster analysis is the statistical concept that group data objects according to the information found in the data describing the objects and their relationships (Selmane et al. 2022; Wu et al. 2022). The greater the similarity within the group and the greater the gap between the groups, the better the clustering effect. Cluster analysis consists of q-mode and r-mode cluster analysis, and the q-mode cluster analysis based on sample classification is applied in this study, whose main idea is to use distance to measure the similarity between samples. The advantages have the following three aspects: (1) it can comprehensively use the information of multiple variables to classify samples; (2) the classification results are intuitive, and the clustering pedigree diagram clearly shows the numerical classification results; and (3) the results obtained by the cluster analysis are more detailed, all-round, and rational than those obtained by traditional classification methods.

Eight models are generated along trace 1 (S2 → S5) of the Cenozoic upper aquifer, as shown in Table 3. CO₂ is in an equilibrium or emission state on the whole, while halite and gypsum minerals are dissolved in the models. Calcite and dolomite minerals remain in the dissolved or precipitation state. Meanwhile, the occurrence of the cation exchange process has been determined by analyzing the molar transfer of NaX, CaX₂, and MgX₂.

PHREEQC generates four models along trace 2 (S10 → S22) in the middle aquifer. The results are presented in Table 4. Calcite mineral shows the precipitation or equilibrium characteristic probably due to the dissolution of CO₂. Halite and gypsum minerals are dissolved in the four models. However, the precipitation of dolomite has been found in most cases. Similarly, the simulated path is also accompanied by the cation exchange process, resulting in a lower concentration of Mg²⁺.

Finally, trace 3 (S32 → S41) is introduced to illustrate the hydrogeochemical process of the Cenozoic lower aquifer. As shown in Table 5, a higher concentration of Na⁺ may be related to the dissolution of halite and NaX, and the cation exchange process. The higher concentration of TDS in the final position is due to the halite and gypsum minerals dissolved in five models.

Through comprehensive analysis of hydro-chemical data using statistical techniques such as correlation analysis and cluster analysis, as well as inverse geochemical modeling based on PHREEQC, the hydrogeochemical processes were elucidated. The following conclusions have been obtained:

• (1) The hydro-chemical type of the upper aquifer is dominated by the K + Na-HCO₃ type, showing a lower concentration of Cl⁻, while the other two aquifers have similar water quality types dominated by the K + Na-Cl type and the K + Na-Cl + SO₄ type. Water samples are divided into five groups by cluster analysis, and the types of water quality are further detailed. The concentrations of Na⁺ + K⁺ and Cl⁻ increase from the upper aquifer to the middle aquifer, demonstrating increasing salinity in groundwater due to the dissolution of evaporated minerals.

• (2) Through geochemistry calculation, the SI of anhydrite, gypsum, halite, and CO₂(g) of aquifers is less than zero, indicating these minerals and gas are unsaturated. While aragonite, calcite, and dolomite of the middle aquifer remain in the unsaturated–saturated state, these minerals in the upper and lower aquifer are in the saturated state.

• (3) After correlation analysis and water–rock interaction analysis, the cation exchange process accelerating the reduction of Ca²⁺ concentration and the increase of concentration widely occurred in three aquifers, and dissolution of calcite and dolomite occurred simultaneously in most cases, except for the middle aquifer.

• (4) Inverse modeling further elucidates the hydrogeochemical processes and water–rock interaction and verifies the preliminary summaries obtained through statistical methods. The dissolution of gypsum and halite minerals widely occurs in three aquifers. Calcite and dolomite minerals also remain in the dissolved state in most cases, except for the middle aquifer. Besides, the cation exchange process releasing Na⁺ and absorbing Ca²⁺/Mg²⁺ is ubiquitous along the traces in three aquifers.

The authors are grateful for the help from the Gubei coal mine with respect to the basic data provided. Also, thanks go to the editors and anonymous reviewers for your constructive comments on the manuscript which greatly improved this paper.

This research was financially supported by the Key Program of the National Natural Science Foundation of China (42372316), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX23_2760), the Fundamental Research Funds for the Central Universities (2023XSCX003), and the Graduate Innovation Program of China University of Mining and Technology (2023WLKXJ003).

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.