Numerical Study on Damage Zones Induced by Excavation and Ventilation in a High-Temperature Tunnel at Depth

Geothermal power is being regarded as depending on techniques derived from hydrocarbon production in worldwide current strategy. However, it has artificially been developed far less than its natural potentials due to technical restrictions. This paper introduces the Enhanced Geothermal System based on Excavation (EGS-E), which is an innovative scheme of geothermal energy extraction. Then, based on cohesion-weakening-friction-strengthening model (CWFS) and literature investigation of granite test at high temperature, the initiation, propagation of excavation damaged zones (EDZs) under unloading and the EDZs scale in EGS-E closed to hydrostatic pressure state is studied. Finally, we have a discussion about the further evolution of surrounding rock stress and EDZs during ventilation is studied by thermal-mechanical coupling. The results show that the influence of high temperature damage on the mechanical parameters of granite should be considered; Lateral pressure coefficient affects the fracture morphology and scale of tunnel surrounding rock, and EDZs area is larger when the lateral pressure coefficient is 1.0 or 1.2; Ventilation of high temperature and high in-situ stress tunnel have a significant effect on the EDZs scale; Additional tensile stress is generated in the shallow of tunnel surrounding rock, and the compressive stress concentration transfers to the deep. EDZs experiences three expansion stages of slow, rapid and deceleration with cooling time, and the thermal insulation layer prolongs the slow growth stage.


Introduction
Geothermal resources have become a key research and development objective among the world's renewable clean energy options because of its source stability, potential high utilization ratio, tremendous reserve and widespread distribution [1]. According to theoretical calculations, the energy reserves in the upper 10 km of the earth's crust are approximately 1.3 × 10 27 J. These geothermal reserves could supply global energy use for approximately 217 million years [2]. Therefore, once large-scale utilization was realized, geothermal energy would act as a key technique to cope with the environmental problems brought from mainly consuming fossil fuels. This promising option can ease the energy crisis, promote and eventually realize sustainability for the prosperity of future world [3,4].
Currently, borehole-based enhanced geothermal system (EGS-D), in which the techniques of Drilling and hydraulic fracturing rooting out of petroleum industry operations are adopted, is the most widely used scheme on extracting heat from hot dry rocks [5]. However, to be frank, EGS-D has to date been failing to be successfully commercialized.
One of the root causes should lie at the drawbacks out of hydraulic fracturing, which is hardly controllable at precise and yields insufficient effectively stimulated volume within quartzite rocks. These two main defects have so far caused many obstacles in extracting heat from Hot Dry Rocks, including limited amount of heat extraction [6,7], difficult to optimal design of hydraulic fracturing [8,9], water loss, ground-water replacement, mixing and contamination [10,11], inducing earthquakes [12,13].
In recent years, deep underground associated excavation technologies are readily available [14][15][16]. Gold mining in South Africa has accessed into hard orebodies being buried over 4 km deep, and intelligent mine robots are rapidly iterating and evolving. Based on these backgrounds, aiming at solution to the bottleneck of EGS-D, a disruptively innovative geothermal extraction solution based on excavation rather than drilling, i.e., Enhanced Geothermal Systems based on Excavation technology (EGS-E), is proposed, as shown in Figure 1. EGS-E overcomes the shortcomings of the boreholes and fracturing based conventional EGS and provides a new large-scale scheme for developing hot dry rocks. The concept of EGS-E has the following connotations and advantages: (1) Excavation of shaft and roadways accessing the hot rock levels and then creating heat reservoirs artificially in a distributive way. The diameter of its vertical shafts is far much larger than that of boreholes, and its provided multi-level flow channels replaces the fractured ones in a far more controllable way, with far larger flow rates and flow amounts; (2) The horizontal roadways are excavated with certain spacing to form a room-and-pillar mining-like complex. Heat is extracted from surrounding rock being stimulated by long drilling, inter-hole fracturing, borehole-blasting or even caving method in which the utmost automatic interplay between stress concentration and gravitational collapse will be induced. Large volume and highly conductive heat reservoir is eventually constructed from roadway networks converging into focal caverns. Only in this way can predetermined volume of cracked zones for enhanced permeability and sufficient area for effective heat exchange be accurately controlled; (3) Heat collected from distributive but systematical roadway branches are transported over to the focused "hot pool" at each level. There focused heat exchange occurs between the colder working-fluid through inlet pipes and the much hotter water under a high pressure. Steam carrying heat was transferred to the relatively shallow level where installed turbines runs to generate electricity [7]. Different patterns of selection are adjusted from roadway complex and then some run and others rest in alternating cycles to ensure the entire system becomes stable and durable. Facing future in the world, engineering structure should be taken as a solution to mining heat from hot rocks which belongs to a vast but scattered resource. In this paper, rock properties under high temperature were first reviewed after investigating a large number of others' existing studies. Then, using numerical modelling, a section of typical roadway with the circular shape at cross section, which is excavated within high-temperature and high-stress hard rock, was analyzed. During the simulation, a sim-  The problem of deep rock excavation engineering, especially the EGS-E engineering, has two main characteristics simultaneously: high in-situ stress and high temperature conditions. For the former, many researchers studied the stress redistribution and damage extent of surrounding rock due to excavation by relying on empirical approach [17], numerical simulation [18,19] and theoretical calculation [20]. For the latter, relevant studies focused on the temperature evolution of surrounding rock [21], effective ventilation distance [22] and heat insulation [23] in tunnel ventilation. However, due to the temperature of the non-geothermal mining tunnel is not high (usually less than 80 • C), previous researches mostly studied the characteristic factors of deep engineering separately and independently. Only a limited number of studies have involved the effect of rock strength and deformation parameters changes caused by high temperature damaged on EDZs. On the other hand, the coupling effect between temperature and stress in ventilation process of damaged surrounding rock after excavation is lack of research, and the damage further evolution of surrounding rock induced by ventilation was not explored.
In this paper, rock properties under high temperature were first reviewed after investigating a large number of others' existing studies. Then, using numerical modelling, a section of typical roadway with the circular shape at cross section, which is excavated within high-temperature and high-stress hard rock, was analyzed. During the simulation, a simple, effective model CWFS was applied to preliminarily estimate the scale and zoning of Excavation Damage Zones in a quantitative way. Finally, the further evolution law of damage caused by tunnel ventilation after surrounding rock damage was studied.

Influence of High Temperature on Granite Properties
Hot Dry Rock geothermal energy resources are reserved in high-temperature depths which usually are radioactive igneous rocks such as granite. In this section, the microstructure of this type of rocks and its influence on its macroscopic properties under the corresponding temperatures which deviate from usual ground or shallow to rather deeply buried conditions are mainly investigated. Adopting conventional methods of thermal treatment, many researchers [24][25][26][27][28][29][30][31][32][33] have undertaken experimental or numerical simulation studies on rock physical and mechanical properties under high temperature.
It has been found that (1) Two microscopic mechanisms were essentially revealed by laboratory examinations when samples had been heated to various temperature levels. They are (a) the extension of existing cracks or the formation of new microcracks which is attributed to the thermal expansion of minerals: The mismatch of thermal expansions of different mineral grains causes microcracks to initiate at mineral contacts, and the anisotropy in the thermal expansion of individual mineral under restricted surrounding leads to the development of intra-granular microcracks; (b) Collapse of mineral structures due to thermally induced stress enhancement: The underlying mechanism is the heterogeneity among thermal deformation of mineral grains and the α-β phase transition of quartz [25]. For instance, granite taken from Quanzhou city, Fujian province of China was tested in the experimental study. As shown in Figure 2, the experimental results reveal that almost all destructive chemical reactions occurred at temperatures in excess of 300-400 Celsius degrees. In the range of from 100 • C to 573 • C, thermal cracks are mainly intergranular ones. While once temperature exceeded 573 • C, due to the phase transition of quartz, intra-granular cracks began to develop [26]. to 573 °C, thermal cracks are mainly intergranular ones. While once temperature exceeded 573 °C, due to the phase transition of quartz, intra-granular cracks began to develop [26]. (2) The transformation of mechanical parameters with temperature is due to the competition between strengthening and degradation effects, which were simultaneously induced by raising temperature [27]. The former is toughness enhancement of biotite under high temperature [28] and inter-crystalline contact compaction due to volume expansion [29]. The latter is intergranular cracks caused by non-uniform thermal expansion of different minerals after the temperature threshold was exceeded [30]. As shown in Figure 3, when temperature was below 200 °C, damage factors only dropped slightly indicating that mechanical properties were first strengthened due to water loss and mineral expansion. Over 200 up to 400 °C, damage factors first dropped shortly and then start rising, mechanical properties immediately deteriorated along with the appearance of intercrystalline cracks [31]. These findings were further confirmed by SEM analysis. After all, it should be concluded that under 300 °C the rock microstructures is only subject to minor change [29].  (2) The transformation of mechanical parameters with temperature is due to the competition between strengthening and degradation effects, which were simultaneously induced by raising temperature [27]. The former is toughness enhancement of biotite under high temperature [28] and inter-crystalline contact compaction due to volume expansion [29]. The latter is intergranular cracks caused by non-uniform thermal expansion of different minerals after the temperature threshold was exceeded [30]. As shown in Figure 3, when temperature was below 200 • C, damage factors only dropped slightly indicating that mechanical properties were first strengthened due to water loss and mineral expansion. Over 200 up to 400 • C, damage factors first dropped shortly and then start rising, mechanical properties immediately deteriorated along with the appearance of intercrystalline cracks [31]. These findings were further confirmed by SEM analysis. After all, it should be concluded that under 300 • C the rock microstructures is only subject to minor change [29]. to 573 °C, thermal cracks are mainly intergranular ones. While once temperature exceeded 573 °C, due to the phase transition of quartz, intra-granular cracks began to develop [26]. (2) The transformation of mechanical parameters with temperature is due to the competition between strengthening and degradation effects, which were simultaneously induced by raising temperature [27]. The former is toughness enhancement of biotite under high temperature [28] and inter-crystalline contact compaction due to volume expansion [29]. The latter is intergranular cracks caused by non-uniform thermal expansion of different minerals after the temperature threshold was exceeded [30]. As shown in Figure 3, when temperature was below 200 °C, damage factors only dropped slightly indicating that mechanical properties were first strengthened due to water loss and mineral expansion. Over 200 up to 400 °C, damage factors first dropped shortly and then start rising, mechanical properties immediately deteriorated along with the appearance of intercrystalline cracks [31]. These findings were further confirmed by SEM analysis. After all, it should be concluded that under 300 °C the rock microstructures is only subject to minor change [29].  (3) The brittleness-ductility transition of granite occurred at temperature as high as 600-700 • C [31]. At 25-600 • C, uniaxial compressive stress in granite dominated the tendency of brittle fracture process, which mainly induced tensile cracks along its axial direction. While at 700-800 • C, thermal stress in granite determined ductile damage process, which turned to result in post-peak ductile deformation [33].
(4) A series of triaxial compression tests were undertaken with granite samples extracted from Hunan Province of China under four different confining pressures, i.e., 10, 30, 60, 90 MPa, and four different temperatures, i.e., Room Temperature, 100, 200, 300 • C. The results [21] showed that within range of up to moderate temperature it would be stress but not temperature that dominantly influenced rock mechanical behavior [27,33]. Although the samples deformed more along with the levels of confining stress and temperature were raised, their stress-strain curves still showed brittle failure feature when the experimental conditions fell into the ranges of confining stress and temperature.
In most geothermal engineering, the rock temperature rarely exceeds 300 • C. Combined with the literature investigation, it ought be concluded that the hard rocks at a temperature that are not exceeding 300 • C still typically belong to brittle ones. Yet their physical characteristics, mechanical behaviors as well as failure process and patterns are roughly similar to those under room temperature.

Numerical Model to Study Tunneling Induced Excavation Damage Zone
Working roadway excavation in EGS-E will confront in-situ rocks with their highstress and high-brittleness problems. The dominant failure mechanisms of excavations in hard, massive rock masses at great depths are governed by the tensile strength of the rock mass and shown as extensile fractures. So-called macroscopic spalling is generated under high-magnitude compressive stresses at low confinement along the excavation boundary [34].
The photograph shown in Figure 4 was taken to view the cross-section of a testing tunnel with radius 1.75 m at depth 420 m within Lac du Bonnet granite. It was completed during Underground Research Laboratory 'Mine-by Experiment project' which was undertaken by Atomic Energy of Canada Limited at southeastern Manitoba. As well, photographs shown in Figure 5 are vividly telling the visible brittle failures surrounding designed open-space during deep tunnelling. It should be pointed out that the micro-crack like shearing damages which were difficult to detect with naked eyes would also concern hydraulic conductivity, especially in leakage and sealing around a tunnel. Its fundamental study is indispensable to any further implementation of permeability enhancement in order to construct an artificial geothermal reservoir in hot dry rocks.
(3) The brittleness-ductility transition of granite occurred at temperature as high as 600-700 °C [31]. At 25-600 °C, uniaxial compressive stress in granite dominated the tendency of brittle fracture process, which mainly induced tensile cracks along its axial direction. While at 700-800 °C, thermal stress in granite determined ductile damage process, which turned to result in post-peak ductile deformation [33].
(4) A series of triaxial compression tests were undertaken with granite samples extracted from Hunan Province of China under four different confining pressures, i.e., 10, 30, 60, 90 MPa, and four different temperatures, i.e., Room Temperature, 100, 200, 300 °C. The results [21] showed that within range of up to moderate temperature it would be stress but not temperature that dominantly influenced rock mechanical behavior [27,33]. Although the samples deformed more along with the levels of confining stress and temperature were raised, their stress-strain curves still showed brittle failure feature when the experimental conditions fell into the ranges of confining stress and temperature.
In most geothermal engineering, the rock temperature rarely exceeds 300 °C. Combined with the literature investigation, it ought be concluded that the hard rocks at a temperature that are not exceeding 300 °C still typically belong to brittle ones. Yet their physical characteristics, mechanical behaviors as well as failure process and patterns are roughly similar to those under room temperature.

Numerical Model to Study Tunneling Induced Excavation Damage Zone
Working roadway excavation in EGS-E will confront in-situ rocks with their highstress and high-brittleness problems. The dominant failure mechanisms of excavations in hard, massive rock masses at great depths are governed by the tensile strength of the rock mass and shown as extensile fractures. So-called macroscopic spalling is generated under high-magnitude compressive stresses at low confinement along the excavation boundary [34].
The photograph shown in Figure 4 was taken to view the cross-section of a testing tunnel with radius 1.75 m at depth 420 m within Lac du Bonnet granite. It was completed during Underground Research Laboratory 'Mine-by Experiment project' which was undertaken by Atomic Energy of Canada Limited at southeastern Manitoba. As well, photographs shown in Figure 5 are vividly telling the visible brittle failures surrounding designed open-space during deep tunnelling. It should be pointed out that the micro-crack like shearing damages which were difficult to detect with naked eyes would also concern hydraulic conductivity, especially in leakage and sealing around a tunnel. Its fundamental study is indispensable to any further implementation of permeability enhancement in order to construct an artificial geothermal reservoir in hot dry rocks.

Damage Initiation-Spalling Limit (DISL) Criterion and CWFS Model
Traditional rock failure criteria, e.g., Hoek Brown strength envelope, was applicable to surface or shallow engineering where rock mass is naturally so jointed that its response to excavation will be controlled by structural issues as its prerequisite condition. At deep underground, however, igneous rock like granite wouldn't be jointed to that extent and its strength yet can be regarded as roughly isotropic. Then there had been unfortunately no methodology available till recent few years. A common criterion on shear yielding couldn't imitate rock brittle failure until it reflected the following two requisite mechanisms: tensile fracture under low constraint and shear rupture under high constraint. The in-situ strength envelop was established by so-called DISL criterion [38], which came from accepted empirical observations in excavating massive to moderately jointed hard rocks and laboratory cases of axial stress-strain measurements to study the thresholds indicating the inception of nonlinear long-term yielding envelop. As illustrated in Figure 6, this criterion is an S-shaped or three-segment straight composite polylines, not a Mohr-Coulomb like linear arc or Hoek-Brown like continuous curve [39]. The reason lies at that rock heterogeneity subjecting to compression leads to local tensile cracks so as to create Griffith type of extensional fractures. Under low constraint conditions, rock mass cannot withstand failing at about one-tenth of its UCS. Much prior to the peak, instead, local high tensile stresses could exceed its corresponding tensile strength, eventually resulting in thin rock slabs referred to as spalling [39]. In laboratory samples, crack propagation is inhibited due to geometry effects, whereas near the excavation wall, cracks can freely propagate. The composite strength envelope for in-situ excavation of brittle rocks as shown by the red solid line was based on conceptual model DISL, which is compared with the corresponding laboratory peak and long-term yield strength envelopes.
The S-type criterion combines the cohesion weak-friction enhancement (CWFS) strength model at low confinement and the shear yield model at high confinement to characterize the lower brittleness at high confining pressures. At room temperature, the fracture and Highly Damaged Zone scale caused by hard rock excavation under high in-situ stress have been verified in some actual field measurements [40][41][42][43][44]. Diederichs et al. adopted this failure criterion to model the stress states corresponding to excavating caused boundary damage as well as surrounding yield. The approach had been testified through matching the failure morphology of Lac du Bonnet granite under the Underground Research Laboratory project within Manitoba province of Canada. Rojat et al. further applied it to investigating spalling phenomena during operate the Swiss Alps tunnel through gneiss, granodiorite and granite formations [40]. So far, it has been applied to simulate fracture open-

Damage Initiation-Spalling Limit (DISL) Criterion and CWFS Model
Traditional rock failure criteria, e.g., Hoek Brown strength envelope, was applicable to surface or shallow engineering where rock mass is naturally so jointed that its response to excavation will be controlled by structural issues as its prerequisite condition. At deep underground, however, igneous rock like granite wouldn't be jointed to that extent and its strength yet can be regarded as roughly isotropic. Then there had been unfortunately no methodology available till recent few years. A common criterion on shear yielding couldn't imitate rock brittle failure until it reflected the following two requisite mechanisms: tensile fracture under low constraint and shear rupture under high constraint. The insitu strength envelop was established by so-called DISL criterion [38], which came from accepted empirical observations in excavating massive to moderately jointed hard rocks and laboratory cases of axial stress-strain measurements to study the thresholds indicating the inception of nonlinear long-term yielding envelop. As illustrated in Figure 6, this criterion is an S-shaped or three-segment straight composite polylines, not a Mohr-Coulomb like linear arc or Hoek-Brown like continuous curve [39]. The reason lies at that rock heterogeneity subjecting to compression leads to local tensile cracks so as to create Griffith type of extensional fractures. Under low constraint conditions, rock mass cannot withstand failing at about one-tenth of its UCS. Much prior to the peak, instead, local high tensile stresses could exceed its corresponding tensile strength, eventually resulting in thin rock slabs referred to as spalling [39]. In laboratory samples, crack propagation is inhibited due to geometry effects, whereas near the excavation wall, cracks can freely propagate. The composite strength envelope for in-situ excavation of brittle rocks as shown by the red solid line was based on conceptual model DISL, which is compared with the corresponding laboratory peak and long-term yield strength envelopes.
ings of several hard brittle rock masses [41,43]. Among them, the Creighton Mine in Sudbury, Canada was worthy of being mentioned as the application reached depth 2400 m and the major principle stress up to 114.1 MPa [44]. Quantitative estimation of risk in spalling damage [45] and sensitivity in predictions using the DISL method has been evaluated [46] based on a large number of field data and empirical formulae [47][48][49][50]. The CWFS model is a strain-dependent yield criterion which accounts for the nonsimultaneous mobilization of friction and cohesion. The changes in these two strength components are controlled by plastic shear strain-a variable which is used as a proxy for damage. The yield criterion of CWFS model can be generalized by Mohr-Coulomb: where c and are cohesion and friction angle, respectively; is the normal stress; and are the characteristic plastic strains that control the loss of cohesion and the enhancement of friction angle, respectively.
According to Figure 7, the CWFS strength model requires six parameters to define: , the initial friction; , the mobilized friction angle; , peak cohesion; , residual cohesion; and are plastic strain limit for cohesion weakening and friction strengthening. In addition, parameters related to tensile strength loss are also used: , peak tensile strength; , residual tensile strength; , plastic strain limit for tensile strength weakening. The programme to obtain accurate CWFS parameters has been explained through performing extensive uniaxial and triaxial compressive tests [51].
In situ strength 0 Figure 6. The S-shaped or three-segment polyline envelope as in-situ rock excavation strength [39].
The S-type criterion combines the cohesion weak-friction enhancement (CWFS) strength model at low confinement and the shear yield model at high confinement to characterize the lower brittleness at high confining pressures. At room temperature, the fracture and Highly Damaged Zone scale caused by hard rock excavation under high in-situ stress have been verified in some actual field measurements [40][41][42][43][44]. Diederichs et al. adopted this failure criterion to model the stress states corresponding to excavating caused boundary damage as well as surrounding yield. The approach had been testified through matching the failure morphology of Lac du Bonnet granite under the Underground Research Laboratory project within Manitoba province of Canada. Rojat et al. further applied it to investigating spalling phenomena during operate the Swiss Alps tunnel through gneiss, granodiorite and granite formations [40]. So far, it has been applied to simulate fracture openings of several hard brittle rock masses [41,43]. Among them, the Creighton Mine in Sudbury, Canada was worthy of being mentioned as the application reached depth 2400 m and the major principle stress up to 114.1 MPa [44]. Quantitative estimation of risk in spalling damage [45] and sensitivity in predictions using the DISL method has been evaluated [46] based on a large number of field data and empirical formulae [47][48][49][50].
The CWFS model is a strain-dependent yield criterion which accounts for the nonsimultaneous mobilization of friction and cohesion. The changes in these two strength components are controlled by plastic shear strain-a variable which is used as a proxy for damage. The yield criterion of CWFS model can be generalized by Mohr-Coulomb: where c and ϕ are cohesion and friction angle, respectively; σ n is the normal stress; ε p c and ε p ϕ are the characteristic plastic strains that control the loss of cohesion and the enhancement of friction angle, respectively.
According to Figure 7, the CWFS strength model requires six parameters to define: ϕ i , the initial friction; ϕ m , the mobilized friction angle; c p , peak cohesion; c r , residual cohesion; ε p c and ε p ϕ are plastic strain limit for cohesion weakening and friction strengthening. In addition, parameters related to tensile strength loss are also used: t p , peak tensile strength; t r , residual tensile strength; ε p t , plastic strain limit for tensile strength weakening. The programme to obtain accurate CWFS parameters has been explained through performing extensive uniaxial and triaxial compressive tests [51].

Numerical Modelling Simulator
The finite-difference simulator FLAC 3D was used to model the mechanical process. This commercial software was chosen because it is convenient to use and widely applicable across mining and Geo-engineering fields. The incremental plastic parameter is defined by FLAC 3D as: where ∆ , ∆ , ∆ represent the major, mean, and minor principal plastic strain increment, respectively [52]. The focus of this modelling study is on those geologic conditions where the behavior of the engineering of interest is governed by continuum failure and damage of intact rock material rather than by deformation in distinct fragmented features. In other words, assuming deep rockmass of EGS-E is initially intact or massively jointed, i.e., its GSI index belongs to more than 75.

Model Geometry, Gridding, Initial and Boundary Conditions
As shown in Figure 8, a numerical model with scale 90 m × 90 m × 1 m was established. The excavation model was composed of 81,600 elements and radially defined as total three categories of unstructured meshing regions. Facing to the cross section, as innermost densified is a circular ring region with radii from 2.0 m to 7.6 m, which was evenly divided into 100 zones along radial direction and 640 zones along circumferential direction. Surrounding that is the transition region from a 7.6 m-radius circle to a square with side 20 m long, which was divided into 20 proportional enlarging zones along radial direction at ratio 1.1 and along circumferential direction 160 even zones. Then the outmost square stripe was divided into 30 zones along radial direction at proportion ratio 1.1 from half-side 10 to 45 m and on the other hand 160 zones along circumferential direction. Along thickness direction, no further discretization was undertaken so there is only one layer of zones. The normal displacements were fixed at the front and the back as well as the bottom, after vertical stress was applied to the top and the corresponding value Kx times of vertical stress to the both lateral sides, the model was run up to the first balance. In term of the inner boundary, the circular wall was at first applied with normal and shearing stresses exactly reaching the initial balance. Next, the tunneling advance which is in actual three-dimensional excavating process was virtually modelled as gradually relaxing both normal and shearing stresses by the same small percent. From their original amount equivalent to the in-situ stress effect, they were stepwise reduced to zero and at each step were applied to the interior boundary, i.e., the wall of a circular opening [38].

Numerical Modelling Simulator
The finite-difference simulator FLAC 3D was used to model the mechanical process. This commercial software was chosen because it is convenient to use and widely applicable across mining and Geo-engineering fields. The incremental plastic parameter is defined by FLAC 3D as: where ∆ε p 1 , ∆ε p m , ∆ε p 3 represent the major, mean, and minor principal plastic strain increment, respectively [52]. The focus of this modelling study is on those geologic conditions where the behavior of the engineering of interest is governed by continuum failure and damage of intact rock material rather than by deformation in distinct fragmented features. In other words, assuming deep rockmass of EGS-E is initially intact or massively jointed, i.e., its GSI index belongs to more than 75.

Model Geometry, Gridding, Initial and Boundary Conditions
As shown in Figure 8, a numerical model with scale 90 m × 90 m × 1 m was established. The excavation model was composed of 81,600 elements and radially defined as total three categories of unstructured meshing regions. Facing to the cross section, as innermost densified is a circular ring region with radii from 2.0 m to 7.6 m, which was evenly divided into 100 zones along radial direction and 640 zones along circumferential direction. Surrounding that is the transition region from a 7.6 m-radius circle to a square with side 20 m long, which was divided into 20 proportional enlarging zones along radial direction at ratio 1.1 and along circumferential direction 160 even zones. Then the outmost square stripe was divided into 30 zones along radial direction at proportion ratio 1.1 from half-side 10 to 45 m and on the other hand 160 zones along circumferential direction. Along thickness direction, no further discretization was undertaken so there is only one layer of zones. The normal displacements were fixed at the front and the back as well as the bottom, after vertical stress was applied to the top and the corresponding value K x times of vertical stress to the both lateral sides, the model was run up to the first balance. In term of the inner boundary, the circular wall was at first applied with normal and shearing stresses exactly reaching the initial balance. Next, the tunneling advance which is in actual threedimensional excavating process was virtually modelled as gradually relaxing both normal and shearing stresses by the same small percent. From their original amount equivalent to the in-situ stress effect, they were stepwise reduced to zero and at each step were applied to the interior boundary, i.e., the wall of a circular opening [38].
It was assumed that the buried depth of tunnels for EGS-E is 3000 m. Accordingly, vertical stress due to pressure from overloading rocks can be calculated as 78 MPa in case that gravitational acceleration is taken as 10 m/s 2 and mass density of rocks is averagely taken as 2600 kg/m 3 . Horizontal stress is calculated through multiplying a lateral pressure coefficient Kx to this vertical stress. Thus compressive 78 MPa were initially set along vertical direction. A few values of Kx timing 78 MPa were initially set as evenly compressive along horizontal direction.

Dependency of Material Mechanical Properties on Temperature Condition
According to the relevant experimental results [24,53], the temperature-dependent physical and mechanical parameters of granite as shown in Table 1 were incorporated into the advanced finite difference modelling workflows using FLAC 3D via an in-built interpolation routine. So, the strength, elastic deformation and plastic failure of Fujian Granite could be predicted over a wide range of temperature conditions. Accordingly, it can be applicable to studying deep geothermal energy extraction.  It was assumed that the buried depth of tunnels for EGS-E is 3000 m. Accordingly, vertical stress due to pressure from overloading rocks can be calculated as 78 MPa in case that gravitational acceleration is taken as 10 m/s 2 and mass density of rocks is averagely taken as 2600 kg/m 3 . Horizontal stress is calculated through multiplying a lateral pressure coefficient K x to this vertical stress. Thus compressive 78 MPa were initially set along vertical direction. A few values of K x timing 78 MPa were initially set as evenly compressive along horizontal direction.

Dependency of Material Mechanical Properties on Temperature Condition
According to the relevant experimental results [24,53], the temperature-dependent physical and mechanical parameters of granite as shown in Table 1 were incorporated into the advanced finite difference modelling workflows using FLAC 3D via an in-built interpolation routine. So, the strength, elastic deformation and plastic failure of Fujian Granite could be predicted over a wide range of temperature conditions. Accordingly, it can be applicable to studying deep geothermal energy extraction. The peak cohesion value can be calculated as [41]: where σ in-situ represents the unconfined strength of the in-situ rock. Typical values of σ in-situ range between 30% and 50% of the laboratory measured unconfined compressive strength (UCS). In the modelling, suggested value 41% was accordingly adopted [41]. In addition to the above parameters, as shown in Table 2, referring to the numerical simulations study on AECL-URL [38], both the friction angle for spalling limit and the friction angle of damage initiation were set unvarying under the various temperature conditions. This is because friction angles are not strikingly different when granite is under 400 • C [27,31]. Both cohesion and tensile strength were set to the minimum value 0.1 kPa. The three plastic strain limits for cohesion weakening, tension weakening and friction strengthening, respectively, came from the numerical simulation of AECL-URL as well.

Studying Tunneling Induced Excavation Damaged Zone (EDZ) in Depths
The partition of rock failure region surrounding a tunnel hole is based on how to switch friction angle to distinguish Highly Damaged Zone (HDZ) from Excavation Damage Zone. When ϕ of damaged rock falls into between ϕ i and ϕ m , inside the space defined by the principal stresses the corresponding state is between the following two envelops: the Damage Initiation (DI) one and the Spalling Limit (SL) one. The situation of material at this position belongs to EDZ. When ϕ attains to equal to ϕ m , the stress state intersects the SL envelope and the surrounding rock there belongs to HDZ.
In order to understand the scale and mechanism of surrounding rock failures caused by tunnel excavation during a EGS-E project, in addition to choosing the material parameters at 200 • C according to Table 2, the numerical modelling under three configurations of initial stresses was first performed. In these configurations, to represent the tectonic stress effect, the vertical stress as both initial condition and top boundary condition is kept constant and the lateral pressure coefficient K x applicable to initial condition as well as the left and right lateral boundaries were defined as 0.8, 1.0 and 1.2, respectively. This is to consider one hydrostatic initial stress case as well as two anisotropic ones, because the in-situ state is close to the hydrostatic pressure state when the depth is greater than 3000 m [54].
As shown in Figure 9, in term of hydrostatic initial stress case, i.e., K x = 1.0, the upper row shows four enlarged views of distribution of HDZ, EDZ and Excavation Influence Zone (EIZ). The former dark and middle blue represent different extent of rock failures and the light grey one represents the rock distinct elastic deformation. These contours were plotted by distinguishing the state of principle stresses at each nodes taking advantage of the above mentioned criterion indicated by fraction angle envelops. The lower row shows a series of enlarged views of quarter contours indicating progressive evolution of the major principal stress. The normal pressure 10 MPa applied to the inner boundary meant that  As the inner confining stresses were reduced due to excavating, the circumferential stress concentration was generated around the tunnel wall. When Mohr's circle for the stresses intersects the Damage Initiation envelope, the rock mass begins to yield, produce plastic strain, which is shown in Figure 9a as damage initiation. Along with the confinement was further reduced, strain localization occurs in some elements in vicinity of the tunnel wall, resulting in HDZ being initiated there. Following that stage, some fractures start mutually interacting. As shown in the fourth diagram within the upper row of Figure  9, prominent log-spiral shearing sliding can be eventually observed at the end of modelling. Besides extensile fractures due to the loss of confinement, shear cracks propagate, and the fractures coalesce and interact mutually. Around the tunnel, the partition rupture appears. Figure 10 shows the phenomenon of intermittent compressive stress distributed along radial horizontal direction within the surrounding rock of a tunnel.  As shown in Figure 11 is the effects of lateral pressure coefficient Kx on EDZ in the aspects of distribution scale, failure pattern, the maximum plastic strain and the maximum principal stress. It can be seen that Kx has distinct effects on the failure issues of surrounding rock. Although the existence of a pair of intersecting shear zones were observed in all cases, they mainly concentrated behand two sidewalls in case Kx = 0.8 and at the backs of roof and floor in case Kx = 1.2. They both illustrated that the damage zones mainly focus along the direction of the minimum principal stress according to their initial stress states. The case of Kx = 0.8 belongs to the typical fracture pattern when hard rock is subjected to high stress in depths, i.e., the two groups of intersecting fractures form a pair of V-shaped failures can be observed at two sides, while relatively intact rock can be sustained at the roof and the floor. At lateral pressure coefficient of Kx = 1.2, greater horizontal stress leads to the shear zones closest to the side walls from where shearing slides continues to prop- As the inner confining stresses were reduced due to excavating, the circumferential stress concentration was generated around the tunnel wall. When Mohr's circle for the stresses intersects the Damage Initiation envelope, the rock mass begins to yield, produce plastic strain, which is shown in Figure 9a as damage initiation. Along with the confinement was further reduced, strain localization occurs in some elements in vicinity of the tunnel wall, resulting in HDZ being initiated there. Following that stage, some fractures start mutually interacting. As shown in the fourth diagram within the upper row of Figure 9, prominent log-spiral shearing sliding can be eventually observed at the end of modelling.
Besides extensile fractures due to the loss of confinement, shear cracks propagate, and the fractures coalesce and interact mutually. Around the tunnel, the partition rupture appears. Figure 10 shows the phenomenon of intermittent compressive stress distributed along radial horizontal direction within the surrounding rock of a tunnel.
As shown in Figure 11 is the effects of lateral pressure coefficient K x on EDZ in the aspects of distribution scale, failure pattern, the maximum plastic strain and the maximum principal stress. It can be seen that K x has distinct effects on the failure issues of surrounding rock. Although the existence of a pair of intersecting shear zones were observed in all cases, they mainly concentrated behand two sidewalls in case K x = 0.8 and at the backs of roof and floor in case K x = 1.2. They both illustrated that the damage zones mainly focus along the direction of the minimum principal stress according to their initial stress states. The case of K x = 0.8 belongs to the typical fracture pattern when hard rock is subjected to high stress in depths, i.e., the two groups of intersecting fractures form a pair of V-shaped failures can be observed at two sides, while relatively intact rock can be sustained at the roof and the floor. At lateral pressure coefficient of K x = 1.2, greater horizontal stress leads to the shear zones closest to the side walls from where shearing slides continues to propagate until intersect each other and then bend up and down, respectively. As shown in Figure 11 is the effects of lateral pressure coefficient Kx on EDZ in the aspects of distribution scale, failure pattern, the maximum plastic strain and the maximum principal stress. It can be seen that Kx has distinct effects on the failure issues of surrounding rock. Although the existence of a pair of intersecting shear zones were observed in all cases, they mainly concentrated behand two sidewalls in case Kx = 0.8 and at the backs of roof and floor in case Kx = 1.2. They both illustrated that the damage zones mainly focus along the direction of the minimum principal stress according to their initial stress states. The case of Kx = 0.8 belongs to the typical fracture pattern when hard rock is subjected to high stress in depths, i.e., the two groups of intersecting fractures form a pair of V-shaped failures can be observed at two sides, while relatively intact rock can be sustained at the roof and the floor. At lateral pressure coefficient of Kx = 1.2, greater horizontal stress leads to the shear zones closest to the side walls from where shearing slides continues to propagate until intersect each other and then bend up and down, respectively. As shown in Figure 12, there had already been theoretical and laboratory modelled studies on the log-spiral shear-failure sliding curves [55][56][57] surrounding a circular tunnel. Tensile/extensional modes were not evident in these physical model studies, which might be due to weak properties and low level of confinement.
Energies 2021, 14, x FOR PEER REVIEW 14 of 20 Figure 11. Effect of various coefficients of lateral pressure on EDZs.
As shown in Figure 12, there had already been theoretical and laboratory modelled studies on the log-spiral shear-failure sliding curves [55][56][57] surrounding a circular tunnel. Tensile/extensional modes were not evident in these physical model studies, which might be due to weak properties and low level of confinement.

Temperature Distribution and Stress Field in EGS-E Tunnel after Ventilation Cooling
Since EGS-E engineering in depths is to be under high temperature environment, for the sake of achieving eligible working temperature for personnel and equipment, it should be necessary to take the initiative in cooling down the surrounding. According to the latest relevant study [22], it was found that at the entrance the heat exchange from surrounding rocks to inside tunnel airflow is the most intense. Due to the largest temperature difference between the tunnel wall and surrounding rock, there is taken as the research objective. Some basic assumptions were made as follows: (1) The thermophysical properties of the surrounding rock as well as the insulation liner are assumed to be homogeneous and isotropic constants; (2) The contribution of rock strain energy to temperature evolution is ignored during cooling process. The coupling between thermal and mechanical simulations is achieved through iterating mechanical calculations till converge after each sufficient small thermal time step, i.e., 60 s in this paper.
In the modelling, the numerical discretization of surrounding rock was kept the same as that in the previous section. As shown in Figure 13, a 10 cm thick thermal insulation liner with its corresponding mesh was appended onto the tunnel wall. It was evenly divided into two blocks along the radial direction and 640 blocks along the circumferential direction.

Temperature Distribution and Stress Field in EGS-E Tunnel after Ventilation Cooling
Since EGS-E engineering in depths is to be under high temperature environment, for the sake of achieving eligible working temperature for personnel and equipment, it should be necessary to take the initiative in cooling down the surrounding. According to the latest relevant study [22], it was found that at the entrance the heat exchange from surrounding rocks to inside tunnel airflow is the most intense. Due to the largest temperature difference between the tunnel wall and surrounding rock, there is taken as the research objective. Some basic assumptions were made as follows: (1) The thermophysical properties of the surrounding rock as well as the insulation liner are assumed to be homogeneous and isotropic constants; (2) The contribution of rock strain energy to temperature evolution is ignored during cooling process. The coupling between thermal and mechanical simulations is achieved through iterating mechanical calculations till converge after each sufficient small thermal time step, i.e., 60 s in this paper.
In the modelling, the numerical discretization of surrounding rock was kept the same as that in the previous section. As shown in Figure 13, a 10 cm thick thermal insulation liner with its corresponding mesh was appended onto the tunnel wall. It was evenly divided into two blocks along the radial direction and 640 blocks along the circumferential direction. In the mechanical aspect, the initial and boundary conditions of this modelling were inherited from the previous simulation in case Kx = 1.0. In the thermal aspect, the temperatures of surrounding rock and the thermal insulation liner were initially set as 200 °C and 30 °C, respectively. The airflow temperature is kept constant as Ta = 30 °C. The outer boundary of the model is Dirichlet condition at fixed 200 °C and the heat exchange of the inner boundary is Robin condition, i.e., between the tunnel and the airflow Equation (2) was satisfied as follows where is thermal conductivity, is the temperature of the surrounding rock at the tunnel wall, h is convective heat transfer coefficient, i.e., 60 W/(m 2 ·K) in this paper.
As shown in Table 3 are the thermophysical parameters of granite as the surrounding rock [58] and the material made of manufacturing thermal insulation liner [22], respectively. On the other hand, the mechanical parameters of CWFS model had been shown in Table 2. It was assumed that the supporting, lining and other constructive implementations are completed during a period of 14 days. In term of modelling this procedure, the simulated magnitude and covering area of EDZs were recorded once per physical hour.  Figure 14, the radial distributions of temperature across surrounding rock are different with and without installing thermal insulation layer (TIL). Due to the influence of low-temperature airflow inside the tunnel, temperature in surrounding rock continues to descend with ventilating time, and the area of low-temperature disturbance gradually expands. In case of no thermal insulation liner, the temperature of surrounding rock decreases up to 35.5 °C when ventilation time reaches 300 h. In case of installing thermal insulation liner, the temperature gradient within thermal insulation liner is much dramatic than that of the surrounding rock in case of bare tunnel. Obviously the rather small heat conductivity of thermal insulation liner makes it become a barrier. After 300 h Figure 13. The meshing of thermal insulation liner which is quarterly shown as its symmetry.
In the mechanical aspect, the initial and boundary conditions of this modelling were inherited from the previous simulation in case K x = 1.0. In the thermal aspect, the temperatures of surrounding rock and the thermal insulation liner were initially set as 200 • C and 30 • C, respectively. The airflow temperature is kept constant as T a = 30 • C. The outer boundary of the model is Dirichlet condition at fixed 200 • C and the heat exchange of the inner boundary is Robin condition, i.e., between the tunnel and the airflow Equation (2) was satisfied as follows where λ is thermal conductivity, T r is the temperature of the surrounding rock at the tunnel wall, h is convective heat transfer coefficient, i.e., 60 W/(m 2 ·K) in this paper. As shown in Table 3 are the thermophysical parameters of granite as the surrounding rock [58] and the material made of manufacturing thermal insulation liner [22], respectively. On the other hand, the mechanical parameters of CWFS model had been shown in Table  2. It was assumed that the supporting, lining and other constructive implementations are completed during a period of 14 days. In term of modelling this procedure, the simulated magnitude and covering area of EDZs were recorded once per physical hour. As shown in Figure 14, the radial distributions of temperature across surrounding rock are different with and without installing thermal insulation layer (TIL). Due to the influence of low-temperature airflow inside the tunnel, temperature in surrounding rock continues to descend with ventilating time, and the area of low-temperature disturbance gradually expands. In case of no thermal insulation liner, the temperature of surrounding rock decreases up to 35.5 • C when ventilation time reaches 300 h. In case of installing thermal insulation liner, the temperature gradient within thermal insulation liner is much dramatic than that of the surrounding rock in case of bare tunnel. Obviously the rather small heat conductivity of thermal insulation liner makes it become a barrier. After 300 h of ventilation, the lowest temperature of surrounding rock is still as high as 153.0 • C, indicating that the insulation layer (TIL) can effectively mitigate the heat loss from surrounding rock during the process of ventilation.  At the end of ventilation cooling, the stress state of surrounding rock along the radial but horizontal direction on the circular tunnel cross-section is exemplified in Figure 15. The falling off of surrounding rock temperature causes additional tensile stress near the tunnel. The reduction of minimum principal stress leads to further rock damages and then descending of the compressive stress which can hold its burden. The rock in the temperature lowered zone shrinks and is unloaded so the concentration of compressive stress transfers to further zone.   At the end of ventilation cooling, the stress state of surrounding rock along the radial but horizontal direction on the circular tunnel cross-section is exemplified in Figure 15. The falling off of surrounding rock temperature causes additional tensile stress near the tunnel. The reduction of minimum principal stress leads to further rock damages and then descending of the compressive stress which can hold its burden. The rock in the temperature lowered zone shrinks and is unloaded so the concentration of compressive stress transfers to further zone. of ventilation, the lowest temperature of surrounding rock is still as high as 153.0 °C, indicating that the insulation layer (TIL) can effectively mitigate the heat loss from surrounding rock during the process of ventilation.
(a) Lined tunnel with thermal insulation layer; ( b) Bare tunnel without thermal insulation layer. At the end of ventilation cooling, the stress state of surrounding rock along the radial but horizontal direction on the circular tunnel cross-section is exemplified in Figure 15. The falling off of surrounding rock temperature causes additional tensile stress near the tunnel. The reduction of minimum principal stress leads to further rock damages and then descending of the compressive stress which can hold its burden. The rock in the temperature lowered zone shrinks and is unloaded so the concentration of compressive stress transfers to further zone.    Figure 16 shows the scales of EDZs and their evolution process along with ventilation time. With the ventilation time passing by, EDZ d , HDZ d , EDZ a and HDZ a gradually grew. Finally cooling as well had a significant impact on the development of fracture zones in surrounding rock. Having undergone ventilation cooling for 14 days with and without TIL, the breadth of EDZ d increased by 0.87 m and 1.87 m, and the coverage of EDZ a increased by 27.7 m 2 and 58.28 m 2 , respectively. The expansion of EDZ scales with cooling time experienced three stages: (1) Slow expansion stage: in the very beginning, i.e., the initial stage of cooling, the rock mass affected by temperature change was concentrated inside HDZ in vicinity of circular boundary. The shrinkage of its volume only caused insignificant stress changes but led to a slow growth of EDZs. Without TIL, this stage corresponds to 0-10 h. While with TIL, this stage was put off to 2 d. (2) Rapid expansion stage: further cooling will reduce temperatures in EDZ as well as EIZ (Excavation Influence Zone) where majority of compressive stresses were withstood, resulting in rock volume shrinkage and additional tensile stress. Thus, the concentration of compressive stress quickly transferred to the further and EDZs rapidly grew into larger scales; (3) Attenuating expansion stage: along with further ventilation and cooling inside the tunnel, the temperature gradient of surrounding rock mitigated, displaying a relatively slow temperature decline. In the same time, non-uniform in the volume shrinkage of EDZs turned to relax. Finally, the increasing rate of EDZs scale slowed down. initial stage of cooling, the rock mass affected by temperature change was concentrated inside HDZ in vicinity of circular boundary. The shrinkage of its volume only caused insignificant stress changes but led to a slow growth of EDZs. Without TIL, this stage corresponds to 0-10 h. While with TIL, this stage was put off to 2 d. (2) Rapid expansion stage: further cooling will reduce temperatures in EDZ as well as EIZ (Excavation Influence Zone) where majority of compressive stresses were withstood, resulting in rock volume shrinkage and additional tensile stress. Thus, the concentration of compressive stress quickly transferred to the further and EDZs rapidly grew into larger scales; (3) Attenuating expansion stage: along with further ventilation and cooling inside the tunnel, the temperature gradient of surrounding rock mitigated, displaying a relatively slow temperature decline. In the same time, non-uniform in the volume shrinkage of EDZs turned to relax. Finally, the increasing rate of EDZs scale slowed down.

Conclusions
In this study, the surrounding rock failure mode and EDZs scale caused by excavation under the conditions of near hydrostatic pressure is analysed, using the numerical simulation method. Further evolution of EDZs at Kx = 1.0 induced by ventilation is studied by coupling temperature field. The following conclu-sions can be drawn: (1) The CWFS model can well represent the progressive failure of surrounding rock from damage to complete loss of strength during unloading, and the initiation and propagation of EDZs. Small confinement at tunnel wall plays an important role in restraining failure. (2) When the tunnel is excavated under near hydrostatic pressure (Kx = 0.8, 1.0 and 1. 2) in hard rock, shear fracture zone appears in all cases. Spalling failure occurs near the tunnel wall and the deep surrounding rock is controlled by shear fracture, which endan-ger the construction safety. Lateral pressure coefficients have a significant effect on the fracture morphology and scale of surrounding rock. When Kx = 0.8, 1.2, damage concentrates in the direction of the minimum principal stress, while the damage distribution is more uniform at Kx = 1.0, and when Kx = 1.0, 1.2, the EDZ scale is larger. (3) Due to the ventilation, the surrounding rock stress significantly changes. Additional tensile stress is generated in the shallow part of the surrounding rock, and the maximum and minimum principal stresses decrease. The compressive stress concentration transfers to the deep, which promotes the further expansion of EDZs, and it experiences slow, rapid and deceleration stages with ventilation time. Therefore, the influence of rock cooling volume shrinkage on EDZs must be considered in a High-

Conclusions
In this study, the surrounding rock failure mode and EDZs scale caused by excavation under the conditions of near hydrostatic pressure is analysed, using the numerical simulation method. Further evolution of EDZs at K x = 1.0 induced by ventilation is studied by coupling temperature field. The following conclu-sions can be drawn: (1) The CWFS model can well represent the progressive failure of surrounding rock from damage to complete loss of strength during unloading, and the initiation and propagation of EDZs. Small confinement at tunnel wall plays an important role in restraining failure. (2) When the tunnel is excavated under near hydrostatic pressure (Kx = 0.8, 1.0 and 1. 2) in hard rock, shear fracture zone appears in all cases. Spalling failure occurs near the tunnel wall and the deep surrounding rock is controlled by shear fracture, which endan-ger the construction safety. Lateral pressure coefficients have a significant effect on the fracture morphology and scale of surrounding rock. When Kx = 0.8, 1.2, damage concentrates in the direction of the minimum principal stress, while the damage distribution is more uniform at Kx = 1.0, and when Kx = 1.0, 1.2, the EDZ scale is larger. (3) Due to the ventilation, the surrounding rock stress significantly changes. Additional tensile stress is generated in the shallow part of the surrounding rock, and the maximum and minimum principal stresses decrease. The compressive stress concentration transfers to the deep, which promotes the further expansion of EDZs, and it experiences slow, rapid and deceleration stages with ventilation time. Therefore, the influence of rock cooling volume shrinkage on EDZs must be considered in a High-Temperature Tunnel. In addition, thermal insulation layer prolongs the slow growth stage.