**Preface to "Geothermal Energy Utilization and Technologies 2020"**

The increasing interest in renewable energy sources is leading the geothermal sector towards playing an important role in meeting the final energy demand of different countries. This can help reduce dependency on energy imports and can vary the energy mix. Possible applications can vary from small scale (residential) to large scale (city). Geothermal energy, unlike other renewable sources, has the significant benefit of being a programmable source, therefore it is not dependent on the time of day or on particular weather conditions.

Geothermal energy could be considered for direct use in space heating and cooling, greenhouse heating, aquaculture, bathing, district heating networks and industrial uses. Where specific conditions are available (high-temperature hydrothermal resources, aquifer systems medium temperatures, and hot dry rock), geothermal energy could be converted into electricity by means of different technologies (flash steam, dry steam, binary). The cogeneration and polygeneration systems could be also considered to exploit geothermal sources. A smart energy community being supported by geothermal sources could be a sustainable and interesting option.

In 2017, the global geothermal power generation was 84.8 TWh, while the cumulative capacity reached 14 GW. The capacity is expected to reach 17 GW by 2023 as estimated by the International Energy Agency. Despite its significant potential, the contribution of geothermal energy to global heating, cooling, and power demand is relatively low. This book aims to advance the contribution towards focusing the attention on technologies to guarantee the faster spread of geothermal energy and on the policies supporting the exploitation of this source.

> **Carlo Roselli, Maurizio Sasso** *Editors*

### *Article* **Assessment of Geothermal Resources in the North Jiangsu Basin, East China, Using Monte Carlo Simulation**

**Yibo Wang 1,2, Lijuan Wang 3, Yang Bai 1,2,4, Zhuting Wang 5, Jie Hu 1,2,4,\*, Di Hu 6, Yaqi Wang 1,2,4 and Shengbiao Hu 1,2,4,\***


**Abstract:** Geothermal energy has been recognized as an important clean renewable energy. Accurate assessment of geothermal resources is an essential foundation for their development and utilization. The North Jiangsu Basin (NJB), located in the Lower Yangtze Craton, is shaped like a wedge block of an ancient plate boundary and large-scale carbonate thermal reservoirs are developed in the deep NJB. moreover, the NJB exhibits a high heat flow background because of its extensive extension since the Late mesozoic. In this study, we used the Monte Carlo method to evaluate the geothermal resources of the main reservoir shallower than 10 km in the NJB. Compared with the volumetric method, the Monte Carlo method takes into account the variation mode and uncertainties of the input parameters. The simulation results show that the geothermal resources of the sandstone thermal reservoir in the shallow NJB are very rich, with capacities of (6.6–12) <sup>×</sup> 1020 J (mean 8.6 <sup>×</sup> <sup>10</sup><sup>20</sup> J), (5.1–16) <sup>×</sup> <sup>10</sup><sup>20</sup> J (mean 9.1 <sup>×</sup> <sup>10</sup><sup>20</sup> J), and (3.2–11) <sup>×</sup> <sup>10</sup><sup>20</sup> J (mean 6.6 <sup>×</sup> <sup>10</sup><sup>20</sup> J) for the Yancheng, Sanduo and Dai'nan sandstone reservoir, respectively. In addition, the capacity of the geothermal resource of the carbonate thermal reservoir in the deep NJB is far greater than the former, reaching (9.9–15) <sup>×</sup> <sup>10</sup><sup>21</sup> J (mean 12 <sup>×</sup> 1021 J). The results indicate capacities of a range value of (1.2–1.7) <sup>×</sup> <sup>10</sup><sup>21</sup> J (mean 1.4 <sup>×</sup> <sup>10</sup><sup>22</sup> J) for the whole NJB (<10 km).

**Keywords:** geothermal resource; Monte Carlo simulation; assessment; thermal reservoir; North Jiangsu Basin

#### **1. Introduction**

In the past two centuries, the extensive use of non-renewable energy (e.g., nuclear energy, coal, and so on) in the world has caused many environmental disasters [1–3]. Among them, the consumption of fossil fuels has led to a sharp increase in the concentration of greenhouse gases in the atmosphere, which has led to an increase in global warming. Therefore, green and renewable energy has been highly regarded by more and more countries. Among those types of energies, geothermal energy is now receiving wide attention because of its unique characteristics: wide distribution, resource-richeness, safety and stability, cleanliness and low-carbon nature.

In the four decades since China's reform and opening up, rapid economic growth and the use of fossil fuels have made China the world's largest carbon emitter (more than

**Citation:** Wang, Y.; Wang, L.; Bai, Y.; Wang, Z.; Hu, J.; Hu, D.; Wang, Y.; Hu, S. Assessment of Geothermal Resources in the North Jiangsu Basin, East China, Using Monte Carlo Simulation. *Energies* **2021**, *14*, 259. https://doi.org/10.3390/en14020259


Received: 4 December 2020 Accepted: 29 December 2020 Published: 6 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

100 million tons of CO2 per year). The direct use of geothermal energy has increased rapidly with the increasing demand for energy efficiency in buildings and the reduction of CO2 emissions, and has maintained its position as the world's number one in recent years [1–3]. The North Jiangsu Basin is located in Jiangsu Province, the most developed province in (Eastern) China, and is a typical 'hot in summer and cold in winter' region. Existing detailed geological data, abundant geothermal resources [4], and an important economic basis are an important basis for the efficient utilization of geothermal resources in the North Jiangsu Basin and, at the same time, provide an important possible model for the future use of geothermal resources in the Yangtze River Delta.

The assessment of geothermal resources is the estimation of the amount of geothermal energy which might be extracted from the inner Earth and economically used in the future. A regional resource assessment can, on the one hand, provide a framework for the government or industry with a long-term energy strategy and policy, and on the other hand, help us to ensure rational planning and exploitation of geothermal resources. A geothermal resources assessment has been carried out using the volumetric method in the NJB [5–7] and the Jiangsu Province [8–11]. However, there are currently three main problems: (1) the low amount of geological and geothermal data constrains the accurate geothermal resources assessment; (2) carbonate thermal reservoirs are currently the most important and favorable reservoirs for geothermal resources, yet they have not attracted enough attention from researchers; (3) the volumetric method does not take into account the uncertainty of parameters involved, but rather assigns a specific value to the thermal storage geometry/physical property parameters.

This paper updates and summarizes high-quality temperature profiles measured in geothermal wells, oil wells, and national hydrological observation wells representative of the major structural units in the NJB. We established the whole density and porosity column of the strata (or rocks) in different geologic ages in the NJB. Specific heats of representative rocks from different thermal reservoirs were also determined. Based on these parameters, we evaluated the uncertainties in parameters such as reservoir area and temperature through monte Carlo simulations and finally calculated the geothermal resource potential of the NJB (less than 10 km). We focused our study on considering parameter uncertainties (Monte Carlo method) and the evaluation of geothermal resources in deep carbonate thermal reservoirs. This research can provide reliable primary data for further study of the geothermal resource planning and development, long-term energy strategy policy, environmental protection, etc.

#### **2. Geologic Setting**

The NJB, situated in the northeastern Lower Yangtze Craton, is a mesozoic-Cenozoic sedimentary basin. The basin is bounded by the Su'nan and the Sulu uplifts to the south and north, respectively, and bounded to the west by the Tancheng-Lujiang Fault Zone, with an area of approximately 35,000 km2. Controlled by NE-NNE faults, the NJB can be divided into the Jianhu uplift, the northern Yanfu depression, and the southern Dongtai depression, and the two depressions contain 22 highs and sags (Figure 1).

#### *2.1. Tectonic Evolution and Stratigraphy*

The basement of the NJB mainly consists of three parts, namely metamorphic rocks in the Proterozoic, the Yangtze marine clastic and carbonates in the Paleozoic-Mesozoic, and volcanic rocks forming in the middle Triassic to Early Cretaceous (Qiu et al., 2006). Carbonate rocks in the NJB are mainly found in the Upper Ordovician-Ordovician (Huangxu Formation, Dengying Formation, Hetang Formation, mufushan Formation, Paotaishan Formation, Guanyinshan Formation, Lunshan Formation, Honghuayuan Formation, Dawan Formation, and Tangtou Formation), Carboniferous (Jinling Formation, Laohudong Formation, Huanglong Formation, and Chuanshan Formation) and Permian (Qixia Formation, Gufeng Formation, Longtan Formation, and Dalong Formation) strata. Extensive drilling

and geophysical data show that the sedimentary thickness since Paleozoic exceeds 11,000 m, with Cenozoic strata over 7000 m thick (Figure 2).

**Figure 1.** Geothermal background in the North Jiangsu Basin. (**A**): Schematic geological map of East Asia (modified after Grimmer, et al. [12]). (**B**): Schematic geological structural map of the North Jiangsu Basin (modified after Wang, et al. [4]).

**Figure 2.** Geological interpretation sections in the North Jiangsu Basin ((**a**,**c**) are modified from Wang, et al. [4] (**b**) is revised from oil company reports, (**d**) is modified from Chen [13], Qiao, et al. [14]). The four sections are shown in Figure 1B. Ar, Archean; Pt2, mesoproterozoic; Z-O, Ediacaran to Ordovician; S-P, Silurian to Permian; K, Cretaceous; E, Paleogene; N-Q, Neogene to Quaternary; γ, magmatic rock.

In the Triassic period, the NJB was uplifted due to the collision between North China Craton and Yangtze Craton, and the possible southward squeeze of the Siberian plate [15–17]. The Lower Yangtze craton thus entered the terrestrial sedimentary environment. Since then, the NJB has been primarily characterized by three tectonic evolution episodes: (1) Yizheng movement with fluvial facies and lacustrine facies, resulting in the Cretaceous Taizhou Formation and the Paleocene Fu'ning Formation; (2) Wubao movement with lake delta fluvial-lacustrine facies, forming the Eocene Dai'nan Formation and Sanduo Formation; and (3) Sanduo movement with the Neocene-Quaternary fluvial facies [18–21].

#### *2.2. Geothermal Background*

Thermal history research in the Lower Yangtze Craton suggests that the last tectonothermal event in the NJB probably occurred in the early Palaeogene and that back-arc extension caused by the subduction of the Izanagi plate is likely to be the main reason. The thermal history inversions from the Yancan-1 and An-1 wells show that the peak heat flow value could be around 85 mW·m−2, after which it gradually decreased and stabilized [22]. The present-day high mantle heat flow (35–43 mW·m−2) indicates that the higher geothermal background in the NJB is probably mainly influenced by mantle activity [23].

A heat flow map is an important representation of the geothermal background of a structural unit and is an important basis for regional geothermal resources assessment. Based on previous work [4], we updated and produced an up-to-date heat flow map of the NJB: a total of 78 high-quality heat flow values from within the basin and 165 published heat flow data from around the basin were used (Figure 1). The distribution density

of heat flow sites and the quality of heat flow data in the NJB are extremely high, both nationally and globally. The measured minimum and maximum heat flow values in the NJB are 46 mW·m−<sup>2</sup> (Yandong-3) and 110 mW·m−<sup>2</sup> (site 004), respectively, with an average value of 67 mW·m−2, which is higher than the mean heat flow in Continental China and line with the average global continental value [24,25]. As shown in Figure 1, the middle and eastern Jianhu Uplift, the northern Jinhu sag, and the eastern Dongtai sag in the NJB are relatively high heat flow regions (>70 mW·m−2). The Jianhu Uplift and Huai'an high have the highest geothermal heat flow, with an average geothermal value of over 72 mW·m−2, followed by the Sujiazui high and Xiaohai high, with average geothermal values of 71 mW·m−2. In general, heat flow values in the Jianghu uplift are higher than those in the depressions, and the heat flow in some highs of the depressions is also higher, with some even exceeding that in the Jianhu Uplift.

#### *2.3. Types of Geothermal Reservoir*

Thermal reservoirs in the NJB contains two primary types: Cenozoic sandstone thermal reservoirs, and Ediacaran-Paleozoic carbonate thermal reservoirs. Cenozoic thermal reservoirs can be divided vertically by lithology, porosity, temperature, etc., into the Palaeogene Dai'nan Formation (E2d), the Sanduo Formation (E3s), and the Neogene Yancheng Formation (Ny). The Dai'nan reservoir with its depth varying from 270 to 2800 m is deeper than the Sanduo reservoir (65–1900 m) and Yancheng reservoir (45–1700 m) and therefore has a higher temperature. The large porosity (>20%) of the Cenozoic sandstone thermal reservoirs results in good water-rich conditions and is an important basis for the efficient use of geothermal resources. The distribution of the Cenozoic sandstone thermal reservoirs in the NJB is shown in Figure 3.

**Figure 3.** Distribution of Cenozoic sandstone thermal reservoirs and carbonate thermal reservoirs in the North Jiangsu Basin. (**a**) Yancheng reservoir; (**b**) Sanduo reservoir; (**c**) Dai'nan reservoir; (**d**) Carbonate reservoir.

Numerous geophysical profiles and sedimentological evidence show that the NJB has a large-scale carbonate thermal reservoir at depth, characterized by a wide range and large scale [4,13,14,26–28]. Carbonate rocks in the NJB are much thicker and hotter than the Cenozoic sandstone thermal reservoirs, suggesting a very rich geothermal resource potential. The carbonate reservoirs in the NJB act as medium-high temperature reservoirs and provide optimum conditions for geothermal resource development. Due to the development of possible karst fissures inside the thermal reservoirs, which provide good access for groundwater activity, the heat energy from the deep can be rapidly transferred in a convective manner to the bottom of the sedimentary cover, giving the thermal reservoirs a high temperature. In addition, the unique water-rich nature and permeability of the carbonate karst zone contribute to the efficient extraction and recharge of groundwater. Although there is some variation in the stratigraphic thickness of the Ediacaran-Ordovician rocks, carbonate rocks develop at depth throughout the North Jiangsu Basin, as shown in Figure 3. Considering the accuracy of the profiles and the difficulty of extracting geothermal resources, we only evaluated carbonate thermal reserves above 10 km.

#### **3. Methodology**

#### *3.1. Assessment method*

There are several methods (e.g., surface heat flux method, heat storage modeling, statistical analysis, analogy method) available for geothermal resources assessment, the most widely used being the volumetric method proposed by muffler and Cataldi [29]. The volumetric method calculates the total geothermal energy in the fluids in the rock masses and pores of a study area, i.e., the geothermal energy accumulation. We use the volumetric method in combination with the monte Carlo method to minimize parameter uncertainty. If the porosity of the reservoir and the thermophysical properties of the water and rock are known, then the total geothermal resources can be calculated by the following three equations (Equations (1)–(3)).

$$\mathbf{Q\_{Total}} = \mathbf{Q\_R} + \mathbf{Q\_W} \tag{1}$$

$$\mathbf{Q\_R} = \mathbf{A} \cdot \mathbf{D} \begin{pmatrix} 1 \ -\ \varphi \end{pmatrix} \not\propto \mathbf{c\_R} \ (\mathbf{T} - \mathbf{T\_0}) \tag{2}$$

$$\mathbf{Q}\_W = \mathbf{A} \cdot \mathbf{D} \text{ } \boldsymbol{\varphi} \text{ } \boldsymbol{\varphi}\_W \text{ } \mathbf{c}\_W \text{ (}\mathbf{T} - \mathbf{T}\_0\text{)}\tag{3}$$

where QTotal, QW, and QR represent the total geothermal energy and that stored in pore water and rock, J; A (m2) and D (m) are the area and thickness of the thermal reservoirs; ρ<sup>R</sup> and <sup>ρ</sup><sup>W</sup> denote the density of the geothermal reservoir and geothermal water, kg·m−3; cR and cW are the specific heat of the geothermal reservoir rock and geothermal water, J·kg−1· ◦C<sup>−</sup>1; ϕ is rock porosity, %; T and T0 represent the average temperature of geothermal reservoir and surface, ◦C.

#### *3.2. Monte Carlo Simulation*

The monte Carlo simulation is a numerical simulation method in which probabilistic phenomena are studied as objects. It is a calculation method for inferring unknown quantities of characteristics by taking statistical values from sample surveys. The main steps involved in carrying out a geothermal resource evaluation using the monte Carlo method are: (1) Defining the input parameter set; (2) Setting different distribution models for the input parameters; (3) Defining the relationship between input and output parameters according to the mathematical model; (4) Defining the output parameters; (5) Setting the number of iterations; (6) Simulating and analyzing the results of the calculations and provide a probability distribution function for each result.

In the monte Carlo simulation process, different distribution models can be selected for the input parameters, e.g., uniform, pert, triangular, and lognormal distributions, with each model having its unique frequency distribution. Taking into account the range and pattern of variation of the different parameters, we give different distribution models for the input parameters. The distribution of the measured porosity in the formations suggests that the triangular distribution may be more consistent with the actual porosity pattern; two-dimensional profile and deep borehole data constrain the variation of the thermal reservoir thickness as the triangular distribution; the overall linear increasing of the temperature in thermal reservoirs indicates that a triangular distribution of reservoir temperature should be chosen; pert distribution models were chosen by the specific heat and rock density based on the skewness of distribution characteristics and continuity of parameter variation. Taking the Jinhu sag as an example, the following models are chosen for the input parameters (Figure 4): the annual average temperature (T0) and water density (ρW) are given a constant uniform model; the geothermal reservoir area (A) is set to a range of uniform distribution model; the variation of geothermal reservoir temperature (T), porosity (ϕ) and thickness (D) are given a triangular distribution model; for the specific heat of the geothermal reservoir rock (cR) and water (cW), and the rock density (ρR), the pert distribution model is chosen. Taking into account the accuracy of the simulation results, the number of iterations in this study was set to be 10,000.

**Figure 4.** Distribution models of input parameters in the Jinhu sag ((**a**) water density; (**b**) the geothermal reservoir area; (**c**) the annual average temperature; (**d**) the geothermal reservoir rock density; (**e**) the specific heat of rocks; (**f**) the specific heat of water; (**g**) geothermal reservoir temperature; (**h**) the geothermal reservoir rock porosity; (**i**) the thickness of the geothermal reservoir)).

#### **4. Database**

#### *4.1. Temperature Logs*

The determination of parameter ranges and distribution models is key for the assessment of geothermal resources. Temperature is one of the most important parameters for geothermal resources assessment and steady-state temperature measurement in boreholes is the most direct and effective way to obtain the true temperature in deeper formations. From 2018 to 2019, temperature profiles of 99 boreholes were acquired [4]. The temperature and depth data of the boreholes was measured using a consecutive logging system of a 5000 m long cable and a Platinum thermal resistance sensor. This field work is the first systematic steady-state temperature measurements in the whole NJB, with the same instruments and researchers. To date, we have obtained high-quality temperature logs from a

total of 110 boreholes (Figure 1) and selected 28 representative temperature logs of each sub-structural unit in the NJB for display, as shown in Figure 5.

**Figure 5.** The representative temperature-depth profiles for each sub-structural unit in the North Jiangsu Basi (**a**) Yanfu Depression; (**b**) Jianhu Uplift; (**c**) Dongtai Depression; (**d**) Tan-Lu Fault Zone. (modified from Wang, et al. [4].

The water level of each borehole (generally less than 40 m) can be determined according to the temperature logs (Figure 5). The temperature increases (nearly) linearly with depth more than 50 m, showing a conductive type [30,31]. Based on the depth and temperature data, the temperature gradients (the portion above (including) the Yancheng Formation) of each sub-structural unit in the NJB were calculated. Regionally, the middleeastern Jianhu Uplift and the eastern Dongtai Depression is a zone of the high temperature gradient, with the average temperature gradient value higher than 35 ◦C·km<sup>−</sup>1. Vertically, the temperature gradient in the NJB increases with depth, and the temperature gradients of Ny and E2s-E2d are generally lower than that of E1f-K2t. Formations of the Neogene and above are characterized by a low temperature gradient, generally between 15–30 ◦C·km<sup>−</sup>1, which may be the result of heat redistribution and groundwater activity in the shallow. The high sand content in the Ny, E2s, and E2d formations, resulting in high thermal conductivity and low temperature gradients; whereas, the Palaeocene and Cretaceous formations have a relatively high mud content and thus a high temperature gradient of 30–40 ◦C·km<sup>−</sup>1.

#### *4.2. Thermal Conductivity*

Thermalphysical properties, particularly thermal conductivity measurement, are fundamental to the study of regional deep geothermal field. We carried out thermal conductivity test on 264 samples (Table 1; sample locations, see Figure 1). moreover, we collected the thermal conductivity data tested by Wang and Shi, and Wang et al. [32,33] (Table 1). According to the data above, we established the thermal conductivity column of different formations (Erathem) in the NJB (Table 1).


**Table 1.** Thermal conductivity of the formations in the North Jiangsu Basin.

\* Previous work is from Wang and Shi, and Wang et al. [32,33].

The results show that the average thermal conductivity of the strata in the region fluctuates considerably, with the smallest value being the clay and sand of the Quaternary with a thermal conductivity of 0.6 W·m−1·K−<sup>1</sup> and the largest value being the Silurian S2m1 with a thermal conductivity of 8.0 W·m−1·K<sup>−</sup>1. In general, the thermal conductivity shows a gradual decrease with the strata from old to new. The thermal conductivity of the Palaeozoic strata is generally high, greater than 3.0 W·m−1·K−1; the mesozoic strata have the second highest thermal conductivity, mostly 2.0–3.0 W·m−1·K−1; and the Cenozoic strata have the lowest thermal conductivity, between 1.5 and 2.5 W·m−1·K−1. The stratigraphic thermal conductivity of the NJB can be roughly divided into two parts, namely the shallow low thermal conductivity section and the deep high thermal conductivity section. The shallow Cretaceous-Quaternary thermal conductivity is low (<2.7 W·m−1·K<sup>−</sup>1), especially the Late Paleozoic-Quaternary thermal conductivity is mostly less than 2.0 W·m−1·K<sup>−</sup>1, which can provide a good cover for heat preservation. The thermal conductivity of the deeper Ediacaran-Jurassic is relatively high, generally greater than 3.0 W·m−1·K<sup>−</sup>1, especially the thermal conductivity of the Silurian-Devonian and Ediacaran-Cambrian can be more than 4.0 W·m−1·K−1, and the thermal conductivity of the Dengying Formation can be more than 6.0 W·m−1·K<sup>−</sup>1, all of which can be used as good thermal storage layers.

#### *4.3. Input Parameters*

The temperature range of the geothermal reservoirs is an important parameter for geothermal resource assessment. In combination with the geophysical profiles, we used Equation (4) to calculate the temperature range corresponding to each thermal reservoir at depth (z), based on the average heat flow (Qs) in the various sub-structural units in the NJB. The temperature at a given depth for heat flow in each sub-structural unit in the NJB with a constant surface heat flow Qs, surface temperature T0, thermal conductivity λ, and radiogenic heat production A is defined as

$$\mathbf{T(z)} = \mathbf{T}\_0 + \mathbf{z} \, \mathbf{Q\_s}/\lambda - \mathbf{z}^2 \, \mathbf{A}/\text{(2\lambda)}\tag{4}$$

where <sup>λ</sup> (W·m−1·K<sup>−</sup>1) and A (μW·m<sup>−</sup>3) represent the measured thermal conductivity and heat production of the rock, respectively.

The main parameters, such as porosity and heat production, are measured data, and detailed explanations are given in Wang, et al. [4]. In addition, we carried out measurements for specific heat, porosity, and density for the main lithologies of each geothermal reservoirs, other parameters were mainly set regarding the literature: the annual average temperature of each sub-structural units [34] and water density (1000 kg·m<sup>−</sup>3) were taken as constant values; the specific heat of the water in the geothermal reservoir varied with temperature [35]; the variation in thickness and area of the geothermal reservoir was mainly based on actual geophysical profiles and sedimentological evidence [4,13,14,28,36]. The parameters for the Yancheng Formation, the Sanduo Formation, the Dai'nai Formation, and the carbonate thermal reservoirs are reported in Tables 2–5, respectively (the most likely values are shown in parentheses).

**Table 2.** Parameters for the Yancheng Formation thermal reservoir.



**Table 2.** *Cont.*



**Table 4.** Parameters for the Dai'nan Formation thermal reservoir.



**Table 5.** Parameters for the carbonate thermal reservoir.

#### **5. Results**

*5.1. Simulation of the Cenozoic Sandstone Thermal Reservoirs*

According to the parameters in Table 2, we used the monte Carlo method to calculate the geothermal resource base for the thermal reservoir of the Neogene Yancheng Formation in the NJB. Figure 6a indicates the geothermal resource base in which the value varies from 6.6 × <sup>10</sup><sup>20</sup> J to 1.2 × 1021 J (mean value (8.6 ± 0.6) × 1020 J), with the highest probability being 0.86 × <sup>10</sup><sup>21</sup> J (probability > 8 %) and a 90 % probability of a range of (0.76–0.96) × <sup>10</sup><sup>21</sup> J. It shows a geothermal resource potential of 3.5 × 1016 <sup>J</sup>·km−<sup>2</sup> by dividing the average value by the Yancheng thermal reservoir area.

**Figure 6.** Simulation results for the primary thermal reservoirs in the North Jiangsu Basin ((**a**) Yancheng thermal reservoir; (**b**) Sanduo thermal reservoir; (**c**) Dai'nan thermal reservoir; (**d**) carbonate thermal reservoir).

Figure 6b shows the potential of the geothermal resource in which the value changes from 5.0 × 1020 J to 1.6 × 1021 J with a mean value of (9.1 ± 0.2) × 1020 J), with nearly 7% probability of the most likely value of 0.86 × 1021 J, and a probability of < 5% that the geothermal resources > 1.2 × 1021 J or < 0.67 × <sup>10</sup><sup>21</sup> J. Dividing by area, it indicates a geothermal potential of 4.5 × <sup>10</sup><sup>16</sup> <sup>J</sup>·km<sup>−</sup>2.

Figure 6c presents a geothermal resource potential for the Dai'nan thermal reservoir of 3.2 × 1020 J to 1.1 × 1021 J with a mean value of (6.6 ± 0.1) × <sup>10</sup><sup>20</sup> J), with a 90 % probability that the energy potential is (0.45–0.90) × <sup>10</sup><sup>21</sup> J. Considering the geothermal resource potential per unit of thermal storage area, it indicates a potential of 6.7 × 1016 <sup>J</sup>·km<sup>−</sup>2.

#### *5.2. Simulation of the Carbonate Thermal Reservoirs*

Figure 6d indicates that the geothermal resource potential for the carbonate thermal reservoir sums to 9.9 × <sup>10</sup><sup>21</sup> J to 1.5 × 1022 J with a mean value of (1.2 ± 0.1) × 1022 J). From Figure 6d, we can get that the geothermal potential is above 1.3 × 1022 J or below 1.1 × <sup>10</sup><sup>22</sup> J in which the probability is less than 5 %, and the most likely value of the whole is 1.2 × <sup>10</sup><sup>22</sup> J (nearly 8%). The geothermal resource potential per unit of the thermal storage area is 3.7 × 1017 <sup>J</sup>·km<sup>−</sup>2, which is much larger than that of the Cenozoic sandstone thermal reservoirs. The geothermal resource base of the carbonate thermal reservoir in the NJB is extremely large, indicating a great potential for geothermal resources.

#### *5.3. Monte Carlo Simulation Results for the Total Geothermal Resources in the North Jiangsu Basin*

The total geothermal resource base of the Cenozoic sandstone and carbonate thermal reservoirs were estimated through running the monte Carlo simulation for iterations, with the results being shown in Figure 7. It indicates the geothermal resource base in which the value ranges from 1.2 × 1022 J to 1.7 × 1022 J (mean value (1.4 ± 0.1) × 1022 J), with the highest probability being 1.4–1.5 × 1022 J (probability nearly 7%), and a 90 % probability of a range of (1.3–1.6) × 1022 J. The Yancheng, Sanduo, Dai'nan, and carbonate thermal reservoirs account for about 6%, 6%, 5%, and 83% of the total geothermal resource base, respectively.

**Figure 7.** Geothermal resource potential simulation results for the sum of the four thermal reservoirs.

#### **6. Conclusions**

This study used the monte Carlo method to perform the detailed assessment of geothermal resources in the North Jiangsu Basin. Geothermal resources within the Cenozoic sandstone and carbonate (<10 km) thermal reservoirs total (1.2–1.7) × <sup>10</sup><sup>22</sup> J, with a mean value of (1.4 ± 0.1) × 1022 J and a most likely value of 1.4 × 1022 J, which is equivalent to 4.9 × 1013 tons of standard coal heat content. moreover, the most promising carbonate thermal reservoir owns an average geothermal resource potential of (1.2 ± 0.1) × 1022 J, about five times that of the Cenozoic thermal reservoirs. The abundant geothermal energy of the North Jiangsu Basin may alleviate energy shortages to a certain extent and promote exemplary geothermal heating in the 'hot in summer and cold in winter' region of the Yangtze River Delta.

**Author Contributions:** Conceptualization, Y.W. (Yibo Wang) and Y.B.; methodology, Z.W.; software, J.H.; validation, J.H. and S.H.; formal analysis, Y.W. (Yibo Wang) and Y.W. (Yaqi Wang); investigation, Y.W. (Yibo Wang), Z.W., D.H. and Y.B.; data curation, L.W. and Z.W.; writing—original draft preparation, Y.W. (Yibo Wang); writing—review and editing, Y.B.; visualization, Z.W.; supervision, S.H.; project administration, S.H.; funding acquisition, L.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was accomplished under the support of the Jiangsu Fund Project (Special Study on the Investigation of Geothermal Field and Crustal Thermal Structure in the North Jiangsu Basin, JITC-1902AW2919; JITC-1802AW2292).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to [Privacy].

**Acknowledgments:** We thank the Geological Survey of Jiangsu Province very much, who provided enormous amount of help on the borehole selection and sample collection that are significant for our temperature logging work and geothermal resources assessment.

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

#### **References**


## *Article* **Techno-Economic Assessment of Mobilized Thermal Energy Storage System Using Geothermal Source in Polish Conditions**

#### **Dominika Matuszewska 1, Marta Kuta <sup>1</sup> and Piotr Olczak 2,\***


Received: 21 May 2020; Accepted: 28 June 2020; Published: 2 July 2020

**Abstract:** The paper considers technical and economic possibilities to provide geothermal heat to individual recipients using a mobile thermal storage system (M-TES) in Polish conditions. The heat availability, temperature and heat cost influence the choice of location—Ba ´nska Ni ˙zna, near Zakopane in the southern part of the Poland. The indirect contact energy storage container was selected with phase change material characterized by a melting temperature of 70 ◦C and a heat storage capacity of 250 kJ/kg, in the amount of 800 kg. The economic profitability of the M-TES system (with a price per warehouse of 6000 EUR, i.e., a total of 12,000 EUR—two containers are needed) can be achieved for a heat demand of 5000 kWh/year with the price of a replaced heat source at the level of 0.21 EUR/kWh and a distance between the charging station and building (heat recipient) of 0.5 km. For the heat demand of 15,000 kWh/year, the price for the replaced heat reached EUR 0.11/kWh, and the same distance. In turn, for a demand of 25,000 kWh/year, the price of the replaced heat source reached 0.085 EUR/kWh. The distance significantly affected the economic profitability of the M-TES system—for the analyzed case, a distance around 3–4 km from the heat source should be considered.

**Keywords:** geothermal energy; mobile thermal energy storage (M-TES); phase change material (PCM); LCOH; heat transport; renewable energy source

#### **1. Introduction**

For several years, a global energy transformation due to the growing awareness of greenhouse gas (GHG) emissions and fossil fuel depletion has been observed. In this context, the use of renewable energy sources, including geothermal energy, became one of the main strategies of this transformation. Geothermal energy is produced mainly as a result of the decay of potassium, uranium and thorium radioisotopes [1] in the planetary interior as thermal energy [2]. Compared with other renewable energy sources, it has many advantages, among others, large reserves of geothermal sources, wide distribution and scope of applications, good stability, and high utilization efficiency of use. Additionally, once stabilized, they can be used for a long period of time [3,4]. In addition, geothermal sources are the only renewable energy source that is not affected by solar radiation and gravitational attraction between the sun and the moon [5]. However, the main problem with geothermal resources is that the speed of geothermal reservoir exploitations can be faster than heat replacement, and nowadays this is considered as the challenge connected with geothermal heat utilization. It largely depends on such factors as geological time scale, geothermal applications and used-heat reinjection methods [6].

The great potential of geothermal resources used for electric power generation has lately been recognized more and more [7]. However, geothermal heat has been widely used directly for centuries by humans and animals, and can be found in various regions of the world [8]. Heat appears on the earth's surface in various forms that affect the way it is used. The wide spectrum of geofluid temperature determines the way it is used, and the range of applications as a temperature function is shown in Figure 1, which is based on the Lindal diagram [9,10]. Currently, geothermal heat can be utilized in two forms: direct use of hot springs for bathing, aquaculture, geothermal heating and indirect in electric power generation [11,12]. Global geothermal power capacity by the end of 2019 totaled 13.93 GW, with annual electricity generation reaching 85.98 TWh in 2017 (most recent data) [13]. Geothermal power plants can be categorized by the technologies exploiting geothermal resources, according to power plant configuration: dry steam plant (23%), single flash (41%), double flash (19%), binary (14%) and advanced geothermal energy conversion systems such as tripe flash, hybrid flash-binary systems, hybrid fossil-geothermal systems, hybrid back pressure system (3%) [14,15]. A two-group classification (binary cycles for lower well enthalpies and steam cycles for higher enthalpies) has been suggested by Valdimarsson [16].

**Figure 1.** Modified Lindal diagram. Source: own study, based on [9,10].

The geothermal heat can be used directly in many applications depending on a temperature range (mostly used reservoirs with temperature between 20 ◦C and 150 ◦C) [17,18]. In 2018, the direct utilization of geothermal energy was identified and reported for 82 countries with an estimated global installed thermal capacity at 70.33 GWth (with 0.265 capacity factor and a growth rate 7.7% annually) [19,20]. Lund and Boyd [21] show that heat energy utilization reaches 163,287 GWh and identify nine utilization categories: geothermal heat pump (90,297 GWh), bathing and swimming (33,147 GWh), space heating (24,493 GWh), greenhouse heating (7347 GWh), aquaculture pond heating (3266 GWh), industrial uses (2939 GWh), cooling and snow melting (653 GWh), agricultural drying (653 GWh) and others (490 GWh). As it can be seen from these categories, the most popular are geothermal heat pumps, which can be used in many applications under different conditions [22,23]. Figure 1 shows that the direct use (DU) of geothermal energy varies in range from high-temperature industrial applications to very low-temperature aquaculture. Additionally, in some ranges, the

geothermal heat utilized in DU can not only meet thermal requirements for some applications, but also can be used in combined heat and power configurations (CHP) [24]. Matuszewska and Olczak [25] show the modified binary system to provide heat and electricity simultaneously from low-enthalpy geothermal reservoirs (normally used for providing district heating in Polish conditions). Carotenuto et al. [26] analyze novel solar–geothermal low-temperature district heating and cooling and domestic hot water systems. The case study has been developed for the district area of Monterusciello in Southern Italy, where the geothermal source obtains 55 ◦C. Hepbasli and Canakci [27] show, in the example of geothermal district heating in Izmir-Balcova, that geothermal energy is cheaper than other energy resources in Turkey (including fossil fuels) and thus it can make a contribution in GHG emissions. Sander's [28] study challenges opportunities of geothermal district heating systems based on the analysis of China, Germany, Iceland, and the U.S. country assessments.

One of the ways to deliver geothermal heat to the recipient in a less conventional way than with the help of standard heat pipelines is to use a mobile heat storage, which will be the subject of the following analysis.

Mobilized thermal energy storage systems are based on the use of heat in a place other than its generation, collection, and storage. The basic assumption is storing waste heat, excess heat or heat produced from renewable energy sources in an amount larger than local demand. Heat transported using mobile heat containers can be stored using phase change materials based on the use of latent heat or using zeolites that use thermochemical transformations. The article cites examples of both possibilities; however, the analysis focuses on mobile thermal energy storage (M-TES) based on phase change materials.

The stored heat can be transported by road, rail, water, or a combination of these. Due to the greatest universality, road transport is the most popular. It can take place using trucks (with large sizes of the container) or cars equipped with a trailer (with smaller container sizes). The heat is transported to recipients, where the M-TES is discharged.

Potential heat sources can be selected from: renewable heat sources, including geothermal heat, industrial plants and biogas plants producing large amounts of waste heat. A number of requirements are set for heat sources: the possibility of effective heat collection and transfer to the phase change material (PCM)-container, an appropriate heat temperature enabling material full phase transformation, constant heat availability, an adequate amount of available heat, high thermal efficiency. Additionally, the heat transportation process from heat source to the recipient location must meet the relevant requirements; the distance cannot be too large. Distance depends on the parameters of the PCM-container and the heat demand at the recipient. The literature shows that it is not advisable for the distance to be more than 30 km [29]. The M-TES needs to be thermally protected to minimize heat losses in transport and during unloading. The M-TES system assumes the use of two heat containers working in parallel and enabling constant access to heat. When the first one is unloaded at the recipient location, the second one is loaded and waits for the delivery to the user.

The place of heat reception should be characterized by heat demand in the right range, suitable heat parameters and operating mode. Potential recipients can be buildings like single and multi-family houses, public buildings, and industry (transport within the same industrial plant).

The mobile thermal energy storage system's work cycle is presented in Figure 2.

Due to the fact that the M-TES parameters can be flexibly adapted to the parameters of the heat source and the recipient, it is possible to select the appropriate PCM filling the heat container, where parameters will be adapted to the temperature of the heat source, and at the same time, will also be adapted to the needs of the recipients. At the stage of M-TES system design, parameters such as the frequency of deliveries or the size of the container dedicated to a given solution should be considered. These parameters are connected to the size and seasonality of the heat demand. The heat demand (at every hour) is effected by a number of parameters, such as: building size and its condition (statistical [30] and methodical approach [31]), number of residents [31], building usage hours and characteristics [32], and requirements for maintaining the temperature at a certain level, related to weather conditions [33].

**Figure 2.** Work cycle of the mobile thermal energy storage system.

In recent years, projects dedicated to mobile heat transport technology have been conducted, as well as the first projects completed with commercial implementations. Selected examples of such solutions are presented below.

Yan Wang et al. [34] designed and tested an M-TES system dedicated to renewable energy and recovered industrial waste heat. The storage unit contained a heat exchanger in the form of compact storage tubes filled with PCM. The system uses 215 kg sodium acetate trihydrate (CH3COONa ·3H2O), which is characterized by: melting temperature: 331 K, density 1330 kg/m3, specific heat 3200 J/(kg·K), thermal conductivity 0.8 W/(m· ◦C) and enthalpy 260 kJ/kg. The tests showed that 1200 s would be enough to melt 215 kg of PCM filling 80 tubes. The tested system has been able to store 125,576 kJ of heat.

Weilong Wang et al. [35] built and tested lab-scale direct and indirect thermal energy storage containers filled with the phase change material erythritol, which is characterized by melting temperature 118 ◦C and heat density 330 kJ/kg. The system has been dedicated for the supply of heat recovered from the industry. Two options of M-TES were compared: direct- and indirect-contact containers, and it was concluded that the direct-contact M-TES can be charged and discharged faster than the indirect-contact container (for tested systems).

Marco Deckert et al. [36] optimized and tested an M-TES system based on the system commercially available. The system consists of two 20-foot-long tanks filled with sodium acetate trihydrate, which melts at 58 ◦C. The heat exchanger is built of 24 tubes. To enhance the charging and discharging process, tubes in one part of the storage were extended with graphite structures. The tested M-TES has been able to store up to 2 MWh of the heat, including 1.3 MWh of latent heat. The authors concluded that the technical and economic feasibility depends on the M-TES storage capacity but also on heat demand of the recipients. The tests indicated areas of improvement and demonstrated the possibility of commercial use for M-TES.

Xuelai Zhang et al. [37] proposed an M-TES system in the form of container with dimensions: 4.1 m × 2 m × 1.25 m (inner tank: 3.6 m × 1.8 m × 1.1 m). The container was filled with stainless steel balls (radius of 60 mm, 80 mm and 100 mm) with new phase change material consisting of 0.4% nanocopper + 99.6% erythritol, which is characterized by phase change temperature 118 ◦C and latent heat 362.2 kJ/kg. The authors concluded that the biggest impact on the results had the temperature of the oil bath (higher speed obtained in lower initial temperature) and the size of the balls (easier melting for smaller diameter, but heat dissipation is faster for larger sizes).

Andreas Krönauer et al. [38] present the results of one-year long tests where the M-TES was charged to 130 ◦C, using extraction steam from a waste incineration plant. The heat was transported to the recipient located 7 km from the heat source and used in the drying processes. The tank was filled with 14 tons of Zeolite and had a storage capacity of 2.3 MWh. The authors conducted analysis and identified areas for improvement to increase CO2 reduction and improve efficiency.

This paper focuses on the economical assessment of M-TES using geothermal heat in Polish conditions. There is a lack of information about the profitability of using M-TES for geothermal reservoir utilizations. Additionally, most research focuses on determination of the operating conditions of the geothermal reservoir, and on the amount of energy that can be extracted, but they do not consider the problems with large dispersion of recipients (which may cause an increase in the costs of building a heating network), as well as topography unfavorable for this type of solution. Considering this, the application of M-TES for geothermal heat is justified primarily in places where it is impossible or unprofitable to build a heating network.

In Poland, the temperature of geothermal water does not exceed 100 ◦C, and most reservoirs have a depth of 3000 m [39]. The most promising and interesting areas from a geothermal point of view are in the southern part of the country, the Podhale Basin. It has sufficient operation conditions for district heating, recreation, balneotherapy and other applications. Additionally, the heating system located there and operated by PEC Geotermia Podhala ´nska S.A. is currently the largest geothermal installation in Poland with a geothermal heating system capacity ca. 40.8 MWth (the total installed capacity of the heating plant is 80.5 MWth) [40]. The analyzed geothermal heating system is among the largest in Europe and is steadily upgraded [41].

The example of PEC Geotermia Podhala ´nska S.A. was considered in further analysis due to the location in a mountainous area, which may prevent the extension of a heating network, or significantly increase construction costs. In addition, the price of heat generation using geothermal sources is competitive in relation to other sources and amounts to approx. 0.0322 EUR/kWh, at the same time, the cost of building a heating network is high. Due to the large dispersion of recipients, as well as the landform, the price of heat from geothermal energy, after inclusion of the heat network cost, is high. PEC Geotermia Podhala ´nska S.A. currently has a heating network with a length of over 100 km, in order to connect new customers not located on the existing network, it is necessary to consider the expansion of the network in terms of the economics of the solution. Due to the high dispersion and differences in the foundation height of buildings, the heating network development costs are high. Therefore, based on global examples, it seems worthwhile to consider the use of M-TES. Figure 3 presents that the price for heat from geothermal energy is lowest compared with other geothermal district heating systems in Poland. This means that, if for PEC Geotermia Podhala ´nska S.A. the solution is economically justified, it can be further considered for facilities with a higher price for geothermal heat.

**Figure 3.** Map of Poland with geothermal sources and heat price (without distribution and charge for power). Source: own study based on [42,43].

Based on the references above, the M-TES for the use of geothermal heat form Ba ´nska Ni˙zna well was analyzed. The economic assessment of these problem was conducted. PEC Geotermia Podhala ´nska S.A. was selected due to the lowest geothermal heat prices (with no additional fee) and highest geothermal temperature from the production well. There are enough studies about M-TES for geothermal heat utilization, especially in Polish conditions. The novelty of this paper is in analyzing the economical profitability of applying such solutions in mentioned conditions. The paper is structured as follows: Section 2 presents used methods and assumptions; in Section 3, the results are presented; and finally, conclusions are provided in Section 4.

#### **2. Methods and Assumptions**

The main purpose of this analysis is to determine the economic profitability of the solution focused on replacing the current heat source in the building with geothermal heat transported by a M-TES. The analysis began by determining the relationship between the M-TES heat capacity and total building heat demand (THD). The next step was to determine the impact of various parameters on Economical Profitability. Finally, the impact of using M-TES transport systems on the reduction of carbon dioxide emission was estimated.

The calculations presented in this article were carried out for M-TES, which uses 800 kg of phase change material for heat storage (maximum amount of PCM that can be safely transported by car with a trailer with a permissible total weight of 3.5 Mg). The material that could potentially be used in this application was selected, taking into account the parameters of the heat source, requirements resulting from the characteristics of the place of heat reception and the phase change materials requirements typical for this type of applications (e.g., appropriate phase transition temperature, high value of heat storage capacity, the highest possible value of the thermal conductivity coefficient, stability of properties after many working cycles, nontoxic, non-harmful for human or environment, nonflammable).

Materials suitable for use in this application include: PureTemp 68 (melting temperature: 68 ◦C, heat storage capacity: 213 kJ/kg, thermal conductivity liquid/solid: 0.15/0.25 W/(m· ◦C) [44]; ATP 70 (melting temperature: 70 ◦C, heat storage capacity: 250 kJ/kg, thermal conductivity 0.2 W/(m· ◦C) [45]; GAIA HS PCM78 (melting temperature: 78 ◦C, heat storage capacity: 260 kJ/kg, thermal conductivity liquid/solid: 0.5/0.98 W/(m· ◦C) [46]; RT 70 HC (melting temperature: 70 ◦C, heat storage capacity: 260 kJ/kg, thermal conductivity 0.2 W/(m· ◦C) [47]. For calculation purposes it has been selected PCM characterized by following properties: temperature: 70 ◦C, heat storage capacity: 250 kJ/kg. The amount of PCM as well as the weight of the entire M-TES was selected in a way which it could be transported by car with a trailer with a permissible total weight of 3.5 Mg. The use of this type of solution was adapted to various transport requirements related to the landform around the heat source (Ba ´nska Nizna).

#### *2.1. Total Heat Demand*

The total heat demand (THD) in the range of 5000–25,000 kWh was adopted for the analysis. The limiting values for the above-mentioned range has been determined on the basis of:

1. Heat demand for households using coal-based heat sources. In the first place, these sources should be replaced by renewable energy, due to carbon dioxide and pollution emission to the air [48]. For coal-based heat sources, the average heat demand per m<sup>2</sup> of usable floor space is 222 kWh/(m2·year) for insulated buildings and 253 kWh/(m2·year) for non-insulated buildings [30]. The average usable floor space of households (in rural areas) is 108.3 m<sup>2</sup> and the total energy demand is 216 kWh/(m2·year) [48]. At the same time, the area up to several km from Ba ´nska Ni ˙zna is characterized by a large diversity of buildings (from small traditional one- and two-room houses to multi-story residential houses) as well as their energy efficiency. DHW demand was determined on the basis of average consumption of DHW per person 35 dm3/(person·day) [31] and the number of households in Lesser Poland Voivodeship [49], as well as on the basis of efficiency of heat production and transport for coal sources. The results are shown in Figure 4.


**Figure 4.** Demand for thermal energy in households in Poland. Data limited to coal sources. Source: own study, based on [30,31,50].

This value applies to heat for central heating (CH) and domestic hot water (DHW). Lower demand values (from 5000 kWh/year) were not taken into account due to the high share of heat losses (calculated according to formula no. 1) [51] from M-TES which is HL(M-TES, year) = 2300 kWh/year to the total amount of heat delivered to the building from M-TES (5000 kWh/year). It was assumed that M-TES thermal insulation will provide a two-time lower heat loss stream than in heat storage solutions typically used inside buildings.

$$\text{HL} \left( \text{M} - \text{TES}, \tau \right) = \frac{1 \text{h}}{1000} \cdot \left( \text{T} \left( \text{M} - \text{TES}, \tau \right) - \text{Ta} \left( \tau \right) \right) \cdot \text{ni} \cdot \text{chl} \cdot \sqrt{\text{V} \left( \text{M} - \text{TES} \right)}, \text{ kWh} \tag{1}$$

where:

HL(M-TES,τ)—heat losses from M-TES at the hour τ, kWh

T(M-TES,τ)—M-TES temperature at the hour τ, ◦C

Ta(τ)—outside temperature, ◦C

ni—number of M-TES for each house: 2 (one close to house, second for loading)

chl—coefficient for heat losses. It has been assumed from (as <sup>1</sup> <sup>2</sup> of value from [51]: 0.08), <sup>W</sup>/(K·dm3/2)

V(M-TES)—volume of M-TES, dm<sup>3</sup>

The calculations assume T(M-TES,τ) at the level of 70 ◦C. For the determination of outside temperature (Ta), data from the Typical Reference Year (Figure 3) [52] were used. On this basis, the annual sum of heat losses was calculated—Equation (2).

$$\text{HL}(\text{M}-\text{TES}, \text{year}) = \sum\_{\tau=1}^{8760} \text{HL}(\text{M}-\text{TES}, \tau) / \frac{\text{kWh}}{\text{year}} \tag{2}$$

The 25,000 is a defined upper safe limit for the use of two heat stores due to the rate of heat distribution at the recipient and its charging at the base station, e.g., in Geotermia Podhala ´nska.

The distribution of total heat demand THD for heating Q(CH) and for domestic hot water preparation Q(DHW) was assumed in a 3:2 ratio [30]. This ratio matters regarding the frequency of M-TES exchange (winter/summer relationship) and for heat losses during the transportation: a greater share of DHW means that, in warmer months, there are more trips (and exchanges) than in the case of a smaller share. However, the analysis does not include heat losses during transport, due to the fact that its time is a maximum of several minutes. In summary, it was considered that the ratio Q(CH)/Q(DHW) according to specific THDs was separated, and has no significant impact on further calculations.

The heat demand for heating, taken from the range of 3000–15,000 kWh, was divided into individual hours of the year according to the hourly temperature data from the Typical Reference Year for Zakopane (10 km from Ba ´nska Ni ˙zna) (Figure 5).

**Figure 5.** Outside temperature in Zakopane according to Typical Reference Year. Source: own study, based on [52].

Data on the heat demand for the preparation of domestic hot water were divided according to hours based on data from work [53].

The change in the amount of heat accumulated in M-TES in each hour was determined from the formula:

$$\mathbf{Q(M-TES,\tau)} = \mathbf{Q(M-TES,\tau-1)} - \text{HL}(\mathbf{M-TES,\tau}) - \mathbf{Q(CH,\tau)} - \mathbf{Q(DHW,\tau)}, \text{ kWh} \tag{3}$$

where:

Q(M-TES,τ)—heat accumulated inside the M-TES at the hour: τ, kWh Q(M-TES,τ − 1)—heat accumulated inside M-TES at the hour: τ − 1, kWh HL(M-TES,τ)—heat losses from M-TES at the hour: τ, kWh Q(CH,τ)—heat needed for the central heating purposes at the hour: τ, kWh Q(DHW,τ)—heat needed for the domestic hot water purposes at the hour: τ, kWh τ—hour (time)

When Q(M-TES, τ) would reach zero in a given hour of computing simulation, it is considered in the calculations that the M-TES was replaced with a fully charged one. In this case, Q(M-TES, τ—1) is replaced by a value equal to the heat capacity M-TES, which is 55 kWh.

#### *2.2. Economical Profitability and NPV*

The following assumptions were made in the analysis:


Fixed prices were assumed in the analysis due to the assumed proportional increase in the prices of heat (geothermal and other), transport costs (including labor costs).

Economical profitability (EP) M-TES solutions were calculated in comparison to other technologies for heat production in buildings on the basis of estimating the costs of their use over 20 years. EP values were presented as the difference between the cost of the current solution (e.g., coal boiler) and planned installation, which is M-TES. This difference was calculated for a 20-year perspective and referred to one year of use—Equation (4). Due to the fact that in both cases (base and M-TES), investment costs refer to the total heat demand, and often the service costs are not linearly dependent on the amount of heat demand, the presentation of the comparison results in a unitary way was considered as too high a simplification.

EP(dist., THD, LCOH.OS, M <sup>−</sup> TES.P) = THD·(LCOH.OS(THD) <sup>−</sup> LCOH.G(dist., M <sup>−</sup> TES.P, THD), EUR year (4)

where:

EP—economical profitability, EUR/year

dist.—distance between building and Geothermal Base Station, km

THD—total heat demand for building, kWh/year

LCOH.OS—Levelized Cost of Heating (for 20 years) for other than Geothermal source, equation 10 (based on [54]), EUR/kWh

LCOH.G—Levelized Cost of Heating (for 20 years) for Geothermal source, Equation (5), EUR/kWh MTES.P—M-TES price, EUR

$$\text{LCOH.G} \left( \text{dist.}, \text{M} - \text{TES.P}, \text{THD} \right) = \frac{\text{n} \cdot \text{(CGH(THD) + TC(dist., THD))} + \text{ni} \cdot \text{M} - \text{TES.P}}{\text{n} \cdot \text{THD}}, \frac{\text{EUR}}{\text{kWh}} \tag{5}$$

where:

LCOH.G—Levelized Cost of Heating (for n years) for Geothermal source, EUR/kWh n—lifetime of M-TES, year

CGH—Cost of Geothermal Heat (calculated by Equation (6)), EUR/year

TC—transport and work cost (calculated by Equation (7)), EUR/year

ni—number of M-TES for each house: 2 (one close to house, second for loading in base station) M-TES.P—M-TES price, EUR

THD—Total heat demand for the building, kWh/year

dist.—distance between building and Geothermal Base Station, km

$$\text{GCH}(\text{THD}) = (\text{THD} + \text{HL}(\text{M} - \text{TES}, \text{year})) \cdot \text{GHP}\_{\text{year}} \frac{\text{EUR}}{\text{year}} \tag{6}$$

where:

CGH—Cost of Geothermal Heat, EUR/year

THD—Total heat demand for the building, kWh/year

HL(M-TES, year)—heat losses from M-TES in whole year, kWh/year

GHP—unit geothermal heat price, EUR/kWh

$$\text{TC(dist., THD)} = \text{M} - \text{TES.} \text{EX}(\text{THD}) \cdot \left( \text{TWC} \cdot \frac{1}{60} \cdot (\text{TCMs} + \text{UTT} \cdot \text{dist.} \cdot 2) + \text{dist.} \cdot 2 \cdot \text{TTC} \right) \cdot \frac{\text{EUR}}{\text{year}} \tag{7}$$

where:

TC—transport and work cost, EUR/year M-TES.EX—number of M-TES exchanging per year TWC—total work cost, EUR/h TfCMs—time for M-TES exchanging, min UTT—unit transport time, min/km

TTC—total transport cost, EUR/km

dist.—distance between building and Geothermal Base Station, km

THD—Total heat demand for the building, kWh/year

To use geothermal heat with the use of M-TES, every household needs to have two heat containers; while one of them is used to supply heat to the building, the other one is charged at the charging station in Geotermia Podhala ´nska S.A. (in the Geothermal Base Station). The cost of a single container is determined from the following relationship:

$$\text{M}-\text{TES.P} = \text{Phase change material cost} + \text{container cost} + \text{trailer cost} \tag{8}$$

LCOH.OS in the calculations was adopted in the interval from 0.05 to 0.25 EUR/kWh. LCOH.OS was not clearly defined because it can differ for every building. Theoretically, the upper limit of LCOH.OS should be slightly higher than the price of electricity in Poland, which is currently around 0.15 EUR/kWh for domestic consumers [30]. However, the adopted range (up to 0.25 EUR/kWh) also considers the possible scenario of an increase in energy prices or the possibility of scaling the proposed solution outside the Poland. An example of the LCOH.OS calculation for a building (THD = 10,000 kWh/year) heated with a coal boiler was made using Equation (9).

$$\text{LCOH.OS} = \frac{\text{n} \cdot (\text{FC.OS}(\text{THD}) + \text{WC}(\text{THD})) + \text{nl-OS.P}}{\text{n-THD}} / \frac{\text{EUR}}{\text{kWh}} \tag{9}$$

where:

LCOH.OS—Levelized Cost of Heating (for n years) for other than Geothermal source, EUR/kWh n—lifetime of M-TES, year

FC.OS—total cost of fuel for other than Geothermal heat source, EUR/year

THD—Total heat demand for the building, kWh/year

WC—work cost, EUR/year

nl—correction coefficient of lifetime for other source = n/lifetime of other source

OS.P—other source price, EUR

The analysis of the example also includes: the cost of fuel 200 EUR/Mg of coal with a calorific value of 24 MJ/kg, efficiency 80% [31]. It displays the cost of 0.0375 EUR/kWh. FC.OS = 375 EUR/year. The cost of service was associated with the working time, and it was estimated as about <sup>1</sup> <sup>2</sup> hours per day in the heating season. In total, it is about 100 h/year for the town of Ba ´nska Ni˙zna. The cost of service, due to the work mainly being the owners (no additional costs), was valued as <sup>1</sup> <sup>2</sup> TWC. The total work cost was obtained: WC <sup>=</sup> 100 h <sup>×</sup> <sup>1</sup> <sup>2</sup> × 5 EUR/h = 250 EUR. For a 10,000 kWh/year demand, a unit cost of 0.025 EUR/kWh is implied. The cost of the heat source, i.e., in this case, a coal boiler, was estimated at EUR 1500 (OS.P). Assuming a lifetime of 10 years and THD = 10,000 kWh/year, the unit investment cost is 0.015 EUR/kWh. The cost of auxiliary energy (consumed in the building) was not included in the analysis since it is needed in virtually every solution, including the application of the M-TES system. The total for the above example with a coal boiler in the conditions of Ba ´nska Ni˙zna, Zakopane and the surrounding area of LCOH.OS is 0.0775 EUR/kWh.

In addition, Net Present Value (NPV) was calculated, according to Formula No. 10, for the following assumptions:


$$\text{NPV} = \left[ \sum\_{t=1}^{n} \frac{\text{CF}\_{t}}{(1+\mathbf{r})^{t}} \right] - \text{I0} \tag{10}$$

*Energies* **2020**, *13*, 3404

where:

CFt is the cash flow in the year t, EUR t—year of the analysis n—lifetime of M-TES, year I0 = ni·M-TES.P initial investment value, EUR r—discount rate, %

$$\text{CF}\_{\text{t}} = \left( \text{CGH}(\text{THD}) + \text{TC}(\text{dist., THD}) \right) - \text{THD-LCOH.OS}\_{\text{year}} \frac{\text{EUR}}{\text{year}} \tag{11}$$

where:

CFt is the cash flow in the year t, EUR

CGH—Cost of Geothermal Heat, EUR/year

THD—Total heat demand for the building, kWh/year

TC—transport and work cost, EUR/year

dist.—distance between building and Geothermal Base Station, km

LCOH.OS—Levelized Cost of Heating (for n = 20 years) for other than Geothermal source, EUR/kWh

#### *2.3. Carbon Dioxide Emission*

The following specific emission values were assumed for the calculations:

	- coal boiler: 0.375 kg CO2/kWh [53,58],
	- gas boiler: 0.200 kg CO2/kWh [58,59],

The reduction of CO2 emissions caused by replacing the current heat source with a geothermal source was calculated according to Formula No. 12:

$$\text{RE.CO2}(\text{THD}, \text{dist.}) = \text{THD} \cdot (\text{EM.CO2} \cdot \text{OS} - \text{EMCO2} \cdot \text{G}) - \text{EMCO2} \cdot \text{T} (\text{M} - \text{TES} \cdot \text{EX}(\text{THD}), \text{dist.}), \frac{\text{k\_f^\*CO}\_2}{\text{year}} \tag{12}$$

where:

RE.CO2—reduction of CO2 emissions, kgCO2/year

THD—Total heat demand for the building, kWh/year

dist.—distance between building and Geothermal Base Station, km

EM.CO2.OS—emission from other than Geothermal source, kg CO2/kWh

EM.CO2.G—Emission connected with Geothermal source, kg CO2/kWh

EM.CO2.T—annual CO2 emissions related to the M-TES transportation, determined from the formula no. 13, kg CO2/year

M-TES.EX—Mobile Thermal Energy Storage exchanged

$$\text{EX.CO2.T(M - TES.EX(THD), dist.) = 2 \cdot dist.M - TES.EX(THD) \cdot EM.CO2.KM} \tag{13}$$

where:

EM.CO2.T—annual CO2 emissions related to the M-TES transportation, kg CO2/year M-TES.EX—Mobile Thermal Energy Storage exchanged THD—Total heat demand for the building, kWh/year

dist.—distance between building and Geothermal Base Station, km

EM.CO2.OS—emission from other than Geothermal source, kg CO2/kWh

EM.CO2.KM—emission from car per km, kg CO2/km

Emissions associated with transport are linearly dependent on the distance. For example, for a distance of 5 km between Geothermal Base Station and Building and M-TES.EX = 200/year, transport related emissions would 480 kg CO2/year. A total of 2000 km would be driven.

#### **3. Results**

#### *3.1. Analysis of M-TES System Operation in Relation to the Total Heat Demand Size for the Building*

For different levels of total heat demand for the building (THD), the hourly heat demand (Q(CH) and Q(DHW)) as well as heat losses were calculated. The results for THD = 25,000 kWh/year in the winter period are shown in Figure 6, and in the summer period in Figure 7. The M-TES charge level (initial Q(M-TES)) 55 kWh in a given hour means that M-TES exchanged occurred in this hour.

**Figure 6.** Heat demand and heat losses and heat accumulated in mobile thermal storage system (M-TES) for each hour in the 1–240-h period of the year (January).

**Figure 7.** Heat demand and heat losses and heat accumulated in M-TES for each hour in the 4000–4240-h period of the year (June).

As can be seen in Figure 6, on the second day of the year, M-TES should be changed twice (for THD = 25,000 kWh/year). In turn, during the warmer period (June), i.e., between 4000–4240 h of the year, the tank should be replaced much less frequently (six times in 10 days; in the example for the winter, the sum was 10), and the heat losses are also lower. On aggregate for months, the number of M-TES exchanges is shown in Figure 8.

For THD = 5000 kWh/year, the maximum number of exchanges per month is 19, for THD = 15,000 it is 48, and for 25,000 it amounts to 78 (December). A feasibility analysis of the M-TES replacement in the context of various types of restrictions was carried out. The results presented below show the example, based on the calculations for the last value (78). The first of the analyzed restrictions is the M-TES charging and discharging rate: 78 exchanges mean almost three exchanges per day—the average charging power is approx. 7 kW. The maximum estimated charging power of 14 kW is realized with a flow temperature of 86 ◦C and a return temperature of 72 ◦C. Under these conditions, the flow rate would be 0.25 dm3/s. For THD = 25,000, the simulation received nine days with four M-TES exchanges per day (this means the average charging power of 9 kW). A similar situation is found for the M-TES discharging process. Further research includes work related to the irregularity of the M-TES maximum charging/discharging rate.

**Figure 8.** Number of M-TES exchanges (M-TES.EX) during each month of the year for different value of total heat demand (THD): 5000, 15,000 and 25,000 kWh/year.

In total, the number of exchanges in the year M-TES (M-TES.EX) is 147, 337 and 533, respectively. The number of exchanges needed to ensure the demand of THD = 15,000 kWh/year is 2.3 times higher than for THD = 5000 kWh/year, while the ratio of demand is three times higher. The reason for this phenomenon is the practically constant level of heat losses (HL (M-TES)) independent of the THD value in the adopted range.

The second possible restriction affecting the number of exchanges is the level of transport complexity in the mountainous area and—resulting from this—the transport time of one of the two M-TES dedicated for the building using this system. The longer the transportation time, the more time should be deducted for the M-TES charging process. The fully charged M-TES is transported to households and replaced with the empty M-TES (Q(M-TES) ≈ 0). Therefore, the discharging process is interrupted only during the M-TES replacements. The remaining time is spent on the M-TES charging process, excluding the transport to Geothermal Base Station and connection to the charging spot.

The use of the proposed M-TES system (with a heat capacity of 55 kWh) for a greater heat demand (more than THD = 25,000 kWh/year) may result in the M-TES dual storage system being insufficient due to the speed of heat use (recipient) or M-TES charging. Therefore, the maximum heat demand assumed for the building was checked. For this purpose, weather data from the last 20 years 2000–2019 were used (Figure 9). ERA5-Land data were analyzed for the location of Ba ´nska Ni ˙zna (Figure 3.) [61,62].

**Figure 9.** Outside temperature (Ta) in Banska Nizna: average temperature in each day hour in each month in 2000–2019.

The lowest average monthly temperature in the period 2000–2019 was recorded in February 2012. For this month, detailed weather data were obtained from ogimet.com. Then, the number of M-TES

exchanges during the following days was determined. The highest accepted value of THD = 25,000 was taken into account. The points of exchanges are shown in Figure 10 in combination with the outside temperature. For the first three days of February 2012, the estimated number of M-TES exchanges would be five per day. This number of exchanges means a minimum charging or discharging time of four hours, which means a necessary charging power of 14 kW.

**Figure 10.** Outside temperature (Ta) in Zakopane and number of M-TES exchanges (×) calculated for February 2012. Source: own study, based on [63].

#### *3.2. Economical Profiltability*

The results in the form of Economical Profitability for THD = 5000, 10,000, 15,000, 20,000 and 25,000 kWh/year are presented in the distance, along with the Levelized Cost of Heating (LCOH) function of the current heat source (LCOH.OS): Figures 11–15: (a) for M-TES price = EUR 2000, (b) for M-TES price = EUR 6000 (c) for M-TES price = EUR 10,000.

**Figure 11.** The Yearly Economical Profitable as a function of Levelized Cost of Heating from non-geothermal source (LCOH.OS) and distance between Geothermal Base Station and House, for Total Heat Demand (THD) = 5000 kWh/year and for: (**a**) M-TES price = 2000 EUR; (**b**) M-TES price = 6000 EUR; (**c**) M-TES price = 10,000 EUR.

**Figure 12.** The Yearly Economical Profitable as a function of Levelized Cost of Heating from non-geothermal source (LCOH.OS) and distance between Geothermal Base Station and House, for Total Heat Demand (THD) = 10,000 kWh/year and for: (**a**) M-TES price = 2000 EUR; (**b**) M-TES price = 6000 EUR; (**c**) M-TES price = 10,000 EUR.

**Figure 13.** The Yearly Economical Profitable as a function of Levelized Cost of Heating from non-geothermal source (LCOH.OS) and distance between Geothermal Base Station and House, for Total Heat Demand (THD) = 15,000 kWh/year and for: (**a**) M-TES price = 2000 EUR; (**b**) M-TES price = 6000 EUR; (**c**) M-TES price = 10,000 EUR.

**Figure 14.** The Yearly Economical Profitable as a function of Levelized Cost of Heating from non-geothermal source (LCOH.OS) and distance between Geothermal Base Station and House, for Total Heat Demand (THD) = 20,000 kWh/year and for: (**a**) M-TES price = 2000 EUR; (**b**) M-TES price = 6000 EUR; (**c**) M-TES price = 10,000 EUR.

**Figure 15.** The Yearly Economical Profitable as a function of Levelized Cost of Heating from non-geothermal source (LCOH.OS) and distance between Geothermal Base Station and House, for Total Heat Demand (THD) = 25,000 kWh/year and for: (**a**) M-TES price = 2000 EUR; (**b**) M-TES price = 6000 EUR; (**c**) M-TES price = 10,000 EUR.

A profitability of at least zero will be achieved at M-TES price = 6000 EUR for the total heat demand (THD):


• 25,000 kWh/year; at the price LCOH.OS = 0.085 EUR/kWh and distance = 0.5 km; at the price 0.18 EUR/kWh and distance = 6 km.

At higher demand values (THD) from the analyzed range, the results in the form of EP assume higher values, because the unit cost of M-TES decreases. This relationship is also influenced by the fact that M-TES price is usually higher than half of the cost of another heat source converted into 20 years of lifetime ( <sup>1</sup> <sup>2</sup> ·nl·OS.P). In addition, the above-mentioned dependence (EP on M-TES.P) is affected by the HL (M-TES)/THD ratio, which is lower at higher THD values, because the HL (M-TES) value is practically constant. HL (M-TES) depends, among other things, on temperature in M-TES and the outside temperature.

An example of the impact of THD on EP is illustrated by a comparison of Figures 11b and 15b (for M-TES.P = EUR 6000). For the same LCOH.OS price equal to 0.15 EUR/kWh, five times more heat demand is added to the increase of EP from −300 EUR/ear to 1500 EUR/year, at distance = 0.5 km. When LCOH.OS is 0.1 EUR/kWh, EP achieves a result of −550 EUR/ear and 200 EUR/year respectively. Additionally, in Figure 15b, there is a higher slope of the isoline compared to Figure 11b. This demonstrates the relatively smaller impact of distance on EP compared to the effect of LCOH.OS on EP. A similar comparison (THD = 5000 and 25,000 kWh/year) was made for LCOH.OS equal to 0.1 EUR/kWh and variable M-TES.P—Figure 16. The impact of M-TES price clearly becomes smaller with higher heat demand. In turn, a greater impact on EP is identified for distance.

**Figure 16.** The Yearly Economical Profitable as a function of M-TES price and distance between Geothermal Base Station and Building, for LCOH.OS = 0.1 EUR/kWh and for: (**a**) Total Heat Demand (THD) = 5000 kWh/year; (**b**) THD = 25,000 kWh/year.

#### *3.3. NPV Analysis*

The Net Present Value (NPV) was calculated (for 20 years and distance = 0.5 km) for M-TES.P = 6000 EUR and THD = 15,000 EUR/year (LCOH.G without M-TES.P = 0.045 EUR/kWh). The results for various discount rates are shown in Figure 17. At a discount rate of 6%, NPV reaches zero for LCOH.OS = 0.115 EUR/kWh, and for a discount rate of 0% at 0.082 EUR/kWh.

**Figure 17.** Net Present Value (NPV) as a function of LCOH.OS for different discount rate value, for distance = 0.5 km.

To verify the impact of changing the value of selected parameters (+/−20%) on the NPV change, a sensitivity analysis was performed—Table 1. The results of the analysis were sorted from the largest to the smallest impact, but this is related to the base level of the parameters—a different base level

would have a different impact on NPV (example distance 0.5 km and 5 km—Table 1). Significantly, the ratio of heat demand for heating and domestic hot water preparation (with the same THD) does not play any role in the NPV value, but only as shown in Section 2 for different frequencies of M-TES replacement in individual months. At the same time, the authors are aware that the results of the sensitivity analysis are affected by the adopted baseline values, e.g., as shown in Table 1, the assumed baseline distance will have a different impact on NPV change.


**Table 1.** NPV sensitivity analysis.

distance\*—only this example is calculated for base value 5 km, others for 0.5 km.

#### *3.4. Carbon Dioxide Emission Reduction*

The results in the form of quantitative CO2 emissions related to M-TES transport are presented in Table 2. The results in the form of a quantitative reduction of CO2 emissions in the event of the replacement of a coal or gas source by a geothermal source by a mobile heat storage are presented in Table 2. The results presented in Tables 2 and 3 were estimated for the assumptions presented in Section 2.3.

**Table 2.** Carbon dioxide emission of M-TES transport as a function of distance between building and M-TES load base, and THD. Source: own study.


Journeys related to transportation M-TES for a year to cover the heat demand for a building THD = 5000 kWh/year for a distance of 0.5 km are associated with CO2 emissions from the car at the level of 38 kg CO2/year, a distance of 10 km will result in 20x greater emissions. At THD = 25,000 kWh/year carbon dioxide emission is about 3.6 higher than in THD = 5000 kWh/year.

For a building with THD = 25,000 kWh/year and a distance of 0.5 km from the M-TES charging station, the use of geothermal source will bring:



**Table 3.** Carbon dioxide emission reduction as a function of distance, THD and replacement heat source. Source: own study.

#### **4. Conclusions**

This article considered, for Polish conditions, the technical and economic possibilities of providing geothermal heat to individual recipients using an M-TES system. For this purpose, heat storage, PCM and the best location regarding geothermal heat availability (temperature and the cost) were selected. The chosen location is Ba ´nska Ni ˙zna near Zakopane (southern part of the Poland). An indirect contact energy storage container filled with PCM characterized by melting temperature 70 ◦C and heat storage capacity 250 kJ/kg was selected, in the amount of 800 kg. The total M-TES heat capacity is 55 kWh.

The discussed PCM-container is satisfactory to supply a building with a total heat demand (central heating and domestic hot water) up to 25,000 kWh/year. This was established based on the data for the month with the lowest average temperature of the last 20 years (February 2012). Above this demand value, consideration should be given to the use of M-TES with larger heat capacity. The economic profitability of the M-TES system (with a price per warehouse of 6,000 EUR, i.e., a total of 12,000 EUR—two containers are needed) can be achieved for a heat demand of 5000 kWh/year with the price of a replaced heat source at the level of 0.21 EUR/kWh and distance between charging station and building (heat recipient) 0.5 km. For the heat demand of 15,000 kWh/year, the price for replaced heat reached EUR 0.11/kWh, and the same distance. In turn, for a demand of 25,000 kWh/year, the price of the replaced heat source reached 0.085 EUR/kWh.

The NPV value for the M-TES system was also determined. With a heat demand of 15,000 kWh/year, the NPV was calculated for 20 years of operation and a 6% discount rate that reached zero for a levelized cost of heating 0.115 EUR/kWh, and for a discount rate of 0% at 0.082 EUR/kWh. The levelized cost of heating in the case of geothermal energy without M-TES costs (2 × 6000 EUR) amounted to 0.045 EUR/kWh.

The economic profitability is significantly affected by the distance. For the adopted assumptions and at current prices, there are no arguments for transporting geothermal heat over distances bigger than 3–4 km away from the heat source, even in the case of the replaced heat price at the level of the current electricity price in Poland. This also concerns heat losses during the transportation, and ecological aspects. Considering the building with 25,000 kWh/year of total heat demand, which is located 0.5 km from the M-TES charging station, the use of a geothermal source can bring 9 Mg CO2 emission reduction in the case of replacing a coal-fueled heat source or 4.8 Mg CO2 emission reduction in the case of replacing a natural gas heat source. By increasing the distance to 10 km, a decrease in emission reduction of up to 6.4 Mg CO2 and 2.2 Mg CO2, respectively, can be observed.

**Author Contributions:** All the authors have contributed toward developing and implementing the ideas and concepts presented in the paper. All the authors have collaborated to obtain the results and have been involved in preparing the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by AGH—University of Science and Technology number 16.16.210.476.

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

#### **Abbreviation**


#### **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* **E**ff**ects of Cyclic Heating and Water Cooling on the Physical Characteristics of Granite**

#### **Xiangchao Shi 1,2,\*, Leiyu Gao 1, Jie Wu 1, Cheng Zhu 3, Shuai Chen <sup>1</sup> and Xiao Zhuo <sup>1</sup>**


Received: 21 January 2020; Accepted: 27 April 2020; Published: 29 April 2020

**Abstract:** This paper aims to study the effect of cyclic heating and flowing-water cooling conditions on the physical properties of granite. Ultrasonic tests, gas measured porosity, permeability, and microscope observations were conducted on granite after thermal treatment. The results showed that the velocity of P- and S-waves decreased as the number of thermal cycles increased. The porosity increased with the number of the thermal cycles attained at 600 ◦C, while no apparent changes were observed at 200 and 400 ◦C. The permeability increased with the increasing number of thermal cycles. Furthermore, microscope observations showed that degradation of the granite after thermal treatment was attributed to a large network of microcracks induced by thermal stress. As the number of thermal cycles increased, the number of transgranular microcracks gradually increased, as well as their length and width. The quantification of microcracks from cast thin section (CTS) images supported the visual observation.

**Keywords:** granite; physical characteristics; cyclic; thermal treatment; water cooling

#### **1. Introduction**

Geothermal energy is an important component of renewable energy, and most of the deep geothermal resource is stored in hot dry rock (HDR) [1]. HDR is defined as a hot and almost waterless geothermal system. Common HDR systems include granite or other crystalline basement rocks. Rock temperature varies from 150 to 500 ◦C at maximum depths of 5–6 km [2,3]. The use of deep HDR resources will contribute to the mitigation of the environmental pollution caused by traditional fossil energy [4]. Enhanced geothermal system (EGS) is an effective engineering method to exploit HDR resources [5,6]. To enhance the permeability of HDR, typical methods include increasing the heat transfer area and improving the efficiency of the hydraulic fracturing system. During the creation of the fracture network and the subsequent thermal energy exploitation, water is injected into the bedrock [7]. Studying the variation of rock physical properties after cyclic thermal treatment with flowing water cooling can provide a theoretical foundation for the EGS system.

A number of previous studies have conducted experiments to investigate the evolution of the physical and mechanical properties of the bedrock under temperature variations. Typical factors to consider include the heating temperature, heating rate, and cooling condition. With increasing temperature, rock experiences more significant deterioration. Hydraulic properties, such as porosity [8,9] and permeability [10,11], increase with a rise of the temperature, whereas mechanical properties, such as the P-wave velocity [11,12], unconfined compressive strength [7,12–15], tensile strength [11,16–18], and Young's modulus [9], decrease with increasing temperature [19]. The different types of microcracking also have been studied using 2D and 3D observations [20–23]. For example, thermal damage in Beishan granite subjected to high temperature treatment (from 100 to 800 ◦C at different heating rates, ranging from 1 to 15 ◦C/min) was studied in order to assess thermal effects on physical and mechanical properties. Results of acoustic emission (AE) monitoring, mechanical and physical properties measurements all indicated that heating rates had a significant impact on thermal damage. 5 ◦C/min was recognized as the critical heating rate for standard samples. Thermal stress induced by temperature gradients plays an important role in governing the damage in samples treated at a heating rate above 5 ◦C/min, at higher heating rates, thermal cracking is dominated by the stress concentrations caused by high thermal gradients [24,25]. Additionally, these physical and mechanical characteristics depend on different cooling conditions. Rock degradation caused by rapid cooling with water cooling or liquid nitrogen cooling was more severe compared with those under slower cooling such as air cooling, or cooling in a furnace [26–28]. Other rock types encountered in EGS projects, such as sandstone, have shown a similar trend [29–32]. The initial permeability of sandstone under certain pressure conditions was found to increased nonlinearly with the increase in temperature. Moreover, unconfined compressive strength and elastic parameters (i.e., elastic modulus and Poisson's ratio) of calcarenite decreased as the temperature was increased from 105 to 600 ◦C [33].

Many previous studies have investigated rock damage after a single heating and cooling treatment, however, few have considered the effect of cyclic treatments. Gräf studied the effects of cyclic thermal-heating treatment on the thermal expansion behavior of granite, however, the heating temperature and the number of treatment cycles were limited [34]. Mahmutoglu investigated the effect of thermal cycles on the mechanical behavior of marble and Buchberger sandstone. The results of unconfined compression, Brazilian and "Continuous Failure State" triaxial tests, pointed out that all of the mechanical parameters decreased gradually with an increasing number of heating cycles [35]. Additionally, Rong et al. studied the effect of thermal cycles on marble and granite subjected to air cooling on P-waves, stress–strain relationships, and acoustic emissions. The results revealed that thermal cyclic loading weakens the mechanical properties of the rock [36]. Furthermore, Wu et al. studied the effects of thermal cycles on the density, permeability, and unconfined compressive strength of granite subjected to liquid nitrogen cooling. Liquid nitrogen cooling was found to have a greater effect on the physical and mechanical properties of granite than air cooling [37]. Previous literature has also revealed that microcracks influence the physical and mechanical characteristics of the rock. Research has also quantified the microscopic responses of rock, particularly microcracking [38], to such processes.

Considering that water is the most common fluid used to extract thermal energy from HDR, in this experimental study, we investigate the effects of cyclic heating and water cooling on the physical and mechanical properties of granite, including a quantitative analysis of the resulting microcracks.

#### **2. Materials and Methodology**

#### *2.1. Rock Samples*

Granite was selected as the experimental material in our study. The samples were fine-grained granite with a grain size ranging from 0.5 to 1.5 mm (Chinese granite G655), which were collected from an outcrop located in Zhangzhou, Fujian, China (117.86◦ E, 24.83◦ N). No fissures were observed in the original rock. Two shapes of granite specimens were used: cylinders with the dimensions of 25 mm in diameter and 50 mm in length and discs with the dimensions of 25 mm in diameter and 10 mm in thickness (Figure 1). The mineralogy of the granite is described in Section 3.4.

**Figure 1.** Photographs of the rock specimens.

#### *2.2. Methodology of Heating and Cooling*

The granite cylinders were used for the measurement of porosity, permeability, P- and S-wave velocity measurements, whereas the granite disc samples were used for microstructural characterization. The P-wave and S-wave velocity of the cylindrical specimens were measured prior to thermal treatment. Thirty-one cylinders and thirty-one discs with similar P-wave velocities were divided into five groups with six specimens in each group (see Table 1). A SX-G04123 box-type electric furnace was used in heating (power 2.5 kW, maximum temperature 1200 ◦C). The process of heating a rock specimen to a pre-determined temperature and then cooling it down to room temperature with water was regarded as one thermal cycle. A heating rate of 1.5 ◦C/min was used to avoid the influence of thermal shock [39]. The pre-determined temperature, once reached, was kept constant for 5 minutes in the furnace to avoid the influence of subsequent heating time at a pre-determined temperature. The specimens were then taken out of the furnace and cooled down to room temperature (20 ◦C) with flowing water (shown in Figure 2). For each group, the pre-determined temperatures were set to 100–600 ◦C to mimic a high-temperature condition of deep bedrocks in an EGS system. The numbers of thermal cycles for each group was either 0 (i.e., no thermal treatment), 1, 2, 4, 8, and 16.


**Table 1.** Thermal treatment conditions of granite employed in this study. (A: disc specimens, B: cylindrical specimens).

**Figure 2.** Schematic diagram of the heating and cooling process.

#### *2.3. Ultrasonic Wave Velocity Measurements*

Before measuring the porosity and permeability, P and S-wave velocities of the rock specimens were measured after thermal treatment for different temperatures and numbers of thermal cycles were measured using an ultrasonic apparatus (HS-YS2A Type). Two transducers fixed by a holder applying the same force were positioned at the ends of the specimen. A pressure gauge and hand wheel were used to apply an equal force with transducers during every test (see Figure 3). Each specimen was tested twice to ensure repeatability.

**Figure 3.** Schematic of the ultrasonic test.

#### *2.4. Porosity and Gas Permeability Measurements*

We measured the porosity and steady-state permeability of the three temperature groups (200, 400, and 600 ◦C) using an SCMS-E high-temperature, high-pressure, multi-parameter core measurement system. The test procedure was performed in accordance with the methods suggested by the Practices for Core Analysis (GB/T 29172-2012) in Chinese. The measurement gas was nitrogen, while the test temperature was room temperature (10–15 ◦C), and the confining pressure was 3.5 MPa (Figure 4).

**Figure 4.** Schematic of the porosity and permeability measurement system.

#### *2.5. Microscopic Observation*

Cast thin sections (CTS), i.e., thin sections of rock impregnated with colored epoxy, were used for highlighting the microcracks. The disc rock samples were saturated with blue epoxy to distinguish pores and fracture from rock matrices. The procedure of obtaining CTS is shown in Figure 5. All disc rock specimens were impregnated with colored epoxy after thermal treatment. After leveling, lapping and polishing, a thin section of size 25 × 0.03 mm was obtained. The impregnated sample was bonded to the surface of a piece of glass for further processing. The thin section image was then photographed under plane polarized light and cross polarized light by using a polarizing microscope (Zeiss Scope A1) with an attached digital camera.

**Figure 5.** Schematic of polarizing microscope observation process.

#### *2.6. Image Processing*

Microcracks appear blue under plane-polarized light because the blue epoxy filled the microcracks. We selected the color of blue epoxy (R: 70, G: 132, B: 171). Usually, it is not consistent in every CTS image. Then, we set the tolerance to 70 to ensure that the blue parts of the image were successfully selected. The results are shown in Figure 6b. Manual interactive thresholding segmentation was used for the segmentation process. The thresholding of 40 was applied to the intermediate image based on the image histogram (see Figure 7). The binarized image result is shown in Figure 6d, which was further used for quantitative analysis (Figure 6c).

**Figure 6.** Schematic view of image processing. (The blue part is the preparation process, while the green part is the measurement process).

**Figure 7.** Segmentation via image thresholding.

We then measured four different parameters to measure: area, size (width), size (length), dendrites—one pixel thick open branches (Figure 8). The number of isolated elements was automatically counted (Figure 6d). It should be noted that the size width was not the actual width of the fracture. Finally, the image statistical data were exported to a worksheet for post processing. The image porosity was calculated according to Equation (1).

$$Porosity = \frac{\text{White pixels}}{\text{All pixels}} \times 100\% \tag{1}$$

**Figure 8.** Three different types of parameters for measurements.

#### **3. Results**

#### *3.1. Porosity*

As shown in Figure 9, high temperatures (600 ◦C) had a greater effect on porosity than low heating temperatures (200 and 400 ◦C). Moreover, at 600 ◦C, porosity increased substantially as the cycle number increased from zero to two. At 200 and 400 ◦C, the influence of the number of thermal cycles on the porosity and permeability is negligible.

**Figure 9.** Effect of the number of thermal cycles and temperature on granite porosity.

#### *3.2. Gas Permeability*

The permeability variations with the number of thermal cycles at different heating are plotted in Figure 10. The trends were similar for all three different heating temperatures. A positive correlation was observed between thermal cycling and permeability increase at 400 and 600 ◦C, while there was an irregular rise against thermal cycles at 200 ◦C. In addition, the permeability of granite at 600 ◦C was significantly higher than that at 400 ◦C and 200 ◦C. At 600 ◦C, the permeability of granite increased from 0.0001 to 4.7770 mD after 16 thermal cycles.

**Figure 10.** Permeability characteristics after thermal treatment. (**a**) Three different temperatures using same y axis. (**b**) The pre-determined temperature was 200 ◦C. (**c**) The pre-determined temperature was 400 ◦C. (**d**) The pre-determined temperature was 600 ◦C.

#### *3.3. P- and S-Wave Velocity*

Ultrasonic wave test is commonly used to detect the interior failure in a rock because of its simple and non-destructive characteristics. The typical P- and S-wave forms recorded and used to determine the velocities are presented in Figure 11 and the velocity results are shown in Figure 12. It appears that both P-wave velocity and S-wave velocity exhibited a similar negative correlation with heating temperature and the number of thermal cycles. The gradient of both P-wave and S-wave velocity decreased as the number of cycles increased. At 600 ◦C and one cycle of the thermal treatment, P- and S-wave velocities decreased by 73.6% and 58.6%, respectively. At 600 ◦C, after 16 cycles of the thermal treatment, the velocity reduced by 84.3% and 82.4%, respectively.

**Figure 11.** Typical P- and S wave forms recorded and used to determine velocities. (**a**) P-wave form. (**b**) S-wave form.

**Figure 12.** Variation of P-wave velocity and S-wave velocity changes with the number of thermal cycles and temperature. (**a**) Variation of P-wave velocity with temperature. (**b**) Variation of S-wave velocity with temperature. (**c**) Variation of P-wave velocity with the number of thermal cycles. (**d**) Variation of S-wave velocity with the numbers of thermal cycles.

#### *3.4. Rock Microstructure*

The microscopic observations of A31 (no thermal treatment) are shown in Figure 13. The granite is mainly composed of feldspar, quartz, and biotite with a small amount of pyroxene and magnetite. Anorthosite is 578 μm in length and 251 μm in width. Quartz is 458 μm in length and 243 μm in width. Clear boundaries were observed between mineral grains (Figure 13b). No blue epoxy was not observed in the plane-polarized image (Figure 13a), indicating that the granite has negligible porosity.

**Figure 13.** Polarized micrographs of granite without thermal treatment. (**a**) Plane polarized image; (**b**) cross polarized image. Qtz—Quartz; An—Anorthite; Bt—Biotite; Px—Pyroxene; Mag—Magnetite.

Figure 14 presents the CTS observation results for the granite after one thermal cycle at 600 ◦C. The mineral grains and the associated microcracks are clearly affected by thermal treatment. Two types of microcracks were observed: grain boundary microcracks, transgranular microcracks (including intracrystalline microcracks) [12]. Grain boundary microcracks describe the microcracks developed at the boundary between different minerals, whereas transgranular microcracks are those developed in the interior of the mineral grains. These microcracks were predominantly confined within minerals, with some passing through multiple grains. Occasionally, transgranular microcracks and grain boundary microcracks developed connectivity at an angle of approximately 90◦. Furthermore, the location of the microcracks was related to mineral species. Transgranular microcracks were typically developed in feldspar, and the development direction was almost perpendicular to the optical twin crystal direction of feldspar. Grain boundary microcracks were typically between feldspar and quartz, and between feldspar and feldspar [25,40]. Almost no transgranular microcracks were observed in biotite.

**Figure 14.** The types and location of granite fracture. (**a**) Plane polarized image; (**b**) cross polarized image. Qt—Quartz; An—Anorthose; Bt—Biotite; Px—Pyroxene. 1—grain boundary microcracks; 2—transgranular microcracks. (600 ◦C, 1 cycle).

Plane-polarized and cross-polarized CTS images of the granite specimens with different heating temperatures (100–600 ◦C, one cycle) are presented in Figures 15 and 16, respectively. Heating temperature affected the microcracks development behaviors. In the range of 100–300 ◦C, minimal blue epoxy was observed in the images. Grain boundary microcracks were sparsely distributed in a specimen heated to 400 ◦C (see Figure 15d), however, as the temperature exceeded 400 ◦C, more grain microcracks could be observed in the granite. For the specimen treated at 500 ◦C heating temperature (see Figures 15e and 16e), transgranular microcracks were observed in feldspar. Grain boundary microcracks also existed but in a smaller proportion than those transgranular microcracks in the granite subjected to 600 ◦C. According to the plane-polarized microscope observation, as shown in Figure 15f (600 ◦C), abundant grain boundary microcracks were found along feldspar and quartz grains. Moreover, the transgranular microcracks and grain boundary microcracks occasionally coalesced.

Thermal cycles also had a significant influence on microcrack evolution, especially that of transgranular. The plane-polarized and cross-polarized photographs of CTS for granite specimens treated at 600 ◦C with different numbers of thermal cycles are presented in Figures 17 and 18. With the increase of thermal cycles, both the length and width of transgranular microcracks increased, as well as the number of grain boundary microcracks and transgranular microcracks. Figures 17f and 18f present sample A30, which was subject to 600 ◦C heating and 16 thermal cycles. The maximal width of transgranular reached approximately 20 μm. Microcracks were well developed and most mineral boundaries had grain boundary microcracks. The width of the transgranular microcracks was larger than that of grain boundary microcracks.

**Figure 15.** Plane polarized images of microcracks in granites subjected to different heating temperatures. (**a**) 100 ◦C; (**b**) 100 ◦C; (**c**) 200 ◦C; (**d**) 400 ◦C; (**e**) 500 ◦C; (**f**) 600 ◦C; 1—grain boundary microcracks; 2—transgranular microcracks.

**Figure 16.** Cross polarized images of microcracks for granites subjected to different heating temperatures. (**a**) 100 ◦C; (**b**) 100 ◦C; (**c**) 200 ◦C; (**d**) 400 ◦C; (**e**) 500 ◦C; (**f**) 600 ◦C. Qtz—Quartz; An—Anorthose; Bt—Biotite; Px—Pyroxene.

**Figure 17.** Plane polarized images of microcracks for granites with 600 ◦C heating temperature subjected to different numbers of the thermal cycle. (**a**) No thermal treatment; (**b**) 1 cycle; (**c**) 2 cycles; (**d**) 4 cycles; (**e**) 8 cycles; (**f**) 16 cycles. 1—grain boundary microcracks; 2—transgranular microcracks.

**Figure 18.** Cross polarized images of microcracks for granites with 600 ◦C heating temperature subjected to different numbers of the thermal cycle. (**a**) No thermal treatment; (**b**) 1 cycle; (**c**) 2 cycles; (**d**) 4 cycles; (**e**) 8 cycles; (**f**) 16 cycles. Qtz—Quartz; An—Anorthose; Bt—Biotite; Px—Pyroxene.

#### *3.5. Microcrack Morphology*

The granite morphology was analyzed through the CTS images (Figure 19), which were sorted according to the microcrack area and arranged vertically. The length of microcracks and the number of inflexion points increased with the number of thermal cycles (see Figure 19). This implies that the development and cross-cutting of grain-boundary microcracks and transgranular microcracks developed and crossed together. Hence, thermal cycles had a substantial influence on high-temperature granite subjected to water cooling.

**Figure 19.** The morphology for granite with 600 ◦C heating temperature and different thermal cycles. (**a**) 1 cycle; (**b**) 2 cycles; (**c**) 4 cycles; (**d**) 8 cycles; (**e**) 16 cycles.

We conducted a statistical analysis of the microcracks area to explore the distribution of pores size. As shown in Figure 20, relative frequency of microcracks descended rapidly from 0.35 to 0.05. 80% of pores contained less than 21 pixels or 17.01 <sup>μ</sup>m2 (1 pixel <sup>≈</sup> 0.81 <sup>μ</sup>m2). Tiny pores were the most numerous among all microcrack sizes.

**Figure 20.** Distribution of microcracks with 600 ◦C heating temperature and 8 thermal cycles.

The results of the CTS image processing are presented in Figure 21. Generally, differences were observed between the porosity measured via gas permeability measurement and microscope observations. The latter increased by 3.5 times from one cycle to 16 cycles, which rose from 1.61% to 5.67%. Conversely, the former only increased from 1.73% to 3.42% (Figure 21a). Measurements of maximum porosity, maximum length, and maximum width tended to increase with a greater number of thermal cycles (see Figure 21b–d). The line charts revealed that the development of microcracks induced by thermal stress increased with the rise of thermal cycles. The number of measured dendrites is shown in Figure 22, which also increased with the number of thermal cycles, by almost 13 times to 63. This clearly illustrates the expansion of microcracks with the number of thermal cycles at high heating temperatures.

**Figure 21.** Results of the cast thin section (CTS) image processing versus the numbers of the thermal cycle. (**a**) Variation of image porosity. (**b**) Variation maximum image porosity. (**c**) Variation maximum fracture length. (**d**) Variation maximum fracture width.

**Figure 22.** Effect of the number of thermal cycles on the number of dendrites.

#### **4. Discussion**

Porosity and crack density are the most important properties, playing a major role in the structural integrity of a rock [30]. They are related to the type and arrangement of mineral grains, internal pores, and microcracks. High temperature thermal treatment could increase the porosity of granite. Some researchers reported that 400 ◦C is a threshold temperature for granite to change its structure [9,30,41]. Below the threshold temperature, no significant relationship was found between the porosity of granite and the number of thermal cycles, suggesting that the increasing number of thermal cycles does not contribute to the propagation of microcracks. However, beyond a temperature of 600 ◦C, the increase of cyclic thermal treatment is associated with the increase of rock porosity.

Permeability is in relation to pores and especially fractures. Generally, fractures play an important role in the percolation capacity. Obviously, the thermal treatments can enhance the permeability of granite. The permeability is significantly improved as heating temperature rises [11,17,42]. Interestingly, the porosity at 400 ◦C was lower than that at 200 ◦C except after two thermal cycles, whereas the permeability exhibited the exact opposite trend. This reflected that the enhanced connectivity between granite microcracks was enhanced. Pore and crack networks facilitated fluid flow through the specimens, leading to the increased permeability. As a result, it demonstrated reversely that propagation of microcracks increased with the number of thermal cycles. The number of thermal cycles also had significant effects on the permeability of granite.

The CTS image processing employed in this study exhibits some limitations. First of all, the images cannot reveal all microcracks in the sample because of their finite resolution and size. Secondly, the images are 2D images rather than 3D images; thus, some microcracks cannot be observed for technical reasons. These two limitations would result in an underestimation of the porosity. Third, due to the inconsistent blue color of epoxy, some noise pollution is inevitable in the CTS image, which would lead to an overestimated porosity. Nevertheless, these images still provide useful and quantitative descriptions of microcracks in the granite.

During the process of heating and cooling, a series of physical and chemical reactions have occurred in granite (see Figure 23). When the heating temperature exceeded 100 ◦C, water inclusions that originally existed in the granite pores due to wettability and capillary force-escaped rock in the form of gas [8]. As a result, the P-wave velocity decreased. The chemical formula for biotite is *K* (Mg<0.67, Fe>0.33)3[AlSi3*O*13](OH)<sup>2</sup> . In the temperature range of 300–500 ◦C, due to the escape of crystal water and dissociation of the H<sup>+</sup> and OH<sup>−</sup> (the existing form of constitution water in the mineral crystal lattice structure), the mineral framework was destroyed and microcracks developed in rock [17,42]. When the heating temperature reaches 573 ◦C, low-temperature α quartz turns into high-temperature β quartz, which is accompanied by a sudden volume expansion [7,43–45]. This transition results in severe derogation of the granite. During the heating process, due to differences in the thermal expansion coefficients of different minerals, thermal stresses accumulate at granular interfaces, even if the temperature field outside is uniform. This expansion mismatch primarily contributes to the development of grain boundary microcracks [24,46]. During the cooling process, granite specimens are placed in the flowing water, which induces the sudden variation of temperature and results in the gradual formation of grain boundary microcracks and transgranular microcracks among rock minerals. Water then invades the connected microcracks, resulting in new chemical reactions in the minerals. Because the transition of quartz to β quartz is a reciprocal reaction, the granite specimens experience repeated damage with an increasing number of thermal cycles.

**Figure 23.** Mechanisms of thermal damage in heating and cooling.

#### **5. Conclusions**

We conducted a set of physical experiments to investigate the effect of cyclic thermal treatment and water cooling on the physical characteristics of granite. The results show that the thermal cycling has a significant influence on the physical characteristics (i.e., porosity, permeability, the seismic velocity). The results contribute to the fundamental understanding of the evolution of porosity and permeability in HDR geothermal systems. Qualitative and quantitative analyses of CTS images led to the following conclusions:


**Author Contributions:** Conceptualization, X.S. and C.Z.; data curation, J.W. and L.G.; formal analysis, L.G. and S.C.; supervision, C.Z. and L.G.; validation, X.Z. and X.S.; writing—original draft, L.G., X.Z. and X.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Science Foundation of China, grant number 51774428, and Sichuan International Science and Technology Innovation Cooperation/Hong Kong, Macao and the Taiwan Science and Technology Innovation Cooperation, grant number 2019YFH0166, and the National Students' Innovation and Entrepreneurship Training Program of China, grant number 20180615013.

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

#### **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*
