分级地下水流系统中地下水年龄分布的数值模拟——以济南为例

 摘要

        了解地下水年龄对含水层管理至关重要。本研究选取济南市一个典型的水文地质断面，建立了地下水流动与年龄输送的耦合模型。模拟结果表明，泰山北麓至齐光断裂发育由嵌套的局部子系统、中间子系统和区域子系统组成的分级地下水流动系统。这些流动子系统的规模主要由地形控制。含水层岩性变化也会影响地下水流动系统。地下水年龄分布与流动系统密切相关。地下水年龄沿流动路径变老，并在内部停滞点附近突然增加，这可以作为不同地下水流动层次之间的特征边界。研究结果不仅有助于了解济南市地下水循环规律，促进地下水资源的合理开发利用，而且为地下水流动系统理论的实际应用提供参考。

1. Introduction

        地下水的年龄作为一种内在性质，承载着地下水循环和演化的重要信息。因此，了解地下水年龄对于水资源管理至关重要。地下水年龄的评估是评估含水层可再生性、补给性和脆弱性的关键，也有助于确定地下水流动路径、估计含水层特性和地下水速度以及补给率[1-3]。
        通常使用两种方法来估计地下水年龄：（i）测量地下水中的天然或人造示踪剂；（ii）通过数学求解输运方程来建模地下水年龄。当使用环境示踪剂来评估地下水系统中的年龄分布时，地下水输送过程的复杂性通常会大大降低。Phillips和Castro、Kazemi等人[4，5]对用于地下水年龄测定的环境示踪剂及其应用和缺陷进行了综述。 
        地下水年龄建模方法通常分为三种方法，即平流模型、溶质运移模型和年龄-质量运移模型[6]。通过比较通过这三种不同建模方法模拟的Carrizo含水层的年龄分布，发现当混合过程不可忽略且地下水速度变化严重时，可以通过年龄-质量传输方法获得最一致的年龄分布[7]。因此，年龄-质量传输方法应该是计算分级地下水流系统中地下水年龄分布的首选方法，其中流速表现出强烈的变化[6]。 
本研究是在中国北方典型的泉水集水区济南进行的。108个泉水分布在市区。随着当地经济的发展，地下水开采量不断增加，但缺乏科学合理的规划和调控。因此，济南的许多泉水都面临着流量衰减甚至干涸的问题。因此，迫切需要全面了解该地区的地下水流动系统。 
济南以往的研究主要涉及地质结构、地下水位和泉水流出动力学[8，9]。王等人根据地下水排泄区的水文和化学信息，将济南的地下水流动系统概念化为局部、中间和区域子系统的层次结构[10]。然而，对济南市地下水流系统的定量研究很少，这对于准确管理和保护地下水资源是必要的。 
本文以济南市地下水分层流动系统为研究对象，进行了现场实例分析。采用地下水流与年龄输运耦合模型对不同的地下水流动层次和年龄分布进行了定量分析。 

2. Geological setting

Jinan is located in a monoclinic geological structure, with topographical elevation gradually declining from south to north. Land view transits successively from mountain terrain in the south, to the piedmont plain, and finally to the alluvial plain in the north (figures 1 and 2).
Metamorphic rocks (Art), exposed in the south, constitute the basement of the aquifer system in Jinan. Limestone and shale are interbedded in the Cambrian strata, overlaying the Art basement. The overlying Ordovician, composed by the thick layers of limestone, constitutes the most productive aquifer in Jinan. Mesozoic igneous rocks, composed by diorite and gabbro, exist in the northern Jinan city and are mostly covered by Quaternary deposits. 

Figure 1. Geological setting of the study area. A–A′ refers to the cross section shown in figure 2
(modified after [10]).

Figure 2. Typical geological cross section in the study area (modified after [10]).

Karst groundwater generally flows from south to north, coinciding with the regional topography and dip direction of karst strata. Karst groundwater mainly discharges in the form of ascending springs in Jinan city as blocked by the igneous rocks. Regional groundwater continues to move northward with deeper circulation depth and longer flow path. The overlying Carboniferous, Permian, Triassic and Neogene constitute the relative aquicludes in the system. 
The cross section A-A’ (figure 2) stretches from the divide of Mount Taishan in the south to Qiguang Fault in the north, covering all the landforms and topographical features in the study area. This cross-section is parallel to the regional groundwater flow direction, and goes through all the aquifer groups in Jinan. Therefore A-A’ is a typical cross section in studying the Jinan groundwater flow system. In this paper, a two-dimensional numerical model is constructed based on the profile A-A’. 

3. Mathematical models

The two-dimensional steady flow assumption has been commonly applied to regional scale subsurface flow and transport model [11,12]. It is often supposed that the groundwater age distribution is constant over many years. So the groundwater flow equation under steady state is used in the present study to calculate the subsurface flow field:

where K is the hydraulic conductivity tensor, h is the hydraulic head. 
In this paper, the age mass transport approach is employed to evaluate groundwater age distribution, given its advantages over other two age modeling methods. The equation for the age mass transport under steady state is [13]: 

where τ is the groundwater age, θ is the effective porosity, u is the pore-water velocity vector, D is the dispersion coefficient tensor. 
The two-dimensional coupled groundwater flow and age mass transport equations under steady state are solved numerically through finite element calculation via FEFLOW software (version 6.2, DHI-WASY). 

4. Numerical model

The profile A-A’ (figure 2) comprises all types of aquifer groups in the study area, including: (1) porous aquifers composed by Quaternary sediments (Q), (2) fissure karst aquifers composed by Ordovician carbonate rocks (O), (3) karstic fissure aquifers composed by Cambrian clastic rocks with interbeds of carbonate rocks (Є), (4) fissure aquifers composed by the metamorphic rocks of Archean Taishan Group (Art). Gabbro (γ) as well as Carboniferous (C), Permian (P), Triassic (T), Jurassic (J), Neogene (N), which are mainly composed by sandstone, mudstone and shale with relatively low permeability, are divided into aquicludes. Based on the hydraulic properties of each strata in the profile A-A’, the model structure can be generalized as shown in the figure 3 with the length of 88 km and altitude difference of approximately 3400 m. 
The model domain is bordered by the average phreatic surface on the top, the divide of Mount Taishan at the south, and the impervious Qiguang Fault at the north (figure 3). The maximum circulation depth of regional groundwater in Jinan is calculated to be 1300-1900 m [14]. Given the rock permeability decreasing with depth, the depth of 2000 m is set as impervious basement of the model. 

Figure 3. Schematic diagram of the generalized model structure and boundary conditions.

The model domain is discretized into triangular meshes, with the maximum element size of 30 m. The mesh has 173960 nodes and 342648 elements. The top boundary for the groundwater flow model is specified hydraulic heads and three other boundaries are impervious. For the age mass transport model, the impervious boundaries correspond to the zero age mass flux, the recharge zones on the top correspond to the zero age mass, and the discharge zones on the top correspond to the zero dispersive age mass flux. 

Figure 4. The relationship between hydraulic conductivity and depth in Ordovician limestone aquifer. The karst development degree is represented by the percentage of boreholes discovering caves in Jinan. The permeability decaying tendency is approximately represented by the attenuation of karst development with depth.

The effective porosity (θ), hydraulic conductivity (K) and dispersity are the principal parameters for the subsurface flow and age mass transport model. The parameters are mainly acquired from numerical simulation and pumping tests as well as tracing experiments conducted in the previous study [9,15,16]. Measured groundwater age and flow paths determined by the hydrochemical and isotopic data are also used to calibrate model parameters. 

Since the permeability of karst aquifer is controlled by the degree of karst development, the attenuation of karst development with depth can approximately represent the permeability decaying tendency. The hydraulic conductivity of Ordovician limestone aquifer exhibits exponential distribution with decay exponent of 0.003 m-1 (figure 4). 
The calibrated horizontal hydraulic conductivity Kx and effective porosity θ are listed in table 1. The anisotropy ratio Kx/Kz is set as 50, and longitudinal dispersity being 80 m, lateral dispersity being 8 m, which make the simulation results more reasonable. The longitudinal dispersity is about 1/1000 of the cross-section length, coinciding with Castro and Goblet and Jiang et al [6,7]. 

5. Simulation results and discussion

5.1. Groundwater flow distributions 

The distributions of subsurface flow system and groundwater age in the profile obtained from the above numerical model are exhibited in figure 5. The primary effect of water-table undulations on basinal water flow is the generation of hierarchically nested flow systems of different orders. Local flow system is distributed in the upper part, underlain by the intermediate and regional flow system. 

Figure 5. The distributions of subsurface flow system and groundwater age in the profile A-A’.

Affected by the topographic relief of the mountain terrains in the south, groundwater mainly discharges at the lower area nearby, forming the local flow system. The local flow system is developed in the shallow part with short flow path and small penetration depth. The intermediate flow system is sourced from the southern piedmont, and then mainly flowing through the carbonate aquifers with longer flow paths. This subsystem is controlled by the macroscopic topography and permeability of  the aquifers. The regional flow system is characterized by the longest flow paths and smallest flow velocities. The regional landform features and large-scale hydrogeological structures dominate the development of this subsystem. The recharge of regional groundwater is from the southern ranges, and then flowing through deep fault and fracture zones with the longest flow paths and largest circulation depths driven by gravity. 
The influence of lithologic changes on the groundwater flow system can be directly exhibited through the groundwater flowline. The flowline is refracted at the stratigraphic interface, and the larger lithologic difference between aquifers, the more obvious the flowline refraction is. The influence of gabbro rocks on the groundwater flow system is also exhibited in figure 5. Local and intermediate groundwater flow are resisted by the gabbro rocks and converged around the contact zone between limestone and igneous rocks. Then it may ascend in the form of springs in fissure-developed areas. Deep regional groundwater is forced to move downward and bypass the gabbro rocks, then entering the northern Jinan geothermal field. 
5.2. Groundwater age distributions 
The groundwater age simulation result is also exhibited in figure 5. Shallow Ordovician karst groundwater ages in discharge area (less than 200 m depth) vary between 10-30 a. Groundwater ages in the deep Ordovician thermal reservoir vary between 13-35 ka. The simulation result coincides with the measured groundwater age based on CFCs and 14C dating method. 
Groundwater age distribution is closely related to the flow system, with younger groundwater near recharge area and becoming older along the flow path. Abrupt change of age distribution can be found in the regional discharge area. Groundwater ages under the regional discharge area are generally larger than 40000 years, and even reaching 99000 years at the basin bottom. Thus it can be inferred that there exists stagnant zone under the regional discharge area. 
An important characteristic of hierarchically nested flow system is the existence of stagnation points or stagnant zones between flow subsystems. The stagnation points can be divided into two types, i.e. basin-bottom stagnation points and internal stagnation points [17]. According to the simulation result (figure 5), three basin-bottom stagnation points SPR1, SPR2, SPR3 are found in the groundwater flow system. SPR1, SPR3 are located at the bottom of regional groundwater divide and discharge area respectively. Groundwater at SPR3 is the oldest in the basin, owing to its longest flow path and smallest velocity. SPR2 is located at the bottom of gabbro rocks, reflecting the influence of gabbro rocks on groundwater flow system.

Figure 6. Age profiles across the internal stagnation points SPL1, SPL2, SPL3.

Internal stagnation points SPL1, SPL2, SPL3 are located inside the basin where the divergence and convergence of four flow subsystems coexist. The shape of groundwater age contours near the internal stagnation points is similar to a ridge pointing to the discharge area. 
Virtual boreholes penetrating stagnation points are used to exhibit the vertical groundwater age variations (figure 6). The profiles of groundwater age distribution exhibit that it reaches a local maximum near the stagnation points. When exceeding the influence zone of stagnation points, groundwater age changes become moderate. Groundwater below SPL1, an internal stagnation point near the recharge area, becomes older with depth. But groundwater ages below SPL2 and SPL3, internal stagnation points near the discharge area, firstly decrease and then increase with depth. This age change tendency is controlled by the regional groundwater flow direction. Regional groundwater in deep Ordovician aquifer is forced to move downward under the gabbro rocks, causing groundwater becoming older with depth. Then groundwater flow passes through the Ordovician aquifer and moves upward to the regional discharge area, which in turn leads to groundwater becoming younger with depth. 
The heterogeneity of aquifer system can also influence the groundwater age profile, exhibiting the inflection point at the interface of different strata (figure 6). But the influence of lithologic change on groundwater age is much smaller than that of the special hydraulic property of stagnation points. 

6. Conclusions

A 2-D coupled subsurface flow and age transport model is constructed to quantitatively analyse the distributions of subsurface flow system and groundwater age in the typical cross section in Jinan. The major findings of this study can be summarized as follows:
A hierarchically nested groundwater flow system, including local, intermediate and regional subsystems, develops in the study area. Topography and geological structure control the scale of groundwater flow hierarchies with different flow path length and circulation depth. 
Groundwater age distribution is closely related to the subsurface flow system. The age of groundwater in different hierarchies is proportional to its flow path length and circulation depth. Local groundwater is much younger than the intermediate and regional groundwater. Abrupt change in groundwater age profile is the characteristic boundary between different groundwater flow systems. 
Hydraulic heterogeneity is an important factor influencing groundwater flow and age distributions. Flowline is refracted at the stratigraphic interface and bypassing the impermeable gabbro rocks, which would further affect the groundwater age distribution. 
Internal stagnation point is an important character of the nested groundwater flow system, and can be used to divide different flow hierarchies. Groundwater age abruptly changes and local maximum age value exists near the stagnation point. 
The results of this study can contribute to the future application of groundwater flow system theory, such as interpreting groundwater tracer-age and hydrochemical evolution. The results of this study can also contribute to the groundwater resources management, especially guaranteeing the perennial outflow of karst springs in Jinan, and protecting groundwater quality. In the future, more hydrogeological work such as groundwater age dating and tracer test is needed to verify the numerical simulation results in this study. 
References 
[1] Tóth J 1999 Groundwater as a geologic agent: an overview of the causes, processes, and
manifestations Hydrogeol. J. 7 1-14 
[2] Tóth J 2009 Gravitational Systems of Groundwater Flow: Theory, Evaluation, Utilization (New
York: Cambridge University Press) 
[3] Gomez J and Wilson J 2013 Age distributions and dynamically changing hydrologic systems:
Exploring topography-driven flow Water Resour. Res. 49 1503-22 
[4] Phillips F and Castro M 2003 Groundwater dating and residence time measurements Treatise on Geochemistry ed J Drever (Oxford: Oxford University Press) pp 451-97 

[5]    Kazemi G, Lehr J and Perrochet P 2006 Groundwater Age (Hoboken: Wiley–Interscience)
[6]    Jiang X, Wan L, Ge S, Cao G, Hou G, Hu F, Wang X, Li H and Liang S 2012 A quantitative study on accumulation of age mass around stagnation points in nested flow systems Water Resour. Res. 48 W12502 

[7] Castro M and Goblet P 2005 Calculation of ground water ages–a comparative analysis Ground
Water 43 368-80 
[8] Qian J, Zhan H, Wu Y, Li F and Wang J 2006 Fractured-karst spring-flow protections: a case
study in Jinan, China Hydrogeol. J. 14 1192-205 
[9] Kang F, Jin M and Qin P 2011 Sustainable yield of a karst aquifer system: a case study of Jinan
springs in northern China Hydrogeol. J. 19 851-63 
[10] Wang J, Jin M, Lu G, Zhang D, Kang F and Jia B 2016 Investigation of discharge-area groundwaters for recharge source characterization on different scales: the case of Jinan in northern China Hydrogeol. J. 24 1723-37 

[11] Castro M, Goblet P, Ledoux E, Violette S and Marsily G 1998 Noble gases as natural tracers of water circulation in the Paris Basin: 2. Calibration of a groundwater flow model using noble gas isotope data Water Resour. Res. 34 2467-83 

[12] Castro M and Goblet P 2003 Calibration of regional groundwater flow models: Working toward a better understanding of site-specific systems Water Resour. Res. 39 1172 

[13]    Goode D 1996 Direct simulation of groundwater age Water Resour. Res. 32 289-96
[14]    Wang J, Jin M, Jia B and Kang F 2015 Hydrochemical characteristics and geothermometry applications of thermal groundwater in northern Jinan, Shandong, China Geothermics 57 185-95 

[15] Cai W, Gao Z, Wang Q, He Q and He D 2013 Research on Connections between JINAN Karst
Waters (Beijing: Geological Publishing House) (in Chinese) 
[16] Li Z, Kang F, Liu G, Xu M and Pan C 2013 Jinan Geothermal and Hot Springs (Beijing:
Geological Publishing House) (in Chinese) 
[17] Jiang X, Wang X, Wan L and Ge S 2011 An analytical study on stagnant points in nested flow systems in basins with depth-decaying hydraulic conductivity Water Resour. Res. 47 W01512