Numerical Simulation and Design Optimization of Large-Scale Geothermal Production Based on a Multiwell Layout in Xianxian Geothermal Field

Yanguang Liu; Lina Guo; Yan Wang; Yuzhong Liao; Gui Zhao; Bisheng Wu

doi:10.2113/2022/2929551

Abstract

Large-scale geothermal production based on a multiwell layout has been promoted by the development of clean energy. However, the theories about the effects of this geothermal production pattern on the temperature field of geothermal reservoirs are yet to be improved, which hinders the design optimization of the deployment schemes of geothermal wells. This study proposed a 3D geothermal multiwell simulation method covering the heat transfer in wellbores by introducing the geothermal well theory of 1D line elements, thus improving the calculation efficiency without reducing the calculation accuracy. By taking the geothermal development using multiple wells in the Xianxian geothermal field, China, as an example, this study firstly analyzed the evolution of the temperature field in geothermal reservoirs during the geothermal exploitation based on a multiwell layout. Then, it revealed the influencing mechanisms of geothermal multiwell effect, reinjection depth, and well spacing on the thermal breakthrough of exploitation wells. Finally, it proposed an optimized design method for the geothermal development based on a multiwell layout. All these will provide scientific support for the large-scale geothermal production based on exploitation-reinjection combination using multiple wells.

1. Introduction

Reinjection in geothermal fields not only serves as a way to treat wastewater in the environment but also provides necessary water supply for geothermal systems to ultimately sustain geothermal exploitation. Since its introduction in the late 1960s, geothermal reinjection has been carried out in more than half of geothermal fields in the world and large-scale geothermal production based on a multiwell layout has been formed in the geothermal heating in cities in China [1–4]. However, the geothermal reinjection without proper management tends to result in adverse effects (e.g., thermal breakthrough) and other challenges affecting the sustainable development of geological reservoirs.

To seek the “optimal” reinjection strategy, early geothermal studies mainly focused on analytical and numerical approaches and productivity monitoring [5–7]. Among them, the analytic approaches generally present geothermal reservoirs as a conceptual model of homogeneous porous media interbedded with a single type of fissures or parallel fissures. For geothermal systems with complex geometric structures, numerical models help to fully understand the water-rock heat transfer mechanisms and major physical processes in the geothermal systems. Moreover, numerical simulation has attracted more attention due to its adaptability to varying geological structures, nonlinear parameters, and complex boundary conditions, especially for high-temperature steam-water systems [8]. As a result, it has gradually become a powerful tool for the assessment of geothermal resources, the design optimization of geothermal exploitation and reinjection schemes, and the thermal breakthrough prediction of geothermal wells [9]. Blöcher et al. employed the open-source software OpenGeoSys (OGS) to compare different exploitation schemes for deep geothermal systems in Groß Schönebeck in the North German Basin [10]. Zhao et al. proposed an equivalent seepage channel model to conduct the inversion of tracer tests and the thermal breakthrough prediction of exploitation wells and put forward an optimized design scheme of geothermal exploitation-reinjection systems [11]. Yin et al. carried out the trace tests under a multiwell layout and researched the effects of fractures on the runoff direction and thermal breakthrough in carbonate geothermal reservoirs [12]. Kim et al. also utilized the commercial software COMSOL in the study of heat flow mechanisms in geothermal doublet systems and predicted the capacity under different operation conditions [13]. Numerical approaches have played an important role in the production and research of deep geothermal resources [14–16].

The exploitation-reinjection strategies for HDR geothermal resources and enhanced geothermal systems have also been studied in recent years. Chen et al. established a discrete fracture model to simulate the fracture networks randomly generated and investigated the characteristics of fluid flow and heat transfer in fissure-type porous media [17]. Hawkins et al. simulated the heat transfer in fractured strata using multiple tracers. Bai et al. conducted an experimental study on the fluid flow and heat transfer in the fractures in granites and developed an analysis formula of overall heat transfer coefficient [18]. Delbar and Chapuis put forward new equations for the advective part considering the distortion of flow lines around wells in the tracer movement in a straight uniform flow [19]. Asai et al. studied various aspects related to the simulation of enhanced geothermal systems, comprehensively investigated the effects of different grid systems on the transient temperature prediction of exploited water, and proposed a practical solution [20].

In-depth studies have been conducted on the dynamic processes of pressure, temperature, hydrochemistry, and stress in geothermal reservoirs equipped with exploitation and reinjection systems. However, they mainly focus on doublet systems consisting of one exploitation well and one reinjection well (on a doublet scale) and emphasize the response regularity of geothermal doublet systems to reinjection [21, 22]. Since geothermal exploitation in geothermal fields tends to involve multiple doublet systems in engineering practices, whether these doublet systems interact with each other has been a high priority in the development and management of geothermal resources [2]. There are tens or even hundreds of geothermal wells in many geothermal fields all over the world. However, the actual modeling process of these geothermal reservoirs involving geothermal fields and multiple geothermal wells encounters the computational bottleneck induced by the significant difference in size between surrounding rocks and well structure. In detail, owing to the difference in size between geothermal wells and geothermal fields (~dm vs. ~km), it is usually required to divide the geometrical structures of the geothermal wells into very small cells, thus inducing a huge amount of computation. Most especially, hundreds of thousands or even millions of cells are required for a 3D geothermal reservoir model, which makes the computational time too long for engineering practices. To overcome such numerical difficulty, Al-Khoury et al. [23] and Saeid et al. [15] established a pseudo-3D wellbore model, in which 1D line elements were adopted to simulate the heat transfer in geothermal wells. Zhao et al. [11] and Wang et al. [24] further extended this model. At present, it is still a cutting-edge topic for the management and sustainable development and utilization of geothermal resources to establish a numerical model of a geothermal multiwell system on the scale of a geothermal field with improved calculation efficiency and satisfactory calculation accuracy.

This paper proposed a 3D geothermal multiwell simulation method covering the heat transfer in wellbores by introducing the theory of simplified 1D line elements to simplify geothermal wells. With this method, the difficulties of the grid division and calculation efficiency of the 3D model of the geothermal reservoirs on a scale of hundreds of square meters were solved. By taking the multiwell system in the Xianxian geothermal field as an example, this study conducted the numerical simulation of multiple exploitation and reinjection wells based on tracer tests and revealed the evolution of the temperature field in the geothermal field. Moreover, it compared and analyzed the effects of different exploitation and reinjection volumetric flow rates on the thermal breakthrough of exploitation wells and proposed the influencing mechanisms of the geothermal multiwell effect, the depth of the reinjection wells, and the distance between the exploitation well and reinjection well in a doublet system (i.e., the well spacing) on the temperature field in geothermal reservoirs. Based on this, it put forward an optimized method for the multiwell exploitation-reinjection design of carbonate geothermal reservoirs. All these will provide scientific support for the thermal breakthrough prediction in large-scale geothermal production.

2. Overview of the Study Area

The Xianxian geothermal field is located in Xianxian County, Hebei Province, China. It stretches across the Cangxian Uplift and the Jizhong Depression, of which the secondary tectonic units include the Raoyang Sag, the Xian Salient, the Fucheng Sag, and the Qingxian Salient (Figure 1). As the most important fault in the Xianxian geothermal field, the Xian Fault is an inherited fault and was active from the Archean to the Cenozoic. Meanwhile, it has a long extension and large fault throw, with a NE-NNE strike and a NW dip. As the boundary between two third-order tectonic units (i.e., the Cangxian Uplift and the Jizhong Depression) in the geothermal field, the Xian Fault affects the formation sedimentation. As a result, the fourth-order tectonic unit Xian Salient lacks Guantao Formation and Lower Tertiary sediments. Porous and karst-fissure-type geothermal reservoirs are developed in the Xianxian geothermal field. Among them, the bedrock geothermal reservoirs mainly consist of the Jixianian Wumishan and Gaoyuzhuang formations of the Mesoproterozoic. They are mainly distributed in the eastern part of Xianxian County, with a burial depth of top boundary of 1100−1500 m, a ratio of effective reservoir thickness to total reservoir thickness of 15%–30%, and a geothermal reservoir thickness of 120−2500 m. In terms of lithology, they are mainly composed of dolomite and flint-bearing dolomite.

A total of 58 geothermal wells are distributed in the Xianxian geothermal field. Among them, 31 geothermal wells were drilled in the bedrock reservoir area (Figure 2 and Table 1). Monitoring well GRY1, exploitation well ZK1, and reinjection well ZK2 are located in the Meizhuangwa Farm in the northeast of Xianxian County. They are more than 5 km away from other geothermal doublets, while the other geothermal doublets are relatively closer to each other, with the smallest distance of only 57 m.

3. Numerical Simulation Method

3.1. Generalization of Geothermal Well Model

All the geothermal wells in the study area are vertical. For instance, the wellbore structure of the typical geothermal well ZK1 is shown in Figure 3, in which the well parts at a depth of 0–512.0 m, 512.0–1337.88 m, and 1337.88–2500.18 m are the pump chamber section, the well wall section, and the open-hole water intake section, respectively. Surface casing and intermediate casing have been placed into the pump chamber section and the well wall section, respectively, and cement mortar is used for cementing. Given that geothermal wells are highly slender, Al-Khoury et al. [23] and Saeid et al. [15] established a pseudo-3D wellbore model using 1D line elements to simulate the heat transfer in geothermal wells. Zhao et al. [11] and Wang et al. [24] further extended this model by applying 1D line elements to simplify the 3D wellbore structure (Figure 4). On this basis, this study established a geothermal multiwell model by introducing the method of simplifying the 3D wellbore structure using 1D line elements, in which the seepage heat transfer along the axial direction of wellbores is covered and the assumptions are as follows:

(1)
The reservoirs are completely saturated with geothermal water, and the gas-liquid phase transition in the reservoirs is not considered
(2)
The cap rocks are completely dry, and the hydraulic connection between different reservoirs is not considered
(3)
The fluids in the wellbore flow along the axial direction of the wellbore, and their radial flow velocities are equal at the same depth
(4)
The correlation between properties of the fluids in the wellbore (i.e., heat conductivity coefficient and specific heat capacity) and the temperature is considered, but the gas-liquid phase transition in the wellbore is not considered

3.2. Mathematical Model of Hydrothermo Coupling

3.2.1. Heat Transfer Equation of Cap Rock

Assuming that the cap rock is completely dry and only considering its heat transfer

\begin{matrix} (1) & ρ_{r} C_{p, r} \frac{\partial T}{\partial t} - \nabla \cdot (k_{r} \nabla T) = 0, \end{matrix}

where

ρ_{r}

is the rock density (kg/m³),

C_{p, r}

is the specific heat capacity of porous media at a constant pressure (J/kg/K),

T

is the temperature (K), and

κ_{r}

is the permeability of porous media (m²).

3.2.2. Hydrothermo Coupling Equation of Geothermal Reservoirs

The saturated pressure flow in deep geothermal reservoirs can be described by Darcy’s law (covering the effects of gravity) assuming that the reservoirs are completely saturated and ignoring the gas-liquid phase transition in the reservoirs:

\begin{matrix} (2) & u_{r} = - \frac{κ_{r}}{μ} (\nabla p + ρ_{f} g \nabla z), \end{matrix}

where

u_{r}

is the Darcy velocity in porous media (m/s),

μ

is the hydrodynamic viscosity (Pa·s),

p

is the fluid pressure (Pa),

ρ_{f}

is the fluid density (kg/m³),

g

is the gravitational acceleration (m/s²), and

z

is the vertical coordinate (m).

The continuity equation of fluids can be expressed as [26]

\begin{matrix} (3) & ρ_{f} S_{r} \frac{\partial p}{\partial t} + \nabla \cdot (ρ_{f} u_{r}) = Q_{m}, \end{matrix}

where

S_{r}

is the water storage coefficient of porous media (1/Pa).

The seepage heat transfer in geothermal reservoirs can be described by the convective heat transfer equation [27, 28]:

\begin{matrix} (4) & {(ρ C_{p})}_{eff} \frac{\partial T}{\partial t} + ρ_{f} C_{p, f} u_{r} \cdot \nabla T - \nabla \cdot k_{eff} \nabla T = Q, \end{matrix}

where

C_{p, f}

is the specific heat capacity of the fluid at a constant pressure (J/kg/K),

k_{eff}

is the effective thermal conductivity of porous media (W/m/K), and

Q

is the heat source (W/m³).

3.2.3. Hydrothermo Coupling Equation of Geothermal Wells

The heat transfer between fluids and surrounding rock was covered using the equivalent heat transfer coefficient, and meanwhile, the effects of the casing and the cement mortar layer were included.

(1) Energy Balance Equation of Fluids in the Wellbores. The energy conservation equation of incompressible fluids in a 1D wellbore is as follows [29]:

\begin{matrix} (5) & ρ_{f} A C_{p, f} \frac{\partial T}{\partial t} + ρ_{f} A C_{p, f} {\bar{u}}_{w} \cdot \nabla T - \nabla A k_{f} \nabla T = f_{D} \frac{ρ_{f} A}{2 d_{i}} {|u|}^{3} + Q_{wall}, \end{matrix}

where

A

is the cross-sectional area of the wellbore (m²),

{\bar{u}}_{w}

is the axial flow velocity (m/s) of the fluids,

f_{D}

is the Darcy friction factor,

d_{i}

is the internal diameter of the wellbore (m), and

Q_{wall}

is the heat exchange between the well wall and its surrounding rock (W/m).

The first term on the right-hand side of equation () is used to describe the energy loss caused by the friction between the fluids and the well wall, whereby the Darcy friction factor

f_{D}

is related to Re and the ratio of

e

to

d_{i}

[30]:

\begin{matrix} (6) & f_{D} = f_{D} (Re, \frac{e}{d_{i}}), \\ Re = \frac{ρ_{f} {\bar{u}}_{w} d_{i}}{μ}, \end{matrix}

where Re is the Reynolds number and

e

is the surface roughness.

For laminar flow with the Reynolds number less than 2000, the Darcy friction factor is independent of the surface roughness of the wellbore and can be calculated by Stokes’ law [29]:

\begin{matrix} (7) & f_{D} = \frac{64}{Re} . \end{matrix}

For turbulent flow with the Reynolds number greater than 2000, the Darcy friction factor can be calculated using the Haaland equation [26].

\begin{matrix} (8) & \sqrt{\frac{1}{f_{D}}} = - 1.8 \log_{10} ({(\frac{e / d_{i}}{3.7})}^{1.11} + (\frac{6.9}{Re})) . \end{matrix}

For example, the Reynolds number Re of the geothermal well ZK1 with a diameter 215.9 mm is $1.6 \times 10^{5} > 2000$ since the fluid density is 1000 kg/m³, the dynamic viscosity of the fluids is 0.001 Pa·s, and the volumetric flow rate is 100 m³/h. Therefore, the fluids are in a turbulent state, and thus, $f_{D}$ should be calculated using equation (8).

(2) Heat Exchange Process between the Fluids in the Wellbores and Surrounding Rock Masses.

{(h z)}_{eff}

and

T_{ext}

yield

Q_{wall}

[25]:

\begin{matrix} (9) & Q_{wall} = {(h Z)}_{eff} (T_{ext} - T_{f}), \end{matrix}

where

{(h z)}_{eff}

is the equivalent heat transfer coefficient of the wellbore wall structure (W/m/K),

T_{ext}

is the temperature of the contact surface between a wellbore and the rock masses (K), and

T_{f}

is the temperature of the fluids in the geothermal well (K).

For the wellbore with a length of

Δ L

(Figure 4), in which the inner diameter and outer diameter of the casing are

r_{0}

and

r_{1}

⁠, respectively, and outer diameter of the cement mortar layer is

r_{2}

⁠, the heat exchange between the fluids in the well and the casing is

\begin{matrix} (10) & Q_{0} = h_{int} A_{Q} \cdot (T_{0} - T_{f}), \\ (11) & A_{Q} = ∆ L 2 π r_{0}, \end{matrix}

where

A_{Q}

is the area of the heat exchange between the fluids and the inner wall of the casing (m²),

h_{int}

is the thermal membrane resistance of the inner wall of the casing (W/(m²·K)), and

T_{0}

is the temperature of the inner surface of the casing (K).

Assuming that energy dissipation only occurs via heat exchange, the same amount of heat will enter the cement mortar layer through the outer wall of the casing:

\begin{matrix} (12) & Q_{1} = ∆ L 2 π r \cdot (- k_{p} \frac{d T}{d r}) . \end{matrix}

The above equation is reorganized and revised into the form of the integral from

r_{0}

to

r_{1}

⁠:

\begin{matrix} (13) & \int_{T_{0}}^{T_{1}} d T = - \int_{r_{0}}^{r_{1}} (\frac{Q_{1}}{∆ L 2 π k_{p}} \cdot \frac{1}{r}) d r . \end{matrix}

The calculation result of the integral is

\begin{matrix} (14) & T_{1} - T_{0} = - \frac{Q_{1}}{∆ L 2 π k_{p}} \cdot \ln (\frac{r_{1}}{r}) . \end{matrix}

It can also be reorganized as

\begin{matrix} (15) & Q_{1} = \frac{∆ L 2 π k_{p}}{\ln (r_{1} / r_{0})} \cdot (T_{1} - T_{0}) . \end{matrix}

Similarly, the energy entering the surrounding rock masses from the cement mortar layer is

\begin{matrix} (16) & Q_{2} = \frac{∆ L 2 π k_{c}}{\ln (r_{2} / r_{1})} \cdot (T_{2} - T_{1}), \end{matrix}

where

k_{p}

is the thermal conductivity of the casing (W/m/K) and

k_{s}

is the thermal conductivity of the cement (W/m/K).

Assuming

Q_{0} = Q_{1} = Q_{2} = Q_{wall} Δ L

and considering the temperatures of the contact surface between the wellbore and the rock masses are equal (i.e.,

T_{2} = T_{ext}

⁠), it can be deduced from equations (), (), and () that

\begin{matrix} (17) & Q_{wall} = \frac{2 π}{(1 / r_{0} h_{int}) + (\ln (r_{1} / r_{0}) / k_{p}) + (\ln (r_{2} / r_{1}) / k_{c})} (T_{ext} - T_{f}) . \end{matrix}

That is,

\begin{matrix} (18) & {(h Z)}_{eff} = \frac{2 π}{(1 / r_{0} h_{int}) + (\ln (r_{1} / r_{0}) / k_{p}) + (\ln (r_{2} / r_{1}) / k_{c})} . \end{matrix}

The thermal membrane resistance

h_{int}

induced by fluid flow in the wellbore can be calculated by

\begin{matrix} (19) & h_{int} = Nu \frac{k_{p}}{d_{i}}, \end{matrix}

where Nu is the Nussle number, which is a constant for laminar flow:

\begin{matrix} (20) & Nu = 3.66 . \end{matrix}

In the case that the fluids inside the wellbore are in a turbulent state (⁠

2000 < Re < 6 \times 10^{6}

⁠,

0.5 < \Pr < 2000

⁠), the Nussle number can be obtained by

\begin{matrix} (21) & Nu = \frac{(f_{D} / 8) (Re - 1000) \Pr}{1 - 12.7 {(f_{D} / 8)}^{1 / 2} (\Pr^{2 / 3} - 1)}, \end{matrix}

where Pr is the Prandtl number:

\begin{matrix} (22) & \Pr = \frac{C_{p, f} μ}{k_{f}} . \end{matrix}

As indicated by the above equations, the main controlling equations (5) and (9) of the heat transfer in the wellbores are retained, although the dimensions of the wellbores are geometrically simplified. Furthermore, the material parameters related to the wellbores and fluids are considered, such as the cross-sectional area of the wellbores, the thermal conductivity of the cement mortar and casing, the thermodynamic properties of the fluids, the volumetric flow rate, and the properties of the contact surface between the wellbores and surrounding rock masses. All of them are incorporated into the equivalent heat transfer coefficient.

In this study, a 3D numerical model of the geothermal reservoir engineering on a scale of well-geothermal field was established based on the finite element platform COMSOL. The numerical model was partitioned using tetrahedron elements, and the wellbores were denoted by 1D line elements. Meanwhile, interface Parallel Direct Sparse Solver in the COMSOL was used to solve the controlling equations (2), (3), and (4) of the hydrothermo coupling in the rock masses of the geothermal reservoirs, the controlling equation (1) of the heat transfer in cap rock, and the controlling equations (5) and (9) of the seepage heat transfer in wellbores. The calculation procedure of the model is shown in Figure 5.

3.3. Numerical Modeling

The region simulated using the numerical model covers an area of about 210 km². It is bounded by the Xian Fault in the west, the Shaohongdian-Gaoguan Fault in the north, the geological boundary of the bedrock in the east, and the administrative boundary of the study area in the north. Its vertical depth varies from 10 m to 4000 m underground, and it can be divided into the Quaternary System, the Neogene Minghuazhen Formation, and the Jixianian Wumishan Formation according to lithological differences. The numerical model was established based on the isolines of bedrock burial depth.

In this study, the grid cell size of the areas surrounding the geothermal wells should not exceed 10 times the diameter of the geothermal wells in principle. That is, only the areas surrounding the geothermal wells were partitioned into small cells, while the areas far away from the geothermal wells were partitioned into large cells, which contributed to the reduction in the grid cell number. In this way, the computational time can be saved. 1D line elements were employed to characterize the 31 geothermal wells listed in Table 1. The grids of the geothermal wells and their surrounding areas were partitioned into small cells, with the maximum cell size of 20 m, while the maximum cell sizes of the grids of the geothermal reservoir area and remaining cap rock were 800 m and 1000 m, respectively. As a result, the numerical model contained a total of 682,284 tetrahedral cells and 4858 edge cells (Figure 6). Parameter sensitivity analysis was conducted for the density of the grid division using the verification method of a single-well geothermal system proposed by Liu (2019). As verified by the analytical results, it is reasonable to geometrically simplify the 3D geothermal wells into 1D line elements in this study [27].

3.4. Initial and Boundary Conditions

The pressure conditions of the model were set as follows. The initial water pressure of the reservoirs was set at the measured hydrostatic pressure. It was assumed that no flow occurred on the top and bottom of the numerical model, and the sides of the model were set as pressure boundaries. The initial water pressure in the wellbores was also set at the measured hydrostatic pressure, and the flow boundary was set as the wellheads after the model started to operate (⁠ $t > 0$ ⁠).

The initial temperature field of the model was set according to the geothermal gradient (⁠

∆ T

⁠). The temperature of the model’s top was set at the fixed ground temperature (⁠

T_{s}

⁠), the model’s bottom was assumed to be of thermal insulation, and the model’s sides were set as open boundaries. The initial temperature in the wellbores was set at the initial temperature of the surrounding rocks, and the wellhead temperature was fixed at

T_{in}

after the model started to operate (⁠

t > 0

⁠).

\begin{matrix} (23) & T = T_{s} + ∆ T \cdot z, \\ T_{in} (z = 0, t) = T_{in}, \end{matrix}

where

T_{s}

is the ground temperature (K),

∆ T

is the geothermal gradient (K/m), and

T_{in}

is the reinjection temperature (K).

The heat flux at the wellheads of the exploitation wells was set as

\begin{matrix} (24) & - k_{f} \frac{d T_{p}}{d z} = 0 . \end{matrix}

Mass flux

M_{t}

was applied to the inner boundaries of the reservoirs at the interfaces between the wellbores and the surrounding rock of the reservoirs

\begin{matrix} (25) & M_{t} = \frac{ρ_{f} q}{l_{open}}, \end{matrix}

where

M_{t}

is the mass flux (kg/m/s),

l_{open}

is the length of the open-hole section of a reinjection well (m), and

q

is the volumetric flow rate (m³/h).

The heat source of the internal boundaries of the reservoirs was set as

\begin{matrix} (26) & Q = C_{p, f} M_{t} (T_{bottom} - T), \end{matrix}

where

T_{bottom}

is the temperature of the open-hole section of a reinjection well (K).

The temperature at the bottom of an exploitation well was set at the average temperature of the well’s open-hole section:

\begin{matrix} (27) & T_{f} (z = H, t) = \frac{1}{l_{open}} \int T (z) d z . \end{matrix}

3.5. Model Recognition and Verification

3.5.1. Temperature Fitting

The geothermal gradients of the strata in the study area vary with lithology. They tend to decrease with an increase in depth under the influence of the factors such as convection [31]. According to the temperature curves of wells GRY1, XXZK1, and XXZk2 in Figure 7, the geothermal gradient was 4.27°C/100 m in the Quaternary (Q) and the Minghuazhen Formation (Nm), 0.30°C/100 m in the upper part of the Wumishan Formation (Jxw), and 1.46°C/100 m from the middle part of the Wumishan Formation to the upper part of the Yangzhuang Formation due to the relatively poor water yield property. Then, the temperature decreased downward and the geothermal gradient was -1.69°C/100 m in the lower part of the Yangzhuang Formation (Jxy) due to the strong aquifers. Moreover, the geothermal gradient increased to 2.45°C/100 m in the Gaoyuzhuang Formation (Jxg) due to the poor water yield property. Besides, the geothermal gradient decreased with a rise in water-bearing capability in the area with karst developing. The initial temperature distribution in the model was simulated (Figure 8) according to the temperature measurement of well GRY1 and the parameters of the model (Table 2).

3.5.2. Fitting of the Burial Depth of Water Level

Seepage was verified through the exploitation-reinjection test of the doublet system consisting of wells ZK1 and ZK2. During the test, the burial depth of the water level was monitored from September 20, 2017, to January 9, 2018, lasting 121 days. Since November 24, 2017, the exploitation volumetric flow rate had been increased to 130 m³/h from 100 m³/h of the first stage and had been completely reinjected from well ZK2. It can be seen in Figure 9 that the decreased amplitude of the water level increased with an increase in the exploitation volumetric flow rate. Given that some data greatly fluctuated due to the effects of pump shutdown, central heating, and pump replacement, these data were removed through data processing. The simulated water level was fitted by adjusting the parameters of the numerical model, such as boundary water level and hydrogeological parameters of geothermal reservoirs (Table 2). By comparison between the simulated and measured water level (Figure 10), the parameters and boundary conditions of the numerical model were verified. Therefore, the model can be used for numerical simulation.

4. Analysis of Results

Given the actual reinjection conditions and the requirements of geothermal development and management in the study area, this study investigated the effects of reinjection rate on the sustainable development performance of geothermal reservoirs under two operation conditions with an average single-well exploitation volumetric flow rate of 100 m³/h. Under the first operation condition, incomplete reinjection using low-temperature geothermal tailwater was adopted and the reinjection rate was 80%, which was the achievable comprehensive reinjection rate of the doublet systems in the study area. Since the pressure of geothermal reservoirs will decrease in the case of insufficient natural recharge, the reinjection rate of 100% should be progressively achieved according to local management requirements for geothermal development. Under the second operation condition, the reinjection rate was 100%. According to the industry requirements, the thermal breakthrough time of an exploitation well was defined as the time when the temperature of the exploitation well decreases by no more than 2°C within 100 years of geothermal exploitation and reinjection.

4.1. The First Operation Condition

The computation under this condition was finished in about 11 hours using the multiwell model. Figure 11 shows the temperature field evolution in the geothermal reservoirs during the 100 years of geothermal exploitation and reinjection. According to Figure 12, the low-temperature zone gradually dispersed outwards and the cold fronts moved from the reinjection wells toward the exploitation wells as the exploitation and reinjection proceeded. As can be seen from Figures 12 and 13, the doublet systems in the study area are not evenly distributed, and seven of them are relatively concentrated. Figure 13 shows that the thermal breakthrough occurred in well P15 in less than one year. This was mainly caused by the unreasonable spacing between well P15 as the exploitation well and well R16 as the reinjection well—only 56 m. The spacing between wells P2 and R3 and wells P34 and R35 is 325 m and 365 m, respectively, and the thermal breakthrough time was about the 30th year and the 36th year, respectively. The remaining exploitation wells are located comparatively far away from their reinjection wells, and no thermal breakthrough occurred in them during the 100 years of operation. Therefore, it is significant for the prevention of thermal breakthrough to appropriately increase well spacing, which is roughly consistent with the research results of exploitation-reinjection systems obtained by Wang et al. [31].

4.2. The Second Operation Condition

Figure 14 illustrates the temperature field evolution in the geothermal reservoirs under the second operation condition. Owing to the increase in reinjection volumetric flow rate, the second operation condition featured quicker migration of cold fronts, earlier thermal breakthrough of some exploitation wells, higher mass flow and corresponding higher heat generation power, and more heat taken away compared with the first operation condition. All these resulted in a quick decrease in the reservoir temperature [32]. Thermal breakthrough occurred in wells P2 and P34 in about the 32nd year and the 46th year, respectively. In contrast, no thermal breakthrough occurred in the remaining exploitation wells with large well spacing, but all of them experienced large amplitude of temperature decrease (Figures 15 and 16). Among them, the exploitation well P48 and the corresponding reinjection well R49 are located 543 m away from each other, and the temperature of well P48 decreased by 1.2°C after 100 years of operation. This indicates that it is feasible to set the well spacing at greater than 500 m under the current operation condition, which is consistent with the geothermal doublet system of Liu et al. [30].

5. Influencing Factors of Temperature Field in Multiwell Geothermal Production

5.1. Geothermal Multiwell Effect

Multiple geothermal doublet systems will interact if they are close enough, which is commonly referred to as the geothermal multiwell effect [27]. A geothermal doublet model generally only covers the effects of the reinjection on the thermal breakthrough of the exploitation well within a doublet system, while it tends to ignore the effects of reinjection on other geothermal wells nearby. As a result, the predicted thermal breakthrough time of an exploitation well may disagree with its actual thermal breakthrough time [17]. Therefore, the research on the geothermal multiwell effect is greatly significant for large-scale geothermal production based on a multiwell layout. There is a doublet system consisting of exploitation well P34 and reinjection well R35 in the study area, and well P34 is only 484 m away from the reinjection well R26 in another doublet system to the west of well P34. As can be seen from Figure 17, a portion of the low-temperature water reinjected from well R26 flowed toward well P34, which inevitably imposed potential effects on the temperature of well P34. Besides, as for the doublet system composed of exploitation well P15 and reinjection well R16, the east side of well P15 is only 672 m away from reinjection well R45 of another doublet system. As shown in Figure 18, a portion of the cold water reinjected from well R45 flowed toward well P15 during geothermal exploitation and reinjection. As indicated by the hydraulic connection between the above doublet systems, the geothermal doublet systems are not kept separate from each other but are hydraulically connected to a certain extent. Therefore, the flow and transport of fluids in the geothermal reservoirs might be driven by pressure gradient and the differences in permeability between different lithologic units. Meanwhile, the temperature of the exploitation well in a doublet system may be affected by the superposition of temperature drop induced by the cold-water reinjection of the doublet systems nearby. However, the superposition effect tends to be ignored in geothermal doublet models, which may result in a low estimate of the thermal breakthrough time of exploitation wells, and thus, the geothermal development and utilization conditions in surrounding areas and their impacts should be comprehensively considered in practices [27]. The multiwell exploitation-reinjection system model proposed in this study covers the interaction among various geothermal doublet systems and thus can more effectively reflect the practical problems in engineering. Therefore, it will provide important scientific support for the optimization of the well deployment in the multiwell geothermal exploitation-reinjection mode consisting of multiple doublet systems.

5.2. Effects of Reinjection Depth on Exploitation Temperature

Previous studies indicate that reinjection depth produces significant effects on the temperature change of geothermal reservoirs. In general, reinjection wells are as deep as or slightly deeper than their corresponding exploitation wells. According to the thermal breakthrough curves of the wellheads of exploitation wells after the reinjection depth increased (Figure 19), the exploitation depth was roughly equal to the reinjection depth for wells P2 and P34 experiencing thermal breakthrough in the Xianxian geothermal field. Herein is the analysis of the sensitivity of the two geothermal wells to the reinjection depth. The depth of their corresponding reinjection wells was both increased by 500 m in the simulation, and the results are shown in Figure 19. It can be seen that thermal breakthrough of wells P2 and P34 occurred in the 49th year and the 64th, respectively, which was postponed for 17 years and 18 years, respectively, compared with the case that the exploitation depth was equal to the reinjection depth. Therefore, the thermal breakthrough was effectively postponed when reinjection wells were deployed more deeply than their corresponding exploitation wells. The main reasons include that deep reinjection wells increased the diagonal distance between the exploitation well bottom and the reinjection well bottom, increased the circulation path length and residence time of the reinjected water in the geothermal reservoirs, and enlarged the contact area between reinjected water and the geothermal reservoirs. As a result, the geothermal water absorbed more heat from the matrix of rock masses, thus postponing the decrease in the temperature of the exploitation wells. Therefore, it is suggested to appropriately deploy the reinjection wells more deeply than the exploitation wells in practical geothermal exploitation and reinjection engineering on the premise of comprehensively considering the overall layout of exploitation and reinjection wells, the technical difficulties in well drilling, and the cost of well completion. This will be greatly significant for sustainable geothermal production.

5.3. Effects of Well Spacing on Exploitation Temperature

The geothermal development and utilization in Xianxian County have a history of more than 30 years. Due to the lack of experience in early geothermal development, some geothermal wells were forced to be shut down as a result of too early thermal breakthrough, which seriously hindered the development and utilization of geothermal resources. With the development of geothermal science, the deployment of the geothermal wells was renovated and doublet systems consisting of an exploitation well and a reinjection well were established. A total of 13 doublet systems have been built up to now. In this way, a multiwell exploitation-reinjection model has been preliminarily formed, with an exploitation volumetric flow rate of about 100 m³/h and a tailwater temperature of about 25°C. All the geothermal tailwater was reinjected. This study analyzed the sensitivity of the exploitation temperature to well spacing based on the actual arrangement of the geothermal wells in Xianxian County, aiming at providing scientific support for the design of the geothermal exploitation and reinjection wells. To this end, the well spacing was set at 300 m, 400 m, and 500 m by referring to previous analytical results [31].

According to Figure 20, with continuous reinjection of the low-temperature geothermal tailwater, low-temperature zones were formed nearby the reinjection wells and cold fronts spread outwards and gradually moved toward the exploitation wells. Furthermore, the arrival time of cold fronts at the exploitation wells was related to well spacing to a certain extent [32]. As indicated by the thermal breakthrough curves of the wellheads of the exploitation wells in the case of different well spacing (Figure 21), the thermal breakthrough time of wells P1, P6, P7, P11, and P22 all correlated closely with the well spacing. In detail, the increase in the well spacing contributed to a large area and long duration of the heat exchange between the reinjected water and the matrix of rock masses. Consequently, the reinjected water absorbed heat from the rock matrix and was then heated up, thus postponing the thermal breakthrough of the exploitation wells. In the case that the mass flow of the exploitation and reinjection wells remained unchanged, the heat extraction ratio decreased with a decrease in the wellhead temperature. However, the heat extraction ratio slightly decreased during the 100 years of geothermal exploitation under the condition of large well spacing for geothermal development, which will not affect the geothermal development and utilization means. According to Table 3, the thermal breakthrough time of the five exploitation wells was less than 50 years in the case that the well spacing was 300 m. This indicates that the well spacing should not be less than 300 m; otherwise, the thermal breakthrough will occur too early. The thermal breakthrough time of well P1 was about 90 years in the case that well spacing was 500 m, and it is suggested to increase the depth of the reinjection wells to postpone the thermal breakthrough of the well. No thermal breakthrough occurred in the other four wells in 100 years. Overall, the increase in well spacing help to postpone the thermal breakthrough of exploitation wells, and it is suggested that the well spacing for geothermal development should not be less than 500 m under the conditions in this study.

6. Conclusions and Outlook

6.1. Conclusions

This paper proposed a 3D geothermal multiwell simulation method covering the heat transfer in wellbores by introducing the theory of simplified 1D line elements and using the heat transfer coefficient method. By taking the multiwell system in the Xianxian geothermal field as an example, this study conducted the numerical simulation of multiple exploitation and reinjection wells and optimized the exploitation-reinjection scheme based on a multiwell layout of carbonate geothermal reservoir. In this manner, the difficulties of the grid division and calculation efficiency of the 3D model of the geothermal reservoirs on a scale of hundreds of square meters were solved, which will provide scientific support for the thermal breakthrough prediction in large-scale geothermal production.

(1)
The 3D geothermal multiwell simulation method based on the theory of 1D line elements was verified using the geothermal multiwell system in the Xianxian geothermal field. As indicated by the verification results, the calculation efficiency was improved without reducing the calculation accuracy. Therefore, this method can serve as a new modeling tool for the design optimization, prediction, assessment, and management of geothermal multiwell systems
(2)
The evolution and influencing mechanisms of the temperature field in the geothermal reservoirs in the Xianxian geothermal field were investigated using the multiwell simulation method covering the heat transfer between geothermal wells and geothermal reservoirs. In a multiwell exploitation and reinjection system consisting of multiple doublet systems, different doublet systems will interact and thus notable multiwell effect will occur. The major factors affecting the thermal breakthrough time of exploitation wells include the exploitation volumetric flow rate, reinjection volumetric flow rate, and well spacing. An ideal well arrangement scheme for the study area should consist of the exploitation volumetric flow rate of 100 m³/h, reinjection volumetric flow rate of 100 m³/h, and well spacing of 500 m. In this case, no thermal breakthrough will occur in general

6.2. Outlook

Although the simulation method proposed in the study can be used to calculate the practical condition of geothermal development using multiple wells, it is yet to cover the factors such as the mechanical deformation of the rock masses in the matrix. Besides, parameter acquisition is to be further studied given the complex geological conditions of geothermal fields, and the details are as follows.

(1)
It is suggested that the temperature-dependent factors of geothermal reinjected fluids should be included in the simulation of geothermal multiwell systems, such as the density, specific heat capacity, dynamic viscosity, and thermal conductivity coefficient. The purpose is to further improve the simulation results of the model
(2)
The geothermal development based on a multiwell layout may change the geological conditions of geothermal resources, such as hydraulic conditions and the temperature of geothermal reservoirs around exploitation and reinjection wells. To determine whether these changes will produce the impacts such as ground deformation, it is necessary to further strengthen the research on the thermo-hydro-mechanical multifield coupling effect in geothermal development
(3)
This study established a numerical model for the geothermal development based on multiwell layout, in which equivalent seepage channels were employed by simplifying geothermal wells using 1D line elements. This model is to be verified by more examples. Most especially, further in-depth studies are required on the effects of secondary porosity features (i.e., faults, fractures, and fissures) in strata on the hydraulic migration in this model for carbonate geothermal reservoirs with relatively developed karst

Data Availability

Some or all data, models, or code generated or used during the study are available from the corresponding author by request (list items).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the S&T Program of Hebei Province, China (No. 20374201D), the foundation of Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences (Nos. SK202104 and SK202008), the Natural Science Foundation of Hebei Province, China (No. D2019330003), the National Key R&D Program of China (Nos. 2018YFC0604300 and 2018YFC0604306), and the Geological Survey Project of China (No. DD20190128).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).

Well number	Use	Depth (m)	Length of the open-hole section (m)
P1	Exploitation well	1751	384
P2	Exploitation well	1700	319
P4	Exploitation well	1698	315
P6	Exploitation well	2260	746
P7	Exploitation well	1853	510
P11	Exploitation well	2050	610
P15	Exploitation well	1750	236
P22	Exploitation well	1981	803
P25	Exploitation well	1848	326
P34	Exploitation well	1767	270
P46	Exploitation well	1697	233
P48	Exploitation well	1870	284
zk1	Exploitation well	2500	1162
R3	Reinjection well	1858	553
R5	Reinjection well	2001	720
R8	Reinjection well	2196	749
R10	Reinjection well	2050	777
R14	Reinjection well	2162	812
R16	Reinjection well	1750	378
R21	Reinjection well	2065	820
R26	Reinjection well	1900	524
R28	Reinjection well	2100	838
R31	Reinjection well	2150	740
R33	Reinjection well	2000	589
R35	Reinjection well	2004	664
R37	Reinjection well	1827	436
R45	Reinjection well	1597	251
R47	Reinjection well	1941	641
R49	Reinjection well	2140	729
R55	Reinjection well	2100	723
zk2	Reinjection well	2004	657
GRY1	Monitoring well	4000	250

Well number	Use	Depth (m)	Length of the open-hole section (m)
P1	Exploitation well	1751	384
P2	Exploitation well	1700	319
P4	Exploitation well	1698	315
P6	Exploitation well	2260	746
P7	Exploitation well	1853	510
P11	Exploitation well	2050	610
P15	Exploitation well	1750	236
P22	Exploitation well	1981	803
P25	Exploitation well	1848	326
P34	Exploitation well	1767	270
P46	Exploitation well	1697	233
P48	Exploitation well	1870	284
zk1	Exploitation well	2500	1162
R3	Reinjection well	1858	553
R5	Reinjection well	2001	720
R8	Reinjection well	2196	749
R10	Reinjection well	2050	777
R14	Reinjection well	2162	812
R16	Reinjection well	1750	378
R21	Reinjection well	2065	820
R26	Reinjection well	1900	524
R28	Reinjection well	2100	838
R31	Reinjection well	2150	740
R33	Reinjection well	2000	589
R35	Reinjection well	2004	664
R37	Reinjection well	1827	436
R45	Reinjection well	1597	251
R47	Reinjection well	1941	641
R49	Reinjection well	2140	729
R55	Reinjection well	2100	723
zk2	Reinjection well	2004	657
GRY1	Monitoring well	4000	250

Type	Name	Symbol	Value	Unit
Geothermal wells	Exploitation volumetric flow rate	$q_{out}$	100	m³/h
	Reinjection volumetric flow rate	$q_{in}$	100	m³/h
	Diameter of geothermal wells	$d_{i}$	215.9	mm
	Reinjection temperature	$T_{in}$	30	°C
	Ground temperature	$T_{s}$	20	°C
	Thickness of casing wall	$Δ w_{1}$	8.05	mm
	Thermal conductivity of casing	$k_{p}$	45	W/(m·K)
	Thickness of cement mortar wall	$Δ w_{2}$	25.55	mm
	Thermal conductivity of cement mortar	$k_{s}$	1.28	W/(m·K)

Changcheng Formation (Chg)	Density	$ρ_{r}$	2900	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	920	J/(kg·K)
	Thermal conductivity	$k_{r}$	5	W/(m·K)
	Porosity	$φ_{r}$	0.25
	Permeability	$κ_{r}$	0.631	m/d

Yangzhuang Formation (Jxy)	Density	$ρ_{r}$	2800	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	920	J/(kg·K)
	Thermal conductivity	$k_{r}$	5	W/(m·K)
	Porosity	$φ_{r}$	0.25
	Permeability	$κ_{r}$	0.631	m/d

Wumishan Formation (Jxw)	Density	$ρ_{r}$	2800	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	920	J/(kg·K)
	Thermal conductivity	$k_{r}$	5	W/(m·K)
	Porosity	$φ_{r}$	0.25
	Permeability	$κ_{r}$	1.95	m/d

Minghua Formation (Nm)	Density	$ρ_{r}$	2500	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	860	J/(kg·K)
	Thermal conductivity	$k_{r}$	2.3	W/(m·K)
Quaternary Formation (Q)	Density	$ρ_{r}$	2300	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	850	J/(kg·K)
	Thermal conductivity	$k_{r}$	2.5	W/(m·K)

Water	Density	$ρ_{f}$	1000	kg/m3
	Thermal conductivity	$k_{f}$	0.68	W/(m·K)
	Specific heat capacity under constant pressure	$C_{p, f}$	4200	J/(kg·K)
	Dynamic viscosity	$μ$	0.001	Pa·s
	Specific heat rate	$γ$	1

Type	Name	Symbol	Value	Unit
Geothermal wells	Exploitation volumetric flow rate	$q_{out}$	100	m³/h
	Reinjection volumetric flow rate	$q_{in}$	100	m³/h
	Diameter of geothermal wells	$d_{i}$	215.9	mm
	Reinjection temperature	$T_{in}$	30	°C
	Ground temperature	$T_{s}$	20	°C
	Thickness of casing wall	$Δ w_{1}$	8.05	mm
	Thermal conductivity of casing	$k_{p}$	45	W/(m·K)
	Thickness of cement mortar wall	$Δ w_{2}$	25.55	mm
	Thermal conductivity of cement mortar	$k_{s}$	1.28	W/(m·K)

Changcheng Formation (Chg)	Density	$ρ_{r}$	2900	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	920	J/(kg·K)
	Thermal conductivity	$k_{r}$	5	W/(m·K)
	Porosity	$φ_{r}$	0.25
	Permeability	$κ_{r}$	0.631	m/d

Yangzhuang Formation (Jxy)	Density	$ρ_{r}$	2800	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	920	J/(kg·K)
	Thermal conductivity	$k_{r}$	5	W/(m·K)
	Porosity	$φ_{r}$	0.25
	Permeability	$κ_{r}$	0.631	m/d

Wumishan Formation (Jxw)	Density	$ρ_{r}$	2800	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	920	J/(kg·K)
	Thermal conductivity	$k_{r}$	5	W/(m·K)
	Porosity	$φ_{r}$	0.25
	Permeability	$κ_{r}$	1.95	m/d

Minghua Formation (Nm)	Density	$ρ_{r}$	2500	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	860	J/(kg·K)
	Thermal conductivity	$k_{r}$	2.3	W/(m·K)
Quaternary Formation (Q)	Density	$ρ_{r}$	2300	kg/m³
	Specific heat capacity under constant pressure	$C_{p, r}$	850	J/(kg·K)
	Thermal conductivity	$k_{r}$	2.5	W/(m·K)

Water	Density	$ρ_{f}$	1000	kg/m3
	Thermal conductivity	$k_{f}$	0.68	W/(m·K)
	Specific heat capacity under constant pressure	$C_{p, f}$	4200	J/(kg·K)
	Dynamic viscosity	$μ$	0.001	Pa·s
	Specific heat rate	$γ$	1

Numerical Simulation and Design Optimization of Large-Scale Geothermal Production Based on a Multiwell Layout in Xianxian Geothermal Field

Abstract

1. Introduction

2. Overview of the Study Area

3. Numerical Simulation Method

3.1. Generalization of Geothermal Well Model