Well-Placement Optimization in an Enhanced Geothermal System Based on the Fracture Continuum Method and 0-1 Programming

Abstract

The well-placement of an enhanced geothermal system (EGS) is significant to its performance and economic viability because of the fractures in the thermal reservoir and the expensive cost of well-drilling. In this work, a numerical simulation and genetic algorithm are combined to search for the optimization of the well-placement for an EGS, considering the uneven distribution of fractures. The fracture continuum method is used to simplify the seepage in the fractured reservoir to reduce the computational expense of a numerical simulation. In order to reduce the potential well-placements, the well-placement optimization problem is regarded as a 0-1 programming problem. A 2-D assumptive thermal reservoir model is used to verify the validity of the optimization method. The results indicate that the well-placement optimization proposed in this paper can improve the performance of an EGS.

1. Introduction

Development and utilization of renewable energy have been a hot topic in society in recent years because of increased energy consumption and pollution [1]. Due to its reproducibility and cleanness, geothermal energy has received extensive attention. Most geothermal energy is preserved in hot dry rock (HDR) with a temperature between 150 °C to 650 °C in a depth range of about 3–10 km [2].

The enhanced geothermal system (EGS) proposed in the 1970s is the representative technology for HDR development [3]. The connected fracture network is formed in HDR through hydraulic fracturing, and the fractured thermal reservoir, called an artificial thermal reservoir, can be injected with cold water to extract the thermal energy [4].

The cost of an EGS in reservoir development and management is expensive, especially in well-drilling. In the process of constructing an EGS, the cost of well-drilling would account for more than 50% of the total cost because the hard reservoir and high temperature in hot dry rock could damage the drill bit quickly [5]. The optimal well-placement and operation are important to the performance of the EGS [6]. Combining numerical simulation with optimization algorithms is an effective method to search for the optimal well-placement.

The seepage in EGS during heat extraction is affected by multiple factors such as multi-field coupling [7], geometrical parameters of porous media [8] and fractures. The effects of fractures on seepage and heat extraction cannot be ignored because the EGS is often developed by hydraulic fracturing [9]. Many methods have been applied to research the fluid flow in the fracture, including the equivalent porosity model (EPM) [10], dual-porosity model (DPM) [11], digital core [12], discrete fracture network (DFN) [13], stochastic continuum (SC) [14], fracture continuum method (FCM) [15] and lattice Boltzmann methods (LBM) [16].

The optimization algorithm has been used in groundwater resource management [17], porous media model building [18] and oil reservoir management [19] for many years. Optimal well-placement [20,21], or well pattern [22,23], has largely improved the performance of the reservoir. In addition, there are several optimization algorithms applied to the well location or well pattern of the reservoir, such as the adjoint gradient algorithm [24,25], genetic algorithm (GA) [26], particle swarm optimization (PSO) [27], and new unconstrained optimization algorithm (NEWUOA) [28].

Some research about the well-placement in an EGS has been proposed. Akin et al. [29] optimized a new injection well-placement using simulated annealing based on the Kizildere geothermal field. A trained artificial neural network replaced the commercial simulators to reduce the processing demand. Chen et al. [30] proposed that a suitable well layout can improve the heat extraction effect. In addition, they designed a five-spot well layout and confirmed its performance using a numerical simulation. Chen et al. [31] used a multivariate adaptive regression spline (MARS) to set a surrogate model to replace the numerical model and used the bound optimization by quadratic approximation (BOBYQA) to optimize the well-placement. Wu et al. [32] studied the relationship between well-placement and heat extraction based on the semi-analytical solution model with a single horizontal fracture. Guo et al. [33] proposed that more production wells are more effective in delaying the breakthrough of the cold front, and the well should be placed at a position with higher rock stiffness.

There is less research on well-placement and all studies used the traditional method to encode the well-placement. On the other hand, they did not fully consider the uneven distribution of lots of fractures. In order to improve the performance of heat extraction, an optimization framework based on 0-1 programming and genetic algorithms is used in EGS well-placement. The purpose of this work is to provide a valid method to determine where the best locations of wells are in an EGS with a complex fracture network. The framework for EGS well-placement optimization consists of two parts: coding the well-placement variable with a 0-1 variable instead of the traditional coordinates of well-placement and reducing the computational cost by the FCM model. The first part is used to decrease the possible well-placements and it also has the potential to do joint optimization for well-placement and the number of wells. The FCM model is used to simplify the fractured reservoir model to reduce the computational costs while preserving the effect of fractures on heat extraction. An assumptive model is used to verify the validity of the method, and GA is used to search the best well-placement in EGS.

2. Method

2.1. Fracture Continuum Method

For a geothermal reservoir with a large number of fractures, the discrete fracture network model needs to discretize each fracture, which is too computationally expensive, making it unsuitable for optimization problems that require multiple iterations. In this work, the FCM is used to describe the flow of fluids in the thermal reservoir, and to minimize the computational cost of numerical simulation while preserving the effect of the fracture network on the fluid flow.

The FCM model can be considered as a stochastic continuum model that preserves the characteristics of the fracture distribution. In the FCM model, the reservoir is divided into several sub-grids, and the permeability of each sub-grid is determined by the fractures passing through that sub-grid.

2.1.1. Backbone Network’s Extraction

Matthäi [34] proposed that the disconnected fractures contribute little to the fluid flowing in the fractured reservoir when the permeability of the matrix is more than six orders of magnitude smaller than the fracture, and the effect of heat convection is much greater than heat conduction during heat extraction. Therefore, the disconnected fractures and the dead-end of fractures are eliminated and the backbone of the fracture network is extracted. Figure 1 shows an original fracture network and the backbone network extracted from the original network.

The permeability mapping in the next section is based on the backbone network, because the fractures that are not connected to each other are able to be connected in permeability mapping due to the dead end or disconnected fracture, which cause a higher permeability. Using the backbone network can largely avoid the higher permeability and retain the effect of the fracture on seepage.

2.1.2. Permeability Mapping Approach

The permeability mapping method is proposed in Ref [35]. An analytical method is used to calculate the permeability of each sub-grid. The permeability of the model in this work is expressed in tensor to preserve the effect of fracture direction on fluid flow. For a fracture with an angle $\theta$ to the x-axis, the permeability tensor of the fracture $k_F$ in the two-dimensional coordinate system can be expressed as:

$$k_F=k_f\left[\begin{array}{cc}{\cos }^2\theta & \sin \theta \cos \theta \\ \sin \theta \cos \theta & {\sin }^2\theta \end{array}\right]$$

(1)

where $\theta$ is the angle between the fracture and the x-axis; $k_f$ is the permeability of the fracture. The contribution of the fracture to the hydraulic conductivity of the sub-grid, which is crossed, can be estimated as $T_f/\Delta$ [36], where $T_f$ is the hydraulic conductivity of the fracture and $\Delta$ is the size of sub-grid. The relationship of the fracture to the permeability of the sub-grid, which the fracture passes through, can be expressed as:

$$k_c=\frac{d_f}{\Delta }{k}_f$$

(2)

where $k_c$ is the permeability contribution of fracture to the sub-grid; and $d_f$ is the width of the fracture. Therefore, the permeability of each sub-grid can be calculated as:

$$k_{i,j}=k_m+\frac{d_f}{\Delta }{k}_f\left[\begin{array}{cc}{\cos }^2\theta & \sin \theta \cos \theta \\ \sin \theta \cos \theta & {\sin }^2\theta \end{array}\right]$$

(3)

where $k_{i,j}$ is the permeability of sub-grid(i, j); and $k_m$ is the permeability of matrix. For cases where multiple fractures pass through the same sub-grid, the permeability sub-grid can be calculated as follows:

$$k_{i,j}=k_m+{\displaystyle \sum _{n=1}^N{k}_c}\left[\begin{array}{cc}{\cos }^2\theta & \sin \theta \cos \theta \\ \sin \theta \cos \theta & {\sin }^2\theta \end{array}\right]$$

(4)

where $N$ is the number of fractures passing through the sub-grid(i, j); and $k_{i,j}$ is the permeability of the sub-grid(i, j). Figure 2 shows the process of fracture permeability mapping.

Considering the error between the FCM model obtained after mapping and the DFN model, the permeability of FCM needs to be corrected as follows:

$$k_{i,j}^{\prime }=C{k}_{i,j}$$

(5)

where $k_{i,j}^{\prime }$ is the corrected permeability of the sub-grid(i, j); C is the permeability correction factor used to correct for the error in flow rate that occurs from mapping. In Ref [36] research C is calculated as $\left|\sin \theta \right|+\left|\cos \theta \right|$. In this work, the correction factor was calculated from the flow ratio between the DFN model and FCM model with uncorrected permeability.

2.2. Governing Equation

The fluid flowing in the thermal reservoir is described by Darcy’s law. The mass balance equation in the porous media is as follows:

$$\frac{\partial \left(\epsilon {\rho }_f\right)}{\partial t}+\nabla \cdot ({\rho }_{f}u)=Q_m$$

(6)

$$\frac{\partial \left(\epsilon {\rho }_f\right)}{\partial t}={\rho }_{f}S\frac{\partial P}{\partial t}$$

(7)

$$u=-\frac{k}{\mu }\nabla P$$

(8)

where $\epsilon$ is the porosity of the rock matrix; ${\rho }_f$ is the density of the fluid; $t$ is the time; $u$ is the Darcy velocity; $S$ is the storage coefficient; $P$ is the pressure; $k$ is the permeability of media; $\mu$ is the fluid dynamic viscosity; $Q_m$ is the source-sink term.

In this work, the local thermal non-equilibrium theory is used to describe the heat exchange between the rock and the fluid flowing in the geothermal reservoir. The energy balance equations are as follows [37,38]:

$$\left(1-\epsilon \right){\rho }_s{C}_{p,s}\frac{\partial T_s}{\partial t}=\nabla \cdot \left[\left(1-\epsilon \right){\lambda }_s\nabla T\right]+q_{sf}\left(T_f-T_s\right)$$

(9)

$$\epsilon {\rho }_f{C}_{p,f}\frac{\partial T_f}{\partial t}+\epsilon {\rho }_f{C}_{p,f}u\nabla T_f=\nabla \cdot \left(\epsilon {\lambda }_f\nabla T\right)+q_{sf}\left(T_f-T_s\right)$$

(10)

where ${\rho }_s$ is the density of the matrix; $T_s$ and $T_f$ are the temperatures of the matrix and fluid respectively; $C_{p,s}$ and $C_{p,f}$ are the specific heat capacities of the matrix and fluid respectively; ${\lambda }_s$ and ${\lambda }_f$ are the matrix and fluid thermal conductivities respectively; $q_{sf}$ is the interstitial convective heat transfer coefficient.

In Section 3 of this paper, the permeability of the FCM model needs to be corrected by using the results of the DFN model. The governing equations in the matrix are the same as that used in the FCM model.

The mass conservation equation in discrete fractures is written as:

$$d_f\frac{\partial {\rho }_f}{\partial t}+{\nabla }_T\cdot d_f{\rho }_f{u}_f=d_f{Q}_m$$

(11)

$$u_f=\frac{k_f}{\mu }{\nabla }_{T}P$$

(12)

where $u_f$ is the Darcy velocity in the fracture. The porosity of fractures is assumed to be 100%, so the temperature of the rock is not considered in the energy balance equation of fractures. The energy balance equation for the fluid in the discrete fractures is written as [39]:

$$Q_{fe}=d_f{\rho }_f{C}_p{}_f\frac{\partial T_{fr}}{\partial t}-d_f{\rho }_f{C}_p{}_f{u}_f{\nabla }_t{T}_{fr}-d_f{\lambda }_f{\nabla }_t{T}_{fr}$$

(13)

where $T_{fr}$ is the temperature of the fluid in fractures; $Q_{fe}$ is a source term to describe the heat transfer between the matrix and fractures, which mainly results from the heat convection.

2.3. Well-Placement Optimization of EGS FCM Model

Generally speaking, there are two principles in EGS well-placement design [30]: longer major flow path and less preferential flow. However, it is difficult to find a long major flow path directly without preferential or short-circuit flow which is a notorious issue annoying EGS researchers and engineers [40]. Combinations of optimization algorithms and numerical simulations provide an idea for solving this problem.

2.3.1. Well-Placement Optimization Problem with 0-1 Programming

When designing an EGS well-placement, all wells, including injection wells and production wells, should pass through fractures because of the low permeability of the matrix. In the FCM model used in this paper, the thermal reservoir is divided into several equal-sized sub-grids, the parameters of each sub-grid represent the fractures’ effect to this sub-grid. Therefore, all wells are located in the fractured sub-grids that have high permeability, which can ensure adequate connectivity between wells and reduce the number of potential well-placements. As shown in Figure 3, there are just 36 potential well-placements in a FCM model with 100 sub-grids.

However, this method of well-placement designing would bring some difficulties to the optimization. Complex fracture distribution in the reservoir results in uneven distribution of high permeability sub-grids, which makes it hard to deal with the constraints of well-placement. Transforming the well-placement optimization problem into a 0-1 programming problem can solve this difficulty.

In this work, the well-placement optimization problem of EGS is considered as a 0-1 programming problem. The one-well injection and multi-well production pattern are applied in this work. All wells are located in the high-permeability sub-grid and only one well at most on each sub-grid. For each grid, the grid without the well is recorded as 0, and the grid where the well is located in is recorded as 1. The coding form is shown in Figure 4, where the yellow grids represent the high-permeability sub-grids and the blue grids represent the low-permeability sub-grids.

This gives variables consisting of 0-1 to indicate the number and location of wells in the thermal reservoir. The high-permeability sub-grids have been numbered and the well-placement would be transferred to a one-dimensional vector consisting of binaries. The injection wells and production wells are not distinguished in coding, but the well closest to the center of the model is used as the injection well, and the other wells are used as production wells. The EGS well-placement problem based on 0-1 programming can be described as follows:

$$\begin{array}{l}\max \begin{array}{cc}& \end{array}f(x)\\ s.t.\begin{array}{ccc}& & {\displaystyle \sum _{i=1}^n{x}_i}\end{array}\le M\\ \begin{array}{cc}\begin{array}{cc}& \begin{array}{cc}& \end{array}\end{array}& {\displaystyle \sum _{i=1}^n{x}_i}\end{array}\ge N\\ \begin{array}{cc}\begin{array}{cc}& \end{array}& \begin{array}{cc}& x_i=\end{array}\end{array}\begin{array}{ccc}0& or& 1\end{array}\end{array}$$

(14)

where $x$ is a vector of 0s and 1s transformed from well-placement, $f(x)$ is the objective function, $M$ is the upper limit of the number of wells, and $N$ is the lower limit of the number of wells.

2.3.2. Genetic Algorithm

In this work, a GA is used to solve the EGS well-placement optimization problem. A GA is an optimization algorithm that searches for the best solution by simulating natural evolution. The iteration of a GA begins with a population of individuals. One individual represents a potential solution and the population represents a potential set of solutions to a problem. The individual is encoded by genes that represent the variables of the problem. The individual is evaluated by the fitness value determined by the user. The fitness value is the result of the objective function in most cases. The process of population regeneration consists of selection, crossover, and mutation. The role of selection is to eliminate individuals with low fitness, and crossover and mutation are used to generate new solutions to keep the diversity of the population. The best individual in the last generation is seen as the approximate optimal solution to the problem.

The strategies in GA adopted in this work are given below:

• Initialization: N individuals are randomly generated before iterations, which is used as the first generation in GA.
• Fitness calculation: the fitness (objective function) of each individual is calculated by a numerical simulation.
• Selection: roulette is used to select parent individuals from the current population, which means that individuals with greater fitness are more likely to be selected, and the selected individuals enter the parents pool.
• Crossover: do the single-point crossover of individuals in the parent pool based on crossover probability.
• Mutation: single-point mutation is employed to make small random changes in the individuals in the parent pool
• Elitist strategy: an elitist strategy is applied in the process of evolution. The individual with the best fitness in the current generation is retained to the next generation without crossover and mutation.
• Stopping criteria: when the number of generations achieves the pre-set value, GA will stop.
• Constraint: the constraint in this work is the number of wells. The first generation is initialized in the feasible region, and the infeasible solution generated in the iteration will be repaired.
• Repair method: the production well closest to the center would be removed if the number of wells is above the upper bound of the number of wells, and the well would be added at random locations if the number of wells is below the lower bound.

The GA is written in MATLAB R2018b which is easy to combine with numerical simulation modules.

3. A Well-Placement Optimization Case

3.1. Computational Model

Natural fractures [41] or hydraulic fracturing fractures [42] in the reservoir can often be obtained from history matching or seismic inversion [43]. In this work, an assumptive model sized 200 m × 200 m is used as the original fracture network and the parameters of fractures are referenced from Ref [35]. Two hundred fractures, which are divided into two sets with different dip angles (i.e., 0° and 90°), are generated, and the fractures’ lengths follow exponent distribution with the maximum length of 200 m, the minimum length of 20 m and exponent of 1. The backbone of the generated network is extracted using the method described in Section 2.1.1. The thickness of the permeable layer is 10 m. Figure 5 shows the original fracture network and the backbone network.

The model is divided into 20 × 20 square sub-grids with a side length of 10 m. The permeability mapping is based on the backbone network shown in Figure 5b. There is no permeability heterogeneity on the x-y and y-x components due to the direction of the fractures. The x-x component and y-y component of the permeability tensor are shown in Figure 6.

The mathematical model of the above governing equations including the matrix and fracture are discretized using the finite element method (the initial and boundary conditions are given in Section 3.2), which is solved using a commercial finite element software COMSOL Multiphysics 5.3 (COMSOL Co. Ltd., Stockholm, Sweden) and the GA written by MATLAB R2018b is linked to the software by LiveLink for MATLAB, which is developed by COMSOL. The parameters of GA are shown in

Table 1. The parameters of genetic algorithm (GA).

Parameters	Value
Population size	400
Max generation	40
Crossover rate	0.6
Mutation rate	0.02
Number of wells	5

.

3.2. Model Parameters

The reservoir is initially saturated with water. All wells, including 4 production wells and 1 injection well work under constant pressure conditions. The working fluid for heat extracting is also water. The parameters are written in

Table 2. Model Parameters.

Parameters	Value
Matrix density (kg/m3)	2700
Matrix porosity	0.01
Matrix permeability (m2)	1 × 10−17
Matrix heat capacity (J/(kg·K))	1000
Matrix heat conductivity (W/m·K)	3
Fracture permeability (m2)	1 × 10−10
Fracture width (m)	0.001
Water density (kg/m3)	1000
Water viscosity (Pa·s)	0.001
Water heat capacity (J/(kg·K))	4200
Water heat conductivity (W/m·K)	0.6
Storage coefficient (1/Pa)	1 × 10−10
Thickness of permeable stratum(m)	10
Correction factor	0.79

, which are referenced from some previous numerical studies [44,45]. The permeability of FCM has been corrected by the flow ratio of the DFN model and uncorrected FCM model.

No-flow and adiabatic boundaries are around the reservoir. The adiabatic boundary is set to better observe the effect of well-placement on heat extraction in temperature distribution. The initial and boundary conditions can be found in

Table 3. Initial and Boundary Conditions.

Conditions	Value
Initial pressure (MPa)	20
Initial temperature (°C)	200
Injection pressure (MPa)	30
Injection temperature (°C)	65
Production pressure (MPa)	20

.

In this work, we just consider a five-spot well-placement pattern with one injection well and four production wells. The injection well and production wells are not distinguished in the 0-1 code. The well closest to the center of the model would be used as an injection well while the other wells would be used as production wells. To facilitate analysis of single well performance, the production well number would be numbered clockwise from the well that is closest to that production well (0,0).

3.3. Objective Function

In this work, accumulative extracted thermal energy is used as the objective function of the well-placement optimization. Considering the adiabatic boundaries, it is equal to the decline in the thermal energy of the reservoir. It is expressed as follow [45]:

$$E={\displaystyle {\iiint }_{V_s}{\rho }_s{c}_{p,s}\left(T_i-T\left(t\right)\right)}dv$$

(15)

where $E$ is the decline in thermal energy in the reservoir; $T_i$ is the initial temperature; and $T(t)$ is the reservoir temperature in time $t$.

Besides $E$, the flow rate ($Q$), the accumulative extracted thermal energy ($\gamma$), the average production temperature ($\overline{T_{out}}$) and output thermal power ($p$) are also used to evaluate the performance of a geothermal reservoir with different well-placement. $Q$, $\gamma$, $\overline{T_{out}}$ and $p$ are defined as:

$$Q={\displaystyle {\int }_{L}u\left(t\right)}dl$$

(16)

$$\gamma ={\displaystyle {\int }_0^{t_s}Q{\rho }_f{C}_{p,f}\left(T_{out}-T_{in}\right)dt}$$

(17)

$$\overline{T_{out}}=\frac{{\displaystyle \int {}_L}T\left(t\right)dl}{L}$$

(18)

$$p=Q{\rho }_f{C}_{p,f}\left(T_{out}-T_{in}\right)$$

(19)

where $L$ is the length of the boundary of the well; $\gamma$ is the accumulative extracted thermal energy; $t_s$ is the simulation time; $T_{out}$ is the average temperature of production water; $T_{in}$ is the temperature of injection water; and $L$ is the length of the outlet boundary.

3.4. Results and Discussion

Genetic algorithms do not require a given initial solution, which is different from traditional optimization algorithms. Therefore, two different five-spot injection/production patterns (named Case 1 and Case 2) are used as two comparisons of the optimization result. As in Section 2.3, the yellow grids represent the high-permeability sub-grids and the blue grids represent the low-permeability sub-grids. In Case 1 the production well 1 and 3 is located on a sub-grid without fractures, while in Case 2 all wells are set in the sub-grids passed through fractures just like the optimization result. Figure 7 shows the comparisons and the optimization result (named Case 3).

Figure 8 and Figure 9 show the best individual of four different generations and the convergence process of the objective function. In Figure 9, the solid blue line indicates the best fitness in each generation, and the orange dotted lines indicate the average fitness of each generation.

The convergence process shows that the best fitness achieved a high value in previous generations, which may be due to the fact that the wells are always placed in a high-permeability sub-grid during the optimization process. The lower average fitness in previous generations illustrates that many low fitness individuals are generated during the population initialization and genetic manipulation of previous generations, and the rapid increase in average fitness indicates that the entire population is evolving, which can prove the validity of 0-1 programming and GA in geothermal well-placement optimization.

The temperature variations of three cases are shown in Figure 10. The distribution of temperature is similar to each other at first because of the near location of the injection well. Gradually the cold front of each case migrates with the high-permeability sub-grid and the placement of production wells. It is obvious that the migrations of the cold front in Case 2 and Case 3 are faster than Case 1. The cold front in Case 1 only migrated to the Pro2 and Pro4 that are located in high-permeability sub-grids caused by the connected fracture. From the temperature variation, it can be observed that the effect of fractures to seepage and heat extraction is preserved in the FCM model.

The difference in heat extraction between Case 2 and Case 3 is not as obvious as between Case 1 and other cases, but it can be found that the temperature of the northeast in Case 3 is lower than Case 2 and the low-temperature region in Case 3 is larger than Case 2 overall. Considering that there is no supply source, the heat extraction in Case 3 is more adequate. The accumulative extracted thermal energy of the three cases are plotted in Figure 11.

As shown in Figure 11, the final $\gamma$ of Case 3 from GA is higher than the two five-spot patterns that set to compare, which can also prove the validity of the well-placement method applied in this work. It also can be found that the heat recovery rate of Case 2 is higher than that of Case 3 in the first 1500 days. Figure 12 shows the change in output thermal power in the three cases, which is consistent with Figure 11. As shown in Figure 12, the power of Case 2 is highest in the first 700 days, but it also has the fastest decline. After 2500 days, the power of Case 2 is the least in all three cases.

At the initial running stage, the higher heat recovery rate indicates a better flow connection between the production well and the injection well, which means there are more fractures connected with the production well and injection well, but it does not mean the final performance will be better. Preferential or short-circuit flow in the thermal reservoir is always a headache for geothermal development and management. High flow velocity may cause a rapid decrease in matrix temperatures beside the connected fractures, which decrease the efficiency of heat convection. The average temperature shown in Figure 13 and the flow rate shown in Figure 14 can prove it.

The flow rate is rapidly stable because of the little storage coefficient. As shown in Figure 13, the average production temperature of Case 2 has a fast drop. It can be inferred from the temperature and flow rate that the preferential flow exists in Case 2.

Figure 15, Figure 16 and Figure 17 show the accumulative energy, average temperature and the flow rate of each production well in the three cases. Consistent with the temperature distribution, the production wells, Pro1 and Pro3 in Case 1 contribute little to heat extraction, and the $\overline{T_{out}}$ of Pro3 shows the preferential flow in Case 2 mainly exists between the injection well and Pro3, and the optimization result has been improved in it.

4. Conclusions and Future Work

A well-placement optimization framework is proposed in this paper. FCM is used to simplify the fractured thermal reservoir model and the GA is used to solve the well-placement optimization problem that was considered as a 0-1 programming problem.

• The developed framework is efficient in the EGS well-placement optimization problem. The extracted thermal energy, which was the objective function, has increased in the convergence process of GA. And the optimization result shows better performance than comparison.
• The FCM model can reflect the effect of fractures on seepage and heat transfer to a certain extent.
• Regarding the well-placement optimization problem as a 0-1 programming problem can reduce the potential well-placements and improve the optimization effect. It also has the potential in joint optimization for well-placement and the number of wells.
• In the well-placement design of EGS, the connectivity between the injection well and production well should be considered as the primary factor. The well in low-permeability contributes little to heat extraction.
• Strong connectivity between wells does not mean better performance. Strong connectivity may lead to preferential flow and early heat breakthrough.

In this study, the framework only includes the 2-D model and the vertical well. In the future, this work will be generalized to the 3-D multi-field coupling model, a horizontal well, and the joint optimization of well-placement and the number of wells. A more advanced algorithm will be applied in EGS well-placement optimization, such as multi-objective optimization [46] and machine learning [47].