Methodology for estimating the ground heat absorption rate of Ground Heat Exchangers

Abstract In Ground Source Heat Pump systems, the heat exchange rate is an important factor with regard to the initial cost of the system. When the Ground Heat Exchanger (GHE) is installed in a lithology with high thermal properties in the presence of groundwater, the heat exchange rates are larger than in the cases with poor thermal response of the ground and no groundwater. This research, hence, focuses on a methodology of measuring and analyzing the thermal properties of the lithologies encountered in an area, which can be used for the prediction of heat injection rates of a GHE, depending on its characteristics, the installation area ground properties and groundwater flow. A tool was created with the use of FlexPDE software, and a study case was chosen in order to validate the results. Twenty-two, 100 m in depth, boreholes located in Lefkosia (Cyprus) were tested through simulation for their geothermal performance over time. Subsequently the estimated heat load for the boreholes, after 24 h of operation in cooling mode, was used with the help of Geographic Information System software for the compilation of a heat load per meter depth map that can be transferred to the ground by a GHE. A review of similar studies and Geographical Information System applications referring to other countries is also presented and their results are compared to the results of this study. The step by step procedure presented in this paper can be used by engineers handling geothermal projects as a useful guide for sizing GHEs and calculating the heat injection rates of any area.


41
Shallow geothermal energy is the heat stored beneath the surface of the earth and keeps the upper 42 parts of the crust at a steady temperature, which is approximately equal to the mean annual 43 atmospheric temperature of the year (Busby, 2015;Rybach and Sanner, 2000;Heath, 1983;Gustav, 44 1840). Technologies developed for the exploitation of this free energy type are getting more and 45 more important as heating and cooling nowadays represent around 50% of EU's total energy 46 consumption. Buildings consume more than 67% of thermal energy in Europe and has become a 47 large untapped potential for renewable and energy efficiency (ReGeoCities Project, 2015).

48
Shallow geothermal energy constitutes a renewable energy source with high energy savings for 49 heating and cooling in residential and commercial buildings. Ground Heat Exchanger (GHE) 50 systems with the help of Heat Pumps can use the ground as a heat source or sink and achieve up to 51 70% energy savings compared to traditional heating/cooling systems (ReGeoCities Project, 2015). 52 In addition to their high efficiency, GHE systems are attractive due to their environmental 53 friendliness.

54
In a Ground Source Heat Pump (GSHP) system during the heating cycle, a fluid circulates through 55 the loop of the GHE absorbing heat from the ground, with the heat pump unit delivering heat into 56 the building. For cooling, the process is simply reversed, with the GHE rejecting heat to the ground.

57
To design the GSHP system a thermal analysis is necessary. Its main objective is to predict the 58 thermal response of the system and size the heat pump to be used. The efficiency of the geothermal 59 project is greatly affected by the correct sizing of tubes (Christodoulides et

65
The current paper discusses cases where it is desired to use GSHP systems to provide heating or 66 cooling. Twenty-two study cases were set up to simulate the geothermal performance of boreholes 67 over time, and to determine the heat injection load for these boreholes. Measurement and analysis 68 of the thermal properties of the ground, the water level and velocity followed for obtaining actual 69 design data. A validation tool was created with the use of FlexPDE software in order to achieve 70 optimization of the system. The actual results and a series of maps drawn using the data can be 71 used by engineers as a useful guide for sizing vertical GHE in their study (design) area, the Greater 72 Lefkosia area. Similar maps can be compiled by engineers handling geothermal projects using the 73 geological data of their area of interest and following the same methodology explained in sections 74 3 and 4.

75
Similar studies, concerning the simulation of the performance of vertical GHE over time and the 76 determination of the thermal load and amount of energy extraction that can be sustained are also 77 available for Germany, Japan, England, Wales, France and Ireland and are presented in the 78 following studies. the spatial correlation of these design and performance factors with the geology, legislation, 82 population issues, and climatic conditions using real data was evaluated. This theoretical study is 83 based on large-scale systems, i.e. system capacity ≥ 175 kW and the determined specific heat 84 extractions used a range between 33 and 63 W m -1 . Although spatial variations of the specific heat 85 extractions are locally observed, the majority of the studied zones have average specific heat 86 extractions between 48 and 50 W m -1 showing no or only a minor correlation between the chosen 87 specific heat extraction and the site-specific geology.

88
For Japan two studies, one for Chikushi Plain and a second one for the capital city, Tokyo, are 89 presented. Nam and Ooka (2011) estimated the effect of the ground's thermal properties on a GHE 90 system performance using a numerical simulation tool in the 23 wards of Tokyo. A method to 91 develop the energy potential for GSHP systems was suggested and its application was performed 92 using geographical information system (GIS) data. Their results showed that the heat exchange 93 rates in cooling season have a very small variation, from 40 to 42 W m -1 . The methodology used in our research is explained in the sections that follow. In Section 2 is 148 presented the geological structure of the Greater Lefkosia Area that was used as a study case. In

173
From the geological point of view, the study area is extensively analyzed by the "Bedrock Geologic 174 Map of the Greater Lefkosia Area" (Harrison et al., 2008). "Geologic" or "geological" map is a 175 general term used by geologists to describe the rock types, usually solid, that underlie the soils or 176 other unconsolidated, superficial material (Neuendorf, 2011). More simply, it shows the lithified 177 rocks that lie under the loose softer material at the surface of the Earth, up to an unspecified depth.

178
The island is divided into four geological terranes, the Kyrenia or Pentadhactylos Range, the 179 Troodos For the needs of our research a 3D geological model was created (Fig 2), in order to visualize the 190 study area, examine its thermal response and the potential of vertical GHE usage in the area. Company. The Lakatameia area was penetrated by the exploration hole KL-1, which was drilled 199 by the Forest Oil Company in the 1960s. This borehole, which reached a depth of 2400 m, is shown 200 on the cross section of the geological model (Fig 2).

202
The 3D model was designed based on "Bedrock Geologic Map of the Greater Lefkosia Area" 203 (Harrison et al., 2008) with the use of the ArcGIS and Adobe Illustrator software. In more detail, 204 ESRI ArcGIS package of software which provides the ability to visualize geological data in 3D 205 was used, for creating the surface of the model. Then, using Adobe Illustrator graphic software, the 206 cross section was "attached" to the 3D geological surface already created in ArcGIS.

207
On the 3D model 9 different geological formations are found as shown below:   In addition to the above, Nicosia Formation is divided into seven geological members: were measured after the samples were dried in an oven at 110±10° C for 24 h (ASTM-C-332).

268
ISOMET 2104 is a portable commercial apparatus for direct measurement of heat transfer . The apparatus is equipped with two types of measurement probes: needle 275 probes for soft materials, and surface probes for hard materials.

276
For conducting measurements with the surface probe of the Isomet 2104 Heat Transfer Analyzer a 277 smooth flat surface was required. In order to create the flat testing surface, a rock cutting machine 278 and a polishing machine were used. At least two flat testing surfaces with a minimum size of 7×7 279 cm 2 were created. The minimum thickness to enable the test was 2.5 cm.

280
In order to calculate the specific heat capacity cp of the samples, their density ρ was also measured 281 in the laboratory (specific heat capacity = volumetric heat capacity / density). Because most of our 282 samples were of irregular size, densities of the collected samples were defined by laboratory tests 283 based on CYS-EN-13383-2:2011. In cases of small in size or soluble in water samples, volumes 284 were measured using the Displacement Method (Archimedes Principle).

287
Laboratory tests for measuring the thermal properties of the collected geological samples were 288 repeated two or three times on each sample. Mean values were calculated based on the geological 289 type of the tested samples (Table 1).

290
Totally 16 samples were collected from the Nicosia Formation, 3 from the Lapatza Formation, 3 291 from the Kythrea Group (excluding Lapatza Formation), 3 from the Apalos Formation and 1 from 292 the Kalavasos Formation.

294
Care was taken in the collection of samples to represent the geology at the corresponding depth in 295 the boreholes. For this reason, surface samples were excluded. Instead samples were collected from 296 nearby road cuts, cliff or sloping sides to make sure that collected samples went the same geological 297 processes as the corresponding borehole strata. 298 299 Table 1 shows that the mean values of the thermal conductivity of all the samples are between 0.56 300 and 1.48 W m −1 K −1 . Calculated values for specific heat capacity are also in a very small range 301 from 0.60 to 1.02 ×10 −3 J K −1 kg −1 .

303
Results are further analyzed in Figs. 3 and 4, where one can see that each geological formation 304 presents a range of values for each thermal property. This is due to the variety of lithologies that 305 are present in each geological formation (see Section 2) and the specific properties of the samples 306 measured. Sample properties like grain size, amount and type of impurities, geological compression 307 when the sample was formed and many others affect greatly the thermal properties of the sample.  (Table 1).

316
The Thermal Conductivity Map in Fig. 5 shows low values for the biggest part of the map and only 317 a small part of the study area has values above the average. On the other hand, the Specific Heat 318 Capacity Map (Fig. 6) shows the biggest part of the study area with an average value of 319 0.6×10 −3 JK −1 Kg −1 . In both cases this is due to the fact that a large part the Bedrock Map is covered 320 with Marl Member of the Nicosia Formation which has poor thermal response. 321 322 All

333
As already said, the main objective of this research is to provide engineers with a method to be 334 followed for sizing vertical GHEs in a study area. This will be achieved through setting up a large 335 number of representative study cases covering the extent of the area of interest. In our example 22 336 representative study cases were chosen in different geological and geographical locations, as can 337 be seen in Fig. 7 (some of them appear in Fig. 2 as well). In each case the influence of ground on 338 the performance of the GHE was investigated so as to calculate the heat load per meter depth that 339 can be transferred to the ground in each borehole of the studied area. The performance of the GHEs was calculated using a validated model in FlexPDE software. The 344 results are tabulated in Table 3 and can be used as a guide for engineers designing GHEs in the 345 Greater Lefkosia Area

346
The FlexPDE software was chosen as the numerical solving tool, because it is a script-driven 347 program that can combine parameters -i.e. equations and boundary conditions -entered by the 348 user and obtain numerical solutions of 3D differential equations, based on the Finite Element 349 Method. For the cases under study, heat distribution over time is described by the general heat 350 transfer equation based on the energy balance. The rate of energy accumulated in the system is 351 equal to the flow of energy entering the system, plus the rate of energy generated within the system, 352 minus the flow of energy leaving the system.

353
The 3D conservation of the transient heat equation is used and applied in FlexPDE software: 354

355
(1) where is the temperature, is time, is the density of the borehole/surrounding soil and rock 358 material, is the heat capacity of the borehole/surrounding soil and rock material at constant 359 pressure, is the density of the heat transferred fluid, is the heat capacity of the heat transferred 360 fluid at constant pressure is the velocity of the groundwater, is the heat source and is Fourier's 361 law heat conduction that describes the relationship between the heat flux vector field and the 362 temperature gradient: where λ is the thermal conductivity of the borehole/soil/rock material. 367 368 In Eq. (1) the first term represents the difference in internal energy, the second term is the part of 369 the heat carried away by the flow of water, the third term represents the net heat conducted as 370 described in Eq.
(2) and the fourth is the heat source.
where is the effective density and volumetric heat capacity of the porous media at The heat conduction in Eq.
(3) can be expressed as: where is the effective thermal conductivity.

393
In our case a weighted geometric mean of the thermal conductivity of both the solid and the fluid 394 materials is needed, and is given by: 395 Note that the numerical model, as described above, has been validated by Florides et. al. (2013).

406
In our study cases we refer to a geothermal system combining a borehole heat exchanger and the 407 surrounding soil mass crossed by an aquifer. A vertical GHE was set up in each borehole consisting 408 of a descending and an ascending leg of polyethylene pipe connected at their ends with a U-joint. 409 Boreholes (BHs) have a diameter of 0.2 m and were filled with thermally enhanced bentonitic clay 410 (Fabien et al., 2011). Bentonitic clay has the ability to expand and completely fill the borehole and 411 hold firmly the GHE in place. Water was used as the heat carrier fluid, circulating in the tubes. In 412 the analysis the area considered was equal to 0.5 m around the borehole and a depth of 100 m was 413 assumed. The tubes used are of 100 m length, 0.0285 m inner diameter and 3.5 mm wall thickness. 414 The distance between the center of the tube and the center of borehole is 0.048 m. The initial 415 temperature of the ground is set to 22° C for the entire study area (based on temperatures measured 416 inside Lakatameia BH) and the temperature of circulating water at 40° C in order to satisfy the 417 requirements of the heat pump. The heat pump type was chosen in accordance to the results of the 418 Technical Requirements Checklist (TRC) test that took place again at Lakatameia BH. Geothermal 419 borehole basic parameters as used in the simulations are presented in Table 2.

420
The knowledge of the geological and hydrological (water level) conditions were necessary for the 421 correct parameter setting. Fig. 8

428
Geological changes for the surface layer (up to 7 m depth) were not taken into consideration, as 429 some boreholes (SHN 1, SHN 2, SHN 16, etc.) were very close to rivers and Alluvial Deposits 430 were found on top. In some other cases (SHN11, SHN 12, etc.) manmade materials were found on 431 the surface layer. In addition, surface layer is affected from seasonal ambient temperature and water 432 flow resulting from rainfall, factors which were neglected in this research.

445
In this section the output results of FlexPDE software are presented. Fig 9 shows  where is the amount of heat energy lost by the circulating water in the tubes, the mass of water 453 that was circulated, the specific heat capacity of water, the final (output) temperature of 454 circulating water and the initial (input, entering) temperature of circulating water. 455 456 457 In Fig. 9 GHEs heat load values are monitored for the first 24 working hours. Values decrease with 458 time, as the difference between input and output temperature of circulating water also decreases. 459 Note that SHN1, SHN12, SHN14 and Aglatzia BH, although they have the same lithology with 460 Lakatameia BH, they exhibit different heat loads. Heat load changes are proportional to the changes 461 of water level in these cases. In the cases where heat exchanger surrounding material and water 462 level in the borehole are the same, density is the key factor for the prediction of the thermal response 463 of the borehole. For example, from our study cases, SHN1, SHN12 and SHN14 have the same 464 lithology and almost the same water level with SHN15, but the heat exchanger surrounding material 465 (despite the same lithology for all cases) have higher density values at the area where SHN15 is 466 located. Heat loss results for each study case are analytically presented in Table 3 after 12, 18 and 467 24 hours of continuous operation.

469
It turns out that SK-5 exhibits a much lower heat loss value compared to all other boreholes located 470 at the north part of the study area (see Fig. 7). This is due to the Gypsum lithology that is present 471 at the drilling location. Gypsum has a unique thermal response compared to other lithologies, as it  The heat loss data used in generating the "Design Load Map of Ground Heat Exchangers for the 479 Greater Lefkosia Area" (Fig. 10) resulted from calculations as explained above. Calculated data for 480 the 22 study cases and for GHE operating for 24 hours (Table 3), were interpolated with ArcGIS 481 software using the Kriging method. This method was chosen as it is an interpolation method based 482 on a statistical model that includes autocorrelation. Because of this, Kriging geostatistical technique 483 not only has the capability of producing a prediction surface but also provides some measure of the 484 certainty or accuracy of the predictions. The digital elevation model presented on the background 485 of the map (Fig. 10) was also created with the use of ArcGIS software supplied by licence of the 486 Cyprus Geological Survey Department.

487
The values that resulted from interpolation and are presented on the map, vary between 26 and 42 488 W m -1 . It is important to note that the heat loss data used for interpolation are not evenly distributed 489 throughout the map. Small areas with dense data surrounded by wide areas containing sparse data 490 are characteristic of the heat load map. As such, the user should keep in mind that the heat loss map 491 is well constrained in areas of dense data and more interpretive in areas of sparse data.
492 Fig. 10 can be used as a tool for visualizing and easily understanding the potential of GHE usage 493 in the Greater Lefkosia Area. From the map it can be visualized that, generally, the heat load 494 transferred to the earth through the GHEs at the north part of the study area is much higher than the 495 transferred heat load at the southern areas. This is due to (a) water levels that are higher at the north 496 part, and (b) the fact that different geological formations are present in the two areas, as they are 497 part of different geological terranes (see section 2, "3D Geological Modeling of the Study Area") 498 499 500 501

503
The following section gives some examples of how similar approaches to evaluate the suitability 504 of a location/area for GSHPs was applied in other studies. Most of them demonstrate a series of 505 maps: heat flow, estimated temperatures, geological maps (geology, geotechnical properties) and 506 underground water level maps in comparison with GHEs suitability maps. 507 508 Maximum and minimum calculated values of similar projects are presented in Table 4 below. Heat 509 exchange load values of GHEs are in a range of 22 to 63 W m -1 and although a similar approach 510 can be adopted in all cases for closed loop systems, the final length of the GHE varies depending 511 on the country/area and the site-specific geology (geology, water level, porosity, etc.). All values 512 can be used only as an indication, as underground site-specific conditions may vary in each case. 513 The large scale web applications developed for France, Ireland, England and Wales do not take 514 into consideration the local special underground conditions and, therefore, they classify the areas 515 only as favorable and less favorable, without giving an actual heat exchange value.

516
A number of projects presenting thermal conductivity maps can also be found in international 517 libraries. Thermal conductivity maps were also proposed as a tool for evaluating the suitability of 518 a place for GHE installations. Thermal conductivity values illustrated in Table 5 clearly show that  519 there is no connection between the thermal properties of rocks that are found in different countries 520 and for each country/area separate measurements of thermal conductivity should be taken. 521 522 523 524 525

527
The current research has focused on a methodology of measuring and analyzing the thermal 528 properties of the lithologies encountered in an area, which can then be used for the prediction of 529 heat injection rates of a GHE in a selected area. The proposed methodology can be applied in any 530 vertical GHE as it describes (a) the full procedure of sampling and thermal testing, (b) the 531 compilation of thermal conductivity and thermal diffusivity maps with the use of GIS, (c) the basic 532 formulas used for calculating the geothermal response of a vertical borehole with respect to the 533 water lever, porosity and the thermal properties of each lithology present in the borehole, and finally 534 (d) the compilation of a heat load map. 535 536 Regarding the area that was used as a study case, the Greater Lefkosia Area in Cyprus is an area 537 with high potential in usage of vertical GHEs due to its geological conditions. The "Design Load 538 Map of Ground Heat Exchangers for the Greater Lefkosia Area" was compiled in order to identify 539 areas favorable for installation of Ground Heat Source Pumps and to provide engineers with a 540 useful guide for sizing vertical GHEs. Values presented on the map are between 26 and 42 W m -1 541 for GHEs of up to 100 m depth, after 24 hours of operation in cooling mode. 542 543 For the compilation of the map, the thermal properties of the ground were required and geological 544 sampling for collecting and measuring samples was carried out. Recorded values are in a range of 545 values for each thermal property of each geological formation. Thermal conductivity λ for the 546 encountered geological formations varies from 0.5 to 1.5 W m -1 K -1 , while specific heat capacity 547 cp from 0.6 to 1.0 J K -1 kg -1 . Each geological formation consists of several lithologies, 548 characterized by different degrees of weathering and fracturing. In the same formation these Values may also vary due to the different types of impurities that each sample may have. Both 551 properties, thermal conductivity λ and specific heat capacity cp are presented in separate maps, 552 which give the geographical aspect of the thermal properties in the area. 553

554
Water level was also considered in calculations, as the flow of underground water has an important 555 effect on the cooling of the vertical heat columns of the heat exchangers. Heat load values that can 556 be transferred to the ground through GHEs are proportionally related (a) to the level of 557 groundwater, and (b) to the density of the underground where the borehole is located. As the level 558 of water and the density of the soil/rock increases, heat exchange rate increases too. The only 559 exception is the borehole where Gypsum is present, with heat loss decreasing as the water level 560 increases.

561
Additionally, the heat load transferred to the earth through the GHEs in the north part of the study 562 area is slightly higher than the amount of heat loss occurring in the southern areas. This is due to 563 the level of water, which is higher in the north part. Also, this change of heat load values from the 564 North to the South that can be observed on the "Design Load Map of Ground Heat Exchangers for 565 the Greater Lefkosia Area" (Fig. 10), can also be interpreted as the boundary between the Kyrenia 566 terrane and the Circum Troodos Sedimentary Succession.

567
Upon reviewing similar studies and GIS applications referring to other countries, it turns out that 568 the heat exchange load values of GHEs vary from 22 to 63 W m -1 . The magnitude of the values is 569 a major factor for deciding the depth of the GHE. It must be stressed though that these values can 570 only be considered as indicative, since the place site-specific conditions may vary in each case. The 571 reviewing process has shown that there is not necessarily any connection between the thermal 572 properties of rocks that are found in different countries, and that for each country/area a separate 573 measurement of thermal conductivity should be performed.          areas are classified as: low (<0.9 Wm −1 K −1 ), medium low (0.9-1.1 Wm −1 K −1 ), medium (1.1-1.3 Wm −1 K −1 ),, medium high (1.3-1.5 Wm −1 K −1 ), and high conductivity (≥1.5 Wm −1 K −1 ) areas Vijdea et al., 2014 831 832 833