# Study on Preliminarily Estimating Performance of Elementary Deep Underground Engineering Structures in Future Large-Scale Heat Mining Projects

Authors: Min Tang; Hong Li; Chun’an Tang
Journal: Geofluids
 (2019)

---

## Abstract

Geothermal energy will become an important part of energy in the future because of its advantages in source stability, sustainability, and potential high utilization ratio. In particular, the development and utilization of deep geothermal energy from HDR have gradually attracted people’s attention. Aiming at solution to the bottleneck of EGS-D, a new EGS-E based on excavation technology is proposed. In this paper, a concise and direct method for estimating the early performance of this disruptive and innovative geothermal development scheme is established as a viable alternative to supercomputing for the subsequent quantitative research of the corresponding relationship between a typical deep engineering structure and its heat extraction efficiency. Firstly, the effects of the fixed temperature at a tunnel wall, the radius of a tunnel, and the rock type on the annual heat extraction rate of the tunnel are studied based on the analytical solution of a one-dimensional radial plane problem of the transient heat conduction through high-temperature surrounding rock to the tunnel wall covering 30 years. Then, three different estimation methods of EGS-E efficiency with comb-shaped and chessboard-shaped underground tunnels, respectively, are studied, and the research ideas for the estimation of the EGS-E system with more complicated cobweb-shaped tunnels are pointed out.


---

## Body

## 1. Introduction

Geothermal energy will become an important part of energy in the future because of its advantages of source stability, sustainability, and potential high utilization ratio. Geothermal energy stored at depths of less than 10 kilometers underground was estimated to be 170 million times the amount of heat released from all coals stored in the earth, based on a study by Pollack and Chapman in 1977. However, the scale of geothermal energy with temperature lower than 150°C at depths shallower than 3 kilometers underground is usually too small to maintain the demand for long-term stable electricity production. This portion of low-temperature hydrothermal resource only accounts for 10% of total geothermal energy stored in the earth [1]. And the official data released by the Ministry of Land and Resources of China shows that the total hot dry rock (HDR) geothermal resources at depths from 3 to 10 kilometers in the mainland are equivalent to 860 trillion tons of coal, which is 260000 times the current annual energy consumption amount in China [2]. Therefore, a broader scheme of the enhanced geothermal system (EGS) which aims at exploiting geothermal energy from HDR at the depths of 3 to 10 kilometers has gradually attracted people’s attention.The current popular heat extraction method from HDR has been developed as EGS based on drilling technology (EGS-D). A high-permeability fracture system (artificial thermal storage) is established by reservoir stimulation technology such as hydraulic fracturing through an injection well. Injected cold water (or other fluids) is heated by the fracture heat exchange structure, and then, the hot water or steam generated by the production well is pumped out to the ground for power generation. The injection well, production well, and underground thermal reservoir form a closed loop system of high-temperature thermal fluids [3]. The keys to the development of HDR include (1) localization of the resource target area. It should be based on the terrestrial heat flow value and combined with the characteristics of the geological structure. Seismic exploration technology, electrical method, and electromagnetic method as well as gravity and magnetic method could be used to carry out detailed exploration in the selected research area [4–6]. (2) Improvement of the reservoir: deep drilling technology under high temperature and pressure conditions is relatively mature. The depths of oil exploration drilling and comprehensive scientific research drilling have exceeded 7000 meters. At present, hydraulic fracturing [7], chemical stimulation [8], and thermal stimulation [9], as well as the combination of these technologies, are commonly used in creating a heat reservoir in EGS. And the fracture monitoring technologies used frequently now include microseismic monitoring, acoustic emission, downhole imaging, and various tracers [10, 11]. (3) Research on THMC (Thermal-Hydraulic-Mechanical-Chemical) multifield coupling of EGS: some scholars have successfully studied the coupling problem of each process. However, the construction of the thermal reservoir model still needs to be developed, which will continue to be the focus of EGS numerical simulation research [12–14]. And in recent years, scholars and relevant institutions have put forward new alternative systems, such as EGS using supercritical carbon dioxide as circulating fluid [15]and EGS for producing supercritical water vapor [16] as well as Radiator-EGS that consists of a family of vertically interconnected vanes produced through sequential horizontal drilling and frac-induced rubblization [17].There were several typical EGS-D development cases in the world such as the Fenton Hill, Hijiori, Ogachi, Cooper Basin, and Soultz. However, none of them has achieved large-scale commercial electricity generation owing to the disadvantages of small scale, low efficiency, geographical restriction, and so on [18]. The biggest challenge it faces is how to economically construct an artificial heat repository. In the HDR development report written by researchers in MIT, it was also pointed out that the key to cost-efficient development of EGS in the next 20 years is to obtain economical and effective multireservoir construction technology to ensure enough heat content (volume>1km3) for long-term geothermal development [19].Past experience in geothermal reservoir construction had shown that the development of an artificial reservoir was mainly caused by shear failure [20] of existing joints, which is different from tensile cracking caused by conventional hydraulic fracturing in oil and gas reservoirs. The shear failure part in the stimulation only occurs where fluid pressure becomes lower than the in situ minimum principal stress of the target reservoir. With the current technology, it is difficult to predict the distribution of the stress field around a well before being drilled, especially the distribution of stress field far from a wellbore. And directions of the in situ principal stresses may change with depth, which as well affects predicting the direction of fracture extension [21, 22]. It is usually difficult to predict the direction of fracture extension in the absence of accurate downhole data measurement. Even if related data had been available, fractures may not develop along the predicted direction. Therefore, it is difficult to achieve the interconnection between wells by hydraulic fracturing. Theoretically, it is better to first have the heat storage developed then find a way (e.g., microseismic monitoring technology and acoustic emission technology) to implement directional drilling secondly [23–25]. However, artificial fracture (caused by perforation as well as fracturing stimulation) only plays a dominant role surrounding the wellbore. The growth and expansion of heat storage during fracturing are mainly controlled by the existing natural fracture system (or joint distribution). The development of heat storage during hydraulic fracturing mainly comes from the activation of existing natural fractures which were controlled by the in situ stress field [26, 27]. Thus, EGS-D has a limitation in geographical selection. Even if the process of heat storage stimulation and directional drilling goes smooth, the cracks formed by hydraulic fracturing are often closed under the action of high in situ stress at depth. This makes the cracks disconnected, thus unable to form a sufficient volume of heat storage [28].In addition, there are many reasons causing the short-circuit effect of fluid flows, which also makes the failure of thermal recovery unavoidable. These include the following: (1) The use of proppant near the wellbore requires a higher injection pressure and flow rate, and excessive injection pressure will cause continuous fracture growth, resulting in increased water loss and/or fluid short circuit. (2) Repeated high-pressure stimulation of existing fractures may lead to more direct connection which results in fluid flow short circuit between injection and production wells. (3) If natural fractures are connected to wells, hydraulic fracturing may increase connectivity and lead to short circuit, especially when the well spacing is small. (4) Any operation of pressurizing reservoirs is irreversible and not always beneficial to developing heat storage. Long-term high-pressure water injection will cause irreversible damage to rocks, resulting in fluid short circuit and excessive water loss into farther regions. (5) Deeper and shallower reservoirs may be communicated through fracture growth in a water injection test or the wellbores penetrating two reservoirs, which will affect the construction of multiple reservoirs [29–33].In a word, from the development status of traditional EGS-D, it can be seen that innovative breakthroughs are urgently needed in the research of the deep geothermal exploitation scheme.


## 2. Materials and Methods

### 2.1. Enhanced Geothermal System Based on Excavation Technology (EGS-E)

Aiming at the disadvantages of EGS-D, such as difficulty on building large-scale stable heat storage, small water flow rate, and easy cause of contamination, a new EGS based on excavation technology (EGS-E) was recently proposed to provide a new scheme for the exploitation of deep HDR heat [34].As shown in Figure1, it includes (1) large heat flow formed through the excavation of a super large shaft, (2) heat source with large volume and high permeability due to crack and fragmentation formation through drilling and blasting which are implemented inside the tunnels, and (3) heat storage with a large capacity and high conductivity formed through the caving method based upon the excavation of an underground hot water reservoir and tunnels radially stretched out. In EGS-E, the horizontal tunnels with 360-degree distribution at different depths of the shaft can be excavated through mining excavation technology, and multiple stable large-scale heat storages can be constructed by the caving method from the tunnels. Compared with EGS-D, it can tremendously extend the scale of stable heat storage, expand the area of heat exchange, and upgrade the magnitude of geothermal transport to achieve large-scale geothermal development.Figure 1
Schematic illustration of the EGS-E concept.At present, the Mponeng Gold Mine in South Africa has mastered the excavation technology at a depth of 4350 meters [35]. As for the high temperature in the process of excavation, it is urgently needed to develop ice-cooled radiation cooling technology and local refrigeration cooling technology, new-type heat-resistant materials (e.g., thermal insulate lining sandwiched membrane), intelligent machinery, remote control robots, and so forth [36]. A R&D organization in South Africa has recently developed a concept robot for drilling holes. The Kiruna Iron Mine in northern Sweden which is the largest underground mine in the world has basically realized “unattended intelligent mining.” These technological breakthroughs are bringing “unmanned mining” closer and closer to us, making EGS-E more and more feasible. It is true that the large-scale production in EGS-E also needs to reduce costs in order to compete with other basic power generation technologies [37]. While attaining a large amount of hot water (or another suitable fluid or fluids for better heat storage, exchange, and transport, e.g., supercritical CO2), EGS-E can be combined with deep mining to reduce the engineering cost of EGS-E and provide an active cooling scheme for the high-temperature environment of deep mining so as to achieve a win-win situation of resources and energy development. It can also be combined with underground quarrying. Underground granite without weathering not only is ideal high-quality stone but also can be developed into building material to compensate for excavation costs. In the meantime, it can protect the ground surface environment. Moreover, the construction of an underground geothermal power plant can also be considered, which can reduce the loss of heat during transportation and save the fee for long distance transportation of fluid.However, it is difficult to realize the numerical simulation of a super large scale of the heat mining system in EGS-E without powerful supercomputing. Therefore, preliminary methods to concisely estimate the early performance of EGS-E concept models are proposed in this paper, and their application scope and error evaluation are also studied, which will provide a viable alternative to supercomputing for the subsequent quantitative research of the corresponding relationship between a typical deep engineering structure and its heat extraction efficiency. In the actual EGS-E model, there are many forced heat convection zones consisting of fracture flow and rock mass around tunnels. After heat transfer between fracture flow and HDR, all hot water is collected to the hot water pool at the bottom of the shaft through the forced convection scheme with circulating pipelines installed on the tunnel walls to realize the optimal control of heat transfer structures. Because the main purpose of this paper is to study the effect of estimation for heat extraction efficiency of EGS-E with different tunnel layouts and the transient simulation to large-scale heat convection of fracture flow in EGS-E requires supercomputing, the mechanism of heat transfer is simplified into the large-scale transient heat conduction to the tunnel wall with an equivalent homogenized thermal conductivity at this stage. The equivalent thermal conductivity here is considered to include the homogenized contribution of heat convection in fracture flow and heat conduction in rock. With this simplification, factors influencing heat extraction of tunnel walls are studied first.


### 2.2. Factors Influencing Heat Extraction of Tunnel Walls in EGS-E

Referring to previous studies [38, 39], the heat conducted in the radial direction of surrounding rock is much larger than that in the axial direction; thus, the latter can be ignored for analytical solution. An axially symmetrical plane problem that the transient heat conduction from high-temperature surrounding rock to a circular tunnel for 30 years is solved in this section to study the factors influencing heat extraction of tunnel walls in EGS-E.As shown in Figure2, considering that the cross sections of the tunnel and rock are both circular. The radii from the inside to the outside are R1 and R2, respectively. Based on the previous study [40] and according to the existing international criteria for the development and utilization of HDR [41], R1 is 3 m and R2 is 100 m; thus, it can be considered that the outer boundary is far enough from the tunnel and its temperature is not affected by heat conduction within 30 years, namely, the outer boundary can be regarded as thermal insulation. It is assumed that the temperature of tunnel wall T1 is 150°C constantly with continuous heat exchange between cold water and the tunnel wall. And the initial temperature of surrounding rock T2 is 250°C. The change of temperature field distribution in surrounding rock of the tunnel is studied. And the annual heat extraction efficiency of the tunnel wall for 30 years can be obtained through data processing with MATLAB.Figure 2
An axially symmetrical plane problem of the transient heat conduction from high-temperature surrounding rock to a circular tunnel.The corresponding analytical solution to the one-dimensional transient heat conduction problem shown in Figure2 is to solve the one-dimensional homogeneous heat conduction differential equation in cylindrical coordinates, which is expressed as follows [42]:
(1)∂2T∂r2+1r∂T∂r=1α∂T∂t,R1≤r≤R2,where α is the thermal diffusivity of surrounding rock and α=λ/cρ. λ, c, and ρ are the thermal conductivity, heat capacity, and density of surrounding rocks, respectively.The temperature of inner boundaryT1 is assumed to be 150°C constantly, which is the first type of boundary condition. The outer boundary is regarded as thermal insulation, namely, the heat flux is 0, which is the second type of boundary condition. Initial temperature of surrounding rock T2 is 250°C. Following analytical solution of the surrounding rock, temperature distribution is obtained through the separation variable method and the orthogonal expansion method [43]:
(2)Tr,t=T1+T2−T1∑m=1∞1Nβme−αβm2tR0βm,r∫R1R2r′R0βm,r′dr′,1Nβm=π22βm2J02R1βmJ02R1βm−J12R2βm,R0βm,r=−J0BmrY1R2βm+J1R2βmY0βmr,where J0 and J1 are the first kind of zero-order and first-order Bessel functions and Y0 and Y1 are the second kind of zero-order and first-order Bessel functions. All the eigenvalues βm above need to be accumulated, and the values of βm are the positive roots of the following equation [44]:
(3)−J0R1βmY1R2βm+J1R2βmY0R1βm=0.Then, the average temperatureT¯ of the model at t=1a,2a,⋯,30a is obtained after processing the temperature of surrounding rock at every point and every moment with MATLAB. And the annual temperature decline ΔT¯ of the model is obtained through subtraction of T¯ year by year within 30 years. According to the formula Q=c⋅m⋅ΔT¯, the analytical solution of the annual heat extraction rate of the one-meter-long tunnel wall varying with time can be obtained.At a certain initial temperature of surrounding rock, the heat extraction efficiency of the tunnel wall varies with the wall temperature, radius of the tunnel, and thermodynamic parameters of rock. Therefore, on the basis of the above formulas, the effects of three factors on the annual heat extraction rate of the tunnel wall are studied, respectively. The corresponding changeable calculation conditions are shown in Table1, and the results are shown in Figures 3–5.Table 1
Different wall temperatures, radii of the tunnel, and rock types.
T1 (°C)
70
90
110
130
150
R1 (m)
1.0
1.5
2.0
2.5
3.0
Rock type (α (m2/d))
Granite (0.1134)
Gneiss (0.0992)
Limestone (0.0827)
Basalt (0.0742)
Dry shale (0.0650)Figure 3
Analytical solution of the annual heat extraction rate varying with time for the one-meter-long tunnel under different wall temperatures when the rock type is granite, andT1 is 250°C, α is 0.1134 m2/d, and R1 is 3 m.Figure 4
Analytical solution of the annual heat extraction rate varying with time for the one-meter-long tunnel under different radii of the tunnel when the rock type is granite, andT1 is 250°C, α is 0.1134 m2/d, and T2 is 150°C.Figure 5
Analytical solution of the annual heat extraction rate varying with time for the one-meter-long tunnel under different rock types whenT1 is 250°C, R1 is 3 m, and T2 is 150°C.From Figures3 to 5, it can be seen directly that at a certain initial temperature of surrounding rock, the annual heat extraction rate of the tunnel wall increases with the lower constant temperature of the tunnel wall, the larger radius of the tunnel, or the higher thermal diffusivity of surrounding rock.


### 2.3. Estimation of EGS-E Performance with Different Tunnel Layouts

It is difficult to realize the numerical simulation of the super-large-scale heat mining system in EGS-E without powerful supercomputing. Does superposition of performance from each individual component give out an acceptable estimate on that of the structure made of them in the early heat mining production? Therefore, on the basis of previous studies, estimation methods of EGS-E efficiency with comb-shaped underground tunnels and chessboard-shaped underground tunnels, respectively, are proposed in this paper. And their application scope and error evaluation are also studied, which will provide the basis for subsequent research. According to the existing technical conditions, the international criteria for the development and utilization of HDR include that the volume of heat storage should be generally larger than 1 km3 so as to have development value. Thus, the size of surrounding rock in the numerical model of this section is 1000m∗1000m∗1000m, and the models with chessboard-shaped and comb-shaped underground tunnels together with their finite element meshes are shown in Figures 6 and 7, respectively. The ranges of chessboard-shaped tunnels and comb-shaped tunnels are 1000m∗1000m and 1000m∗200m, respectively. The cross section of the tunnel is still circular, and the radius of the tunnel is 3 m. The initial temperature of surrounding rock is 250°C. It is assumed that the temperature of the tunnel wall is 150°C constantly with continuous heat exchange between cold water and the tunnel wall. The surrounding boundary is set as thermal insulation, and the type of surrounding rock is granite. The equivalent thermal conductivity should be set as 175 W/(m·K) according to half of the available thermal energy of the cubic geothermal field to be mined in 30 years with chessboard-shaped tunnels (namely, the average temperature of the cubic geothermal field is reduced from 250°C to 200°C), which is reasonable and achievable under the condition of forced convection. Other thermodynamic parameters required for simulation are shown in Table 2, which are derived from the built-in material library of COMSOL Multiphysics.Figure 6
Numerical discretization of the simplified chessboard-shaped tunnel model of EGS-E and its local enlarged graph.Figure 7
Numerical discretization of the simplified comb-shaped tunnel model of EGS-E and its local enlarged graph.Table 2
Required thermodynamic parameters for simulation.
Rock type
Heat capacity (J/(kg·K))
Density (kg/m3)
Equivalent thermal conductivity (W/(m·K))
Granite
850
2600
175The core idea of the estimation is to turn the simulation of a large-scale structure into the superposition of performance from each individual component. As for the simplified chessboard-shaped tunnel model and comb-shaped tunnel model of EGS-E shown in Figures6 and 7, they can both be decomposed into several cross tunnels. Based on the following three models shown in Figures 8(a), 8(b), and 8(c), respectively, different superposition methods for estimation of EGS-E efficiency with chessboard-shaped and comb-shaped underground tunnels are studied.Figure 8
Three superposition models for estimation of EGS-E efficiency. (a) The ranges of tunnels and surrounding rock are200m∗200m and 1000m∗1000m∗1000m, respectively. (b) The ranges of tunnels and surrounding rock are 1000m∗1000m and 1000m∗1000m∗1000m, respectively. (c) The ranges of tunnels and surrounding rock are 200m∗200m and 200m∗200m∗1000m, respectively.
(a)
(b)
(c)As shown in Figure8(a), the range of surrounding rock in the model is 1000m∗1000m∗1000m. The range of cross tunnels is 200m∗200m, which is the same as the cross tunnels marked by the red dotted frame in Figures 6 and 7. Namely, the range of chessboard-shaped tunnels in Figure 6 is equivalent to 25 cross tunnels whose range is 200m∗200m, and the range of comb-shaped tunnels in Figure 7 is equivalent to 5 cross tunnels. These models have boundary effect in simulation calculation, but the tunnels are far from reaching the edge of the geothermal field in the actual situation. Therefore, considering the heat extraction efficiency of the whole chessboard-shaped or comb-shaped tunnels is as several times as that of their central cross tunnels to reduce the influence of the model boundary on estimation. The initial conditions, boundary conditions, and parameter settings required for calculation are the same as those of the models in Figures 6 and 7. The heat extraction rates of cross tunnel walls for those three models in Figure 8 at any time within 30 years are obtained and compared with those of chessboard-shaped and comb-shaped tunnels to study the application scope and error evaluation of three different models in Figure 8.Figure8(b) shows another estimation scheme. The range of surrounding rock in the model is 1000m∗1000m∗1000m. The range of cross tunnels is 1000m∗1000m. The initial conditions, boundary conditions, and parameter settings required for calculation are the same as those of the models in Figures 6 and 7, but only the heat extraction rate of cross tunnels marked by the red dotted frame in Figure 8(b) is taken for comparison. And the range of this part is also 200m∗200m.In the estimation scheme shown in Figure8(c), the range of surrounding rock is 200m∗200m∗1000m, and the range of cross tunnels is 200m∗200m. The initial conditions, boundary conditions, and parameter settings required for calculation are the same as those of the models in Figures 6 and 7. The heat extraction rate of the cross tunnel walls at any time within 30 years is obtained and compared with those of central cross tunnels marked in Figures 6 and 7.The heat extraction rates of cross tunnels with a range of200m∗200m at any time within 30 years based on the above three estimation schemes are compared with those of central cross tunnels with the same range marked in chessboard-shaped and comb-shaped models, respectively, as shown in Figure 9. It can be seen directly from the figure that the relative error of scheme (c) for chessboard-shaped tunnels is the smallest in superposition estimate of the heat extraction rate covering 30 years. This is because of the fact that for the chessboard-shaped tunnel model in Figure 6, its tunnel part and rock part can be exactly divided into 25 models as shown in Figure 8(c) with the same layout and size. However, there is still a tiny relative error that is less than 1% within 30 years. This is because the thermal field part distributed evenly to the central cross tunnels will begin to be affected by heat extraction of other surrounding tunnels at a certain point within 30 years with the equivalent thermal conductivity, and this part of heat loss cannot be completely compensated for by that itself “plunders” from other thermal field parts around it due to the boundary effect in simulation.Figure 9
Relative error in the heat extraction rate estimated by different superposition methods varying with time for EGS-E with different tunnel layouts.For comb-shaped tunnels, the relative error is the smallest when choosing scheme (b), which falls in less than 20% in 30 years. Analysis can also be given from the perspective of average division of the total thermal field volume to the total tunnel length. Scheme (b) is relatively the best for efficiency estimation of comb-shaped tunnels because its total length of tunnels and range of thermal field are the same with those of the model shown in Figure7. However, except for the central cross tunnels marked by the red dotted frame in Figure 8(b), the layout of other tunnels is quite different from that of comb-shaped tunnels, so the interaction of central cross tunnels with surrounding tunnels in heat extraction is quite different, which is the main source of its relative error of less than 20%.And for both chessboard-shaped and comb-shaped tunnels, the largest relative error of the above three estimation schemes during the first year is only about 20%. This conclusion is also of practical significance to estimate the early performance of EGS-E.Generally speaking, as explained above, it is advisable to replace it with the superposition estimation schemes studied in this paper to some extent when there is no powerful supercomputing to realize the numerical simulation of the super-large-scale heat mining system in EGS-E. And the analysis of relative error can also provide inspiration for further improvement of the estimation scheme.


## 2.1. Enhanced Geothermal System Based on Excavation Technology (EGS-E)

Aiming at the disadvantages of EGS-D, such as difficulty on building large-scale stable heat storage, small water flow rate, and easy cause of contamination, a new EGS based on excavation technology (EGS-E) was recently proposed to provide a new scheme for the exploitation of deep HDR heat [34].As shown in Figure1, it includes (1) large heat flow formed through the excavation of a super large shaft, (2) heat source with large volume and high permeability due to crack and fragmentation formation through drilling and blasting which are implemented inside the tunnels, and (3) heat storage with a large capacity and high conductivity formed through the caving method based upon the excavation of an underground hot water reservoir and tunnels radially stretched out. In EGS-E, the horizontal tunnels with 360-degree distribution at different depths of the shaft can be excavated through mining excavation technology, and multiple stable large-scale heat storages can be constructed by the caving method from the tunnels. Compared with EGS-D, it can tremendously extend the scale of stable heat storage, expand the area of heat exchange, and upgrade the magnitude of geothermal transport to achieve large-scale geothermal development.Figure 1
Schematic illustration of the EGS-E concept.At present, the Mponeng Gold Mine in South Africa has mastered the excavation technology at a depth of 4350 meters [35]. As for the high temperature in the process of excavation, it is urgently needed to develop ice-cooled radiation cooling technology and local refrigeration cooling technology, new-type heat-resistant materials (e.g., thermal insulate lining sandwiched membrane), intelligent machinery, remote control robots, and so forth [36]. A R&D organization in South Africa has recently developed a concept robot for drilling holes. The Kiruna Iron Mine in northern Sweden which is the largest underground mine in the world has basically realized “unattended intelligent mining.” These technological breakthroughs are bringing “unmanned mining” closer and closer to us, making EGS-E more and more feasible. It is true that the large-scale production in EGS-E also needs to reduce costs in order to compete with other basic power generation technologies [37]. While attaining a large amount of hot water (or another suitable fluid or fluids for better heat storage, exchange, and transport, e.g., supercritical CO2), EGS-E can be combined with deep mining to reduce the engineering cost of EGS-E and provide an active cooling scheme for the high-temperature environment of deep mining so as to achieve a win-win situation of resources and energy development. It can also be combined with underground quarrying. Underground granite without weathering not only is ideal high-quality stone but also can be developed into building material to compensate for excavation costs. In the meantime, it can protect the ground surface environment. Moreover, the construction of an underground geothermal power plant can also be considered, which can reduce the loss of heat during transportation and save the fee for long distance transportation of fluid.However, it is difficult to realize the numerical simulation of a super large scale of the heat mining system in EGS-E without powerful supercomputing. Therefore, preliminary methods to concisely estimate the early performance of EGS-E concept models are proposed in this paper, and their application scope and error evaluation are also studied, which will provide a viable alternative to supercomputing for the subsequent quantitative research of the corresponding relationship between a typical deep engineering structure and its heat extraction efficiency. In the actual EGS-E model, there are many forced heat convection zones consisting of fracture flow and rock mass around tunnels. After heat transfer between fracture flow and HDR, all hot water is collected to the hot water pool at the bottom of the shaft through the forced convection scheme with circulating pipelines installed on the tunnel walls to realize the optimal control of heat transfer structures. Because the main purpose of this paper is to study the effect of estimation for heat extraction efficiency of EGS-E with different tunnel layouts and the transient simulation to large-scale heat convection of fracture flow in EGS-E requires supercomputing, the mechanism of heat transfer is simplified into the large-scale transient heat conduction to the tunnel wall with an equivalent homogenized thermal conductivity at this stage. The equivalent thermal conductivity here is considered to include the homogenized contribution of heat convection in fracture flow and heat conduction in rock. With this simplification, factors influencing heat extraction of tunnel walls are studied first.


## 2.2. Factors Influencing Heat Extraction of Tunnel Walls in EGS-E

Referring to previous studies [38, 39], the heat conducted in the radial direction of surrounding rock is much larger than that in the axial direction; thus, the latter can be ignored for analytical solution. An axially symmetrical plane problem that the transient heat conduction from high-temperature surrounding rock to a circular tunnel for 30 years is solved in this section to study the factors influencing heat extraction of tunnel walls in EGS-E.As shown in Figure2, considering that the cross sections of the tunnel and rock are both circular. The radii from the inside to the outside are R1 and R2, respectively. Based on the previous study [40] and according to the existing international criteria for the development and utilization of HDR [41], R1 is 3 m and R2 is 100 m; thus, it can be considered that the outer boundary is far enough from the tunnel and its temperature is not affected by heat conduction within 30 years, namely, the outer boundary can be regarded as thermal insulation. It is assumed that the temperature of tunnel wall T1 is 150°C constantly with continuous heat exchange between cold water and the tunnel wall. And the initial temperature of surrounding rock T2 is 250°C. The change of temperature field distribution in surrounding rock of the tunnel is studied. And the annual heat extraction efficiency of the tunnel wall for 30 years can be obtained through data processing with MATLAB.Figure 2
An axially symmetrical plane problem of the transient heat conduction from high-temperature surrounding rock to a circular tunnel.The corresponding analytical solution to the one-dimensional transient heat conduction problem shown in Figure2 is to solve the one-dimensional homogeneous heat conduction differential equation in cylindrical coordinates, which is expressed as follows [42]:
(1)∂2T∂r2+1r∂T∂r=1α∂T∂t,R1≤r≤R2,where α is the thermal diffusivity of surrounding rock and α=λ/cρ. λ, c, and ρ are the thermal conductivity, heat capacity, and density of surrounding rocks, respectively.The temperature of inner boundaryT1 is assumed to be 150°C constantly, which is the first type of boundary condition. The outer boundary is regarded as thermal insulation, namely, the heat flux is 0, which is the second type of boundary condition. Initial temperature of surrounding rock T2 is 250°C. Following analytical solution of the surrounding rock, temperature distribution is obtained through the separation variable method and the orthogonal expansion method [43]:
(2)Tr,t=T1+T2−T1∑m=1∞1Nβme−αβm2tR0βm,r∫R1R2r′R0βm,r′dr′,1Nβm=π22βm2J02R1βmJ02R1βm−J12R2βm,R0βm,r=−J0BmrY1R2βm+J1R2βmY0βmr,where J0 and J1 are the first kind of zero-order and first-order Bessel functions and Y0 and Y1 are the second kind of zero-order and first-order Bessel functions. All the eigenvalues βm above need to be accumulated, and the values of βm are the positive roots of the following equation [44]:
(3)−J0R1βmY1R2βm+J1R2βmY0R1βm=0.Then, the average temperatureT¯ of the model at t=1a,2a,⋯,30a is obtained after processing the temperature of surrounding rock at every point and every moment with MATLAB. And the annual temperature decline ΔT¯ of the model is obtained through subtraction of T¯ year by year within 30 years. According to the formula Q=c⋅m⋅ΔT¯, the analytical solution of the annual heat extraction rate of the one-meter-long tunnel wall varying with time can be obtained.At a certain initial temperature of surrounding rock, the heat extraction efficiency of the tunnel wall varies with the wall temperature, radius of the tunnel, and thermodynamic parameters of rock. Therefore, on the basis of the above formulas, the effects of three factors on the annual heat extraction rate of the tunnel wall are studied, respectively. The corresponding changeable calculation conditions are shown in Table1, and the results are shown in Figures 3–5.Table 1
Different wall temperatures, radii of the tunnel, and rock types.
T1 (°C)
70
90
110
130
150
R1 (m)
1.0
1.5
2.0
2.5
3.0
Rock type (α (m2/d))
Granite (0.1134)
Gneiss (0.0992)
Limestone (0.0827)
Basalt (0.0742)
Dry shale (0.0650)Figure 3
Analytical solution of the annual heat extraction rate varying with time for the one-meter-long tunnel under different wall temperatures when the rock type is granite, andT1 is 250°C, α is 0.1134 m2/d, and R1 is 3 m.Figure 4
Analytical solution of the annual heat extraction rate varying with time for the one-meter-long tunnel under different radii of the tunnel when the rock type is granite, andT1 is 250°C, α is 0.1134 m2/d, and T2 is 150°C.Figure 5
Analytical solution of the annual heat extraction rate varying with time for the one-meter-long tunnel under different rock types whenT1 is 250°C, R1 is 3 m, and T2 is 150°C.From Figures3 to 5, it can be seen directly that at a certain initial temperature of surrounding rock, the annual heat extraction rate of the tunnel wall increases with the lower constant temperature of the tunnel wall, the larger radius of the tunnel, or the higher thermal diffusivity of surrounding rock.


## 2.3. Estimation of EGS-E Performance with Different Tunnel Layouts

It is difficult to realize the numerical simulation of the super-large-scale heat mining system in EGS-E without powerful supercomputing. Does superposition of performance from each individual component give out an acceptable estimate on that of the structure made of them in the early heat mining production? Therefore, on the basis of previous studies, estimation methods of EGS-E efficiency with comb-shaped underground tunnels and chessboard-shaped underground tunnels, respectively, are proposed in this paper. And their application scope and error evaluation are also studied, which will provide the basis for subsequent research. According to the existing technical conditions, the international criteria for the development and utilization of HDR include that the volume of heat storage should be generally larger than 1 km3 so as to have development value. Thus, the size of surrounding rock in the numerical model of this section is 1000m∗1000m∗1000m, and the models with chessboard-shaped and comb-shaped underground tunnels together with their finite element meshes are shown in Figures 6 and 7, respectively. The ranges of chessboard-shaped tunnels and comb-shaped tunnels are 1000m∗1000m and 1000m∗200m, respectively. The cross section of the tunnel is still circular, and the radius of the tunnel is 3 m. The initial temperature of surrounding rock is 250°C. It is assumed that the temperature of the tunnel wall is 150°C constantly with continuous heat exchange between cold water and the tunnel wall. The surrounding boundary is set as thermal insulation, and the type of surrounding rock is granite. The equivalent thermal conductivity should be set as 175 W/(m·K) according to half of the available thermal energy of the cubic geothermal field to be mined in 30 years with chessboard-shaped tunnels (namely, the average temperature of the cubic geothermal field is reduced from 250°C to 200°C), which is reasonable and achievable under the condition of forced convection. Other thermodynamic parameters required for simulation are shown in Table 2, which are derived from the built-in material library of COMSOL Multiphysics.Figure 6
Numerical discretization of the simplified chessboard-shaped tunnel model of EGS-E and its local enlarged graph.Figure 7
Numerical discretization of the simplified comb-shaped tunnel model of EGS-E and its local enlarged graph.Table 2
Required thermodynamic parameters for simulation.
Rock type
Heat capacity (J/(kg·K))
Density (kg/m3)
Equivalent thermal conductivity (W/(m·K))
Granite
850
2600
175The core idea of the estimation is to turn the simulation of a large-scale structure into the superposition of performance from each individual component. As for the simplified chessboard-shaped tunnel model and comb-shaped tunnel model of EGS-E shown in Figures6 and 7, they can both be decomposed into several cross tunnels. Based on the following three models shown in Figures 8(a), 8(b), and 8(c), respectively, different superposition methods for estimation of EGS-E efficiency with chessboard-shaped and comb-shaped underground tunnels are studied.Figure 8
Three superposition models for estimation of EGS-E efficiency. (a) The ranges of tunnels and surrounding rock are200m∗200m and 1000m∗1000m∗1000m, respectively. (b) The ranges of tunnels and surrounding rock are 1000m∗1000m and 1000m∗1000m∗1000m, respectively. (c) The ranges of tunnels and surrounding rock are 200m∗200m and 200m∗200m∗1000m, respectively.
(a)
(b)
(c)As shown in Figure8(a), the range of surrounding rock in the model is 1000m∗1000m∗1000m. The range of cross tunnels is 200m∗200m, which is the same as the cross tunnels marked by the red dotted frame in Figures 6 and 7. Namely, the range of chessboard-shaped tunnels in Figure 6 is equivalent to 25 cross tunnels whose range is 200m∗200m, and the range of comb-shaped tunnels in Figure 7 is equivalent to 5 cross tunnels. These models have boundary effect in simulation calculation, but the tunnels are far from reaching the edge of the geothermal field in the actual situation. Therefore, considering the heat extraction efficiency of the whole chessboard-shaped or comb-shaped tunnels is as several times as that of their central cross tunnels to reduce the influence of the model boundary on estimation. The initial conditions, boundary conditions, and parameter settings required for calculation are the same as those of the models in Figures 6 and 7. The heat extraction rates of cross tunnel walls for those three models in Figure 8 at any time within 30 years are obtained and compared with those of chessboard-shaped and comb-shaped tunnels to study the application scope and error evaluation of three different models in Figure 8.Figure8(b) shows another estimation scheme. The range of surrounding rock in the model is 1000m∗1000m∗1000m. The range of cross tunnels is 1000m∗1000m. The initial conditions, boundary conditions, and parameter settings required for calculation are the same as those of the models in Figures 6 and 7, but only the heat extraction rate of cross tunnels marked by the red dotted frame in Figure 8(b) is taken for comparison. And the range of this part is also 200m∗200m.In the estimation scheme shown in Figure8(c), the range of surrounding rock is 200m∗200m∗1000m, and the range of cross tunnels is 200m∗200m. The initial conditions, boundary conditions, and parameter settings required for calculation are the same as those of the models in Figures 6 and 7. The heat extraction rate of the cross tunnel walls at any time within 30 years is obtained and compared with those of central cross tunnels marked in Figures 6 and 7.The heat extraction rates of cross tunnels with a range of200m∗200m at any time within 30 years based on the above three estimation schemes are compared with those of central cross tunnels with the same range marked in chessboard-shaped and comb-shaped models, respectively, as shown in Figure 9. It can be seen directly from the figure that the relative error of scheme (c) for chessboard-shaped tunnels is the smallest in superposition estimate of the heat extraction rate covering 30 years. This is because of the fact that for the chessboard-shaped tunnel model in Figure 6, its tunnel part and rock part can be exactly divided into 25 models as shown in Figure 8(c) with the same layout and size. However, there is still a tiny relative error that is less than 1% within 30 years. This is because the thermal field part distributed evenly to the central cross tunnels will begin to be affected by heat extraction of other surrounding tunnels at a certain point within 30 years with the equivalent thermal conductivity, and this part of heat loss cannot be completely compensated for by that itself “plunders” from other thermal field parts around it due to the boundary effect in simulation.Figure 9
Relative error in the heat extraction rate estimated by different superposition methods varying with time for EGS-E with different tunnel layouts.For comb-shaped tunnels, the relative error is the smallest when choosing scheme (b), which falls in less than 20% in 30 years. Analysis can also be given from the perspective of average division of the total thermal field volume to the total tunnel length. Scheme (b) is relatively the best for efficiency estimation of comb-shaped tunnels because its total length of tunnels and range of thermal field are the same with those of the model shown in Figure7. However, except for the central cross tunnels marked by the red dotted frame in Figure 8(b), the layout of other tunnels is quite different from that of comb-shaped tunnels, so the interaction of central cross tunnels with surrounding tunnels in heat extraction is quite different, which is the main source of its relative error of less than 20%.And for both chessboard-shaped and comb-shaped tunnels, the largest relative error of the above three estimation schemes during the first year is only about 20%. This conclusion is also of practical significance to estimate the early performance of EGS-E.Generally speaking, as explained above, it is advisable to replace it with the superposition estimation schemes studied in this paper to some extent when there is no powerful supercomputing to realize the numerical simulation of the super-large-scale heat mining system in EGS-E. And the analysis of relative error can also provide inspiration for further improvement of the estimation scheme.


## 3. Discussion

(1)
The paper [40] has inspired another type of underground tunnel layout (cobweb), which is also convenient for the outward stretching of the underground tunnel structure, thus further increasing the volume of heat storage, as shown in Figure 10. Similarly, the heat extraction efficiency of cobweb-shaped tunnels can also be estimated. The difference lies in the need for numerical simulation of heat transfer of several central cross tunnels with different sizes (as marked by the red dotted frame in Figure 10) and several single cylindrical tunnels with different lengths (namely, the intervals between different cross tunnels). Then, the estimation results of overall heat extraction efficiency of the model with cobweb-shaped tunnels can be accumulated
(2)
The main purpose of this paper is to study the effect of estimation for heat extraction efficiency of EGS-E with different tunnel layouts, and the transient simulation to large-scale heat convection of fracture flow in EGS-E requires supercomputing; thus, the mechanism of heat transfer is simplified into the large-scale transient heat conduction to the tunnel wall with an equivalent homogenized thermal conductivity for the time beingFigure 10
Simplified cobweb-shaped tunnel model of EGS-E.The accuracy of the 3-D model simulation is verified by the following two aspects: (i) The results of simulations for models with different discretizations are nearly the same. (ii) Relative error in the heat extraction rate estimated by the superposition method shown in Figure8(c) for EGS-E with chessboard-shaped tunnels is less than 1% within 30 years.This simplification should also consider and involve necessary limitations. That is, the modelling is mainly based on great simplification of a comprehensive homogenized thermal conductivity which already contains effects of the heat convection mechanism through fracture networks within surrounding rock. This simplification would surely depend on the sufficiency of enhanced fracturing, either by alike caving method based on blasting and natural collapse or by hydraulic fracturing stimulation through drilling boreholes. And the sufficiency will be dependent on the total contact area of fracture walls which are interfaces between rock and flows, total length and spacing of fracture and rock lumpiness, etc.In addition, the influence factors for heat convection in each individual fracture include not only the rock properties and contact area but also the fluid properties (thermal conductivity, viscosity, heat capacity and density, etc.) and the heat transfer surface properties (shape, size and roughness, etc.). It is also closely related to the phase transition. Bidirectional coupling between temperature field and flow field is an important factor that must be considered as well.In general, compared with heat conduction, the heat convection in fracture flow is a much more complicated process that affected on many factors, and it is necessary to conduct related research according to the classification of influence factors. Subsequently, as a complementary part, a representative local detailed numerical model should be studied based on the heat transfer between fracture flow and HDR, namely, selecting a small fractured area for numerical simulation of bidirectional coupling between the temperature field and flow field to comprehensively study the factors such as flow rate, flow flux, and heat extraction efficiency assuming that the heat storage is designed to consist of several relatively small fractured areas in parallel.


## 4. Conclusions

Geothermal energy will become an important energy component in the future because of its advantages of stability, sustainability, and efficient utilization. In particular, the development and utilization of deep geothermal energy from HDR have gradually attracted people’s attention. Aiming at mitigating the bottleneck of EGS-D, a new EGS-E based on excavation technology was proposed. In this paper, a simple and direct method for estimating the early performance of the large-scale deep geothermal heat mining is studied and established for its applicability in the subversive and innovative scheme, i.e., EGS-E large-scale heat mining, in the near future. A preliminary exploration is made to quantitatively study the corresponding relationship between the deep engineering structure and its heat extraction efficiency. The relevant researches and conclusions are as follows:(1)
The major characteristics (large heat flow, heat source with large volume and high permeability, and heat storage with large capacity and high conductivity) and advantages of EGS-E are introduced. The breakthroughs and prospects about key technologies involved in the construction and operation of the EGS-E system are expounded as well. Moreover, the innovative schemes to reduce the costs of EGS-E are also put forward
(2)
The effects of the tunnel wall temperature, tunnel radius, and rock type on the annual heat extraction rate are studied based on the analytical solution of a one-dimensional radial plane problem of the transient heat conduction through high-temperature surrounding rock to the tunnel wall covering 30 years. The results show that at a certain initial field temperature, the annual heat extraction rate from the tunnel wall increases with the lower inner boundary fixed temperature, the longer radius of the tunnel, or the greater thermal diffusivity of surrounding rock
(3)
Through undertaking numerical simulations with COMSOL Multiphysics, three different estimation methods of EGS-E efficiency with comb-shaped and chessboard-shaped underground tunnels, respectively, are proposed, and the research ideas for the estimation of the EGS-E system with more complicated cobweb-shaped tunnels are pointed out. Relative optimum estimation schemes for comb-shaped and chessboard-shaped underground tunnels are obtained, respectively. The relative error of scheme (c) is the smallest for the superposition estimation of the heat extraction rate in 30 years for chessboard-shaped tunnels, which has been less than 1% within 30 years. For comb-shaped tunnels, the relative error is the smallest when choosing scheme (b), which has been less than 20% in 30 years. And for chessboard-shaped and comb-shaped tunnels, the largest relative error of these three estimation schemes in the first year is only about 20%. Thus, generally speaking, it is advisable to take advantage of convenience and effectiveness of the superposition estimation schemes studied in this paper to some extent when there is no powerful supercomputing to realize the numerical simulation of the super-large-scale heat mining system in EGS-E



---

*Source: 3456307-2019-06-13.xml*