# The Numerical Simulation of Heat and Mass Transfer on Geothermal System-A Case Study in Laoling Area, Shandong, China

Authors: Cheng Gao; Le Zhang; Chuanxiang Sun; Jiayuan He
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Geothermal resources have become increasingly attractive and promising because of their abundant resource base and environmental protection, especially for hydrothermal resources, which are widely developed for heating in winter. It is indicated that numerical simulation is an important tool for high efficiency geothermal system development. Compared with other methods, numerical simulation is a more comprehensive, scientific, and effective method in the evaluation of recoverable resources and the formulation of development plans. However, there are some problems in the existing software, such as deviation in physical property calculation, incomplete multifields coupling, and poor applicability for low-permeability reservoir and large-scale models. Therefore, based on the optimized multiphysics coupling mathematical model and MPI architecture, a simulator for a multiphysics-coupling geothermal system (SMG) was developed by Sinopec. A case study in the Laoling area, Shandong, China was primarily conducted by using geothermal software SMG, where the effect of well spacing and injection fluid temperature on the production was illustrated. Moreover, a geothermal recoverable resources evaluation method based on numerical simulation is proposed and used for the evaluation of regional geothermal dynamic recoverable resources in Laoling.


---

## Body

## 1. Introduction

In recent years, the high efficiency development of clean renewable energy has become a significant method to meet the demand of increasing energy and resolve critical environmental problems such as the global warming issue. Geothermal energy is an important existing form of renewable energy. It is widely stored in the rock of the Earth and subsurface natural fluids. Geothermal energy has its own advantages compared to other types of renewable energy, such as independence from weather conditions and seasons [1]. We should pay more attention on geothermal energy to relieve the pressure of energy demand and decrease the amount of CO2 emissions.


### 1.1. Current Situation of Geothermal Exploration and Development

The total geothermal reserves are 1.45 × 1026 J on Earth, which amounts to 4.948 × 1015 t standard coal. Geothermal resources provide energy in the form of electricity and direct heating and cooling, totaling an estimated 147 TW in 2012 [2]. As shown in Figure 1, high-temperature geothermal energy (150°C above) can be used to produce electricity economically; medium temperature geothermal energy (90∼150°C and low-temperature geothermal energy (25∼90°C) can be used mainly on heating, drying, industry, medicine, and other things related to human’s daily life. Based on the data collected by the World Geothermal Congress 2012 [3], it indicates that 24 countries produce electricity by using geothermal energy, the total capacity is 10751 MW and the average coefficient of usage is 72%. The top 10 countries which have the largest amount of electricity produced by geothermal energy are the United States, Philippines, Indonesia, Turkey, New Zealand, Mexico, Italy, Iceland, Kenya, and Japan. Generally, a lot of countries all over the world have paid much attention on the geothermal development.Figure 1
Cascade utilization of geothermal energy.China has primarily formed a pattern of direct usage of geothermal energy, such as heating, bathing, agriculture, and so on, as well as generating electricity by using high-temperature geothermal systems. The total geothermal reserves in China are 2.5 × 1022 J (equivalent to 8.532 × 1011 t standard coal), which are primarily distributed in 12 basins [4]. In addition, the total amount of recoverable resources is 7.5 × 1021 J which is equivalent to 2.56 × 1011 t standard coal. There are no technical limitations on using shallow geothermal resources in medicine, tourism, entertainment, and so on. For electricity-generating projects using high-temperature geothermal resources in China, there are 5 geothermal power plants with a total capacity of about 27.28 MW [5].The development of hydrothermal systems for heating houses in winter can effectively resolve the haze problem. However, the utilization of geothermal resources only accounts for 5% of the renewable resources in eastern China [6]. There are still some critical problems, such as low resource recovery and difficult reinjection [7]. During the recharge process, water temperature, seepage velocity, and aquifer environmental characteristics affect the formation mechanism and development of clogging, which exists in physical, chemical, and biological forms, making it difficult to diagnose the causes of clogging [8]. Therefore, numerical simulations of geothermal reservoirs are a good tool for recoverable resources assessment and production scheme optimization to increase the operating lifetime of the system.For the geothermal resources evaluation, the usual method is a steady-state calculation method such as the surface heat flow method, reservoir volume method, and so on. It is very important to evaluate geothermal resources by considering the effects of geothermal development and the utilization process. For the geothermal developing process, the current method to assign the well spacing is over simplified as an empirical method is used [9]. In addition, the key operating parameters such as production rate, injection rate, and so on, all depend on the energy demand. However, the temperature distribution, fluid flow, stress field, and chemical process have a significant dynamic effect on the operating scheme at different stages. Numerical simulation technology is highly recommended to deal with these important problems to make the best use of geothermal resources. The numerical simulation method is a dynamic geothermal resource evaluation method that can be tracked and optimized [10, 11]. As a geothermal management tool, it can dynamically calculate the amount of geothermal recoverable resources in the process of development. It also can consider the heat supplement inside the geothermal system and the effect of reinjection on the geothermal recovery process.Numerical simulation technology of geothermal system can accurately describe the thermal recovery process and can accurately evaluate the geological resources, recoverable resources, economically recoverable resources, and green economically recoverable resources in the target block of geothermal resources, which provides guidance for project evaluation. Moreover, the production potential of geothermal wells can be predicted under different injection and production schemes, which gives the satisfied well layout and production scheme [12, 13]. To sum up, the development of integrated numerical simulation technology of geothermal system has very important strategic and practical significance for guiding the current development of hydrothermal geothermal resources.


### 1.2. The Status Analysis of the Numerical Simulation Method on Geothermal System

A lot of universities and institutes have developed simulation technologies on geothermal development, and the main commercial softwares are TOUGH2, OpenGeoSys, STARS, and FEHM [14–21]. Table 1 shows major numerical simulation projects conducted by some well-known companies in recent years.Table 1
Recent research on geothermal system simulation.
Institute and timeName of paper/ReportDescriptionNorth California power agency, 2010Integrated geothermal system simulator for the enhancement of geothermal steam production at the geysers [22]Tough2, stars and FEHM are used in reservoir; Thermo Flex is used in well-flowReykjavík university, landsvirkjun co. ltd., 2016Numerical modeling of the hágöngur geothermal reservoir in central Iceland [23]Coupled simulation by using TOUGH2 and PetraSimSupreme energy, 2016A reservoir simulation of the Maura Laboh field, Indonesia [24]Reservoir and well simulation by using TOUGH2 and PetraSimCarotenuto et al. conducted a literature review study on hydrothermal system simulation technology [25]. The results showed that TOUGH2 coupled with well simulation, could effectively deal with geothermal system simulation. Many coupled models were built and associated modeling was done by the software TOUGH2 to analyze hydrothermal or enhanced geothermal systems [26–29]. Fowler et al. compared CSMP-GEMS, OpenGeoSys-GEM, and TOUGHREACT by simulating a geothermal reservoir [30]; they concluded that all these three kinds of software can correctly simulate the process of fluid flow, heat transfer, and chemical reaction in a geothermal reservoir. However, TOUGHREACT is limited to solve the complex geological model because it cannot deal with unstructured grid models. On the other hand, OpenGeoSys-GEM could only simulate quite simple incompressible flow processes because of its dependence on the groundwater seepage model. CSMP-GEMS could handle incompressible flow problems with unstructured grids model. Butler and Enedy [31] conducted an integrated simulation by using TETRAD/TOUGH2/STAR coupled with ThermoFlex/AspenPlus/PEPSE in support of the Geysers geothermal project. Wang et al. analyzed the advantages and disadvantages of numerical simulation technology on Enhanced Geothermal Systems (EGS) and Hydro Geothermal System [32] and described the development significance of EGS. Zhang et al. used TOUGH2 software to study a geothermal system in the Subei area, China [33]. Ganguly et al. reviewed the numerical simulation method on an enhanced geothermal system [34]. The results showed that water/rock interaction or two-phase flow may be an important issue, which should be taken into consideration. The thermal behavior of rock during geothermal energy development has also been investigated by several researchers [35, 36].Through the review on geothermal simulation of the above domestic and foreign software, the calculation of physical parameters of working fluid and heat reservoir is complex and inaccurate, and the operability for engineering application is poor for reinjection optimization and simulation speed. Therefore, based on the thermal-hydraulic-mechanical-chemical coupling mechanism and the MPI parallel computing platform, this study developed and formed a numerical simulation software of geothermal resources, which realizes five functions: thermal dynamic prediction of thermal recovery process, optimization of injection and production system, recoverable resources assessment, optimization of development scheme design, and formulation of development technology policy.


## 1.1. Current Situation of Geothermal Exploration and Development

The total geothermal reserves are 1.45 × 1026 J on Earth, which amounts to 4.948 × 1015 t standard coal. Geothermal resources provide energy in the form of electricity and direct heating and cooling, totaling an estimated 147 TW in 2012 [2]. As shown in Figure 1, high-temperature geothermal energy (150°C above) can be used to produce electricity economically; medium temperature geothermal energy (90∼150°C and low-temperature geothermal energy (25∼90°C) can be used mainly on heating, drying, industry, medicine, and other things related to human’s daily life. Based on the data collected by the World Geothermal Congress 2012 [3], it indicates that 24 countries produce electricity by using geothermal energy, the total capacity is 10751 MW and the average coefficient of usage is 72%. The top 10 countries which have the largest amount of electricity produced by geothermal energy are the United States, Philippines, Indonesia, Turkey, New Zealand, Mexico, Italy, Iceland, Kenya, and Japan. Generally, a lot of countries all over the world have paid much attention on the geothermal development.Figure 1
Cascade utilization of geothermal energy.China has primarily formed a pattern of direct usage of geothermal energy, such as heating, bathing, agriculture, and so on, as well as generating electricity by using high-temperature geothermal systems. The total geothermal reserves in China are 2.5 × 1022 J (equivalent to 8.532 × 1011 t standard coal), which are primarily distributed in 12 basins [4]. In addition, the total amount of recoverable resources is 7.5 × 1021 J which is equivalent to 2.56 × 1011 t standard coal. There are no technical limitations on using shallow geothermal resources in medicine, tourism, entertainment, and so on. For electricity-generating projects using high-temperature geothermal resources in China, there are 5 geothermal power plants with a total capacity of about 27.28 MW [5].The development of hydrothermal systems for heating houses in winter can effectively resolve the haze problem. However, the utilization of geothermal resources only accounts for 5% of the renewable resources in eastern China [6]. There are still some critical problems, such as low resource recovery and difficult reinjection [7]. During the recharge process, water temperature, seepage velocity, and aquifer environmental characteristics affect the formation mechanism and development of clogging, which exists in physical, chemical, and biological forms, making it difficult to diagnose the causes of clogging [8]. Therefore, numerical simulations of geothermal reservoirs are a good tool for recoverable resources assessment and production scheme optimization to increase the operating lifetime of the system.For the geothermal resources evaluation, the usual method is a steady-state calculation method such as the surface heat flow method, reservoir volume method, and so on. It is very important to evaluate geothermal resources by considering the effects of geothermal development and the utilization process. For the geothermal developing process, the current method to assign the well spacing is over simplified as an empirical method is used [9]. In addition, the key operating parameters such as production rate, injection rate, and so on, all depend on the energy demand. However, the temperature distribution, fluid flow, stress field, and chemical process have a significant dynamic effect on the operating scheme at different stages. Numerical simulation technology is highly recommended to deal with these important problems to make the best use of geothermal resources. The numerical simulation method is a dynamic geothermal resource evaluation method that can be tracked and optimized [10, 11]. As a geothermal management tool, it can dynamically calculate the amount of geothermal recoverable resources in the process of development. It also can consider the heat supplement inside the geothermal system and the effect of reinjection on the geothermal recovery process.Numerical simulation technology of geothermal system can accurately describe the thermal recovery process and can accurately evaluate the geological resources, recoverable resources, economically recoverable resources, and green economically recoverable resources in the target block of geothermal resources, which provides guidance for project evaluation. Moreover, the production potential of geothermal wells can be predicted under different injection and production schemes, which gives the satisfied well layout and production scheme [12, 13]. To sum up, the development of integrated numerical simulation technology of geothermal system has very important strategic and practical significance for guiding the current development of hydrothermal geothermal resources.


## 1.2. The Status Analysis of the Numerical Simulation Method on Geothermal System

A lot of universities and institutes have developed simulation technologies on geothermal development, and the main commercial softwares are TOUGH2, OpenGeoSys, STARS, and FEHM [14–21]. Table 1 shows major numerical simulation projects conducted by some well-known companies in recent years.Table 1
Recent research on geothermal system simulation.
Institute and timeName of paper/ReportDescriptionNorth California power agency, 2010Integrated geothermal system simulator for the enhancement of geothermal steam production at the geysers [22]Tough2, stars and FEHM are used in reservoir; Thermo Flex is used in well-flowReykjavík university, landsvirkjun co. ltd., 2016Numerical modeling of the hágöngur geothermal reservoir in central Iceland [23]Coupled simulation by using TOUGH2 and PetraSimSupreme energy, 2016A reservoir simulation of the Maura Laboh field, Indonesia [24]Reservoir and well simulation by using TOUGH2 and PetraSimCarotenuto et al. conducted a literature review study on hydrothermal system simulation technology [25]. The results showed that TOUGH2 coupled with well simulation, could effectively deal with geothermal system simulation. Many coupled models were built and associated modeling was done by the software TOUGH2 to analyze hydrothermal or enhanced geothermal systems [26–29]. Fowler et al. compared CSMP-GEMS, OpenGeoSys-GEM, and TOUGHREACT by simulating a geothermal reservoir [30]; they concluded that all these three kinds of software can correctly simulate the process of fluid flow, heat transfer, and chemical reaction in a geothermal reservoir. However, TOUGHREACT is limited to solve the complex geological model because it cannot deal with unstructured grid models. On the other hand, OpenGeoSys-GEM could only simulate quite simple incompressible flow processes because of its dependence on the groundwater seepage model. CSMP-GEMS could handle incompressible flow problems with unstructured grids model. Butler and Enedy [31] conducted an integrated simulation by using TETRAD/TOUGH2/STAR coupled with ThermoFlex/AspenPlus/PEPSE in support of the Geysers geothermal project. Wang et al. analyzed the advantages and disadvantages of numerical simulation technology on Enhanced Geothermal Systems (EGS) and Hydro Geothermal System [32] and described the development significance of EGS. Zhang et al. used TOUGH2 software to study a geothermal system in the Subei area, China [33]. Ganguly et al. reviewed the numerical simulation method on an enhanced geothermal system [34]. The results showed that water/rock interaction or two-phase flow may be an important issue, which should be taken into consideration. The thermal behavior of rock during geothermal energy development has also been investigated by several researchers [35, 36].Through the review on geothermal simulation of the above domestic and foreign software, the calculation of physical parameters of working fluid and heat reservoir is complex and inaccurate, and the operability for engineering application is poor for reinjection optimization and simulation speed. Therefore, based on the thermal-hydraulic-mechanical-chemical coupling mechanism and the MPI parallel computing platform, this study developed and formed a numerical simulation software of geothermal resources, which realizes five functions: thermal dynamic prediction of thermal recovery process, optimization of injection and production system, recoverable resources assessment, optimization of development scheme design, and formulation of development technology policy.


## 2. Software Illustration and Geothermal System Modeling

For most geothermal reservoirs, our model allows consideration and tracking of the components H2O, NaCl (as a full compound, not as Na+ and Cl− ions), and any individual (or mixtures of more than one) from among the following 11 gases: CH4, C2H6, C3H8, n-C4H10, H2S, CO2, N2, O2, H2, i-C4H10, and CH3CH2OH. Obviously, consideration of a mixture of more than one constituent gas as a pseudo-component presupposes an assumption of constant gas composition. Halite phase, a solid phase created by the precipitation of NaCl, can also be considered. A complete description of the phase diagram of H2O and brine that cover the temperature range from 0°C to 2000°C and the pressure range up to 100 MPa. The physical chemical processes associated with the partitioning of mass components between phases. These include dissolution and exsolution of gases, dissolution and precipitation of salt, as well as the phase changes (vaporization and condensation) encountered in geothermal systems during production process.


### 2.1. Governing Equations

The mathematical and numerical framework starts from the integral form of mass and energy balance equations describing fluid and heat flow in a multiphase, multicomponent system:(1)ddt∫VnMκdV=∫ΓnFκ⋅ndA˜+∫VnqκdV,where Vn is an arbitrary subdomain for integration bounded by the close surface Гn; the quantity M, F, and q denote the accumulation term (mass or energy per volume), the flux term, and sink/source term, respectively. n is a normal vector on the subdomain surface pointing inward into the element; the superscript k represents one component.The general form of the mass accumulation term is(2)Mκ=∑ϕSβρβXβκ,where ϕ is the effective porosity, ρβ is the density of phase β, Sβ is the phase saturation, and Xβk is the mole fraction of component k in phase ß.The heat accumulation term consists of contributions from the rock matrix and liquid phases. It is given by the following equation;(3)Mκ=1−ϕρRCRT+ϕ∑Sβρβuβκ,where ρR and CR are the grain density and specific heat of the host rock, respectively, T is the temperature and uβ is the specific internal energy in phase ß.In the mass flux term,Fk is the advective mass flux summing over phases:(4)Fκ=∑βXβκFβ+Jβκ,where Fβ can be expressed by the multiphase extensions of the Darcy’s law:(5)Fβ=−kokrβρβμβ∇Pβ−ρβg,where ko is the absolute permeability, krβ is the relative permeability to phase β, µβ is the viscosity, Pβ is the pressure in the ß phase, and g is the gravity.In highly permeable porous media, the gas phase often exhibits flow that is turbulent and cannot be accurately described by Darcy’s law. In this regime, turbulent flow is described by the Forchheimer equation. For a low-permeable geothermal reservoir, it is necessary to account for inertial and gas slippage effects that are caused by gas flow in low-permeability media, in which the Brownian motion of the gas molecules is affected and inhibited by the very small pores of such media. The mass flux of the gas phase can be more complicated and is given by the relationship in following equation:(6)FG=−k1+bPGkr,GρGμG∇PG−ρGg.The termJβκ is the diffusive flux of component κ in the phase β, which can be estimated from(7)Jβκ=−ϕτ0τβρβDβκ∇Xβκ,where τ0 is intrinsic (medium-related) tortuosity, τβ is the mobile phase-related component of tortuosity, and Dβκ is multicomponent molecular diffusion coefficient of component κ in the phase η in the absence of a porous medium.The heat flux term accounts for the conduction and advection heat transfer and is given by the following equation:(8)Fκ=−λeff∇T+∑βhβFβ+εσΔT4,where hβ is the specific enthalpy of phase β, ε is the radiance emittance factor of the radiating surface, σ is the Stefan-Boltzmann constant, and λeff is the effective thermal conductivity.Heat conduction has a significant impact on the geothermal recharge process, so the calculation of the thermal conductivity of reservoir needs to be accurate. The previous calculation methods do not consider the influence of water saturation on thermal conductivity and do not comprehensively consider the thermal conductivity of adjacent grid points. The greater the heterogeneity of heat storage, the greater the error of the calculation results adopted by other software compared with the improved physical property calculation method.The harmonic average is introduced to calculate the effective thermal conductivity considering saturation, as shown in Figure2.(9)λ1=λR+ϕ1S1,AquλAqu,λ2=λR+ϕ2S2,AquλAqu,λ¯=λ1λ2d1λ1+d2λ2.Figure 2
Calculation diagram of thermal conductivity of adjacent grids.The phase diagram of H2O and its associated thermophysical properties in the vapor and liquid states are estimated from the fast parametric relationships of the IAPWS97 formulation and its more recent updates. These correlations provide the fastest (computationally) and most accurate estimates of the boundaries of the various regions of the H2O phase diagram, as well as all relevant properties of H2O: vapor pressure, density, viscosity, enthalpy, and thermal conductivity.The effect of temperature and pressure change on the porosity of the matrix is defined by the following standard exponential equation:(10)ϕ=ϕ0expαTT−T0+αpp−p0,where αT is the thermal expansivity of porous media (1/K) and αp is the pore compressibility (1/Pa), which can be a fixed number or a function of pressure. The permeability and porosity relationship is defined by the following general expression of Rutqvist and Tsang (2002) [26]:(11)k/k0=expγϕ/ϕ0−1,where γ is an empirical permeability reduction factor that ranges between 5 and 29 (5 is for soft unconsolidated media; 29 is for lithified, highly consolidated media). It should be noted that the above equations are simple and convenient to apply to matrix ϕ and κ changes when the changes in pressure and temperature are not very large.


### 2.2. Description of Geothermal System Simulation Software-SMG

The integral finite difference (IFD) method is used to discretize the continuum (1) in space, the advantage of IFD is not limited to regular grids, but is directly applicable to any irregular grids as long as the line connecting the centers of gravity between two grid blocks is perpendicular to their common areas. The resulting discretized equations are expressed in residual form, then are linearized by the Newton-Raphson method and are solved fully implicitly.The simulator for the multiphysics-coupling geothermal system (SMG) is developed by the Exploration and Production Research Institute, Sinopec. SMG is a fully implicit numerical code for the simulation of multiphase, multicomponent flow, and transport of mass and heat through complex porous and fractured media in geothermal systems, and represents the state of the art of geothermal simulators as it incorporates the most recent advances in the physics, numeric, and computing. SMG is programmed in the C++ language is based on the MPI architecture, and allows seamless application on a wide range of multiprocessor platforms running a wide variety of systems. It can be run on PCs, Workstations, and supercomputers under Windows/Unix/MacOS systems. The flow chart of the SMG simulation is shown in Figure3. It employs dynamic memory allocation, thus minimizing storage requirements. It follows the tenets of Object-Oriented Programming (OOP) and involves entirely new data structures that describe the objects upon which the code is based. Minimization of global variables and practically complete reliance on local ones minimizes inter-processor communication (the bane of parallel computing) within the MPI framework and, thus, results in a very fast multiprocessor application.Figure 3
Flow chart of SMG software.SMG can deal with multiple coupled process problems in support of geothermal systems. The software includes fluid flow, heat transfer, mechanical engineering, and chemical reaction. The temperature variations lead to thermal stress, which will also affect the reservoir temperature distribution and the values of reservoir permeability, then affect the fluid flow. On the other hand, fluid flow in geothermal reservoir will cause heat transfer, which will affect temperature distribution. In conclusion, the fully coupled geothermal system will change the permeability and conductivity of the reservoir during the long-term development, then affect the thermal energy output. On the one hand, by accurately describing the mechanism of flow and heat transfer in high-temperature fractures, this technology realizes the numerical simulation of geothermal system. On the other hand, the mechanism of multicomponent flow and mineral dissolution and precipitation is fully considered, which provides a means for the formulation of the reinjection strategy of hydrothermal storage and the prediction of scaling. In addition, the technology realizes MPI parallel computing under Linux and Windows systems and meets the needs of refined simulation of complex heat storage and large-scale development of geothermal engineering projects.


### 2.3. Model Verification

The international geothermal standard example is used for model validation. Table2 summarizes parameters used in the simulation. The problem concerns fluid flow and heat transfer along a one-dimensional horizontal rock with the length of 1 m. Mechanical and precipitation and gravity effect are unconsidered. The system contains single-phase water at initial temperature of 30°C, and then a constant water mass injection rate 0.1 kg/s with a colder temperature of 20°C is imposed at the inlet. The boundary of outlet is outflow and adiabatic.Table 2
Parameters and values for model verification.
ParameterValue and unitLength1 mRock thermal conductivity3.0 W/(m∙°C)Permeability2.96 × 10−2 DInitial temperature30°CInjection rate0.1 kg/sRock density2.6 × 103 kg/m3Rock heat capacity1000 J/(kg∙°C)Porosity0.3Initial pressure3.0 × 106 PaInjection temperature20°CFigures4 and 5 present simulation results of the temperature and pressure at different time vs position from outlet. It is observed that simulation results of SMG have a good agreement with analytical solutions.Figure 4
Temperature distribution at different times.Figure 5
Pressure distribution at different times.


## 2.1. Governing Equations

The mathematical and numerical framework starts from the integral form of mass and energy balance equations describing fluid and heat flow in a multiphase, multicomponent system:(1)ddt∫VnMκdV=∫ΓnFκ⋅ndA˜+∫VnqκdV,where Vn is an arbitrary subdomain for integration bounded by the close surface Гn; the quantity M, F, and q denote the accumulation term (mass or energy per volume), the flux term, and sink/source term, respectively. n is a normal vector on the subdomain surface pointing inward into the element; the superscript k represents one component.The general form of the mass accumulation term is(2)Mκ=∑ϕSβρβXβκ,where ϕ is the effective porosity, ρβ is the density of phase β, Sβ is the phase saturation, and Xβk is the mole fraction of component k in phase ß.The heat accumulation term consists of contributions from the rock matrix and liquid phases. It is given by the following equation;(3)Mκ=1−ϕρRCRT+ϕ∑Sβρβuβκ,where ρR and CR are the grain density and specific heat of the host rock, respectively, T is the temperature and uβ is the specific internal energy in phase ß.In the mass flux term,Fk is the advective mass flux summing over phases:(4)Fκ=∑βXβκFβ+Jβκ,where Fβ can be expressed by the multiphase extensions of the Darcy’s law:(5)Fβ=−kokrβρβμβ∇Pβ−ρβg,where ko is the absolute permeability, krβ is the relative permeability to phase β, µβ is the viscosity, Pβ is the pressure in the ß phase, and g is the gravity.In highly permeable porous media, the gas phase often exhibits flow that is turbulent and cannot be accurately described by Darcy’s law. In this regime, turbulent flow is described by the Forchheimer equation. For a low-permeable geothermal reservoir, it is necessary to account for inertial and gas slippage effects that are caused by gas flow in low-permeability media, in which the Brownian motion of the gas molecules is affected and inhibited by the very small pores of such media. The mass flux of the gas phase can be more complicated and is given by the relationship in following equation:(6)FG=−k1+bPGkr,GρGμG∇PG−ρGg.The termJβκ is the diffusive flux of component κ in the phase β, which can be estimated from(7)Jβκ=−ϕτ0τβρβDβκ∇Xβκ,where τ0 is intrinsic (medium-related) tortuosity, τβ is the mobile phase-related component of tortuosity, and Dβκ is multicomponent molecular diffusion coefficient of component κ in the phase η in the absence of a porous medium.The heat flux term accounts for the conduction and advection heat transfer and is given by the following equation:(8)Fκ=−λeff∇T+∑βhβFβ+εσΔT4,where hβ is the specific enthalpy of phase β, ε is the radiance emittance factor of the radiating surface, σ is the Stefan-Boltzmann constant, and λeff is the effective thermal conductivity.Heat conduction has a significant impact on the geothermal recharge process, so the calculation of the thermal conductivity of reservoir needs to be accurate. The previous calculation methods do not consider the influence of water saturation on thermal conductivity and do not comprehensively consider the thermal conductivity of adjacent grid points. The greater the heterogeneity of heat storage, the greater the error of the calculation results adopted by other software compared with the improved physical property calculation method.The harmonic average is introduced to calculate the effective thermal conductivity considering saturation, as shown in Figure2.(9)λ1=λR+ϕ1S1,AquλAqu,λ2=λR+ϕ2S2,AquλAqu,λ¯=λ1λ2d1λ1+d2λ2.Figure 2
Calculation diagram of thermal conductivity of adjacent grids.The phase diagram of H2O and its associated thermophysical properties in the vapor and liquid states are estimated from the fast parametric relationships of the IAPWS97 formulation and its more recent updates. These correlations provide the fastest (computationally) and most accurate estimates of the boundaries of the various regions of the H2O phase diagram, as well as all relevant properties of H2O: vapor pressure, density, viscosity, enthalpy, and thermal conductivity.The effect of temperature and pressure change on the porosity of the matrix is defined by the following standard exponential equation:(10)ϕ=ϕ0expαTT−T0+αpp−p0,where αT is the thermal expansivity of porous media (1/K) and αp is the pore compressibility (1/Pa), which can be a fixed number or a function of pressure. The permeability and porosity relationship is defined by the following general expression of Rutqvist and Tsang (2002) [26]:(11)k/k0=expγϕ/ϕ0−1,where γ is an empirical permeability reduction factor that ranges between 5 and 29 (5 is for soft unconsolidated media; 29 is for lithified, highly consolidated media). It should be noted that the above equations are simple and convenient to apply to matrix ϕ and κ changes when the changes in pressure and temperature are not very large.


## 2.2. Description of Geothermal System Simulation Software-SMG

The integral finite difference (IFD) method is used to discretize the continuum (1) in space, the advantage of IFD is not limited to regular grids, but is directly applicable to any irregular grids as long as the line connecting the centers of gravity between two grid blocks is perpendicular to their common areas. The resulting discretized equations are expressed in residual form, then are linearized by the Newton-Raphson method and are solved fully implicitly.The simulator for the multiphysics-coupling geothermal system (SMG) is developed by the Exploration and Production Research Institute, Sinopec. SMG is a fully implicit numerical code for the simulation of multiphase, multicomponent flow, and transport of mass and heat through complex porous and fractured media in geothermal systems, and represents the state of the art of geothermal simulators as it incorporates the most recent advances in the physics, numeric, and computing. SMG is programmed in the C++ language is based on the MPI architecture, and allows seamless application on a wide range of multiprocessor platforms running a wide variety of systems. It can be run on PCs, Workstations, and supercomputers under Windows/Unix/MacOS systems. The flow chart of the SMG simulation is shown in Figure3. It employs dynamic memory allocation, thus minimizing storage requirements. It follows the tenets of Object-Oriented Programming (OOP) and involves entirely new data structures that describe the objects upon which the code is based. Minimization of global variables and practically complete reliance on local ones minimizes inter-processor communication (the bane of parallel computing) within the MPI framework and, thus, results in a very fast multiprocessor application.Figure 3
Flow chart of SMG software.SMG can deal with multiple coupled process problems in support of geothermal systems. The software includes fluid flow, heat transfer, mechanical engineering, and chemical reaction. The temperature variations lead to thermal stress, which will also affect the reservoir temperature distribution and the values of reservoir permeability, then affect the fluid flow. On the other hand, fluid flow in geothermal reservoir will cause heat transfer, which will affect temperature distribution. In conclusion, the fully coupled geothermal system will change the permeability and conductivity of the reservoir during the long-term development, then affect the thermal energy output. On the one hand, by accurately describing the mechanism of flow and heat transfer in high-temperature fractures, this technology realizes the numerical simulation of geothermal system. On the other hand, the mechanism of multicomponent flow and mineral dissolution and precipitation is fully considered, which provides a means for the formulation of the reinjection strategy of hydrothermal storage and the prediction of scaling. In addition, the technology realizes MPI parallel computing under Linux and Windows systems and meets the needs of refined simulation of complex heat storage and large-scale development of geothermal engineering projects.


## 2.3. Model Verification

The international geothermal standard example is used for model validation. Table2 summarizes parameters used in the simulation. The problem concerns fluid flow and heat transfer along a one-dimensional horizontal rock with the length of 1 m. Mechanical and precipitation and gravity effect are unconsidered. The system contains single-phase water at initial temperature of 30°C, and then a constant water mass injection rate 0.1 kg/s with a colder temperature of 20°C is imposed at the inlet. The boundary of outlet is outflow and adiabatic.Table 2
Parameters and values for model verification.
ParameterValue and unitLength1 mRock thermal conductivity3.0 W/(m∙°C)Permeability2.96 × 10−2 DInitial temperature30°CInjection rate0.1 kg/sRock density2.6 × 103 kg/m3Rock heat capacity1000 J/(kg∙°C)Porosity0.3Initial pressure3.0 × 106 PaInjection temperature20°CFigures4 and 5 present simulation results of the temperature and pressure at different time vs position from outlet. It is observed that simulation results of SMG have a good agreement with analytical solutions.Figure 4
Temperature distribution at different times.Figure 5
Pressure distribution at different times.


## 3. Applications

### 3.1. Simulation Model in Laoling Area, Shandong Province, China

Based on the geothermal system in the Laoling area, Shandong Province, China, the influences of well spacing and injection temperature on reservoir temperature distribution were investigated. The numerical simulations of the three-dimensional, unsteady flow, heat transfer, and precipitation process in the computational domain were performed using the software SMG.The conceptual model is shown in Figure6 which includes the gap layer and geothermal reservoir. The related data are shown in Table 3. The depth of gap layer is 3000 m, and the thickness of geothermal reservoir is 180 m. The permeability of gap layer along X, Y, and Z directions are 5 × 10−5 mD, 5 × 10−5 mD, and 1.25 × 10−5 mD, respectively. The permeability of the geothermal reservoir along X, Y, and Z directions are 0.2 D, 0.2 D, and 0.02 mD, respectively. The porosities of the gap layer and geothermal reservoir are 1% and 28.9%. All the above data are collected from field. In order to reduce the computation consumption, half of the physical model is simulated. The reservoir model is large enough, so the boundary conditions around the model are impermeable and adiabatic. The top boundary of the reservoir is impermeable and adiabatic and the bottom boundary is set as constant heat flux, 72 W/m/K, which considers the heat supply from the bottom rock. The mass flow rate boundary condition is applied at the reservoir inlet and the outflow boundary condition is used at the reservoir outlet. All other boundaries are treated as no-flow and adiabatic boundaries.Figure 6
Schematic of conceptual model.Table 3
Parameters for Laoling area geothermal reservoir simulation.
ParameterValue and unitReservoir length3200 mReservoir width1000 mReservoir depth180 mRock thermal conductivity2.6 W/(m∙°C)Production/Injection rate104.6 m3/hRock density2.6 × 103 kg/m3Rock heat capacity880 J/(kg∙°C)Salinity6398 mg/LInitial average pressure31.2 MPaInitial rock temperature48 ∼52°CFor transient simulations, the initial conditions have an important effect on the initial transient. The temperature and pressure differences in the reservoir increase with depth, so the analysis of the transient behavior with geothermal reservoir requires appropriate initial reservoir conditions. Thus, a steady-state simulation of the reservoir without any water injection or production is used to determine the initial stage. The reservoir is filled with water with the pressure set at 31.2 MPa and the temperature at the top surface fixed at 48°C, the salinity 6398 mg/L. The other boundary conditions are as described in the above. The initial stationary reservoir conditions are obtained with the residuals all converged to 1 × 10−12.In order to capture the detailed fluid flow and heat transfer process near the well, meshes near the local well is refined, as shown in Figure7. Mesh independence tests are conducted with the numbers of grids are 204600, 336400, and 560000, respectively. The number of grids is set to be 336400 based on the simulation results comparison.Figure 7
Mesh of simulation model (meshes near the wells are refined).


### 3.2. Effect of Well Distance on Reservoir Temperature Distribution

Figure8 shows the relationship between temperature at the production well and the time under different well spacing conditions. As time increases, the temperature of production decreases because of hot water production and cold water (∼30°C) injection. At the primary stage of production, the water temperature decreases quickly, then the temperature of production maintains constant due to the stable heat transfer area and fluid flow in the reservoir. It is obvious that the temperature of production decreases with the well spacing increases.Figure 8
Relationship of water temperature at the production well and time.The temperature distribution of the reservoir, which has been developed for 50 years, as shown in Figure9. The temperature influence area of production well increases with the increment of well spacing. According to Figures 8 and 9, the total amount of heat energy production for 50 years is increasing with the well spacing increment.Figure 9
Temperature distribution for 50 years of development.


### 3.3. Effect of Injection Temperature on Reservoir Temperature Distribution

Under the condition of 800 m well spacing, this paper investigated the effect of injection temperature (28°C, 32°C, and 40°C) on the temperature of production water. As shown in Figure10, during the 45 years of production, the lower the injection temperature, the lower the production temperature. According to Figure 11, during 120 days of production, the production temperature is the same under different injection temperatures.Figure 10
Relationship of production temperature at the bottom hole and time.Figure 11
Relationship of production temperature at the bottom hole and time (120 days).


### 3.4. Evaluation of Regional Geothermal Dynamic Recoverable Resources in Laoling

As noted, geothermal dynamic recoverable resources are more useful for geothermal process assessment. Therefore, a geothermal recoverable resources evaluation method based on numerical simulation is proposed, as described in Figure12.Figure 12
Evaluation method of regional geothermal dynamic recoverable resources.Figure 13
Temperature distribution of geothermal reservoir of the Guantao formation in Laoling after 100 years heating.Most geothermal well groups of Guantao Formation in Laoling have one injection and one production well, with a spacing of 500 m. The hydrothermal geothermal resources in this area are usually used for heating, which is used for heating for 4 months in one year and well shutdown for recovery in 8 months. Generally, the evaluation period of geothermal recoverable resources is set at 100 years. Therefore, this paper uses SMG software and the above model to evaluate the geothermal dynamic recoverable resources of the Guantao Formation in the Laoling area.It can be seen from the figure that after 100 years of heating, the influence range of reinjection wells in the Guantao Formation reservoir is 550 m × 760 m. Therefore, when the well spacing is 500 m, the control range of one injection and one injection well group is 550 m × 760 m. Based on the above evaluation method, the results of the reasonable well groups and total geothermal dynamic recoverable resources in this area, and the available heating area are listed in Table4.Table 4
Regional geothermal dynamic recoverable resources of the Guantao Formation in Laoling.
ParameterValue and unitReservoir area38.59 × 106 m2Influence area of a well group0.418 × 106 m2Recoverable resources of well group0.328 × 1016JHeat extraction rate of well group3.13 MWRequired heat flux for house heating40 W/m2Influence scope550 m × 760 mNumber of well groups92Regional recoverable resources29.256 × 1016 JRegional heat extraction rate287.61 MWRegional heating area719.03 × 106 m2


## 3.1. Simulation Model in Laoling Area, Shandong Province, China

Based on the geothermal system in the Laoling area, Shandong Province, China, the influences of well spacing and injection temperature on reservoir temperature distribution were investigated. The numerical simulations of the three-dimensional, unsteady flow, heat transfer, and precipitation process in the computational domain were performed using the software SMG.The conceptual model is shown in Figure6 which includes the gap layer and geothermal reservoir. The related data are shown in Table 3. The depth of gap layer is 3000 m, and the thickness of geothermal reservoir is 180 m. The permeability of gap layer along X, Y, and Z directions are 5 × 10−5 mD, 5 × 10−5 mD, and 1.25 × 10−5 mD, respectively. The permeability of the geothermal reservoir along X, Y, and Z directions are 0.2 D, 0.2 D, and 0.02 mD, respectively. The porosities of the gap layer and geothermal reservoir are 1% and 28.9%. All the above data are collected from field. In order to reduce the computation consumption, half of the physical model is simulated. The reservoir model is large enough, so the boundary conditions around the model are impermeable and adiabatic. The top boundary of the reservoir is impermeable and adiabatic and the bottom boundary is set as constant heat flux, 72 W/m/K, which considers the heat supply from the bottom rock. The mass flow rate boundary condition is applied at the reservoir inlet and the outflow boundary condition is used at the reservoir outlet. All other boundaries are treated as no-flow and adiabatic boundaries.Figure 6
Schematic of conceptual model.Table 3
Parameters for Laoling area geothermal reservoir simulation.
ParameterValue and unitReservoir length3200 mReservoir width1000 mReservoir depth180 mRock thermal conductivity2.6 W/(m∙°C)Production/Injection rate104.6 m3/hRock density2.6 × 103 kg/m3Rock heat capacity880 J/(kg∙°C)Salinity6398 mg/LInitial average pressure31.2 MPaInitial rock temperature48 ∼52°CFor transient simulations, the initial conditions have an important effect on the initial transient. The temperature and pressure differences in the reservoir increase with depth, so the analysis of the transient behavior with geothermal reservoir requires appropriate initial reservoir conditions. Thus, a steady-state simulation of the reservoir without any water injection or production is used to determine the initial stage. The reservoir is filled with water with the pressure set at 31.2 MPa and the temperature at the top surface fixed at 48°C, the salinity 6398 mg/L. The other boundary conditions are as described in the above. The initial stationary reservoir conditions are obtained with the residuals all converged to 1 × 10−12.In order to capture the detailed fluid flow and heat transfer process near the well, meshes near the local well is refined, as shown in Figure7. Mesh independence tests are conducted with the numbers of grids are 204600, 336400, and 560000, respectively. The number of grids is set to be 336400 based on the simulation results comparison.Figure 7
Mesh of simulation model (meshes near the wells are refined).


## 3.2. Effect of Well Distance on Reservoir Temperature Distribution

Figure8 shows the relationship between temperature at the production well and the time under different well spacing conditions. As time increases, the temperature of production decreases because of hot water production and cold water (∼30°C) injection. At the primary stage of production, the water temperature decreases quickly, then the temperature of production maintains constant due to the stable heat transfer area and fluid flow in the reservoir. It is obvious that the temperature of production decreases with the well spacing increases.Figure 8
Relationship of water temperature at the production well and time.The temperature distribution of the reservoir, which has been developed for 50 years, as shown in Figure9. The temperature influence area of production well increases with the increment of well spacing. According to Figures 8 and 9, the total amount of heat energy production for 50 years is increasing with the well spacing increment.Figure 9
Temperature distribution for 50 years of development.


## 3.3. Effect of Injection Temperature on Reservoir Temperature Distribution

Under the condition of 800 m well spacing, this paper investigated the effect of injection temperature (28°C, 32°C, and 40°C) on the temperature of production water. As shown in Figure10, during the 45 years of production, the lower the injection temperature, the lower the production temperature. According to Figure 11, during 120 days of production, the production temperature is the same under different injection temperatures.Figure 10
Relationship of production temperature at the bottom hole and time.Figure 11
Relationship of production temperature at the bottom hole and time (120 days).


## 3.4. Evaluation of Regional Geothermal Dynamic Recoverable Resources in Laoling

As noted, geothermal dynamic recoverable resources are more useful for geothermal process assessment. Therefore, a geothermal recoverable resources evaluation method based on numerical simulation is proposed, as described in Figure12.Figure 12
Evaluation method of regional geothermal dynamic recoverable resources.Figure 13
Temperature distribution of geothermal reservoir of the Guantao formation in Laoling after 100 years heating.Most geothermal well groups of Guantao Formation in Laoling have one injection and one production well, with a spacing of 500 m. The hydrothermal geothermal resources in this area are usually used for heating, which is used for heating for 4 months in one year and well shutdown for recovery in 8 months. Generally, the evaluation period of geothermal recoverable resources is set at 100 years. Therefore, this paper uses SMG software and the above model to evaluate the geothermal dynamic recoverable resources of the Guantao Formation in the Laoling area.It can be seen from the figure that after 100 years of heating, the influence range of reinjection wells in the Guantao Formation reservoir is 550 m × 760 m. Therefore, when the well spacing is 500 m, the control range of one injection and one injection well group is 550 m × 760 m. Based on the above evaluation method, the results of the reasonable well groups and total geothermal dynamic recoverable resources in this area, and the available heating area are listed in Table4.Table 4
Regional geothermal dynamic recoverable resources of the Guantao Formation in Laoling.
ParameterValue and unitReservoir area38.59 × 106 m2Influence area of a well group0.418 × 106 m2Recoverable resources of well group0.328 × 1016JHeat extraction rate of well group3.13 MWRequired heat flux for house heating40 W/m2Influence scope550 m × 760 mNumber of well groups92Regional recoverable resources29.256 × 1016 JRegional heat extraction rate287.61 MWRegional heating area719.03 × 106 m2


## 4. Conclusions

An efficient modeling approach, named SMG which can simulate fluid flow and heat transfer in geothermal reservoirs was developed based on an optimized multiphysics coupling mathematical model and MPI architecture. Compared with other geothermal modeling approaches, it considers comprehensive THMC coupling process and increases computational efficiency by simplifying Jacobian matrix size based on a parallel method while keeping sufficient accuracy. Our numerical results show that this software captures all the key patterns of fluid flow and heat transfer in geothermal reservoirs. Moreover, a geothermal recoverable resources evaluation method based on numerical simulation is proposed. The geothermal numerical simulator SMG and the proposed evaluation method are good tool for geothermal resources evaluation and optimization of the geothermal development scheme for projects, as shown in the example of Guantao formation in Laoling area in this paper.



---

*Source: 3398965-2022-05-25.xml*