A combined method using Lattice Boltzmann Method (LBM) and Finite Volume Method (FVM) to simulate geothermal reservoirs in Enhanced Geothermal System (EGS)

Abstract

With the development of industrial activities, global warming has accelerated due to excessive emission of CO₂. Enhanced Geothermal System (EGS) utilizes deep geothermal heat for power generation. Although porous medium theory is commonly employed to model geothermal reservoirs in EGS, Hot Dry Rock (HDR) presents a challenge as it consists of impermeable granite with zero porosity, potentially distorting the physical interpretation. To address this, the Lattice Boltzmann Method (LBM) is employed to simulate CO₂ flow within geothermal reservoirs and the Finite Volume Method (FVM) to solve the energy conservation equation for temperature distribution. This combined method of LBM and FVM is implemented using MATLAB. The results showed that the Reynolds numbers (Re) of 3,000 and 8,000 lead to higher heat extraction rates from geothermal reservoirs. However, higher Re values may accelerate thermal breakthrough, posing challenges to EGS operation. Meanwhile, non-equilibrium of density in fractures becomes more pronounced during the system's life cycle, with non-Darcy's law becoming significant at Re values of 3,000 and 8,000. Density stratification due to buoyancy effects significantly impacts temperature distribution within geothermal reservoirs, with buoyancy effects at Re=100 under gravitational influence being noteworthy. Larger Re values (3,000 and 8,000) induce stronger forced convection, leading to more uniform density distribution. The addition of proppant negatively affects heat transfer performance in geothermal reservoirs, especially in single fractures. Practical engineering considerations should determine the quantity of proppant through detailed numerical simulations.

Graphical abstract

Keywords

Lattice boltzmann method / Finite volume method / Enhanced geothermal system / Geothermal reservoir / Proppant / Re / Heat extraction rate

Introduction

With the surge in industrial activities, the acceleration of global warming due to excessive CO₂ emission has become a pressing concern (Sun et al. 2022). Among various sources, the combustion of fossil fuels for power generation accounts for a substantial portion of CO₂ emission (Luo et al. 2021). In response to this challenge, renewable energy resources have gained significant attention (Ding et al. 2024). Geothermal energy, in particular, has emerged as a promising option for its abundant reserves, minimal pollution emissions, and consistent power generation capabilities (Zhang et al. 2022; Rohit et al. 2023). Enhanced Geothermal System (EGS) represents a key approach for utilizing deep geothermal heat (Li et al. 2022), making them applicable across diverse underground environments.

An EGS typically consists of injection wells, production wells, a geothermal reservoir, and power generation infrastructure (Avanthi Isaka et al. 2019). The power generation system in EGS often employs the Organic Rankine cycle, based on the range of geothermal gradients. Moreover, the geothermal reservoirs typically feature multi-fractures induced by hydraulic fracturing, forming a complex fracture network comprising natural and artificial fractures (Dokmak et al. 2024). These fractures serve as the primary conduits for geothermal fluid. Notably, the pioneering Fenton Hill project in the United States, located in Valles Caldera, demonstrated the potential of EGS, utilizing impermeable granite formations known as Hot Dry Rock (HDR) (Olasolo et al. 2016). As a result, the impermeability of most HDR formations necessitates hydraulic fracturing to enhance EGS performance.

Considering the costly and time-consuming nature of experimental tests, numerical simulations offer an effective method to assess ESG performance, including heat transfer and seepage flow within geothermal reservoir. Conventionally, geothermal reservoirs are treated as porous medium, allowing for the modelling of coupled thermo-hydraulic-mechanical-chemical mechanism.

Mukuhira et al. (2020) introduced a novel method for examining the entire flow path system within a geothermal reservoir, successfully delineating the main permeable zone and fluid flow tendencies. In most cases, geothermal fluid was observed to flow through the lower center of the production well. Hu et al. (2022) utilized 3DEC software to develop a numerical model of the FORGE site, incorporating relevant geological structures and natural fracture networks. Subsequently, PFC2D software was employed to analyze the variation in natural fractures. The findings indicated that cyclic injection facilitated the expansion of fracture opening in the HDR matrix, while gradient injection was more conducive to activating existing natural fractures. Chen et al. (2019) developed a numerical model for naturally occurring geothermal reservoirs, integrating a coupled thermo-poroelastic model, a stochastically generated fracture network model, and a shear dilation based geomechanical model to capture the coupled Thermal-Hydraulic-Mechanical (THM) processes within the reservoir. The results highlighted injection-induced state changes around the injection wellbore, explaining differences between the pre- and post-simulated stress states. Within the framework of FLAC3D-TOUGH2MP, Liao et al. (2023) established a 3D THM numerical model embedding the DFN model to investigate the long-term performance of the CO₂-EGS in complex fracture networks. The results revealed that the temperature decreases led to variations in thermal stress and increased fracture width, subsequently reducing injection pressure. Sun et al. (2017) developed a numerical model for fractured EGS reservoir, accounting for local thermal non-equilibrium and Thermo-Hydraulic-Mechanical (THM) coupling effects. The model treated HDR matrix blocks and discrete fractures as differential porous media, demonstrating improved coupling between EGS and heat exchanger through water as heat transfer medium. Yao et al. (2018) further explored the heat production performance and economic benefits of EGS using a 3D stochastically generated fracture model, based on a case from the Desert Peak geothermal field. Results indicated sensitivity of average outlet temperature to factors such as distance between injection and production wells, fracture characteristics, injection temperature and HDR characteristics. Zhou et al. (2022) employed a coupled Thermo-Hydraulic (TH) model with a discrete fracture network, accounting for variability in fracture aperture and associated DFN properties. The study revealed that increased mean fracture aperture resulted in higher heat production. Li et al. (2021) utilized an embedded discrete fracture model and extended finite element method to simulate the Thermal-Hydraulic-Mechanical (THM) coupling process in a hypothetical EGS, analyzing the evolution of pressure, temperature and displacement. Zhang et al. (2023) developed a Thermo-Hydraulic-Mechanical-Chemical (THMC) model with an integrated discrete break model, evaluating the environmental benefits of EGS and investigating the effects of operational and physical parameters on production performance, including the injection temperature, the injection flow rate, the initial reservoir temperature and the initial reservoir permeability. The study indicated efficient temperature production and power generation coefficient over the first 5 years, with subsequent decreases by 7.8% and 15.4%, respectively, over 15 years. This work has valuable insights for both coupled THMC modeling and practical engineering of EGS.

In the field of EGS, there has been limited research on employing the Lattice Boltzmann Method (LBM) to solve for temperature, velocity and pressure distributions. Jiao et al. (2021) introduced a pore-scale model that integrates LBM with the Discrete Element Method (DEM). This comprehensive approach considers Thermo-Hydraulic-Mechanical coupling processes and accurately analyzes convective heat transfer, fluid flow and thermal strain, particularly in fractures. Similarly, Al Balushi et al. (2023) explored the concept of fracture conductive tuning using temperature-sensitive proppants to enhance heat extraction and delay thermal breakthrough from the particle scale to the reservoir scale. Their LBM-based model demonstrated a 47.8% improvement in heat extraction compared to typical proppants, offering a robust method for determining the properties of such temperature-sensitive proppants to achieve uniform thermal gradients along different subsurface flow paths autonomously.

Despite these advancements, most EGS models have traditionally relied on porous medium theory to describe the fractured reservoir at macro-dimensions, neglecting detailed fracture distributions. Incorporating such details necessitates substantial computational resources. Additionally, geothermal reservoirs are typically treated as porous media, overlooking the impermeable nature of HDR formations with zero porosity. This oversight may lead to deviations in the physical interpretation of HDR. Based on the work of Jiao et al. (2021) and Al Balushi et al. (2023), LBM demonstrates high computational effectiveness in identifying fractures and modelling geothermal fluid flow, particularly in complex fracture networks. By simplifying meshes, LBM efficiently addresses the complexities of fracture distribution in geothermal reservoir, enabling accurate modelling of heat transfer and flow in fractures. Conversely, the Finite Volume Method (FVM) is commonly employed in EGS research due to its effectiveness in simulating heat transfer and fluid flow in porous media. However, FVM often neglects specific fracture distributions within geothermal reservoirs, requiring additional computational resources to establish governing equations for heat transfer and fluid flow in fractures.

To overcome the limitations of porous medium theory, this paper employs LBM to model CO₂ fluid flow in geothermal reservoirs, while FVM is utilized to solve the energy conservation equation for temperature distribution. Implemented using MATLAB, this combined LBM-FVM approach enables the analysis of velocity, density and temperature fields in fractured reservoirs with complex fracture networks. Furthermore, the Reynolds number (Re) is used to evaluate the characteristic of heat transfer and CO₂ flow in geothermal reservoirs. Subsequently, the effects of proppant in complex and single fractures are compared, and the buoyancy effect of CO₂ in single fractures is investigated.

1 Mathematical model

1.1 Model description

In this study, a two-dimensional, micro-compressible and non-Darcy's law flow in a geothermal reservoir is investigated. The reservoir is conceptualized as a porous medium with artificial fractures, and CO₂ serves as the geothermal fluid for heat transfer within the fractured reservoir. CO₂ enters the reservoir at the inlet area with a uniform velocity and temperature. The top and bottom layers of the reservoir act as impermeable boundaries with constant temperatures, determined by data from the Gonghe site (Zhong et al. 2022). The boundaries at the top and bottom of the reservoir exhibit no-slip conditions. The operational process of the EGS is shown in Fig. 1, involves injecting CO₂ into the reservoir via an injection well, facilitating heat transfer between CO₂ and HDR within fractures. Thus, fracture is the only channel for CO₂ flow in geothermal reservoir. Subsequently, CO₂ is extracted from the reservoir via a production well, providing high temperature CO₂ for power generation. Thus, understanding the characteristic of heat transfer and fluid flow is crucial for optimizing power generation efficiency in EGS.

The objective of this study is to examine the impact of different inlet temperatures and pressures on temperature and pressure distributions within the geothermal reservoir, aiming to identify the optimal configuration for heat transfer efficiency throughout the system's life cycle.

1.2 Governing equation of seepage flow

The Lattice Boltzmann Method (LBM) is employed to simulate the flow and heat transfer in this study using the Bhatnager-Gross-Ktook model (Qian et al. 1992). The BGK model is widely used in LBM to describe particle collisions. The D2Q9 model is used to solve the convection-diffusion problem within the fractured reservoir, where nine velocity directions represent the motion of particles in a two-dimensional space.

The discrete distribution function f_i with status parameter of discrete velocity c_i and the evolution equation are expressed as follows (Babamahmoudi et al. 2023):

(1)$ {{f}}_{{i}}\text{(}\text{x}{+\Delta }{x}\text{,}{t}{+\Delta }{t})-{{f}}_{{i}}\text{(}{x}\text{,}{t}\text{)=}{{ \varOmega }}_{{i}}\text{(}{f}\text{)} $

Where: Ω_i represents the discrete collision operator. The BGK collision operator, associated with a Single Relaxation Time (SRT) model, can be written as:

(2)$ {{ \varOmega }}_{{i}}\text{(}{f}\text{)=}\frac{\text{1}}{\text{τ}}\text{(}{{f}}_{{i}}^{\text{eq}}-{{f}}_{{i}}\text{)} $

Where: τ is the relaxation time and f_i^eq is the equilibrium distribution function.

1.2.1 LB equation of the velocity field

According to Equation (1) and Equation (2), the discrete distribution function evolution equation under SRT-LBM with the BGK approximation can be expressed as:

(3)$ {{f}}_{{i}}\text{(}{x}{+\Delta }{x}\text{,}{t}{+\Delta }{t}\text{)-}{{f}}_{{i}}\text{(}{x}\text{,}{t}\text{)=}-\frac{1}{\tau }\text{[}{{f}}_{{i}}\text{(}{x}\text{,}{t})-{{f}}_{{i}}^{\text{eq}}\text{(}{x}\text{,}{t}\text{)]} $

Where: Δx = c_i×Δt, τ is the single relaxation time, f_i^eq represents the equilibrium distribution function and can be written as:

(4)$ {{f}}_{{i}}^{\text{eq}}\text{(}{x}\text{,}{t}\text{)=}\rho \omega_{{i}}\left[\text{1+}\frac{{\overrightarrow{{c}}}_{{i}}\cdot \overrightarrow{{u}}}{{{c}}_{{s}}^{\text{2}}}\text+\frac{{\text{(}\overrightarrow{{c}}}_{{i}}\cdot \overrightarrow{{u}}{\text{)}}^{\text{2}}}{\text{2}{{c}}_{{s}}^{\text{4}}}-\frac{\text{|}\overrightarrow{{u}}{\text{|}}^{\text{2}}}{\text{2}{{c}}_{{s}}^{\text{2}}}\right] $

(5)$ {\overrightarrow{{c}}}_{{i}}=\left\{\begin{array}{l}\text{(}\text{0,0}\text{)}\quad {i}\text=\text{0}\\ \text{(}\text{c}\text{o}\text{s}\text{[}\text{(}{i}\text{-1)π}/\text{2}\text{]}\text{,}\;\text{s}\text{i}\text{n}\text{[}\text{(}{i}\text{-1)π}/\text{2}\text{]}\text{)}\quad {i}\text=\text{1,2}\text{,}\text{3,4}\\ \text{(}\text{c}\text{o}\text{s}\text{[}\text{(2}{i}\text{-9)π}/\text{4}\text{]}\text{,}\;\text{s}{i}\text{n}\text{[}\text{(2}{i}\text{-9)π}/\text{4}\text{]}\text{)}\sqrt{\text{2}}\quad {i}\text=\text{5,6}\text{,}\text{7,8}\end{array}\right. $

Where: c_s denotes the lattice sound speed and can be calculated as ${{c}}_{{s}}\text{=1/}\sqrt{\text{3}}$. Meanwhile, ω_i is the weight coefficient of the velocity vectors and can be written as:

(6)$ {{ \omega }}_{{i}}=\left\{\begin{array}{l}\text{4}\text{/}\text{9}\quad {i}\text=\text{0}\\\text{1}\text{/}\text{9}\quad {i}\text=\text{1,2}\text{,}\text{3,4}\\\text{1}\text{/}\text{36}\quad {i}\text=\text{5,6}\text{,}\text{7,8}\end{array}\right. $

The macro-scale parameters including density ρ, velocity u and pressure P are calculated as:

(7)$ \rho \text=\sum _{\text{i}}{{f}}_{{i}} $

(8)$ \rho {u}=\sum _{\text{i}}{{f}}_{\text{i}}\cdot{\overrightarrow{{c}}}_{\text{i}} $

(9)$ {P}\text={{c}}_{{s}}^{\text{2}}\rho $

The kinematic viscosity v of CO₂ is related to the relaxation time τ and the lattice sound speed c_s.

(10)$ {v}=\left(\text{τ}-\frac{\text{1}}{\text{2}}\right){{c}}_{{s}}^{\text{2}}{\Delta }{t} $

1.3 Governing equation of heat transfer

In the field of heat transfer, the Lattice Boltzmann Method (LBM) has evolved into Multi-Speed (MS) and Double Distribution Function (DDF) models (Moradi and D'Orazio, 2023). However, both MS and DDF models encounter difficulties in describing heat transfer under centain boundary conditions. To address this issue, in this study, the Finite Volume Method (FVM) is employed to depict heat transfer in the geothermal reservoir. FVM offers advantages on multiple fronts. Firstly, it allows us to conceptualize the geothermal reservoir as the porosity medium, enabling effective resolution of heat transfer over large scales. Additionally, FVM facilitates a clearer understanding of the physical processes involved, as the discrete equations stem from control volumes representing mass and energy exchanges. In the context of EGS, forced convective heat transfer typically prevails within geothermal reservoir. Consequently, the governing equation for heat transfer can be formulated as follows:

(11)$ \begin{split} &\rho_{\text {eff }} \mathcal{C}_{\mathrm{eff}} \frac{\partial T}{\partial t}+\rho_{\mathrm{eff}} \mathcal{C}_{\mathrm{eff}} u \frac{\partial T}{\partial x}+\rho_{\mathrm{eff}} \mathcal{C}_{\mathrm{eff}} v \frac{\partial T}{\partial y}=\\&\frac{\partial}{\partial x}\left(\lambda_{\mathrm{eff}} \frac{\partial T}{\partial x}\right)+\frac{\partial}{\partial y}\left(\lambda_{\mathrm{eff}} \frac{\partial T}{\partial y}\right)\end{split} $

(12)$ \lambda_{\text {eff }}=\phi \lambda+(1-\phi) \lambda_{\mathrm{s}} $

(13)$ \alpha=\frac{\lambda_{\mathrm{eff}}}{\rho c-(1-\phi) \rho_{\mathrm{s}} c_{\mathrm{s}}} $

Where: λ_s is the thermal conductivity coefficient of HDR, c_s the specific heat capacity of HDR, ρ_s the density of HDR, x the positive in the direction from injection well to production well, y the positive in the upward direction, u the seepage velocity in the x direction, v the seepage velocity in the y direction, T the temperature of geothermal reservoir, is the porosity of geothermal reservoir, λ_eff the effective thermal conductivity coefficient of geothermal reservoir, c_eff the effective specific heat capacity of geothermal reservoir, ρ_eff the effective density of geothermal reservoir, and α the effective thermal diffusion coefficient.

2 Numerical solution method

To successfully simulate seepage flow and heat transfer in a fractured reservoir using LBM, it's crucial to establish appropriate boundary conditions. The corresponding boundary conditions for f_i and energy conservation equation are illustrated in this section.

2.1 Initial and boundary conditions

The inlet boundary condition is set as a specific velocity to simulate the outlet condition of the injection well. Based on Fig. 1, when CO₂ flows uniformly into the fractured reservoir, the distribution functions f₁, f₅ and f₈ are unknown. According to Zou and He (1997), f₁, f₅ and f₈ can be solved as follows:

(14)$ f_1=f_3+\frac{2}{3} \rho u_0 $

(15)$ f_5=f_7-\frac{1}{2}\left(f_2-f_4\right)+\frac{1}{6} \rho u_0 $

(16)$ f_8=f_6-\frac{1}{2}\left(f_2-f_4\right)+\frac{1}{6} \rho u_0 $

(17)$ \rho=\frac{1}{1-u_0}\left[f_0+f_2+f_4+2\left(f_3+f_6+f_7\right)\right] $

Where: u₀ is the inlet velocity in geothermal reservoir, f is the distribution function, the subscript of f is corresponding to the direction in Fig. 2.

The outlet velocity of outlet boundary condition is unknown and it is assumed as a developed flow in the fractured reservoir. Guo et al. (2002) method was used for the developed flow boundary conditions in lattice Boltzmann method. According to Fig. 1, the unknown out velocity in outlet boundary conditions can be solved by the second order extrapolation as follows:

(18)$ f_3=2 f_{(3, n x-2)}-f_{(3, n x-1)} $

(19)$ f_6=2 f_{(6, n x-2)}-f_{(6, n x-1)} $

(20)$ f_7=2 f_{(7, n x-2)}-f_{(7, n x-1)} $

For the upper and bottom boundaries of the fractured reservoir, It is assumed that obstacles boundary conditions due to the impermeable nature of the HDR without hydraulic fracturing. Non-slip boundary conditions are applied to the upper and bottom boundaries of the fractured reservoir. The bounce-back scheme (He et al. 1997) is employed to define these boundary conditions as follows:

(21)$ f_i=f_{ {inv }(i)} $

Where: f_i is the unknown density function in direction i and f_inv(i) is the density distribution function in the opposite of direction i for the fluid.

For energy conservation equation of the fractured reservoir, the boundary conditions can be written as:

(22)$ \left.\bar{T}\right|_{x=0}=\left.T_1\right|_{z=b}, 0 \leqslant z \leqslant b, 0 \leqslant x \leqslant l, 0 \leqslant y \leqslant h $

Where: b represents the depth of wellbore, l denotes the length of the geothermal reservoir, h signifies the height of geothermal reservoir. The subscript 1 indicates the parameter associated with the injection well. A constant temperature boundary condition is applied to both the upper and bottom boundaries of the fractured reservoir, in accordance with the geothermal gradient.

2.2 Solution process

In this study, CO₂ is injected into the fractured reservoir under high pressure from the injection well. Consequently, forced convective heat transfer occurs within the fractured reservoir. The velocity and density fields influence only the temperature field. There exists an uncoupled relationship between the governing equations of fluid flow and heat transfer. In addition, the solution necessitates discrete equations for energy conservation. These discrete equations can be expressed as follows:

(23)$ \begin{split} & \frac{T_{\mathrm{i}, \mathrm{j}}^{\mathrm{n}+1}-T_{\mathrm{i}, \mathrm{j}}^{\mathrm{n}}}{\Delta t}+u_{\mathrm{i}, \mathrm{j}}^{\mathrm{n}} \frac{T_{\mathrm{i}, \mathrm{j}}^{\mathrm{n}+1}-T_{\mathrm{i}-1, \mathrm{j}}^{\mathrm{n}+1}}{\Delta x}+v_{\mathrm{i}, \mathrm{j}}^{\mathrm{n}} \frac{T_{\mathrm{i}, \mathrm{j}}^{\mathrm{n}+1}-T_{\mathrm{i}, \mathrm{j}-1}^{\mathrm{n}+1}}{\Delta y}=\\& \alpha_{\mathrm{eff}, \mathrm{i}}^{\mathrm{n}} \frac{T_{\mathrm{i}+1, \mathrm{j}}^{\mathrm{n}+1}-2 T_{\mathrm{i}, \mathrm{j}}^{\mathrm{n}+1}+T_{\mathrm{i}-1, \mathrm{j}}^{\mathrm{n}+1}}{(\Delta x)^2}+\alpha_{\mathrm{eff}, \mathrm{i}}^{\mathrm{n}} \frac{T_{\mathrm{i}, \mathrm{j}+1}^{\mathrm{n}+1}-2 T_{\mathrm{i}, \mathrm{j}}^{\mathrm{n}+1}+T_{\mathrm{i}, \mathrm{j}-1}^{\mathrm{n}+1}}{(\Delta y)^2}\end{split} $

Where: Index i and j represent x and y direction, n represents the current time step. The convective term is discretized using first-order windward format, and the diffusion term is discretized using the first-order center differential format.

The solution process for the combined method with LBM and FVM is as follows:

(1) Define the variables T, u and ρ and input their initial value;

(2) Input the weight coefficients of the velocity vectors and the initial equilibrium distribution function f_i^eq;

(3) Set the residual criteria for LBM and FDM;

(4) Calculate the inlet and outlet discrete LB equation and calculate the process of collision and propagation.

(5) Output the velocity and density distribution in the geothermal reservoir;

(6) Update the thermophysical properties of CO₂ and calculate the discrete equation of the energy conservation equation.

(7) Output the temperature distribution in the geothermal reservoirs using Gauss-Seidel iteration.

(8) Repeat steps (4) to (7) until all time steps have been completed.

(9) Output the velocity, density and temperature distribution in the geothermal reservoir.

2.3 Validation

2.3.1 Grid independency

In this paper, grid independence validation was conducted using grid numbers of 50×100, 100×200, 150×300, 200×400 and 250×500. To compare the results of the velocity distribution, $\bar u/u_{\rm{max}} $ is used as the criterion to evaluate the grid independency. Fig. 3 illustrates this criterion for different grid number. According to the comparison, it was determined that the grid with dimensions of 150×300 showed acceptable accuracy, with a maximum error of 0.5%.

2.3.2 Comparison between numerical results and analytical solution

Currently, there is only one real EGS project in commercial operation, making it difficult to obtain experimental results for validation. However, validation of the numerical model using the LBM method is necessary. To verify the numerical method, a lid-driven cavity with fluid flow and heat transfer is employed for comparison. In this case, only the upper boundary condition is considered as a slip boundary condition. The numerical results using FDM and LBM are illustrated in Fig. 4, with the velocity u used as a criterion to evaluate the difference between FDM and LBM. The dimension of the lid-driven cavity is 200 m × 200 m and the Re in the cavity is 150 during fluid flow and heat transfer. As shown in Fig. 3, there is good agreement between the two methods, with a maximum error of only 0.57%.

Next, an analytical solution of heat transfer and fluid flow in a single fracture is employed to verify the numerical model. The heat transfer and fluid flow in a single fracture are shown in Fig. 5. The analytical results are presented as the following equation (Ghassemi et al. 2008):

(24)$ T_{{D}}(x, y, t)=\operatorname{erfc}\left[\frac{x}{u_{\mathrm{i}} c_{\mathrm{f}} \rho_{\mathrm{f}}} \sqrt{\frac{\lambda_{\mathrm{r}} \rho_{\mathrm{r}} c_{\mathrm{r}}}{t}}+\frac{y}{2} \sqrt{\frac{\rho_{\mathrm{r}} c_{\mathrm{r}}}{\lambda_{\mathrm{r}} t}}\right] $

Where: T_D represents the dimensionless temperature, which can be written as T_D = (T_r0-T)/(T_r0-T_f0). T_r0 denotes the initial temperature in the rock region at a specific point. T_f0 represents the initial temperature of the geothermal fluid in the single fracture. u_i signifies the inlet temperature of the geothermal fluid in the single fracture, and t denotes time. The suffix r indicates the rock region, and the suffix f signifies the geothermal fluid in the single fracture. Additionally, ρ represents density, and c denotes the constant pressure specific heat capacity.

The initial and boundary conditions for the validation case are presented in Table 1. Fig. 6 depicts both the analytical and numerical solutions in the single fracture. Both results indicate that the production temperature stabilizes when the time reaches approximately 20 years. This stabilization occurs because the production temperature decreases due to the injection of cooling fluid. When the heat extraction rate of the injection fluid matches the heat addition rate from the rock temperature, the production temperature tends to stabilize. The results demonstrate that the maximum error of production temperature is 0.98%. Therefore, the numerical model exhibits high accuracy in solving the heat transfer and fluid flow in the geothermal reservoir.

3 Results and discussion

To illustrate the impact of CO₂ flow on performance of EGS, a 30-year operation time is employed to depict the variations in velocity, density and temperature fields. To demonstrate the characteristic of heat transfer and fluid flow in the geothermal reservoir, we focus on a small-scale region within the reservoir, as depicted in Fig. 7. Fig. 7 show cases the scale of complex fracture distributions in the geothermal reservoir, where Fig.7(a) shows complex fracture distributions, Fig. 7(b) displays complex fracture distributions with proppant, Fig. 7(c) represents a single fracture, and Fig. 7(d) illustrates a single fracture with proppant. The proppants are simplified as spherical with diameters of 10 mm in complex fractures distributions and 5 mm in single fractures. Additionally, obtaining accurate distributions of proppants in fractures is challenging. Hence, it is simplified that proppants follow a Monte Carlo distribution within the fractures. The initial and boundary conditions for this section are outlined in Table 2.

3.1 Effect of Re on the performance of geothermal reservoirs

In this section, Re values of 100, 3,000, and 8,000 are used to compare the effect of CO₂ flow on heat transfer in the geothermal reservoirs. Fig. 8 displays the fractures distribution in the geothermal reservoirs while Fig. 9 illustrates the temperature distribution in the reservoirs over a period of 30 years. Over the 30-year period, the low-temperature area expands with increasing Re. This suggests that heat transfer is enhanced with higher Re, leading to a higher rate of heat extraction during the 30-year operation. In the initial stages of operation, the low-temperature area is larger with Re of 100 compared to Re values of 3,000 and 8,000. This indicates that the slower velocity of CO₂ at lower Re facilitates a more fully developed convective heat transfer compared with that of the larger Re of 3,000 and 8,000. The slower velocity allows CO₂ to more effectively transmit through fractures, resulting in a higher rate of heat extraction at the early stages of operation. However, as the operation progresses, a higher rate of heat extraction is achieved with increasing Re. Fractures with Re values of 3,000 and 8,000 lead to lager low-temperature areas along the fractures compared to the small Re of 100. If this low-temperature area extends to the production well, it may lead to thermal breakthrough in the geothermal reservoirs, indicating diminished potential for sustained power generation. Therefore, lower Re values contribute to improved heat extraction rates in the early stages of operation. However, as the injection process continues, higher Re values lead to higher rates of heat extraction in the geothermal reservoirs. Larger Re values reduce the time required to reach thermal breakthrough, which is not favorable for the operation of Enhanced Geothermal System (EGS).

Fig. 10 illustrates the density distribution over a period of 30 years in the geothermal reservoirs. In the initial stages of operation, the density of CO₂ is the highest gradually decreasing over the 30-year period of EGS operation? Moreover, higher Re values of 3,000 and 8,000 result in higher densities of CO₂ within fractures. Consequently, higher initial pressures are observed in the early stages of operation with these higher Re values. However, these higher initial pressures negatively impact heat transfer performance, necessitating higher injection pressures and consuming more power generation. As the heat transfers progresses, fractures with lower Re values (e.g. 100) exhibit a disequilibrium density distribution within the complex fracture network. This is due to weaker forced convective heat transfer within fractures with lower mass flow rate of CO₂, making them more susceptible to the high temperatures of HDR. This leads to a larger density difference across different areas of fractures over the 30-year period, with significant non-equilibrium density observed towards the end of EGS operation. Additionally, the density of CO₂ within fractures exhibits a higher rate of decrease over the period of 30 years for higher Re values. For instance, at the beginning of the 30-year period, the maximum density of CO₂ with Re of 100 is 962.0 kg/m³, decreasing to 834.8 kg/m³ by the end of the period. In contrast, for CO₂ with Re of 8,000, the maximum density starts at 1,034.0 kg/m³ and decreases to 885.5 kg/m³ over the same duration. The density drops by 4.24 kg/(m³·a) for CO₂ with Re of 100, and by 4.95 kg/(m³·a) for CO₂ with Re of 8,000. Higher Re values of CO₂ in geothermal reservoirs result in increased seepage velocity, enhancing convective heat transfer intensity within fractured reservoirs. Then, the temperature drops rates in geothermal reservoirs with Re of 3,000 and 8,000 are higher compared to those with Re of 100. As a result, according to the Equation of State (EOS) for CO₂, the minimum density of CO₂ increases for all Re values. However, the trend is oppsite for the maximum density of CO₂. CO₂ with higher Re values such as 3,000 and 8,000 experiences higher temperature or lower pressures. This leads to larger pore pressures within fractures for CO₂ with higher Re values compared to those with lower Re values. Consequently, increasing injection pressure with higher Re values only improves the maximum density of CO₂ within fractures without deceasing the maximum density of CO₂. The higher temperature of CO₂ within fractures is the primary reason for the larger drop rate in density observed over the 30-year period for CO₂ with higher Re values.

On the other hand, when Re is 100 for CO₂ flow in geothermal reservoirs, low-density regions concentrate on the upper surface of the fracture, while high-density regions concentrate on the lower surface. It is evident that weakened convective heat transfer occurs with the small Re of 100, leading to a non-equilibrium density distribution within fractures. In CO₂ flow with non-Darcy's law in fractures, the gravitational effect is more pronounced with the small Re of 100. Additionally, weakened forced convection with a small Re of 100 enhances density stratification. Both gravitational effects and weakened forced convection enhance the buoyancy effect of CO₂ in fractures, resulting in significant stratification of density distribution. In addition, significant density stratification occurs at the end of the EGS operation, indicating the generation of laminar flow. Therefore, when temperature distribution stabilizes, the flow pattern tends towards laminar flow for CO₂ flow with Re of 100. In contrast, this phenomenon is not observed for CO₂ flow with higher Re values of 3,000 and 8,000. Turbulent flow dominates over the 30-year period, resulting in a more uniform density distribution. The influence of non-Darcy's law is more significant for CO₂ flow with Re values larger than 100 in fractures.

3.2 Effect of proppant on the performance of geothermal reservoirs

The addition of proppant is a common method in the oil or gas industry to prevent fractures from closing and maintain production capacity. Previous research has confirmed its importance in this regard. However, there have been few studies on the effect of proppant on heat transfer in fractured reservoirs. In this section, the same fracture distribution from section 4.1 is employed and Re is set at 500. Fig. 11 illustrates the fracture distribution with proppant in geothermal reservoirs, where the proppant used is low-density ceramic granule.

Fig. 12 depicts the temperature distribution with or without proppant over the period of 30 years. The results reveal that the low-temperature regions are larger without proppant compared to those with proppant. Additionally, the distance between the low-temperature region along the fractures and the injection well is greater without proppant than with proppant. This suggests that the temperature drop rate is higher with proppant than without it. Consequently, the addition of proppant reduces the rate of heat extraction of CO₂ in fractures. On one hand, the addition of proppant increases the flow resistance in fractures, leading to a reduction in the velocity of CO₂ within them. Consequently, convective heat transfer is weakened due to the addition of proppant. On the other hand, the addition of proppant does not improve the heat transfer area between CO₂ and HDR; it only alters the flow field within the fractures. Therefore, the addition of proppant has a negative effect on heat transfer for CO₂ in geothermal reservoirs. In essence, the addition of proppant reduces the connectivity of geothermal reservoirs, which is the primary reason for its negative impact on heat transfer. Therefore, when fractures in a geothermal reservoir tend to close, adding proppant is not the optimal method to maintain the fracture structure. Instead, considerations should focus on maintaining the rate of heat extraction, such as improving pore pressure in fractures, to sustain fracture structure and ensure effective heat transfer. An indiscriminate addition of proppant is counterproductive to heat transfer in geothermal reservoirs.

Fig. 13 presents the density distribution over a 30-year period in geothermal reservoirs. With the presence of proppant, the average density gradually increases during the operation of EGS. Meanwhile, there is a noticeable reduction in fracture space compared to when no proppant is present. In the absence of proppant, the density distribution exhibits similar variations to that with proppant, albeit with a lower average density. Under identical boundary conditions with a Re value of 500, the addition of proppant reduces the available fracture space in the fractured reservoir. Consequently, this leads to higher flow resistance within the fractures when proppant is introduced. Subsequently, pore pressure increases with the addition of proppant in fractures. Given the dominance of forced convection in fractures under these conditions, alteration in the flow field due to proppant leads to significant variations in temperature distribution alongside changes in velocity field. Ultimately, the addition of proppant results in a reduction in the rate of heat extraction within fractures, contributing to larger low temperature regions compared to those without proppant. Overall, the inclusion of proppant in fractures negatively impacts CO₂ flow and heat transfer in geothermal reservoirs.

Furthermore, considering the temperature distribution depicted in Fig. 12 and the density distribution shown in Fig. 13, it becomes evident that the low-temperature regions extend along the direction of fractures. Meanwhile, areas with higher density tend to accumulate at the bottom of the geothermal reservoirs. In these regions, the rate of heat extraction significantly surpasses that of other heat extract rates. Consequently, the density stratification induced by buoyancy effects plays a crucial role in shaping the temperature distribution within the fractured reservoir. The influence of buoyancy force, under the influence of gravity, should not be underestimated, unless forced convection exhibits sufficient intensity to uniformly distribute density. Therefore, the quantity of proppant added should be determined after rigorous numerical simulations.

3.3 Effect of buoyancy force on the performance of single fracture

To better elucidate the influence of buoyancy force on heat transfer and fluid flow, a single fracture is introduced into the geothermal reservoir to evaluate the performance of EGS. The single fracture is depicted in Fig. 14, while Fig. 15 illustrates the single fracture with proppant in the geothermal reservoirs. Fig. 16 displays the temperature distribution in the fractured reservoir with a single fracture over a span of 30 years. Both fractures (with and without proppant) exhibit enhanced heat transfer over the 30-year period. However, the fracture without proppant facilitates a higher rate of heat extraction in the fractured reservoir compared to the fracture with proppant. Thus, whether it's a complex fracture network or a single fracture, the addition of proppant tends to have a negative impact on heat transfer in the fractured reservoir. Especially, the effect on heat transfer is more pronounced in the case of the single fracture. In complex fracture distributions, the coefficient of convective heat transfer is enhanced, thereby mitigating the weakening effect of proppant to some extent.

From an engineering perspective, if hydraulic fracturing results in a less connected fractured reservoir, the addition of proppant should not be considered as the method to increase production capacity since it may impair heat extraction performance. Conversely, if a fractured reservoir exhibits sufficient connectivity and a high coefficient of convective heat transfer, the addition of proppant could be considered to maintain fracture space and achieve high connectivity. However, given the complex structure of real fracture networks governing CO₂ flow and heat transfer, the quantity of proppant should be determined through detailed numerical simulations.

4 Conclusions

This paper presents a combined method to address CO₂ flow and heat transfer in fracture reservoirs for EGS. The approach integrates the Lattice Boltzmann Method (LBM) for solving the continuity equation and momentum equations of CO₂ flow, and the Finite Volume Method (FVM) for handling energy conservation equations. The study investigates the effects of Re in geothermal reservoirs and the impact of proppant on the performance of geothermal reservoirs. The key conclusions are summarized as follows:

(1) Higher Re values of 3,000 and 8,000 result in increased heat extraction rates in geothermal reservoirs. However, elevated Re levels may accelerate the time to reach thermal breakthrough, which is detrimental to EGS operation. In the initial stages, a smaller Re of 100 exhibits a higher heat extraction rate in fractures.

(2) In geothermal reservoirs with a small Re of 100, non-equilibrium of density in fractures becomes more pronounced towards the end of EGS operation. Conversely, as Re increases to 3,000 and 8,000, density distribution gradually becomes more uniform, and buoyancy effects diminish. Non-Darcy flow becomes more significant for CO₂ flow at higher Re values such as 3,000 and 8,000.

(3) Density stratification induced by buoyancy effects significantly influences temperature distribution in fractured reservoirs. Buoyancy force impact should be considered, especially if forced convection is insufficient to ensure uniform density distribution.

(4) The addition of proppant negatively impacts geothermal reservoirs, particularly in single fractures where heat transfer is significantly weakened. From an engineering perspective, the quantity of proppant should be determined through detailed numerical simulations to optimize reservoir performance.

References

[1] Al Balushi F, Zhang Q, Dahi Taleghani A. 2023. Improving enhanced geothermal systems performance using adaptive fracture conductivity. Applied Thermal Engineering, 233: 121206. DOI:10.1016/j.applthermaleng.2023.121206.

[2] Avanthi Isaka BL, Ranjith PG, Rathnaweera TD. 2019. The use of super-critical carbon dioxide as the working fluid in enhanced geothermal systems (EGSs): A review study. Sustainable Energy Technologies and Assessments, 36: 100547. DOI:10.1016/j.seta.2019.100547.

[3] Babamahmoudi S, Saeedi Dehaghani AH, Hosseini Moghadam A. 2023. Absolute permeability assessment of porous structures under different boundary conditions using lattice Boltzmann method. Geoenergy Science and Engineering, 221: 211357. DOI:10.1016/j.geoen.2022.211357.

[4] Cheng Q, Wang X, Ghassemi A. 2019. Numerical simulation of reservoir stimulation with reference to the Newberry EGS. Geothermics, 77: 327−343. DOI:10.1016/j.geothermics.2018.09.011.

[5] Ding Q, Huang J, Chen J, et al. 2024. Climate warming, renewable energy consumption and rare earth market: Evidence from the United States. Energy, 290: 130276. DOI:10.1016/j.energy.2024.130276.

[6] Dokmak H, Faraj K, Faraj J, et al. 2024. Geothermal systems classification, coupling, and hybridization: A recent comprehensive review. Energy and Built Environment, 1-20.

[7] Ghassemi A, Nygren A, Cheng A. 2008. Effects of heat extraction on fracture aperture: A poro–thermoelastic analysis. Geothermics, 37: 525−539. DOI:10.1016/j.geothermics.2008.06.001.

[8] Guo Z, Zheng C, Shi B. 2002. An extrapolation method for boundary conditions in lattice Boltzmann method. Physics of Fluids, 14: 2007−2010. DOI:10.1063/1.1471914.

[9] He X, Zou Q, Luo LS, et al. 1997. Analytic solutions of simple flows and analysis of nonslip boundary conditions for the Lattice Boltzmann BGK model. Journal of Statistical Physics, 87: 115−136. DOI:10.1007/BF02181482.

[10] Hu Z, Xu T, Moore J, et al. 2022. Investigation of the effect of different injection schemes on fracture network patterns in hot dry rocks - A numerical case study of the FORGE EGS site in Utah. Journal of Natural Gas Science and Engineering, 97: 104346. DOI:10.1016/j.jngse.2021.104346.

[11] Jiao K, Han D, Li J, et al. 2021. A novel LBM-DEM based pore-scale thermal-hydro-mechanical model for the fracture propagation process. Computer and Geotechnics, 139: 104418. DOI:10.1016/j.compgeo.2021.104418.

[12] Krüger T, Kusumaatmaja H, Kuzmin A, et al. 2016. The Lattice Boltzmann Method-Principles and practice. Springer International Publish.

[13] Li S, Wang S, Tang H. 2022. Stimulation mechanism and design of enhanced geothermal systems: A comprehensive review. Renewable and Sustainable Energy Reviews, 155: 111914. DOI:10.1016/j.rser.2021.111914.

[14] Li T, Han D, Yang F, et al. 2021. Modeling study of the thermal-hydraulic-mechanical coupling process for EGS based on the framework of EDFM and XFEM. Geothermics, 89: 101953. DOI:10.1016/j.geothermics.2020.101953.

[15] Liao J, Hu K, Mehmood F, et al. 2023. Embedded discrete fracture network method for numerical estimation of long-term performance of CO₂-EGS under THM coupled framework. Energy, 285: 128734. DOI:10.1016/j.energy.2023.128734.

[16] Luo Z, Wu Y, Zhou L, et al. 2021. Trade-off between vegetation CO₂ sequestration and fossil fuel-related CO₂ emissions: A case study of the Guangdong–Hong Kong–Macao Greater Bay area of China. Sustain Cities and Society, 74: 103195. DOI:10.1016/j.scs.2021.103195.

[17] Moradi I, D'Orazio A. 2023. Lattice Boltzmann Method Pore-scale simulation of fluid flow and heat transfer in porous media: Effect of size and arrangement of obstacles into a channel. Engineering Analysis with Boundary Elements, 152: 83−103. DOI:10.1016/j.enganabound.2023.04.007.

[18] Mukuhira Y, Ito T, Asanuma H, et al. 2020. Evaluation of flow paths during stimulation in an EGS reservoir using microseismic information. Geothermics, 87: 101843. DOI:10.1016/j.geothermics.2020.101843.

[19] Olasolo P, Juárez MC, Morale MP, et al. 2016. Enhanced geothermal systems (EGS): A review. Renewable and Sustainable Energy Reviews, 56: 133−144. DOI:10.1016/j.rser.2015.11.031.

[20] Rohit RV, Kiplangat DC, Veena R, et al. 2023. Tracing the evolution and charting the future of geothermal energy research and development. Renewable and Sustainable Energy Reviews, 184: 113531. DOI:10.1016/j.rser.2023.113531.

[21] Sun R, Shen J, Grasby SE, et al. 2022. CO₂ buildup drove global warming, the Marinoan deglaciation, and the genesis of the Ediacaran cap carbonates. Precambrian Research, 383: 106891. DOI:10.1016/j.precamres.2022.106891.

[22] Sun ZX, Zhang X, Xu Y, et al. 2017. Numerical simulation of the heat extraction in EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model. Energy, 120: 20−33. DOI:10.1016/j.energy.2016.10.046.

[23] Qian YH, D'Humières D, Lallemand P. 1992. Lattice BGK Models for Navier-Stokes Equation. Europhys Lett, 17: 479. DOI:10.1209/0295-5075/17/6/001.

[24] Yao J, Zhang X, Sun Z, et al. 2018. Numerical simulation of the heat extraction in 3D-EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model. Geothermics, 74: 19−34. DOI:10.1016/j.geothermics.2017.12.005.

[25] Zhang EY, Wen DG, Wang GL, et al. 2022. The first power generation test of hot dry rock resources exploration and production demonstration project in the Gonghe Basin, Qinghai Province, China. China Geology, 5(3): 372−382. DOI:10.31035/cg2022038.

[26] Zhang W, Han D, Wang B, et al. 2023. Thermal-hydraulic-mechanical-chemical modeling and simulation of an enhanced geothermal system based on the framework of extended finite element methods - Embedded discrete fracture model. Journal of Cleaner Production, 415: 137630. DOI:10.1016/j.jclepro.2023.137630.

[27] Zhong C, Xu T, Yuan Y, et al. 2022. The feasibility of clean power generation from a novel dual-vertical-well enhanced geothermal system (EGS): A case study in the Gonghe Basin, China. Journal of Cleaner Production, 344: 131109. DOI:10.1016/j.jclepro.2022.131109.

[28] Zhou D, Tatomir A, Niemi A, et al. 2022. Study on the influence of randomly distributed fracture aperture in a fracture network on heat production from an enhanced geothermal system (EGS). Energy, 250: 123781. DOI:10.1016/j.energy.2022.123781.

[29] Zou Q, He X. 1997. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Physics of Fluids, 9: 1591−1598. DOI:10.1063/1.869307.

RIGHTS & PERMISSIONS

Journal of Groundwater Science and Engineering Editorial Office