Numerical Investigation of Techno-Economic Multiobjective Optimization of Geothermal Water Reservoir Development: A Case Study of China

: Renewable geothermal utilization is a signiﬁcant approach for residential heating with two principal modes, direct geothermal district heating systems (DGDHSs) and indirect geothermal district heating systems (IGDHSs). The key principle of geothermal development design is to prevent premature thermal breakthrough, which could result in low e ﬃ ciency of geothermal heating systems. In this paper, a new approach considering building heating demand, geothermal water resource protection, and optimal economic beneﬁts is presented systematically. The results simulated by OGS software show that well spacing, reinjection temperature, and production rate are the most signiﬁcant parameters a ﬀ ecting thermal breakthrough in geothermal reservoirs. In addition, production rate and reinjection temperature have a huge e ﬀ ect on the payback period of investment. Comparing IGDHS to DGDHS, the investment in construction of geothermal wells and the annual water consumption decrease by up to 10% and 50%, respectively. Additionally, electricity costs increase by 5% to 30%. The indirect geothermal district heating system with a well spacing of 300 m, a production rate of 100 m 3 / h, and a reinjection temperature of 301.15 K is much better for this case, both technically and economically. The systematic calculation approach can be reasonably applied to other regions with geothermal energy utilization.


Introduction
With the development of the economy and improvement of people's living standards, primary energy consumption has maintained a rapid growth. In northern China, energy consumption of building heating and space cooling account for 20% and 14%, respectively [1]. One measure for energy conservation and emission reduction is to improve building insulation levels and the air tightness of building envelopes. Another measure is to use geothermal resources. Geothermal energy is a stable resource from the earth's interior [2,3]. Theoretically, it can be utilized continuously [4]. As long as the upper limit of geothermal utilization is properly controlled, the service life of a geothermal reservoir can be extended [5,6]. Besides, geothermal energy is generally unaffected by geography, climate, and season, unlike solar, water, wind, biomass, and ocean energy [7,8].
Geothermal heating technology is an effective way to fully develop and utilize geothermal energy [9,10]. It has been applied in some countries because of its high operational efficiency and significant environmental benefits [11]. However, the high initial investment and long payback period have hindered the popularization of this technology [12,13]. Geothermal well infrastructure investment can exceed 50% of system construction costs for all types of geothermal utilization systems [14,15]. Statistically, the expense of drilling a geothermal well is about $300 per meter [16]. The cost increases with deepening of the wells [17]. In addition, all extracted groundwater needs to be reinjected to maintain the pressure of the geothermal reservoir. The reinjection rate is far lower than the production rate, resulting in a waste of ground water [18]. Cold reinjection water also leads premature thermal breakthrough, which affects the system's service behavior and operating efficiency.
Researchers have studied the energy production of geothermal temperature fields. Zhang et al. developed a numerical model to evaluate geothermal reservoirs and optimize the production performance of geothermal systems [19]. Held et al. derived a reservoir mathematical model that considered fluid-rock interaction in the utilization of geothermal resources [20]. Nam carried out 3D numerical simulations and verified the reliability by comparing results with experimental results [21]. Salimzadeh et al. set up a fully coupled thermal-hydraulic-mechanical (THM) finite element model [22]. Pandey et al. found that thermochemical and thermomechanical effects were the most important consequences of heterogeneity [23]. Nikhil et al. simplified the crack model by a statistical linear congruent generator method, and carried out sensitivity analysis of fracture aperture size and fluid flow [24].
The previous efforts offered many insights into modeling heat extraction from geothermal reservoirs. Rock deformation should be emphasized, which may lead to unwanted sand components in the water flow [25]. However, previous studies ignored the economic benefits of geothermal projects. The disadvantages of geothermal projects without economic benefit analysis were huge. It not only increased the cost, but also caused wasted geothermal resources. High investment and operating cost hindered the development of geothermal heating technology.
In this study, based on a coupled mathematical and multiobjective optimization model, numerical investigation is carried out to evaluate the heat extraction and economic efficiency of a geothermal project in Xinji County. First, coupling of the thermal-hydraulic-mechanical model and multiobjective function of the geothermal heating system is established. Second, the main factors causing thermal breakthrough are analyzed based on data of a geothermal reservoir in Xinji City. Third, production parameters of direct and indirect geothermal district heating systems are optimized based on the ground heating load and system operation strategies. Finally, the initial investment and operation cost of direct and indirect geothermal district heating systems are analyzed, and the optimal method and production parameters are thus obtained. This optimization and economic evaluation method can be applied to geothermal heating systems in other regions, which is beneficial to the further development of geothermal heating technology.

Background
The area studied in this paper is located in Xinji City, Hebei, China. The heat-providing area is 250,100 m 2 , which includes residential and commercial areas. The designed heating load is 9040 kW. After this project is put into operation, the annual heating consumption will be 65,600 GJ.
The sedimentary Miocene Guantao Formation is the most favorable for geothermal heating systems [26,27]. According to a geological investigation conducted by China Petroleum and Chemical Corporation, the geothermal resources of Xinji is regarded as consisting of homogeneous porous media with high permeability. The reservoir thickness is 512 m and the average porosity of this geothermal reservoir is 30%. The depth of the bottom boundary varies from 1600 to 2100 m. The single-well water rate is between 78 and 120 m 3 /h, and the production temperature is in the range of 333.15-3337.15 K. Laboratory tests show that the reservoir rock permeability is around 160 mD. The average reservoir temperature at target layers around 1300 m deep is about 333.15 K. The ground water of this geothermal reservoir can be directly used for building heating due to the low relative temperature [26].
To ensure the stable operation of a geothermal system, the geological conditions and well test data need to be considered in the design of the heating system [27]. Reinjection operation data of Xinji oil production well no. 5 are shown in Figure 1. Except during the filter element replacement period, the reinjection rate is mostly maintained at 120 m 3 /h, which ensures production water being totally reinjected underground. water of this geothermal reservoir can be directly used for building heating due to the low relative temperature [26].
To ensure the stable operation of a geothermal system, the geological conditions and well test data need to be considered in the design of the heating system [27]. Reinjection operation data of Xinji oil production well no. 5 are shown in Figure 1. Except during the filter element replacement period, the reinjection rate is mostly maintained at 120 m 3 /h, which ensures production water being totally reinjected underground.

Geothermal Heating Strategy
There are two kinds of geothermal heating modes designed for Xinji. One is a direct geothermal district heating system (DGDHS). As shown in Figure 2a, ground water extracted from the geothermal production well is transported to the heat exchanger. The heating circulation water absorbs heat from the geothermal water. Then, the heat is piped to residential areas. The other one is an indirect geothermal district heating system (IGDHS), shown in Figure 2b. The heat of geothermal water is partially transferred to the heating circulation water at the primary heat exchanger with higher returning water temperature, and the heat is piped to some residential areas. After the initial heat transfer, the geothermal water is transported to the sequential secondary heat exchanger, where the heat is transferred from the geothermal water to the intermediate circulation water. Then, the heat pump extracts the heat from the intermediate circulation water and discharges it into the heating circulation water. Finally, the heat is transported to the remaining residential areas.

Geothermal Heating Strategy
There are two kinds of geothermal heating modes designed for Xinji. One is a direct geothermal district heating system (DGDHS). As shown in Figure 2a, ground water extracted from the geothermal production well is transported to the heat exchanger. The heating circulation water absorbs heat from the geothermal water. Then, the heat is piped to residential areas. The other one is an indirect geothermal district heating system (IGDHS), shown in Figure 2b. The heat of geothermal water is partially transferred to the heating circulation water at the primary heat exchanger with higher returning water temperature, and the heat is piped to some residential areas. After the initial heat transfer, the geothermal water is transported to the sequential secondary heat exchanger, where the heat is transferred from the geothermal water to the intermediate circulation water. Then, the heat pump extracts the heat from the intermediate circulation water and discharges it into the heating circulation water. Finally, the heat is transported to the remaining residential areas. water of this geothermal reservoir can be directly used for building heating due to the low relative temperature [26].
To ensure the stable operation of a geothermal system, the geological conditions and well test data need to be considered in the design of the heating system [27]. Reinjection operation data of Xinji oil production well no. 5 are shown in Figure 1. Except during the filter element replacement period, the reinjection rate is mostly maintained at 120 m 3 /h, which ensures production water being totally reinjected underground.

Geothermal Heating Strategy
There are two kinds of geothermal heating modes designed for Xinji. One is a direct geothermal district heating system (DGDHS). As shown in Figure 2a, ground water extracted from the geothermal production well is transported to the heat exchanger. The heating circulation water absorbs heat from the geothermal water. Then, the heat is piped to residential areas. The other one is an indirect geothermal district heating system (IGDHS), shown in Figure 2b. The heat of geothermal water is partially transferred to the heating circulation water at the primary heat exchanger with higher returning water temperature, and the heat is piped to some residential areas. After the initial heat transfer, the geothermal water is transported to the sequential secondary heat exchanger, where the heat is transferred from the geothermal water to the intermediate circulation water. Then, the heat pump extracts the heat from the intermediate circulation water and discharges it into the heating circulation water. Finally, the heat is transported to the remaining residential areas.

Thermal-Hydraulic-Mechanical Model
In this paper, the problems of water flow, heat transport, and rock deformation in the Xinji geothermal reservoir are studied jointly. The process of groundwater extraction is considered to be a flowing process of liquid in porous media consisting of a liquid phase and a solid phase [28]. The process of modeling and theoretical analysis consists of the following assumptions [28][29][30][31][32][33]: (1) As the microstructures are well connected, thus, the hydraulic and transport characteristics of the rock matrix can be described by averaged quantities. (2) Local thermodynamic equilibrium is assumed between liquid phase and solid phase.
(3) There is diffusion, convective and conductive heat transfer in porous media. Radiation heat transfer, which has little effect on geothermal water, is ignored. (4) As total fluid balance (not transported species) and yields are considered in the model, there is no water loss in the geothermal reservoir.
The mass balance equation, energy balance equation, and rock deformation in porous media are shown in Equations (1)-(3) [31,32]: where φ is rock porosity; ρ f is water density (kg/m 3 ); t is time (s); Q f is the mass flow of sink/source item (kg/(m 3 ·s)); ρc p is the heat capacity of porous media (J/(m 3 ·K)); T is temperature (K); λ is the thermal conductivity of porous media (W/(m·K)); σ is the typical Cauchy stress tensor; and u f s is fluid flow velocity between rock and water flow (m/s). Both the heat capacity and thermal conductivity of porous media can be estimated as an arithmetic mean of each phase property weighted by its volume fraction. Then, the flow velocity u is divided into the rock velocity u s and the relative velocity of fluid u f s . Huo et al. and Salimzadeh et al. proved that u s is associated with the expansion coefficient of rock and fluid ∂ε ∂t , and ε is calculated from the displacement solved by Equation (3) [32,33]. The term u s is used to couple the water flow with rock deformation. Equation (1) is rewritten into Equation (4): where p is fluid pressure (Pa); c t is total compressibility (Pa −1 ), c t = c f + c s ; c f is water compressibility (Pa −1 ); and c s is rock compressibility (Pa −1 ). In Equation (4), the effect of rock deformation caused by geothermal water extraction is considered. Thus, the flow velocity u f s + u s depicts the movement of water in porous media. This is a very important step. Although u s is usually a very small value, it cannot be neglected, as it accounts for the rock deformation mechanism.
As shown in Equation (5), u fs further expressed by Darcy's law that shows the heat-fluid coupling process: where g is gravitational acceleration (m/s 2 ); z is the depth of the geothermal reservoir (m); and k is the intrinsic permeability tensor (m 2 ). In order to specify the solution for the above equations, one needs to prescribe appropriate boundary conditions. The details are shown in [29,30].
(1) Prescribed pressure (Dirichlet condition) for imposing hydrostatic conditions at domain boundaries or pressure at wellbores.
(2) Prescribed fluid flux (Neumann condition) with the normal vector n for imposing inflow or outflow rate at wellbores.
(3) Prescribed temperature (Dirichlet condition) for imposing natural thermal gradients at model boundaries and the temperature of the injected fluid. The finite element method is used in Open-GeoSys (OGS) to provide numerical solutions to the aforementioned coupled formulation. More information can obtained from [29][30][31][32][33].

Calculation of Well Number based on Heating Load
Based on the heating load of residential buildings, the number of geothermal production wells is calculated, which is used to calculate construction costs. The building heat load is estimated based on local weather conditions, the function of the building, and characteristics of the protected exterior construction [34]. The calculation formula is shown in Equation (9): where W c is heat load per unit area of residential building (W/m 2 ); W b is heat load per unit area of commercial building (W/m 2 ); S c is residential heat-providing area (m 2 ); S b is commercial heat-providing area (m 2 ); and Q l is total heat load (W). The number of geothermal production wells is determined by the production rate and heating mode, shown in Equations (10) and (11): where n 1 is the number of production wells for DGDHS (rounded up); α is the heat loss coefficient of the pipe network (5%-10%); c w is the specific heat of water (J/(kg·K)); t g is production temperature (K); t h is reinjection temperature (K); and q p is production rate (m 3 /h); and where n 2 is the number of production wells for IGDHS (rounded up); t h1 is water temperature for the first return (K); t h2 is water temperature for the second return (K); and COP is the coefficient of performance.

Multiobjective Optimization Model for Heating Systems
In this part, a multiobjective optimization model for geothermal heating systems is established. The project's profitability is analyzed by coupling the initial investment and operating cost of the system. The initial investment, annual operating cost, and revenue are shown in Figure 3

Multiobjective Optimization Model for Heating Systems
In this part, a multiobjective optimization model for geothermal heating systems is established. The project's profitability is analyzed by coupling the initial investment and operating cost of the system. The initial investment, annual operating cost, and revenue are shown in Figure 3 [35][36][37]. Geothermal well construction costs account for the main part of all construction costs which are based on well number, type, and construction depth. Due to limited drilling area, a vertical well is adopted for geothermal well construction, and the rest are directional wells, shown in Figure 4. The construction depth of a directional well is related to well spacing and directional angle, as in Equation (12), and the construction cost of the project can be estimated according to Equation (13): where 1 h is the depth of the vertical well (m); 2 h is the depth of the directional well (m); L is the depth of the first section (m); and  is directional angle (°); and Geothermal well construction costs account for the main part of all construction costs which are based on well number, type, and construction depth. Due to limited drilling area, a vertical well is adopted for geothermal well construction, and the rest are directional wells, shown in Figure 4. The construction depth of a directional well is related to well spacing and directional angle, as in Equation (12), and the construction cost of the project can be estimated according to Equation (13): where h 1 is the depth of the vertical well (m); h 2 is the depth of the directional well (m); L is the depth of the first section (m); and θ is directional angle ( • ); and where C z is construction cost of a geothermal well (k USD); P 1 is construction cost per meter of vertical well (k USD); P 2 is construction cost per meter of directional well (k USD); C f is cost of attached ground facility (k USD); and n is total number of geothermal heating wells. The payback period is the time required for the return on investment to "repay" the total investment of the geothermal heating system [38], as expressed in Equation (14), and is calculated from the beginning year of construction [39,40]: where Pt is static payback time (year); CI is cash inflow; CO is cash outflow; and (CI − CO) t is net cash flow for the year, t.
where Pt is static payback time (year); CI is cash inflow; CO is cash outflow; and ( ) t CI CO  is net cash flow for the year, t. There are two main objectives of geothermal heating system optimization: one is to have the shortest payback period, and the other is to have minimum well spacing. In order to improve the profitability of the system, its production parameters must be optimized. However, three factors restrict the optimization of a geothermal heating system. The first is the area of equipment, which must be smaller than the project area. This factor is becoming more significant due to the limited space of urban areas. Second, a geothermal system is affected by groundwater resources, which is related to the hydrogeological conditions in the specific region. Finally, the system must meet the people's requirements for thermal comfort in the built environment. The multiobjective optimization functions are shown in Equation (15).
Optimization function: There are two main objectives of geothermal heating system optimization: one is to have the shortest payback period, and the other is to have minimum well spacing. In order to improve the profitability of the system, its production parameters must be optimized. However, three factors restrict the optimization of a geothermal heating system. The first is the area of equipment, which must be smaller than the project area. This factor is becoming more significant due to the limited space of urban areas. Second, a geothermal system is affected by groundwater resources, which is related to the hydrogeological conditions in the specific region. Finally, the system must meet the people's requirements for thermal comfort in the built environment. The multiobjective optimization functions are shown in Equation (15).
Optimization function: minl(x 1 , x 2 , x 3 . . . x n ) Constraint conditions: Scale constraint: where l is well spacing (m); S g is the area of equipment (m 2 ); S is project area (m 2 ); q x is the design value of production rate (m 3 /h); q pmax is the maximum of production rate (m 3 /h); Q x is the design value of system load (W); and x 1 , x 2 , x 3 . . . x n represent different working conditions. If hydrogeological conditions in a specific region allow the use of geothermal heating systems, these multiobjective optimization functions can be used for assessment.

Simulation Procedure
According to the established thermal-hydraulic-mechanical and multiobjective optimization models, using simulation software, the calculation process is shown in Figure 5. A geothermal heating system can be optimized through the entire calculation process. First, formation parameters are input into the model. The movement of the cold water thermal front is calculated. By comparing the temperature of the production well (T P0 is the initial temperature of the production well, and T P50 is the temperature at the 50th year), the optimal spacing of wells is selected. Finally, the initial investment and operating cost are calculated. Unprofitable geothermal heating systems are simply discarded. Though the multiobjective optimization model, the optimal production parameters can be obtained.

Sensitivity Analysis
This part studies the influence mechanism of various factors on the formation of thermal breakthrough. The basic conditions of the Xinji geothermal reservoir provided by China Petroleum and Chemical Corporation are listed in Table 1. These parameters, such as well spacing, production rate, and reinjection temperature, affect the formation of thermal breakthrough. They also determine the heat extraction efficiency of a geothermal heating system. In addition, periodic heating is considered in the modeling: a geothermal heating system works for 5 months each winter. Geothermal wells are closed during the remainder of the year. Results are obtained at the final time

Sensitivity Analysis
This part studies the influence mechanism of various factors on the formation of thermal breakthrough. The basic conditions of the Xinji geothermal reservoir provided by China Petroleum and Chemical Corporation are listed in Table 1. These parameters, such as well spacing, production rate, Water 2019, 11, 2323 9 of 17 and reinjection temperature, affect the formation of thermal breakthrough. They also determine the heat extraction efficiency of a geothermal heating system. In addition, periodic heating is considered in the modeling: a geothermal heating system works for 5 months each winter. Geothermal wells are closed during the remainder of the year. Results are obtained at the final time step at the end of the 50th year. Well spacing is the distance between a heat production well and a reinjection well, which determines the placement of geothermal wells. In this section, different well spacings (200, 400, and 600 m) are set into the mathematical model. The remaining parameters are as follows: production rate is 100 m 3 /h and reinjection temperature is 311.15 K. Temporal evolutions of the water temperature and pressure at the production well are plotted in Figure 6. Well spacing is the distance between a heat production well and a reinjection well, which determines the placement of geothermal wells. In this section, different well spacings (200, 400, and 600 m) are set into the mathematical model. The remaining parameters are as follows: production rate is 100 m 3 /h and reinjection temperature is 311.15 K. Temporal evolutions of the water temperature and pressure at the production well are plotted in Figure 6.
In general, large well spacing can slow down the formation of thermal breakthrough. It can be observed from Figure 6a that when the well spacing is 200 m, the temperature of the production well begins to decline at the 25th year. The reason is that the movement of cold water thermal front reaches the production wells. This temperature change would reduce the heat extraction efficiency of a geothermal heating system. Eventually, ground heating demands are not met. For the other well spacings, the temperature of production wells is constant at 333.15 K. Note that there is no thermal breakthrough. As shown in Figure 6b, the pressure of the production well decreases with increased well spacing. This is because the reinjected water is better able to repressurize the depleted area around the production well with smaller spacing. It means that large well spacing can effectively avoid premature thermal breakthrough. The movement of the cold water thermal front takes a long time to arrive at the production well. However, it loses the significance of reinjection to maintain water pressure. As commercial, land use, economic, and policy factors should be considered in the selection of geothermal well spacing, optimizing the well spacing is a challenging task. In general, large well spacing can slow down the formation of thermal breakthrough. It can be observed from Figure 6a that when the well spacing is 200 m, the temperature of the production well begins to decline at the 25th year. The reason is that the movement of cold water thermal front reaches the production wells. This temperature change would reduce the heat extraction efficiency of a geothermal heating system. Eventually, ground heating demands are not met. For the other well spacings, the temperature of production wells is constant at 333.15 K. Note that there is no thermal breakthrough.
As shown in Figure 6b, the pressure of the production well decreases with increased well spacing. This is because the reinjected water is better able to repressurize the depleted area around the production well with smaller spacing. It means that large well spacing can effectively avoid premature thermal breakthrough. The movement of the cold water thermal front takes a long time to arrive at the production well. However, it loses the significance of reinjection to maintain water pressure. As commercial, land use, economic, and policy factors should be considered in the selection of geothermal well spacing, optimizing the well spacing is a challenging task.

Effect of Reinjection Temperature
The results of first part show that there is an apparent decline in production temperature of a geothermal heating system with well spacing of 200 m and a production rate of 100 m 3 /h. These parameters are helpful to visualize the influence of reinjection temperature on thermal breakthrough. In this case, four values of the parameters (288.15, 298.15, 308.15, and 318.15 K) are considered. Temporal evolution of the water temperature and pressure at a geothermal production well is plotted in Figure 7. The results of first part show that there is an apparent decline in production temperature of a geothermal heating system with well spacing of 200 m and a production rate of 100 m 3 /h. These parameters are helpful to visualize the influence of reinjection temperature on thermal breakthrough. In this case, four values of the parameters (288.15, 298.15, 308.15, and 318.15 K) are considered. Temporal evolution of the water temperature and pressure at a geothermal production well is plotted in Figure 7. The temperature of reinjection water has a great influence on the formation of thermal breakthrough and the temperature of the production well. As is observed in Figure 7, thermal breakthrough is formed first when the reinjection temperature is 288.15 K. The lower the reinjection temperature is, the earlier the thermal breakthrough is formed, and the more the production water temperature drops. This is because the temperature decrease caused by the reinjected water results in the movement of the cold water thermal front originating from the reinjection well.
Although lower reinjection temperature leads to premature thermal breakthrough, it also points to a high utilization ratio of geothermal resources. Therefore, reinjection temperature must be controlled reasonably in geothermal heating projects.

Effect of Production Rate
The production rate is the water production rate at the production well. All geothermal water must be reinjected underground, so the reinjection rate is always equal to the production rate. According to the previous analysis of well spacing and reinjection temperature, when the well spacing is 200 m and the reinjection temperature is 308.15 K, the temperature of the production well changes markedly. Therefore, those parameters help to observe the sensitivity of temperature to production rates. Different production rates (80, 100, and 120 m 3 /h) are input into the model. The results are shown in Figure 8. The temperature of reinjection water has a great influence on the formation of thermal breakthrough and the temperature of the production well. As is observed in Figure 7, thermal breakthrough is formed first when the reinjection temperature is 288.15 K. The lower the reinjection temperature is, the earlier the thermal breakthrough is formed, and the more the production water temperature drops. This is because the temperature decrease caused by the reinjected water results in the movement of the cold water thermal front originating from the reinjection well.
Although lower reinjection temperature leads to premature thermal breakthrough, it also points to a high utilization ratio of geothermal resources. Therefore, reinjection temperature must be controlled reasonably in geothermal heating projects.

Effect of Production Rate
The production rate is the water production rate at the production well. All geothermal water must be reinjected underground, so the reinjection rate is always equal to the production rate. According to the previous analysis of well spacing and reinjection temperature, when the well spacing is 200 m and the reinjection temperature is 308.15 K, the temperature of the production well changes markedly. Therefore, those parameters help to observe the sensitivity of temperature to production rates. Different production rates (80, 100, and 120 m 3 /h) are input into the model. The results are shown in Figure 8. A large production rate promotes the movement of the cold water thermal front toward the production well. As observed Figure 8a, with decreased production rate, the temperature decrease at the production well is delayed, and the drop in production temperature is not obvious. This is because the decreased well rate slows down the spread of the temperature decrease initiated by the movement of the cold water thermal front.
According to Figure 8b, the greater the production rate, the lower the pressure. This is because a lower well rate results in a smaller pressure gradient between the reinjection well and the production well.
A small production rate delays the temperature drop of production water. However, it also results in poor heating extraction of the geothermal heating system. Therefore, a reasonable range of production rates should be controlled on the basis of local test data.

Parameter Optimization
Taking IGDHS with a production rate of 100 m 3 /h and reinjection temperature of 301.15 K as an example, results are shown in Figure 9. Figure 9a shows that production temperature does not change when well spacing is over 300 m. As can be seen from Figure 9b, with increased well spacing, the pressure of the production well keeps dropping. This is similar to the first part of the simulation results, proving that the results are correct. A large production rate promotes the movement of the cold water thermal front toward the production well. As observed Figure 8a, with decreased production rate, the temperature decrease at the production well is delayed, and the drop in production temperature is not obvious. This is because the decreased well rate slows down the spread of the temperature decrease initiated by the movement of the cold water thermal front.
According to Figure 8b, the greater the production rate, the lower the pressure. This is because a lower well rate results in a smaller pressure gradient between the reinjection well and the production well.
A small production rate delays the temperature drop of production water. However, it also results in poor heating extraction of the geothermal heating system. Therefore, a reasonable range of production rates should be controlled on the basis of local test data.

Parameter Optimization
Taking IGDHS with a production rate of 100 m 3 /h and reinjection temperature of 301.15 K as an example, results are shown in Figure 9. Figure 9a shows that production temperature does not change when well spacing is over 300 m. As can be seen from Figure 9b, with increased well spacing, the pressure of the production well keeps dropping. This is similar to the first part of the simulation results, proving that the results are correct. A large production rate promotes the movement of the cold water thermal front toward the production well. As observed Figure 8a, with decreased production rate, the temperature decrease at the production well is delayed, and the drop in production temperature is not obvious. This is because the decreased well rate slows down the spread of the temperature decrease initiated by the movement of the cold water thermal front.
According to Figure 8b, the greater the production rate, the lower the pressure. This is because a lower well rate results in a smaller pressure gradient between the reinjection well and the production well.
A small production rate delays the temperature drop of production water. However, it also results in poor heating extraction of the geothermal heating system. Therefore, a reasonable range of production rates should be controlled on the basis of local test data.

Parameter Optimization
Taking IGDHS with a production rate of 100 m 3 /h and reinjection temperature of 301.15 K as an example, results are shown in Figure 9. Figure 9a shows that production temperature does not change when well spacing is over 300 m. As can be seen from Figure 9b, with increased well spacing, the pressure of the production well keeps dropping. This is similar to the first part of the simulation results, proving that the results are correct. Temperature distribution nephograms of geothermal reservoirs under different situations are shown in Figures 10-12. The results show that the radius of the cold water thermal front is directly proportional to the production rate, and the radius of IGDHS is larger than that of DGDHS. In addition, the radius of the cold water thermal front increases with decreased reinjection temperature. As shown in Figure 11b, the radius of the cold water thermal front is about 260 m. Considering that the uneven distribution of formation and the instability of groundwater accelerate the flow of low-temperature groundwater, it is recommended that well spacing should be over 300 m in this case. The production parameters of each geothermal heating system are shown in Table 2.
Water 2019, 11, 2323 13 of 18 Figure 9. Time evolution of a production well for indirect geothermal district heating systems (IGDHSs) with production rate of 80 m 3 /h: (a) temperature at production well; (b) pressure at production well.
Temperature distribution nephograms of geothermal reservoirs under different situations are shown in Figures 10-12. The results show that the radius of the cold water thermal front is directly proportional to the production rate, and the radius of IGDHS is larger than that of DGDHS. In addition, the radius of the cold water thermal front increases with decreased reinjection temperature. As shown in Figure 11b, the radius of the cold water thermal front is about 260 m. Considering that the uneven distribution of formation and the instability of groundwater accelerate the flow of lowtemperature groundwater, it is recommended that well spacing should be over 300 m in this case. The production parameters of each geothermal heating system are shown in Table 2.      Figure 9. Time evolution of a production well for indirect geothermal district heating systems (IGDHSs) with production rate of 80 m 3 /h: (a) temperature at production well; (b) pressure at production well.
Temperature distribution nephograms of geothermal reservoirs under different situations are shown in Figures 10-12. The results show that the radius of the cold water thermal front is directly proportional to the production rate, and the radius of IGDHS is larger than that of DGDHS. In addition, the radius of the cold water thermal front increases with decreased reinjection temperature. As shown in Figure 11b, the radius of the cold water thermal front is about 260 m. Considering that the uneven distribution of formation and the instability of groundwater accelerate the flow of lowtemperature groundwater, it is recommended that well spacing should be over 300 m in this case. The production parameters of each geothermal heating system are shown in Table 2.      Figure 9. Time evolution of a production well for indirect geothermal district heating systems (IGDHSs) with production rate of 80 m 3 /h: (a) temperature at production well; (b) pressure at production well.
Temperature distribution nephograms of geothermal reservoirs under different situations are shown in Figures 10-12. The results show that the radius of the cold water thermal front is directly proportional to the production rate, and the radius of IGDHS is larger than that of DGDHS. In addition, the radius of the cold water thermal front increases with decreased reinjection temperature. As shown in Figure 11b, the radius of the cold water thermal front is about 260 m. Considering that the uneven distribution of formation and the instability of groundwater accelerate the flow of lowtemperature groundwater, it is recommended that well spacing should be over 300 m in this case. The production parameters of each geothermal heating system are shown in Table 2.      Projects using groundwater for heating can be optimized with the above-mentioned method. The input parameters are adjusted by the geological characteristics of each area and heating modes. The optimal well spacing of the geothermal heating system can be obtained. This methodology provides a good reference for cities developing geothermal heating technology.

Economic Optimization of Geothermal System
In this paper, a static economic evaluation method is used to calculate the initial investment and operation cost of DGDHS and IGDHS. The basic conditions of geothermal systems provided by China Petroleum and Chemical Corporation are listed in Table 3 (1 USD = 7 RMB). Through the multiobjective optimization model, a geothermal system with high investment return and short investment payback period is selected for the case in Xinji.  Figure 13a shows the relationship between well price and production rate. The well price increases with increased production rate, and the same result can be achieved by lowering the reinjection temperature. The reason is that large well spacing is required for a geothermal heating system with large production rates and low reinjection temperatures to prevent thermal breakthrough. There is a great difference in electricity charges and initial investments of various geothermal heating systems. As shown in Figure 14a, the electric energy consumed by IGDHS is nearly 1.5 times that of DGDHS. The power consumption of IGDHS presents an irregular change with increased well production rate. The reason is that the decreased reinjection temperature leads to increased secondary heat transfer and heat load of the heat pump, which causes increased power consumption. However, the trend in water consumption is reversed, shown in Figure 14b. Water consumption of IGDHS is much less than that of DGDHS, reduced by 10%-50%. This is because IGDHS generates heat by electricity and has better environmental adaptability than DGDHS. Even in extreme weather, IGDHS meets residential heating needs. It can be clearly seen from Figure 13b that the initial construction cost of DGDHS is the highest. Although the construction cost of DGDHS decreases with increased production rate, it is still higher than that of IGDHS. For IGDHS, increasing the production rate or reducing the reinjection temperature can reduce the initial construction cost, and the construction investment of geothermal wells will decrease by up to 30%.
There is a great difference in electricity charges and initial investments of various geothermal heating systems. As shown in Figure 14a, the electric energy consumed by IGDHS is nearly 1.5 times that of DGDHS. The power consumption of IGDHS presents an irregular change with increased well production rate. The reason is that the decreased reinjection temperature leads to increased secondary heat transfer and heat load of the heat pump, which causes increased power consumption. However, the trend in water consumption is reversed, shown in Figure 14b. Water consumption of IGDHS is much less than that of DGDHS, reduced by 10%-50%. This is because IGDHS generates heat by electricity and has better environmental adaptability than DGDHS. Even in extreme weather, IGDHS meets residential heating needs. There is a great difference in electricity charges and initial investments of various geothermal heating systems. As shown in Figure 14a, the electric energy consumed by IGDHS is nearly 1.5 times that of DGDHS. The power consumption of IGDHS presents an irregular change with increased well production rate. The reason is that the decreased reinjection temperature leads to increased secondary heat transfer and heat load of the heat pump, which causes increased power consumption. However, the trend in water consumption is reversed, shown in Figure 14b. Water consumption of IGDHS is much less than that of DGDHS, reduced by 10%-50%. This is because IGDHS generates heat by electricity and has better environmental adaptability than DGDHS. Even in extreme weather, IGDHS meets residential heating needs. As shown in Figure 15, DGDHS has a longer payback period than IGDHS. Actually, annual power consumption of DGDHS is lower, because there are no heat pumps, but the high cost of construction investment and equipment depreciation causes a longer payback period. On the contrary, IGDHS consumes large amounts of electricity during operation, but the initial investment for a geothermal heating system is lower, which greatly reduces the payback period. As reinjection temperature is lower, the payback period becomes longer. That is because the annual power consumption of a geothermal heating system with low reinjection temperature is significantly higher than other systems. An indirect geothermal district heating system for building heating in the case of Xinji, China, is the best heating mode. The optimal production rate, reinjection temperature, and well spacing for a 50-year heating period are 100 m 3 /h, 301.15 K, and 300 m, respectively. As shown in Figure 15, DGDHS has a longer payback period than IGDHS. Actually, annual power consumption of DGDHS is lower, because there are no heat pumps, but the high cost of construction investment and equipment depreciation causes a longer payback period. On the contrary, IGDHS consumes large amounts of electricity during operation, but the initial investment for a geothermal heating system is lower, which greatly reduces the payback period. As reinjection temperature is lower, the payback period becomes longer. That is because the annual power consumption of a geothermal heating system with low reinjection temperature is significantly higher than other systems. An indirect geothermal district heating system for building heating in the case of Xinji, China, is the best heating mode. The optimal production rate, reinjection temperature, and well spacing for a 50-year heating period are 100 m 3 /h, 301.15 K, and 300 m, respectively. In fact, the income and operating expenses of geothermal heating systems in Xinji fit nicely with economic calculations. It is shown that the calculation results by the static technical economic evaluation are correct. Due to different hydrological conditions and economic levels of cities, all items in this paper cannot be directly used. However, the mathematical and multiobjective optimization models proposed in this paper can be used for assessment, and calculation results can be referenced for the choice of geothermal heating systems.

Conclusions
Renewable geothermal utilization is a significant approach for residential heating in geothermal enrichment regions. The development of geothermal reservoirs has to meet building heating demands with optimal economic benefits, and the geothermal water should be protected. In this paper, a new approach with practical engineering requirements and environmental and economic considerations is established systematically. In fact, the income and operating expenses of geothermal heating systems in Xinji fit nicely with economic calculations. It is shown that the calculation results by the static technical economic evaluation are correct. Due to different hydrological conditions and economic levels of cities, all items in this paper cannot be directly used. However, the mathematical and multiobjective optimization models proposed in this paper can be used for assessment, and calculation results can be referenced for the choice of geothermal heating systems.

Conclusions
Renewable geothermal utilization is a significant approach for residential heating in geothermal enrichment regions. The development of geothermal reservoirs has to meet building heating demands with optimal economic benefits, and the geothermal water should be protected. In this paper, a new approach with practical engineering requirements and environmental and economic considerations is established systematically.

1.
Mathematical models of geothermal development with direct and indirect geothermal district heating systems are established. Furthermore, the optimal well number, well spacing, production rate, and reinjection temperature can be obtained to meet engineering, environmental, and economic criteria.

2.
Well spacing, reinjection temperature, and production rate are the most significant parameters affecting thermal breakthrough in geothermal reservoirs. With decreased well spacing and reinjection temperature or increased production rate, premature thermal breakthrough would occur, resulting in low efficiency of geothermal heating systems.

3.
For indirect geothermal district heating systems, the construction investment of geothermal wells is reduced by up to 20%, and annual water consumption is reduced by up to 50%, but electricity consumption costs increase by 5% to 30%. 4.
Well number and reinjection temperature have a huge effect on the payback period of investment for geothermal development. The payback period of direct geothermal district heating systems is longer than that of indirect systems. Indirect systems have a good environmental adaptability and profitability which is more suitable for building heating.

5.
In the case of Xinji, China, an indirect geothermal district heating system is much better than a direct geothermal district heating system, both technically and economically. The optimal production rate, reinjection temperature, and well spacing for 50 years of building heating are 100 m 3 /h, 301.15 K, and 300 m, respectively. This optimal production parameters can be the reference for the design of geothermal heating systems in other regions. The systematical calculation approach can be reasonably applied to the selection and optimization of other geothermal systems.