Characteristics of the water cycle of Jinan karst spring in northern China

Karst is widely distributed all over the world. Rich water resources could guarantee the development of our society. However, under the influence of human activities, the vulnerability of karst ecosystem is prominent. The sharp decrease of karst spring flow and the deterioration of water quality have become serious problems all over the world (Chu et al. 2017; Sun et al. 2017; Liang et al. 2018). The karst area in northern China is 992,000 km², which is distributed in 13 provinces (Liang & Zhao 2018) (districts and cities). In recent years, the karst ecological environment has also been plagued by many problems such as pollution and destruction (Xing et al. 2010; Yang et al. 2016). Over the years, great investment has been made in the restoration of famous springs in Shanxi (Yuan 2010; Wang et al. 2020), but the restoration effect is not significant. Taking Jinan as an example, to restore the continuous gushing of famous springs, a number of spring conservation measures have been taken successively, such as mining outside and supplementing inside, diverting the Yellow River from the reservoir to protect the spring, shutting down mining, and a surface water to groundwater artificial recharge project (Xing 2006; Hu et al. 2011; Xing et al. 2018). However, the threat of spring water cut-off in the dry season still cannot be eliminated. The fundamental reason is the lack of in-depth understanding of the cyclic evolution mechanism of different levels of flow system, and there is a certain blindness in formulating spring protection countermeasures.

The regional groundwater flow theory (Toth 1963) proposed by the hydrogeologist Toth describes the law of regional groundwater circulation better, and believes that the groundwater circulation does not completely depend on the geological structure, but is also controlled by the topographic potential energy (Liang et al. 2012, 2022). The local groundwater flow system has short runoff, shallow depth, large circulation depth, and long retention time. Therefore, the hierarchical groundwater flow system has different hydrodynamic, physical and chemical characteristics and manifestations. Taking the groundwater flow system of Ordos Cretaceous Basin as an example, Jiang (2011) deduced the analytical solutions of water head and flow function on the two-dimensional section of the basin, and directly simulated the age distribution in the groundwater flow system by using the groundwater age control equation. Wang (2015) proposed three methods of geophysical prospecting, conductivity, and residence time to divide the level of flow system, established the ‘late peak’ method based on the retention time distribution of groundwater in the basin, and verified the effectiveness of this method, which has valuable reference significance for the study of the karst spring flow system. Wang (2016) revealed the relationship between karst hot water and regional groundwater flow in carbonate heat reservoir and the circulation characteristics of groundwater at different scales. Through the retention time and age distribution of groundwater, they identified the nested structure of multistage groundwater flow system in the basin, and reexamined the circulation and evolution process of groundwater from the perspective of time domain. Xiao et al. (2022) used numerical simulation and an indoor sandbox experiment to explore the interaction between seepage field and temperature field under different groundwater flow system development modes, and confirmed that temperature difference is an important driving force of a groundwater flow system in geothermal-rich areas. Gao et al. (2022) demonstrated the application of groundwater layered exploration technology in the study of a groundwater flow system, which laid a foundation for the determination of spatial movement direction of groundwater flow and the delimitation of groundwater flow system. In terms of practical application, Wang et al. (2021) used isotopes and hydrogeochemistry to reveal the groundwater flow pattern and recharge source in the aquifer system of arid inland basins in north-west China, and inferred the hydrogeochemical process in the process of groundwater circulation. Selene et al. (2020) studied the water rock interaction and mixing process of karst water system under the mining conditions of Mexico City, divided the groundwater flow system involved in mining into four parts, and confirmed the existence of water exchange among groundwater flow systems through reverse hydrogeochemical simulation. Hu et al. (2015) used hydrogen and oxygen stable isotopes to study the seasonal recharge of surface karst springs and the average residence time of surface karst water. Agoubi (2021) took the geothermal water resources in Tunisia as the research object, determined the heat source and flow pattern in different levels of groundwater flow system, and analyzed that the groundwater temperature and ion source in nonvolcanic areas depend on the retention time and circulation depth of groundwater.

The existing research results show that the regional groundwater flow system theory provides a theoretical basis for the division of a groundwater flow system and the in-depth exploration of the dynamic law of karst water. Only by correctly understanding and reasonably dividing the levels of groundwater flow system can we fully reveal the integrity, hierarchical structure, correlation, and systematic differences of karst springs, further excavate the unique hydrogeological information hidden in the recharge and discharge of spring water in Jinan and the common characteristics of karst springs in northern China. Although predecessors have studied the karst water system in Jinan, the lack of classification of the level of spring water flow system is not conducive to accurately depicting the dynamic changes of spring water, nor can it achieve the purpose of accurate spring conservation.

Therefore, in this study, the groundwater in Jinan spring karst aquifer is characterized by fixed depth stratified sampling monitoring, isotope and hydrochemical data, so as to more accurately characterize the characteristics of different levels of ground-water flow system. The purpose of this study is to explore the water cycle evolution model of karst water system, and correctly identify the hierarchical karst groundwater flow system combined with multiple factors, which can not only enrich the theory of karst water system, but also provide a theoretical basis for the protection of spring water in Jinan.

A total of 71 sampling points were arranged in this study (Figure 1), including 3 recharge sampling points, 14 spring water sampling points, and 54 karst water sampling points of Cambrian and Ordovician limestone. Water samples were collected and stored in accordance with the operating procedures of the ‘Standard Test Method for drinking water’ (GB5750-2006). The conventional ion components are tested by hydrochemical titration and flame method. The test indexes include K⁺, Na⁺, Ca²⁺, Mg²⁺, ⁠, ⁠, Cl⁻, and Sr²⁺ are detected by inductively coupled plasma emission spectrometer. Isotopic ²H and ¹⁸O ratios are detected by Xi'an Guolian Quality Detection Technology Co., Ltd with Picarro L2130-1 water isotope analyzer. The detection basis is BKR-FF-01-2016 ‘Determination of hydrogen and oxygen stable isotopes in water laser cavity ring down spectroscopy’. Isotopic ³H samples were sent to the Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences for testing. Aqua TROLL 600 multiparameter water quality detector produced by In-Situ company was used to monitor the conductivity of spring and well water at fixed depth in real time in the field, and the conductivity resolution was 0.01 μS/cm, the monitoring frequency was 1–2 times per month, and the observation time was 2017a–2021a.

The highly heterogeneous and anisotropic development of karst aquifer media makes the law of groundwater migration extremely complex (Dong 2020), and there must be hydrodynamic dispersion (Liu et al. 2003). In this study, the inverse statistical model (Rose 1993) (ISM) is used to calculate the age of groundwater. Taking the dispersion parameters in the process of groundwater mixing as the standard deviation, the following algorithm can be used to calculate the average residence time range of groundwater level under different standard deviations.

is the tritium decay correction concentration in year i (⁠⁠); is the tritium input concentration in a given year i (⁠⁠); is the decay constant, generally taken as 0.0558; t is the input year.

The percentage weight of all mixed water is allocated to the average decay correction value during adjacent years, that is, each value obtained from Equation (2).

The total tritium concentration in the period with n-year data is the sum of each weighted average value in that period.

Hydrogen and oxygen isotopes are relatively stable and can be used as tracers to study the interaction between different water bodies (Sappa et al. 2012; Khalil et al. 2015). In this study, the IsoSource linear model (Phillips & Gregg 2003) is used to calculate the contribution rate of hydrogen and oxygen isotopes in different groundwater recharge sources to the concentration of hydrogen and oxygen isotopes in spring water.

represents the number of combinations, i represents the increment range, and s represents the number of sources. The weighted average value of each combination is compared with the isotope value actually measured by the mixture to obtain a feasible solution. In all feasible solutions, the frequency of each resource contribution percentage is analyzed to obtain the contribution rate of different recharge sources.

The calculation process is executed by the Visual Basic program of IsoSource. All feasible source combinations can be obtained by providing the source of spring water recharge, the hydrogen and oxygen isotope values of the samples of the four large spring groups, as well as the required source increment and mass balance tolerance, and the IsoSource program can be downloaded and used on the official website of the U.S. Environmental Protection Agency (Phillips & Gregg 2003).

According to the data of precipitation tritium concentration in Jinan spring area and formula (1) in this paper, the decay correction value of precipitation tritium concentration in Jinan spring area in each year is obtained (Table 1). Then calculate according to formulas (2) and (3). Based on the standard deviation of 5, 10, 15,…, 95, 100, the mean transit time of groundwater from 5a to 100a in Jinan spring domain is calculated (Table 2).

Table 2

Calculated tritium concentration corresponding to different average residence time and standard deviation in Jinan area

View Large

The test results indicated that the tritium concentrations of Baotu Spring, Heihu Spring, Tanxi Spring, and Wangfuchizi Spring were 8.5 TU, 6.9 TU, 8.6 TU, and 8.4 TU, respectively. Among them, the tritium concentration of Heihu Spring was significantly smaller than that of other springs. According to the water quality monitoring results (Table 3), the Cl⁻ concentration of the Heihu Spring group was greater than 70 mg/L, which was substantially larger than that of other spring groups by about 40 mg/L–60 mg/L. The ³H concentration value was similar to that of domestic water, which was obviously smaller than that of atmospheric precipitation (Ren et al. 2000). We considered the Heihu Spring as a spring greatly affected by human activities, thus, the tritium concentration of the Heihu Spring was used as the basis for evaluating the age of spring water. Based on the calculation results of the ISM model, different tritium concentrations had three age distribution ranges: less than 10a, 30a–40a, and 50a–70a. Combined with the piston flow model and full mixing model (Yang 2008), the age distribution range of the four large spring groups was estimated to be 6a–35a. Combined with the hydrogeological conditions and genetic model of the four large spring groups in Jinan, the circulation depth of Baotu Spring and Heihu Spring was assessed to be less than that of Tanxi Spring and Wangfuchizi Spring (Sun et al. 2021). Moreover, the retention time of groundwater in the local groundwater flow system participating in the shallow cycle is considerably smaller than that in the regional groundwater flow system participating in the deep cycle (Jiang 2011). Thus, the spring age of Baotu Spring and Heihu Spring should also be less than that of Tanxi Spring and Wangfuchizi Spring (Fang et al. 1989; Zhang et al. 2012). Therefore, the groundwater age of Baotu Spring and Heihu Spring was corrected to 6a–10a, and the groundwater age of Tanxi Spring and Wangfuchizi Spring was corrected to 30a–35a. The age of groundwater is the embodiment of the difference between groundwater circulation and replacement rate, suggesting that the recharge sources and circulation depth of the four large spring groups are different.

Table 3

Ion content of spring water in four large spring groups

Ion content (mg/L) .	Heihu Spring .		Baotu Spring .		Tanxi Spring .		Wangfuchizi .
	2019a .	2020a .	2019a .	2020a .	2019a .	2020a .	2019a .	2020a .
K+	1.18	0.85	1.2	1.67	2.27	2.42	–	1.02
Na+	34.02	30.60	24.59	36.50	25.33	28.5	–	17.50
Ca2+	130.26	137.00	120.84	111.00	106.81	111	–	114.00
Mg2+	23.09	22.20	22	21.70	21.15	22.1	–	20.70
Cl−	71.97	70.30	60.62	65.80	55.66	59.4	–	47.60
	119.74	119.00	104.08	116.00	100.14	105	–	86.60
	318.43	299.36	297.2	262.70	272.94	262.7	–	268.81

Ion content (mg/L) .	Heihu Spring .		Baotu Spring .		Tanxi Spring .		Wangfuchizi .
	2019a .	2020a .	2019a .	2020a .	2019a .	2020a .	2019a .	2020a .
K+	1.18	0.85	1.2	1.67	2.27	2.42	–	1.02
Na+	34.02	30.60	24.59	36.50	25.33	28.5	–	17.50
Ca2+	130.26	137.00	120.84	111.00	106.81	111	–	114.00
Mg2+	23.09	22.20	22	21.70	21.15	22.1	–	20.70
Cl−	71.97	70.30	60.62	65.80	55.66	59.4	–	47.60
	119.74	119.00	104.08	116.00	100.14	105	–	86.60
	318.43	299.36	297.2	262.70	272.94	262.7	–	268.81

View Large

According to the test results of ²H and ¹⁸O sampling analysis (Table 4), the percentage of spring water recharge sources in each division is calculated by using IsoSource linear model (Table 5). The calculation results of 10 recharge sources are summarized according to the type of water-bearing medium, recharge direction and recharge type, and the mixing ratio of spring water sources of four large spring groups is obtained (Table 6).

According to the correlation between different aquifers and Baotu Spring water level, the correlation between Zhangxia Formation aquifer water level and Baotu Spring water level in recharge area Z1 was poor, and R² was only 0.26. This indicated that the regional groundwater flow system participated in by the Zhangxia Formation aquifer had little contribution to Baotu Spring. The correlation between the water level of the aquifer of the Ordovician Yeli and Liangjiashan formations in the runoff area Z2 and the water level of Baotu Spring was general, and R² was 0.63. Although the aquifers of Z3 and Z4 located in Jingshi Road and close to the four spring groups were the same, there was a great difference in correlation with the water level of Baotu Spring. Z3 had a good correlation with the water level of Baotu Spring, R² = 0.84, while Z4 had a poor correlation with the water level of Baotu Spring, R² = 0.45. Following the tracer tests at Z3 and Z4 (Chi 2019; Li 2021), the tracer injected by Z3 was not detected only in the Heihu Spring group, and different groundwater flow directions were identified according to the variation of tracer concentration at different depths. The tracer injected by Z4 was not detected in the four large spring groups. This demonstrates that there are multiple groundwater flow systems at different levels in the two places.

The surface lithology of the exposed position of Xiaoling Spring is a quaternary system, the underlying Ordovician stratum has been missed, and the surrounding mountains are exposed by the Cambrian Fengshan Formation. After precipitation infiltration, the head pressure of the Cambrian Fengshan Formation aquifer increases. Due to the relative water isolation of the underlying Cambrian Gushan and Changshan Formation, it discharges and flows in Xiaoling Spring and participates in the local groundwater flow system greatly affected by seasonal changes.

In addition, the weathering layer of Archaean metamorphic rocks, Cambrian Mantou Formation, Cambrian Zhangxia Formation, and Ordovician Yeli and Liangjiashan Formation aquifers in the southern mountainous area all have seasonal springs participating in the local groundwater flow system, such as Qingsong Spring, Daling Spring, Heiyu Spring, Ruizhu Spring, etc. The outflow time is from July to October with large rainfall, and the rest of the time is cut off. This shows that the local groundwater flow system is widely distributed and basically covers the whole Jinan spring domain under the influence of topography.

Field observations showed that borehole X1 experienced artesian flow many times after heavy rainfall in August 2019 and July 2021. By analyzing the fixed depth conductivity value combined with the formation data, we found that the conductivity curve presented a ‘three-stage’ structure affected by three different groundwater flow systems (Figure 9(a)). Due to the different recharge sources and runoff paths of different groundwater flow systems, the conductivity of groundwater varied greatly. The analysis revealed that the artesian flow of the X1 borehole in 2019 and 2021 belonged to the artesian flow of different groundwater flow systems. Since August 10, 2019, the duration of heavy rainfall is three days. The statistical rainfall of the rainfall station nearest to the borehole was 323.8 mm. The conductivity test results after gravity flow manifested that the buried depth was 170 m, the shallow conductivity value was the fixed value of 860 μS/cm, and there was a sudden change in the conductivity value at the buried depth of 100 m–130 m. We inferred that there should be a large groundwater seepage channel here, which is discharged along the borehole under the influence of direct precipitation recharge. Since July 27, 2021, the duration of heavy rainfall is three days, and the rainfall is 116 mm. The conductivity test results after artesian flow suggested that the buried depth was 110 m, and the shallow conductivity value remained 585 μS/cm. The comprehensive analysis demonstrated that the conductivity value after gravity flow in 2019 was 860 μS/cm which was meaningfully bigger than the conductivity value of 585 μS/cm after gravity flow in 2021. The conductivity value in the depth of 170 m after the artesian flow was equal to that before the artesian flow, and the conductivity value in the depth of 170 m after the artesian flow was less different from that in the second section of the ‘three-stage’ structure. Therefore, we deduced that the artesian depth in 2019 was greater than that in 2021, belonging to the artesian flow of the intermediate groundwater flow system participated by the Fengshan Formation. The conductivity value after gravity flow in 2021 was 585 μS/cm, which was the same as the conductivity value of section I-I’ in the ‘three-stage’ structure before artesian flow. Because of the influence of sealing and water stop during drilling, the fixed depth curve did not reflect the conductivity change of the local groundwater flow system participated by the Majiagou Formation. Therefore, the artesian flow in 2021 belonged to the confined artesian flow of the intermediate groundwater flow system participated by the Yeli and Liangjiashan formations. Artesian well is a typical hydrogeological phenomenon of groundwater circulation. Thus, studying the development mechanism and relevant characteristics of an artesian well can reveal the hydrodynamic characteristics of the groundwater flow system and can provide a reference basis for the identification of a multistage groundwater flow system (Hou et al. 2008; Wang 2015).

Boreholes X1, X2, X3, and X4 were located in the direct recharge area of the spring area, and boreholes X5 and X6 were located in the runoff discharge area. By comparing the conductivity values of boreholes, we found that although the water intake sections of boreholes X1 and X5 were Ordovician Yeli and Liangjiashan formation aquifers, the conductivity of X5 was about 600 μS/cm greater than that of X1. Combined with the geological profile, the groundwater runoff recharge path of borehole X1 was much smaller than that of borehole X5, so the conductivity value was too large. The shallow Ordovician Majiagou Formation aquifer has a high degree of karst development, small groundwater circulation depth, and fast flow velocity. Therefore, the conductivity value of borehole X6 was large, and the conductivity value was 1,300 μS/cm–1,400 μS/cm.

To sum up, according to the conductivity value of groundwater in different aquifers, combined with the length of groundwater runoff path and the influence of karst development degree, it can be identified that the Ordovician Majiagou Formation is mainly involved in the local groundwater flow system, the Ordovician Yeli, Liangjiashan Formation and Cambrian Fengshan Formation are mainly involved in the intermediate groundwater flow system, and the Cambrian Zhangxia Formation is mainly involved in the regional groundwater flow system, and the groundwater flow system at all levels is affected by the development degree of karst at different depths, and many sub-flow systems are developed.

The study area can be divided into two areas. Region I is mainly the recharge runoff area, which is recharged by atmospheric precipitation. The alternating intensity of groundwater circulation is large and the runoff path is short, resulting in a short time of water rock interaction, making [Sr²⁺]/[Ca²⁺] small, and TDS increase gradually with the runoff path. Region II is mainly the runoff and discharge area of spring groups. The runoff path of karst water from the exposed limestone area to the discharge area is long, the water rock interaction time is long, the Sr²⁺ concentration is large, and there is a large [Sr²⁺]/[Ca²⁺]. The [Sr²⁺]/[Ca²⁺] in the sampling and test results of the four large spring groups are greater than 0.0025, distributed between 0.0025 and 0.00425. The distribution ranges of [Sr²⁺]/[Ca²⁺] of the Zhenzhu Spring group, Wulongtan Spring group, and Baotu Spring group are slightly different, which are 0.0035–0.0040, 0.0032–0.0041 and 0.0036–0.0039 respectively, while the [Sr²⁺]/[Ca²⁺] of Heihu Spring group is significantly smaller, which is 0.0026–0.0031, indicating that the four large spring groups are supplied by different groundwater flow systems and different runoff paths. Combined with the analysis of tritium isotope dating and the calculation results of the mixing ratio of stable isotope recharge sources, the age of Tanxi Spring and Wangfuchizi is 30a–35a. The karst water of the Cambrian Zhangxia Formation, which is mainly involved in the regional groundwater flow system, has the largest recharge proportion to the Zhenzhu Spring group, with a recharge proportion of 36%, followed by Wulongtan Spring group, with a recharge proportion of 23%. Therefore, the average residence time of Wangfuchizi is greater than Tanxi Spring. The age of Baotu Spring and Heihu Spring is 6a–10a. The recharge proportion of karst water from Ordovician Yeli Formation, Liangjiashan Formation and Cambrian Fengshan Formation, which mainly participate in the local groundwater flow system and intermediate groundwater flow system, to Baotu Spring group and Heihu Spring group is 28 and 54% respectively. Therefore, the horizontal residence time of Baotu Spring is greater than Heihu Spring. To sum up, Pearl Spring group and Wulongtan Spring group are replenished by the regional groundwater flow system with a long runoff path, and the spring water circulation renewal rate is slow, while Baotu Spring group and Heihu Spring group are replenished by the intermediate groundwater flow system and local groundwater flow system with relatively short runoff path, and the spring water circulation renewal rate is fast.

The metamorphic rocks of the Neoarchean Mount Tai Group in the spring domain have large thicknesses and low fracture development. For this reason, it is difficult to form deep underground runoff. Usually, seasonal fissure springs with a small flow are formed, which converge to the valley rivers with the surface base flow and participate in the local groundwater flow system. The aquifer of Cambrian Mantou and Xuzhuang formations has general karst development and limited water storage space, so it is hard to form deep runoff. It is usually discharged into springs in the valley. This aquifer group that supplies the Zhangxia Formation aquifer locally is affected by structural faults and has little contribution to the four spring groups. Hence, this aquifer is also mainly involved in the local groundwater flow system. The recharge elevation of the Cambrian Zhangxia Formation aquifer in the indirect recharge area of four large spring groups is higher than that in the direct recharge area (Wang 2016). The surface karst and deep fractures are developed well. They receive the recharge of atmospheric precipitation to form underground runoff and then migrate to the deep. Some of them recharge the Ordovician aquifer of the Cambrian Fengshan Formation through faults and participate in the intermediate groundwater flow system. Furthermore, some of them expose to springs that are in the contact fracture zone of igneous rocks and participate in the regional groundwater flow system. The remaining migrate to the deep underground under the action of gravity to replenish the groundwater north of the Yellow River (Li 2021).

The comprehensive analysis shows that the Cambrian Zhangxia Formation aquifer in the study area mainly participates in the regional groundwater flow system. The Cambrian Fengshan Formation and Yeli and Liangjiashan Formation aquifers are mainly involved in the intermediate groundwater flow system. In the area with deep fissures, the groundwater will also carry out deep circulation under the action of gravity and participate in the regional groundwater flow system. The aquifer of the Ordovician Majiagou Formation mainly participates in the local groundwater flow system with shallow burial. In the fracture development area, there is a recharge effect of gravity on the aquifers of the Yeli, Liangjiashan Cambrian, and Fengshan Formation and participates in the intermediate groundwater flow system. Affected by the terrain, the Cambrian Zhangxia and Fengshan Formation, Ordovician Yeli and Liangjiashan Formation may participate in the local groundwater flow system and are mostly discharged in the form of descending springs.

The contents of δ²H and δ¹⁸O in the four large spring groups are different, indicating that the four spring groups receive mixed recharge from hierarchical groundwater flow systems participated in by different aquifers of Ordovician and Cambrian at the same time, and the recharge proportion is different. The calculation results of IsoSource model show that the aquifer of Ordovician Cambrian Fengshan Formation, which mainly participates in the local and intermediate groundwater flow system, accounts for the largest proportion of recharge to the four large spring groups, and Heihu Spring group > Wulongtan Spring group > Baotu Spring group > Pearl Spring group, accounting for 54, 33, 28, and 19% respectively. The Cambrian Zhangxia Formation aquifer, which is mainly involved in the regional groundwater flow system, accounts for a relatively small proportion of the recharge to the four large spring groups, and the recharge proportion of Pearl Spring group > Wulongtan Spring group > Baotu Spring group > Heihu Spring group accounts for 36, 23, 21, and 5% respectively. In addition, the recharge proportion of artificial reinjection water to the four large spring groups is also large, and the recharge proportion of Baotu Spring group > Pearl Spring group > Heihu Spring group > Wulongtan Spring group accounts for 42, 36, 33, and 32% respectively. The direct recharge area of spring water of the four large spring groups is the main recharge area of local groundwater flow system and intermediate groundwater flow system. Therefore, from the perspective of spring protection, under the condition of protecting the direct recharge area of spring water, the strata of Ordovician and Cambrian Fengshan Formation should be the best choice for artificial recharge.

It should be noted that the cyclic evolution of groundwater is an extremely complex hydrogeological process, and the recharge and mixing of groundwater in different periods are different. The recharge and mixing ratio of the multistage groundwater flow system to the four large spring groups were calculated by the samples taken in this study only to represent the static value at a certain time.

• (1)

  According to the ISM model, the groundwater ages of Baotu Spring and Heihu Spring are 6a–10a, and the groundwater ages of Tanxi Spring and Wangfuchizi are 30a–35a. ISM model fully considers the hydrodynamic dispersion of groundwater in a karst aquifer medium, and the calculation is simpler. Combined with hydrogeological conditions and traditional isotope mathematical model, a more accurate groundwater age distribution range can be obtained. Therefore, ISM model can be used as a basic method to determine the age of groundwater.

• (2)

  Through the identification of the correlation between the groundwater level and the water level of Baotu Spring, the characteristics of hydrogen and oxygen isotopes of groundwater, groundwater spillover zone, artesian well, fixed depth conductivity curve, ion concentration ratio, and other factors, it is considered that the recharge proportion of the four large spring groups by different levels of groundwater flow system is different. The local groundwater flow system and intermediate groundwater flow system mainly participated by the aquifers of the Ordovician and Cambrian Fengshan Formation replenish the Heihu Spring group and Wulongtan Spring group in a large proportion, accounting for 54 and 33% respectively, and the Baotu Spring group and Pearl Spring group in a relatively small proportion, accounting for 28 and 19% respectively. The regional groundwater flow system mainly participated in by the Cambrian Zhangxia Formation aquifer has a large recharge proportion to the Pearl Spring group and Wulongtan Spring group, accounting for 36 and 23% respectively, and a relatively small recharge proportion to the Baotu Spring group and Heihu Spring group, accounting for 21 and 5% respectively.

• (3)

  In the karst water system of the Jinan monoclinal structure, the hierarchical groundwater flow system is divided according to different aquifer structures. The shallow Ordovician Majiagou Formation aquifer is mainly involved in the local groundwater flow system, the Ordovician Yeli, Liangjiashan Formation aquifer and Cambrian Fengshan Formation aquifer are mainly involved in the intermediate groundwater flow system, and the Cambrian Zhangxia Formation aquifer is mainly involved in the regional groundwater flow system, which is affected by the terrain, each aquifer may participate in the local groundwater flow system. Among them, the Cambrian Fengshan Formation Ordovician aquifer participates in the intermediate groundwater flow system with the most complex characteristics and the largest contribution rate to the four large spring groups, which has great scientific research and practical significance.

Thanks to the Shandong Engineering Technology Institute for Groundwater Numerical Simulation and Contamination Control and Shandong Provincial Geo-mineral Engineering Exploration Institute for their technical and data support for this study.

This research was funded by the National Natural Science Foundation of China (41772257) and the innovation team project of colleges and universities in Shandong Province (2018GXRCO12).

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

The authors declare there is no conflict.