Numerical Study on Application Conditions of Equivalent Continuum Method for Modeling Heat Transfer in Fractured Geothermal Reservoirs

: The equivalent continuum method an effective approach for modeling heat transfer in fractured geothermal reservoirs. However, presently there is a lack of systematical and profound study on application conditions of the equivalent porous media (EPM) method. In this study, we numerically investigated the application conditions of the EPM method based on geological data of Yangbajing geothermal ﬁeld. The results indicate that when fracture spacing is within 3–25 m, the results of the EPM method are basically in the same levels as those of the MINC method. However, when the fracture spacing is within 25–300 m, differences of the EPM method from the MINC method increase with the fracture spacing, so when the fracture spacing is within 25–300 m, it is unreasonable to adopt the EPM method to simulate the fractured reservoirs. With the fracture spacing increasing within 25–300 m, the system production temperature and electric power will gradually decrease; the injection pressure, reservoir impedance and pump power will gradually increase; and the energy efﬁciency will gradually decrease.


Background
Actual enhanced geothermal system (EGS) reservoirs are three-dimensional and interconnected fracture networks and are not a single fracture or parallel multiple fracture sets [1][2][3][4][5][6]. Presently there are mainly two different methods to describe the fracture system: equivalent continuum method (ECM) and discrete fracture network (DFN) method. The ECM treats the actual and discrete fracture network as an equivalent continuous porous medium, and this method is mainly used to model highly fractured systems. The ECMs mainly include the equivalent porous media (EPM) method or effective continuum method, the dual porosity media (DPM) method and the multiple interacting continua (MINC) method [7][8][9][10]. Because only few fracture distribution data are required and the computation efficiency is very high, ECMs have been widely used in fractured geothermal reservoir simulation, with the EPM method being most used [7][8][9][10]. The DFN method considers fracture spacing, fracture aperture and fracture orientation characteristics of the fracture system and establishes a discrete fracture network model based on the statistical properties. The DFN method requires many fracture distribution data, and its calculation efficiency is very low; in fact, although it has rarely been used, its use has recently been increasing in geothermal reservoir modeling because of the advancement of computational technology [4,5]. Though the DFN method takes more time, it usually provides more information on the response of the fracture system.
The three methods based on the equivalent continuum assumption have been widely used in geothermal engineering. Birdsell et al. used the EPM method to study performance characteristics of Fenton Hill EGS considering the coupling of water flow, heat transfer and tracer transport, and the results indicate that the reservoir impedance rises with time [11]. McDermott et al. employed the EPM to investigate the impacts of the coupling interaction of hydrologic-thermal-mechanical-chemical processes on the heat production performance of EGSs [12]. Watanabe et al. used the EPM to analyze the uncertainty of coupled hydrologic-thermal-mechanical processes in an EGS reservoir [13]. Zeng et al. used the EPM to investigate the heat production potential of water circulating through two horizontal wells in Desert Peak geothermal field [14]. Sanyal et al. used the DPM to study the electricity generation potential at Desert Peak geothermal field, and the results indicate that net generation profile versus time and reservoir heat recovery factors are the most appropriate criteria and that improving permeability without improving matrix-to-fracture heat transfer area has little benefit in heat recovery or net generation [15]. Taron et al. used the DPM to study the performance characteristics of EGS reservoirs considering the coupling of water flow, heat transfer, solid deformation and chemical reaction [16,17]. Gelet et al. used the DPM to investigate the performance characteristics of EGS reservoirs under the condition of thermodynamic nonequilibrium effect, and the results demonstrate that the difference between DPM and EPM decreases with the decrease in fracture spacing; the effect of thermodynamic nonequilibrium on the production temperature gradually decreases with the decrease in fracture spacing [18,19]. Pruess et al. used the MINC method to investigate the performance characteristics of EGS adopting CO 2 as working fluid, and the results show that CO 2 has lower viscosity and better heat mining performance than water [20,21]. Spycher et al. used the MINC method to investigate regularities of phase distribution in CO 2 -brine mixture [22]. Borgia et al. employed the MINC method to study the salt precipitation process in the fractures of a CO 2 EGS [23]. Sun et al. used the DFN method to investigate the joint influence of in situ stress and fracture network geometry on heat transfer in fractured geothermal reservoirs [24,25]. Davarpanah et al. performed an experimental investigation and mathematical modeling of gas diffusivity by carbon dioxide in fractured reservoirs [26] and used hybrid fuzzy approaches to analyze the hydraulic fracturing technique [27]. Hu et al. analyzed the thermodynamic effects of cycling carbon dioxide injectivity in fractured shale reservoirs [28].
Though the equivalent continuum methods have been widely used to model heat transfer in fractured geothermal reservoirs, there is still a lack of understanding of the application conditions of the EPM method and the DPM method [6]. Wu et al. researched the application conditions of the EPM method in detail and found that the EPM method is reliable only under the condition of local thermodynamic equilibrium between fluid and rock [7]. Long et al. researched the effects of fracture spacing, fracture aperture and fracture orientation on reservoir permeability and found that when the fracture is highly dense, apertures are constant rather than distributed, orientations are distributed rather than constant and the fracture systems behave more like porous media [8]. Robinson et al. analyzed the influence of fracture spacing on heat conduction and convection and found that when the fracture spacing is lower than 2-3 m, the assumption of local thermodynamic equilibrium between fluid and rock is reasonable [2,29]. Gelet et al. compared single porosity model with double porosity model and found that the effect of local thermodynamic nonequilibrium on heat production is gradually enhanced with the increase in the fracture spacing; when the fracture spacing is small, the results of the single porosity model are very close to those of the double porosity model, while with the increase in the fracture spacing, the difference between the single porosity model and the double porosity model gradually increases, thus under the condition of large fracture spacing, the local thermodynamic nonequilibrium between the fluid and rock must be taken into consideration [18,19]. With the increase in fracture spacing, the temperature gradient between the water and rock gradually increases, and thus the influence of the distance between water and rock on heat conduction must be taken into consideration, and the MINC method should be used to divide the rock matrix [20][21][22][23]29].
In order to investigate the influence of fracture spacing on heat conduction and convection in fractured geothermal reservoirs and determine the application conditions of the EPM method and the DPM method, in this study, we numerically investigated the application conditions of the EPM method and the DPM method with TOUGH2 codes based on the geological data from Yangbajing geothermal field. The Yangbajing geothermal field is the first high-temperature hydrothermal convective geothermal field in China [30][31][32]. In [30], detailed geological information on the Yangbajing geothermal field is reported. This will establish a solid foundation for choosing a reasonable fracture model to simulate the fractured geothermal reservoirs.
The studies presented in [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29] have two aspects of drawbacks. On one hand, the method of determining the water production rate of an EGS is not considered, and previous studies all neglected that the injection pressure must be lower than the minimum reservoir principal stress. On the other hand, these studies did not consider the appropriate method of selection for choosing between the EPM and DPM methods, and they neglected that fracture spacing has a significant impact on the simulation results. Thus, there are two main novel points in this study: first, we considered that the injection pressure must be lower than an upper limit, thus allowing the optimal water production rate to be determined; second, we considered that the fracture spacing has obvious influence on the heat transfer processes, and thus the MINC method was used as a standard to determine the application ranges of the EPM method.

Research Objectives
This work mainly aims to investigate the applicable fracture spacing of the EPM method and the DPM method based on the MINC method and to establish a foundation for fractured reservoir representation and simulation.

The EPM Method
In the EPM method, the fractures in the medium are assumed to be highly random and interconnected, so that it is possible to define average properties at each point in the sense of statistics. This method does not consider the physical structure of every single fracture, and the whole fractured medium is treated as an equivalent continuous porous medium. The advantages of this method are that we can directly use mature seepage theory in porous medium and the computation efficiency is very high. The disadvantages of this method include three aspects: (1) the application conditions are usually limited because not all the fractured medium can be treated as a continuous porous medium; (2) the size of the representative elementary volume and the equivalent hydraulic parameters are very difficult to determine; (3) treating the fractured rocks as a continuous porous medium does not well describe the special water conductivity ability of the fractures [33]. Unlike in the rock matrix, in the fractures, the water can freely flow and rapidly be conducted from one site to the other site. The conservation equations of mass, momentum and energy of the EPM method are listed in [33,34].

The DPM Method
In the DPM method, the fractures are considered as the main water and heat flowing path, and the permeability of fractures is much higher than that of the rock matrix; the rock matrix is considered as the main water and heat storage space, while the porosity of the rock matrix is much higher than that of the fractures. Thus, the fractures and rock matrix form two different hydrodynamic systems that are both independent and interrelated [33]. Both hydrodynamic systems are treated as a continuum, and each system has its own porosity and permeability. Therefore, in each space point, we can define two porosities and two permeabilities in the DPM method, and thus there are two pressures and two velocities in each space point. The advantages of the DPM method lie in that the interchange of water and heat between the two hydrodynamic systems is taken into consideration, and the disadvantages of the DPM method lie in that the fractured rock matrix is assumed to have the same size and shape, which is too simplified and is very easy for the discretization of the problem domain [33]. The conservation equations of mass, momentum and energy of the DPM method are listed in [33,34].

The MINC Method
The transient periods for interporosity flow for coupled fluid and heat flows or multiphase flows can be very long. It is necessary to resolve the driving temperature, pressure and mass fraction gradients at the matrix-fracture interface in order to accurately describe such flows. Resolution of these gradients is achieved by appropriate subgridding of the matrix blocks in the MINC method [34]. For detailed information about the MINC method, the reader can refer to [35,36].

Governing Equations
The numerical studies on heat production and electricity generation through vertical wells or horizontal wells at Yangbajing geothermal field show that the performance of using horizontal wells to mine heat from the fractured granite reservoir is more excellent [30,[37][38][39]. Therefore, in this study, we take the horizontal well system as an example to investigate the applicable fracture spacing of the EPM method and DPM method at Yangbajing geothermal field. One injection well is located at bottom of the reservoir, and two adjacent production wells are located at top of the reservoir. The injected cold water enters into the fractured reservoir from the injection well, is heated by the fractured rock matrix and then is produced out from the production well, as shown in Figure 1. A suction pump is installed on the ground to draw the heated thermal water out from the reservoir, and the bottomhole production pressure is maintained as constant. We installed an injection pump on the ground directly connected with the injection well, and the water injection rate was maintained as constant [14,38]. This kind of production scheme can effectively lower reservoir pressure and flowing impedance and reduce water loss [14,38]. Because of the thermal contraction of the rock, during the production time, the fracture aperture will increase with time. However, degradation of fractures due to chemical scaling will cause a decrease in the fracture aperture. Previous studies have proved that these two opposite effects are equally possible [15]. Therefore, in this work, we have assumed that after stimulation the fracture characteristics remain unchanged over the heat production time, and thus the fracture propagation and fracture fluid effects are neglected in this study [15]. For simplification, we neglected the water loss in this study and assumed that the water production rate is equal to the water injection rate: q = q inj = q pro [14]. The governing equation for mass conservation is as follows: where ρ is water density, φ is reservoir porosity and V is velocity vector. The governing equation for momentum conservation in porous media is the famous Darcy law, which is as follows: where V x , V y and V z are velocity components along x, y and z directions, respectively. K x , K y and K z are reservoir permeability along x, y and z directions. µ is water dynamic viscosity, p is pressure and g = 9.80 m/s 2 is gravity.
The governing equation for energy conservation in porous media is as follows: where T is temperature, (ρc p ) s is the product of solid density and solid specific heat capacity and (ρc p ) f is the product of water density and water specific heat capacity. c p is specific heat capacity, the k s is solid heat conductivity and k f is water heat conductivity.

Domain, Grid and Parameters
In this study, we used the TOUGH2-EOS1 codes to simulate the heat production and electricity generation process [34]. TOUGH2 software is a numerical simulator for nonisothermal flows of multicomponent, multiphase fluids in one-, two-and three-dimensional porous and fractured media. This software used the integral finite difference method to solve the conservation equations of mass, momentum and energy. The accuracy of the software has been proved by comparison with many different analytical and numerical solutions, and thus it has been widely used in geothermal reservoir engineering [34]. The study domain is located in a depth of 950-1350 m, corresponding to the fractured granite reservoir at Yangbajing geothermal field. Water supply from the deep formation and the weak heat conduction between the impermeable cap rock and the permeable reservoir are neglected in this study [23,30]. There are two production wells around the injection well, and the horizontal distance between the injection well and the production well is 500 m, as shown in Figure 2. Because of symmetry, only half of the domain in Figure 2 needs to be simulated, shown as the gray area in Figure 2. Assuming the whole length of the horizontal wells is 500 m, only a 10 m length along the well needs to be simulated, and thus the water production rate of the simulated domain is 50 times that of the 10 m length domain. The length of the simulated domain is 500 m along the x-direction, 10 m along the y-direction and 400 m along the z-direction. If the water production rate of the 10 m length domain is q, then the whole water production rate of the three horizontal well system is Q: Q = q × 50 × 2 = 100 q. For the EPM method, it is uniformly divided into 25 grid blocks along the x-direction, and the width of every grid block is 20 m; there is only one grid block along the y-direction, and the width of this grid block is 10 m; it is uniformly divided into 20 grid blocks along the z-direction, and height of every grid block is 20 m. Therefore the whole two-dimensional system in Figure 3 comprises 25 × 1 × 20 = 500 grid blocks. When the fracture spacing is less than 2-3 m, the assumption of local thermodynamic equilibrium between the water and rock is reasonable, and under this condition, the computation results of the EPM method are in the same levels as those of the MINC method [14]. For the MINC method, we divided each grid block of the EPM method into 5 continua with volume fractions of 0.0060, 0.0663, 0.1325, 0.2651 and 0.5301, respectively. Therefore, the whole system for the MINC method has 25 × 1 × 20 × 5 = 2500 grid blocks. In this study, for the MINC method, we investigated the production performance and efficiency under the condition of fracture spacing of 3, 5, 10, 25, 50, 100, 200 and 300 m [14]. Figure 3 shows the simulated domain and the corresponding grid used in the simulation as stated above. The injection well is located at depth of z = −1330 m, and the production wells are located at depth of z = −970 m. Because the radius of the horizontal well is only r w = 0.10 m, far smaller than the size of the grid block, the well is represented by the grid blocks in which the well is located in the form of source term [34].    Table 1 [30], and other parameters of the rock matrix and fractures can be found in [23]. For the EPM method, there is only one porosity and one permeability; thus, the porosity φ and permeability k of the fractured reservoir are determined by Equations (1) and (2), where φ m and k m are the porosity and permeability of the rock matrix, respectively, and φ f and k f is the porosity and permeability of the fractures, respectively.

Boundary and Initial Conditions
The saturated vapor pressure corresponding to 248 • C is 3.84 MPa, which is far lower than the bottomhole pressure at the production well, which is 5.00 Mpa, so water in the wells and reservoir is all kept as liquid [30]. The initial reservoir temperature is 248 • C, and the initial reservoir pressure is P = −0.0089z − 0.4444 (MPa). As stated above, the slight conductive heat exchange between the impermeable cap rock or base rock and the permeable reservoir is neglected, so the boundaries at the top and bottom of the reservoir in Figure 3 are no-flow for mass and heat. Due to symmetry, the left and right boundaries in Figure 3 are also no-flow for mass and heat.

The Determination of Water Injection Rate
The injection of cold water through the injection well into the fractured reservoir will increase the pressure at the wells and the reservoir; too high pressure will cause fracture extension and reservoir damage, so the bottomhole injection pressure P inj must be kept lower than a limit P max . Based on the study of [14], in this work we assume the following: where P max = fP w0 ; f = 1.2 is the safety factor, which is determined by the actual geological conditions of the reservoir [14]; P w0 is the initial pressure of the injection wellbore. In this work, the injection well is located at z = −1330 m, and the corresponding initial pressure is P w0 = 11.39 MPa; thus, P max = 13.67 MPa. The production wells are located at z = −970 m, and the corresponding initial pressure is 8.19 MPa; the bottomhole production pressure is kept at P 0 = 5.00 MPa to maintain production. The injection temperature is T inj = 60 • C [14], and the corresponding injection specific enthalpy is about h inj = 261.20 kJ/kg. Based on (3), the maximum water injection rate of the simulated domain is q = 2.0 kg/s, and the whole water injection rate of the three horizontal well system is Q = 100q = 200.0 kg/s. Figure 4 shows the evolution of production temperature T pro for the EPM method and the MINC method. For the MINC method, the investigated fracture spacings are 3, 5, 10, 25, 50, 100, 200 and 300 m, while the other conditions are the same as those for the EPM method. The corresponding relationships of the porosity and permeability between the EPM method and MINC method are determined by Equations (1) and (2). Figure 4 indicates that when using the EPM method, the T pro during the 20 years is within 248.00-156.96 • C. Using the MINC method, when the fracture spacing is 3, 5, 10 and 25 m, the T pro is within 248.00-153.51 • C, 248.00-153.51 • C, 248.00-153.48 • C and 248.00-152.33 • C, respectively. Therefore, when the fracture spacing D is within 3-25 m, the T pro values found with the EPM method are basically in the same levels as those of the MINC method. When D increases to be 50, 100, 200 and 300 m, the corresponding T pro during the 20 years is within 248.00-147.82 • C, 248.00-136.99 • C, 248.00-112.82 • C and 248.00-99.11 • C, respectively. Therefore, when fracture spacing is within 25-300 m, with the increase in the fracture spacing, the difference in T pro between the EPM method and the MINC method gradually increases; thus, when D is larger than 25 m, we should use the MINC method to simulate the fractured reservoirs, and the EPM method is not reasonable under this condition. Besides, Figure 4 shows that with the increase in D, the duration of the stable stage gradually decreases, the duration of the declining stage gradually increases and the declining rate of T pro during the declining stage gradually increases. This is because higher D effectively reduces the heat transfer area between the water and the rock, and this will obviously decrease the thermal power according to the heat transfer formula of convection and thus decrease the T pro . This is consistent with results from Sanyal et al. using the DPM method [15].

Electric Power
The electric power W e is calculated by the equation in [14], and the heat rejection temperature of the Yangbajing geothermal power plant is T 0 = 9 • C = 282.15 K [25]. Figure 5 shows the evolution of the electric power W e of the EPM method and the MINC method. During the 20 years, using the EPM method, the W e is within 33.61-12.50 MW. Using the MINC method, when the fracture spacing is 3, 5, 10 and 25 m, the W e is within 33.62-11.86 MW, 33.62-11.86 MW, 33.63-11.85 MW and 33.6-11.64 MW, respectively. Therefore when D is within 3-25 m, the difference in W e between the EPM method and the MINC method is very small; thus, under this condition, we can use the EPM method to simulate the fractured geothermal reservoirs. When D increases to be 50, 100, 200 and 300 m, the corresponding W e during the 20 years is within 33.63-10.82 MW, 33.63-8.93 MW, 33.63-5.22 MW and 33.63-3.44 MW, respectively. Therefore, when D is within 25-300 m, with the increase in D, the difference in W e between the EPM method and the MINC method gradually increases; thus, when D is within 25-300 m, it is unreasonable to use the EPM method to model the fractured geothermal reservoirs. With the increase in fracture spacing, the heat conduction in the rock mass becomes predominant and the heat convection in the fracture water becomes secondary. The MINC method can carefully consider the heat conduction in the rock mass; in this method, the pressure and temperature changes in the matrix are controlled locally by the distance from the fractures. Thus, for greater fracture spacing, the MINC method is more suitable. For smaller fracture spacing, the heat convection in the fracture water is predominant and the heat conduction is secondary; thus, both the EPM method and the MINC method are suitable. Besides, Figure 5 also shows that during the stable stage, the W e basically remains unchanged with the increase in D; during the declining stage, the rate of decline in W e gradually increases with the increase in D. As stated above, this is because higher D obviously reduces the heat transfer area between the rock and water, and this will obviously decrease the thermal power and the electric power [14].  Figure 6 shows the evolution of the injection pressure P inj of the EPM method and the MINC method. Using the EPM method, during the 20 years, the injection pressure P inj is within 9.81-13.42 MPa. Using the MINC method, when the fracture spacing is 3, 5, 10 and 25 m, the P inj is within 9.78-13.52 MPa, 9.80-13.52 MPa, 9.85-13.51 MPa and 10.03-13.45 MPa, respectively. Therefore, when D is within 3-35 m, the results of the EPM method are in the same levels as those of the MINC method, and under this condition, we can use the EPM method to simulate the fractured geothermal reservoirs. However, when D increases to be 50, 100, 200 and 300 m, the corresponding P inj during the 20 years is within 10.27-13.36 MPa, 10.42-13.26 MPa, 10.36-13.55 MPa and 10.15-13.84 MPa, respectively. Therefore, when D is within 25-300 m, the difference in P inj between the EPM method and the MINC method increases with the increase in D; thus, under this condition, it is unreasonable to use the EPM method to model the fractured geothermal reservoirs. Besides, Figure 6 shows that the P inj gradually increases with time. This is because as reservoir temperature gradually decreases, the water temperature in the reservoir also decreases with time, and this causes a gradual increase in water viscosity and an increase in reservoir impedance. This effect has been clearly stated in a previous study from Zeng et al. [14].

Reservoir Impedance
The reservoir impedance I R is calculated by the equation in [14]. I R represents the power consumption of the unit production rate for penetrating the fractured reservoir, P pro = P 0 = 5.00 MPa is the bottomhole production pressure and q 0 = 50q = 100.0 kg/s is the water production rate between the injection well and the production well. Figure 7 shows the evolution of the reservoir impedance I R of the EPM method and the MINC method. Using the EPM method, during the 20 years, the I R is within 0.048-0.084 MPa/(kg/s). Using the MINC method, when the fracture spacing is 3, 5, 10 and 25 m, the I R is within 0.048-0.085 MPa/(kg/s), 0.048-0.085 MPa/(kg/s), 0.049-0.085 MPa/(kg/s) and 0.050-0.085 MPa/(kg/s), respectively. Therefore, when D is within 3-25 m, the results of the EPM method are basically in the same levels as those of the MINC method, and under this condition, we can use the EPM method to model the fractured geothermal reservoirs. However, when D increases to be 50, 100, 200 and 300 m, the corresponding I R during the 20 years is within 0.053-0.084 MPa/(kg/s), 0.05-0.083 MPa/(kg/s), 0.054-0.085 MPa/(kg/s) and 0.051-0.088 MPa/(kg/s), respectively. Therefore when D is within 25-300 m, with the increase in D, the difference in I R between the EPM method and the MINC method gradually increases, and under this condition, it is unreasonable to use the EPM method to model the fractured geothermal reservoirs. Besides, as mentioned above, Figure 7 also shows that during the 20 years the I R gradually increases with time, which is caused by the increase in water viscosity caused by the decline in reservoir temperature [14].

Pump Power
The internal energy consumption of the pump power can be calculated by the equation in [14]. The water density ρ varies with pressure and temperature, and this will significantly impact the calculation of Wp [14]. The average value of the minimum density and maximum density is used for computations in this study [14]. The maximum density is 989.0 kg/m 3 , and the minimum density is 804.0 kg/m 3 ; thus, the average value of ρ = 896.5 kg/m 3 is adopted for calculations. Figure 8 shows the evolution of the pump power Wp of the EPM method and the MINC method. Using the EPM method, during the 20 years, the Wp is within 0.46-1.47 MW. Using the MINC method, when the fracture spacing is 3, 5, 10 and 25 m, the Wp is within 0.45-1.49 MW, 0.46-1.49 MW, 0.47-1.49 MW and 0.52-1.47 MW, respectively. Therefore, when D is within 3-25 m, the results of the EPM method are in the same levels as those of the MINC method; thus, under this condition, the EPM method can be used to model the fractured geothermal reservoirs. However, when D increases to be 50, 100, 200 and 300 m, the corresponding Wp during the 20 years is within 0.59-1.45 MW, 0.63-1.42 MW, 0.61-1.50 MW and 0.56-1.58 MW, respectively. Therefore, when D is within 25-300 m, with the increase in D, the difference in Wp between the EPM method and the MINC method gradually increases, and under this condition, it is unreasonable to use the EPM method to model the fractured geothermal reservoirs.
Besides, Figure 8 also shows that during the 20 years the Wp gradually increases with time. Based on the equation of Wp in [14], this is caused by the gradual increase in P inj , and this is in accordance with previous studies [14].

Energy Efficiency
The energy efficiency η of the system can be calculated by the equation in [14]. As stated above, the average value of the water density of ρ = 896.5 kg/m 3 is adopted in the calculation of η. Figure 9 shows the evolution of the energy efficiency η of the EPM method and the MINC method. Using the EPM method, during the 20 years, the η is within 73. 28-8.53. Using the MINC method, when the fracture spacing is 3, 5, 10 and 25 m, the η is within 74.30-7.94, 73.59-7.94, 71.33-7.95 and 64.61-7.89, respectively. Therefore, when D is within 3-25 m, the results of the EMP method are in accordance with those of the MINC method; under this condition, we can use the EPM method to simulate the fractured geothermal reservoirs. However, when D increases to be 50, 100, 200 and 300 m, the corresponding η during the 20 years is within 57.08-7.46, 53.33-6.29, 54.80-3.48 and 60.82-2.17, respectively. Therefore, when D is within 25-300 m, with the increase in D, the difference in η between the EPM method and the MINC method gradually increases; under this condition, it is unreasonable to use the EPM method to model the fractured geothermal reservoirs. Besides, Figure 9 also shows that during 20 years the η gradually decreases. Based on the equation of η in [14], this is caused by the gradual decline in T pro and the gradual increase in P inj , and this is in accordance with a previous study from Zeng et al. [14].

Applicable Fracture Spacing of the DPM Method
Because the DPM method considers fracture spacing, fracture porosity and permeability and rock matrix porosity and permeability, the DPM method can be used to model the fluid flow and heat transfer process in a fractured porous medium under all kinds of fracture spacing, just like the MINC method. Sanyal et al. used the DPM method to investigate production performance and efficiency characteristics of the EGS at Desert Peak geothermal field and analyzed the influence of fracture spacing on electricity generation for fracture spacing within 0.33-300 m [15]. The results show that increasing permeability without increasing matrix-to-fracture heat transfer area has little benefit in net generation or heat recovery; an increase in stimulation volume increases generation level without obviously impacting the shape of the generation profile; and the recovery factor can be reasonably correlated with stimulation volume, irrespectively of well geometry, fracture spacing and fracture domain permeability. Because Sanyal et al. discussed the influence of fracture spacing on the heat production performance and efficiency of EGS reservoir in [15], we directly refer to the corresponding conclusions in this study. For more information about the applicable fracture spacing in the DPM method, please refer to [15].

Limitation of the Model
There are three main limitations of the model in this work. First, the water loss is neglected. In actual EGS reservoirs, the water loss is ever-present, and neglecting water loss will improve the production performance of the EGS. Second, the rock deformation due to pore elasticity and thermal contraction of rocks is neglected. Third, the chemical dissolution and precipitation of rock minerals are neglected. The second and third effects on the production performance of EGSs need further study in the future.

Conclusions
In this study, we mainly investigated the influence of the fracture spacing on the electricity generation performance and efficiency for fracturing spacing within 3-300 m by means of numerical simulation, and we discussed the applicable fracture spacing for the EPM method based on the geological data from Yangbajing geothermal field. According to the simulation results, we can draw the following conclusions: (1) Under the reference conditions, the applicable fracture spacing for the EPM method is within 3-25 m. Within this range, with the increase in the fracture spacing, the production temperature, electric power, injection pressure, reservoir impedance, pump power and energy efficiency are only very slightly influenced, and the results of the EPM method are in accordance with those of the MINC method. (2) When the fracture spacing is larger than 25 m, with the increase in the fracture spacing, the difference in the simulated results between the EPM method and the MINC method gradually increases; under this condition, it is unreasonable to use the EPM method to model the fractured geothermal reservoirs, and the DPM method or the MINC method should be used. (3) When the fracture spacing is within 25-300 m, with the increase in the fracture spacing, the heat exchange area between the rock and water gradually decreases; the thermal power and electric power gradually decline; the injection pressure, reservoir impedance and pump power gradually increase; and the energy efficiency gradually decreases. (4) The above conclusions are only valid for the current setup, and although it may be possible to transfer the conclusions to other fields, any such generalization should be done with caution. Institute of Marine Geology, China Geological Survey, which significantly improved the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
C p specific heat capacity, J/(kg· • C) D fracture spacing, m g gravity, 9.80 m/s 2 h well depth, m h 1 depth of injection well, m h 2 depth of production well, m h inj injection specific enthalpy, kJ/kg hpro production specific enthalpy, kJ/kg I R reservoir impedance, MPa/(kg/s) k reservoir permeability, m 2 k f fracture permeability, m 2 k m matrix permeability, m 2 k x intrinsic permeability along x, m 2 k y intrinsic permeability along y, m 2 k z intrinsic permeability along z, m 2 P pressure, MPa P max critical pressure, MPa Pinj injection pressure, MPa Ppro production pressure, MPa P 0 bottomhole production pressure, MPa q water production rate, kg/s Q total water production rate, kg/s T temperature, • C T 0 mean heat rejection temperature, 282.15 K Tpro production temperature, • C V velocity vector, m/s Wp electric power of pump, MW W e electric power, MW x, y, z Cartesian coordinates, m φ reservoir porosity η energy efficiency η p pump efficiency, 80% ρ water density, kg/m 3 µ water dynamic viscosity, Pa·s