Simulation of thermal breakthrough factors affecting carbonate geothermalto-well systems

Jia-xing Sun, Gao-fan Yue, Wei Zhang

1 Langfang Comprehensive Survey Center of Natural Resources, China Geology Survey, Langfang 065000, Hebei Province, China.2 Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences, Shijiazhuang 050061, China.3 Technology Innovation Center of Geothermal & Hot Dry Rock Exploration and Development, Ministry of Natural Resources, Shijiazhuang 050061, China.4 China University of Petroleum (Beijing), Beijing 102249, China.

Abstract: Fractures play a pivotal role in carbonate thermal storage systems, serving as primary hydraulic conductivity channels that significantly influence thermal breakthrough times and heat extraction efficiency in geothermal-to-well systems.Their impact is critical for well placement and system life prediction.This paper focuses on a geothermal-to-well system within the carbonate reservoir of the Wumishan formation in the Rongcheng geothermal field, Xiong'an new area.It employs a combination of field tests and numerical simulations to determine the permeability of the reservoir and the evolution of fractures between wells.It also examines the influence of fracture width and roughness coefficient on the seepage and temperature fields under various injection scenarios and predicts thermal breakthrough times for production wells.The results show: Higher permeability is observed near well D16 compared to well D22 within the studied geothermal-to-well systems.Wider fractures between wells result in faster temperature decline in production wells.Lower injection flow rates lead to slower temperature reduction in injection wells.The use of roughness coefficients minimizes temperature variations in production wells.This study not only offers guidance for the development and utilization of the geothermal well system, but also contributes to a deeper understanding of the groundwater seepage and heat transfer process influenced by fractures.

Keywords: Geothermal recharge; Influencing factor; Thermal breakthrough; Seepage field; Temperature field

Introduction

Geothermal energy is a clean and sustainable energy sources known for its minimal greenhouse gas emissions and the enhancement of energy efficiency.Its applications include heating, bathing,agriculture, industry and aquaculture crop drying and mineral water production (Moya et al.2018;Zhu et al.2015).Giving these advantages, the development of geothermal energy carries significant importance.

In China, particularly in regions like Xiong'an,carbonate rock thermal storage is a prominent component of geothermal energy (Wang et al.2019).The geological characteristics of formations like the Wumishan Formation within the Jixian system offer shallow burial, high temperatures, ample water resources, and ease of extraction and irrigation (Wang et al.2019).With robust policy support and technical expertise, efficient large-scale geothermal harvesting and irrigation have become a reality (Pang et al.2017; Long et al.2021).To ensure the sustainable development of geothermal resources, geothermal-to-well mining technology has been widely utilized.Understanding the factors that influence thermal breakthrough in geothermal-to-well systems is crucial for the sustainable development of geothermal resources.

Previous research in the Xiong'an new area has covered various phases of geothermal system development.It has evolved from investigating the causal mechanism of geothermal fields to examining extraction and irrigation during the exploitation phase and applying numerical simulation to predict the potential of well systems (Chen et al.1990; He and Li, 2005).Cheng et al.(2011) applied numerical simulation methods to model the temperature fields of the extraction and irrigation geothermal-to-well and determine optimal extraction and irrigation parameters such as distances and recharge temperatures.The transition from twodimensional to three-dimensional model has allowed for a more comprehensive analysis of groundwater flow, temperature and chemical field evolution within geothermal systems, providing insights into the impact of mining and irrigation on the lifespan of production wells (Luo, 2018).In terms of geothermal fluid flow-heat coupling,domestic researchers have developed heat transfer coupling models for fractured rock bodies, in which the heat conduction equation of the rock bodies are discretized by using the Galerkin method with weighted residuals, and the temperature field distribution of fissured rock bodies can be determined by three-dimensional numerical simulations (Zhang and Xu, 2011).Shan et al.(2019) employed a two-dimensional random fracture model to simulate geothermal-to-wells systems, integrating temperature and seepage fields.Their analysis delved into the effects on system capacity and lifetime.Mirmahdi et al.(2019) developed a three-dimensional thermal storage model using Tough2 to replicate geothermal fluid flow and heat transfer in Sabaran geothermal system, located in Northwest Iran.Furthermore,Wang et al.(2019) examined coupled fluid flow and heat transfer processes within geothermal reservoirs of geothermal fields.Their findings support that one-dimensional geothermal well models can be simplified for resource assessment and optimization of injection and extraction strategies.Despite these advancements in geothermal well simulation research, there remains a gap in studies on well development in the Xiong'an new area,with relatively simple parameters.The recent adoption of karst fracture hot storage reinjection technology (Dong et al.2021) has introduced a new dimension to the field.Fractures play a vital role in facilitating the conduction and convection of geothermal fluids.Carbonate thermal storage permeability is influenced by both rock medium and fracture characteristics, however, previous well-to-well studies have often neglected the influence of fractures, especially inter-well fractures.

In this paper, we investigate the synthesis of dual pore-fracture medium within carbonate thermal storage.Our study focuses on the geothermalto-well system, specially wells D16 and D22,located in Xiong'an new area.The objective is to determine the permeability of the reservoir's pore medium through on-site pumping tests.Subsequently, a multi-field coupled model is established including inter-well fractures.By combining the data from on-site pumping tests to determine the permeability of the reservoir's pore medium, we create a multi-field coupling model that accounts for the influence of inter-well fractures.Our analysis includes an examination of changes in the seepage field and temperature field in the presence of fractures, enabling us to predict the thermal breakthrough of the production wells.This research seeks to identify the key factors that impact the thermal breakthrough in geothermal-to-well systems within carbonate thermal reservoirs.

1 Study area

The study area is situated in the northern region of the Jizhong depression (III level) within the Bohai North China Basin (II level), located on the Sino Korean Paraplatform (I level).On the structural units, the study area encompasses several geographical features, including the Gaoyang uplift,Rongcheng uplift, Niutuozhen uplift, Baxian Sag,Raoyang Sag, Baoding Sag and Xushui Sag.These features collectively create a complex landscape with a pattern of "convex, broken, concave", featuring multi-convex, multi-concave, and concaveconvex interspersed distribution (Chen et al.1990;Lu et al.2019).The North China Plain exhibits a geothermal crust structure typical of a cold crust overlying a hot mantle (Lin et al.2022).The study area is also surrounded by significant fractures,including the Rongcheng Fracture, Xushui Fracture, Niu Nan Fracture and others.These fractures,in conjunction with tectonic deformation, serve as crucial conduits for thermal storage and transportation (Dai et al.2019).In terms of geological strata,the study area's sequence from youngest to oldest comprises the Quaternary, Neoproterozoic, Paleoproterozoic, Jixian, Changcheng, and Paleoproterozoic basement layers.The primary focus of this paper is on the Jixian thermal reservoir.The tectonic context of the study area is illustrated in Fig.1.

Fig.1 Tectonic background map of Xiong'an new area

The study site is located in the Rongcheng uplift and encompasses wells D16 and D22, which are approximately 200 meters apart.Well D16 functions as a production well, while well D22 serves as an injection well.Based on drilling information,the Rongcheng Fracture is located 7.5 kilometers to the east of drill hole D16, while the Xushui Fracture is found 3.5 kilometers to the southwest.The Rongcheng Fault exhibits a NNE strike, eastward inclination, and an approximately 45° dips.It acts as a growth fault controlling the development of the Neoproterozoic System.The Xushui Fault is a significant, deep, and long-term active fault with nearly east-west strike, southward inclination, and dips of around 45°.The heat transfer channels are influenced by regional active faults and concealed deep fractures, and large tectonic fractures impact tectonic cleavage (Tang, 2020; Hong et al.2019;Wei et al.2022).The site is characterized by the presence of regional active fractures, necessitating an investigation into the heat transfer characteristics of these fractures.Well D16's drilling encountered various strata, including the Quaternary,Neoproterozoic Minghuazhen Formation, Wumishan Mountain Formation, Yangzhuang Formation and Gaoyuzhuang Formation (undisclosed).Among these, the Wumishan Formation is found at a burial depth of 1,986.0 m and has a layer thickness of 1,006 m.It primarily consists of white,grayish-white, dark-gray, grayish-black thickbedded dolomite, muddy dolomite, and dolomite with flint strips.The upper part is characterized by a high degree of weathering and a prevalence of fractures.The wellhead water temperature measures around 69°C, with a water flow rate of gushing water 201.7 m3/h.Well D22 reaches a depth of 3,517.18 m, with the Wumishan Formation at a depth of 1,961 m.The lithology is mainly gray,grayish-white, and brownish-red dolomite, with well-developed fractures.Details regarding the Gaoyuzhuang Formation are not provided.The outlet water temperature from well D22 is 68°C.Notably, these two wells are positioned 200 m apart, and specific information about them is provided in Table 1.

The Jixian Formation thermal reservoir is classified as a deep-circulation tectonic rift type, charac-terized by its richness in thermal storage resources(Ma et al.2020).The area contains substantial recoverable resources of geothermal fluid, estimated at 401.77×106m3/a (Zhu et al.2022).

Table 1 Basic information on geothermal wells

Within the study area, fractures are classified into three main types: Rock-forming fractures,tectonic fractures and weathering fractures.Among these, tectonic fractures are the most widely distributed.A combination of data and logging information reveals that fractures within the site predominantly strike in the east-west and near north-east directions.At greater depths, the dominant fracture orientations are in the range of 125°-180° and 234°-300°, with the dip angles varying between 50° and 90° (Tang, 2020; Ma et al.2020).

2 Model construction

2.1 Conceptual model

A percolation-heat transfer model is developed based on the following assumptions:

(1) Simplification of the wellbore into a onedimensional line cell.

(2) Simplification of the reservoir as a porous,continuous medium with non-homogeneous anisotropy.

(3) Exclusion of the effects of thermal radiation.

(4) Omission of considerations regarding the effect of gravity.

Considering that the model's scope covers the thermal reservoir of the Wumishan Formation,with a layer thickness of 1,000 m in this area, an expansion of the model's height to 1,300 m is implemented to minimize boundary effects on simulation results.The model is configured as a rectangle shape with dimensions of 800 m in length,800 m in width, and 1,300 m in height.Wells are spaced 200 m apart.The model is initially divided into left and right sections based on logging data from production and injection wells, which facilitate subsequent parameter inversion analysis.For the left part corresponding to well D16, the model is further subdivided into two layers with thicknesses of 613.5 m and 686.5 m respectively.The permeability for the right side, represented by well D22, is interpolated accordingly.It's important to note that this stratification is introduced solely for computational and assignment convenience; in reality, the strata on both sides of the model are continuous.

Initial values for various parameters are determined based on well logging data and pumping test data.The initial pressure values are calculated using the surface pressure and the relevant formula:

Where: P0represents the atmospheric pressure and has a value is 1.013×105Pa,ρdenotes the density of water and has a value of 1,000 kg/m3，z is the grid elevation value, and h represents the water level and depth, with a value of -980 m.

Since the boundary area is significantly larger than the flow influence area, the upper, lower and surrounding boundaries are set as no-flow and thermally insulated boundaries.The model's visual representation is shown in Fig.2.

Fig.2 Schematic diagram of the model

2.2 Mathematical model and coupling relationship

(1) Flow equation Fluid flow in a thermal reservoir can be described by utilizing Darcy's law (Sun et al.2020), with the mathematical expression:

Where: εprepresents the porosity of the rock matrix, dimensionless;ρis the density of water,kg/m3;uis the velocity of fluid percolation within the matrix, m/s;Qmrepresents sources and sinks of seepage;kstands for permeability of bedrock, m2.

The conservation equation for fluid flow in the fracture is:

Where:dfrepresents the fracture opening, mm;kfis the fracture permeability, m2;efsignifies the volumetric strain at the fracture surface.Qfrepresents exchange terms between matrix and fracture,m3/s; τ denotes derivation along the tangential direction of the fracture.

(2) Heat transfer equation

Heat transfer in thermal storage follows the law of conservation of energy, including heat transfer in the rock matrix and heat transfer processes in fractures, the heat transfer equation in the rock matrix is given in Equation 5, the heat transfer equation in the fracture for the transient case is as follows (Equation 6).

Where:Csis the specific heat capacity of the matrix rock, J/(m3·K); ρsrepresents density of matrix rock, kg/m3;Tsis the substrate rock temperature, K; λsstands for thermal conductivity of matrix rock, W/(m·K);Wsignifies heat source term, W/m2.

Where:Tfrepresents fluid temperature in fracture, K; ρfis the fluid density, kg/m3;Cfsignifies the specific heat capacity of a fluid, J/ (m3·K);λfstands for the thermal conductivity of fluid, W/(m·K);hdenotes the heat transfer coefficient, W/(m2·K);ufrepresents the fluid velocity inside the fracture, m/s.

(3) Coupling relationship

This paper adopts a bidirectional coupling of flow and heat transfer: (1) The flow rates in the energy equation are derived from calculations of the Darcy equation; (2) Fluid properties are affected by temperature, including dynamic viscosity,heat capacity, density, and thermal conductivity of water (Qu et al.2017; Zhang et al.2018) (Equations 7-10).

The dynamic viscosity of water can be written as:

The heat capacity of water can be written as:

The density of water can be written as:

The thermal conductivity of water can be written as:

2.3 Numerical model

In this study, COMSOL Multiphysics software is used to solve the coupled equations described above, utilizing Darcy's law and the porous media heat transfer interface.

The model was discretized using a free tetrahedral mesh.Grid refinement was applied to ensure accuracy, particularly in the calibration of well D16, resulting in a well section dissected into a fixed cell number of 60.This process generated a total of 16,482 tetrahedral cells.Similarly, for well D22, grid refinement was implemented with a fixed cell number of 60, resulting in a total of 16,590 tetrahedral cells.The study involved transient analysis.

The matrix heat transfer of the model is calculated using Equation 5, the heat transfer in the fracture is calculated using Equation 6, and the fluid flow in the fracture is calculated using Equation 4.

The geological characteristics of the site indicate a rifted structure with predominant sub-eastwest and sub-north-east orientations, trending southeast.The dips within this area vary from 50°to 90°.Based on the logging data that provide information on fracture strike, inclination and dip,the interwell fracture surface was defined with a NE-SW strike, an inclination of SE and a dip of 60°.The injection well has a diameter of 24.45 cm,while the production well has a diameter of 17.78 cm.Refined grid discretization was applied to the production and injection well, resulted in a total of 17,013 tetrahedral cells, 2,956 triangular cells, 359 edge cells, and 45 vertex cells (Fig.3).

The output step size is set to 1 year, and the total computational step size is 50 years.The specific values of each parameter utilized in the model are detailed in Table 2.

2.4 Research program

Fig.3 Mesh sectioning case diagram

(1) Inversion of thermal storage permeability The model is constructed with parameters as described above, with well D16 as the production well and well D22 as the injection well.The placement of injection wells was considered in line with the effect of natural flow direction (Sun et al.2022).In the study area, the natural hydraulic gradient between production well and injection well is not significant, thus, upstream and downstream impacts are not considered in this case.The initial permeability was interpolated from logging data obtained from production well.Additionally,production well flow parameters were interpolated using data from pumping tests.The actual pumping flow is depicted in Fig.4.

The model was run to generate simulation results, which were subsequently fitted to the actual pumping test data.The fitting results are presented in Figs 5 and 6.By manually adjusting the permeability parameter, the fitting results consistently converged with the pumping test data.As demonstrated in Fig.7, a scatter plot of the actual water level drawdown and the simulated water level drawdown, the R2value of the linear fit between the two is 0.9188.This significant correlation indicates a high degree of linear fit, confirming that the simulation results align closely with the pumping test results.This high correlation underscores the reliability of the inversion model.

(2) Impact of fracture and production programs on geothermal systems

During the actual production process, various factors influence the temperature and seepage fields of geothermal systems.These factors includefracture characteristics (width, shape, size, and roughness), injection fluid velocity, injection temperature, well distance, and production time, as supported by previous research (Shan et al.2020;Xiao et al.2021; Li et al.2016; Sanaz et al.2015).In this study, we specifically focused on three sensitive parameters: Fracture width (df), roughness coefficient (ff), and injection flow rate.

Table 2 Model parameter values

Fig.4 Pumping flow conditions

Fig.5 Well D16 fitting results

Fig.6 Well D22 fitting results

Fig.7 Linear fit of actual water level drop to simulated water level drop (42 data sets)

To assess the impact of each parameter on the temperature of the production well, and to conduct a sensitivity analysis of the relevant parameters, we assigned specific values to each sensitive parameter, as shown in Table 3.

Table 3 Values for parameter sensitivity analysis

The main purpose of this study is to analyze variations in the temperature field and seepage field in production wells, with a specific focus on understanding the thermal breakthrough phenomena.In this context, thermal breakthrough is defined as a 2°C drop in production well temperature (Liu et al.2019).Furthermore, the study aims to investigate how the above-mentioned sensitive parameters, such as fracture characteristics and production parameters, influence the temperature variations in production wells.

3 Results

3.1 Permeability distribution

The permeability results obtained from the inversion fitting analysis in the first part of section 2.4 are as follows: On the left side of the model, where well D16 is located, the permeability ranges from 3.14e-16m2to 8.07e-13m2.On the right side of the model, where D22 is located, the permeability ranges from 4e-17m2to 1.16e-12m2.

Comparing these results with the logging data from well D16 reveals that the overall average permeability to the left of well D16 exceeds that of well D16 itself.The lithology analysis of well D16 indicates the presence of dolomite, tuff and dolomitic sandstone, suggesting a high overall dolomite content and extensive fracture development in the D16 area.The interpretation of seismic profiles in the area (He et al.2018) further supports these findings.The Rongcheng uplift exhibits a significant presence of dolomites, and the area's intense tectonic activity has resulted in numerous large faults, providing a conductive environment for fracture development.Analyzing the rock cores from well D16 and D22 affirms these observations.The core of well D16 displays increased fragmentation, high weathering, and well-developed fractures, reinforcing the conclusion that the left area is more permeable than the well section.Conversely,well D22 exhibits less developed fractures and slows signs of amphibolite intrusion.

3.2 Impact of fractures on the distribution of seepage and temperature fields

For the simulation, a fracture width of 3 mm, fracture roughness 1.6, and an injection flow rate of 165 m3/h were considered.The simulation results indicate that the injection of low-temperature water into the production well initiates a leftward movement of the cold frontal surface, forming a fanshaped pattern from the end of the injection well to the end of the production well.Analyzing the changes in temperature and seepage fields over 1 year, 5 years, 20 years, and 50 years (Fig.8), the seepage rate between the production well on the left and the injection well on the right generally accelerates over time, particularly in the vicinity of the two wells.Although the narrow fracture surface might not exhibit an obvious change in the seepage field, there are clear seepage channels between the two wells.

Regarding temperature, a gradual decrease is observed in the production well for about 10 years after the start of production.Subsequently, the temperature decreases at a faster rate, leading to thermal breakthrough, defined as a 2°C temperature drop, occurring in the production well in year 16.Fig.9 illustrates the temperature changes in the production well over a 50-year period.

3.3 Sensitivity analysis of parameters

The results, illustrated in Figs 10 to 12, indicate that fracture width, roughness coefficient, and injection flow rate exert varying influence on the temperature of the production well.Notably, the impact of the roughness coefficient on the temperature of the production well is marginal, making it largely negligible.Fracture width emerges as the most influential parameter, showing a substantial effect on temperature in production well.A wider fracture corresponds to a more rapid decline in production well temperature, accelerating the onset of thermal breakthrough and consequently reducing the lifespan of the production well.Additionally, a higher injection flow rate correlates with an increased likelihood of achieving thermal breakthrough and amplifying temperature change in the production well.

Fig.8 Variation of seepage and temperature fields after 1 year (a), 5 years (b), 20 years (c), 50 years (d) of production

Fig.9 Plot of temperature changes over 50 years in production wells

In practical applications, optimizing the injection flow rate becomes crucial.Furthermore, using well logging data, combined with strategies such as fracturing for production enhancement, allows for targeted selection of areas characterized by larger fracture widths.

4 Conclusion

The lithology surrounding wells D16 and D22 is predominantly dolomite, exhibiting relatively high degree of weathering and extensive fracture development.The permeability analysis reveals that well D16's perimeter has a greater permeability range (3.14e-16m2to 8.07e-13m2) than compared to well D22's perimeter (4e-17m2to 1.16e-12m2).The average permeability around well D16 surpasses that of well D22, forming a discernible seepage channel within the well-to-well system formed by the two wells.

Formation parameter such as fracture width,fracture roughness, and fluid injection rate significantly influence the temperature and seepage fields of the geothermal system.Fracture width and injection rate are sensitive factors affecting thermal breakthrough, while the roughness coefficient exhibits negligible sensitivity.Wider fractures and higher injection rate correlate with accelerated temperature drops in the production well and faster achievement of thermal breakthrough.

Over time and with the injection of cold water recharge, fluid seepage in the hot storage fissure undergoes dynamic changes.The cold front within the hot storage layer advances ahead of the fracture, shortening the heat exchange duration between fracture water and the rock matrix.

Fig.10 Effect of different fracture widths on production well temperature

Fig.11 Effect of different roughness coefficients on production well temperature

Fig.12 Effect of different injection flow rates on temperature of production well

Acknowledgements

This study was supported by basic research project of Chinese Academy of Geological Sciences (No.YK202309).Special thanks are due to the reviewers and editors of this journal for their valuable suggestions and revisions of the manuscript.