Next Article in Journal
Natural and Human-Induced Flow and Sediment Transport within Tidal Creek Networks Influenced by Ocean-Bay Tides
Previous Article in Journal
Climate Variability and Climate Change Impacts on Land Surface, Hydrological Processes and Water Management
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Simulation of Deep Thermal Groundwater Exploitation in the Beijing Plain Area

1
School of Water Resources and Environment, China University of Geosciences (Beijing), Beijing 100083, China
2
MOE Key Laboratory of Groundwater Circulation and Environmental Evolution, China University of Geosciences (Beijing), Beijing 100083, China
3
School of Science, China University of Geosciences (Beijing), Beijing 100083, China
*
Author to whom correspondence should be addressed.
Water 2019, 11(7), 1494; https://doi.org/10.3390/w11071494
Submission received: 3 June 2019 / Revised: 4 July 2019 / Accepted: 13 July 2019 / Published: 18 July 2019
(This article belongs to the Section Water Resources Management, Policy and Governance)

Abstract

:
Thermal groundwater is relatively abundant in the deep-seated bedrock underlying the Beijing plain area. The main geothermal reservoir is composed of dolomites of the Wumishan Group of the Meso–Neoproterozic Jixian System. The thermal groundwater has been developed and utilized since the 1970s and significant declines in groundwater levels were observed. A 3D unsteady flow model of an anisotropic karst-fissure aquifer based on the equivalent continuum is used to describe the flow of thermal groundwater and heat transport. The heat transportation is described by the governing equation including convection and dispersion. The simulation of this paper aims to solve such problems as uneven distribution and thinness of the aquifer, insufficient initial monitoring data, and poor knowledge of the properties of the horizontal boundary. They are solved by considering vertical stratification of the aquifer with equal thickness, replacing initial water level data by surface elevation, and choosing natural boundary far away from the exploitation areas. Through a trial–error procedure, the simulated and measured groundwater level and temperature in the simulation period are well fitted. Three exploitation schemes are proposed to predict the spatial and temporal changes in groundwater level and temperature of the thermal groundwater in the study area. The prediction results show that the reinjection can effectively slow the decline in the thermal groundwater levels. Except for the Dongnanchengqu, Xiaotangshan, and Liangxiang subgeothermal fields, the other six subgeothermal fields have the potential for further development of thermal groundwater.

1. Introduction

Located in the northern part of the North China Plain Basin (Figure 1), Beijing is one of the cities in China underlain by relatively abundant geothermal resources [1]. Thermal groundwater in the deep-seated carbonates in the Beijing plain area is of sedimentary basin type of low to medium temperature. The geothermal reservoir in the study area is composed of the Proterozoic siliceous dolomite and Cambrian–Ordovician limestones, with a total thickness of more than 5000 m. Especially, the carbonates of the Wumishan Group of the Jixian System of the Meso–Neoproterozoic with a thickness of more than 2000 m (m is short for meter) are widely spread, accounting for over 90% of the basement rocks in the plain area. The Wumishan carbonates is characterized by highly-developed karst landform and thus become the most important geothermal reservoir in the study area. The burial depth of the top of the geothermal reservoir varies greatly from about 150 m to 4000 m or more (at the center of Beijing depression). The caprocks of the geothermal reservoir are mainly composed of Mesozoic and Cenozoic clastic rocks and range in thickness from tens of meters to several kilometers. The well head temperature of the thermal groundwater ranges from 48 °C to 95 °C.
Exploration of thermal groundwater in the bedrock aquifers in the Beijing plain area can be dated back to the year 1971 [2]. Since then, number of production wells and total exploitation keep on increasing. Due to the over-exploitation of thermal groundwater, the groundwater level of the geothermal reservoir continues to decline. By monitoring the amount of production, water level, and hydrochemistry of the thermal groundwater, it is found that there is a positive correlation between the amount of production and the decreasing water level. Simulation and prediction of the temporal and spatial evolution of groundwater level and temperature in the geothermal system are of important significance, even though the present data of geological exploration and routine observation of groundwater level and temperature at some geothermal wells are insufficient.
Thermal groundwater in deep bedrock in a sedimentary basin is a special type of geothermal water. It is characterized by deep burial, great change of temperature and heat exchange with surrounding rock [3,4]. In the study of such a complex thermal groundwater system, numerical simulation techniques can systematically reveal the distribution, migration, and formation mechanism of thermal groundwater; analyze the resource potential of the whole geothermal system; and help to formulate reasonable exploitation planning. Numerous mathematical models describing thermal groundwater flow in a geothermal system have been developed during the past decades. Voss [5] created a complex model of energy transport in the subsurface system, which was two dimensions and could accurately reproduce the flow and energy with proper discretization. Eckstein and Maurath [6] performed a simulation in the Momotombo geothermal area of Nicaragua in Central America with a time-dependent model and the results of the simulation revealed the origin and velocity of the hot fluids. Zhang [7,8] studied the function and rule of natural convection in heat transport, and the flow equation and heat transfer equation describing natural convection in the heat transfer of aquifers were derived. The natural convection was well demonstrated by simulating the energy storage of a group of wells in Shanghai. Ondrak and Moller [9] used a numerical model coupling heat and mass transport to simulate the temporal and spatial evolution of temperature of a shear zone and its vicinity in Germany, and obtained a better understanding of the time scales at which the hydrothermal vein deposits form. Dilsiz [10] used a 3D hydrogeological conceptual model to explain the occurrence of the Pamukkale hot spring in the southwestern Turkey. Xue [11] discussed the problem of heat transport in the saturated zone and energy storage in aquifers, and the numerical solution of the corresponding model was also given. Ostermann [12] developed a 3D model to simulate the heat transport in hydrothermal systems which was based on a transient advection–diffusion-equation for a porous medium and proved that Galerkin schemes could be available for the numerical realization. Diersch [13] discussed the finite-element computation of heat transport in porous media and found that it could be applied to the nonisothermal porous-medium processes, for instance, exploitation of geothermal reservoirs. The movement of thermal groundwater involves the flow and heat migration of hot water, and the 3D unsteady flow equations are needed to describe both the hot water flow and heat migration in a geothermal reservoir. However, such models are difficult to calibrate because they require field observations that are relative expensive to conduct on a regular basis over an extended period of time. In this simulation, a 3D unsteady groundwater flow and heat transport model in carbonates based on the equivalent continuum [14] is used to simulate the thermal groundwater system in the study area. Field observation data from some of the geothermal wells are used for a first calibration. Predictions of future changes in hydraulic head and temperature are also performed. The results may provide some insight into the geothermal structure of the North China Plain, and have implications for the development of geothermal water resources in similar sedimentary basins.

2. Description of the Study Area

At present, geothermal resources which can be exploited and utilized in Beijing mainly occur in the plain area, with an area of over 6000 km2. Tectonically, the Beijing plain area is located in the north of the China–North Korea Paraplatform. It is a faulted basin developed from the north China Platform and has undergone several tectonic changes and multi-stage and multi-cycle geological structural evolution [15]. The structural framework of Beijing formed during the Mesozoic Yanshan Movement and sedimentation basins with different ranges, different sediments and different subsidence periods formed due to the later differential movements. The basins are covered by Eogene and Quaternary unconsolidated sediments which range in thickness from tens of meters to nearly thousands of meters. The Jurassic, Cretaceous, Paleogene, Neogene, and Quaternary sediments with the maximum thickness of over 4000 m were deposited in the sedimentation center, namely, the central part of the Beijing depression.
A series of tectonic movements led to mainly two sets [15] of faults in the region (Figure 1). One set trends nearly in the NE direction, including the Babaoshan fault (F2), the Huangzhuang–Gaoliying fault (F1), etc. Another set is the NW direction, represented by the Nankou–Sunhe fault (F13) and the Yongdinghe fault (F11). The occurrence of geothermal anomalies are closely related to the active faults. Deep faults, such as the Huangzhuang–Gaoliying fault, can be regarded as deep flow and migration channels for thermal groundwater and also as boundaries of the geothermal fields. The faults concentrated in the middle of the Beijing geothermal field play a crucial role in increasing the development of karst-fissures in the carbonate aquifers, communicating deep heat sources, and improving geothermal gradients.
The strata from top to bottom in the geothermal system in Beijing include: the Ordovician System (O), Cambrian System (Cm), Qingbakou System (Qn), Tieling Group of the Jixian System (Jxt), Hongshuizhuang Group of the Jixian System (Jxh), Wumishan Group of the Jixian System (Jxw) [16]. The Wumishan Group of the Jixian System occurs all over the study area, while other strata are discontinuous in the horizontal direction and have a relative small thickness (Figure 2 and Figure 3). Borehole data [2,15,16] show that the dolomite of the Wumishan Group of the Jixian System has a large amount of water produced by single wells and is the most important aquifer in this area. Other strata also developed a certain degree of fracture and karst phenomenon, and the specific well-yield is smaller than that of the Jixian Wumishan Group (Table 1).
Isotopic data [17,18,19,20,21,22] show that the thermal groundwater in the study area is of meteoric origin. The 14C dating data [18,19,21,22] which show that the age of thermal groundwater in the bedrock aquifers increases gradually from northwest to southeast, indicating that the flow direction of the hot water is from northwest to southeast. As shown in Figure 4, precipitation infiltrates in the northern and western mountainous areas and groundwater flows to the deep aquifers of the Beijing plain area along the fractures and fault zones. During the migration process, it is heated by the heat flow from below and obtains heat. The huge Jixian System has well-developed crisscrossed joints and karst fissures, and provides favorable conditions for the flow and storage of thermal groundwater. The Cenozoic and Mesozoic clastic rocks above the geothermal reservoir form a good and thick insulating cover. The thermal groundwater enriches in the fault zones and is available for exploitation.
The Beijing plain area is located in the north of North China Basin and the characteristics of the geothermal field are consistent with those of the North China Basin in both the horizontal and vertical directions. By analyzing the drilling data, ground temperature logging data and the well head temperature of the geothermal wells [2,15,16] in Beijing, the contours of temperature at the elevation of −2000 m was obtained (Figure 5). The horizontal temperature of the geothermal field in the study area is characterized by a rising trend from northwest to southeast. The main vertical temperature characteristic is that the geotemperature gradient of the caprocks is higher than that of the geothermal reservoir (Table 2).
In the Beijing plain area, the area of the reservoir at the burial depth of less than 3500 m and with temperature of more than 50 °C is approximately 2760 km2. The Beijing geothermal field was previously divided into 10 subgeothermal fields [23] that are somewhat related to each other (Figure 6). In recent years, with the development of economic construction, the number of exploitation wells of geothermal resources is increasing and the exploitation amount of the hot water is also increasingly growing. At present, there are two mainly subgeothermal fields in which the exploration of thermal groundwater is relatively strong. One is the Dongnanchengqu subgeothermal field with an area of 117 km2 in the central city. The other is the Xiaotangshan subgeothermal field with an area of 20 km2 to the north of the city. By the end of 2016, there were about 200 exploitation wells of thermal groundwater in Beijing plain area, about half of which were concentrated in the the Dongnanchengqu subgeothermal field and Xiaotangshan subgeothermal field. The accumulative exploitation amount of thermal groundwater is over 1.3 × 108 m3 (m3 is short for cubic meter) (Table 3). As a result, the groundwater levels of the geothermal reservoir continues to substantially decline with an average drop of more than 2.5 m per year.

3. Mathematical Model

The carbonate geothermal reservoir in the study area is an underdeveloped karst landform and is of certain anisotropy and heterogeneity. However, the thermal groundwater in the reservoir flows slowly and has a uniform potentiometric surface. The groundwater flow is thought to obey the Darcy’s law. The thermal groundwater flow is affected by the viscosity and density of hot water, which are dependent on temperature. In this simulation, a 3D unsteady groundwater flow and heat transport model [24] in an anisotropic carbonate aquifer based on the equivalent continuous medium is established and the influence of groundwater density and viscosity on the simulated groundwater levels is taken into account.
The groundwater flow equation that describes the change of density with temperature is
x i [ K i j ( h x j + ρ r e j ) ] = S s h t + ϕ ρ r T × T t ( 1 + ρ r ) q
where K i j is the hydraulic conductivity and K i j = ρ g k i j / μ ; S s is the specific storage; ρ is density of water; μ is the viscosity coefficient; k i j is the permeability; h is the hydraulic head; t is the time; T is the fluid temperature; ρ r = { ρ ( P , T ) / ρ 0 1 } is the density variation of water; ϕ is the porosity of aquifers; e j is the j th component of the unit vector of the gravitational direction; q is the source and sink.
K i j ρ r e j is the convection caused by density difference and φ ( ρ r / T × T / t ) indicates the change in liquid mass caused by the density change in Equation (1). If the density is constant, Equation (1) becomes the common groundwater flow equation.
Heat transport is mainly controlled by convection and diffusion. The following mathematical equation including convection and diffusion is used to describe heat transport
c e T t + c ρ v i T x j = x i ( λ i j T x j ) + c ρ ( T T ) q
where c e is the heat capacity of aquifers; v i is the Darcy’s velocity; c is the heat capacity of water; λ i j is the thermal conductivity of aquifers; T is the temperature of the source and sink.
Combined with the initial conditions and boundary conditions of groundwater flow and temperature variation (mentioned in Section 4.3 and Section 4.4), Equation (1) and Equation (2) constitute the mathematical model of this simulation. FEFlOW, a groundwater numerical simulation software, is used for calculation.

4. Model Structure

Due to the complexity of geological structure in the Beijing plain area, the following difficulties exist in the process of establishing geological model and determining the initial conditions. First of all, there is no waterproof structure in the study area that can be treated as the ideal boundary (impermeable boundary) of the model. Secondly, the great difference in the roof elevation of the thermal reservoir and some discontinuous strata in the Beijing plain area make it difficult to layer the model in the vertical direction. Finally, geothermal water level data at the initial time of simulation are not sufficient to determine the initial distribution of geothermal water level. This simulation will focus on solving these problems.

4.1. Model Scope

The karst-fissured geothermal reservoir is of sedimentary basin type of low-medium temperature [25]. Currently, exploitation of thermal groundwater mainly concentrates in three subgeothermal fields including the Dongnanchengqu, Xiaotangshan, and Liangxiang. Exploitation wells in these areas are densely distributed, and the Wumishan Group of Jixian System is used as the main aquifer for exploitation. No geological structure with obvious water resistance property was found in the study area which can be treated as the simulation boundary, so the hydrogeological property of the simulation boundary that we choose is not clear. In order to decrease the influence of the boundary properties on the simulation results of the dense exploitation areas, the boundary of the simulated area should be as far away from these areas as possible. To better simulate the movement of the deep geothermal water and predict spatial and temporal changes in groundwater levels and temperature, the model area is kept away from the downtown areas and extends to the natural boundaries. Faults far from the exploitation centers of geothermal water can be chosen as the boundary of the calculation area. In the southeast, the Huairou–Laishui deep fault and intrusive rock boundary in the Xincheng area are chosen to be the southwestern boundary. In the south, the Gu’an–Changli fault (F14) is selected as the boundary. In the east, the Xianghe Fault (F15), Nankou–Sunhe Fault (F13), and Xiaqian Fault (F10) constitute the eastern boundary of the model. In the north, the plain–mountain boundary is taken as the model boundary, and it bypasses two intrusive rock areas in the plain areas in Huairou and Changping County (Figure 1).
The roof of the model (Figure 7) is determined according to the elevation of the roof of the geothermal reservoir. The bottom of the model is determined according to the thickness of the Jixian Wumishan Group. The deepest elevation of the roof of the reservoir is −4100 m in the model, and Jixian Wumishan Group is about 2000 m thick. To make the model contain the complete Wumishan Group, a thickness of 2000 m is extended downward and the elevation plane of −6100 m is chosen to be the model baseplate.

4.2. Geologic Model Structure

Elevation of the geothermal reservoir roof in the Beijing plain area is of significant difference, varying from −150 m to −4000 m. In the central part of Beijing, the roof elevation reaches the lowest level of about −4000 m. Except the Jixian Wumishan Group, other strata occur discontinuously in the horizontal direction and have a relative small thickness. Due to the discontinuous strata and the difference of roof elevation of the model, it is impossible to layer the model according to the lithology of strata in the vertical direction. The horizontal plane with the elevation of −4100 m is chosen to be an interface. There are 10 layers below −4100 m. All the 10 layers are horizontal layers, and the thickness of each layer is 200 m. The region above −4100 m is divided into 10 layers and at the same horizontal position, the thickness of each layer is the same. In this way, a 20-layered model is obtained (Figure 8). The upper 10 layers contain all the strata in the geothermal reservoir, and the elevation in each layer is different. Each layer near the central part of Beijing depression has a small thickness of about 10 m. It reaches the maximum thickness of about 400 m in the central part of the Shuangqiao subgeothermal field. The lower 10 layers contain mainly the Jixian Wumishan Group, and each layer has the same thickness of 200 m.

4.3. Boundary Conditions

Boundary conditions of the carbonate geothermal reservoir in the Beijing plain area include the groundwater flow boundary and heat migration boundary. In this simulation, the boundary of the simulation area is generalized and treated as follows.

4.3.1. Groundwater Flow Boundary

The faults to the south-east of the plain and the boundary of Beijing plain together form the boundary of the simulation area. Recharge of the deep geothermal water in the Beijing plain area comes from the lateral flow, which come from the northern and western mountainous areas where groundwater receives infiltration of precipitation. The flow direction of the geothermal water is from northwest to southeast, so the geothermal water in the Beijing plain flows out of the study area in the southeastern boundary as the form of lateral runoff. Therefore, the boundary between the mountain and plain in the north and west of the study area is set to be boundary of fixed flow, and the intrusive rock boundary strip is set to be the no-flow boundary. Since the fault boundaries in the south and east of the study area is not completely impermeable faults, they are set to be boundary of fixed flow. The flow value of the boundary is according to the Darcy’s law, and is evaluated and corrected several times. The flow values of the specific boundaries are shown in Figure 9 and Table 4.
The horizontal plane at elevation of −6000 m is taken as the bottom boundary of the model, where karst and fissure development are with poor maturity. Therefore, the bottom boundary of the model is set to be no-flow boundary. There is a weak hydraulic connection between the model roof and the overlying strata. The hot water exchange capacity is insignificant compared with the geothermal water exploitation quantity in the study area. Therefore, the model roof is also defined as the no-flow boundary.

4.3.2. Heat Migration Boundary

Geothermal survey shows that there is no significant change in temperature in the geothermal reservoir. In addition, geothermal reinjection only exerts a small influence on the temperature. Therefore, the lateral boundary of the model is set to be constant temperature boundary. The specific temperature is determined by the initial condition of the temperature field. Although the roof and baseplate of the model are no-flow boundaries, there exists heat transport from the underlying strata to the geothermal reservoir and from the reservoir to the overlying caprocks. Therefore, the upper and lower boundaries of the model are defined as heat flux boundaries. For the bottom boundary, the average heat flow value (0.065 W/m2) in North China [26] is selected as the flow value of heat flux boundary and for the top, it is 0.005 W/m2.

4.4. Initial Conditions

Extensive monitoring of the deep thermal groundwater levels began in the early 1980s [2]. In the earlier simulation of the 1970s, there were few geothermal water level data available to analyze and obtain the initial condition of the model. In the early 1970s, some geothermal wells slightly flowed thermal groundwater, and the groundwater levels were a little higher than the land surface elevation, while groundwater levels at other geothermal wells were slightly lower than the land surface. Thus, the land surface of the Beijing plain area can be considered to approximately stand for the initial groundwater surface. The deep geothermal water in the plain area flows from northwest to southeast before exploitation. The topography in Beijing is high in the northwest and low in the southeast, which is consistent with the initial distribution of thermal groundwater levels. Therefore, initial groundwater surface can be placed by the land surface approximately in this simulation. Therefore, the elevation data of the land surface are set to be initial groundwater levels. Combined with the boundary conditions in the Section 4.3.1, the natural state (no exploitation state) is simulated, and finally the approximately stable state is obtained as the initial groundwater level at the beginning of the simulation (Figure 10).
The initial temperature at the elevation of −2000 m is obtained by interpolating the temperature field data (Figure 4). Since the geothermal gradient of each strata in the thermal reservoir are nearly equal to each other, the average geothermal gradient (1.7 °C/100 m) of Cambrian System, Tieling Group, and Wumishan Group of Jixian System can be selected as the average geothermal gradient of the thermal reservoir in the study area (Table 2). The following equation is used to calculate the temperature of the other locations of the geothermal reservoir
T ( x , y , z ) = T ( x , y , 2000 ) z + 2000 100 × 1.7
where T ( x , y , z ) is the temperature at any point in the thermal reservoir and T ( x , y , 2000 ) is the temperature at the corresponding point on the plane of −2000 m. At last, the initial temperature field from layer 1 to 20 in the model are obtained (Figure 11).

4.5. Source and Sink

Thermal groundwater in the geothermal reservoir of the Beijing plain area receives recharge from infiltration of precipitation, and discharges through springs, artificial exploitation, and lateral run-off. During the simulation period, the springs in the study area is zero-flow, and in the 21st century, the thermal groundwater is artificially reinjected into the reservoir. The lateral run-off inflow from the mountainous area and the lateral runoff inflow from the southwest boundary of the study area can be simulated as the boundary of fixed flow (mentioned in the Section 4.3.1). Exploitation wells, reinjection wells, and observation wells are distributed in the study area and concentrated in the main subgeothermal fields. The location of exploitation wells, reinjection wells, and observation wells were taken to coincide as much as possible with the nodes of the discretization mesh. During the modeling period, 100 production wells were operated at withdraw rates ranging from 40 to 1000 m3/d. Daily averaged values of pumping rates were distributed over the corresponding nodes of the pumping wells. The temperature of the reinjection wells are set to be 35 °C, an average temperature of reinjection water [2]. The reinjection amount is given with the actual amount and the observation wells are set in the various positions of the subgeothermal fields.

5. Model Calibration and Results

The existing hydrogeological materials suggest that the bedrock geothermal reservoir in the Beijing plain area is a confined aquifer containing geothermal water. The numerical model describes the unsteady flow of geothermal water and heat migration occurring simultaneously in the system. It considers the influence of temperature on water viscosity and density, and the influence of concentration is ignored.

5.1. Element Subdivision in the Computational Domain and Time Discretization

The triangle algorithm is adopted to a conduct triangle mesh generation in a plane in the study area. To improve the accuracy, secondary mesh encryption is conducted in the main fault zone and near the main exploitation wells in the model. The plane grid (Figure 12) consists of 10,422 nodes in a single plane and 20,338 elements and there are 20 model layers in total.
According to long-term monitoring materials of thermal groundwater levels and well head temperature data in the study area, January 1971 to December 1999 is chosen to be the simulation period. In this period, there is just thermal groundwater exploitation without geothermal reinjection. January 2000 to December 2016 is chosen to be the verification period. In this period, geothermal exploitation and reinjection of the geothermal water are designed. The simulated time step is 30 days.

5.2. Paremeters

3D mathematic models describing groundwater flow and heat migration mainly contain such aquifer parameters as the coefficient of permeability, thermal coefficient of conductivity, porosity and specific storage. Parameter division (Figure 13) of the model domain are based on the geological characteristics of the Beijing plain area, the roof fluctuation of the geothermal reservoir, the rock types and the water yielding of different areas. The coefficient of permeability of the geothermal reservoir is calculated based on the pumping test data in the Beijing plain area. The coefficient of heat conductivity of the rocks is lacking in the study area, and it is estimated by “heat conductivity areal map of various rocks” [26] and “heat conductivity of different minerals” [27,28]. There are few data about porosity and specific storage in the study area, and these parameters are evaluated with empirical values.
The subdivision elements, initial conditions, boundary conditions, parameter partition, and initial values in the study area discussed above are used for numerical computation to get the groundwater level and temperature changes in each node in the simulation period. The groundwater level and temperature in the nodes of the observation wells are fitted with the actual measured values. Trial and error method is applied to repeatedly adjust various relevant parameters until satisfactory fitting results are obtained. Parameters used in the numerical model are listed in Table 5. The results show that the values of all the aquifer parameters are within a reasonable range. As shown in Table 5, the hydraulic conductivity in the study area ranges from 0.003 to 0.06 m/d. The values of hydraulic conductivity show a relatively low permeability in the deep geothermal reservoir in the Beijing plain area, and indicate a slow geothermal water flow and slow renewal. The reason why the permeability coefficient is small is that the entire karst geothermal reservoir is generalized as an equivalent porous media model, and the coefficients in the fractured aquifer are replaced by the equivalent parameters. Another reason may be that part of the exploitation amount in the research area has not been counted and the real exploitation amount is greater than that input in the model. In order to achieve the same decrease in depth of the groundwater level as the measured value, the parameters related to water yielding in the model are a bit less than the actual values.
In the fault zones, the fracture and karst developed well with broken bedrocks and sound water yielding, which can provide a good condition for geothermal water flow and storage. In this simulation, subdivisions 9, 10, 11, 12, 15, 29, 30, and 31 (Figure 13) are bedrock fractured zones with a relatively larger coefficient of permeability and better well yielding than other areas.

5.3. Simulation Result and Discussion

Since the 1970s, geothermal water exploitation in Beijing keeps on increasing. The geothermal water level in various subgeothermal fields dropped in varying degrees [2]. The simulated changes of thermal groundwater level are shown in Figure 14. The mean of residual error (M), root mean square deviation (RMS), and linear correlation coefficient (r) between the simulated water level and the observed water level in each geothermal field are calculated to evaluate the fitting between the simulated and observed values.
M = 1 N i = 1 N ( c a l i o b s i )
R M S = [ 1 N i = 1 N ( c a l i o b s i ) 2 ] 1 / 2
r = i = 1 N ( c a l i c a l _ _ _ _ ) ( o b s i o b s _ _ _ _ ) ( c a l i c a l _ _ _ _ ) 2 ( o b s i o b s _ _ _ _ ) 2
where N is the number of observed values; c a l i is the calculated values of groundwater level; o b s i is the observed values of groundwater level; c a l _ _ _ _ and o b s _ _ _ _ mean the average values of calculated values and observed values, respectively.
The calculated results are shown in Table 6. The values of M and RMS are less than 10% of the drop of thermal groundwater level during the whole simulation, indicating that the error between the calculated and the observed values is small. The linear correlation coefficient is close to 1, which indicates that the calculated values are highly correlated with the observed values. A satisfactory fitting is achieved; in this simulation, the model can reasonably reflect the process of geothermal water level decline.
The groundwater level in the Liangxiang subgeothermal field declines the most. In the overall simulation period (from 1971 to 2016), it drops by 94.17 m, with an annual average drop of 2.05 m. The exploitation quantity in the Dongnanchengqu subgeothermal field and the Xiaotangshan subgeothermal field is far greater than those in other subgeothermal fields. In the simulation period, their groundwater levels decrease by 80.98 m and 51.71 m, respectively, with an annual average drop of 1.76 m and 1.12 m. The exploitation of the Liangxiang subgeothermal field ranks third among the nine subgeothermal fields in the Beijing plain area. Although it cannot compete with the Dongnanchengqu subgeothermal field and the Xiaotangshan subgeothermal field, its exploitation is still greater than those of the other six subgeothermal fields. Fracture and karst development in this area is not mature, and water diversion and storage capacity are not good enough. Relatively large exploitation quantity and poor well yielding make the Liangxiang subgeothermal field have the fastest decrease in thermal groundwater level in the study area. Different drop rates of thermal groundwater levels in various positions indicate that thermal groundwater flow in the Beijing plain area changes from the native northwest–southeast flow to the cones of depression in groundwater level centering on the Liangxiang and Dongnanchengqu subgeothermal fields (Figure 15). In the simulation period, the temperature field slightly changes (Figure 16), and the temperature change in each observation point does not exceed 2 °C. Reinjection with constant temperature can just reduce the geothermal reservoir temperature surrounding the reinjection wells to a certain degree. The results show that the maximum scope of the influence of the temperature of reinjection wells is 2 km, which exerts only a slight impact on the temperature field.
According to the calculation data of the model, the mass balance of geothermal water in the past 45 years is shown in Table 7. The average recharge of geothermal water in the study area is 5021.1 m3·d−1 and the average discharge is 24,598.9 m3·d−1, leading to the negative equilibrium of 19,577.8 m3·d−1. The average exploitation amount is 21,803.4 m3·d−1, which is the main reason for the continuous decline of groundwater water level in Beijing.

5.4. Sensitivity Analysis

In the sensitivity analysis, the single parameter and flow value of boundaries are analyzed. The horizontal hydraulic conductivity (Kx and Ky), the specific storage (Ss), and flow value of boundaries are selected for uncertainty analysis. Based on the parameters in Table 5, each parameter and flow value are increased or reduced by 10% respectively, and the variation of thermal groundwater water level at the end of the simulation period is calculated (Table 8).
As shown in Table 8, when the specific storage is changed, the variation range of thermal groundwater level is the largest, indicating that the specific storage has the greatest influence on the simulation results than other parameters. The Dongnanchengqu subgeothermal field is less affected by the boundary conditions than other subgeothermal fields because it is located in the center of the simulation area and far from the simulation boundary. The effect of horizontal hydraulic conductivity on the simulation results of the Dongnanchengqu subgeothermal field and the Liangxiang subgeothermal field is larger than the boundary condition, while in the Xiaotangshan subgeothermal field and the Lisui subgeothermal field it is slightly smaller than the boundary condition.

6. Prediction and Analysis

With the increase in exploitation quantity of the thermal groundwater, excessive exploitation has become an issue in the development of geothermal resources in Beijing. After entering the 21st century, geothermal reinjection has notably increased, especially in the areas where thermal groundwater is being extensively exploited. Reinjection and exploitation are conducted at the same time to relieve the decline trend of the groundwater levels and good effect has been achieved. However, there is an increasing demand for geothermal resources, and exploitation of geothermal water will be ever-growing. Reasonable increase in exploitation and reinjection allocation must be considered now. The numerical model established above can be used to forecast the production performance of geothermal water in the thermal groundwater system.
Schemes with exploitation and reinjection are added to this simulation. A total of three prediction schemes are designed. Scheme 1 is to keep the existing exploitation intensity unchanged, and forecast the trend of groundwater level and temperature in the geothermal water system in the Beijing plain area without reinjection. Scheme 2 refers to keep the existing exploitation intensity and added reinjection intensity unchanged, and judge the influence of current exploitation conditions on the geothermal water system in the Beijing plain area. Scheme 3 involves to keep existing reinjection unchanged and adjust the exploitation quantity of the main exploitation wells in various subgeothermal fields to seek for the optimal exploitation scheme and keep the minimum drawdown of thermal groundwater level at the same time. All three schemes will last for 10 years from January 2017 to December 2026.
The boundary conditions and parameters in the three schemes are the same as those in the simulation period. The data of groundwater level and temperature at the end of the simulation period are used as the initial conditions. In Scheme 1, the exploitation quantity is the average value from 2006 to 2016 with no reinjection. In Scheme 2, a reinjection amount is added based on Scheme 1, which is the average value from 2006 to 2016. In Scheme 3, 30% of the total exploitation quantity of subgeothermal fields in the Dongnanchengqu, Xiaotangshan, and Liangxiang were cut and allocated to the rest of the six subgeothermal fields on average. The total exploitation and reinjection amount is the same as that in Scheme 2 (Table 9).
Simulation of the prediction model (Figure 17 and Figure 18) shows that in Scheme 1, the geothermal water level in the Dongnanchengqu, Xiaotangshan, and Liangxiang subgeothermal fields drop by 26.17, 20.21, and 23.51 m, respectively. In Scheme 2, 13.76, 11.78, and 13.8 m are obtained respectively, with a reduction of 12.41, 8.43, and 9.71 m than those in Scheme 1. In Scheme 2, reinjection amount and exploitation quantity in the Jixian System in this three subgeothermal fields account for 23.9%, 39.1%, and 10.9%, respectively. The results show that the reinjection can effectively slow down the descending speed of the geothermal water levels.
In Scheme 3, the geothermal water levels in this 3 subgeothermal fields drops by 11.87, 5.63, and 3.6 m, respectively, which are 1.89, 6.15, and 10.2 m less than those in Scheme 2. In the simulation of this scheme, geothermal water level in the Dongnanchengqu, Xiaotangshan, and Liangxiang subgeothermal fields drops the least of the three schemes. The cones of depression of thermal groundwater levels have a tendency to slow down. In Scheme 3, the ratio of reinjection amount to exploitation quantity in these three subgeothermal fields are 33.0%, 54.1%, and 15.6%, respectively. The results show that this exploitation scheme can effectively alleviate the speed of drop in thermal groundwater levels in these fields. Meanwhile, the total exploitation amount in Scheme 3 remains unchanged, suggesting that the other six subgeothermal fields are of geothermal potential to be further developed.

7. Conclusions

The Beijing plain area is underlain by low to medium temperature thermal groundwater of sedimentary basin type. The well head temperature of geothermal water ranges from 48 °C to 95 °C. Due to deep burial and relatively high temperature, the density of the thermal groundwater varies with temperature, leading to complicated changes in the hydraulic field and temperature field in such a geothermal system. When the system is numerically modeled, it is necessary to simultaneously describe thermal groundwater flow and heat migration. This article establishes a 3D unsteady flow model to simulate the spatial and temporal changes in thermal groundwater level and temperature under the artificial exploitation and reinjection in the study area from 1971 to 2016. Compared with common groundwater flow equation, the flow equation in this simulation contains the natural convection caused by the density variation and heat transport is mainly controlled by convection and diffusion. In the process of establishing the model, many problems have been encountered, such as there is no waterproof structure in the study area that can be treated as the simulated boundary, the great difference of roof elevation of thermal reservoir, the discontinuous strata, the lack of water level data in the initial stage of simulation, and so on. We solved these problems by choosing faults far away from the dense exploitation areas as the simulated boundaries, layering the model with the same thickness in the vertical direction and using the elevation data of the land surface as the initial groundwater levels. In this model, the elevation of −6100 m is chosen to be the model baseplate and the elevation of −4100 m water level is taken as a boundary to divide the upper and lower regions and each has 10 layers on average. The hydrological properties of different parts of the geothermal reservoir is shown by various parameters. Constant flow boundary and temperature boundary are adopted as the boundary conditions of the model. The boundary flow is calculated according to Darcy’s law, and the boundary temperature is obtained by analyzing temperature field data of the geothermal reservoir. In the numerical simulation, parameters of the geothermal reservoir are repeatedly adjusted until the geothermal water level and temperature are in reasonable agreement with the observed values. The simulation results of parameters show that the coefficient of permeability in the study area varies from 0.003 to 0.07 m/d and a larger permeability coefficient and greater heat conductivity are in the intensive structural development zones. The prediction calculations shows that (1) reinjection can effectively reduce the drop of thermal groundwater levels; (2) geothermal water level in the Dongnanchengqu, Xiaotangshan, and Liangxiang subgeothermal fields will drop by 13.76, 11.78, and 13.8 m, respectively, if the current exploitation quantity and reinjection amount are maintained (Scheme 2) for 10 years; (3) part of the exploitation can be allocated to six other subgeothermal fields with exploitation potential to slow down the cone of depression acceleration of groundwater levels. Since the 3D model has a high requirement of field observation data, which are hard to collect in practice, more elaborate numerical simulation are expected to be conducted in the future study.

Author Contributions

Conceptualization, Z.X. and X.Z.; methodology, Z.X. and X.Z.; software, Z.X., R.C. and X.Z.; validation, Z.X. and X.Z.; formal analysis, Z.X. and X.Z.; investigation, X.Z., X.Z., Y.S. Z.S. and K.H.; resources, X.Z.; data curation, Z.X. and X.Z.; writing—original draft preparation, Z.X. and X.Z.; writing—review and editing, Z.X. and X.Z.; visualization, Z.X.; supervision, X.Z.; project administration, X.Z.; funding acquisition, X.Z. and Y.S.

Funding

This work was supported by the Natural Science Foundation of Beijing (8152026) and the Fundamental Research Funds for the Central Universities of China (2652017168).

Acknowledgments

The authors want to thank Jialan Zhu, Qingxiao Zhang, and Xiaowei Shen for their help in the field work.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhang, X. Development of Geothermal Resources in China: A Review. J. Earth Sci. 2018, 29, 452–467. [Google Scholar] [CrossRef]
  2. Beijing geological engineering survey institute. Report on Dynamic Monitoring of Geothermal Resources in Beijing; Beijing geological engineering survey institute: Beijing, China, 2011. [Google Scholar]
  3. Zhou, X.; Zhao, W. Calculation of the rising value of dynamic water level in deep underground hot water drilling wells. Surv. Sci. Technol. 2000, 5, 33–35. [Google Scholar]
  4. Zhou, X.; Chen, M.; Zhao, W.; Li, M. Modeling of a deep-seated geothermal system near Tianjin, China. Ground Water 2002, 39, 443–448. [Google Scholar]
  5. Voss, C.I. A finite-element simulation model for saturated-unsaturated, fluid-density-dependent ground-water flow with energy transport or chemically-reactive single-species solute transport. Water Resour. Investig. Rep. 1984, 84, 4369. [Google Scholar]
  6. Eckstein, Y.; Maurath, G.; Ferry, R.A. Modeling the thermal evolution of an active geothermal system. J. Geodyn. 1985, 4, 149–163. [Google Scholar] [CrossRef]
  7. Zhang, Z.; Xue, Y.; Wu, J. Research on natural convection in underground hot water transport. Hydrogeol. Eng. Geol. 1995, 22, 16–117. [Google Scholar]
  8. Zhang, Z.; Xue, Y. Numerical simulation of natural convection during heat transfer in aquifer. Adv. Water Sci. 1996, 7, 99–104. [Google Scholar]
  9. Ondrak, R.; Möller, P. Modeling coupled heat and mass transport applied to the hydrothermal system of the Upper Harz Mountains (Germany). Chem. Geol. 1999, 155, 171–185. [Google Scholar] [CrossRef]
  10. Dilsiz, C. Conceptual hydrodynamic model of the Pamukkale hydrothermal field, southwestern Turkey, based on hydrochemical and isotopic data. Hydrogeol. J. 2006, 14, 562–572. [Google Scholar] [CrossRef]
  11. Xue, Y.; Xie, C. Numerical Simulation of Groundwater; Science Press: Beijing, China, 2007. [Google Scholar]
  12. Ostermann, I. Three-dimensional modeling of heat transport in deep hydrothermal reservoirs. GEM Int. J. Geomath. 2011, 2, 37–68. [Google Scholar] [CrossRef]
  13. Diersch, H.J. Heat Transport in Porous Media; Springer: Berlin/Heidelberg, Germany, 1972. [Google Scholar]
  14. Bear, J. Dynamics of Fluids in Porous Media; American Elsevier Publishing Co.: New York, NY, USA, 1972. [Google Scholar]
  15. Bureau of Geology and Mineral Resources of the Municipality of Beijing. Regional Geology of Beijing; Geological Publishing House: Beijing, China, 1991. [Google Scholar]
  16. Bureau of Geology and Mineral Resources of the Municipality of Beijing. Stratigraphy of the Municipality of Beijing; China University of Geosciences Press: Wuhan, China, 2008. [Google Scholar]
  17. Xin, Y. Discussion on the cause of formation of Beijing xiaotangshan geothermal field. Geophys. Geochem. Explor. 1989, 13, 473–475. [Google Scholar]
  18. Pan, X.; Wang, Z. Geochemical characteristics of hot water in Xiaotangshan geothermal field. Beijing Geol. 1999, 4, 7–15. [Google Scholar]
  19. Liu, J.; Pan, X.; Yang, Y.; Wang, X.; Liu, Z.; Zhang, L.; Zheng, K. Geochemistry of hot water in a geothermal well in Beijing urban area. Mod. Geol. 2002, 16, 318–321. [Google Scholar]
  20. Zhao, X. Distribution and migration characteristics of deep geothermal geology in Beijing Olympic park. Hebei Coal 2005, 6, 1–3. [Google Scholar]
  21. Lu, J.; Che, Y.; Wang, J.; Liu, Z.; Liu, C.; Zheng, G. Hydrogeochemical characteristics of hot water and genetic model of geothermal system in north Beijing area. Seismology 2006, 28, 419–429. [Google Scholar]
  22. Wang, X.; Xie, Z.; Li, S.; Cui, Y.; Shao, J. Environmental isotopic studies of groundwater near suburbs of Beijing city. Geosci. Lead. Edge 2006, 13, 205–210. [Google Scholar]
  23. Territory Resources Bureau of Beijing. Sustainable Utilization Plan of Geothermal Resources of Beijing from 2006 to 2020; Territory Resources Bureau of Beijing: Beijing, China, 2006. [Google Scholar]
  24. Diersch, H.G. Variable-Density Flow, Mass and Heat Transport in Porous Media; Springer: Berlin/Heidelberg, Germany, 2014. [Google Scholar]
  25. Chen, C. Geothermics of North China; Science Press: Beijing, China, 1988. [Google Scholar]
  26. Kappelmeyer, O.; Haenel, R. Geothermics with Special Reference to Application; Berlin Gebrueder Borntraeger Geoexploration Monographs: Berlin, Germany, 1974. [Google Scholar]
  27. Birch, F.; Clark, H. The thermal conductivity of rocks and its dependence upon temperature and composition; Part II. Am. J. Sci. 1940, 238, 529–558. [Google Scholar] [CrossRef]
  28. Beck, A.E.; Beck, J.M. On the measurement of the thermal conductivities of rocks by observations on a divided bar apparatus. Eos Trans. Am. Geophys. Union 1958, 39, 1111–1123. [Google Scholar] [CrossRef]
Figure 1. Map of the Beijing plain area, showing boundaries and faults in the study area. 1—Boundary between mountainous area and plain area in Beijing; 2—Fault; 3—Administrative boundary of Beijing; 4—Boundary of model area; F1—Huangzhuang—Gaoliying Fault; F2—Babaoshan Fault; F3—Chegongzhuang Fault; F4—Lianhuachi Fault; F5—Liangxiang–Qianmen Fault; F6—Chongwenmen Fault; F7—Nanyuan–Tongxian Fault; F8—Shunyi Fault; F9—Nangong Fault; F10—Xiaqian–Mafang Fault; F11—Yongdinghe Fault; F12—Taiyanggong Fault; F13—Nankou–Sunhe Fault; F14—Guan–Changli Fault; F15—Xianghe Fault.
Figure 1. Map of the Beijing plain area, showing boundaries and faults in the study area. 1—Boundary between mountainous area and plain area in Beijing; 2—Fault; 3—Administrative boundary of Beijing; 4—Boundary of model area; F1—Huangzhuang—Gaoliying Fault; F2—Babaoshan Fault; F3—Chegongzhuang Fault; F4—Lianhuachi Fault; F5—Liangxiang–Qianmen Fault; F6—Chongwenmen Fault; F7—Nanyuan–Tongxian Fault; F8—Shunyi Fault; F9—Nangong Fault; F10—Xiaqian–Mafang Fault; F11—Yongdinghe Fault; F12—Taiyanggong Fault; F13—Nankou–Sunhe Fault; F14—Guan–Changli Fault; F15—Xianghe Fault.
Water 11 01494 g001
Figure 2. Geological sketch map of the Beijing plain area, showing strata and their distribution. 1—Administrative boundary of Beijing; 2—Faults; 3—Boundary between mountainous area and plain area in Beijing; 4—Stratigraphic boundary; 5—Cretaceous System; 6—Jurassic System and Cretaceous System; 7—Jurassic System; 8—Carboniferous System and Permian System; 9— Ordovician System; 10—Cambrian System; 11—Qingbaikou System; 12—Jixian System; 13—Changcheng System; 14—Archean Erathem; 15—Intrusive rock (granite); 16—Intrusive rock (diorite).
Figure 2. Geological sketch map of the Beijing plain area, showing strata and their distribution. 1—Administrative boundary of Beijing; 2—Faults; 3—Boundary between mountainous area and plain area in Beijing; 4—Stratigraphic boundary; 5—Cretaceous System; 6—Jurassic System and Cretaceous System; 7—Jurassic System; 8—Carboniferous System and Permian System; 9— Ordovician System; 10—Cambrian System; 11—Qingbaikou System; 12—Jixian System; 13—Changcheng System; 14—Archean Erathem; 15—Intrusive rock (granite); 16—Intrusive rock (diorite).
Water 11 01494 g002
Figure 3. Profile showing the Beijing depression, section location is indicated in Figure 2. 1—Quaternary System; 2—Neogene System; 3—Palaeogene System; 4—Cretaceous System; 5—Jurassic System; 6—Carboniferous System and Permian System; 7—Ordovician System; 8—Cambrian System; 9—Qingbaikou System; 10—Tieling Group of Jixian System; 11—Hongshuizhuang Group of Jixian System; 12—Wumishan Group of Jixian System.
Figure 3. Profile showing the Beijing depression, section location is indicated in Figure 2. 1—Quaternary System; 2—Neogene System; 3—Palaeogene System; 4—Cretaceous System; 5—Jurassic System; 6—Carboniferous System and Permian System; 7—Ordovician System; 8—Cambrian System; 9—Qingbaikou System; 10—Tieling Group of Jixian System; 11—Hongshuizhuang Group of Jixian System; 12—Wumishan Group of Jixian System.
Water 11 01494 g003
Figure 4. Sketch profile showing formation of the thermal groundwater, section location is indicated in Figure 4. 1—Fault; 2—Rainfall; 3—Flowing direction of thermal groundwater; 4—Terrestrial heat flow; 5—Temperature contour; 6—Formation symbol; Q—Quaternary System; N—Neogene System; E—Palaeogene System; K—Cretaceous System; J—Jurassic System; C—Carboniferous System; P—Permian System; O—Ordovician System; Cm—Cambrian System; Jx—Jixian System.
Figure 4. Sketch profile showing formation of the thermal groundwater, section location is indicated in Figure 4. 1—Fault; 2—Rainfall; 3—Flowing direction of thermal groundwater; 4—Terrestrial heat flow; 5—Temperature contour; 6—Formation symbol; Q—Quaternary System; N—Neogene System; E—Palaeogene System; K—Cretaceous System; J—Jurassic System; C—Carboniferous System; P—Permian System; O—Ordovician System; Cm—Cambrian System; Jx—Jixian System.
Water 11 01494 g004
Figure 5. Contour map showing the temperature at the elevation of −2000 m in the Beijing plain area. 1—Mountainous area; 2—Main roads in Beijing; 3—Fault; 4—Temperature contour; 5—Horizontal position of A–A’ section.
Figure 5. Contour map showing the temperature at the elevation of −2000 m in the Beijing plain area. 1—Mountainous area; 2—Main roads in Beijing; 3—Fault; 4—Temperature contour; 5—Horizontal position of A–A’ section.
Water 11 01494 g005
Figure 6. Location of the subgeothermal fields in Beijing. 1—Administrative boundary of Beijing; 2—Boundary between mountainous area and plain area in Beijing; 3— Boundary of each county in Beijing; 4—Boundary of each subgeothermal field; 5—Mountainous area; 6—Plain area; 7—Administration center of each county in Beijing; 8—Yanqing subgeothermal field; 9—Jingxibei subgeothermal field; 10—Xiaotangshan subgeothermal field; 11—Houshayu subgeothermal field; 12—Lisui subgeothermal field; 13—Liangxiang subgeothermal field; 14—Tianzhu subgeothermal field; 15—Shuangqiao subgeothermal field; 16—Dongnanchengqu subgeothermal field; 17—Fengheying subgeothermal field.
Figure 6. Location of the subgeothermal fields in Beijing. 1—Administrative boundary of Beijing; 2—Boundary between mountainous area and plain area in Beijing; 3— Boundary of each county in Beijing; 4—Boundary of each subgeothermal field; 5—Mountainous area; 6—Plain area; 7—Administration center of each county in Beijing; 8—Yanqing subgeothermal field; 9—Jingxibei subgeothermal field; 10—Xiaotangshan subgeothermal field; 11—Houshayu subgeothermal field; 12—Lisui subgeothermal field; 13—Liangxiang subgeothermal field; 14—Tianzhu subgeothermal field; 15—Shuangqiao subgeothermal field; 16—Dongnanchengqu subgeothermal field; 17—Fengheying subgeothermal field.
Water 11 01494 g006
Figure 7. Contour of the elevation of the geothermal reservoir roof.
Figure 7. Contour of the elevation of the geothermal reservoir roof.
Water 11 01494 g007
Figure 8. Sections showing layers of the geologic model, section locations are indicated in Figure 7.
Figure 8. Sections showing layers of the geologic model, section locations are indicated in Figure 7.
Water 11 01494 g008
Figure 9. Locations of the fluid flux boundaries.
Figure 9. Locations of the fluid flux boundaries.
Water 11 01494 g009aWater 11 01494 g009b
Figure 10. Distribution of hydraulic heads of groundwater at the beginning of the simulated period.
Figure 10. Distribution of hydraulic heads of groundwater at the beginning of the simulated period.
Water 11 01494 g010
Figure 11. Distribution of temperature of groundwater at the beginning of the simulated period.
Figure 11. Distribution of temperature of groundwater at the beginning of the simulated period.
Water 11 01494 g011
Figure 12. Discretization in a plane of the simulation area.
Figure 12. Discretization in a plane of the simulation area.
Water 11 01494 g012
Figure 13. Parameter division of the simulation area. 1-32 represent the number of partition for parameter.
Figure 13. Parameter division of the simulation area. 1-32 represent the number of partition for parameter.
Water 11 01494 g013
Figure 14. Curves of the observed and simulated hydraulic heads during the simulated period, locations of the observed points are indicated in Figure 10.
Figure 14. Curves of the observed and simulated hydraulic heads during the simulated period, locations of the observed points are indicated in Figure 10.
Water 11 01494 g014
Figure 15. Distribution of hydraulic heads of groundwater at the end of the simulated period.
Figure 15. Distribution of hydraulic heads of groundwater at the end of the simulated period.
Water 11 01494 g015
Figure 16. Simulated temperature of different observation points in Layer 1 during the simulated period.
Figure 16. Simulated temperature of different observation points in Layer 1 during the simulated period.
Water 11 01494 g016
Figure 17. Curves of simulated hydraulic heads during the predicted period.
Figure 17. Curves of simulated hydraulic heads during the predicted period.
Water 11 01494 g017
Figure 18. Simulated temperature of different observation points in Layer 1 during the predicted period.
Figure 18. Simulated temperature of different observation points in Layer 1 during the predicted period.
Water 11 01494 g018
Table 1. Characteristics of the basement aquifers in the Beijing plain area.
Table 1. Characteristics of the basement aquifers in the Beijing plain area.
AgeFormation SymbolLithologyBurial Depth to Top (m)Thickness (m)Temperature (°C)Hydrogeologic Unit
Ordovician SystemODolomitic limestone200–200085060–78Aquifer
Cambrian SystemCmDolomitic limestone500–1500<50040–60Aquitard
Qingbaikou SystemQnShale and quartz sandstone400–1000200-500/Aquifer
Jixian SystemJxSiliceous dolostone50–4000>260038–88Aquifer
Table 2. Geotemperature gradient of different strata.
Table 2. Geotemperature gradient of different strata.
StrataFormation SymbolAverage Geothermal Gradient (°C/100 m)
Quaternary SystemQ4.18
Tertiary SystemT3.98
Jurassic SystemJ3.53
Ordovician SystemO1.85
Cambrian SystemCm1.78
Qingbaikou SystemQn3.42
Jixian System Tieling GroupJxt1.80
Jixian System Hongshuizhuang GroupJxh4.80
Jixian System Wumishan GroupJxw1.52
Table 3. Statistical table of thermal ground water exploitation in the study area (unit: 10,000 m3).
Table 3. Statistical table of thermal ground water exploitation in the study area (unit: 10,000 m3).
Year197119721973197419751976197719781979198019811982198319841985
Amount1530781051331872232332303047858909571024776
Year198619871988198919901991199219931994199519961997199819992000
Amount766827843827910844867963885837877881980921799
Year200120022003200420052006200720082009201020112012201320142015
Amount8718908818047727189009381289143614171420135013801350
Table 4. Flow value of different boundaries.
Table 4. Flow value of different boundaries.
Boundary NumberFlow Value (m/d) Boundary NumberFlow Value (m/d)
1−4.32 × 10−661.3 × 10−6
2−4.12 × 10−67−2.5 × 10−6
39.6 × 10−88−2.88 × 10−6
42.04 × 10−69−3.02 × 10−6
51.7 × 10−6
d = day.
Table 5. Values of parameters used with the model
Table 5. Values of parameters used with the model
Number of ZoneKx (m/d)Ky (m/d)Kz (m/d)Ss (1/m)φΛ (J/m/s/K)
10.0040.0040.00082 × 10−70.023.5
20.0060.0060.00122 × 10−70.023.5
30.0080.0080.00162 × 10−70.023.5
40.0140.0140.00282 × 10−70.023.5
50.0040.0040.00082 × 10−70.023.5
60.0050.0050.0012 × 10−70.023.5
70.0030.0030.00062 × 10−70.023.5
80.020.020.0042 × 10−70.023.5
90.0310.0310.0314 × 10−70.045
100.0470.0470.0474 × 10−70.045
110.0560.0560.0564 × 10−70.045
120.060.060.064 × 10−70.045
130.0070.0070.00142 × 10−70.023.5
140.0060.0060.00122 × 10−70.023.5
150.010.010.014 × 10−70.045
160.0060.0060.00122 × 10−70.023.5
170.0280.0280.00562 × 10−70.023.5
180.020.020.0042 × 10−70.023.5
190.010.010.0022 × 10−70.023.5
200.0250.0250.0052 × 10−70.023.5
210.0060.0060.00122 × 10−70.023.5
220.0030.0030.00062 × 10−70.023.5
230.0040.0040.00082 × 10−70.023.5
240.0020.0020.00042 × 10−70.023.5
250.0130.0130.00262 × 10−70.023.5
260.0080.0080.00162 × 10−70.023.5
270.0260.0260.00522 × 10−70.023.5
280.0040.0040.00082 × 10−70.023.5
290.0360.0360.0364 × 10−70.045
300.0280.0280.0284 × 10−70.045
310.0480.0480.0484 × 10−70.045
Note that Kx, Ky, Kz are coefficient of permeability of x, y, z direction; Ss is specific storage; φ is porosity; λ is thermal conductivity. J =Joule; s = second; and K = Kelvin).
Table 6. Fitting error of thermal groundwater level in simulation period.
Table 6. Fitting error of thermal groundwater level in simulation period.
Geothermal FieldM(m)RMS(m)rDrop of Thermal Groundwater Level during the Whole Simulation (m)
Dongnanchengqu−2.765.890.9880.98
Xiaotangshan−4.024.60.9851.71
Liangxiang0.875.290.8794.17
Lisui−2.183.680.9652.95
Table 7. Table of average water balance of the study area from 1971 to 2016.
Table 7. Table of average water balance of the study area from 1971 to 2016.
RechargeAmount/m3·d−1Ratio/%
Inflow of lateral boundary2920.858.17
Reinjection amount2100.341.83
Total5021.1100
Difference between recharge and discharge−19577.8
DischargeAmount/m3·d−1Ratio/%
Outflow of lateral boundary2795.511.36
Exploitation amount21803.488.64
Total24598.9100
Difference between recharge and discharge−19577.8
Table 8. Fitting error of thermal groundwater level in simulation period.
Table 8. Fitting error of thermal groundwater level in simulation period.
Subgeothermal FieldKxy ↑ 10%Kxy ↓ 10%Ss ↑ 10%Ss↓ 10%Flow value of Boundaries ↑ 10%Flow Value of Boundaries ↓ 10%
Xiaotangshan0.068−0.4024.071−5.0850.432−0.403
Dongnanchengqu2.299−3.0725.197−6.3720.0220.049
Liangxiang3.198−4.1775.308−6.5130.073−0.008
Lisui0.222−0.5434.955−6.111−0.5040.590
Table 9. Exploitation and reinjection amount under different schemes in the predicted period.
Table 9. Exploitation and reinjection amount under different schemes in the predicted period.
SchemeSubgeothermal FieldExploited StrataExploitation Amount (× 104 m3/y)Reinjection Amount (× 104 m3/y)Net Exploitation Amount (× 104 m3/y)Ratio of Reinjection Amount to Exploitation Amount
1DongnanchengquJx272/272/
XiaotangshanJx370/370/
O-Cm30/30/
LiangxiangJx110/110/
JingxibeiJx180/180/
LisuiJx77/77/
TianzhuJx95/95/
ShuangqiaoJx30/30/
HoushayuJx21/21/
FengheyingO2/2/
2DongnanchengquJx2726520723.9%
XiaotangshanJx37014522539.1%
O-Cm30/30/
LiangxiangJx110129810.9%
JingxibeiJx1803015016.7%
LisuiJx77532468.8%
TianzhuJx95207521.1%
ShuangqiaoJx3092130.0%
HoushayuJx21/21/
FengheyingO2/2/
3DongnanchengquJx1976513233.0%
XiaotangshanJx26814512354.1%
O-Cm30/30/
LiangxiangJx77126515.6%
JingxibeiJx2733024311.0%
LisuiJx117536445.3%
TianzhuJx1452012513.8%
ShuangqiaoJx4593620.0%
HoushayuJx32/26/
FengheyingO3/3/
y = year. Scheme 1: Keeping the existing exploitation without reinjection; Scheme 2: Keeping the existing exploitation and reinjection; Scheme 3: Keeping the existing reinjection and adjust the exploitation.

Share and Cite

MDPI and ACS Style

Xu, Z.; Zhou, X.; Chen, R.; Shen, Y.; Shang, Z.; Hai, K. Numerical Simulation of Deep Thermal Groundwater Exploitation in the Beijing Plain Area. Water 2019, 11, 1494. https://doi.org/10.3390/w11071494

AMA Style

Xu Z, Zhou X, Chen R, Shen Y, Shang Z, Hai K. Numerical Simulation of Deep Thermal Groundwater Exploitation in the Beijing Plain Area. Water. 2019; 11(7):1494. https://doi.org/10.3390/w11071494

Chicago/Turabian Style

Xu, Zhongping, Xun Zhou, Ruige Chen, Ye Shen, Ziqi Shang, and Kuo Hai. 2019. "Numerical Simulation of Deep Thermal Groundwater Exploitation in the Beijing Plain Area" Water 11, no. 7: 1494. https://doi.org/10.3390/w11071494

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop