Fast calculation of the technical shallow geothermal energy potential of large areas with a steady-state solution of the finite line source

,


Introduction
The heating and cooling sector is one of the most energy intensive sectors in the European Union (30 % of all energy consumed) and accounts for the majority of the final energy use per household (79 % in EU households, Fleiter et al., 2017).In 2018, 75 % of the energy used in the space heating and cooling sector was generated from fossil fuels.In order to meet climate targets set out in the Paris Agreement and the European Green Deal (European Parliament, 2020;UNFCCC, 2015) and independence from fossil fuel imports as set out in the REPower EU Action (European Commission, 2022), this sector is in a high need for decarbonisation and a reduction of energy consumption.As part of the decarbonisation drive, municipalities within the EU are required to develop heating and cooling plans in the next decade (European Parliament, 2023), with tighter regulations already in place for some parts of the EU such as the German state of Baden-Wuerttemberg, where larger municipalities need to have a heating and cooling plan in place by the end of 2023 (Baden-Württemberg, 2023).Commonly discussed options for the heating and cooling sector in renewable energy systems include the decentralisation of energy systems (Orehounig et al., 2015), district heating (Inayat and Raza, 2019;Lund et al., 2010), and the widespread use of heat pumps (Gaur et al., 2021;Lund, 2007).From an efficiency point of view, ground source heat pumps (GHSPs) are preferable to air-water heat pumps as the ground has less temperature fluctuations than ambient air and can provide a more stable temperature throughout the year resulting in a higher source temperature in winter.Thus, an extensive utilization of shallow geothermal energy could help to strongly reduce carbon emissions from the heating and cooling sector (Lund and Boyd, 2016).
The most commonly installed shallow geothermal energy systems use vertical borehole heat exchangers (BHEs) through which a heat carrier fluid is circulated in a closed tubing system.If two BHEs are placed close to each other (i.e.within a radius less or equal to the depth of a BHE), thermal interference between the boreholes occurs, degrading the heat extraction potential of both BHEs (Meng et al., 2019;Vienken et al., 2015).High BHE densities may also result in a long-term decrease of ground temperatures and lead to a decreasing GSHP efficiency over time (Patton et al., 2020).The thermal interference between BHEs is thus a crucial element of any potential analysis of shallow geothermal energy.
Commonly, the feasibility and potential of GSHPs are explored on a local or city district level, and various approaches have been developed in the past (Bayer et al., 2019;Casasso and Sethi, 2016;Heim et al., 2022;Noorollahi et al., 2017).While such studies can support urban planners and policy makers in individual cases, regional studies are needed to truly enable the utilisation of the shallow geothermal energy potential on the scale needed for a successful energy transition.The majority of studies that assess the regional-scale potential of shallow geothermal energy estimate the theoretical potential (Bertermann et al., 2015;Majorowicz et al., 2009), which is the physically available energy in a given ground volume (Bayer et al., 2019).However, the technical potential, which is the technically extractable amount of heat with consideration of the built environment and the thermal interference, is more useful for evaluating the potential of the implementation of GSHPs on a large scale.Some studies quantify the regional technical potential but do not take the thermal interference of boreholes into account (Casasso and Sethi, 2016;Galgaro et al., 2015).This is partly because the calculation of the thermal interference for large borehole arrays is computationally intensive: studies which take the thermal interference into account either rely on the distribution of boreholes along a set grid and precomputed interferences (Walch et al., 2021), or estimate the interference based on BHE density (Miocic and Krecher, 2022).
The effective temperature variation at the borehole walls in a BHE field due to a constant heat extraction rate from a BHE field can be calculated using thermal response factors, or g-functions (Eskilson, 1987).These g-functions can be used to superimpose heating (and cooling) loads to predict variations of borehole wall temperatures throughout time (Cimmino and Cook, 2022).Initially developed via numerical finite difference models (Eskilon, 1987), analytical solutions of the finite line source have been proposed to evaluate g-functions (Cimmino, 2018;Lamarche and Beauchamp, 2007).Following a series of simplifications and algorithmic improvements to speed-up the evaluation of g-functions, the continuously developed open-source pygfunction package (Cimmino, 2018) is often seen as state of the art when it comes to the evaluation of g-functions of arbitrarily positioned boreholes.However, for the g-function evaluation of a large number (> 100.000) of BHE fields with a large number of BHEs with varying lengths/depths for regional studies the computational time is a limitation.
Here, we argue that the temperature variation at the borehole walls caused by thermal interference by neighbouring boreholes within a borehole field in a steady state, which commonly only commences after years of operation, are of importance when a regional potential is analysed, and that g-functions prior to the steady state can be ignored for considering this effect.An analytical solution of the steady-state finite line source solution is proposed which enables the rapid evaluation of steady-state g-functions of large borehole fields with BHEs of varying depths.This methodology is then applied to calculate a regional technical shallow geothermal potential for the German state of Baden-Württemberg, a study area of 35.000 km 2 and more than 8.5 million BHEs.

Materials and methods
The regional technical shallow geothermal potential can be calculated using the following steps: (1) Identify suitable BHE locations in the study area; (2) Assess the interference between the placed BHEs; (3) Calculate the maximum technical possible heat extraction of all BHEs taking the thermal interference between individual BHEs into account and determine the technical shallow geothermal potential of each parcel in the study area.

BHE placement & length 2.1.1. BHE placement
Suitable parcels for placing BHEs were identified by filtering the cadastral land register data (ALKIS) of the federal state of Baden-Württemberg, Germany, for all parcels with residential building using QGIS (v.3.24).To ensure a minimum distance between BHEs and buildings, a buffer of 2 m was placed around all buildings located on these parcels.Subtracting the building and buffer area from the parcels gave the area where theoretically BHEs could be placed, not considering obstacles such as large trees or garden sheds not registered in the cadastral land register.Within this available area, BHEs were randomly placed using a minimum distance between BHEs of 10 m (globally).This approach was selected as it allows for a higher BHE density in build-up areas.To evaluate the technical shallow geothermal potential two cases are studied: A low potential, with one BHE per parcel, and a high potential, where up to 20 BHEs were placed per parcel.Note that these potentials can be seen as minimal potentials as cooling is not taken into consideration.
The placement of BHEs in the study area is partly prohibited due to groundwater protection areas or mineral water production areas.Therefore, these areas were excluded.One exception are groundwater protection areas of Zone IIIB which are considered, but BHEs within these areas are only allowed to operate with pure water as heat carrier as opposed to a glycol-water mixture in all other areas.In certain regions with a complex geology, the placement of BHEs has to be assessed on a case-by-case basis by the state office for geology, resources, and mining (LGBR).These regions were therefore also excluded from the study.

BHE length
Additionally, the geological setting of the study area, particularities in tectonically strongly affected regions and in areas where swellable rocks (e.g., anhydrite) are present within the shallow subsurface leads to further restrictions of BHE placement, particularly the maximum allowed drilling depth.The State Geological Service of the federal state of Baden-Württemberg (LGRB) provided maximum allowed depth data which takes these geological restrictions into account.The depth data is based on a regional 3D geological model, and is available with a limited resolution online (https://isong.lgrb-bw.de/).Based on the geological setting, placed BHEs were assigned the maximal allowed drillable depth or 100 m, whichever restriction applies first.The 100 m was chosen due to stricter state regulations at deeper depths.BHEs with a depth of less than 10 m were disregarded considering economic aspects and because such very shallow systems are similar to shallow horizontal heat collectors which need a different modelling approach.

Defining BHE fields
As neighbouring BHEs interfere with each other, borehole fields were

Distance between borehole heat exchangers (m). H Height of the boreholes (m). r,z
Cylindrical coordinates (m).z ′ Coordinate, where the heat source is located (m).

Δθ
Change of the dimensionless temperature response.

Subscripts rec
Referring to the receiver borehole heat exchangers.trans Referring to the transmitter borehole heat exchangers.

GSHP
Ground source heat pump.
J.M. Miocic et al. defined to calculate the impact of BHE interference within these fields (Fig. 1).The borehole fields were defined based on the individual parcels by placing a buffer of the maximal borehole depth within the parcel around the parcel, and all BHEs within this buffer (and within the parcel) were defined as a borehole field.The buffering was carried out in QGIS, while the intersection between the BHE locations (points) and the buffer (polygon) was implemented using the R package sf (v.1.0.10,(Pebesma et al., 2021).For each of the resulting borehole fields the thermal interference between boreholes was then calculated (see below).

BHE interference 2.2.1. Analytical solution of the steady-state finite line source
BHEs are commonly dimensioned regarding their heat extraction using the finite line source model (Eq.( 1)).This model is based on the instantaneous point source, which is then further integrated over the length H of the borehole and mirrored at the ground surface using the method of images to bound the infinite solid (Carslaw and Jaeger, 1986).Eq. ( 1) can be used to calculate the dimensionless temperature response θ of a BHE with length H and radius r at a considered time t: Where α is the ground thermal diffusivity, z' is the location of the source and z is the location at which the equation is evaluated.
There are two common approaches to evaluate this model concerning its height.Either it is evaluated at the middle of the borehole depth (z = H/2) or averaged by integrating over the entire length H of the BHE and then dividing by the length H (Eskilson, 1987).In the direction r perpendicular to the line source, the evaluation is usually carried out at the borehole wall, i.e. at a distance r = r b from the line source.
As thermal interference occurs for neighbouring BHEs, this effect has to be considered while determining the geothermal potential.Therefore, a distinction must be made between the transmitter BHE, from which the thermal interference originates, and a receiver BHE, which is affected by the interference.If the interference of BHE 1 on BHE 2 is taken into account, BHE 1 is called the transmitter BHE and BHE 2 the receiver BHE.This designation changes if the interference effect of BHE 2 on BHE 1 is to be determined.By using the line source methodology, the line source always represents the transmitter BHE.If the line source is not evaluated at the borehole wall but at a distance r = B, where the receiver BHE is located, the thermal influence of the transmitter BHE (Δθ) with length H trans on the receiver BHE with another length H rec can be determined (Fig. 2).
This can be done in a dimensionless way by defining the mean dimensionless temperature change caused by the transmitter BHE on the receiver BHE, which subsequently can be superposed by the mean dimensionless temperature response at the receiver BHE itself.This is mathematically done by integrating Eq. ( 1) once more over the height H rec and then dividing by H rec to get the mean value, as shown in Eq. (2).
For steady-state conditions, which are of primary interest when it comes to calculating the interference of BHEs at a certain distance for dimensioning of a BHE field, time becomes infinite and Eq. ( 2) simplifies to: The double integral of Eq. ( 3) can be analytically solved to remove the need for numerical integration: Fig. 1.Definition of borehole fields based on parcel buffering.The radius (r) of the buffer around a parcel equals the length (H) of the deepest borehole in the parcel considered.For the resulting borehole field, the steady-state g-function is calculated using the methodology described below.This is carried out for every parcel in the study area.The mean dimensionless temperature change calculated with Eq. ( 4) is only due to the interference of one transmitter BHE on one receiver BHE.If various boreholes cause an interference, the temperature changes of all surrounding transmitter BHEs have to be calculated and added up.Subsequently, all dimensionless temperature responses due to the neighbouring BHEs can be superposed by the dimensionless temperature change at the considered (receiver) borehole wall due to its own heat extraction or injection.This allows the determination of the overall temperature change at the borehole wall.

Superposition technique
The solution of the interference of one BHE on other BHEs discussed above has been implemented in a Python code (https://github.com/johomio/GSHP_potential) which allows the calculation of the mutual influence of large number of BHEs.The code uses the following steps to derive field-wide g-functions, with the x,y-coordinates as well as the depth of the individual BHEs as input parameters.

Distance calculation.
The distance between two BHEs is calculated using the Pythagoras theorem and all distances are fed into a n x n matrix (Table 1).The main diagonal is always 0, because the distance of a BHE to itself is 0. Further, the matrix is symmetrical to the main diagonal, since B ij = B ji . (5)

Assignment of transmitter and receiver height.
There is an associated height for each BHE.Since each BHE acts both as a transmitter and receiver, the individual heights have to be assigned accordingly.
A list with the respective heights of all BHEs is used to derive two n x n matrices of the BHE height, one for the transmitters and one for the receivers (Fig. 3).For the receiver matrix, the rows are expanded until the size n x n is reached (Fig. 3c), while for the transmitter matrix the initial list of BHE heights is transposed (Fig. 3b) and then expanded (Fig. 3d).Note that both resulting matrices have equal main diagonals.
The advantage of this representation is that all input parameters are in the shape of an n x n matrix with the same size.This allows the values that are in the same cell in all matrices to be used in the algorithm one after the other, which minimizes the computing time.

Calculation of the dimensionless temperature response including thermal interference
. By utilising the matrices described above, Eq. ( 4) can be used to calculate the interference between every pair of boreholes within the field.The distance matrix serves as input parameter for r, the transmitter matrix for H trans and the receiver matrix for H rec .This results in a matrix with every single interference between two boreholes within the field.Next, the diagonal of this matrix is replaced with the steadystate value of the g-function of the individual BHE under consideration.Specifically, a steady-state g-value of 6.6 is used for a single BHE with a predetermined radius and length, which is adjusted to the ratio of r b /H using Eskilson's radius correction (Eskilson, 1987).In order to get the overall dimensionless temperature response of a specific borehole, i. e., including the interference of all surrounding boreholes, the resulting matrix is summed up row-wise.Finally, the mean dimensionless temperature response for the whole BHE field is caluclated by summing up all values within the matrix and dividing by the number of BHEs.The results of the presented method to calculate the steady-state dimensionless temperature response including thermal interference is compared with the established method implemented in pygfunctions (Cimmino and Cook, 2022).The comparison yields almost identical results, as can be seen from percentage deviation of the dimensionless temperature responses in Fig. 4.

Computational efficiency
For testing the computational efficiency, steady-state g-functions for BHE fields with 1000, 3000, 5000, and 10,000 randomly distributed BHEs were calculated with both our python code and the pygfunction tool (Cimmino and Cook, 2022).As seen in Fig. 4, the calculation times increase exponentially with the number of BHEs within the borehole

Table 1
Distance matrix between all BHEs within a borehole field.field.However, our proposed calculation procedure for the g-function is significantly faster than the pygfunction tool.

Determination of the technical geothermal potential
The geothermal potential of the BHE fields is calculated using the well accepted approach that the heat extraction rate of each BHE can be divided into three main components: a stationary component, which includes the impact of heat extraction over long periods of time as well as the interaction between BHEs, an annual periodic component, which takes heating seasonality into account, and a peak load component, which considers peak loads over small periods of time, i.e. operation intervals of the heat pump (HBC, 2022;Koenigsdorff, 2011).For a detailed description of the calculation methodology of the technical geothermal potential the reader is referred to Miocic & Krecher (2022).The methodology needs several input parameters which include the annual average surface temperature, terrestrial heat flow, thermal conductivity of the subsurface, borehole related parameters (radius, thermal resistivity, g-function), operational parameters (annual operation time, peak load time), as well as regulatory restrictions including the maximal temperature spread between inflow and outflow of the heat pump.Table 2 lists these input parameters.The geothermal potential calculations have been implemented in a R code (https://github.com/johomio/GSHP_potential).The output is given per parcel and includes the maximum heat extraction rate per meter borehole (W/m), geothermal power (kW), and energy yield (kWh) of the GSHP system within the parcel considered.The regional technical geothermal potential can then be calculated by summing up output from the parcels of interest.

BHE placement
Within the study area, 1.509 million parcels suitable for the placement of BHEs were identified.The resulting amount of BHEs that can be placed within these parcels while accommodating spatial restrictions ~ 1.5 million in case of 1 BHE per parcel, and 8.60 million when up to 20 BHEs are placed per parcel.In the latter case, the median number of BHEs per parcel is 4.45, with the majority of parcels having less than 10 BHEs (Fig. 5).The parameters used to randomly place BHEs do have a small impact on the amount of BHEs that can be placed, however the random placement does allow for a better use of space then using a gridlike placing approach (Fig. 7).Although the placement of the BHEs takes the built environment into account, the total number of BHEs is likely overestimated as placement restrictions on a per parcel level additional to the building such as garden sheds or large trees are not accounted for.

BHE thermal interference
The interference between neighbouring BHEs was calculated using the analytical steady-state finite line source solution (see 2.2.).Two different potential cases based on the number of BHEs placed per parcel, called low potential and high potential, were calculated.The impact of the thermal interference on the calculated heating power is illustrated in Fig. 7.For the low potential (one BHE per parcel) the mean steady-state value of the g-function of the borehole fields is 11.4,while for the high potential (up to 20 BHEs per parcel) it is 28.8 on average.The low potential steady-state values of the g-functions are close to normally distributed while the high potential values show a bimodal distribution (Fig. 6).This is likely a reflection of parcel size and building structure,  with higher g-functions values occurring in areas with a high parcel density, while low g-function values occur in individual rural parcels (Fig. 7).Similarly, it can also represent the variable geology of the study area, as there are many areas with restrictions of the maximal BHE depth.Given the minimum separation distance of 10 m between BHEs, shallow BHEs interfere less strongly than deeper BHEs.The way the borehole fields are defined strongly influences the interference between BHEs calculated.Here, the borehole fields are defined by a buffer around the individual parcels, with the buffer radius equal to the maximum borehole depth within the parcel.The steady-state values of the gfunctions from the interference calculation represent the field-wide average interference between individual BHEs.However, the spatial distribution of BHEs within the borehole field, as well as BHEs located outside the borehole field, will influence the true interference field.Thus, the calculated average g-function value for the borehole field may either be higher or lower than the true value for the BHEs within the parcel for which the geothermal potential is calculated.A potential solution could be to calculate the interference for the borehole field, but   (f) Technical geothermal potential with no thermal interference (g-function of each BHE = 6.6).This case results in much higher calculated heating powers for all parcels and highlights the need for taking the thermal interference into account when calculating a technical geothermal potential.
only take the average of the g-function value of the BHEs located within the parcel.However, as the placement of the BHEs is strongly influenced by the local planning, this is not necessary for a regional potential assessment, but should be taken into account when calculating the potential of individual borehole fields.

Technical geothermal potential
The overall technical shallow geothermal potential of the whole study area ranges from 9.3 TWh/a for the low potential study with one BHE per parcel to more than 34 TWh/a when the higher amount of BHEs is distributed.For the lower geothermal potential, the heating power supplied by BHEs per parcel are close to normally distributed and average at 3.27 kW (Fig. 8a).The corresponding specific heat extraction rates average at 32.6 W/m (Fig 8c).For the high potential, the heating power per parcel has an average of about 9 kW and shows a distribution with a positive skewness (Fig. 8b).The heat extraction rates for the high potential have an average of 20.4 W/m (Fig. 8d), which is significantly lower than for the lower potential and results from the strong thermal interference between the more closely placed neighbouring boreholes.The potential varies within the study area, predominantly due to differences in soil thermal conductivities and restrictions in BHE depth (Fig. 9).
The wide range of the heating power occurring in the high potential case is the result of the highly variable number of BHEs that can be placed per parcel.Large single parcels, as commonly observed with farms in rural areas, can have very high power potential, while small parcels within densely built-up areas will have only a small number of BHEs which are subject to strong thermal interference.
While the number of BHEs was prioritised during BHE placement, the results of the maximum potential case show that for optimising the shallow geothermal potential it may be advantageous to take thermal interference into account: If BHEs are optimally placed (in x, y direction) with respect to the thermal interference while considering the optimal borehole depth, a lower number of BHEs may be able to extract a similar amount of heat as the high potential calculated here.
The calculated technical geothermal potential per parcel has a median ranging from 11 to 28 kWh/m 2 /a for the low and high geothermal potential respectively (Fig. 8e&f).This is in the same range as a previously calculated technical geothermal potential for the same study area which had a mean of 25.7 kWh/m 2 /a (Miocic and Krecher, 2022).For a neighbouring region in Northern Switzerland Walch et al. (2021) have estimated an average technical geothermal potential of 16.4 kWh/m 2 /a.Note that, while the low potential has a lower technical potential, the difference is not as high as the difference in BHEs per parcel would suggest.This highlights the optimisation potential and that for placing BHEs the thermal interference has to be taken into account.

Practical implications
For urban and rural planning, policy making, and the development of regulations, regional scale estimations of renewable energy potential, including the technical shallow geothermal potential, are required (Gormally et al., 2012;Walch et al., 2021).Heating and cooling plans will be required for all larger (> 45.000 inhabitants) EU municipalities in the future, and for the German state of Baden-Wuerttemberg by the end of 2023.For German municipalities the most recent proposal for a federal law on heating planning and decarbonisation of the heating sector requires heating and cooling plans by 2028.The presented methodology and results are available for all municipalities within the study area of Baden-Württemberg (Germany) (through the state energy agency KEA-BW) and are used to develop local municipal heating plans which will form the base for climate-neutral heating supply by 2050.The proposed methodology, combined with building heat consumption data and heat demand assessment, will thus support the implementation of the energy transition.
While the large-scale regional technical shallow geothermal potential presented here can be used for potential evaluation on a regional scale, it is associated with uncertainties arising from the choice of input parameters such as the COP of the heat pump or the thermal conductivity of the soil and more detailed local models need to be applied before BHEs are installed.These models should take the exact BHE locations, heat pump specifications and operational parameters, the geological situation including the exact thermal conductivity, the thermal interference, the heat transfer from ground water flow, and the urban heat island effect into account.However, the methodology to calculate the steady-state thermal interference of large BHE fields with varying distances between BHEs and variable BHE depths presented here will allow the rapid calculation of regional scale shallow geothermal potential elsewhere.It can also be used for the local specific planning of individual areas or parcels with arbitrarily arranged BHE fields.
Technical potential assessments, such as the one outlined in this study, typically concentrate on a single energy resource, aiming to determine its maximum achievable potential.However, this approach exhibits several limitations.Employing a diverse array of resources, such as a combination of district heating, shallow geothermal systems, and other renewable heat sources, is likely to yield significantly higher efficiency, e.g., when thermal regeneration of BHE fields by heat input is taken into account.In the case of the shallow geothermal potential presented in this study, the thermal interference between adjacent BHEs results in a calculated potential that is lower for most parcels than what would be observed in reality, as not all parcels would utilize BHEs.Nevertheless, technical potential assessments can still be valuable for identifying areas best suited for the evaluated technology.

Summary and conclusion
Ground source heat pumps could strongly reduce carbon emissions from the space heating sector.For urban and rural planning, understanding the technical potential of these systems is crucial.Regionalscale studies, which are needed for a high-level planning, need to take the thermal interference between individual BHEs into account.Here, a new methodology for the rapid calculation of the thermal interference of large borehole fields with varying borehole depths, based on the analytical solution of the steady-state finite line source is presented.

Fig. 2 .
Fig. 2. Visualisation of the influence of a transmitter BHE on the receiver BHE.See text for details.

4 )Fig. 3 .
Fig. 3. Visualisation of the height assignment procedure on the basis of the tables: a) list of the heights for each BHE b) transposed heights-matrix c) heights-matrix of the receiver-BHEs d) heights-matrix of the transmitter-BHEs.

Fig. 4 .
Fig. 4. Calculation times (on the left-hand axis) and the deviation of the values of the steady-state solution (on the right-hand axis) of a BHE field g-function for different borehole field sizes, for both the newly proposed methodology and the pygfunction toolbox (Cimmino and Cook, 2022).

Fig. 5 .
Fig. 5. Histogram of number of BHEs per parcel when BHE placement is maximized.

Fig. 6 .
Fig.6.Histograms and kernel density plots of the steady-state g-function per parcel for the low potential (left) and high potential (right).Due to the higher number of BHEs within the borehole fields interference is much higher in case of the high potential.

Fig. 7 .
Fig. 7.Maps illustrating the steady-state g-functions of individual parcels for the high potential (up to 20 BHEs per parcel) and the impact of BHE placement as well as the impact of thermal interference for a small subsection of the study area: (a) BHE placement along a grid, (b) random placement A which represents the algorithm that was used for the whole study area, (b) random placement B which uses a different seed within the placement algorithm, which results in a different placement, (d) Half-violin-point plot showing the distribution of steady-state g-functions for the different placements.Note that for each placement a different number of BHEs can be placed, with the random placements achieving much higher BHE densities.Also note the correlation between the built environment and gfunctions, with strongly build up areas with large residential parcels having higher steady-state function values than individual rural parcels.(e) Technical geothermal potential taking thermal interference between BHEs into account (based on placement and thermal interference illustrated in (b).(f) Technical geothermal potential with no thermal interference (g-function of each BHE = 6.6).This case results in much higher calculated heating powers for all parcels and highlights the need for taking the thermal interference into account when calculating a technical geothermal potential.

Fig. 8 .
Fig. 8. Distributions of heating power, heat extraction, and geothermal potential per parcel for the low potential (left), and the high potential (right) models.

Fig. 9 .
Fig. 9. Map of the study area showing the high potential case upscaled to a 100×100 m raster.

Table 2
Input parameters for the determination of the geothermal potential of each parcel.