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

: 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 ﬂow 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) × 10 20 J (mean 8.6 × 10 20 J), (5.1–16) × 10 20 J (mean 9.1 × 10 20 J), and (3.2–11) × 10 20 J (mean 6.6 × 10 20 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) × 10 21 J (mean 12 × 10 21 J). The results indicate capacities of a range value of (1.2–1.7) × 10 21 J (mean 1.4 × 10 22 J) for the whole NJB (<10 km).


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][2][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 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][16][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][19][20][21].  [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; Pt 2 , 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][16][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][19][20][21].

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 −2 (Yandong-3) and 110 mW·m −2 (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.

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 (E 2 d), the Sanduo Formation (E 3 s), 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.

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  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][27][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.

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)).
where Q Total , Q W , and Q R represent the total geothermal energy and that stored in pore water and rock, J; A (m 2 ) and D (m) are the area and thickness of the thermal reservoirs; ρ R and ρ W denote the density of the geothermal reservoir and geothermal water, kg·m −3 ; c R and c W are the specific heat of the geothermal reservoir rock and geothermal water, J·kg −1 · • C −1 ; ϕ is rock porosity, %; T and T 0 represent the average temperature of geothermal reservoir and surface, • C.

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 (T 0 ) 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 (c R ) and water (c W ), 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.
cording 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.

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. ment 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. 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 −1 . Vertically, the temperature gradient in the NJB increases with depth, and the temperature gradients 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 −1 . Vertically, the temperature gradient in the NJB increases with depth, and the temperature gradients of Ny and E 2 s-E 2 d are generally lower than that of E 1 f-K 2 t. Formations of the Neogene and above are characterized by a low temperature gradient, generally between 15-30 • C·km −1 , which may be the result of heat redistribution and groundwater activity in the shallow. The high sand content in the Ny, E 2 s, and E 2 d 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 −1 .

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). 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 −1 and the largest value being the Silurian S 2 m 1 with a thermal conductivity of 8.0 W·m −1 ·K −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 −1 ), especially the Late Paleozoic-Quaternary thermal conductivity is mostly less than 2.0 W·m −1 ·K −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 −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 −1 , all of which can be used as good thermal storage layers.

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 (Q s ) 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 Q s , surface temperature T 0 , thermal conductivity λ, and radiogenic heat production A is defined as T(z) = T 0 + z Q s /λ − z 2 A/(2λ) (4) where λ (W·m −1 ·K −1 ) and A (µW·m −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 −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).      Figure 6b shows the potential of the geothermal resource in which the value changes from 5.0 × 10 20 J to 1.6 × 10 21 J with a mean value of (9.1 ± 0.2) × 10 20 J), with nearly 7% probability of the most likely value of 0.86 × 10 21 J, and a probability of < 5% that the geothermal resources > 1.2 × 10 21 J or < 0.67 × 10 21 J. Dividing by area, it indicates a geothermal potential of 4.5 × 10 16 J·km −2 . Figure 6c presents a geothermal resource potential for the Dai'nan thermal reservoir of 3.2 × 10 20 J to 1.1 × 10 21 J with a mean value of (6.6 ± 0.1) × 10 20 J), with a 90 % probability that the energy potential is (0.45-0.90) × 10 21 J. Considering the geothermal resource potential per unit of thermal storage area, it indicates a potential of 6.7 × 10 16 J·km −2 . Figure 6d indicates that the geothermal resource potential for the carbonate thermal reservoir sums to 9.9 × 10 21 J to 1.5 × 10 22 J with a mean value of (1.2 ± 0.1) × 10 22 J). From Figure 6d, we can get that the geothermal potential is above 1.3 × 10 22 J or below 1.1 × 10 22 J in which the probability is less than 5 %, and the most likely value of the whole is 1.2 × 10 22 J (nearly 8%). The geothermal resource potential per unit of the thermal storage area is 3.7 × 10 17 J·km −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.

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 × 10 22 J to 1.7 × 10 22 J (mean value (1.4 ± 0.1) × 10 22 J), with the highest probability being 1.4-1.5 × 10 22 J (probability nearly 7%), and a 90 % probability of a range of (1.3-1.6) × 10 22 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 6b shows the potential of the geothermal resource in which the value changes from 5.0 × 10 20 J to 1.6 × 10 21 J with a mean value of (9.1 ± 0.2) × 10 20 J), with nearly 7% probability of the most likely value of 0.86 × 10 21 J, and a probability of < 5% that the geothermal resources > 1.2 × 10 21 J or < 0.67 × 10 21 J. Dividing by area, it indicates a geothermal potential of 4.5 × 10 16 J·km −2 . Figure 6c presents a geothermal resource potential for the Dai'nan thermal reservoir of 3.2 × 10 20 J to 1.1 × 10 21 J with a mean value of (6.6 ± 0.1) × 10 20 J), with a 90 % probability that the energy potential is (0.45-0.90) × 10 21 J. Considering the geothermal resource potential per unit of thermal storage area, it indicates a potential of 6.7 × 10 16 J·km −2 . Figure 6d indicates that the geothermal resource potential for the carbonate thermal reservoir sums to 9.9 × 10 21 J to 1.5 × 10 22 J with a mean value of (1.2 ± 0.1) × 10 22 J). From Figure 6d, we can get that the geothermal potential is above 1.3 × 10 22 J or below 1.1 × 10 22 J in which the probability is less than 5 %, and the most likely value of the whole is 1.2 × 10 22 J (nearly 8%). The geothermal resource potential per unit of the thermal storage area is 3.7 × 10 17 J·km −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.

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 × 10 22 J to 1.7 × 10 22 J (mean value (1.4 ± 0.1) × 10 22 J), with the highest probability being 1.4-1.5 × 10 22 J (probability nearly 7%), and a 90 % probability of a range of (1.3-1.6) × 10 22 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. Energies 2021, 14, x FOR PEER REVIEW 16 of 18 Figure 7. Geothermal resource potential simulation results for the sum of the four thermal reservoirs.

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) × 10 22 J, with a mean value of (1.4 ± 0.1) × 10 22 J and a most likely value of 1.4 × 10 22 J, which is equivalent to 4.9 × 10 13 tons of standard coal heat content. Moreover, the most promising carbonate thermal reservoir owns an average geothermal resource potential of (1.2 ± 0.1) × 10 22 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.

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].