1 Introduction

Hot dry rock (HDR), as one of geothermal energy, has high exploitation valuableness and prospect due to its high quality of thermal storage. Enhanced geothermal systems or engineering geothermal systems (EGS) are the exploitation technology of HDR (shown in Fig. 1). Its fundament thought is to obtain high-temperature water or steam by injecting cold water into the artificial heat storage system formed by the hydraulic fracturing technology. The thermal extraction processes mainly include the following four steps. ① The low-temperature working fluid is injected into the deep storage system through the injection well; ② After fully exchanging heat between the working fluid and high-temperature rock, the temperature of working fluid increases and flows to the bottom of production well; ③ The high-temperature working fluid is extracted through the production well; ④ The working fluid is transported to the ground power generation system. Due to the complexity of heat recovery process, the exploitation of EGS involves the complex spatial-temporal evolution processes of multi-physics fields such as fluid flow, convective heat transfer and mechanical deformation in fractured rock mass [1, 2]. Numerical simulation is considered one of the effective ways to study this process [3,4,5,6,7,8,9,10,11]. Because the artificial thermal reservoir contains the fracture network caused by the hydraulic fracturing, the difficulty of numerical simulation lies in how to characterize these fractures and establish efficient and reliable multi-physics field coupling model.

Fig. 1.
figure 1

Schematic of the HDR exploitation process

Full size image

The characterization of fracture network has made a lot of progress in the simulation of petroleum exploitation. A widely used fracture model is the discrete fracture model (DFM) [12, 13]. Based on the DFM, researchers have established some geothermal exploitation models. Sun et al. [14] and Yao et al. [15] used COMSOL Multiphysics to establish a two-dimensional (2D) and three-dimensional (3D) thermal-hydraulic-mechanical (THM) coupling models of the HDR exploitation process, respectively. However, the DFM has a significant disadvantage when meshing. The fracture must be aligned with the matrix grid, and a large number of unstructured grids must be used for the simulation of complex fracture distribution, which seriously affects the efficiency of numerical calculation. This is the reason why the DFM cannot be directly applied to the actual engineering applications. To overcome the shortcomings of DFM, Lee et al. [16] developed an embedded discrete fracture model (EDFM) (shown in Fig. 2), which divides the computational domain into two kinds of media, termed matrix and fracture. The matrix can be divided by the structured grids, and the fractures are directly embedded in the grid system of the matrix, thus avoiding the complicated unstructured grids of the traditional DFM. Based on the EDFM, Gunnar et al. [2] established an efficient method to investigate the role of thermal stress during hydraulic simulation over short and long periods. Yan et al. [17] developed a hydro-mechanical coupling model of the fractured rock mass. The results of the numerical calculation are not much different from those of COMSOL Multiphysics, and the computational speed is faster. Karvounis [18] improved the EDFM and proposed an adaptive hierarchical fracture model to simulate the operational process of EGS. In EDFM, the dimension of fracture is reduced. For the 2D problem, the fracture is regarded as a 1D line segment, while for the 3D problem, the fracture is considered as a 2D plane.

Fig. 2.
figure 2

Schematic of the embedded discrete fracture model (EDFM)

Full size image

Therefore, in view of the advantages of the EDFM, the exploitation model of HDR is established by using this fracture model in this paper. Two sets of flow and energy conservation equations are used to describe the flow and heat transfer processes in the matrix and the fracture, respectively. The accuracy of the proposed model is verified by comparing with COMSOL Multiphysics. Finally, combined with fractal theory, an artificial thermal storage system with a fractal tree is generated to characterize the artificial fracture network in the exploitation process of HDR. The spatial-temporal evolution processes of pressure and temperature fields are analyzed.

2 Governing Equations

2.1 Seepage Equations

As shown in Fig. 2, the EDFM regards the fractured rock mass as two sets of independent media. The coupling between the matrix and the fracture is implemented through the transfer flow function. The mass conservation equation of the matrix media is as follow.

$$ \rho_{f} c_{t} \phi^{m} \frac{{\partial p^{m} }}{\partial t} + \rho_{f} \nabla \cdot ({\mathbf{v}}^{m} ) = Q^{fr - m} + Q_{w} $$
(1)

For the fracture media, we have:

$$ \rho_{f} c_{t} \phi^{fr} \frac{{\partial p^{fr} }}{\partial t} + \rho_{f} \nabla \cdot ({\mathbf{v}}^{fr} ) = Q^{m - fr} + Q^{fri - frj} $$
(2)

In Eqs. (1)–(2), ρf is the fluid density; p and T represent the pressure and temperature field, respectively; t is the time; ϕ is the porosity; ct is the coefficient of compressibility; Qw is the source or sink term, described by the Peaceman formula [19, 20].

Velocity fields in Eqs. (1)–(2) are calculated by Darcy’s equation.

$$ v^{i} = - \frac{{{\mathbf{k}}^{i} }}{{\mu_{f} }}(\nabla p^{i} - \rho_{f} {\mathbf{g}}),\,i = m\,{\text{or}}\,fr $$
(3)

where μf is the dynamic viscosity; k is the permeability tensor; g is the gravity acceleration.

The coupling between the matrix and the fracture is described by the Eq. (4).

$$ Q^{m - fr} = {\text{CI(}}P^{m} - P^{fr} ) $$
(4)

The exchange flow between intersecting fractures is governed by the Eq. (5).

$$ Q^{fri - frj} = {\text{PI(}}P_{i}^{fr} - P_{j}^{fr} ) $$
(5)

In Eqs. (4)–(5), \( {\text{CI = }}A\Xi /\bar{d} \) is the exchange coefficient; A represents the area of matrix grids occupied by the fracture segment; \( \bar{d} \) denotes the average distance from all points in the matrix grid to the fracture segment; \( \Xi \) is the average flow coefficient between the matrix and the fracture. \( {\text{TI}} = \alpha_{i} \cdot \alpha_{j} /(\alpha_{i} + \alpha_{j} ) \) denotes the exchange coefficient between intersecting fractures, in which \( \alpha_{i} = 2b{}_{i}k_{i} /(\mu L_{i} ) \), \( b{}_{i} \) indicates the aperture of the fracture segment. Li represents the length of the ith fracture segment.

2.2 Energy Conservation Equations

According to the idea of the EDFM, two sets of energy equations are still used to describe the matrix and the fracture media, respectively, and the exchange energy between them should be considered. For the matrix media:

$$ (\rho c_{p} )_{eff} \frac{{\partial T^{m} }}{\partial t} + (\rho c_{p} )_{f} {\mathbf{v}}^{m} \cdot \nabla T^{m} = \nabla \cdot (k_{eff} \nabla T^{m} ) + E^{fr - m} + Q_{E} $$
(6)

where \( E^{fr - m} = {\text{H}}({\mathbf{V}}^{fr - m} ){\mathbf{V}}^{fr - m} h^{fr} + {\text{H}}( - {\mathbf{V}}^{fr - m} ){\mathbf{V}}^{fr - m} h^{m} + k_{eff} {\text{CI(}}T^{fr} - T^{m} ) \) represents the energy exchange between the matrix and the fracture; \( {\mathbf{V}}^{fr - m} \) represents the relative velocity between the matrix and the fracture; \( H( \cdot ) \) is the Heaviside function; \( Q_{E} \) is the energy term of well, including the injected and produced energy.

The energy conservation equation of the fracture media is described by the Eq. (7).

$$ (\rho c_{p} )_{eff} \frac{{\partial T^{fr} }}{\partial t} + (\rho c_{p} )_{f} {\mathbf{v}}^{fr} \cdot \nabla T^{fr} = \nabla \cdot (k_{eff} \nabla T^{fr} ) + E^{m - fr} + E^{fri - frj} $$
(7)

where \( E^{fri - frj} = ({\text{H}}({\mathbf{V}}^{fri - frj} ){\mathbf{V}}^{fri - frj} h^{fri} ) + {\text{H}}( - {\mathbf{V}}^{fri - frj} ){\mathbf{V}}^{fri - frj} h^{frj} + k_{eff} {\text{TI(}}T^{fri} - T^{frj} ) \) is the energy exchange between different fractures; \( {\mathbf{V}}^{fri - frj} \) represents the relative velocity between the ith fracture and the jth fracture; \( (\rho c_{p} )_{eff} = \phi^{i} (\rho c_{p} )_{f} + (1 - \phi^{i} )(\rho c_{p} )_{s} \) denotes the effective physical parameters; \( k_{eff} = \phi^{i} k_{f} + (1 - \phi^{i} )k_{s} \) represents the effective thermal conductivity;

2.3 Evolutions of the Fluid Properties

In geothermal exploitation, the working fluid is usually injected into the target position at a low temperature. Due to the high temperature and pressure in the thermal reservoir, the thermal properties of the working fluid will be greatly affected [21]. The equation proposed by Wagner and Kurse [22] is used to describe the changing relationship of density and viscosity with pressure and temperature. For the change of density, we have:

$$ \rho (\psi ,\tau ) = p/(RT\psi \lambda_{psi} ) $$
(8)

where \( \psi = p/p^{*} \); the reference pressure \( p^{*} = 16.53{\text{Mpa}} \); \( \tau = T^{*} /T \); the reference temperature \( T^{*} = 1386K \); the specific gas constant of water \( R = 461.526{\text{J}}/(kg \cdot {\text{K}}) \); \( \lambda_{psi} = \sum\limits_{i = 1}^{34} { - n_{i} l_{i} (7.1 - \psi )^{{l_{i} - 1}} (\tau - 1.222)^{{J_{i} }} } \).

For the viscosity, we have:

$$ \mu = \mu^{*} (\tau^{0.5} \sum\limits_{i = 0}^{3} {\tau^{i} } )^{ - 1} \exp (\delta \sum\limits_{i = 1}^{19} {n_{i} (\delta - 1)^{{l_{i} }} (\tau - 1)^{{J_{i} }} } ) $$
(9)

where \( \mu^{*} = 5.5071 \times 10^{ - 5} {\text{Pa}} \cdot {\text{s}} \) is the reference viscosity; \( \tau = 647.226/T \); \( \delta = \rho /317.763{\text{kg/m}}^{ 3} \).

Figure 3 shows the variations of density and viscosity with pressure and temperature. It can be seen from Fig. 3b, the change of viscosity with temperature is significant. From Eq. (3), it is obvious that the change of viscosity will affect the calculation of Darcy’s velocity. Therefore, it is not accurate to set the viscosity as a fixed value in the actual simulation.

Fig. 3.
figure 3

Variations of the fluid property with pressure and temperature

Full size image

3 Model Solution and Verification

3.1 Model Solution

In this paper, the finite volume method (FVM) is used to discretize the flow and heat transfer equations. The seepage equations are discretized by the two-point flux approximation scheme (TPFA). For the energy equations, the convection terms are discretized by the QUICK scheme, the diffusion terms are discretized by the second order central difference scheme, and the unsteady term is discretized by the first-order backward Euler scheme.

3.2 Verification-Five Discrete Fractures

Consider the 2D square matrix with a side length of 100 m as shown in Fig. 2a. Five discrete fractures are predefined in the computational domain. The left boundary simulates the injection well with a pressure of 106 Pa. The right boundary simulates the production well with a pressure of 105 Pa. The initial pressure is 105 Pa. The impermeable boundary conditions are applied on the top and bottom boundaries. The injection temperature at the left boundary is 20 ℃ and the initial temperature is 180 ℃. The remaining boundaries are adiabatic. Other main parameters are shown in Table 1. The time step is one year and the total simulation time is 40 years.

Table 1. Main model parameters
Full size table

To verify the correctness of the model, the results are compared with the commercial software-COMSOL. Figure 4a shows the grid used in COMSOL. The total number of free triangular grids is 7,019. Figure 4b is the computational grid adopted in the proposed model. The number of matrix grids is 2,500. The number of fracture grids is 121 and the total computational grid is 2,621. The computing time of COMSOL is 17 s and that of the proposed model is 8 s. Figure 5 shows the changes in pressure and temperature fields after 40 years. The results of COMSOL are displayed in the left and the right is the results of the proposed model. There is no apparent difference between these two groups of results. To quantitatively verify the computational accuracy of the pressure and temperature fields, the horizontal fracture in Fig. 4 is selected. With COMSOL as the reference solutions, the computational results are shown in Fig. 6. It can be seen that both the pressure and temperature field are in good agreement, which verifies the correctness of the model and algorithm in this paper.

Fig. 4.
figure 4

Computational grids for the COMSOL and the proposed model.

Full size image
Fig. 5.
figure 5

Results comparison of COMSOL (left) and the proposed model (right).

Full size image
Fig. 6.
figure 6

Comparisons of pressure and temperature in the horizontal fracture with COMSOL.

Full size image

4 Model Application

The EGS can be exploited after artificial hydraulic fracturing. The fracture network will affect the whole thermal extraction process. However, due to the complex structure of deep thermal reservoirs, it is difficult to get all the information of fracture distributions. Some studies have shown that the fracture network after artificial fracturing has fractal characteristics to a certain extent [23]. Therefore, based on fractal theory, fractal trees with 63 fractures are generated to represent the fracture network of deep thermal reservoir. As shown in Fig. 7b, the start-point and end-point coordinates of the first trunk are (20 m, 20 m) and (40 m, 40 m). The angle between the branch and the horizontal and vertical directions is 30°, and the number of iterations of the fractal tree is 4 times. The computational grids of fractal fracture tree are shown in Fig. 8. Each matrix grid contains only one fracture segment. The number of matrix grids is 10,000. The number of fracture grids is 409 and the total number of grids is 10,409. The time step is one year and the total simulation time is 40 years. The specific model parameters are shown in Table 2. The boundary conditions of the model are as follows.

Fig. 7.
figure 7

Fractal fracture tree including 63 fractures

Full size image
Fig. 8.
figure 8

Computational grids of fractal fracture tree

Full size image
Table 2. Parameters of the fractal tree model
Full size table
  1. (1)

    Pressure boundary: all boundaries are subjected to an impermeable boundary condition. The injection well is located at the bottom left corner, and the injected flow rate is 130 m3/day. The produced pressure of the production well at the top right corner is 105 Pa.

  2. (2)

    Temperature boundary: The matrix is surrounded by adiabatic boundary conditions. The initial reservoir temperature is 180 ℃ and the injected temperature is 20 ℃.

4.1 Temporal and Spatial Evolutions of Pressure Field

Figure 9a, b and c show the spatial distributions of pressure field in thermal reservoirs at different years. The existence of fractures makes the seepage field exhibits anisotropy. Because of the high permeability of fractured media, it provides the main pathways of fluid. Therefore, the pressure in fractured media decreases rapidly. Due to an artificial fracture system existing between the injection well and the production well, the hydraulic exchange between the matrix and the fracture is significant in this area, which leads to a rapid rise in pressure of the surrounding matrix.

Fig. 9.
figure 9

Pressure field distributions at various time

Full size image

As the time of exploitation goes on, the pressure near the injection well increases gradually, while the pressure near the production well decreases gradually. This is mainly due to the obvious temperature difference between wells. In the early stage of exploitation, with injecting the cold water, the temperature near the injection well decreases rapidly. The density and viscosity become large, which causes the high flow resistance in the area, whereas the trend near production well is the opposite. However, after 20 years of exploitation, the temperature near the production well began to decrease. Therefore, the density and viscosity of fluid increase, which leads to the increasing pressure gradient around the production well, while that around the injection well is the opposite. If we do not consider the changes in density and viscosity with pressure and temperature in the simulation, the evolutions of the pressure field around the wells cannot be observed. This shows that the influences of fluid properties on the pressure field cannot be ignored.

4.2 Temporal and Spatial Evolutions of Temperature Field

Figure 10a, b and c show the spatial distributions of temperature field in thermal reservoirs for 5, 20 and 40 years. It can be seen that the range of temperature decrease gradually expands to the production well with the growing time of exploitation. After 5 years of exploitation, the cold front has advanced to the root zone of fracture tree. In the early stage of exploitation, the heat exchange between low-temperature fluid and high-temperature rock mass is mainly through the heat conduction. The temperature change near the injection well is roughly an arc-shaped distribution.

Fig. 10.
figure 10

Temperature field distributions at various time

Full size image

The leading edge of temperature after 5 years extends to a distance of about 56 m from the injection well and the trailing edge, which is consistent with the bottom hole temperature, is about 10 m. When the exploitation time is 20 years, the cold front has been pushed into most areas of the fracture trees. At this time, the effect of fracture conductivity become prominent and the low-temperature zone in the reservoir is also enlarged. Meanwhile, the temperature distributions show anisotropy. Because of the high permeability of fractured media, the internal velocity is faster. Therefore, the convective heat transfer is obviously enhanced. When the time of exploitation reaches 40 years, the recovery temperature is lower than the initial temperature. The thermal breakthrough has gradually formed and the heat recovery efficiency of the EGS began to decline.

5 Conclusion

This paper provides a way to simulate the process of geothermal energy exploitation. According to the idea of differentiating matrix and fracture in the EDFM, two sets of governing equations are adopted to describe the fluid flow and heat transfer processes of matrix and fracture, respectively. The proposed model also considers the influence of fluid properties. The spatial-temporal evolutions of pressure and temperature fields are analyzed for a hypothetical EGS case. The main conclusions are as follows.

  1. (1)

    The accuracy of the proposed model and in-house code is proved by comparing with commercial software, COMSOL Multiphysics. For numerical simulation of the flow and heat transfer in EGS with fractures, the finite volume method based on the embedded discrete fracture model can be a good alternative. The matrix can be divided by the structured grids, which is not restricted by the fractures and avoids using the complicated unstructured grids.

  2. (2)

    The distributions of the pressure field are mainly affected by the fracture distributions. The pressure of the matrix will increase rapidly in areas with dense fractures. The interwell pressure gradient is mainly influenced by the thermal reservoir permeability and the fluid properties. The distribution of leading edge of the temperature field is obviously affected by the fracture conductivity, while the trailing edge is mainly affected by the injected water temperature. The interwell temperature distributions are greatly influenced by the through fractures.