**Dynamic Simulation and Energy Economic Analysis of a Household Hybrid Ground-Solar-Wind System Using TRNSYS Software**

#### **Rafał Figaj, Maciej Zoł ˛ ˙ adek \* and Wojciech Goryl**

Department of Sustainable Energy Development, Faculty of Energy and Fuels, AGH University of Science and Technology, 30059 Krakow, Poland; figaj@agh.edu.pl (R.F.); wgoryl@agh.edu.pl (W.G.)

**\*** Correspondence: mzoladek@agh.edu.pl

Received: 18 June 2020; Accepted: 5 July 2020; Published: 8 July 2020

**Abstract:** The adoption of micro-scale renewable energy systems in the residential sector has started to be increasingly diffused in recent years. Among the possible systems, ground heat exchangers coupled with reversible heat pumps are an interesting solution for providing space heating and cooling to households. In this context, a possible hybridization of this technology with other renewable sources may lead to significant benefits in terms of energy performance and reduction of the dependency on conventional energy sources. However, the investigation of hybrid systems is not frequently addressed in the literature. The present paper presents a technical, energy, and economic analysis of a hybrid ground-solar-wind system, proving space heating/cooling, domestic hot water, and electrical energy for a household. The system includes vertical ground heat exchangers, a water–water reversible heat pump, photovoltaic/thermal collectors, and a wind turbine. The system with the building is modeled and dynamically simulated in the Transient System Simulation (TRNSYS) software. Daily dynamic operation of the system and the monthly and yearly results are analyzed. In addition, a parametric analysis is performed varying the solar field area and wind turbine power. The yearly results point out that the hybrid system, compared to a conventional system with natural gas boiler and electrical chiller, allows one to reduce the consumption of primary energy of 66.6%, and the production of electrical energy matches 68.6% of the user demand on a yearly basis. On the other hand, the economic results show that that system is not competitive with the conventional solution, because the simple pay back period is 21.6 years, due to the cost of the system components.

**Keywords:** household hybrid system; dynamic simulation; economic analysis; TRNSYS software

#### **1. Introduction**

Over the last few years, the main goal of worldwide energy policy is decreasing the use of fossil fuels, and, consequently, greenhouse gas emissions. Many changes are taking place to make possible the achievement of such a goal, especially in the form of policy of incentives for investments in renewable energy, scientific grants for actions in developing innovative energy systems, or requirements for energy-efficient buildings [1]. In this context, positive effects in the energy sectors may be achieved. In fact, due to policy, research, and incentive strategies, European total energy consumption in the last decade decreased by about 10%, to 1105 Mtoe in 2017. In general, all the energy-consumption sectors, such as the transport, industry and residential, are involved in the developed strategies. In the framework of energy use, about 25% of European final energy demand, 283 Mtoe, is used in the residential sector, especially for heating and cooling purposes [1].

Further decreasing demand in this sector is possible due to the development of energy-efficient buildings, especially nearly zero-energy buildings (NZEB) [2–4]. The most common energy sources for such buildings are photovoltaic (PV) and photovoltaic-thermal collector (PVT) systems [5–7], wind turbines [8,9] and heat pumps [10,11]. Nevertheless, efforts in reducing energy consumption lead also to the development of highly efficient hybrid and polygeneration systems [12].

Geothermal energy technologies, involving earth-air heat exchangers or electric heat pumps (EHP), can be coupled with other renewable energy devices to improve the efficiency of the whole system [11]. Li et al. [13] presented a parametric study of a standalone installation with a 7 kW wind turbine, a set of batteries, and a heat pump for a single-family house in Sweden. Dynamic simulations of such a system proved that wind energy, due to its intermittent characteristics, cannot fully satisfy the needs of the heat pump.

Roselli et al. [14] provided an analysis of ground source heat pump supported by a 5 kW wind turbine with battery storage. The model was prepared for 200 m<sup>2</sup> office and analyzed for two locations in Italy—Cagliari and Naples. The dynamic model developed in the Transient System Simulation (TRNSYS) software showed that the coefficient of performance (COP) and energy efficiency ratio (EER) for Naples were consequently 3.97 and 4.59. The primary reduction per kWh of final energy demand was 0.8 for Naples and 1.24 for Cagliari. The use of batteries allowed to lower the value of energy exported to the grid in the range of 27% to 63%, depending on the capacity. The fraction of the energy demand met by renewables was about 25% for Naples and 48% for Cagliari.

PV-based systems also are unsuitable to fully cover the electrical energy needs of heat pumps. Kemmler [15] provided a simulation of 4 reference buildings in Germany with a PV-EHP system, controlled by an algorithm that was maximizing the self-consumption of PV energy. The authors proved that 25.3–41.0% of the electricity demand may be covered with a PV system. This value increases in the case of adding batteries, but such a solution increases the payback time.

Psimopoulos et al. [16] evaluated the impact of advanced control strategies on energy and economic performance of a residential heat pump system with a photovoltaic field and electrical energy storage operating under a wide range of climate conditions by means of TRNSYS software. Presented results show that forecast-based control systems lead to greater final energy savings. The presented strategy allowed to reduce final energy use up to 842 kWh per year, which leads to about 175 EUR savings per year in Spain and Germany and about half of this value for Italy, France, and the United Kingdom.

Another trend in the hybridization of EHP systems consists in connecting them with phase-changing material thermal storage. The study reported in ref. [17] showed that in Italy, due to the use of PCM storage it was possible to reduce electric energy consumption by 18%, thanks to improved COP in terms of heating (increase from 3.5 to 4.13) and cooling (increase from 4.0 to 5.9).

Another interesting concept that may be used with EHP installation is the installation of PVT collectors. Such systems provide electricity for operating the heat pump and also generating heat with high efficiency, allowing to lower the electrical energy consumption of ECH needed to match the thermal demand. Research presented in reference [18] showed that the installation of PVT collectors in Italy leads to a reduction of the non-renewable primary energy requirements of buildings, nonetheless, it is connected with higher investment costs.

Most of the literature studies describe EHP driven by renewable energy provided only by one source, like wind turbines or PV fields. However, those energy resources are mostly available in different seasons—solar energy abounds in summer months, and wind energy may be available during the whole year or seasonally, depending on local wind resource conditions. Connecting those technologies in one hybrid installation may provide benefits in terms of matching a significant part of electrical energy demand with renewables.

Ozgener provided an experimental study of a solar-assisted geothermal heat pump connected with a 1.5 kW wind turbine to meet the heating demand of a greenhouse in Izmir, Turkey. Only about 3% of energy consumption was provided by the turbine. The author concluded that such a system could be economically viable in areas with higher wind resources, and the results of the presented study could be extended to the residential sector [19].

Vanhoudt et al. [20] described an installation based on EHP for a residential building located in Belgium. The system consisted of 7.7 kW of PV panels and a 5.8 kW wind turbine. Authors proved that a dynamic control strategy allows us to decrease the one percent peak power demand by up to 5%, and the average power demand about 17%. The self-consumption of renewable energy increased by about 20%. Due to frequent switching, energy consumption increased by 8–12%.

Stanek et al. [21] proposed an energy and thermo-ecological analysis to evaluate the impact of a 3.5 kW PV field and small (3 kW) wind turbine on the operation of EHP in a 163 m2 residential building in Katowice, Poland. Authors considered a combination of operation of both systems, a wind-turbine (WT) system only and three configurations of PV fields. The average exergy efficiency of wind-turbine and PV panels was equal to 25.9% and 12.94% respectively, while the value for system power plants was equal to 33.67%. Average thermo-ecological costs for EHP installation supported by a wind turbine and grid, PV field and grid and wind turbine, PV field, and grid were, consequently, 2.14, 2.16, 1.24.

Li et al. [22] presented an analysis of a residential building HP system with a 5 kW wind turbine and flat-plate solar collectors. Collectors were used for heating of domestic hot water and increasing heat pump evaporation temperature for room heating. Wind power was used for meeting the heat pump power demand. Wind turbines provided about 7.5% of the yearly power demand of heat pump to satisfy thermal load of 198 m2 residential building located in Beijing.

Rivoire et al. [23] analyzed the application of a ground-coupled heat pump (GCHP) in different buildings and climates. The highest profitability of such systems was found for poorly insulated buildings in cold climates, because of its high demands (8.6–9.9 years). However, a balance between heating and cooling needs was found in temperate climate zones for highly insulated buildings. This parameter is important for a ground heat pump (GHP), because of reducing the thermal imbalance of the ground. Authors proved that HP with a capacity of 60% of peak load meets 82–96% of the annual demand. The reduction of primary energy consumption was 33–75% and CO2 emission was reduced by 27–56%.

A numerical study taken for the tropical environment in Thailand [24] proved, that using a ground-source heating pump (GSHP) only for cooling purposes allows us to achieve about 40% energy savings. Analysis was performed for three types of administration buildings (about 9 working hours per day, 200, 40, and 25 m2). Electricity consumption for air conditioning in such buildings was 29,835, 4874, 3372 kWh/year respectively. The authors concluded, that GSHP may be feasible in Thailand's condition, however, there are problems with installing the required amount of borehole heat exchangers. As a solution, the authors proposed installing them in car parks.

Another study for Thailand described the installation of GSHP with horizontal heat exchangers [25]. The authors compared experimental data from two months of operation for the GSHP and air-source heat pump (ASHP). It was found that the GSHP consumed about 19% less electricity than the ASHP system. Comparing those two systems it was proved that the reduction of CO2 emission in the case of GSHP was about 3000 kg. Because of the high initial cost, such installation was found to be unprofitable, however, the growth of electricity costs may significantly increase this parameter.

The literature analysis shows that existing works are mainly dedicated to the description and investigation of systems with one renewable energy source connected to electric ground source heat pumps. Applications with hybridization of energy sources in GHP systems are not so numerous. This paper presents the numerical analysis of a hybrid HP system driven by renewable electrical energy provided by PVT collectors and a wind turbine. To the authors' best knowledge, this is the first time in literature that a system like the one considered in this paper is comprehensively investigated. In fact, there are no papers available in literature dealing with the proposed system configuration, detailed simulation of the system operation and adoption of realistic user demand. The presented control strategy leads to a maximum reduction of non-renewable primary energy requirements. Analysis was provided for a year's operation of the system, considering the energy and economic aspects of the presented system. Furthermore, a parametric analysis of the system is performed, where solar field extension and wind turbine power are varied to investigate the effect of these parameters on the system performance.

#### **2. Methodology**

The proposed hybrid geothermal-solar-wind system was modelled and dynamically simulated adopting TRNSYS software [26]. The selected tool is widely used in commercial and scientific applications to perform advanced multipurpose analyses of complex, hybrid and novel energy systems based on close to the reality simulations [27,28]. The software comprises a vast built-in library of components, spacing between different categories of equipment (heating cooling and ventilation, solar devices, etc.), models of which are validated experimentally and/or are intrinsically validated because they are based on manufacturer data [26]. Therefore, the simulation of systems developed within TRNSYS software carries out reliable results.

In the software used, selected components and connections among them were used to develop the system structure in a user-friendly graphical environment. In particular, the modelling and the simulation of the system were performed using both built-in library components (e.g., ground heat exchanger, photovoltaic thermal collectors, wind turbine, etc.) and user-defined components (control system, energy and economic model). As regards the building, Type 56 was used to simulate its thermal behavior, while the SketchUp TRNSYS3d plug-in [29] was used to implement the geometric structure of the building (thermal model).

The list of the built-in components used to develop the model of the system is shown in Table 1. In this section, only a brief introduction of the main system components models was provided, since the detailed presentation of the models is available in the software reference [26]. The model of the fan-coil operating in heating model was taken from the literature [30], as for the heat pump, manufacturer data (Aermec WRL 026/161 [31]) was used based on the approach reported in [32].


**Table 1.** List of the built-in components used to develop the hybrid system.

#### *2.1. Layout and Operation Strategy of the System*

The proposed system includes vertical ground heat exchangers coupled with a reversible heat pump for space heating and cooling purposes, and a photovoltaic/thermal solar collector field assisting the ground heat exchanger in the heating operation during winter and preheating the domestic hot water all year long. The solar field and the micro wind turbine are used for the generation of electrical energy, matching in part the demand of the user.

The basic layout of the system, including all the main components and loops, is shown in Figure 1. The loops of the system were designed to manage adequately the thermal energy flows in order to match the thermal energy demand of the user. The loops consist of:



**Figure 1.** Layout of the system with main loops and components.

The devices that were integrated in the system layout are the following:


The system layout included also several devices in order to manage properly the fluid loops, as diverters (D), mixers (M), and single (P2, P3, P4, and P5) and variable speed (P1) pumps. It is worth noting that the layout presented in Figure 1, is a simplified version of the system implemented in the simulation environment, since this one includes also other components, like controllers, weather database and components generating the results of the simulation as plotters, variable integrators and printers. The controller components were coupled together in order to develop an efficient control strategy of the system, described as follows.

In the heating season, the ground heat exchanger heated the working fluid provided by the bottom of the storage tank. The heated fluid was supplied to TK1 by means of a stratification supply system that allows one to supply the fluid to the tank to the node with the temperature closest to the entering fluid. In this way, the fluid with a higher temperature was supplied to the top part of the tank while the fluid with a lower temperature was supplied to the bottom part, ensuring thermal stratification of the storage. From top of TK1, the mixture of glycol and water was supplied to the source side of the reversible heat pump, which operated in order to keep the temperature of the buffer tank TK3 within a fixed dead band. RHP was activated when the temperature inside the tank drops to 35 ◦C while it was turned off when the temperature increased to 38 ◦C, in order to avoid an unnecessarily heating of TK2. Moreover, during the heating operation, the heating power of RHP was modulated by means of a proportional controller as a function of the tank temperature. The control strategy was developed to set the power of RHP to 100% when a temperature of 35 ◦C inside the TK3 tank was reached, while it was set to 20% for a temperature of 37 ◦C. It is worth noting that the activation of the pump dedicated to GHX was activated only when RHP was turned on in order to avoid an unnecessarily electrical energy consumption, and this control strategy was performed during both heating and cooling periods.

RHP ensured a proper operation of TK3 in terms of temperature required by the fan-coil system for the space heating, with the last one activated in order to ensure a temperature of the air inside the building rooms between 20 and 22 ◦C.

During the cooling season, the RHP was activated in a chilling mode in order the keep the TK3 temperature between 10 and 7 ◦C, and similarly in the winter period the unit cooling power was modulated proportionally the TK3 temperature. In fact, the load factor of RHP was set to 100% when a temperature of 10 ◦C inside the TK3 tank was reached, while it was set to 20% for a temperature of 8 ◦C. During the operation of RHP in cooling mode, the heat rejected by the condenser of the unit was supplied by means of GF to the top of TK1, and from there was dissipated by GHX.

The operation of the loops at the load and source side of TK1 were managed by a diverter and mixing valve system, consisting of D3–D6 diverters and M3–M6 mixers. In particular, the valves were used to supply GHX with the fluid stored in the bottom part of TK1 (lowest temperature) and to supply RHP with the fluid stored in the top part of TK1 (highest temperature) during winter. The same system was used in summer in order to supply from the bottom part of TK1 tank the fluid entering the condenser of RHP and to supply GHX with the fluid stored in the top part. The adoption of this valve system allowed one to operate both GHX and RHP with in a proper way from the level of temperature point of view.

The photovoltaic/thermal collectors were used all year long to heat the water stored in the bottom part of TK3, in order to produce DHW, while during winter they were used also to heat TK1 aiding the space heating operation performed by RHP. The pump of the solar collectors was activated when the solar radiation increases above 100 W/m2 in order to avoid a possible cooling of the fluid flowing within the collectors. Moreover, in order to avoid a possible cooling of TK1, during winter, and TK2, all year long, a by-pass consisting of D1 and M1 valves was adopted. On the basis of which the tank was supplied by solar collector, the control strategy checked continuously that the fluid temperature exiting the collectors was at least 2 ◦C higher than the one inside the tanks before supplying it to the storage. If such a condition was not met, SF was recirculated within the collectors and M1 and D1 valves. In detail, in this control strategy, the solar collectors' outlet temperature was compared to the top temperature of TK1 and to the water temperature inside TK2 at the level of the internal

heat exchanger supplied by the solar loop. The setpoint temperature of the solar collectors was set by means of the operation of the variable speed pump and depended on the tank to be loaded: 30 ◦C for TK1 and 60 ◦C for TK2.

The heating operation of TK2 by the solar system was set to be performed until the water near the bottom internal heat exchanger of the tank reached the setpoint value of 55 ◦C. In case the DHW demand is not met by solar thermal energy, GB is activated. In particular, then the TK2 top temperature drops to 50 ◦C GB is turned on in order to heat the water up to 55 ◦C. The auxiliary operation of GB allows one to match the user DHW independently from the availability of solar energy.

The solar loop operation during winter was designed to heat primarily TK3 and secondarily TK1. The heating of TK1 was allowed only when the temperature of the fluid exiting the solar collectors was not enough high to heat the water near the lower internal exchanger of TK2 (SF temperature 2 ◦C higher than that of the water). Once the temperature of the fluid exiting PVT increases above the threshold relative to the temperature of TK2, the heating of the DHW is allowed.

Finally, the electrical energy produced by the PVT and WT is supplied to the user or to the grid depending on the demand. The grid allows one to store virtually the energy produced in excess and use it part when the demand overcomes the production.

#### *2.2. Ground Heat Exchanger*

The component has been implemented using the Type557 model, simulating a vertical heat exchanger that transfers heat with the ground. The adopted model allowed one to simulate both U-tube ground heat exchanger or a concentric tube ground heat exchanger by means of a subroutine. Nevertheless, for the purposes of this study, the U-tube sub-model was adopted.

The component simulated both heating and cooling operation of the heat exchanger, since the working fluid may absorb heat from or reject it to the ground as a function of the temperatures of the fluid and the ground surrounding the heat exchanger.

The adopted model was based on the following assumptions:


The detailed description of the model is here omitted fake of brevity. Its comprehensive presentation is available in reference [33].

#### *2.3. Phototoltaic*/*Thermal Collectors*

The solar field is equipped with flat-plate photovoltaic/thermal solar collectors, consisting of a conventional solar thermal device with an absorber covered by a PV cell layer encapsulated within a transparent protective layer, typical of PV panels. The model of the PVT collector was based on the flat-plate thermal collector mathematical model modified in order to take into account the presence of the PV cell panel. In detail, the TRNSYS model Type50 was used for the simulation of the solar field, where constant values for the overall energy loss coefficient, the glass transmittance, and the absorbance of the absorber were assumed. The model was based on several equations, and the most important ones have been presented below.

The modified Hottel–Willier–Bliss method [34] and an energy balance was adopted in order to calculate the PV cell temperature (*Tcell*). The overall thermal loss coefficient of the collector per unit area (*Ul*\*), overall collector heat removal efficiency factor (*FR*) and the PVT thermal power output (*QPVT*) were calculated with the following equations:

$$
\mathcal{U}I\_1^\* = \mathcal{U}I\_1 - \pi\_\mathcal{g} I\_{\text{tot}} \eta\_{PVT} \left( T\_{\text{cell}} - T\_{\text{amb}} \right) \tag{1}
$$

$$FR = m\_f c\_f \frac{1 - e^{\frac{f \eta L\_f^\* A\_{PVT}}{m\_f c\_f}}}{\mathcal{U}\_l^\* A\_{PVT}} \tag{2}$$

$$Q\_{PVT} = A\_{PVT} I\_{tot} FR \pi\_{\mathcal{g}} \alpha\_{PVT} \mathcal{U}\_l^\* \left( T\_{f,in} - T\_{amb} \right) \tag{3}$$

where *Itot* is the solar radiation, τ*<sup>g</sup>* is the glass transmission coefficient, *fp* is the efficiency factor, *mf* is the mass flow rate, *cf* is the fluid specific heat, *APVT* is the area of PVT, α*PVT* is the PVT absorption coefficient *Tf,in* is the fluid inlet temperature and *Tamb* is the ambient temperature. The PVT electrical efficiency (η*PVT*) was calculated on the basis of the temperature of the photovoltaic cell:

$$
\eta\_{PVT} = \eta\_{PVT,ref} \left[ 1 - \beta\_{PVT} \left( T\_{cell} - T\_{ref} \right) \right] \tag{4}
$$

where η*PVT,ref* is the reference PVT electrical efficiency, β*PVT* is the efficiency temperature coefficient, *Tcell* is the photovoltaic cell temperature and *Tref* is the reference temperature. In addition, energy balances of the collector were used to calculate the electric power and the outlet temperature of the fluid.

#### *2.4. Wind Turbine*

The built-in library TRNSYS component Type90 was used to simulate the micro-wind turbine. The adopted model calculated the power produced (*PWT*) as a function of air density (ρ), power coefficient (*cp*), rotor area (*Ar*) and wind speed (*v*), as follows:

$$P\_{\rm WT} = \frac{1}{2} \rho\_{\rm air} c\_p A\_{\rm r} v^3 \tag{5}$$

where *cp* is dependent on the axial induction factor (*a*) according to the following equation:

$$c\_p = 4a(1-a)^2\tag{6}$$

When *a* reaches the value of 1/3, *cp* reaches the maximum theoretical value of 0.593, according to the Betz limit [35].

The effect of the variation of the height above the ground of the site of WT installation was taken into account by the model in terms of the air density variation and the increase of wind speed. The density of air as a function of the height (*z*) was calculated with the ideal gas law where temperature (*Tz*) and pressure (*pz*) are considered:

$$
\rho\_z = \frac{p\_z}{RT\_z} \tag{7}
$$

*Tz* is determined on the basis of the vertical gradient of temperature ("lapse rate" [36]), as follows:

$$T\_z = T\_0 - B\,z\tag{8}$$

where *T*<sup>0</sup> was assumed constant to 15 ◦C and *B* was equal to 6.5 ◦C/km [37] Moreover, the pressure *pz* was calculated as follows:

$$p\_z = p\_0 \left(1 - \frac{Bz}{T\_0}\right)^{5.26} \tag{9}$$

where the reference pressure (*p*0) is provided on the basis of dynamic local conditions. The wind speed (*v*) variation with the elevation was evaluated on the basis of Von Karman analysis [38]:

$$\frac{v\_1}{v\_2} = \left(\frac{z\_1}{z\_2}\right)^a \tag{10}$$

where α is the wind shear exponent, dependent on the site location and its dynamic conditions like atmospheric conditions, presence of objects obstructing the airflow, roughness of the ground surface, etc. In this study, the shear exponent was fixed to 0.14, a value that characterized ideal boundary layer conditions. During the simulation, Equation (10) was used to calculate the wind speed at a certain elevation taking into account time-dependent wind velocities provided at ground level by the Meteonorm weather database included in TRNSYS.

Besides the mathematical formulation of the model, the component in order to calculate the power produced by WT required an external file, where the geometry of the wind turbine as well as the characteristic curve of power as a function of the wind speed were provided. In order to implement such manufacturer data, the data from the datasheet of an ENAIR 30PRO unit [39] were considered. In order to scale the wind-power characteristic, the manufacturer data was normalized with respect to the nominal power.

#### *2.5. Energy and Economic Model*

In order to assess the global energy and economic performance of the hybrid system (HS), a comparison was performed assuming a conventional system (CS) accordingly to the methodology reported in reference [27,28]. This approach was adopted since without the adoption of a reference, the energy and economic performance of the hybrid could only have been evaluated in absolute terms, with are not exhaustive for a comprehensive analysis of the system. Using this approach, it was assumed that both HS and CS must supply the same amount of final energy to the user, in terms of space heating and cooling, domestic hot water and electrical energy.

In detail, the reference system consists of:


HS energy performance was evaluated on the basis of a primary energy (*PE*) consumption comparison on a yearly basis with respect to CS. *PE* consumption of CS was calculated as follows:

$$PE\_{\rm CS} = \frac{E\_{\rm th,heating} + E\_{\rm th,DHIV}}{\eta\_{\rm GB,ref}} + \frac{E\_{\rm th,cooling}}{\rm COP\_{\rm cool,ref} \eta\_{\rm el,ref}} + \frac{E\_{\rm el,user}}{\eta\_{\rm el,ref}} \tag{11}$$

where *Eth,heating* and *Eth,cooling* was the yearly thermal energy provided for space heating and cooling, respectively, and *Eel,user* was the yearly electrical energy demand of the user.

The equation used to calculate the *PE* consumption of HS was:

$$PE\_{HS} = \frac{E\_{th,GB,DHW}}{\eta\_{GB,ref}} + \frac{E\_{el,grid,net} - E\_{el,mem,overall,grid}}{\eta\_{el,ref}} \tag{12}$$

where *EthGB,DHW* is the auxiliary thermal energy supplied by GB for the production of DHW and *Eel,grid,net* is the energy electrical energy provided by the grid, excluding the electrical energy supplied to and recovered from the grid (net metering) and the *Eel,unrecovered,grid* is the electrical energy supplied to the grid and not recovered by the user. On the basis of primary energy consumption, the primary energy saving ratio (*PESr*) was evaluated as follows:

$$PESr = \frac{PE\_{CS} - PE\_{HS}}{PE\_{CS}} \tag{13}$$

The economic performance of HS was evaluated taking into account its investment costs and the operating costs of both HS and CS. In particular, the literature and market cost data were used to determine the capital cost of proposed system components, in accordance with the procedure adopted in reference [28]. In particular, the following cost function were adopted:

$$\mathcal{L}\_{\rm GHX} = 20 \cdot L\_{\rm GHX} \tag{14}$$

$$C\_{PVT} = 400 \cdot A\_{PVT} \tag{15}$$

$$\mathcal{C}\_{WT} = \mathcal{c}\_{WT} P\_{WT} = 4298.75 P\_{WT}^{-0.141} \cdot P\_{WT} \tag{16}$$

$$\mathcal{L}\_{\rm RHP} = 4.7108 Q\_{\rm muon,heat,RHP}^2 + 139.69 Q\_{\rm muon,heat,RHP} + 3845.7 \tag{17}$$

$$C\_{TK1} = 494.9 + 808V\_{TK1} \tag{18}$$

where *LGHX* is the length of the vertical ground heat exchangers, *APVT* is the area of the PVT collector field, *PWT* is the nominal power of WT, *Qnom,heat,RHP* is the nominal thermal power of RHP and *VTK1* is the volume of TK1. It is worth noting that the costs of the other tanks (TK2 and TK3) were not taken into account since the devices were assumed to be already present in CS.

The total cost of HS was evaluated summing the cost evaluated by Equations (13)–(17) and adding 10% in order to take into account for the balance of plant components like valves, pumps, etc.

The operation costs for both HS and CS were calculated assuming constant electrical energy tariff (*cel*) of 0.133 €/kWh and natural gas price (*CNG*) of 0.50 €/m3 characterized by a reference low heating value (*LHVNG*) of 9.59 kWh/m3 [40,42].

For HS, a net metering system, thus the possibility to store and utilize on demand the electrical energy produced in excess, was implemented by means of a bidirectional connection with the grid. In particular, it is assumed that in case of a system with nominal power up to 10 kW it is possible to recover freely up to 80% of the energy supplied to the grid, while the eventual quote exceeding 80% is paid with a mean price for electrical energy. In the case of a nominal power between 10 and 40 kW, the limit is lowered to 70%. This net metering method is adopted in Poland in order to increase the share of small-scale renewable systems in the electrical energy market [43].

Therefore, the operation cost of HS (*Cop,HS*) and CS (*Cop,CS*) and the savings (Δ*Cop,HS*) were calculated as follows:

$$\mathbf{C}\_{\rm op,HS} = \frac{E\_{\rm th,GB,DHW}}{\eta\_{\rm GB,ref} LHV\_{\rmNG}} \mathbf{c}\_{\rm NG} + \left( E\_{\rm cl,grid,net} + E\_{\rm cl,vacuum,lim,grid} \right) \mathbf{c}\_{\rm el} \tag{19}$$

$$\mathbf{C}\_{\rm op,CS} = \frac{E\_{\rm th,h\,heating} + E\_{\rm th,DHW}}{\eta\_{\rm GB,ref} LHV\_{\rmNG}} \mathbf{c}\_{\rm NG} + \left(\frac{E\_{\rm th,cooling}}{COP\_{\rm cool,ref}} + E\_{\rm el,user}\right) \mathbf{c}\_{\rm el} \tag{20}$$

$$
\Delta C\_{op,HS} = C\_{op,CS} - C\_{op,HS} \tag{21}
$$

where *Eel,recovered,grid* was the electrical energy produced by the system and recovered from the grid above the limit established by the net metering prosumer contract (70 or 80% of the energy sent to the grid)

Finally, the simple pay back (SPB) index was introduced in the economic model in order to evaluate in a simple way the performance of the system in terms of profitability.

#### **3. Case Study**

The proposed system was investigated under a case study consisting of a single-family household with a ground floor, an attic and a sloped roof. This building as case study was used in a previous paper of the authors [28]. However, here it was adopted under different climatic conditions from those used in the previous paper, and with a dynamic electrical energy demand of the user and the simulation of the heating load. The building structure is reported in Figure 2.

The ground floor area was 100 m2 with a floor height of 2.70 m, while the attic had the same height as the ground floor and a useful floor area of 75 m2. The roof slope was 30◦. The ground floor consisted of two rooms with an area of 25 m<sup>2</sup> and a room of 50 m2. The components of the building envelope, like walls, roof, and floor, were implemented assuming several series of layers with a global transmittance reported in Table 2. The detailed structure of the envelope elements has not been reported here for reasons of brevity.

**Figure 2.** Structure of the case study building.



The climatic conditions of Gdansk, Northern Poland, were selected in order to simulate both the system and building operation. For this purpose, the Meteonorm weather database was used, implementing in the simulation a TMY2 file where the mean yearly weather conditions determined on the basis of a 10-year long period were stored. The mean hourly air temperature, horizontal solar radiation, and wind speed at ground level are shown in Figure 3.

**Figure 3.** The mean hourly air temperature, horizontal solar radiation, and wind speed at ground level for Gdansk locality.

The air-conditioning system of the building integrates an independent fan coil unit for each zone of the house providing space heating and cooling. The space heating period was set from 15 November to 31 March, and the cooling one from 1 May to 15 October. The heating operation was allowed 24 h per day in order to keep the room temperature between 20 and 22 ◦C, while space cooling was performed between 8:00 am to 10:00 pm in order to maintain the air temperature between 24 and 26 ◦C.

The model of the building was completed by adding typical loads of a residential house, in order to simulate reliable thermal gains and losses. In particular, the following loads were considered: human activity (5 house inhabitants—sensible heat 75 W, latent heat 75 W), equipment (3.3 W/m2), lights (5.0 W/m2), fresh air infiltration (0.25 Vol/h). In addition, a dynamic DHW demand profile was implemented assuming that each person consumes 60 dm<sup>3</sup> of hot water consumption a day.

The electrical energy demand of the user was normalized with respect to the daily demand, and it is reported in Figure 4 for different days and periods of the year. The profiles were implemented for different days, workdays, Saturdays, Sundays/Holidays, for two periods in the year, period 1 and period 2 assumed from November to March and from April to October, respectively. In particular, the profiles were developed on the basis of standard electrical energy consumption data available for users as the one here investigated [45] performing a normalization of the hourly data on the basis of daily energy consumption (sum of the hourly consumption). Therefore, the profiles developed were characterized by an integral over 24 h equal to 1.0, thus, in this way it was possible to scale the normalized profiles in order to achieve a fixed annual electrical energy consumption for the user. The yearly electrical energy consumption of the building was set to 5.0 MWh, according to common energy demand levels for this kind of utility [45].

**Figure 4.** Normalized daily electrical energy demand of the user for different year periods and types of day.

The main design and operating parameters of the proposed system components are reported in Table 3. Such parameters were selected on the basis of manufacturer data in order to allows the system a proper operation in terms of the capability of energy transfer, number of activation-deactivation cycles of the components, and to match the user thermal demand.



**Table 3.** *Cont.*

#### **4. Results and Discussion**

The dynamic simulation of the system was carried out for a one-year basis, from 0 to 8760 h. with a 0.04 h time step, calculating temperature and power trends and integrated variables. Under these simulation parameters, a large amount of dynamic data was generated by the simulation, thus for reasons of brevity, only the most important results were reported in this paper. In detail, the dynamic operation of the system is shown for a typical winter operation day in terms of temperature and power trends, while the behavior of the system over the yearly operation was presented on monthly basis. Finally, the yearly results are presented showing the global energy and economic performance. The system is also investigated changing the area of the PVT collector and the nominal power of the wind turbine, remaining the other parameters unchanged. In fact, the dimensioning of the ground heat pump system was performed on the basis of space heating and cooling demand of the user, which is constant under the adopted assumptions.

#### *4.1. Daily Operation of the System in Winter*

The dynamic operation of the system was analyzed for all the days of year, nonetheless due to the high amount of data this analysis cannot be presented here; thus, only the results for a winter operation day have been reported. The selected day was 18 February, falling between 1176th and 1200th hour of the year, which presented the full operation characteristic of the system operation from the point of view devices activation.

The most important temperatures of the system have been shown in Figures 5 and 6. During the firsts hours of the day, both inlet and outlet temperature of GHX (GHX, in and GHX, out) slowly decreases because the heat pump is in operation and this determines a decrease of the temperature at the bottom of TK1, which supplies GHX. The decrease is also due to the fact that thermal energy produced by the PVT during the previous day slightly increased the temperature of the GHX loop. From 0:00 am to about 8:00 am the heat pump is supplied with GF at a temperature about 3 ◦C (TK1, out), while the temperature exiting the evaporator of RHP (RHP, source, out) oscillates slightly above 0 ◦C, due to the supplied heat to run the heat pump. After 8:00 am, PVT solar collectors start to supply heat to TK1 and the heat pump starts to reduce the heat output, determining an increase of its top temperature and as a consequence an increase of the temperatures in the GHX-RHP loop. The temperature at the outlet reaches a temperature of about 9 ◦C just after 12:00 am, and after then starts to decrease because TK2 starts to be heated with solar thermal energy. However, when this operation stops, the heating of TK1 is again performed but the temperature keeps to decreasing due to the operation of RHP.

It is worth noting that during the central hours of the day the temperature at the outlet temperature of RHP condenser slightly decreases while the heated water (HCW) at the outlet of TK3 (TK3, out) increases. This situation is achieved since the heat pump thermal output decreases as a function of a reduction of the space heating load of the building.

**Figure 5.** Main temperatures of GHX, TK1, RHP and TK3, daily analysis, winter day.

As concerns the PVT and TK2 loos, the results point out a decrease of the inlet temperature of PVT (PVT, in) in the night hours due to the heat losses on the pipes connected to the solar field. At the same time, the water temperature at the top of TK2 (TK2, out) decreases being the heat transferred to the lower parts of the tank where the temperature is slightly lower. In fact, the temperature of the water inside the tank near the solar heat exchanger (TK2, HX1) increases during the night by about 2 ◦C. After 8:00, am the temperature at the outlet of PVT collectors (PVT, out) increases above the inlet one allowing to start the heating of TK1. However, this operation is performed only when the outlet temperature is 2 ◦C higher than the top of TK1, according to the developed control strategy. The TK1 heating operation ends after 12:00 am when the PVT outlet temperature overcomes by 2 ◦C the temperature of the water inside TK2 near HX1, which is the condition that triggers the heating of TK2. The control system sets the PVT setpoint temperature to 60 ◦C from the previous value of 30 ◦C, thus the outlet temperature increases due to the reduction of the mass flow rate of SF operated by P1 and the proportional controller. In the result, the solar heat increases the water temperature at the bottom of TK2 to about 28 ◦C, allowing the preheating of DHW, and this temperature level is maintained until occurs the next DHW request near 7:00 pm. During the selected day the heating operation of DHW is completed by GB, which maintains the temperature at the top of TK2 (TK2, out) between 48 and 55 ◦C.

**Figure 6.** Main temperatures of PVT and TK2, daily analysis, winter day.

The thermal and electrical powers for the selected day are reported in Figures 7 and 8, respectively. The thermal power of the ground heat exchanger (GHX) varies as a function of the temperature of the tank (see Figure 5) and the thermal power supplied by the heat pump to the user (RHP, load). The thermal demand for space heating is higher, and the heat extracted from the ground is higher. During the day, the reduction of the heating demand along with the heating of TK1 by the solar field implies a decrease of the thermal extracted by GHX to zero just after 12:00 am, since the increase of the temperature of the working fluid of GHX is not possible due to the thermal conditions of the ground coupled with a relatively high bottom temperature of TK1. The solar field thermal power (PVT) increases up to 4 kW at noon, then the heat starts to be transferred to TK2 with lower power (TK2, HX1). This is due to a higher water temperature inside TK2 with respect to that present at the bottom of TK1, which limits the solar heat transfer rate to the tank. The heating of TK2 by renewable thermal energy is not adequate to match the DHW demand, thus the activation of the GB and the heating of the top part of TK2 by HX2 is mandatory (TK2, HX2).

**Figure 7.** Main thermal powers of the system, daily analysis, winter day.

The electrical power trends highlight that during the night hours the output produced by the wind turbine (WT) is lower than the total power demand of the user, consisting of the building equipment, heat pump, and system auxiliaries operation. Therefore, a significant part (about 1/3) of the user demand (user) must be matched by the energy supplied by the grid. This mainly consists of the grid

energy that is not taken into account by the net-metering system (net grid) since the contribution of the energy recovered by the grid (from the grid) is marginal. Only for about the first hour of the day, the user demand is matched by the energy that was produced during the previous day and supplied to the grid.

The system operates in the central hour of the day denotes that a certain amount of energy produced by the system is supplied to the grid (to grid), due to the output of the PVT field. The energy virtually stored by the grid is immediately supplied again to the user once the power generated by the system (PVT + WT) decreases below the demand just before 3:00 pm.

**Figure 8.** Main electrical powers of the system, daily analysis, winter day.

#### *4.2. Monthly Results*

The results were analysed by performing an integration of the powers on a monthly basis, thus the oscillations typical of the dynamic system operation were agglomerated in the integral results. The main thermal energy flows in the system are shown in Figure 9. Here, the effect of the user space conditioning demand on the system behavior is clearly shown, since the thermal energies flows inside the system are significantly higher during the winter heating season compared to the cooling one.

**Figure 9.** Main thermal energies of the system, monthly analysis.

The energy extracted by the ground heat exchanger (GHX) rages between 0.78 MWh to 1.41 MWh assuming a trend dependent on the thermal energy supplied by the load side of RHP (RHP, load). It is worth noting that as expected the major part of the thermal demand of the RHP evaporator is provided by GHX (at least 75.3% for all the months), while a relatively small part is provided by PVT (between 6.5% and 26.5%). From the point of view of summer operation, the heat rejected by GHX is from 11 to 37% higher than the one transferred by the condenser of the heat pump (RHP, source), and this is due to the fact that the pump of GHX (P3) operates independently from that of RHP (P4).

Moreover, it is important to note that the monthly thermal energy produced by PVT (PVT) outside the heating season is relatively small compared to the solar energy available (solar). This occurs because the solar field area is oversized with respect to DHW demand present during summer, and this determines a decrease in the thermal performance of the solar collector.

The monthly electrical energy flows are shown in Figure 10. The solar field electrical energy production (PVT) following the trend of the solar energy availability reported in Figure 9 ranges between 0.048 and 0.30 MWh. The energy production of the PVT field is lower than the one achieved for the wind turbine (WT) during winter, while the opposite occurs during summer. In particular, the system receives from WT between 0.11 and 0.40 MWh of electrical energy which is used to match in part the user demand.

The results also point out that the operation of the heat pump for space heating significantly affects the user's electrical energy demand during winter, since in summer the thermal energy removed by RHP from the user in the form of space cooling is relatively small, as outlined in Figure 9. Under these conditions, the production of electrical energy during the year by both PVT and WT determines the that during winter the amount of energy withdrawn from the grid excluding the net-metering (net grid) is relatively high, being between 41.1% and 58.8% of the demand (user). Conversely, it is null during the summer due to favorable conditions of demand and energy production. In addition to this, it noticeable that the energy produced in excess by the system (to grid) is always above zero, and increases significantly during the central months of the year, as well as the energy recovered from the grid in the frame of the net-metering contract. However, is it interesting to observe that the electrical energy stored virtually in the grid during April and May can be recovered in October and November when the PVT energy production drops.

**Figure 10.** Main electrical energies of the system, monthly analysis.

The main efficiency parameters of the hybrid system are reported in Figure 11. As preliminarily mentioned, the thermal efficiency of the PVT field (η*, th, PVT*) in winter is relatively higher than the one achieved in summer, since the temperature of the solar loop increases when only DHW is produced. In fact, the aperture area of the collector allows one to match the major part of the DHW demand in the summer months (more than 70% of the demand). It is worth noting that it is not possible to achieve a higher value due to the time of the DHW demand, occurring in the morning and late evening hours when the solar energy availability is relatively lower.

The decrease of the PVT thermal performance also negatively affects the electrical efficiency (η*, el, PVT*), being the PV cell affected by the operation temperature of the solar collectors. The efficiency decreases of about 15.9% with respect to the maximum value of 0.177 achieved in January and December. The wind turbine performance is presented in Figure 11 as a function of the normalized equivalent number of operation hours (H, WT) defined as the ratio between the equivalent hours of the wind turbine and the total number of hours in the considered time interval. As expected, the performance of WT decreases in summer, operating in terms of equivalent hours less than 30% of the time due to a lower availability of the wind energy source characteristic of the selected location.

Finally, the performance of the heat pump in terms of COP for space heating (COP, heat) and cooling (COP, cool) highlights that the device operates under stable conditions in terms of source and load temperatures. The COP in winter varies less than 1.5% compared to the maximum value of 4.57, while during summer the variation is 2.8% with respect to the maximum value of 8.17.

**Figure 11.** Energy efficiency parameters of the system, monthly analysis.

#### *4.3. Global Energy and Economic Results*

As done for the monthly analysis, the yearly results are analysed performing the integration of the variables from 0 to 8760 h. The main thermal and electrical energies of the system in the yearly time scale are reported in Table 4. Analysis of the results shows that despite a significant availability of solar energy (*Esol,PVT*) the production of thermal energy by the PVT field (*Eth,PVT*) is less more than 50% of the thermal energy produced by GHX in total (*Eth,GHX,heat* and *Eth,GHX,cool*). As shown by the monthly analysis, the ground heat exchanger is significantly more exploited in winter compared to summer due to different space conditioning demands, as outlined by the fact that the winter energy extracted 6.4 times higher than the heat rejected in summer.

It is interesting to note that the yearly DHW demand (*Eth,DHW*) is matched by solar energy only in 38.7%, and this occurs because the contribution of the solar energy to the production of DHW (*Eth,HX1*) during winter is almost null, while during summer it not exceeds 75.3%.

From the point of view of electrical energies, among the generation devices, WT produces more energy. In fact, the yield of PVT (*Eel,PVT*) is 0.87 MWh less than the one achieved by WT (*Eel,WT*), due to the comparatively better availability of wind energy compared to the solar one for the selected locality. The results also show that the incidence of the system auxiliary devices (*Eel,AUX*) on the system electrical energy balance is negligible, the energy demand of such equipment being only 5.6% of the total amount (*Eel,demand*), while the yearly energy demand of the heat pump is significant, since its operation determines 25.4% of the annual energy bill.

Moreover, the analysis of the energy flows among the system and the grid points out that the electrical energy supplied to the grid (*Eel,to grid*) and recovered (*Eel,from grid*) is the same, and thus the full bidirectional net-metering balance is achieved. In particular, the user sends to the grid 24.5% of the electrical energy produced, while the grid supplies the user with net-metering-free energy (*Eel,net grid*) in order to cover the demand in 31.4%.


**Table 4.** Main thermal and electrical energies of the system in MWh, yearly results.

The energy and economic indexes of the system are presented in Table 5. The savings of HS in terms of primary energy with respect to CS (*PESr*) are remarkable (66.6%), nonetheless, some primary energy must be consumed by HS due to the grid intervention to meet the electrical energy demand and to the operation of GB, especially during winter.

The scarce thermal efficiency of the PVT field affects the total efficiency of the collector units (33.1%), while the wind conditions of the selected locality allows one to achieve a good performance of WT, the normalized number of equivalent operation hour being equal to 0.333. This value underlines that the amount of electrical energy produced by WT is one-third of the maximum achievable for such a unit under ideal operational wind conditions (wind speed always above 11 m/s).

The yearly performance of the heat pump for heating and cooling overcomes the nominal conditions, indeed both COPs are higher than the values reported in Table 3. Such a result is achieved because: (i) during the winter the operation of GHX and PVT allows the evaporator of RHW to operate at a relatively higher temperature, and the outlet condenser temperature is lower than the nominal one (45 ◦C), (ii) in summer GHX allows one decrease the condenser temperature below the nominal operation conditions (30 ◦C).

On the basis of the comparison in terms of operational cost between CS and HS, the proposed innovative system allows one to save 0.900 k€/year which is 65.4% of the operational cost of the conventional system. Nevertheless, taking into account the fact that the investment for the whole hybrid system is relatively high (almost 20 k€), the economic performance of the investigated solution is scare. SPB is slightly over 21 years, which implies a feasibility of the system outside reliable criteria for economic investments. However, assuming that both PVT and WT are financed by 70% from incentive policies, SPB value decreases to 11.6 years, which is a value that starts to be economically viable. It is worth noting, that such incentive strategies, based on capital incentives, are increasingly common across Europe, and thus its applicability is reasonable. Moreover, it is worth noting that the proposed economic analysis of the system is performed under the worst scenario, where that total cost of the investment is considered. In case of the necessity of the installation of a new system providing DHW, heating and cooling, or in case of the substitution of the existing one, the cost that should be considered in the investigation would be only the difference of cost between CS and HS. Under such conditions the profitability of the system would be better compared to the present analysis. An interesting possibility of incentive for the proposed system could be that based on the consumption of produced renewable thermal and electrical energy, thus an additional saving may be generated by the hybrid system. In the invested case study, if an incentive of 0.05 €/kWh is adopted for both thermal

(heating and cooling) and electrical energy produced and consumed by the user, SBP decreases to 13.1 years.


**Table 5.** Main energy and economic parameters of the system, yearly results.

#### *4.4. Parametric Analysis*

The energy and economic performance of the system was also investigated changing the extension of the PVT field and the nominal power of WT. This analysis was performed with the aim of determining the effect of such design parameters on the global results when all the other ones are assumed to be constant. The PVT area was varied from 3.0 to 15.0 m2, with a step to 3 m2, whereas WT power was increased from 0.3 to 1.5 kW, with a step of 0.3 kW.

The results of the parametric analysis as a function of the PVT field area is presented from the point of view of energy in Figure 12, and from that of the efficiency and economical viewpoint in Figure 13.

**Figure 12.** Main thermal and electrical energies of the system vs. variation of the PVT area, parametric analysis.

The variation of PVT area affects almost linearly the electrical energy producible by the solar field, thus the negative effect of a mean higher solar loop temperature achievable for a higher collector extension is limited. As a consequence, the electrical energy supplied by the grid and the net of the net-metering contract decreases linearly. The opposite occurs for the thermal energy produced by PVT collectors, increasing less when the area is gradually augmented. In particular, the increase of the thermal energy produced by the field configuration of 3 to 6 m<sup>2</sup> is 0.606 MWh, whereas from 12 to 15 m2 it is only 0.288 MWh. The increase of the production of the thermal energy by PVT implies a decrease of the thermal energy extracted by GHX, due to a higher mean temperature of TK1, and an increase of the heat supplied to TK2 by means of HX1. Nevertheless, the yield of GHX is scarcely affected by such variation, because the additional thermal input of the PVT field to TK1 achieved in case of a higher area is marginal, and a similar situation is achieved for HX1. In this case, the solar field reaches the almost the maximum capability of providing heat to TK2, due to the DHW demand and dynamics of solar thermal energy production during the day, for and PVT area above 6 m2. Additionally, it can be observed that the increase of the collector area does not affect the thermal energy used by the load side of RHP (evaporator) to match the used space heating demand, the latter being constant. In addition to this, also the COP of RHP in heating model remains practically constant, showing that the contribution of the solar field to the heating of TK1 is negligible.

The thermal efficiency of PVT decreases more than two times when the area increases from 3 to 15 m2, due to the previously discussed reasons, highlighting that the increase of the PVT field extension significantly affects the operation of the solar system in terms of thermal efficiency. On the other hand, a limited decrease of the electrical efficiency of PVT is achieved along the range of investigated area values (only 5.1%), because the increase electrical losses of PVT due to a higher PV cells temperature is marginal.

The savings of primary energy achievable by HS increase as a function of the PVT area on the basis of a higher thermal and electrical energy produced and used in the system. However, the main reason for the increase of *PESr* is due to the production of electrical energy rather than the production of heat, being the effect of the last one on the energy balance of the system relatively marginal. Therefore, in terms of primary energy savings, it is convenient to increase the PVT area, and the same condition is achieved from the economic viewpoint, since *SPB* decreases for a higher area. However, the decrease is not sufficient to achieve the system economic profitability because the value remains over 21 years.

**Figure 13.** Efficiency, primary energy and economic parameters of the system vs. variation of the PVT area, parametric analysis.

The effect of WT power on the system performance is shown in Figures 14 and 15.

The capacity of WT affects only the electrical energy flows in the system as the thermal and electrical parts of the system are decoupled. The increase of WT nominal power, keeping constant the hub height and the wind conditions, determines a proportional increase of the generated electrical energy, which is equal to 2.92 MWh per kW. In view of this trend, the electrical energy supplied to the grid and withdrawn from the grid increases as well for a larger wind turbine, although the increase, in this case, is not linear due to the characteristic curve of the device. In fact, the increase of power from 0.3 to 0.6 kW implies a variation of energy supplied to the grid of 0.196 MWh, although the increase of 0.544 MWh for the variation from 1.2 to 1.5 kW. Moreover, it is important to note that beyond 1.2 kW of WT power, it is not possible to recover all the electrical energy supplied to the grid, since due to the intrinsic characteristic of the user demand and dynamic production of electrical energy by PVT and WT. Therefore, even when a turbine with a larger diameter is installed, the electrical energy needed

from the grid is not null for the WT capacities considered. In practice, the increase of the wind turbine yield per unit of power determines the same reduction of the electrical energy produced within the grid and required by the user is needed to match the demand.

Finally, the analysis of primary energy savings shows that the trend of energy production achieved increasing the WT capacity directly affects *PESr*, and this is achieved because a higher production of electrical energy means a higher saving of non-renewable energy provided by the grid to the user. In the investigated case, the numerical data shows that increasing WT power 5 times implies an increase of *PESr* of 82.3% compared to the value achieved for a 0.3 kW unit.

**Figure 14.** Main electrical energies of the system vs. variation of the WT power, parametric analysis.

As concerns the economic savings of HS with respect to CS, they increase of 67.0% passing from a wind turbine with 0.3 kW to a one with 1.5 kW, and on the other side in the same range of power, the total cost of the system (*C, tot*) increases of 30.6%. Therefore, coupling these two economic conditions, it is clear why the trend of SPB is decreasing as a function of WT power. In particular, *SPB* trend passed from 25.4 years for a 0.3 kW unit to 19.8 years for a power 5 times higher, consequently, the effect of higher savings overcomes the effect of the system cost increase when larger wind turbines are considered.

**Figure 15.** Primary Energy and economic parameters of the system vs. variation of the WT power, parametric analysis.

#### **5. Conclusions**

In the paper, a novel micro-scale hybrid geothermal-solar-wind system for a residential user has been investigated through dynamic simulation performed using a complex model developed in TRNSYS software [26]. The layout of the system included vertical ground heat exchangers, photovoltaic/thermal solar collectors, wind turbine, net-metered bidirectional connection with the grid, and a water-to-water reversible heat pump. The system was used to match a part of the user electrical energy demand, and to provide space heating and cooling along with DHW. The dynamic operation of the system and the building was simulated simultaneously in order to take into account their mutual interaction. In order to investigate the system, a case study of a single-family house under Polish northern climatic conditions, and a conventional reference system for the user, were considered.

The operation of the system was investigated from the viewpoints of temperature and thermal and electrical power trends for a selected winter day. Monthly and yearly analyses were carried out in order to assess the energy and economic system performance during the year and in terms of global performance. The study was concluded performing a parametric analysis of the system, where the solar collector field area and micro wind turbine power were varied.

The main results of the study were the following:


The present study outlined that the hybrid system under consideration achieves a satisfactory performance in terms of energy, whereas its economic profitability is not adequate to permit its application as a valid alternative compared to an existing reference system based on a natural gas boiler from domestic hot water production and heating, an electrical chiller for space cooling and the grid for matching user demand. However, considering the installation of a new system matching the user demand, or in the case of the substitution of the existing one, only the extra cost for the hybrid system with respect to the reference one must to be considered. Thus, the profitability of the hybrid system will increase and the gap between the reference system in terms of economic competitiveness will be significantly reduced.

The investigations of the proposed system will be expanded by performing analyses regarding the modification of the system layout, comprehensive parametric analysis of the system, and optimizations aiming at determining the characteristics of the system depending on the design parameters and user type and behavior.

**Author Contributions:** Conceptualization, R.F. and M.Z.; methodology, R.F.; software, R.F.; formal analysis, M. ˙ Z.; ˙ investigation, R.F.; data curation, R.F.; writing—original draft preparation, R.F. and M.Z; writing—review and ˙ editing, R.F. and W.G.; visualization, R.F.; supervision, R.F. and W.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** Part of this work was funded by the Polish Ministry of Higher Education on the basis of the decision number 0086/DIA/2019/48.

**Acknowledgments:** This work was carried out under Subvention and Subvention for Young Scientists, Faculty of Energy and Fuels, AGH University of Science and Technology, Krakow, Poland. The authors acknowledge for the use of the infrastructure of the Center of Energy, AGH UST in Krakow.

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

#### **Nomenclature**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Geothermal Power Production from Abandoned Oil Reservoirs Using In Situ Combustion Technology**

#### **Yuhao Zhu 1, Kewen Li 1,2,3,\*, Changwei Liu 4,\* and Mahlalela Bhekumuzi Mgijimi <sup>1</sup>**


Received: 11 October 2019; Accepted: 20 November 2019; Published: 24 November 2019

**Abstract:** Development of geothermal resources on abandoned oil reservoirs is considered environmentally friendly. This method could reduce the rate of energy consumption from oil fields. In this study, the feasibility of geothermal energy recovery based on a deep borehole heat exchanger modified from abandoned oil reservoirs using in situ combustion technology is investigated. This system could produce a large amount of heat compensated by in situ combustion in oil reservoir without directly contacting the formation fluid and affecting the oil production. A coupling strategy between the heat exchange system and the oil reservoir was developed to help avoid the high computational cost while ensuring computational accuracy. Several computational scenarios were performed, and results were obtained and analyzed. The computational results showed that an optimal water injection velocity of 0.06 m/s provides a highest outlet temperature of (165.8 ◦C) and the greatest power output of (164.6 kW) for a single well in all the performed scenarios. Based on the findings of this study, a geothermal energy production system associated with in situ combustion is proposed, specifically for economic reasons, because it can rapidly shorten the payback period of the upfront costs. Modeling was also performed, and based on the modeling data, the proposed technology has a very short payback period of about 4.5 years and a final cumulative net cash flow of about \$4.94 million. In conclusion, the present study demonstrates that utilizing geothermal resources or thermal energy in oilfields by adopting in situ combustion technology for enhanced oil recovery is of great significance and has great economic benefits.

**Keywords:** geothermal power production; abandoned oil reservoirs; in situ combustion

#### **1. Introduction**

Generating geothermal power from abandoned oil reservoirs could be the answer to the energy crisis facing today's world. However, the upfront cost in a geothermal project is a barrier. Increasing the heat production and shortening the payback period of the upfront costs are two main objectives in a geothermal project. High outlet temperature and larger amount of geothermal fluid extraction could help to enhance the total heat production, whereas low drilling costs are conducive to reduce the payback period of the geothermal project. In this regard, most oilfields possess great potential for producing geothermal energy, because they are associated with production of geothermal fluid from oil reservoirs. Furthermore, abandoned oil reservoirs are good potential areas for harvesting and producing geothermal energy, and are most favorable because they are easily retrofitted into geothermal wells. However, the outlet temperature of geothermal wells in oil fields is usually lower, thus resulting in low heat production and long payback period.

Over recent years, the petroleum industry has seen and experienced rapid development due to the high demand of energy worldwide. In the midst of this rapid development, a larger number of oil and gas wells has been left abandoned due to either technical problems or economic benefit limitations, particularly heavy oil reservoirs. To enhance heavy oil recovery, thermal recovery technology such as in situ combustion method is often used. According to the latest update on drilled oil wells data, there exist about 20–30 million abandoned oil wells around the world [1]. Abandoned oil reservoirs may not hold oil productivity after a long-term recovery process, but such wells can produce brine fluids with high temperature. Thus, geothermal energy from such high-temperature fluids could be directly or indirectly recovered back to the surface for further utilization. Li et al. [2] and Zhang et al. [3] proposed a geothermal system combined with power generation which extracts geothermal heat from hot water in oilfields. Zheng et al. [4] proposed a concept to produce geothermal energy from abandoned oil and gas reservoirs by oxidizing the residual oil with injected air. This system, however, lacks a detailed heat exchanger model describing the heat transfer from reservoir to geothermal well, which is the essential part of evaluating the efficiency of a geothermal system. Zheng et al. [4] proposed a concept which lacks a detailed heat exchanger model, Davis and Michaelides [5] investigated the feasibility of the heat exchanger using abandoned wells for power generation by ORC (Organic Rankine Cycle). Their proposed geothermal system can provide an output power of 3.4 MW when at optimal condition.

Wight and Bennett [6] pointed out that the closed loop system (i.e., heat exchanger) retrofitted from abandoned oil well has the advantage of protecting both the equipment and the environment from performance or safety risks caused by minerals and contaminants commonly present in naturally occurring geothermal fluid. This is because the system can produce geothermal energy without extracting geothermal fluid from the formation. Another advantage of utilizing abandoned oil reservoirs is that there is little or no drilling cost at all, which is generally the main upfront investment in any geothermal project [7]. Moreover, the cost of retrofitting abandoned oil well into geothermal well using work-over techniques is equivalent to one tenth the cost of drilling a new geothermal well, and the maximum expense is even less than one seventh of the cost of drilling a new well [8]. Therefore, with regard to the above, it is clear that, if abandoned wells could be re-developed and utilized, large amount of drilling cost could be cut down, and more economic gains could be achieved.

Numerous mathematical models have been developed to both calculate the outlet temperature and evaluate the performance of the heat exchanger. Bu et al. [8] considered a two-dimensional model for numerically simulating the heat transfer between the surrounding rocks and the heat exchanger. This model neglects the variation of temperature in the tube wall of the injection annulus and the extraction well, which decreases the accuracy of the temperature distribution and overestimates the heat extracted from the surrounding rocks. Furthermore, Nian and Cheng [9] pointed out that it is unsound to use the assumption of vanishing well radius as a boundary condition in the model, because this could overestimate the heat flux from the surrounding rocks to working fluid in the heat exchanger.

Noorollahiet et al. [10] simulated a geothermal system which extracted geothermal energy by using two abandoned oil wells in Ahwaz oil field in Southern Iran. They developed a 3D numerical model using finite element method software ANSYS. They found that the sensitive parameter analysis indicates that casing geometry will strongly affect the heat transfer between the geothermal well and the surrounding rocks. Mokhtari et al. [11] found that when evaluating the pressure drop and thermal efficiency of the heat exchanger, the diameter ratio of the inner to outer pipe is very important. Their results showed that the optimal diameter ratios for pressure drop minimization and thermal efficiency maximization are 0.675 and 0.353, respectively.

In general, there are three grades of temperatures for extractive geothermal energy recognized by the industry namely: high temperature (>180 ◦C), intermediate temperature (100 to 180 ◦C) and low temperature (30 to 100 ◦C), respectively [12]. High-temperature geothermal resources can be used for power generation, while the low and medium temperature resources could be/are normally used directly for district heating and geothermal source heat pump (GSHP) [13,14]. However, a large number of geothermal reservoirs have low outlet temperatures. Low-temperature geothermal resources are economically not viable for commercial plants because they cannot generate sufficient heat. Thus, they cannot recover the upfront investment for retrofitting the abandoned wells and the construction of the power plants. To extract more heat from subsurface reservoir, low boiling point organic fluids with high heat-to-electricity efficiency are often used as working fluid mostly in ORC plants. Mokhtari et al. [11] investigated four organic fluids (R123, R134a, R245fa, R22) often used as working fluids in geothermal power plants. They found that R123 has the most desirable characteristics in comparison with other working fluids. Noorollahi et al. [15] further investigated two other types of working fluid (isobutene and ammonia) for power generation. A binary power cycle model was applied for power generation in their study to calculate the total electric power output. The two fluids (ammonia and isobutene) were chosen as secondary fluids in their power cycle, and their power outputs were compared. They found that for similar flow rates, the total net power for isobutene is greater than that of ammonia.

Another method for enhancing geothermal system is extracting geothermal energy from oil reservoirs with in situ combustion. Zheng et al. [4] proposed a geothermal system enhanced by oil reservoir using in situ combustion technology. They used a retrofitted abandoned well for their geothermal system. Similarly, Cinar [16] proposed a geothermal model using a perforated well to extract heated fluid from the oil reservoir where wet in situ combustion was applied. However, as pointed out in his study, wet in situ combustion may extinguish the combustion process, which causes reduction in the mass production of heated water resulting into less sufficient heated water to sustain a commercial scale power plant. Additionally, this may also result in the large depletion of the formation fluid. Cheng et al. [17] investigated the enhancement effect of geothermal power generation from abandoned oil reservoirs with thermal reservoirs. They found that a geothermal well with thermal reservoirs could produce a heat and electric power output of about four times the amount that can be produced by a geothermal well without thermal reservoirs. Although the cost of drilling can be reduced or even eliminated, the cost of the power plant installation accounts for the major portion of the total investment of the geothermal project, which further prolongs the payback period. Therefore, in this regard, more heat extraction from the subsurface is needed to sustain a commercial-sized geothermal project/plant.

In this study, a new model for recovering and re-utilizing oil for geothermal energy from abandoned wells is proposed. For this model, an in situ combustion technique is applied in the reservoir to recover heavy oil by injecting air into the formations. The air injection helps to oxidize the oil and heat up the reservoir. A schematic representation of the newly proposed model can be seen in Figure 1.

Figure 1 is a schematic representative summary of the development and utilization of the resources in abandoned oil reservoirs by in situ combustion. As can be seen from the figure, during the recovery process, a large amount of heat is generated from the in situ combustion. The abandoned well becomes a heat exchanger (i.e., geothermal well). A schematic representation of the heat exchange system is shown in Figure 2.

As can be seen from Figure 2, water is injected through the injection annulus and extracted through the extraction well. The oil reservoir with in situ combustion provides large amounts of heat at the bottom of the geothermal system. The system mainly includes heat transfer in the geothermal well (i.e., extraction well and injection annulus), oil reservoir and the surrounding rocks (i.e., strata). In general, the geothermal fluid extraction process may cause some problems in the formation, including groundwater recession, corrosion and scaling problem, the high cost of geothermal, and drilling of the re-injection well [18]. However, the system shown in Figure 2 has no direct contact in the formation fluid and no effect on oil production. It only extracts the heat from the formation and the oil reservoir by recycling the working fluid in a closed loop concentric tube, which is combined with in situ combustion by air injection in the oil reservoir. The detailed methodology adopted in this study is shown in Figure 3.

**Figure 1.** A new method of geothermal production utilizing abandoned wells and in situ combustion.

**Figure 2.** Schematic of heat exchange system.

In this study, a numerical model describing retrofitting abandoned oil wells into useful geothermal wells with in situ combustion for the purpose of recovering geothermal resources was developed. This model provides a new coupling strategy between the heat exchange system and the oil reservoir. The main idea is that the reservoir temperature with in situ combustion is expressed as a function of time, and then imported into the heat exchange system as a boundary condition. This model takes advantage of the commercial reservoir simulator CMG for simulating in situ combustion and helps to avoid large amounts of computational cost when calculating the heat transfer from the surrounding rocks and oil reservoir to the heat exchanger. Several parameters known to affect the system performance

were analyzed using the proposed model in this study. Moreover, it should be noted that turbulent flow, which increases the computational efficiency without losing much computational accuracy is not considered in this model. For the findings of this study, the computational results indicate that the outlet temperature increases remarkably after the combustion front reaches the area near the heat exchange wellbore. Furthermore, a geothermal system using advanced in situ combustion was proposed to extract more heat from the reservoir right from the beginning, thereby significantly shortening the payback period of the upfront costs.

**Figure 3.** Flowchart followed when investigating the proposed geothermal system.

#### **2. Model Description of In Situ Combustion**

The process of in situ combustion is basically the injection of oxidized gas or oxygen-enriched air to generate heat by burning a portion of residual oil, where most of the oil is driven toward the producers by a combination of gas-drive, steam or water-drive. The main purpose of applying in situ combustion is to generate a large amount of heat for the geothermal system. For this study, the modeling was performed using STARS by CMG. The in situ combustion model assumes that the reservoir has: uniform porosity, isotropic permeability, and closed upper and lower layers. The model utilizes a five spots pattern, as shown in Figure 4.

**Figure 4.** Five spots pattern applied in oil reservoir model.

The numerical model domain, which is the dark area shown in Figure 4, takes advantage of the symmetric flow, with injection in one quarter and a production well in two corners. To simulate the process of in situ combustion in a general form, a direct conversion cracking kinetics scheme for three components of oil was chosen in the simulation, as it does not depend upon the stoichiometry of the products, and thus reduces the degree of uncertainty in the simulation results, as the number of unknowns is reduced [19]. Here we applied the chemical reaction data of the template in STARS of CMG [20] to obtain a typical in situ combustion scenario. Seven components were considered in the process of in situ combustion, and these were: water (H2O), heavy oil component (HC), light oil component (LC), inert gas (IG), oxygen (O2), carbon dioxide (CO2), and coke, respectively. Three chemical reactions were introduced to describe the components change in reservoirs and the heat generated in the process. The enthalpy and activation energy of these three reactions were as shown in Table 1.

**Table 1.** Combustion reactions and their respective kinetic parameters.


Their relative permeability curves and the reservoir properties are shown in Table 2 and Figure 5.


**Table 2.** STARS input parameters for in situ combustion model.

**Figure 5.** Relative permeability curves used in model simulation.

#### **3. Model Description of Geothermal Well**

As seen in Figure 2, an abandoned well can be retrofitted into geothermal well by coating the inner tubing (extraction well) with insulation and sealing the bottom of the well. Water is then used as working fluid, and is injected through the annulus space and extracted from the inner tubing to the surface [21]. When the water flows downward along the annulus, it is heated by the surrounding rocks and by the reservoir, where there is a large amount of heat generated by in situ combustion. The heated water is then extracted back to the surface. This is a concentric tube heat exchanger, and the working fluid is not in direct contact with the surrounding rocks. Based on the data of the main oil fields in China, the geothermal gradient is generally about 0.03 K/m [22–26]. The model developed in this study utilized a concentric tube heat exchanger that was designed to retrofit an abandoned well with a typical casing outer diameter of 19.6 cm and an inner diameter of 15 cm. The inner diameter of the extraction well was 4 cm, the thickness of insulation was 2 cm, the well depth was 4000 m, and the thickness of the reservoir was 10 m. The whole reservoir was fully penetrated by the well. This system was modeled with the finite element modeling software COMSOL Multiphysics.

#### *3.1. Governing Equations*

A two-dimensional axi-symmetric cylindrical model was applied to describe the whole system. Heat exchange takes place between the working fluid and the rocks simultaneously [21,27], see Figure 6.

**Figure 6.** Dimensions of concentric pipe heat exchange system.

The dashed vertical line on the left side of Figure 6 shows the axis of symmetry, which is parallel to the direction of the *z* axis. The total depth of the heat exchange system is 4000 m from the surface. There is a buffer zone of 5 m at the bottom of the heat exchange system. The thickness of the oil reservoir is 10 m. The time-dependent governing equation corresponds to the convection-diffusion equation, which contains additional contributions of heat flux and no other heat source [28]. The heat flux describes the heat transfer from the rock to the injection annulus and from the injection annulus to the extraction well. Therefore, the expressional equation is

$$
\rho \mathbf{C}\_p \frac{\partial T}{\partial t} + \rho \mathbf{C}\_p \overrightarrow{\boldsymbol{\mu}} \cdot \nabla T + \nabla \overrightarrow{q} = 0 \tag{1}
$$

where <sup>ρ</sup> is the density (kg/m3), *Cp* is the specific heat capacity (J/(kg·K)), *<sup>T</sup>* is the absolute temperature (K), <sup>→</sup> *u* is the velocity vector (m/s) and <sup>→</sup> *q* is the heat flux by conduction (W/m2), which can be described using the Fourier's three-dimensional diffusion law

$$
\overrightarrow{q'} = -k\nabla T\tag{2}
$$

where *k* is thermal conductivity (W/(m·K)). Some parameters used in the numerical simulation are listed in Table 3. This heat transfer equation is already built into COMSOL, and can be called up in the software directly.

**Table 3.** The parameters used in the numerical simulation of the concentric tube heat exchanger.


For the water, the values of heat capacity, density and thermal conductivity depend on temperature and are already built into the COMSOL materials database, being expressed by the following formulas:

$$C\_{p\\_water} = 12010.1 - 80.4 \times T + 0.3 \times T^2 - 5.4 \times 10^{-4} \times T^3 + 3.6 \times 10^{-7} \times T^{-4} \tag{3}$$

$$
\rho\_{\text{\\_water}} = 838.5 + 1.4 \times T - 0.003 \times T^2 + 3.7 \times 10^{-7} \times T^3 \tag{4}
$$

$$k\_{\text{\\_water}} = -0.9 + 0.009 \times T - 1.610^{-5} \times T^2 + 8.010^{-9} \times T^3 \tag{5}$$

#### *3.2. Initial Conditions*

The velocity and the initial temperature of the injected water are both uniform along the well, and their values are 0.03 m/s and 30 ◦C, respectively. It should be noted that the velocity of water in the extraction well should be a corresponding value in order to retain a constant mass flow rate in the numerical model. The surface temperature is 15 ◦C. The initial temperature of the rocks is given by the following equation

$$T\_{\mathbb{R},0}(z) = T\_{srf} + \mathbf{G} \cdot \mathbf{z} \tag{6}$$

where *TR,0* is the initial rock temperature (K), *Tsrf* is the temperature at surface (K), *G* is the geothermal gradient (K/m), and *z* is the depth of rock (m).

#### *3.3. Boundary Conditions*

The radius of influence may be at a modest distance from the well over the time frame modeled in this study. Therefore, the approximation that the temperature is constant at a 100-m radius is justified. The chosen distance of 100 m is explained in the following section. The boundary condition of the whole system, including the rock, is given by

$$\begin{array}{c} \left. T\_{R,b} \right|\_{\begin{array}{c} r = R \\ z = Z \end{array}} = \left. T\_{R,0}(z) \right|\_{\begin{array}{c} \\ z = Z \end{array}} \tag{7} \\ \begin{array}{c} \left. T\_{R,0} \right|\_{\begin{array}{c} \\ z = Z \end{array}} \end{array} \tag{7}$$

where *TR,b* is the rock temperature at a constant temperature boundary (K), *R* and *Z* are distances at constant temperature boundary in *r* and *z* direction (m).

However, if in situ combustion is applied in the reservoir, the boundary condition at the depth of the reservoir corresponds to the temperature influenced by the in situ combustion. Hence, the boundary condition at the depth of reservoir is given by

$$\left.T\_{R,b}\right|\_{r=R} = T\_{\dot{\mathbf{x}}}(z,t) \tag{8}$$

where *Tic* is the temperature of reservoir with in situ combustion and it is the function of depth *z* and time *t*. Please note that the distance of the constant temperature boundary is still *R*. This is reasonable, given that the influence of the in situ combustion on the heat exchange system is dominant considering the whole region of the oil field, while the reduction of the reservoir temperature and the heat extraction caused by geothermal well is limited.

#### **4. Results and Discussion**

A two-dimensional axi-symmetric cylindrical geothermic model was applied to simulate the temperature distribution of the geothermal well system as well as the surrounding rock system. To obtain more accurate temperature data, the two systems were coupled by using COMSOL Multiphysics simulator.

#### *4.1. Mesh Independence Study*

Before applying the heat exchanger model to calculate the outlet temperature, a mesh independence study was conducted to obtain stable temperature data. Triangular meshing was applied to the numerical model of the heat exchanger. The outlet temperatures after 5 years of operation are shown in Figure 7. When the number of triangular elements is small, the outlet temperature data oscillate and then begin to converge at certain values as the number of triangular elements becomes larger than 60,000. The numerical model of the heat exchanger in this study had 64,426 triangular elements, which guaranteed the mesh independence while using a relatively small computational time.

**Figure 7.** The outlet temperature at different numbers of triangular elements in the numerical model.

#### *4.2. Validation and Comparison of the Model*

The proposed model was compared and validated with the results from Bu et al. [8] and Templeton et al. [29]. The assumptions made by Bu et al. [8] include neglecting the variation of the tube wall temperature of the injection annulus and the extraction well, and also the Dittus-Boelter relation in the convection model. These assumptions decrease the accuracy of the temperature distribution and will overestimate the heat extracted from the surrounding rocks in their system. Furthermore, the temperature transfer is very sensitive to the properties of the extraction well and the insulation. Again, the Dittus-Boelter relation is not suitable for parts with larger temperature differences, especially near the top of well, and is also not suitable for the annular injection well [30]. For the model used by Templeton et al. [29], the whole simulation domain does not consider the rock below the heat exchanger. It should be noted that both studies neglect the changes in the thermal properties of water, which are a function of temperature.

To achieve fair comparison with the newly proposed model in this study, we applied parameters similar to those applied in the two studies done by Bu et al. [8] and Templeton et al. [29]. To be more precise, the inner diameter of the extraction and the injection wells are 10 cm and 30 cm; the thickness of the insulation is 1 cm; the thickness of casing is 2 cm; the thermal conductivity of the surrounding rocks and the insulation are 2.1 W/(m·K) and 0.027 W/(m·K); the density and the heat capacity of the surrounding rocks are 2730 kg/m3 and 1098 J/kg/K; and the velocity of the injected water and the geothermal gradient are 0.03 m/s and 0.045 K/m, respectively. Figure 8 shows the outlet temperatures calculated by these three models over ten years.

**Figure 8.** Verification and comparison of the results calculated by Bu et al. and Templeton et al. with the results of the proposed model under similar conditions.

As seen from Figure 8, the temperature result calculated by Bu et al. [8] model is about 26% to 44% higher than the results of both the model of Templeton et al. [29] and the proposed model in this study under similar conditions. The overestimation of the model proposed by Bu et al. [8] is because of the assumptions applied in the model, as mentioned in the previous section. The results of both the model of Templeton et al. [29] and the proposed model in this study have a similar shape of trend throughout the operational time/period. However, the negligence of the surrounding rocks below the heat exchange system causes the temperature result of the model of Templeton et al. [29] to be slightly smaller than the result of the proposed model in this study. Figure 9 shows the temperature profile of the injection annulus and the extraction well with 0.03 m/s injection velocity after a period of two months in operation.

**Figure 9.** Comparison of the temperature profile between the proposed model of this study and the model of Bu et al.

As seen from Figure 9, the solid lines and the dashed lines show the temperature profiles of the injection annulus and extraction well, respectively. The highest temperature is found at the bottom of the well at about 4000 m. The variation of the fluid temperature curves of this study (red line) show that there is a slight decrease in the temperature near the top of the injection annulus. The reason for this is that the temperature of the injected water is 30 ◦C, and this is higher than that of the surface temperature, which is 15 ◦C. This causes the injected water to lose some heat in some regions near the wellhead. However, we cannot find any such phenomenon in the results from the model of Bu et al. [8]; thus, this overestimation reduces the accuracy of the simulation results of the model of Bu et al. [8].

#### *4.3. Analysis of Heat Transfer Characteristics*

The parameters applied in the calculation are listed in Tables 2 and 3, and Figures 5 and 6. The temperature distribution in the surrounding rocks and the heat exchange system with in situ combustion are shown in Figure 10.

**Figure 10.** The temperature distribution of the heat exchange system after one year of operation (left figure), and the temperature distribution of the surrounding rocks after 50 years of operation (right figure).

As seen in Figure 10, there is a temperature drawdown near the heat exchange system. The high temperature zone is strongly marked at the depth of the in situ combustion reservoir at about 4000 m. The isothermal contours show that the initial geothermal temperature is not affected by the heat exchange system when the radius is greater than 80 m, even after a period of 50 years of operation. This could help to explain why the 100 m distance (the radius of influence in boundary conditions) from the heat exchange system to the rock in the *r*-direction, and the 500 m distance in the *z*-direction, were chosen as the boundary of the surrounding rocks for the model proposed in this study.

#### *4.4. Heat Production Analysis*

Considering the high upfront investment required for the initial development of the geothermal power plant and retrofitting of the abandoned well, the outlet temperatures of the conventional geothermal well project may not be high enough to generate enough electric energy and recover the upfront costs. Therefore, in order to obtain more heat from the reservoir, in situ combustion was applied to supply the heat extraction from the heat exchange system. A three-dimensional model of the reservoir with in situ combustion was simulated by using the STARS simulator in CMG. The modeling grid numbers are: 20, 20 and 5 in the i, j and k directions, respectively. The temperature distributions in the oil reservoir during the 6th year and 11th year of the total operation period of 50 years are shown in Figure 11. The measurement unit of the temperatures of the reservoir is degrees Celsius.

**Figure 11.** The temperature distribution of the oil reservoir at the 6th year (left figure) and the 11th year (right figure).

In Figure 11, the air is injected from the upper left injector. The lower right producer will be retrofitted into the heat exchange system. The extreme high temperature zone in this figure is where the combustion front is located. The temperature near the wellbore (block 18, 20, 3) and the average temperature of the whole reservoir are shown in Figure 12.

**Figure 12.** The average reservoir temperature (red dashed line) and the temperature near the heat exchange system wellbore (blue solid line).

As seen in Figure 12, when the combustion front approaches the wellbore (heat exchange system), the temperature around the wellbore increases sharply until reaching a peak point at about 350 ◦C. With continuous injection of air, the temperature decreases, and this is followed by a long stage of declination. This indicates that a large amount of heat is generated in the reservoir when in situ combustion is applied. Considerable amounts of heat can be extracted by the heat exchanger system from both the surrounding rocks and the combustion reservoir, especially when the peak point of combustion front has been reached.

The coupling between the oil reservoir model and the heat exchanger model is at the depth of the oil reservoir. If the average reservoir temperature is applied (red dashed line in Figure 12) as the boundary condition of the heat exchanger model, this may strongly underestimate the heat production of the heat exchanger. Therefore, in this study, we applied the temperature profile data of the block (18, 20, 3) into the heat exchanger model as a time-dependent boundary condition at the depth of the oil reservoir. This method reduces the computing cost while also preserving the computational validity of the result. This is reasonable, because as can be seen from Figure 10, the affected region in the surrounding rocks is limited to the place very close to the heat exchange wellbore. After applying the temperature data of the block (18, 20, 3) into the proposed heat exchanger model, the outlet temperature was strongly enhanced by the combustion reservoir. The temperature curves are as shown in Figure 13.

**Figure 13.** The effect of in situ combustion on performance.

It should be noted that in Figure 13, there is a bump up of the outlet temperature when the combustion front reaches the wellbore. The highest outlet temperature of the extraction well is about 150 ◦C. An enhancement effect is strongly marked when the geothermal heat is compensated by the in situ combustion technique. The corresponding temperature profiles of the injection annulus and the extraction well are shown in Figure 14.

In Figure 14, the solid lines and dashed lines represent the temperature profiles of the injection annulus and the extraction well, respectively. The lines with different colors show the temperature profiles of the heat exchange system at different operational times. The temperature bump up is also marked in the temperature profiles at the depth of reservoir (4000 m). As seen in Figure 14, before the combustion front reaches the wellbore (in less than 10 years), the temperature of both the injection annulus and the extraction well at the depth of (4000 m) are almost the same, or at least very similar. However, after the combustion front reaches the wellbore (after 10 years), the temperature of the extraction well is significantly higher than that of the injection annulus. This indicates that the in situ combustion has greatly enhanced the temperature of the water in the extraction well.

**Figure 14.** The temperature profiles of different operational time with in situ combustion.

According to Cheng et al. [31], the power provided by the heat exchange system can be simply given by the following equation:

$$P = M(T\_{out} - T\_{in})C\_p \eta\_{ri} \eta\_{m} \eta\_{\mathcal{S}} / 1000 \tag{9}$$

where *P* is the actual generated power (kW), *M* is the mass flow rate (kg/s), *Tout* is the outlet temperature of extraction well (K), *Tin* is the inlet temperature of injection annulus (K), *Cp* is the specific heat capacity of water (J/(kg·K)), η*ri* is the relative internal efficiency of steam turbine (0.8), η*<sup>m</sup>* is the mechanical efficiency of steam turbine (0.97), η*<sup>g</sup>* is the generator efficiency (0.98) [31]. For this study, the mass flow rate *M* is calculated by multiplying the injection velocity 0.03 m/s and the cross-sectional area of the injection well (see Figure 6).

Figure 15 demonstrates the variation of the outlet temperature and the actual gained power generated by the heat exchange system as functions of thermal conductivity of the insulation, inlet temperature and the velocity of the injected water in 30 years with in situ combustion.

(**a**) Thermal conductivity of the insulation.

(**c**) Velocity of the injected water.

**Figure 15.** The variation of the outlet temperature and the actual power of the heat exchange system on different parameters.

It can be seen from Figure 15a,b that when the thermal conductivity of the insulation decreases and the inlet temperature increases, the outlet temperature increases. The outlet temperature is much more sensitive to the value of the thermal conductivity of the insulation. This indicates that insulation with better thermal resistance properties is essential to reducing the heat loss from the extraction well to the injection annulus, and to obtaining higher outlet temperature. This confirms that it helps to extract more heat from the formation. Although there is a higher outlet temperature due to higher inlet temperature, this does not guarantee a higher actual power. The reason for this is that a high inlet temperature does not provide a large temperature difference between the heat exchange system and the surrounding rocks. This reduces the efficiency of the heat transfer from the surrounding rocks to the heat exchange system. In Figure 15c, because of the increase in the velocity of the injected water, both the outlet temperature and the actual power increased first, and then decreased. This is because the water injection velocity has an optimal value of 0.06 m/s. This explains that there is a large amount of heat loss from the extraction well to the injection annulus when the water injection velocity is small. For situations in which the water injection velocity is larger, the efficiency of the heat extraction from the surrounding rocks to heat exchange system will be very low.

With respect to the temperature near the wellbore of the heat exchange system in the early stage (about 10 years in Figures 13 and 14), it is still not significantly enhanced, even after applying in situ combustion in the reservoir. This means that in situ combustion needs to be applied in advance to avoid the low-temperature stage in order to recover the upfront costs. This is called the advanced in situ combustion method, and comes from a similar concept of advanced water injection technology. In other words, the operation of retrofitting the abandoned well will be executed only when the temperature of the bottom hole increases drastically (i.e., when the combustion front approaches the wellbore).

As can be seen from Figure 16, the outlet temperature of the scenarios without in situ combustion, with in situ combustion, and with advanced in situ combustion are compared. It can be clearly seen that the outlet temperatures of the scenarios without in situ combustion (blue dashed line) and with in situ combustion (red solid line) are relatively low at the beginning. If the heat exchange system starts to operate when the combustion front approaches the area near the wellbore, the outlet temperature will be enhanced significantly at a very early stage, and after that, decrease slowly (green dashed line). To find out how much energy could be extracted by the heat exchange system, the corresponding electricity data were calculated, and are shown in Figure 17.

**Figure 16.** The effect of advanced in situ combustion on performance.

The curves demonstrate that the total electricity generated by the geothermal power plant with compensation for in situ combustion (green dashed line and red solid line) is much more than the electricity generated without in situ combustion (blue dashed line). The cumulative electric energy of the scenario with in situ combustion (red solid line) after a period of 50 years' operation is 50.3 × <sup>10</sup><sup>6</sup> kW·h. The scenario with advanced in situ combustion (green dashed line) has a higher cumulative electric energy of 51.4 <sup>×</sup> <sup>10</sup><sup>6</sup> kW·h.

**Figure 17.** Cumulative electricity generated by the proposed geothermal power system.

#### *4.5. Economic Appraisal*

To determine the duration of the payback period of the upfront investment and the cumulative net cash flow after a period of 50 years' operation, an economic appraisal was performed by considering the following information: the geothermal development engineering and management engineering; the requirements of economic appraisal methods and parameters promulgated by China National Development and Reform Commission; the current fiscal and taxation system and the pricing system in China; and the current status of geothermal energy development.

A typical cost for a vertical geothermal well is about \$2.3 million, based on the data of the geothermal project in Menengai, Kenya [32]. As stated previously, the cost of work over techniques in such wells is equivalent to about one tenth of the total cost of drilling a new geothermal well [33]. Hence, we assumed that the cost of retrofitting an abandoned well would be \$0.23 million for a single well. For the investment of a power plant [34], the cost per installed kW comes to about \$1500/kW (power peak is 159.2 kW during the 13th year in these scenarios). Based on the data given by Yambajan geothermal power generation [35], the electricity sale price is \$0.14 (¥0.93) per kW·h. The management expense is 1% of the annual sales. The corporate income tax rate is 25%. Therefore, for a period of 50 years in operation, the cumulative net cash flow (NCF) curves are shown in Figure 18.

As shown in Figure 18, the scenario with advanced in situ combustion (green dashed line) has the shortest payback period (about 4.5 years) and the largest final cumulative NCF (\$4.94 million). The final cumulative NCF of the two scenarios with in situ combustion (green dashed line and red solid line) are both higher than the scenario without in situ combustion (blue line). Although the scenario with simultaneous in situ combustion (red solid line) has relatively high final cumulative NCF, its payback period is very close to the scenario without in situ combustion (blue dashed line) and it is longer comarable to the scenario with advanced in situ combustion (green dashed line).

**Figure 18.** Cumulative net cash flow curves for the proposed geothermal power system.

In regard to the information presented above, a geothermal project with advanced in situ combustion is recommended. Please note that when evaluating geothermal projects, upfront investment, power output, fiscal policy, etc., should be fully considered.

#### *4.6. Discussion*

Note that this discussion leads to the conclusions made in this study, and was based upon a correlation of the study model with other previous models. The main theme of this paper is the application of in situ combustion for geothermal development with feasible economic gains. Although Bu et al. [8] and Templeton et al. [29] used different assumptions for their models, their purposes were similar to that of this study: the development of an economically friendly project. What makes the model of this study unique is the introduction of in situ combustion, and the fact that this model considers all of the parameters neglected by the other models. For example, the model of Bu et al. [8] neglects temperature variation in the tube wall of the injection annulus and the extraction well, and further neglects the Dittus-Boelter relation in convection model. The model of Templeton et al. [29] does not consider the rock below the heat exchanger. In short, both studies neglect the changes in the thermal properties of water, which are functions of temperature. Such negligence decreases the accuracy of the temperature distribution and may overestimate the heat extracted from the surrounding rocks in the system. However, since temperature is a major key in geothermal projects, the model proposed in this study considered all these aspects, especially temperature transfer in the extraction well and insulation. The results of this model are fair and reliable, because even with different assumptions from the previous models, the parameters applied are similar. These include the inner diameter of the extraction and the injection wells, the thickness of insulation and casing, the thermal conductivity of the surrounding rocks, the velocity of injected water, geothermal gradient, etc.

Furthermore, during the geothermal fluid extraction process, most models may cause some problems in the formation, such as groundwater recession, corrosion and scaling problems, the high cost of geothermal, and the drilling of the re-injection well, but the newly proposed model has no direct contact with the formation fluid and no effect on oil production; thus, it only extracts the heat from the formation and the oil reservoir by recycling the working fluid in a closed-loop concentric tube, which is combined with in situ combustion by air injection in the oil reservoir. The coupling between the oil reservoir model and the heat exchanger model is done at the depth of the oil reservoir. This depth may vary for each reservoir, so as not to underestimate the heat production of the heat exchanger. Again, this phenomenon is neglected in the results of the model of Bu et al. [8], and this reduces the accuracy of the simulation results. Note that the experimental data used in this study are similar to those of the two models described above, but the final results, such as temperature, are all different, which is due to the differences in the parameters that are assumed and neglected between the models. However, the results of this study model have more economic benefits than the other two models. A final point to note is that more attention should be paid to the figures, as more detail is apparent in the figures.

#### **5. Conclusions**

Based on the calculation and analysis conducted in this study, the following conclusions can be drawn:


**Author Contributions:** Conceptualization, C.L. and K.L.; Data curation, C.L.; Formal analysis, Y.Z., C.L. and K.L. Investigation, C.L. Methodology, Y.Z. Software, Y.Z.; Supervision, K.L.; Validation, C.L.; Writing—original draft, Y.Z.; Writing—review & editing, C.L. and K.L.; M.B.M. revised and polished the language of the final version.

**Funding:** This research received no external funding.

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

#### **Nomenclature**


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
