Effect of Vertical Permeability Heterogeneity in Stratiﬁed Formation on Electricity Generation Performance of Enhanced Geothermal System

: Because geologic sedimentation and hydrofracturing processes are not homogeneous, the reservoirs of enhanced geothermal systems (EGSs) are also heterogeneous; this has a signiﬁcant inﬂuence on the electricity generation performance of EGS. Presently, there are a lack of systematic and profound studies on the effect of vertical permeability heterogeneity in stratiﬁed formation on the electricity generation performance of EGS. In order to uncover the effect of vertical permeability heterogeneity on electricity generation performance of EGS, in this work we analyzed the inﬂuence of vertical permeability heterogeneity on electricity generation performance of EGS through a numerical method based on geological data at the Yangbajing geothermal ﬁeld. The results indicate that when the average permeability of stratiﬁed formations is constant for a homogeneous reservoir, the system attains maximum water production rate, maximum electric power, minimum reservoir impedance and maximum pump power; with the increasing of the vertical permeability heterogeneity, the water production rate gradually decreases, the electric power gradually declines, the reservoir impedance gradually increases and the pump power gradually declines. When the average permeability of stratiﬁed formations is constant, with the increasing of the vertical permeability heterogeneity, the injection pressure and energy efﬁciency only changes very slightly; this indicates that the vertical permeability heterogeneity is not the main factor affecting the system injection pressure and energy efﬁciency. MPa, thus the minimum pressure at the bottomhole of the production well is about (8.23 − 3.40) MPa = 4.83 MPa. In this study, we assume that the production well is kept at P 0 = 5.0 MPa to maintain the heat production. Based on (3), the maximum water production rate can be determined. The injection temperature is T inj = 60 ◦ C [4,5]; for an initial wellbore pressure of P W0 = 11.35 MPa, the injection speciﬁc enthalpy is about h inj = 260.66 kJ/kg.


Background
The EGS geothermal resource is the main body of the whole geothermal resource. The total EGS resource within a 10 km depth that can be exploited occupies more than 90% of the geothermal resource [1], so only the development and utilization of the EGS resource opens the inventory of the geothermal resource. Compared with other renewable energy sources, the EGS resource is more concentrated and stable and very suitable for generating base-load electric power, with a high utilization efficiency and nearly no pollution emission [1,2]. Therefore, the development of EGS resources has received worldwide attention. The total EGS resource reserve in China within a 3-10 km depth amounts to 20.90 M EJ; if the recoverable fraction is taken as 2%, the recoverable EGS resource amounts to 4400 times the total annual energy consumption in 2010 in China [3]. It is predicted that humans will utilize EGS on a large scale to generate electricity by 2030, and EGS will provide about 100,000 MW electric power by 2050 in the U.S.A., occupying about 10% of the total electricity generating capacity [1].
Because numerical studies of EGS reservoir are timesaving, inexpensive and relatively easy, numerical simulations are of great significance for analyzing the characteristics of the long-term electricity generation performance of EGS [1,2]. Zeng et al. investigated the performance and efficiency characteristics of EGSs through vertical wells and horizontal wells at the Desert Peak geothermal field, and found that reservoir impedance through a single vertical fracture is within 1.14-1.93 MPa/(kg/s), which is far higher than the commercial criterion of 0.1 MPa/(kg/s); the reservoir impedance through two horizontal wells is within 0.10-0.13 MPa/(kg/s), which is very close to the commercial criterion, so adopting two horizontal wells can obtain better heat production and electricity generation performance [4,5]. Based on the geological data at Yangbajing geothermal field, Zeng et al. numerically studied the electricity generation performance of EGSs through a single horizontal well, a single vertical well, multi vertical wells and multi horizontal wells, respectively, and found that the heat recovery ratio of a single well system is very low and adopting multi wells can obtain a higher heat recovery ratio [2,[6][7][8][9]. Zhang et al. researched the electricity generation performance and efficiency characteristics of EGS through two horizontal wells at the Xujiaweizi area at the Daqing oilfield, and found that during the past 20 years, the reservoir impedance has been within 0.17-1.49 MPa/(kg/s), which is far higher than the commercial criterion of 0.1 MPa/(kg/s), so the horizontal well technology still needs to be improved [10]. Zhang et al. numerically investigated the electricity generation and heating potential of the enhanced geothermal system in Songliao Basin, and the results indicate that Yingcheng formation is suitable for the development of HDR (hot dry rock) resources in this region [11]. Zhang et al. numerically evaluated the potential of developing an enhanced geothermal heating system in northeast China, and the results demonstrate that thermal energy production from the stimulated reservoir is satisfactory for heating purposes, but the water flow impedance maintains at a relatively high level during the production period [12]. Huang et al. forecasted the heat extraction and electricity production of a prospective enhanced geothermal system site in Songliao Basin, and the results suggest that the desirable thermal efficiency and gross power output could be obtained initially, whereas the decrease in production arising from thermal depletion of the reservoir is significant at later stages of the operation [13]. Huang et al. conducted a parametric study of an enhanced geothermal system based on thermo-hydro-mechanical modeling of a prospective site in Songliao Basin, and the results show that the properties that are able to strengthen the heat convection in the reservoir are able to influence the power output significantly, and a wellbore arrangement that takes advantages of the buoyancy flow is also a valuable way to increase the production [14].
Besides simulations in the geothermal field scale as stated above, the theory for simulating fractured EGS reservoirs have made great progress in recent years [15][16][17][18]. First, besides the equivalent porous media (EPM) method to represent the fractured reservoirs, the double porosity method (DPM) [19][20][21][22][23], the multiple interacting continua (MINC) method [24][25][26] and the discrete fracture network (DFN) method [27][28][29][30][31][32] are increasingly used. Sanyal et al. used the DPM method to investigate the productivity of typical EGS, and found that fracture spacing has a significant influence on the electric power [19]. Taron et al. established a thermal-hydrologic-mechanical-chemical model based on the DPM method, and this greatly improved the simulation accuracy for EGS reservoirs [20,21]. Gelet et al. employed the DPM method to study the heat production performance of EGS and found that the single porosity approach overestimates the thermal contraction of fractured reservoirs [22,23]. Pruess et al. used the MINC method to study the production behavior of enhanced geothermal systems with CO 2 as the working fluid, and the results showed strong effects of gravity on the mass flow and heat extraction due to the large contrast of CO 2 density between cold injection and hot production conditions [24,25]. Zeng et al. investigated the electricity generation potential at the Yangbajing geothermal field through the MINC method, and found that the system attained an electric power of 22.2-19.3 MW and an energy efficiency of 23.92-10.82 during 30 years through five vertical wells [26]. Because the DFN method can accurately reflect the heterogeneity and anisotropy of the fractured reservoirs, it has become a hot spot for EGS reservoir modeling, such as the studies reported in [27][28][29][30][31][32]. Second, traditional EGS modeling mainly considers the coupling hydrological-thermal effect, while in recent years, the three coupling effects of the hydrological-thermal-mechanical model or four coupling effects of the hydrological-thermal-mechanical-chemical model are considered more and more in simulations [15][16][17][18][27][28][29][30][31]. Pandey et al. summarized the latest advances in the hydrological-thermal-mechanical-chemical model for EGS, and found that the poroelastic effects are important for short time periods (hours to days), the thermoelastic effects are important for intermediate time periods (days to months) and the chemical effects are important for long time period (months to years) [18]. Salimzadeh et al. presented a coupled hydraulic-thermal-mechanical-chemical model for EGS reservoirs and studied the aperture change under in situ stresses during heat production [27]. Shi et al. investigated heat extraction performance of a multilateral-well EGS based on a hydraulic-thermal-mechanical model and found that reservoir properties have significant effects on rock deformation [28]. Wang et al. proposed a strongly coupled hydro-thermomechanical model to simulate the heat production process in Habanero EGS and the results demonstrate that the proposed model is robust and effective [29]. Sun et al. studied heat extraction in EGS with the hydraulic-thermal-mechanical model based on the discrete fracture model, and the results indicate that it is important to take into account the HTM coupling effects [30]. Yao et al. used a hydraulic-thermal-mechanical model to simulate the Desert Peak geothermal reservoir and the results show that the electric power can meet the commercial standard [31]. Third, different from the traditional thermal model which assumes that water and rock are in local thermodynamic equilibrium, recently, the local thermodynamic non-equilibrium model isapplied more and more [33][34][35]. Jiang considered the actual existence of a local thermal non-equilibrium between the rock matrix and flowing fluid in the porous heat reservoir, and established a three-dimensional transient model for the EGS subsurface thermo-hydraulic process [33,34]. Shaik et al. simulated fluid-rock coupling heat transfer in a naturally fractured geothermal system that considered the local thermal non-equilibrium effect, and the results show that the heat transfer between rock and fluid and the fracture connectivity have a profound effect on the economic potential of a geothermal reservoir [35].
Presently, in most EGS reservoir models, the reservoir is assumed as homogeneous, and the effect of the permeability heterogeneity of the layered reservoirs on the heat production and electricity generation is neglected [1]. In fact, because the geologic sedimentation and hydrofracturing processes are not homogeneous and the EGS reservoirs are also layered and heterogeneous, the influence of the heterogeneity of reservoirs on production performance should be taken into consideration in simulation studies [1]. In order to study the effect of the vertical permeability heterogeneity of the EGS reservoirs on the heat production and electricity generation performance, in this work, we investigated the influence of the permeability heterogeneity on the production performance through a numerical simulation based on the geological data at the Yangbajing geothermal field. The geologic structure, and hydrological and geothermal features of the Yangbajing geothermal field are reported in references [6,36,37], in detail. This will lay a solid foundation for future research and the development of EGS [1].

Research Objectives
The research objectives of this work are to investigate the effect of the vertical permeability heterogeneity of the EGS reservoirs on heat production and electricity generation performance, and to find out the optimum reservoir permeability conditions for the electricity generation of EGS reservoirs. This will provide guidance for the future site selection and performance optimization of EGS.

Physical-Mathematical Model
According to experiences of numerical simulation studies of heat mining, in this work, we used the five-spot well configuration, as shown in Figure 1 [8]; one injection well is located in the center and four production wells are located at the corners. The distance between two adjacent production wells is 1000 m. Because of symmetry, only one fourth of the model domain needs to be simulated, as the grey fraction shown in Figure 1. Based on the geological data at the Yangbajing geothermal field, the deep fractured granite reservoir is located within 950-1350 m, and thickness of the reservoir is 400 m. Assuming, after stimulation, the average fracture spacing is lower than 2-3 m, the heat transfer between the cap rock or base rock and the reservoir can be neglected [6]. Because the average fracture spacing is lower than 2-3 m, the local thermodynamic equilibrium between the fractured rock and circulating water is reasonable [4], thus the fractured granite reservoir can be represented as an equivalent porous medium [2,[4][5][6][7][8][9]. The cold water enters into the reservoir through the perforated interval of 50 m length at the bottom position of the injection well and leaves the reservoir through the perforated interval of 50 m length at the top position of the production well. This kind of production scheme can greatly prolong the reservoir lifetime, reduce the reservoir impedance and improve the production efficiency [2]. In this study, we used the TOUGH2-EOS1 codes to carry out the simulations, which was developed at the Lawrence Berkeley National Laboratory [38]. The TOUGH2 codes employed the integral finite difference method to solve the conservation equations of mass, momentum and energy, and the comparison between the computation through TOUGH2 codes and the experiment has proved that the calculation results are highly reliable [38]. The governing equations of fluid flow and heat transfer based on the equivalent porous medium method can be found in [38].
Processes 2021, 9, x FOR PEER REVIEW 5 of 17   In order to investigate the influence of vertical permeability heterogeneity in stratified formation on the system heat production, we equally divided the 950-1350 m reservoir into 8 layers with different permeabilities; the thickness of every layer is 50 m and the horizontal and vertical permeability in each layer are different. In this study, we neglected the dependence of permeability on the reservoir porosity and assumed that all porosities in each layer are 0.10 [39]. The average horizontal permeability of the 8 layers was kept equal at 50 × 10 −15 m 2 . We assumed that the ratio of the horizontal permeability k h to the vertical permeability k v is constant as k h /k v = α [39]. The degree of reservoir vertical heterogeneity is given by the non-dimensional parameter β as Equation (1) (2). The higher the reservoir permeability heterogeneity, the greater the β; for the homogenous reservoir, β = 0. In this study, we researched, under the three conditions of α = 2, α = 5 and α = 10 and under the four conditions of β = 0 (for homogenous reservoir), β = 0.3, β = 0.5 and β = 0.8, the influence of the permeability heterogeneity on the system electricity generation performance. In order to compare the performance of the heterogeneous reservoir with the homogenous reservoir, we calculated and analyzed the production performance of homogenous reservoirs of β = 0 under the conditions of α = 2, α = 5 and α = 10, respectively. The reservoir permeability under different cases are listed in Table 1. Due to geologic sedimentation and compaction, the horizontal permeability decreases with increasing depth under the same material composition, thus in this study, the horizontal permeability declines with increasing depth, and this more agrees with the actual situation.

Domain, Grid and Parameters
In order to investigate the influence of vertical permeability heterogeneity on the heat production performance of EGS reservoirs, we first researched the production performance of a homogenous EGS reservoir, then we compared the production performance of the heterogeneous reservoirs with that of the homogenous reservoirs. Figure 2 shows the grid system used in this simulation. There are 25 gridblocks in the x direction and in the y direction, and every length of the gridblocks is 20 m. In the z direction, there are 8 gridblocks, and every height of the gridblocks is 50 m. Thus, the three-dimensional grid system in Figure 2 comprises a total of 25 × 25 × 8 = 5000 gridblocks. A mesh independence check was done to make sure the computation results were not influenced by the mesh division. Table 2 shows the properties and conditions of the simulated domain of the 950-1350 m fractured reservoir [2,6]. According to previous studies, the bottomhole injection pressure P inj must be lower than an upper limit P max [2,6]: where P max = P W0 + P b [2,6]; P W0 is the initial pressure of the wellbore and P b is the maximum pressure buildup. In this study, the perforated interval of the injection well is located at −1350-−1300 m and the perforated interval of the production well is located at −1000-−950 m. The intermediate depth of the injection perforated interval is z = −1325 m, corresponding to an initial pressure of P W0 = 11.35 MPa. Based on the engineering data at the Desert Peak geothermal field, P b can be assigned as 6.90 MPa [5]. So P max = (11.35 + 6.90) MPa = 18.25 MPa. The intermediate depth of the production perforated interval is z = −975 m, corresponding to an initial pressure of 8.23 MPa. According to the engineering data at the Desert Peak geothermal field, the maximum pressure drawdown is 3.40 MPa, thus the minimum pressure at the bottomhole of the production well is about (8.23 − 3.40) MPa = 4.83 MPa. In this study, we assume that the production well is kept at P 0 = 5.0 MPa to maintain the heat production. Based on (3), the maximum water production rate can be determined. The injection temperature is T inj = 60 • C [4,5]; for an initial wellbore pressure of P W0 = 11.35 MPa, the injection specific enthalpy is about h inj = 260.66 kJ/kg.

Domain, Grid and Parameters
In order to investigate the influence of vertical permeability heterogeneity on the heat production performance of EGS reservoirs, we first researched the production performance of a homogenous EGS reservoir, then we compared the production performance of the heterogeneous reservoirs with that of the homogenous reservoirs. Figure 2 shows the grid system used in this simulation. There are 25 gridblocks in the x direction and in the y direction, and every length of the gridblocks is 20 m. In the z direction, there are 8 gridblocks, and every height of the gridblocks is 50 m. Thus, the three-dimensional grid system in Figure 2 comprises a total of 25 × 25 × 8 = 5000 gridblocks. A mesh independence check was done to make sure the computation results were not influenced by the mesh division. Table 2 shows the properties and conditions of the simulated domain of the 950-1350 m fractured reservoir [2,6]. According to previous studies, the bottomhole injection pressure Pinj must be lower than an upper limit Pmax [2,6]: where Pmax = PW0 + Pb [2,6]; PW0 is the initial pressure of the wellbore and Pb is the maximum pressure buildup. In this study, the perforated interval of the injection well is located at In this study, we assume that the production well is kept at P0 = 5.0 MPa to maintain the heat production. Based on (3), the maximum water production rate can be determined. The injection temperature is Tinj = 60 °C [4,5]; for an initial wellbore pressure of PW0 = 11.35 MPa, the injection specific enthalpy is about hinj = 260.66 kJ/kg.

Boundary and Initial Conditions
The saturated vapor pressure corresponding to 248 • C is 3.84 MPa, far lower than the bottomhole pressure at the production well of 5.00 MPa, so the water in the reservoir and wells is all maintained as a liquid. The initial reservoir pressure is P = −0.0089z − 0.4444 (MPa), and the initial reservoir temperature is 248 • C [2,6]. Due to symmetry, the lateral boundaries in Figure 2 are no-flow for mass and heat. As stated above, because the slight conductive heat exchange between the permeable reservoir and the impermeable cap rock or base rock is neglected, the boundaries at the top and bottom of the reservoir in Figure 2 are also no-flow for mass and heat. The relative convergence error is assigned as 1.0 × 10 −5 and the absolute convergence error is assigned as 1.0 for every conservation equation. Table 3 shows the maximum water production rate Q = 4q under the various conditions of permeability distribution shown in Table 1. According to (3), P inj will increase with increasing q, and the maximum q can be obtained only when the P inj is closing to the P max . When the α = 2, comparing case 1 with case 2, case 3 and case 4, we can find the water production rate q decreases with the increasing of the degree of reservoir vertical heterogeneity β; under the condition of the homogenous reservoir of β = 0, the system obtains a maximum water production rate. When the α = 5 or α = 10, comparing case 5 with case 6, case 7 and case 8, or comparing case 9 with case 10, case 11 and case 12, we can find that only under the condition of the homogenous reservoir of β = 0, the system obtains a maximum water production rate; with the increasing of the degree of reservoir vertical heterogeneity, the β increases and the water production rate q declines. When the β = 0, comparing case 1 with case 5 and case 9, we can find that the vertical permeability k v has an important influence on the water production rate q, and higher vertical permeability will obviously increase the q. Similarly, comparing the group of β = 0.3, the group of β = 0.5 or the group of β = 0.8 in Table 1, we can find similar conclusions. Overall, when the average value of the horizontal permeability of the layered reservoir is constant, the water production rate q will decrease with the increasing of the degree of reservoir vertical heterogeneity β. Table 3. The water production rate Q = 4q under various conditions of permeability distribution.  Figure 3 shows the evolution of the production temperature under various conditions of permeability distribution, shown in Table 1, over 50 years of heat production [2,6]. When α = 2, comparing case 1 with case 2, case 3 and case 4, we can easily find that with the increasing of the degree of reservoir vertical heterogeneity β, the duration of the stable stage gradually increases, the duration of the declining stage gradually decreases, and the production temperature during the declining stage slowly increases. When α = 5 or α = 10, comparing case 5 with case 6, case 7 and case 8, or comparing case 9 with case 10, case 11 and case 12, we can obtain similar findings. When β = 0, comparing case 1 and case 5, case 9, we can find that with the increasing of α, the duration of the stable stage slowly increases, the duration of the declining stage slowly decreases, and the production temperature during the declining stage gradually increases. When the β = 0.3, β = 0.5 or β = 0.8, through comparison we can obtain similar findings. Overall, when the average value of the horizontal permeability of the layered reservoir is constant, with the increasing of the degree of reservoir vertical heterogeneity β, the duration of the declining stage gradually decreases, and the production temperature during the declining stage slowly increases. This is in good agreement with the fact that the q decreases with the increasing of the degree of reservoir vertical heterogeneity β and a lower q obtains a higher production temperature during the declining stage.

Production Temperature
Processes 2021, 9, x FOR PEER REVIEW 8 Figure 3 shows the evolution of the production temperature under various c tions of permeability distribution, shown in Table 1, over 50 years of heat production When α = 2, comparing case 1 with case 2, case 3 and case 4, we can easily find that the increasing of the degree of reservoir vertical heterogeneity β, the duration of the s stage gradually increases, the duration of the declining stage gradually decreases, an production temperature during the declining stage slowly increases. When α = 5 o 10, comparing case 5 with case 6, case 7 and case 8, or comparing case 9 with case 10 11 and case 12, we can obtain similar findings. When β = 0, comparing case 1 and c case 9, we can find that with the increasing of α, the duration of the stable stage s increases, the duration of the declining stage slowly decreases, and the production perature during the declining stage gradually increases. When the β = 0.3, β = 0.5 o 0.8, through comparison we can obtain similar findings. Overall, when the average of the horizontal permeability of the layered reservoir is constant, with the increasi the degree of reservoir vertical heterogeneity β, the duration of the declining stage ually decreases, and the production temperature during the declining stage slow creases. This is in good agreement with the fact that the q decreases with the increas the degree of reservoir vertical heterogeneity β and a lower q obtains a higher produ temperature during the declining stage.

Electric Power
The electric power We is calculated as Equation (4) in reference [4], where the inje specific enthalpy inj   Table 1 over 50 years of heat production. When α = 2, comparing 1 with case 2, case 3 and case 4, we can easily find that with the increasing of the d

Electric Power
The electric power W e is calculated as Equation (4) in reference [4], where the injection specific enthalpy h inj is h inj = 260.66 kJ/kg. The Yangbajing geothermal power plant comprises a single-flash power plant and double-flash power plant, which are open loop systems. The heat rejection temperature of the Yangbajing geothermal power plant is T 0 = 9 • C = 282.15 K [2,6]. Figure 4 shows the evolution of electric power under various conditions of permeability distribution in Table 1 over 50 years of heat production. When α = 2, comparing case 1 with case 2, case 3 and case 4, we can easily find that with the increasing of the degree of reservoir vertical heterogeneity β, the electric power gradually decreases. Higher β obviously decreases the q, which prolongs the duration of the stable stage, decreases the duration of the declining stage, and, meanwhile, obviously declines the electric power over the declining stage. When α = 5 or α = 10, we can obtain similar findings through comparisons. When β = 0, comparing case 1 with case 5, case 9, we can find that with the increasing of α, the duration of the stable stage gradually increases, and the electric power during the stable stage also gradually decreases; the duration of the declining stage slowly decreases, and the electric power over the declining stage also slowly declines. When β = 0.3, or β = 0.5 or β = 0.8, we can obtain similar findings through comparisons. This is mainly because the reservoir vertical heterogeneity obviously influences the water production rate, and according to Equation (4), this will significantly affect the electric power. Overall, when the average value of the horizontal permeability of the layered reservoir is constant, for a homogenous reservoir, the system will obtain maximum electric power; with the increasing of the reservoir vertical heterogeneity, the electric power will gradually decline.
Processes 2021, 9, x FOR PEER REVIEW 9 of reservoir vertical heterogeneity β, the electric power gradually decreases. Higher viously decreases the q, which prolongs the duration of the stable stage, decreas duration of the declining stage, and, meanwhile, obviously declines the electric p over the declining stage. When α = 5 or α = 10, we can obtain similar findings thr comparisons. When β = 0, comparing case 1 with case 5, case 9, we can find that wi increasing of α, the duration of the stable stage gradually increases, and the electric p during the stable stage also gradually decreases; the duration of the declining stage s decreases, and the electric power over the declining stage also slowly declines. Wh 0.3, or β = 0.5 or β = 0.8, we can obtain similar findings through comparisons. This is m because the reservoir vertical heterogeneity obviously influences the water produ rate, and according to Equation (4), this will significantly affect the electric power. Ov when the average value of the horizontal permeability of the layered reservoir is con for a homogenous reservoir, the system will obtain maximum electric power; wi increasing of the reservoir vertical heterogeneity, the electric power will gradually de  Figure 5 shows the evolution of the injection pressure under the various cond of permeability distribution, shown in Table 1, over 50 years of heat production. W = 2, for β = 0, β = 0.3, β = 0.5 and β = 0.8, the injection pressure is 11.85-17.11 MPa, 1 17.92 MPa, 11.78-18.30 MPa and 11.53-18.43 MPa, respectively. So, with the increas the reservoir vertical heterogeneity β, the Pinj has only a very slight change. When α α = 10, through comparisons, we can easily obtain similar findings. When β = 0, for (case 1), α = 5 (case 5) and α = 10 (case 9), the corresponding Pinj is 11.85-17.11 MPa, 1 17.90 MPa and 11.76-17.78 MPa, respectively. So, with the increasing of α, the Pinj ha a very slight change. For β = 0.3, β = 0.5 and β = 0.8, through comparisons, we can obtain similar findings. Overall, the reservoir vertical heterogeneity β has only a slight influence on the Pinj and it is not the main factor affecting the injection pre Figure 6 also displays the Pinj slowly increasing during heat production (see Section This is because the reservoir impedance is gradually increasing; when the produ  Figure 5 shows the evolution of the injection pressure under the various conditions of permeability distribution, shown in Table 1, over 50 years of heat production. When α = 2, for β = 0, β = 0.3, β = 0.5 and β = 0.8, the injection pressure is 11.85-17.11 MPa, 11.85-17.92 MPa, 11.78-18.30 MPa and 11.53-18.43 MPa, respectively. So, with the increasing of the reservoir vertical heterogeneity β, the P inj has only a very slight change. When α = 5 or α = 10, through comparisons, we can easily obtain similar findings. When β = 0, for α = 2 (case 1), α = 5 (case 5) and α = 10 (case 9), the corresponding P inj is 11.85-17.11 MPa, 11.82-17.90 MPa and 11.76-17.78 MPa, respectively. So, with the increasing of α, the P inj has also a very slight change. For β = 0.3, β = 0.5 and β = 0.8, through comparisons, we can easily obtain similar findings. Overall, the reservoir vertical heterogeneity β has only a very slight influence on the P inj and it is not the main factor affecting the injection pressure. Figure 6 also displays the P inj slowly increasing during heat production (see Section 3.5). This is because the reservoir impedance is gradually increasing; when the production pressure is constant, based on Equation (5), the P inj will gradually increase.

Reservoir Impedance
The reservoir impedance R I is calculated as Equation (5) [4,5]: where R I represents the power consumption of the unit production rate for penetra the fractured reservoir, the pro P = P0 = 5.00 MPa is the bottomhole production pres and q is the water production rate between the injection well and the production Figure 6 displays the evolution of the reservoir impedance under various condition permeability distribution, shown in Table 1 I . Figure 6 also shows that during heat production, the R I gradu increases. This is because with the heat production, the temperature of the reservoir and fluid gradually declines, causing the slow increase in water viscosity, and this fin

Pump Power
The internal energy consumption p W of both the injection pump and produ pump can be calculated as Equation (6) [4]: h is the depth of the injection well and 2 h is the depth of the production

Reservoir Impedance
The reservoir impedance I R is calculated as Equation (5) [4,5]: where I R represents the power consumption of the unit production rate for penetrating the fractured reservoir, the P pro = P 0 = 5.00 MPa is the bottomhole production pressure, and q is the water production rate between the injection well and the production well. Figure 6 displays the evolution of the reservoir impedance under various conditions of permeability distribution, shown in Table 1 respectively. So, with the increasing of α, the system I R gradually increases. For β = 0.3, β = 0.5 and β = 0.8, through comparisons, we can easily get similar findings. As stated in Section 3.1, this is because a higher reservoir heterogeneity obviously decreases the water production rate. When the P inj basically maintains at a constant, according to Equation (5), this will obviously increase the I R . Overall, when the average value of the horizontal permeability of the layered reservoir is constant, for the homogenous reservoir, the system will obtain a minimum I R ; with the increasing of reservoir vertical heterogeneity β, the I R gradually increases. This proves that the stimulation of an EGS reservoir should be controlled as homogeneously as possible to obviously reduce the I R . Figure 6 also shows that during heat production, the I R gradually increases. This is because with the heat production, the temperature of the reservoir rock and fluid gradually declines, causing the slow increase in water viscosity, and this finally causes increasing of the I R . This has been reported in detail in previous studies from Zeng et al. [2,4,6].

Pump Power
The internal energy consumption W p of both the injection pump and production pump can be calculated as Equation (6) [4]: where h 1 is the depth of the injection well and h 2 is the depth of the production well. In this work, h 1 = 1325.0 m, h 2 = 975.0 m and h 1 − h 2 = 350.0 m. The water density ρ changes versus the temperature and pressure, and this will obviously influence the calculation of W p [4,5]. In this study, the average values of the maximum density and the minimum density are adopted for calculations in Equations (6) and (7), and this treatment method has proved to be effective in previous studies [4,5]. The maximum density is 998.0 kg/m 3 and the minimum density is 804.0 kg/m 3 , thus the average value of ρ = 901.0 kg/m 3 is adopted for calculations in Equations (6) and (7). Figure 7 displays the evolution of the pump power under the various conditions of permeability distribution in Table 1 over 50 years of heat production. When α = 2, for β = 0, β = 0.3, β = 0.5 and β = 0.8, the corresponding W p is 0.52-1.25 MW, 0.42-1.09 MW, 0.31-0.85 MW and 0.09-0.26 MW, respectively. It can be found that for the homogenous reservoir (β = 0), the system obtains maximum W p ; with the increasing of reservoir vertical heterogeneity β, the W p gradually declines. This is mainly because with the increasing of reservoir vertical heterogeneity β, the q gradually decreases, according to Equation (6), and the W p gradually declines. For α = 5 or α = 10, we can get similar findings through comparisons. When β = 0, for α = 2, α = 5 and α = 10, the corresponding W p is 0.52-1.25 MW, 0.48-1.25 MW and 0.41-1.08 MW, respectively. So, with the increasing of α, the q gradually declines and the W p gradually decreases. For β = 0.3, β = 0.5 and β = 0.8, we can obtain similar findings through comparisons.
Overall, the reservoir vertical heterogeneity has a significant influence on the W p when the average value of the horizontal permeability of the layered reservoir is constant and the W p gradually declines with the increasing of the reservoir vertical heterogeneity; for a homogenous reservoir, the system obtains a maximum W p . Figure 7 also shows that the W p gradually increases with heat production, and this is because the P inj gradually increases with time; these are in good agreement with previous studies from Zeng et al. [2,[4][5][6].

Energy Efficiency
The energy efficiency η of the system can be calculated as Equation (7) Table 1 ov years of heat production. When α = 2, for β = 0, β = 0.3, β = 0.5 and β = 0.8, the corresp ing η is 32.29-10.70, 32.27-11.36, 32.83-11.65 and 35.27-11.72, respectively. It can be fo that the reservoir vertical heterogeneity β has only a very slight influence on the η. W α = 5 or α = 10, we can get similar findings through comparisons. When β = 0, for α = = 5 and α = 10, the corresponding η is 32.29-10.70, 32.53-10.44 and 33.04-11.41, re tively; it can be found that the α also has only a very slight influence on the η. Whe 0.3, β = 0.5 or β = 0.8, we can get similar findings through comparisons. Overall, whe average value of the horizontal permeability of the layered reservoir is constant, the ervoir vertical heterogeneity has only a very slight influence on the η. According to E tion (7), this is because the reservoir vertical heterogeneity has only a very slight influ on the Pinj as stated in Section 3.4; meanwhile, it has also a very slight influence on the , thus the reservoir vertical heterogeneity has only very slight influence on the η. Fig  also shows that the η will decline during the heat production. This is because the Pin gradually increase and pro h and Tpro will gradually decline; this is in good accord with previous studies from Zeng et al. [2,4,5].

Energy Efficiency
The energy efficiency η of the system can be calculated as Equation (7) based on [4,5]: As stated above, the average value of the water density of ρ = 901.0 kg/m 3 is adopted in the calculation of Equation (7). Figure 8 displays the evolution of the energy efficiency under the various conditions of permeability distribution in Table 1 over 50 years of heat production. When α = 2, for β = 0, β = 0.3, β = 0.5 and β = 0.8, the corresponding η is 32.29-10.70, 32.27-11.36, 32.83-11.65 and 35.27-11.72, respectively. It can be found that the reservoir vertical heterogeneity β has only a very slight influence on the η. When α = 5 or α = 10, we can get similar findings through comparisons. When β = 0, for α = 2, α = 5 and α = 10, the corresponding η is 32.29-10.70, 32.53-10.44 and 33.04-11.41, respectively; it can be found that the α also has only a very slight influence on the η. When β = 0.3, β = 0.5 or β = 0.8, we can get similar findings through comparisons. Overall, when the average value of the horizontal permeability of the layered reservoir is constant, the reservoir vertical heterogeneity has only a very slight influence on the η. According to Equation (7), this is because the reservoir vertical heterogeneity has only a very slight influence on the P inj as stated in Section 3.4; meanwhile, it has also a very slight influence on the h pro , thus the reservoir vertical heterogeneity has only very slight influence on the η. Figure 8 also shows that the η will decline during the heat production. This is because the P inj will gradually increase and h pro and T pro will gradually decline; this is in good accordance with previous studies from Zeng et al. [2,4,5].

Model Validation
Strictly, the numerical results should be compared with the experimental da make sure that the computation is accurate and reliable. However, because the ex mental modeling of geothermal production is complex and time-consuming, until experimental data about heat production at the Yangbajing geothermal field have scarce. For simplification, we compared the results of Case 1 in this work with a prev study in reference [8], where most of the model conditions and properties are basi the same. Figure 9 shows the comparison of the production temperatures between C in this work and the basic case in [8]. In Case 1, the production temperature drops 248.0 °C to 221.8 °C over 50 years, while in the basic case in [8], the production tempera drops from 248.0 °C to 232.5 °C over 50 years. The relative error of the production perature between the two cases increases with time, and finally, it reaches 4.8%, whi much smaller than 10%. Because the results of the two simulated cases are very close ative error is smaller than 5%), we can regard the simulations in this work as reliable

Model Validation
Strictly, the numerical results should be compared with the experimental data to make sure that the computation is accurate and reliable. However, because the experimental modeling of geothermal production is complex and time-consuming, until now, experimental data about heat production at the Yangbajing geothermal field have been scarce. For simplification, we compared the results of Case 1 in this work with a previous study in reference [8], where most of the model conditions and properties are basically the same. Figure 9 shows the comparison of the production temperatures between Case 1 in this work and the basic case in [8]. In Case 1, the production temperature drops from 248.0 • C to 221.8 • C over 50 years, while in the basic case in [8], the production temperature drops from 248.0 • C to 232.5 • C over 50 years. The relative error of the production temperature between the two cases increases with time, and finally, it reaches 4.8%, which is much smaller than 10%. Because the results of the two simulated cases are very close (relative error is smaller than 5%), we can regard the simulations in this work as reliable.

Model Validation
Strictly, the numerical results should be compared with the experimental data to make sure that the computation is accurate and reliable. However, because the experimental modeling of geothermal production is complex and time-consuming, until now, experimental data about heat production at the Yangbajing geothermal field have been scarce. For simplification, we compared the results of Case 1 in this work with a previous study in reference [8], where most of the model conditions and properties are basically the same. Figure 9 shows the comparison of the production temperatures between Case 1 in this work and the basic case in [8]. In Case 1, the production temperature drops from 248.0 °C to 221.8 °C over 50 years, while in the basic case in [8], the production temperature drops from 248.0 °C to 232.5 °C over 50 years. The relative error of the production temperature between the two cases increases with time, and finally, it reaches 4.8%, which is much smaller than 10%. Because the results of the two simulated cases are very close (relative error is smaller than 5%), we can regard the simulations in this work as reliable.

Conclusions
In this work, we numerically investigated the effect of vertical permeability heterogeneity on the electricity generation performance of EGS based on the geological data at the Yangbajing geothermal field. Based on the study, the following conclusions are made.
(1) When the average value of the horizontal permeability of the layered reservoir is constant, for a homogenous reservoir, the system obtains a maximum water production rate. (2) When the average value of the horizontal permeability of the layered reservoir is constant, with the increasing of the reservoir permeability heterogeneity, the duration of the declining stage gradually declines and the production temperature during the declining stage gradually increases. (3) When the average value of the horizontal permeability of the layered reservoir is constant, for a homogenous reservoir, the system obtains a maximum electric power. (4) When the average value of the horizontal permeability of the layered reservoir is constant, for a homogenous reservoir, the system obtains a minimum reservoir impedance. (5) When the average value of the horizontal permeability of the layered reservoir is constant, for a homogenous reservoir, the system obtains a maximum pump power. (6) When the average value of the horizontal permeability of the layered reservoir is constant, with the increasing of the vertical permeability heterogeneity, the water production rate gradually decreases, the electric power gradually declines, the reservoir impedance gradually increases, and the pump power gradually declines. (7) When the average value of the horizontal permeability of the layered reservoir is constant, with the increasing of the reservoir permeability heterogeneity, the injection pressure and energy efficiency have only very slight changes. This proves that the reservoir permeability heterogeneity is not the main factor affecting the injection pressure and energy efficiency. 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 h pro 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 P inj injection pressure, MPa P pro 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 T pro production temperature,°C W p 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