Predictions of Rock Temperature Evolution at the Lahendong Geothermal Field by Coupled Numerical Model with Discrete Fracture Model Scheme

: The comprehensive exploitation of geothermal fields has an impact on the productivity of the reservoirs. To realize sustainable steam production, changes in the rock temperature need to be predicted and controlled. A coupled thermo ‐ hydro ‐ mechanical (THM) model employing COMSOL Multiphysics was proposed to study the characteristics of heat transfer, fluid flow, and solid deformation at the Lahendong geothermal field, in North Sulawesi, Indonesia. The numerical results were compared with analytical and measured data in order to validate the numerical simulation. Based on the results, the predicted temperatures of the production wells showed significant decrease with the production time. In addition, a reduction in the reservoir temperature leads to lower specific gross electrical power within the production well, which should significantly reduce the sustainability of the power plant.


Introduction
The growing use of renewable energy sources around the world is the effect of increased efforts to contribute to a decrease in the dependence on gas and other fossil fuels. As a country located in the ring of fire, Indonesia has many active volcanoes with enormous geothermal potential [1]. With a total installed capacity of 1948 MWe, as of the end of 2018, Indonesia has become the second-largest nation in the world in terms of employing geothermal energy for electricity [2]. Although geothermal energy is renewable energy, the comprehensive exploitation of geothermal fields has an impact on the geothermal reservoir evolution. A large amount of geothermal exploitation will cause dryness in the reservoirs, pressure drawdown, thermal drawdown, silica scaling, and other issues [3,4]. To realize sustainable steam production, changes in the rock temperature at the reservoirs need to be predicted and controlled. Therefore, further field monitoring, laboratory experiments, theoretical research, and numerical simulations at geothermal fields are still needed.
In this work, a three-dimensional (3D) numerical model was developed to study the characteristics of heat transfer, fluid flow, and deformation at a geothermal field by using a multiphysics FEM solver, COMSOL Multiphysics (COMSOL Multiphysics ® v. 5.5, 2019). The Lahendong geothermal field, located in Tomohon, North Sulawesi, Indonesia (Figure 1), was examined for the present case study. The Lahendong geothermal field has four units power plant with a total production capacity of 80 MWe [5]. Unit 1, the first 20-MWe power plant, has been producing electricity since June 2001; Unit 2, the second 20-MWe power plant, has been producing electricity since 2007; and Units 3 and 4 (2 × 20 MWe) have been producing electricity since 2009. It is important to understand the evolution of the reservoir characteristics in order to evaluate the sustainability of the reservoir. The reservoir temperature is divided into two different reservoir categories. There are very high temperatures in the southern zone (350 °C) and modest temperatures in the northern zone (250 °C). The southern and northern zones show dryness ranging from about 80% to 100% and 30% to 50%, respectively [6].  [5] and Koestono [7]).
To date, several numerical simulations have been conducted to simulate the Lahendong geothermal field. Yani [8] applied TOUGH2 to simulate the natural state condition of the Lahendong geothermal reservoir and validated the predictions using measured pressure and temperature data. He also predicted the changes in reservoir pressure and temperature after production had been started at Unit 1, including the scenario of increased production with the addition of Units 2 and 3 during the 30-year operation period. However, the scenario of production have not yet been validated because of the lack of measured pressure drawdown data [8]. In 2015, Sumantoro et al. [9] improved Yani [8] reservoir modeling by adjusting the permeability, enthalpy, and flow rate parameters. Moreover, they added the actual topography, extended the boundary conditions, and used the airwater equation of state to model the shallow unsaturated zone. Brehme et al. [10] conducted a thermal-hydraulic simulation to determine the permeability of the Lahendong geothermal reservoir by using the finite element software FEFLOW. Another study carried out an inverse modeling technique to achieve the optimal model parameters and simulated the natural state condition in the Lahendong geothermal field [11]. Nevertheless, the above studies only focused on the simulation of the initial reservoir condition in the Lahendong geothermal field. In addition, predictions on future production have not yet been conducted. Therefore, accurate and reliable predictions of the initial conditions and the future production at this location are still required.
The most important task in many rock engineering projects is to model the complex interaction among the fluid migration, heat transfer, rock mass deformation, species transport, and chemical reactions [12][13][14][15]. However, heat transfer and fluid migration are the most important processes for estimating the economic efficiency and production period in geothermal reservoirs [16]. In this study, the heat transfer, fluid migration, and rock mass deformation processes have been considered; they are generally referred to as thermal, hydraulic, and mechanical (THM) coupling.
THM coupling has been widely used to simulate geothermal systems especially in enhanced geothermal systems (EGSs). The number of authors has indicated the ability of EGS as a possible means of providing renewable energy [17]. Thermoporoelastic deformation of the rock matrix has been extensively studied with the THM coupling simulations for single fracture [18,19]. Together with increased fluid pressure in the fracture, heat contraction of the rock matrix reduces contact stress on the fracture surfaces and enhances the opening of fracture. The effect of thermoelastic between multiple fractures by using THM coupling has confirmed the enhancement of the flow locally and reduces the efficiency of the system in the EGS [20]. In order to better understand the efficiency of EGS, Zhao et al. [21] simulated a 3D THM coupled model of geothermal energy extraction plan at the Tengchong geothermal field, in China. Numerical results showed an exponential increase in the temperature with the distance between well injection and well production. Sun et al. [22] utilized a two-dimensional (2D) THM coupling to study the heat extraction process in EGS at the Cooper Basin, Australia. Yao et al. [23] and Sun et al. [24] proposed a 3D THM coupling to evaluate the EGS heat production performance. Nevertheless, THM coupling has not often been adopted for examining hydrothermal geothermal resources.
In hydrothermal geothermal resources, a natural state simulation of the reservoir is required before simulating the actual production. The predicted results should be in good agreement with the measured data, such as the pressure and temperature of the wells. Magnenet et al. [25] developed a two-dimensional (2D) THM model to predict the natural hydrothermal circulation at Soultz-sous-Forêts geothermal field. However, no faults or fracture networks are explicitly considered in the study. Alternatively, by considering the large-scale faults, the coupled THM model was used to simulate the hydrothermal circulation at Rittershoffen geothermal field [26]. Nonetheless, these approaches were not readily applicable to predict production performance.
The main purpose of this study was to better understand the natural state and the production evolution of a geothermal reservoir in a hydrothermal geothermal resource. First, a steady state fluid migration, heat transfer, and rock mass deformation were simulated. Subsequently, the THM model was analyzed comprehensively to understand the production evolution of geothermal reservoir. A coupled THM model was developed with a discrete fracture model (DFM). The reservoir was divided by the DFM into the rock matrix and discrete fractures. For each fracture, flexible Delaunay triangulation was used to simulate the complex fracture network [27,28]. Finally, the developed 3D model was applied to the Lahendong geothermal field and then validated by comparing the predictions with the measured data.

Conceptual Model
The conceptual model was developed based on analyses of geological, geophysical, geochemical, and well-log data, which are necessary for developing an appropriate numerical model for the target geothermal field. Several conceptual models have been proposed for the Lahendong geothermal field [7,29,30]. In this chapter, some interpolation and interpretation data for making the conceptual model are presented. Figure 2 shows the conceptual model of the Lahendong geothermal field reported by Brehme et al. [5].

Geological
Geological mapping of the Lahendong field was started by P.T. Pertamina Geothermal Energy in 1982 and adopted by many researchers [7,8,30,32,33]. Mount Lengkoan and Mount Tampusu are two volcanoes dominating this location, and a crater lake, Lake Linau, is located in the center of the geothermal field [5]. A geological and topographic map of the Lahendong geothermal field is shown in Figure 3.

Geophysical
Gravity, resistivity, magnetotelluric (MT), and Mise-á-la-masse methods have been conducted in the Lahendong geothermal field [7,33,35]. Initially, the focus of the geophysical investigations was to find the center of the study area, where the up-flow zone is located. Regarding the gravity investigation, the up-flow zone of the study area was in Lake Linau [33]. Raharjo et al. [35] reported a 3D MT inversion and showed an up-dome resistivity structure around Lake Linau. This up-dome resistivity structure was characterized by a resistivity value ranging from 20-60 Ωm and covered by a conductor resistivity structure with a value of <10 Ωm as the clay cap ( Figure 4).

Well-Log Analysis
As reported by Sumintadireja et al. [33], the Lahendong geothermal field has the following formations found at the following depths. The Post Tondano/Pangolombian formation is located at a depth of 0 m-850 m, the Tondano formation is located at a depth of 350 m-1100 m, and the Pre-Tondano formation is located at a depth of 1100 m-1600 m. Andesite, basalt, tuff, volcanic breccia, and pyroclasts are the main rock types in this location. Based on the measured data, the pore pressure in this area is close to hydrostatic pressure. This means that the reservoir is fluid-dominated and that the fluid flows through the southern zone to the northern zone in this area [33]. The temperature of the reservoir is controlled by the convection regime from that fluid flow.

Methodology
The following assumptions are made to formulate the mathematical model in this study: 1. The rock mass is treated as a 3D fractured porous media consisting of the rock matrix and discrete fractures. The discrete fracture model is applied to the fault zone in this study area. 2. The fractured porous geothermal reservoir has a single-phase fluid saturation. Therefore, the water and fluid flow both in the rock matrix block and fractures comply with the Darcy's Law. 3. The model ignores the variation of fracture aperture. 4. The diameter of the well is small, so that storage is negligible. 5. As shown in the published paper [22], the water density and dynamic viscosity are not constant but a function of pressure and temperature. 6. The density, porosity, permeability, specific heat, and thermal conductivity of the fractured porous media are assumed to be constant.

Discrete Fracture Model
In the discrete fracture model (DFM), the rock properties are divided into the rock matrix and discrete fractures [27,36]. The matrix to matrix, fracture to fracture, and matrix to fracture relationships are represented by grid nodes. Figure 5 shows the physical domain and the grid domain in the fractured porous medium. The physical domain represents the fractured porous medium in the real condition, and the discretization from that physical domain is the grid domain. Although the fracture thickness is not denoted in the grid domain, it is considered in the computational domain for the flow-rate estimation. These simplified grid domains are useful for calculating the fractures implicitly. Furthermore, elimination of fracture intersections could improve the simulator's stability over a longer time-step and improve the simulator's calculation capacity [27,37].

Governing Equations
The governing equations of the model describe the coupled THM processes among fluid migration, heat transfer, and rock mass deformation. The differential equations are solved by COMSOL Multiphysics (COMSOL Multiphysics ® v. 5.5, 2019). The differential equations used in this study are described in [23,24]; they are divided into the rock matrix and discrete fractures as explained below.

Fluid Migration
The permeability in the rock matrix is less than that in the discrete fractures. Consequently, the seepage in the rock matrix follows Darcy's law, whereas the seepage in the fractures is calculated by the modified Darcy's law. The equations for the rock matrix and the discrete fracture flow can be expressed as

Rock Mass Temperature
The rock mass temperature in the rock matrix can be calculated by where ρ is the rock density [kg/m 3 ], Cs is the specific heat capacity of the rock matrix [J/kg/K], Ts is the rock temperature [°C], λ is the thermal conductivity of the rock matrix [W/m/K], and q is the heat source [W/m 3 ]. Energy conservation in the discrete fractures is defined as where In the injection well, the diameter of the well is small and can be represented by a line. Therefore, the line heat source can be calculated by where M0 is the mass flow rate [kg/s], l is the well length [m], and Tin is the injection temperature [°C].

Rock Mass Stress Field
The equilibrium equation, used to evaluate the mechanical process in the rock structure, is given as where σij,j is the Cauchy stress tensor and Fi is the body force per unit volume in the x, y, z coordinate in 3D. The deformation equation of the rock matrix is expressed as where E is the elastic modulus [Pa], v is the Poisson's ratio [-], u is the displacement [mm], αB is the Biot-Willis coefficient [-], and αT is the thermal expansion coefficient The deformation equations for the fractures are given as where uf is the displacement of the fractures [mm], σ' is the effective stress of the fractures [Pa], and k is the fracture stiffness [Pa/m]. Subscripts n and s refer to the normal and tangential directions, respectively. The governing equations for calculating the natural state condition are in a stationary condition. This means that the reservoir variable does not change over time to reach the steady flow and pressure field.

Validation of the TH and THM Coupling Model
The purposed model is validated by analytical solution of TH and THM coupling model. A twodimensional (2D) single-fracture model is studied to examine the accuracy of the proposed model using COMSOL Multiphysics. A single fracture is modeled to represent the heat transfer in the complex fracture networks in the reservoir. Furthermore, an analytical solution created by Lauwerier [38] is used to validate the numerical model. Lauwerier [38] solves the heat transfer in saturated porous media by assuming that the heat transfer in the discrete fractures (x direction) is transported only by convection. Meanwhile, the heat transfer in the rock matrix is restrained in the y direction. Based on the Lauwerier problem, a schematic of the 2D single-fracture model is shown in Figure 6. The analytical solution for calculating the temperature from the heat transfer in a single fracture can be written as follows: where T0 is the initial temperature [°C], Tin is the injection temperature [°C], erfc is the complementary error function, and U is a unit step function.  Table 1 presents a list of the material parameters used in the 2D single-fracture model. In the numerical simulation, the geometry is modeled in 100 m × 100 m with the aperture thickness (2df) of 1 mm. Figure 7 shows the distribution of temperature contours at different times. It shows that the temperature change in the fracture decreases faster than that in the rock matrix. Figure 8 provides more detailed information on the temperature distribution along the fracture. Figure 8a shows the temperature change in the fracture at different times (t = 10 days, 100 days, and 1000 days), while Figure 8b shows the temperature at three different positions in the fracture (x = 10 m, 50 m, and 100 m), respectively. It can be seen that the numerical results are in excellent agreement with the analytical solutions.
A THM coupling analysis is validated using the thermal consolidation model proposed by Bai [39]. Figure 9 illustrates an example of the soil thermal consolidation problem. The size of the model is 1.0 m × 0.2 m, with the initial pore pressure and temperature that are 0.1 MPa and 10 °C, respectively. In addition to an external pressure of 0.1 MPa, 0 MPa pore pressure and 60 °C temperature are applied on the top boundary. The lateral and bottom boundaries are considered impermeable and insulated where displacements are normally restricted. The material parameters used in the thermal consolidation problem are presented in Table 2, and the analytical solutions are obtained from the literature [39]. Figures 10 and 11 show the comparison results between the analytical and numerical solutions. It can be seen that the numerical solutions of temperature, pore pressure, and displacement agree well with the analytical solutions. Therefore, the numerical model is proved to be accurate for a THM coupling process.

Natural State Calibration
The natural state or the pre-production calibration was studied to determine the initial condition of the reservoir before starting the production. The natural state condition is obtained when the model is run until it reaches a steady result. In order to reach the results close to the measured data, the specific area, location, and total amounts of heat source and flow rate were adjusted based on the conceptual model and a trial-and-error process. Furthermore, the initially measured pressure and temperature data were used to calibrate and validate the numerical results.

Computational Model
The geometry of the research area is displayed in Figures 12 and 13. The study area covers 35 km 2 with a domain of 5000 m × 7000 m × 2800 m. Three geological layers are set, and the fault zone is applied as the discrete fracture model. The finite element mesh system used in the calculation is illustrated in Figure 14. The mesh consists of 416,850 tetrahedral elements and 65,195 triangular elements. A finer mesh is generated for all elements in the model.

Initial and Boundary Conditions
For the initial temperature on the surface, constant temperatures of 28 °C was considered. There is a cooling pluton beneath Lake Linau with temperatures of 350 °C [10]. The recharge water source was recharged by precipitation [40]. Table 3 shows the boundary conditions of the model. The temperature gradient and hydraulic head were assumed based on the measured data of the wellbores that were located closely to the boundary. Whereas, the heat flux parameter was defined based on the numerical study published by Brehme et al. [10]. The original topography on top of the model was also considered based on the digital elevation model (DEM). The DEM data were collected from a 1 Arc-second DEM United States geological survey (USGS) map. For the evaluation of the reservoir production temperature, five production wells in Cluster-4 (LHD-8, LHD-10, LHD-11, LHD-12, and LHD-15) and two production wells in Cluster-13 (LHD-17 and LHD-18) were examined. These production wells have been providing steam for the turbine, namely, Unit 1 since June 2001 and Unit 2 since June 2007 [41]. For the injection wells, the injection well in Cluster-7 (LHD-7) was examined. Cold water, at the temperature of 40 °C, was injected in the injection well, and a production rate of 12 kg/s was assumed for each production well. Table 4 shows the parameters used in the numerical simulation. The rock parameters were chosen based on the available data from the literature [8,10,29,42]. However, the parameters that were not available from the literature were determined by referring to the existing numerical simulations and making assumptions based on them [24,43].  3 Brehme et al. [10], 4 Silitonga et al. [42]. Figure 15 shows the results of the rock temperature under the natural state condition. Based on the results shown in the figure, the reservoir temperature in the southern zone is higher than that in the northern zone. Therefore, the wells in the southern zone became locations for the production wells in this geothermal system. In the northern zone, the brine and condensation from the production wells were reinjected into the reservoir, as explained by Prabowo et al. [6]. It was confirmed that the current model can replicate the natural state conditions obtained from the literature [8,10]. Comparisons between the numerical simulation results and the measured data are shown in Figures 16 and 17. Six of the horizontal wells (LHD-1, LHD-3, LHD-4, LHD-5, LHD-7, and LHD-13) were explained as being representative of the other wells. This is because horizontal or straight wells are easier to define than deviated wells, that cannot be well-identified for their locations. Figure 16 shows the natural state of the pore pressure based on our predictions compared to the actual measurement data and numerical predictions reported by Yani [8]. It can be seen that the initial pore pressure in this area is hydrostatic. Although several wells, such as wells LHD-4, LHD-5, and LHD-13, have a lower pore pressure than the actual measured data, almost all of the predicted pore pressure values match closely with the actual measured data.

Depth [m]
Temperature

Depth [m]
Temperature The results of the temperature predictions under the natural state condition are shown in Figure  17. These results were compared to the actual measured data and the numerical model reported by Yani [8]. As shown in Figure 17, the temperature in well LHD-1 is lower than the measured data. However, our model predictions show a relatively good agreement with the measured data at the bottom of the well. The predicted temperature in well LHD-3 matches the measured temperature data. The temperature in this well is similar to the thermal gradient, because well LHD-3 is located outside of the main geothermal reservoir area.
Regarding the well temperature results for well LHD-4, the temperature began at 28 °C and then rose steadily to reach 329 °C by the well depth of 2300 m bgl. The results for the well temperature in well LHD-4 do not show a similar trend to that of the measured data. However, the temperature at the bottom of this well has a similar value to that of the measured data. For well LHD-5, the results are lower than those of the measured data. This is because cold and hot water from the northeast and southwest zones, respectively, are mixed in this well [8]. The predicted temperature in well LHD-7 is higher than the measured data until the well depth of 1300 m bgl. On the contrary, the temperature becomes constant after the well depth of 1300 m bgl. In well LHD-13, the temperature is relatively similar to the measured data. Based on the natural state simulations, it can be concluded that fluid migration and heat transfer should be important processes for these simulations. Figure 18 shows a quantitative comparison of the pressure and temperature at the bottom of the wells under the natural state condition between the predictions by the current model and the measured data reported by Yani [8]. The figure shows that the simulated values are relatively congruent with the measured values. The results that are located above the line represent higher pressure or higher temperature values in the predictions than those in the measured data.

Depth [m]
Temperature

Evolution of the Reservoir Temperature
The evolution of the reservoir temperature for long-term production (36 years) has been simulated and is shown in Figure 19. Although all the wells have an initial temperature more or less than the initial condition data, the modeled values are still in good agreement. The temperature of injection well LHD-7 began to decline after the cold-water injection was started. Production wells LHD-11 and LHD-17 are representatives of all the production wells in Cluster-4 and Cluster-13, respectively. The reservoir temperature in LHD-11 was seen to decrease from 275 °C to 248 °C over the 36-year period, while the temperature in production well LHD-17 was seen to decrease from 282 °C to 238 °C in 30 years. The time required for the process of the decreasing reservoir temperature is longer for representative well LDH-11 because almost all of the well positions in Cluster-4 are farther from the injection well. The wellhead pressure data from the seven production wells (LHD-8, LHD-10, LHD-11, LHD-12, LHD-15, LHD-17, and LHD-18) at the Lahendong geothermal field are used for comparison [44]. The data were collected in June 2008 and converted to temperatures using the table of properties for saturated water [45].   Figure 19. Comparison of reservoir temperatures among initial condition, wellhead temperature, and prediction by the current model.
Except well LHD-17, compared to the wellhead production temperature data, that reservoir has temperatures higher than those of the production data. This is because to obtain the maximum mass flow from the reservoir to the power plant, the optimum pressure in the wellhead should be determined.
The existing geothermometer and enthalpy monitoring data were also used to validate the simulation results. A previous study [46] showed the changes in chemical characterization, dryness, and temperature at the Lahendong geothermal reservoir based on the chemical geothermometer and enthalpy monitoring data. These data were used to monitor the geothermal reservoir after it had been exploited for more than 12 years. The enthalpy monitoring data were converted to temperatures using the table of properties for saturated water [46]. A comparison of the changes in reservoir temperature among the geothermometer, enthalpy monitoring data [46], and the results obtained by the simulation is shown in Figure 20. Production well LHD-11 was used as the representative well for this comparison. As shown in Figure 20, the simulation results showed a trend that was close to that of the gas geothermometer, wherein the current model results match the H2S/H2 gas geothermometer calculations well. The temperature of the reservoir fluctuated at around 270 °C over the given period. This means that there was no significant change in temperature in this reservoir. Compared with EGS, high temperature and high pressure of injection well are not considered in the hydrothermal geothermal resources [17]. The mechanical effect may not be significant even when changing the reservoir temperature in this location. Unlike the geothermal reservoirs in a non-  volcanic area such as an enhanced geothermal system (EGS) or hot dry rock (HDR), the geothermal reservoirs in the volcanic area has a high thermal gradient and high permeability [47]. The naturally fractured rock formations have a sufficient supply of steam and hot water to generate electricity. Therefore, the hydrofracturing for enhancing the permeability of the rock matrix and fractures may not be necessary in this present study. However, the mechanical model is still required to make the realistic models of the reservoir, especially in real conditions. The rock deformation and fracturing in the volcanic geothermal reservoirs may occur in the active volcanic area caused by volcanic activity and collapse-resurgence processes [48]. Moreover, the reinjection of fluids into the reservoir may bring about induced seismicity [49]. The specific mechanical effects needs to be studied in the future.

Specific Gross Electrical Power Prediction
To identify the influence of the reservoir temperature evolution on power plant capacity, a thermodynamic analysis was conducted to estimate the gross electrical power of a geothermal resource. This analysis employs the properties of the fluid such as the specific enthalpy of saturated liquid (hf, kJ/kg), specific enthalpy of saturated vapor (hg, kJ/kg), specific entropy of saturated liquid (sf, kJ/kg/K), and specific entropy of saturated vapor (sg, kJ/kg/K) before and after the expansion process. After the temperature or pressure of the reservoir was predicted, the thermodynamic properties of fluid can be defined by using the table of properties for saturated water [45].
The Lahendong geothermal power plant has a single-flash geothermal system [41]. Therefore, the single-flash scheme and T-s cycle diagram ( Figure 21) were considered in determining the properties of the fluid [50]. Because there is no data of the fluid mass flow rate of each production well, the electrical power predictions can be calculated only by the specific gross electrical power (kW/kg/s) for each production well. More details of the gross electrical power calculation procedure at the Lahendong geothermal power plant are given in the literature [44]. Figure 22 shows the specific gross electrical power prediction at the production wells over 36 years. It can be seen that the effects of the reservoir temperature evolution on the specific gross electrical power are significant. The specific gross electrical power decreases with decrease of the temperature. The production well LHD-12 has specific gross electrical power higher than the others. However, after this production well starts to produce the geothermal fluid, the specific gross electrical power starts to decrease and has relative similar productivity with the production well LHD-10 which has lower specific gross electrical power than LHD-12 at the initial condition.

Conclusions
The evolution of the reservoir temperature at the Lahendong geothermal field has been evaluated through the use of a thermo-hydro-mechanical (THM)-coupled model. The computational model was defined based on the previous conceptual model, numerical study, and published data. Based on the results of the natural state condition, a good agreement was obtained between the numerical simulation results and the actual measured data. The reservoir temperature in the southern zone was seen to be higher than that in the northern zone. This means that the southern zone was suitable for the production well locations, while the northern zone was good for the reinjection well locations. The evolution of the reservoir temperature for long-term production (i.e., 36 years) was simulated, and the significant changes were observed. Moreover, the estimated temperatures of the reservoir were consistent with the thermal interpretations of H2S/H2 gas geothermometer studies. This suggests that the cold-water flow from the injection well to the production wells was not the main problem.
The main problem at this location may have been the decline in enthalpy, which may have caused the cease of fluid production at one production well (LHD-10) [46]. Therefore, as a vapordominated reservoir, the two-phase flow in the simulation should be considered in order to solve that problem. However, the reservoir at this geothermal field is still capable of supplying geothermal fluid during its lifetime. Although the mechanical effect does not significantly influence the natural state simulation, fluid migration and heat transfer became essential matters for this simulation. Furthermore, the effects of the chemical processes were not considered in this numerical simulation, but should be considered in future work in order to obtain more reasonable estimations of the energy production.