Indirect Methods for Validating Shallow Geothermal Potential Using Advanced Laboratory Measurements from a Regional to Local Scale—A Case Study from Poland

: This paper presents a broad overview of laboratory methods for measuring thermal properties and petrophysical parameters of carbonate rocks, and analytical methods for interpreting the obtained data. The investigation was conducted on carbonate rock samples from the Krak ó w region of Poland in the context of shallow geothermal potential assessment. The measurement techniques used included standard macroscopic examinations; petrophysical investigations (porosity, density); analysis of mineral composition thermal conductivity (TC) and speciﬁc heat measurements; and advanced investigations with the use of computed tomography (CT). Various mathematical models, such as layer model, geometric mean, and spherical and non-spherical inclusion models, were applied to calculate thermal conductivity based on mineralogy and porosity. The aim of this paper was to indicate the optimal set of laboratory measurements of carbonate rock samples ensuring su ﬃ cient characterization of petrophysical and thermal rock properties. This concerns both the parameters directly characterizing the geothermal potential (thermal conductivity) and other petrophysical parameters, e.g., porosity and mineral composition. Determining the quantitative relationship between these parameters can be of key importance in the case of a shortage of archival thermal conductivity data, which, unlike other petrophysical measurements, are not commonly collected. The results clearly show that the best correlations between calculated and measured TC values exist for the subgroup of samples of porosity higher than 4%. TC evaluation based on porosity and mineral composition correlation models gives satisfactory results compared with direct TC measurements. The methods and results can be used to update the existing 3D parametric models and geothermal potential maps, and for the preliminary assessment of geothermal potential in the surrounding area.


Introduction
Energy policies and climate protection are, in addition to safety issues, some of the most important global challenges at the present time. One of the main elements of climate protection is energy transformation, particularly in the fast-developing countries of Asia (China, India, etc.) and East Europe, where the share of fossil fuel, comprised mainly of coal and natural gas, remains dominant in particular industries and heating.
Here, a wide range of investigations and laboratory measurements were conducted, ranging from standard mineralogical analyses and petrophysical measurements of density and porosity, and thermal conductivity and heat capacity, to advanced methods using computed tomography techniques. The obtained data were applied to calculate thermal conductivity with the use of different mathematical models. Determining the thermal conductivity based on other parameters can be of key importance in the case of a shortage of archival thermal conductivity analyses, which, unlike other petrophysical measurements, are not commonly performed.
The conducted investigations and measurements will both extend the scant base of the petrophysical and thermal parameters of the rocks of the Kraków region, and broaden the knowledge on applicable measurement methods and eventual limitations in their application. A second aim was to indicate the appropriate methods for performing laboratory measurements and data analysis that can be applied to carbonate rocks in other locations.

Theoretical Background
Thermal conductivity depends on petrophysical parameters such as mineral composition, porosity, grain size, degree of cementation, size and shape of pores, the presence of fractures and cavities, and pressure and temperature. These parameters should be considered in the analysis of the obtained laboratory thermal conductivity values. The most important parameters, with a brief description of their possible impact on thermal conductivity, are listed below.
Energies 2020, 13, 5515 5 of 32 Table 1. Thermal conductivity values of chosen minerals according to [27,28].  [27] Porosity-the thermal conductivity coefficient value decreases with an increase in porosity because both the air and media saturating the pore space are characterized by a considerably lower thermal conductivity than the rock framework. The thermal conductivity of water is 0.61 W/m·K; oil, 0.14 W/m·K; and air, 0.026 W/m·K [12,29].

Mineral
Fractures and cavities-these can be of great significance in the case of igneous, metamorphic, and carbonate rocks. The thermal conductivity value is influenced by the number of fractures in addition to their geometry, and by the properties of the filling substances [12]. The presence of fractures is also connected with a higher vulnerability of rock to pressure changes. With an increase in pressure, the fractures close and the thermal conductivity values significantly increase [12].
Grain size-the thermal conductivity value increases with the size of the grain. This is due to the decreased contact between grains and, subsequently, lower resistance between grain contacts [30,31].
Rock structure-anisotropy influences the heat flow rate; thus, in the case of rocks with an orientated structure, the thermal conductivity value depends on the measurement direction [12,13,29,32]. The anisotropy of thermal conductivity is determined by the coefficient Kλ, which is defined as the ratio of the λ value measured parallel to rock structure (λŞ) to λ values measured perpendicularly (λ⊥). The highest Kλ values are characteristic of shales, schists, and gneisses, and the lowest for carbonate rocks [29,32].
Pressure-better heat flow in grain contacts, tightening fractures, and decreasing porosity caused by increasing pressure result in higher thermal conductivity [12]. Experiments conducted on mudstone and dolomite have shown that the thermal conductivity increases with an increase in pressure up to circa 80 MPa and then stabilizes [33].
Temperature-thermal conductivity depends on temperature. This relationship is connected to the structure of the material. With increasing temperature, the thermal conductivity of crystalline materials (such as minerals) decreases, and increases in the case of amorphous materials. For the majority of rocks, the thermal conductivity decreases as temperature increases [12].
To summarize, the thermal conductivity value of a rock is a function of both the content and the thermal conductivity of the rock's minerals and media saturating the pore space. It also depends on the pore space structure. Thus, when assessing this parameter, different models can be applied based on a simplification of the rock's internal geometry (structure), allowing for calculation of the rock thermal conductivity based on the properties of its components [12][13][14][15][16][17][18][19][20][21][22].
To determine thermal conductivity, various mathematical models have been applied, from simple layer models, to more sophisticated spherical and non-spherical inclusion models.

Layer Models
The layer models in which the heat flow is parallel or perpendicular to the boundary between components is assumed to determine the extreme values, i.e., limits within which the true thermal conductivity values are placed [12,17,20].
The thermal conductivity for the heat flow parallel to the boundary between components is defined by the arithmetic mean (λ_ aryt ), thus determining the upper limit of the investigated value [12,17,20]: where: λ m -thermal conductivity of the grain framework [W/m·K], λ f -thermal conductivity of the pore media [W/m·K], φporosity [%]. The thermal conductivity for the heat flow perpendicular to the boundary between components is determined by the harmonic mean (λ_ harm ), and defines the lower limit of the investigated value [12,17,20]: (2) where: λ m -thermal conductivity of the grain framework [W/m·K], λ f -thermal conductivity of the pore media [W/m·K], φporosity [%].

Spherical Inclusion Models
Clausius-Mossotti spherical inclusion models [12] for rocks consisting of spherical inclusions of thermal conductivity λ 1 , dispersed in host material of thermal conductivity λ 2 , are determined by the following equation [12]: where: V-volume fraction of the inclusions, λ-thermal conductivity of the whole rock [W/m·K], λ 1 -thermal conductivity of the inclusion material [W/m·K], λ 2 -thermal conductivity of the host material [W/m·K]. When the basic component is the grain framework of thermal conductivity λ m , and pore solutions of thermal conductivity λ p appear as spherical inclusions, the equation takes on the following form [12]: where: η = λ m /λ f , λ m -thermal conductivity of the grain framework, λ f -thermal conductivity of the pore media [W/m·K], φ-porosity.
For rock consisting of spherical grains suspended in a fluid, the following equation was derived [12]: where: Energies 2020, 13, 5515 7 of 32 λ f -thermal conductivity of the pore media, φporosity.
To determine the thermal conductivity of a rock based on the thermal conductivity of individual components, empirical models using the geometric mean have often been applied [16][17][18]: where: n-number of the components, λ i -thermal conductivity of the i-th constituent of the rock, V i -fractional volume of the i-th constituent of the rock.

Non-Spherical Inclusion Models
Non-spherical inclusion models describing rocks consisting of a grain matrix with spheroidal, non-intersecting pores are used for the simulation of mechanical and acoustic features [20]. A spheroid (rotating ellipsoid) is characterized by two equal axes. The parameter defining the shape of the rotating ellipsoid is the aspect ratio α, which determines the ratio of the length of the unequal axis to that of one of the equal axes. In borderline cases, the rotating ellipsoid can assume the shape of either a needle (α→∞), a sphere (α→1), or an extensive, flattened disk (α→0) [12,20]. In this work, a non-spherical inclusion model based on the generalization of a spherical inclusion Clausius-Mossotti model [12] was presented: where: λ m -thermal conductivity of the grain framework [W/m·K], λ f -thermal conductivity of the pore media [W/m·K], φporosity [%]. Parameter R mi , which represents the function of depolarization coefficients along the ellipsoid's axes L a , L b , and L c is expressed by the following equation [34]: The values of both the depolarization coefficients and the parameter R mi for chosen pore shapes (borderline cases) [12,34] are presented in Table 2. Table 2. Depolarization coefficients and parameter R mi values for pores of the needle, the sphere, and the disk shape [12,34].

Pore Shape Depolarization Coefficients Parameter R mi
Sphere

Geological Setting of the Pilot Area
Kraków is situated in the southern part of Poland, in the Vistula valley, at a distance of several kilometers from the Carpathian thrust and about 100 km from the Tatra Mountains, at 219 m above sea level. From a geological perspective, the Kraków region is situated on the border of three extensive Energies 2020, 13, 5515 8 of 32 geological units: the Silesia-Kraków Monocline, Nida Basin, and Carpathian Foredeep, and its geological structure is complex. The border between the Silesia-Kraków Monocline and the Nida Basin runs within the so-called Ojców plate and conventionally has been assumed to be in accordance with the course of the outcrops of Cretaceous formations. From the south, the Ojców plate border runs along an extensive zone of tectonic horsts. This area-to the south of the Ojców plate up to the Carpathian thrust-is part of the Carpathian Foredeep [35].
Two large rock complexes, so called tectonic structural stages can be distinguished in the Kraków region. The older complex includes Devonian, Carboniferous, and older deposits that were tectonically deformed during Variscan orogenesis [36]. The younger complex, built from Permian, Triassic, Jurassic, and Cretaceous deposits, is the so-called Mesozoic-Permian stage. The rocks of both of the complexes became inclined towards the north east-probably during the Laramic phase, between the Cretaceous and the Tertiary-which caused the formation of a monoclinal structure [35,37]. Tectonic modeling of the Silesia-Kraków Monocline took place in stages (mainly in the Neogene) and caused the emergence of horsts and depressions.
The geological structure of the town center is of similar character, and horsts creating distinct elevations are visible in the land morphology. Historically, some of these were exploited to situate various edifices.
Neogene tectonics in the Kraków region are connected mainly with phases of rocks forming movements in the nearby Carpathians. Here, faults span several generations and their detailed dating is not always possible [35,38]. Rocks appearing within the faults are mainly Upper Jurassic limestone and, less frequently, Upper Cretaceous, whereas depressions are filled with Miocene clay deposits including, locally, evaporates (clays with gypsum, gypsum rocks). These formations belong to the Carpathian Foredeep, with the smallest width in Poland (ca. 10-15 km) in this area, and its northern border shows an erosive character [39].
The youngest formations in the geological profile of the Kraków region are represented by the Quaternary formations (Pleistocene and Holocene), which mainly fill the paleovalleys of the Vistula and its tributaries, in addition to other morphological depressions. Pleistocene formations connected with glaciations are represented by glacial tills, fluvioglacial and alluvial sands, and gravels and loesses. In contrast, younger Holocene formations create a series of terraces, mainly in the Vistula and Rudawa valleys, and are represented by sands, gravels, and alluvial soils [40]. In the youngest Quaternary, so-called Anthropocene formations appear mainly as banks connected to human settlements, and-in the center of Kraków-to the town's historical layers [40,41]. In Kraków, deposits of minerals have been found, e.g., natural aggregates (sand and gravels), clays for building ceramics, limestones and marls for the limestone industry, and deposits of mineral and healing waters.
Kraków is located on the upper course of the Vistula river (the Baltic Sea reception basin). In the Kraków region, underground waters are connected with rocks of the Paleozoic (Devonian and Carboniferous), Jurassic, Cretaceous, Miocene, Eocene, and Quaternary stages [40,41]. The location of the investigated boreholes on a geological map of the Kraków region without Quaternary and terrestrial Tertiary deposits is shown in Figure 2.

Materials and Methods
The investigations were conducted on 22 rock samples collected from four boreholes in the Kraków region: Kopiec-1G, Kopiec-4G, Opatkowice-1, and Trojanowice-2 (Figures 2 and 3).   Rock samples were collected during the realization of the GeoPLASMA-CE "Shallow Geothermal Energy Planning, Assessment and Mapping Strategies in Central Europe" project carried out in 2016-2019. The project aimed at encouraging shallow geothermal use in heating and cooling strategies in central Europe. Within the GeoPLASMA-CE project, only thermal conductivity measurements of dry samples were conducted. The methods and the range of the measurements undertaken in the current study exceeded the scope of the above-mentioned project. The thermal conductivity values for both dry and saturated samples, and the quantitative mineral composition and porosity, were determined for all rock samples. To address the large heterogeneity of the pore space (unequal distribution of vugs and fractures), computed tomography was used for the selection of samples suitable for particular investigations ( Figure 4). Rock samples were collected during the realization of the GeoPLASMA-CE "Shallow Geothermal Energy Planning, Assessment and Mapping Strategies in Central Europe" project carried out in 2016-2019. The project aimed at encouraging shallow geothermal use in heating and cooling strategies in central Europe. Within the GeoPLASMA-CE project, only thermal conductivity measurements of dry samples were conducted. The methods and the range of the measurements undertaken in the current study exceeded the scope of the above-mentioned project. The thermal conductivity values for both dry and saturated samples, and the quantitative mineral composition and porosity, were determined for all rock samples. To address the large heterogeneity of the pore space (unequal distribution of vugs and fractures), computed tomography was used for the selection of samples suitable for particular investigations (Figure 4).  Several laboratory methods allow for the determination of thermal conductivity values. These can be divided into two groups: steady-state methods and transient techniques. Selection of a laboratory method for thermal conductivity determination should be based on the characteristics of the tested material, namely thermophysical properties of the sample, and the quantity and form of the supplied material, i.e., powdery or solid sample. Limitations of each of the methods should also be taken into consideration. Transient methods make it possible to obtain results in a shorter time, and steady-state methods require a longer measurement period and a large sample [43]. Due to frequent unequal distribution of vugs and fractures in the pore space, it is difficult to determine the thermal conductivity of carbonate rocks. Hence, if the pore space heterogeneity is to be considered, the measurements of such rocks should be conducted using samples that are as large as possible. The measurements of the investigated carbonates were conducted with the steady-state method by determining the size of the heat flow through the sample using a FOX 50 LaserComp apparatus. The tests were carried out on samples in the shape of slices of 5 cm diameter and 1.5 cm thickness. The mean temperature was 25 • C, and the difference between the heating and the cooling plates was 20 • C, with 5% accuracy. The measurements were conducted on dry samples (samples dried for 12 h in 105 • C) and saturated (water saturation) samples.
Grain density was determined using the helium pycnometry method with a Micrometrics (USA) AccuPyc 1330 apparatus. Bulk density was determined using the mercury displacement method. The porosity value was determined based on density measurements. The mineral content analysis was undertaken with the quantitative X-ray method based on the Rietveld technique [44] using a Panalytical (United Kingdom) X'Pert Pro X-ray diffractometer. The specific heat measurements were made using the differential scanning calorimetry method (TG-DSC) with a Netzsch (Germany) STA 449 F3 Jupiter thermal analyzer. Computed tomography investigations were conducted using a Geotek (United Kingdom) RXCT (Rotating CT system). The mercury injection capillary pressure (MICP) data were collected on a Micromeritics AutoPore IV 9520 device (USA) MICP analysis allows measurement of pore sizes from about 3.5 to 500 µm and provides a number of petrophysical parameters, e.g., pore size distribution, porosity, permeability, skeletal and apparent density, and specific surface area of a sample.
Organodetritical limestones from Kopiec-1G and Kopiec-4G ( Figure 5A,B) have beige, grey-beige, and grey hues with visible microbial and sponge structures also containing fragments of bivalvia, brachiopoda, and bryozoa. In several samples of stylolites, selective dolomitization and recrystallization processes were observed.
The samples from Opatkowice-OB1 and Trojanowice-2 ( Figure 5C,D) were represented by pelitic limestones with grain elements of beige and grey-beige hues, and sponges and bryozoa were visible in some places.
The mineral composition of limestones was minimally diversified (Table 3); the dominant mineral was calcite, and its content generally exceeded 99%. Quartz appeared in small amounts (below 1-2%) in most rocks. Only in samples 6 and 18 were high amounts of quartz-95.6% and 35.2%, respectively-found due to the presence of silificated sponges. In several rocks, dolomite admixtures (from 1.5% to 8%) were present, most likely connected with selective dolomitization processes.
The investigated rocks were characterized by highly diversified porosity of 1.63% to 12.14% ( Table 3). The lowest porosity value (below 2%) was found in sample 3 from Kopiec-1G and in samples 16 and 17 from Kopiec-4G, in which recrystallization processes were observed. Samples 11-13 from Kopiec-4G, and 21 and 22 from Trojanowice-2, possessed the highest (above 10%) porosity values. In these rocks, large, usually connected to sponges ( Figure 5B,D), vugs were macroscopically visible. Such vugs are described as primary porosity [46,47].   The results of investigations conducted with the MICP method verified the fractured and vuggy character of the pore space. Figure 6 shows both the cumulative intrusion curves obtained from MICP tests and the pore size distribution. The incremental intrusion curves ( Figure 6) represent the pore volume accessed through pore throats of a given size [48,49]. Assessment of the pore space geometry in the studied carbonates indicates the porous-fracture type of the reservoir space, which is mainly unimodal in nature. Pore size classification is based on the grading proposed by Hartmann and Beamont [50]. The dominant pore systems for the studied samples occur in the pore size diameter range of 0.02-5 µm. This indicates that the pore space corresponds mainly to micro-and nanopores ( Figure 6, Table 4). The relationship between the pore space formation, porosity, and thermal conductivity is most visible in the carbonate rock samples from Trojanowice-2. These are samples with dominant micropores (Figure 6), followed by high porosities obtained both by MICP and helium methods (Tables 3 and 4). The share proportion of particular pore types is the main factor conditioning the obtained different results of thermal conductivity of the analyzed samples. The differences in the porosity results obtained by helium and MICP methods are due to different measurement ranges of these methods [51], in addition to the measurement methodology. In the case of MICP, measurement is carried out using crushed samples, thus eliminating big pores and fractures, thus decreasing the porosity. The results of investigations conducted with the MICP method verified the fractured and vuggy character of the pore space. Figure 6 shows both the cumulative intrusion curves obtained from MICP tests and the pore size distribution. The incremental intrusion curves ( Figure 6) represent the pore volume accessed through pore throats of a given size [48,49]. Assessment of the pore space geometry in the studied carbonates indicates the porous-fracture type of the reservoir space, which is mainly unimodal in nature. Pore size classification is based on the grading proposed by Hartmann and Beamont [50]. The dominant pore systems for the studied samples occur in the pore size diameter range of 0.02-5 μm. This indicates that the pore space corresponds mainly to micro-and nanopores ( Figure 6, Table 4). The relationship between the pore space formation, porosity, and thermal conductivity is most visible in the carbonate rock samples from Trojanowice-2. These are samples with dominant micropores (Figure 6), followed by high porosities obtained both by MICP and helium methods (Tables 3 and 4). The share proportion of particular pore types is the main factor conditioning the obtained different results of thermal conductivity of the analyzed samples. The differences in the porosity results obtained by helium and MICP methods are due to different measurement ranges of these methods [51], in addition to the measurement methodology. In the case of MICP, measurement is carried out using crushed samples, thus eliminating big pores and fractures, thus decreasing the porosity.    Figure 7). Low TC values are generally reflected by high porosity (Table 4). A distinctive decrease in thermal conductivity values in line with increasing porosity (Table 3, Figures 8 and 9) is observed in regards to both the dry and saturated samples.
Heat capacity values (specific heat) range from 0.889 to 0.924 J/g K (Table 5). Table 4. Results of mercury injection capillary pressure (MICP) and thermal conductivity measurements for samples from the four borehole locations.

Sample
No.  Figure 7). Low TC values are generally reflected by high porosity (Table 4). A distinctive decrease in thermal conductivity values in line with increasing porosity (Table 3, Figures 8 and 9) is observed in regards to both the dry and saturated samples. Heat capacity values (specific heat) range from 0.889 to 0.924 J/g K (Table 5).

Mathematical Models
The mathematical models noted previously were used to compute the thermal conductivity of the investigated samples. Calculations were carried out for 22 carbonate rocks from the Kraków region. All analyses were performed on saturated samples.
The analyses were carried out using the following procedure: 1. Determination of the thermal conductivity of all samples using arithmetic, harmonic, and geometric means and comparing them with the laboratory measured conductivities ( Figure 10); 2. Removal of two outstanding samples of different mineral composition (high content of quartz) from the applied models ( Figure 11) and correlation analysis of obtained results ( Figure 12); 3. Applying the models to samples divided into two groups: of porosity higher and lower than 4%

Mathematical Models
The mathematical models noted previously were used to compute the thermal conductivity of the investigated samples. Calculations were carried out for 22 carbonate rocks from the Kraków region. All analyses were performed on saturated samples.
The analyses were carried out using the following procedure: 1. Determination of the thermal conductivity of all samples using arithmetic, harmonic, and geometric means and comparing them with the laboratory measured conductivities ( Figure 10); 2. Removal of two outstanding samples of different mineral composition (high content of quartz) from the applied models ( Figure 11) and correlation analysis of obtained results ( Figure 12); 3. Applying the models to samples divided into two groups: of porosity higher and lower than 4%

Mathematical Models
The mathematical models noted previously were used to compute the thermal conductivity of the investigated samples. Calculations were carried out for 22 carbonate rocks from the Kraków region. All analyses were performed on saturated samples.
The analyses were carried out using the following procedure: 1.
Determination of the thermal conductivity of all samples using arithmetic, harmonic, and geometric means and comparing them with the laboratory measured conductivities ( Figure 10); 2.
Removal of two outstanding samples of different mineral composition (high content of quartz) from the applied models ( Figure 11) and correlation analysis of obtained results ( Figure 16);

3.
Applying the models to samples divided into two groups: of porosity higher and lower than 4% (Figure 12), and correlation analysis of obtained results ( Figure 17); 4.
Determining the thermal conductivity using spherical inclusion models (Clausius-Mossotti) (Figures 13 and 18), taking into consideration the division of samples into two groups of different porosity ( Figure 21) and correlation analysis of obtained results ( Figure 24); 5.
Thermal conductivity determination of all samples using non-spherical inclusion models (Figures 14 and 19) taking into consideration the division of samples into two groups of different porosity ( Figure 22) and correlation analysis of obtained results ( Figure 25); 6.
Introducing a correction allowing for approximation of the modelled values ( Figure 15) to the laboratory measured values (Figures 20, 23 and 26-31).
Energies 2020, 13, x FOR PEER REVIEW 16 of 31 4. Determining the thermal conductivity using spherical inclusion models (Clausius-Mossotti) (Figures 15 and 16), taking into consideration the division of samples into two groups of different porosity ( Figure 17) and correlation analysis of obtained results ( Figure 18); 5. Thermal conductivity determination of all samples using non-spherical inclusion models (Figures 19 and 20) taking into consideration the division of samples into two groups of different porosity ( Figure 21) and correlation analysis of obtained results ( Figure 22); 6. Introducing a correction allowing for approximation of the modelled values ( Figure 23) to the laboratory measured values (Figures 24-31).
The thermal conductivity of the rock matrix used in the applied models was calculated based on the content and thermal conductivity of minerals ( Table 1). The content (percentage by weight) of minerals was obtained using XRD analysis (Table 3) and converted into volume percent.

Layer Models and Geometric Mean
The layer (λ_arytm and λ_harm) and geometric mean (λ_geom) models were compared with laboratory measurements for all investigated samples ( Figure 10). A distinct influence of mineral composition on matching calculated to measured thermal conductivity values was observed. Two rocks (samples 6 and 18, indicated with a red ellipse in Figure  10) are characterized by overestimated values that are distinct from others. These samples contain significantly higher amounts of quartz (35% and 95%) than other rocks (maximum 2.5%). Regarding different mineralogy, they were excluded from further consideration.
Comparison of λ_arytm and λ_harm ( Figure 11) with laboratory measurements verified the theoretic assumption: the highest values are indicated by λ_arytm (the model assuming the heat flowing parallel to the boundary between components) and the lowest by λ_harm (the model assuming the heat flowing perpendicular to the boundary between components). The values obtained based on the λ_geom model are similar to those of λ_harm ( Figure 11). Each of the three models overestimated results in regard to the laboratory measured values and displayed medium correlation factors (Figure 12a-c). Then, the set of data was divided into two subgroups: samples of porosity above 4% (Figure 13a) and samples of porosity below 4% (Figure 13b). Two samples characterized by calculated thermal conductivity values that differed most from laboratory measurements (samples 8 and 10, Figures 11 and 13a) were excluded from the group of higher porosity samples. This procedure allowed considerably better correlations for rocks of porosity higher than 4% to be obtained compared to the whole set of samples ( Figure 14). The best correlation coefficient (R 2 = 0.72) is characteristic for λ_geom (Figure 14b) and the calculated values are most similar to the measured values for λ_harm (Figure 14c). The values obtained by each of the three models for samples of porosity lower than 4% are significantly overestimated (Figure 13b) and do not correlate with the measured values. , λ_harm (c), and the laboratory measured values λ_lab for samples of porosity above 4% (see Figure  13a for samples numbers).

Clausius-Mossotti Spherical Inclusions Model
In the next step, the Clausius-Mossotti spherical inclusion models were tested [12]. Both of the applied models λ_sph_m (rock composed of a grain matrix and pores in the shape of spherical inclusions) and λ_sph_f (rock composed of spherical grains dispersed in a pore solution) show that thermal conductivity values are overestimated in relation to the measured values ( Figure 15). Similarly, as in the case of the layer and geometric mean models (Figures 12 and 14), the correlations carried out for the whole collection of samples ( Figure 16) improve after distinguishing a subgroup of samples of porosity above 4% (Figures 17 and 18). In that group, the λ_sph_f model yields values most similar to the measurements. Model λ_nsph_d gave slightly lower values that were more similar to the measured values. As in the case of other models, significantly better correlations were obtained for the distinguished subset of samples of porosity above 4% (Figures 21 and 22) than for all of the rocks (Figure 20).  above 4% than for the whole set of samples (R 2 = 0.55-0.59). The overestimated values obtained in the result of calculations are probably due to understated porosity values caused by the appearance of isolated pores, which are not reflected in the laboratory measurements. The influence of lithology on the calculated thermal conductivity values was also observed for siliciclastic rocks [21,22].
Comparison of all of the models calculated for the subgroup of rocks of porosity over 4% ( Figure  23) allowed the models that best fit the laboratory measurements to be identified.  Figure  13a for samples numbers).

Clausius-Mossotti Spherical Inclusions Model
In the next step, the Clausius-Mossotti spherical inclusion models were tested [12]. Both of the applied models λ_sph_m (rock composed of a grain matrix and pores in the shape of spherical inclusions) and λ_sph_f (rock composed of spherical grains dispersed in a pore solution) show that thermal conductivity values are overestimated in relation to the measured values ( Figure 15).  Model λ_nsph_d gave slightly lower values that were more similar to the measured values. As in the case of other models, significantly better correlations were obtained for the distinguished subset of samples of porosity above 4% (Figures 21 and 22) than for all of the rocks (Figure 20).     Figure 17a for samples numbers).

Non-Spherical Inclusion Models
For models assuming pores of ellipsoid shapes, extreme cases of elongated, needle shaped pores (λ_nsph_n), pores of spherical shape (λ_nsph_s), and disk-shaped pores (λ_nsph_d) were analyzed ( Figure 19). Similar results were obtained for λ_nsph_s (counterpart of spherical inclusion model λ_sph_m) and λ_nsph_n models ( Figure 19). Model λ_nsph_d gave slightly lower values that were more similar to the measured values. As in the case of other models, significantly better correlations were obtained for the distinguished subset of samples of porosity above 4% (Figures 21 and 22) than for all of the rocks (Figure 20).

Analysis of the Obtained Results
The comparison of the laboratory measured thermal conductivity values with those obtained by mathematical models clearly indicates the impact of both the mineralogy and petrophysical properties of the investigated rocks on the models tested. The influence of mineralogy is visible in the case of two samples (samples 6 and 18) with outstandingly high quartz content (Figure 10). These samples are characterized by the distinctively different from the other ones, overestimated thermal conductivity values obtained from layer and geometric mean models. As mentioned previously, these samples were excluded from further considerations. Another highly important factor is     The thermal conductivity of the rock matrix used in the applied models was calculated based on the content and thermal conductivity of minerals ( Table 1). The content (percentage by weight) of minerals was obtained using XRD analysis (Table 3) and converted into volume percent.

Summary and Conclusions
A method using mineralogical, petrophysical, and thermal data to determine the thermal potential of carbonate rocks from the Kraków region was introduced. Analyses of mineral composition, porosity, and thermal conductivity were conducted, and the thermal conductivity coefficient was assessed using mathematical models. The proceedings followed the investigation scheme presented in Figure 32.
The obtained model allows the thermal conductivity database to be extended in cases where it is not possible to measure this parameter due to a sparse amount of rock material.

Layer Models and Geometric Mean
The layer (λ _arytm and λ _harm ) and geometric mean (λ _geom ) models were compared with laboratory measurements for all investigated samples ( Figure 10).
A distinct influence of mineral composition on matching calculated to measured thermal conductivity values was observed. Two rocks (samples 6 and 18, indicated with a red ellipse in Figure 10) are characterized by overestimated values that are distinct from others. These samples contain significantly higher amounts of quartz (35% and 95%) than other rocks (maximum 2.5%). Regarding different mineralogy, they were excluded from further consideration.
Comparison of λ _arytm and λ _harm ( Figure 11) with laboratory measurements verified the theoretic assumption: the highest values are indicated by λ _arytm (the model assuming the heat flowing parallel to the boundary between components) and the lowest by λ _harm (the model assuming the heat flowing perpendicular to the boundary between components). The values obtained based on the λ _geom model are similar to those of λ _harm (Figure 11). Each of the three models overestimated results in regard to the laboratory measured values and displayed medium correlation factors (Figure 16a-c). Then, the set of data was divided into two subgroups: samples of porosity above 4% (Figure 12a) and samples of porosity below 4% (Figure 12b). Two samples characterized by calculated thermal conductivity values that differed most from laboratory measurements (samples 8 and 10, Figures 11 and 12a) were excluded from the group of higher porosity samples. This procedure allowed considerably better correlations for rocks of porosity higher than 4% to be obtained compared to the whole set of samples ( Figure 17). The best correlation coefficient (R 2 = 0.72) is characteristic for λ _geom (Figure 17b) and the calculated values are most similar to the measured values for λ _harm (Figure 17c). The values obtained by each of the three models for samples of porosity lower than 4% are significantly overestimated (Figure 12b) and do not correlate with the measured values.

Clausius-Mossotti Spherical Inclusions Model
In the next step, the Clausius-Mossotti spherical inclusion models were tested [12]. Both of the applied models λ _sph_m (rock composed of a grain matrix and pores in the shape of spherical inclusions) Energies 2020, 13, 5515 25 of 32 and λ _sph_f (rock composed of spherical grains dispersed in a pore solution) show that thermal conductivity values are overestimated in relation to the measured values ( Figure 13). Similarly, as in the case of the layer and geometric mean models (Figures 16 and 17), the correlations carried out for the whole collection of samples ( Figure 18) improve after distinguishing a subgroup of samples of porosity above 4% (Figures 21 and 24). In that group, the λ _sph_f model yields values most similar to the measurements.
Model λ _nsph_d gave slightly lower values that were more similar to the measured values. As in the case of other models, significantly better correlations were obtained for the distinguished subset of samples of porosity above 4% (Figures 22 and 25) than for all of the rocks ( Figure 19).

Analysis of the Obtained Results
The comparison of the laboratory measured thermal conductivity values with those obtained by mathematical models clearly indicates the impact of both the mineralogy and petrophysical properties of the investigated rocks on the models tested. The influence of mineralogy is visible in the case of two samples (samples 6 and 18) with outstandingly high quartz content ( Figure 10). These samples are characterized by the distinctively different from the other ones, overestimated thermal conductivity values obtained from layer and geometric mean models. As mentioned previously, these samples were excluded from further considerations. Another highly important factor is porosity. Two groups of samples of different pore space geometry-compact rocks of low porosity and vuggy samples of porosity above 4%-were distinguished in the set of the investigated rocks. Considerably better correlation coefficients (R 2 = 0.69-0.74) were obtained for rocks of porosities above 4% than for the whole set of samples (R 2 = 0.55-0.59). The overestimated values obtained in the result of calculations are probably due to understated porosity values caused by the appearance of isolated pores, which are not reflected in the laboratory measurements. The influence of lithology on the calculated thermal conductivity values was also observed for siliciclastic rocks [21,22].
Comparison of all of the models calculated for the subgroup of rocks of porosity over 4% ( Figure 15) allowed the models that best fit the laboratory measurements to be identified.
As shown in Figure 15, the best conformity of results was obtained for models λ _harm , λ _sph_f , and λ _nsph_d . Each of the three models assume a structure of the pore space in which the contacts between the framework grains are limited. As described above, model λ _sph_f is that of rock composed of spherical grains dispersed in the pore solution, model λ _nsph_d is characterized by the presence of extended fractures, whereas model λ _harm assumes the heat flow perpendicular to the layers consisting of respective components of the rock (including pore solutions). It appears that these theoretic assumptions can be applied to rocks where porosity is connected to unequally distributed and sometimes large vugs filled with gas or brines constituting "a barrier" against the heat flow.
A correction allowing for the approximation of values obtained by the models to the laboratory measured values was applied. Such a correction can be introduced when the correlation between the calculated and measured values is of sufficient quality. The obtained linear equations showing relationships between the calculated models and the laboratory measurements were used, and calculations were made for each of the models (Table 6).

Summary and Conclusions
A method using mineralogical, petrophysical, and thermal data to determine the thermal potential of carbonate rocks from the Kraków region was introduced. Analyses of mineral composition, porosity, and thermal conductivity were conducted, and the thermal conductivity coefficient was assessed using mathematical models. The proceedings followed the investigation scheme presented in Figure 32.
The obtained model allows the thermal conductivity database to be extended in cases where it is not possible to measure this parameter due to a sparse amount of rock material.
Due to computed tomography measurements, it was possible to depict the inner structure and to trace the distribution of fractures and vugs in the investigated limestones. These also allowed representative core fragments to be chosen for thermal conductivity and porosity measurements. It is particularly important in the case of rocks with high variability of the pore space, such as the investigated limestones. The tomography image clearly shows that porosity is related to unevenly distributed, sometimes large vugs and fractures which constitute a "barrier" against the heat flow. The best match of the measured and calculated results was obtained by the models that assumed such a structure of the pore space, in which the grain contacts of the framework are limited by either the advantage of pore solutions over the framework grains (λ _sph_f ) or the presence of extended voids (λ _harm and λ _nsph_d ). Due to computed tomography measurements, it was possible to depict the inner structure and to trace the distribution of fractures and vugs in the investigated limestones. These also allowed representative core fragments to be chosen for thermal conductivity and porosity measurements. It is particularly important in the case of rocks with high variability of the pore space, such as the investigated limestones. The tomography image clearly shows that porosity is related to unevenly distributed, sometimes large vugs and fractures which constitute a "barrier" against the heat flow. The best match of the measured and calculated results was obtained by the models that assumed such a structure of the pore space, in which the grain contacts of the framework are limited by either the advantage of pore solutions over the framework grains (λ_sph_f) or the presence of extended voids (λ_harm and λ_nsph_d).
Mineral composition investigations were essential in the construction of mathematical models, taking into consideration the volume of particular minerals with appropriate thermal conductivity coefficients. These measurements also made it possible to analyze the influence of mineralogy on the calculated models. As a result of the analysis, two samples (samples 6 and 18) of different mineral composition (high amounts of quartz) were excluded from the input data.
An analysis of porosity measurements proved the dominant influence of porosity on the thermal conductivity value and allowed the measurement data to be arranged to obtain the best fitting of the tested models. In most cases, the calculated thermal conductivity values were overestimated compared to the laboratory measured values. This was probably due to the presence of isolated pores that was not reflected in the laboratory measurements of porosity. This was confirmed by the considerably better conformity of results for all investigated models obtained for a selected subgroup of samples of porosity exceeding 4% (R 2 from 0.69 to 0.74) than for all samples (R 2 from 0.55 to 0.59).
Using the good correlation coefficients obtained for the group of rocks of porosities above 4%, a correction was made so that the calculated values approximated the measured thermal conductivities. The calculation of thermal conductivity based on both mineral composition and porosity ( Figure 33) is essential for completing the thermal database of carbonate rocks in the investigated region. As mentioned, the amount of material from a number of boreholes is too sparse Mineral composition investigations were essential in the construction of mathematical models, taking into consideration the volume of particular minerals with appropriate thermal conductivity coefficients. These measurements also made it possible to analyze the influence of mineralogy on the calculated models. As a result of the analysis, two samples (samples 6 and 18) of different mineral composition (high amounts of quartz) were excluded from the input data.
An analysis of porosity measurements proved the dominant influence of porosity on the thermal conductivity value and allowed the measurement data to be arranged to obtain the best fitting of the tested models. In most cases, the calculated thermal conductivity values were overestimated compared to the laboratory measured values. This was probably due to the presence of isolated pores that was not reflected in the laboratory measurements of porosity. This was confirmed by the considerably better conformity of results for all investigated models obtained for a selected subgroup of samples of porosity exceeding 4% (R 2 from 0.69 to 0.74) than for all samples (R 2 from 0.55 to 0.59).
Using the good correlation coefficients obtained for the group of rocks of porosities above 4%, a correction was made so that the calculated values approximated the measured thermal conductivities. The calculation of thermal conductivity based on both mineral composition and porosity ( Figure 33) is essential for completing the thermal database of carbonate rocks in the investigated region. As mentioned, the amount of material from a number of boreholes is too sparse to conduct laboratory measurements of thermal conductivity, and samples from outcrops are often weathered. Thermophysical rock properties are key parameters for the assessment of shallow geothermal resources, future design of geological work, effective heat/cold extraction, and sustainable resource management. From the study of the thermo-and petrophysical properties of the carbonate rocks in the Kraków area, the following conclusions can be drawn:


The studied Upper Jurassic rocks are a heterogenous formation of limestones, which can be differentiated in structure and fabric. These differences affect the thermo-and petrophysical properties of the rocks and show related trends;  Thermal conductivity was assessed by means of mathematical models based on a simplification of the rock's structure, allowing for calculation of the rock thermal conductivity based on the properties of its components;  The best correlations between calculated and measured thermal conductivity values (R 2 from 0.69 to 0.74) were obtained for the subgroup samples of porosity higher than 4%. The subgroup of samples of porosity lower than 4% did not show a satisfactory fit of data and was not taken into consideration in the subsequent analyses;  The models which yielded values most approximate to the laboratory measurements were λ_harm, λ_sph_f, and λ_nsph_d;  A correction based on the obtained linear equations showing relationships between the respective calculated models to the laboratory measurements was introduced. The corrected thermal conductivity values were considerably closer to the laboratory measurements;  The developed measurement workflow allows for the use of mineralogical and petrophysical analysis to identify the geothermal potential of carbonate rocks;  Calculated mathematical models of thermal conductivity provide information on the main factors controlling this property. This will enable selection of optimal rock parameters in the design phase of a shallow geothermal installation;  The proposed methodology appears to be justified for the case in which a reliable set of data is used to verify mathematical models of thermal conductivity of carbonate rocks. In particular, the discussed methods and models can be used to refine existing geothermal potential maps or calibration of 3D parametric models used to design prognostic calculations of technical installations for shallow geothermal energy use.
The preparatory research for this article, combined with the established methodology for estimating the thermal conductivity of carbonate rocks in the Kraków area, is the first stage of our work, which will most likely be continued by our team in the future. Thermophysical rock properties are key parameters for the assessment of shallow geothermal resources, future design of geological work, effective heat/cold extraction, and sustainable resource management. From the study of the thermo-and petrophysical properties of the carbonate rocks in the Kraków area, the following conclusions can be drawn:

•
The studied Upper Jurassic rocks are a heterogenous formation of limestones, which can be differentiated in structure and fabric. These differences affect the thermo-and petrophysical properties of the rocks and show related trends; • Thermal conductivity was assessed by means of mathematical models based on a simplification of the rock's structure, allowing for calculation of the rock thermal conductivity based on the properties of its components; • The best correlations between calculated and measured thermal conductivity values (R 2 from 0.69 to 0.74) were obtained for the subgroup samples of porosity higher than 4%. The subgroup of samples of porosity lower than 4% did not show a satisfactory fit of data and was not taken into consideration in the subsequent analyses; • The models which yielded values most approximate to the laboratory measurements were λ _harm , λ _sph_f , and λ _nsph_d ; • A correction based on the obtained linear equations showing relationships between the respective calculated models to the laboratory measurements was introduced. The corrected thermal conductivity values were considerably closer to the laboratory measurements; • The developed measurement workflow allows for the use of mineralogical and petrophysical analysis to identify the geothermal potential of carbonate rocks; • Calculated mathematical models of thermal conductivity provide information on the main factors controlling this property. This will enable selection of optimal rock parameters in the design phase of a shallow geothermal installation; • The proposed methodology appears to be justified for the case in which a reliable set of data is used to verify mathematical models of thermal conductivity of carbonate rocks. In particular, the discussed methods and models can be used to refine existing geothermal potential maps or calibration of 3D parametric models used to design prognostic calculations of technical installations for shallow geothermal energy use.
The preparatory research for this article, combined with the established methodology for estimating the thermal conductivity of carbonate rocks in the Kraków area, is the first stage of our work, which will most likely be continued by our team in the future.

Acknowledgments:
We would like to thank Andrzej Urbaniec for his help in the macroscopic description of the investigated limestone samples.

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

Terminology and Abbreviations
TC thermal conductivity CT computed tomography BHE borehole heat exchanger TRT thermal response test GSHP ground source heat pump λ m thermal conductivity of the grain framework λ f thermal conductivity of the pore media φ porosity V volume fraction of the inclusions λ thermal conductivity of the whole rock λ 1 thermal conductivity of the inclusion material λ 2 thermal conductivity of the host material λ _lab laboratory measured thermal conductivity of saturated sample λ _aryt thermal conductivity calculated by applying arithmetic mean (Equation (1)) λ _harm thermal conductivity calculated by applying harmonic mean (Equation (2)) λ _geom thermal conductivity calculated by applying geometric mean (Equation (6)) λ _sph_m thermal conductivity calculated by applying spherical inclusion models, for rock consisting of grain framework and pore solutions appearing as spherical inclusions (Equation (4)) λ _sph_f thermal conductivity calculated by applying spherical inclusion models, for rock consisting of spherical grains suspended in a fluid (Equation (5)) λ _sph_av thermal conductivity calculated by applying spherical inclusion models, average value λ _nsph_s thermal conductivity calculated by applying non-spherical inclusion models, case α→1, pores of the sphere shape (Equation (7)) λ _nsph_n thermal conductivity calculated by applying non-spherical inclusion models, case α→∞, pores of the needle shape (Equation (7)) λ _nsph_d thermal conductivity calculated by applying non-spherical inclusion models, case α→0, pores of the disk shape (Equation (7)) λ _aryt-cor corrected thermal conductivity value, calculated by applying arithmetic mean λ _harm-cor corrected thermal conductivity value, calculated by applying harmonic mean λ _geom-cor corrected thermal conductivity value, calculated by applying geometric mean λ _sph_m-cor corrected thermal conductivity value, calculated by applying spherical inclusion models, for rock consisting of grain framework and pore solutions appearing as spherical inclusions λ _sph_f-cor corrected thermal conductivity value calculated by applying spherical inclusion models, for rock consisting of spherical grains suspended in a fluid λ _nsph_n-cor corrected thermal conductivity value calculated by applying non-spherical inclusion models, case α→∞, pores of the needle shape λ _nsph_d-cor corrected thermal conductivity value calculated by applying non-spherical inclusion models, case α→0, pores of the disk shape