Next Article in Journal
Blue-Light Levels Emitted from Portable Electronic Devices Compared to Sunlight
Next Article in Special Issue
Explicit Multipole Formula for the Local Thermal Resistance in an Energy Pile—The Line-Source Approximation
Previous Article in Journal
Performance Comparison of PD Data Acquisition Techniques for Condition Monitoring of Medium Voltage Cables
Previous Article in Special Issue
Effects of the Circuit Arrangement on the Thermal Performance of Double U-Tube Ground Heat Exchangers
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

Reviewing the Modeling Aspects and Practices of Shallow Geothermal Energy Systems

1
Faculty of Engineering and Technology, Cyprus University of Technology, 31 Arch. Kyprianou, P.O. Box 50329, 3603 Limassol, Cyprus
2
Geotechnics Department-National Laboratory for Civil Engineering, 1700-066 Lisbon, Portugal
3
Slovenian National Building and Civil Engineering Institute, 1000 Ljubljana, Slovenia
4
EKIT Department, University of Plovdiv “Paisii Hilendarski”, 4000 Plovdiv, Bulgaria
5
Department of Mechanics, Technical University of Sofia, Plovdiv Branch, 4000 Plovdiv, Bulgaria
6
Department of Electrical Engineering, Computer Engineering and Informatics, Cyprus University of Technology, 31 Arch. Kyprianou, P.O. Box 50329, 3603 Limassol, Cyprus
*
Author to whom correspondence should be addressed.
Energies 2020, 13(16), 4273; https://doi.org/10.3390/en13164273
Submission received: 11 July 2020 / Revised: 5 August 2020 / Accepted: 14 August 2020 / Published: 18 August 2020
(This article belongs to the Special Issue Computational Geothermal Energy Applications)

Abstract

:
Shallow geothermal energy systems (SGES) may take different forms and have recently taken considerable attention due to energy geo-structures (EGS) resulting from the integration of heat exchange elements in geotechnical structures. Still, there is a lack of systematic design guidelines of SGES. Hence, in order to contribute towards that direction, the current study aims at reviewing the available SGES modeling options along with their various aspects and practices. This is done by first presenting the main analytical and numerical models and methods related to the thermal behavior of SGES. Then, the most important supplementary factors affecting such modeling are discussed. These include: (i) the boundary conditions, in the form of temperature variation or heat flow, that majorly affect the predicted thermal behavior of SGES; (ii) the spatial dimensions that may be crucial when relaxing the infinite length assumption for short heat exchangers such as energy piles (EP); (iii) the determination of SGES parameters that may need employing specific techniques to overcome practical difficulties; (iv) a short-term vs. long-term analysis depending on the thermal storage characteristics of GHE of different sizes; (v) the influence of groundwater that can have a moderating effect on fluid temperatures in both heating and cooling modes. Subsequently, thermo-mechanical interactions modeling issues are addressed that may be crucial in EGS that exhibit a dual functioning of heat exchangers and structural elements. Finally, a quite lengthy overview of the main software tools related to thermal and thermo-hydro-mechanical analysis of SGES that may be useful for practical applications is given. A unified software package incorporating all related features of all SGES may be a future aim.

1. Introduction

Despite the potentially huge financial savings and reduction of fossil fuels, the use of geothermal systems, either for electricity production in power plants (see for example [1] and related devices) or for the exploitation of the more widely used shallow geothermal energy, has been suffering from a lack of proper official international standards. Shallow geothermal energy system (SGES may take a number of different forms such as vertical and horizontal GHEs (with different pipe configurations), surface water, standing columns, but also as thermo-active structure (TAS) systems or energy geo-structures (EGS). The latter find applications like diaphragm walls, retaining walls, embankments, tunnel linings, barrette piles, shallow foundations and energy piles (EPs) [2,3] (Figure 1). Some national professional associations have produced guidance documents [4,5,6], while a general guidance proposing principles of SGES design within the framework of existing geotechnical standards (e.g., Eurocodes) has yet to be proposed. In particular, the design of heat exchange elements integration within various geotechnical structures, i.e., EGS, which are needed to support buildings or ground, has not been sufficiently addressed. When designing such foundations, one has to examine the construction processes from the initial phase through all construction stages up to the loading phase, with consequently induced stress changes in the ground. Another consideration, when thermal changes are applied in EGS during the use phase, are the induced additional stress changes in the ground [7]. Hence, to help the engineer dealing with the design of such structures, it is very useful, if not essential, to provide proper tools such as software, in line with appropriate guidelines.
Although, the use of foundation elements for building energy has been around for more than three decades, according to the authors knowledge, only Swiss guideline SIA D0190 [8] considers thermally-activated piles in detail and provides design guidance for thermal effects on differential movements in pile groups and on consequent additional pile loads. Thus, SIA D0190 views EGS from only the thermo-active piles’ perspective. SIA D0190 refers to the numerical tool PILESIM used for the calculation of the thermal performance of energy piles [9].
However, besides thermo-active piles, whose use began in the 1990s [10], base slabs, diaphragm walls, tunnels [11,12,13] and even sewers [14] are nowadays utilized for heating and cooling purposes, while borehole fields have the potential of being used for seasonal thermal storage [15]. The European standard describing the design of heat pump (HP) heating systems (EN 15,450 [16]) does not include design guidelines for ground-source heat pumps (GSHPs) coupled with EP or borehole heat exchanger (BHE) fields. Among national EGS design guidance documents one can refer to the German [17], the Austrian [18], the Swiss [8]), the Italian [19], or the French [20] guidelines.
There also exist other kinds of recommendations to help on the proper exploitation of the ground for geothermal purposes without adversely affecting the ground, groundwater, and interacted structures, by the operation of the geothermal system. Most of them focus on the proper provision of energy from below the earth’s surface via various kinds of heat exchangers [21,22], while some of them consider the mechanical impact of temperature variations on the EGS as well [23,24].
Generally, GSHP system design and related environmental considerations are covered within such guidance documents [17], where the distribution of responsibilities among participating system designers are presented [23]. In addition, the design of thermo-active piles is covered within some of the existing documents [6], while other types of EGS are not included. Regardless of the ultimate/serviceability limit state analysis with pile thermal effects taken into account [23], more useful design methods are still needed to guide the general geotechnical design of EGS in practice [25]. The current geotechnical design of EP is generally dominated by empirical considerations, presenting over-conservative approaches through the enlargement of the factor of safety (FS), leading to substantial extra costs and the subsequent use of non-standard construction methods (e.g., [25,26]). Amatya et al. [27] assume that the FS used for EP is at least twice as large as that used for traditional piles without a heat exchanger, and hence a probabilistic geotechnical analysis of EP [28] has been proposed instead.
COST Action GABI (Geothermal energy Applications in Buildings and Infrastructure) aimed to compile different kinds of thermo-active structure design procedures and prepare a common European framework for design of shallow geothermal energy system [29] integrated into foundation structures.
Motivated by the lack of complete guidelines leading to the design of SGES, the current study aims at reviewing the modeling aspects of such systems. This is done through the presentation of: (i) the main analytical and numerical models, including thermal energy testing (such as TRT), related to the thermal behavior of SGES; (ii) auxiliary factors affecting the modeling such as boundary conditions, spatial dimensions, determining SGES parameters, short-term vs. long-term analysis and presence of groundwater; (iv) very important related information on the thermo-mechanical interactions modeling; (v) very useful practical applications through the use of the main software tools related to thermal and thermo-mechanical analysis of SGES.
Now, to move to GSHP design, it is important to determine the thermal characteristics of the GHE. Although there are several methods developed for the determination of ground thermal characteristics, the most commonly used is the TRT, where the temperature perturbation, resulting by the heating, with an electric element, of the heat transfer fluid flowing in the GHE, is monitored (e.g., [30,31]).
There are various recommendations in literature for the TRT duration, namely for a minimum of 12–40 h, 36–48 h, 48 h with data to be collected for a minimum of 10 minutes’ steps, 50 h, and 60 h. Clearly, the longer the TRT duration the more expensive the TRT. However, short durations may underestimate the thermal conductivity and/or overestimate the borehole length [32]. The parameters or characteristics to be computed through the TRT are, usually, (i) the ground thermal conductivity, (ii) the effective borehole thermal resistivity, (iii) the ground specific heat capacity, (iv) the grout conductivity, (v) the grout specific heat capacity, (vi) the shank spacing [33]. As technology has progressed, the original TRT test [34], has undergone variations such as the enhanced TRT, where supplementary information with regard to thermal conductivity distribution and the accuracy of the thermocouples [35], or the distributed TRT, where the thermal conductivity is measured at multiple depths of the borehole with the aid of a fiber optic cable [36].
Sharma et al. [37] were the first to suggest and field test the use the optical fiber sensors (optrode) on geothermal wells. The authors used the dependence of Raman scattering light strength on temperature and measured the temperature distribution along the optical fiber. The method was then further developed and more systems were tested [38,39,40,41], with the advantage of using having measurements at certain depths of a GHE, also at different times [42], with an interpretation of the data to be made either by using the ILSM [43] or the infinite CSM [41]. Similar to distributed TRT, enhanced TRT [44,45] rely on heat injection at constant power to the GHE and gather data along the depth of the GHE with the use of optical fibers using copper heating cables incorporated into one hybrid cable [46].
For modeling purposes, both analytical and numerical approaches exist in literature. All the models are governed by the Fourier’s law of heat conduction [47,48,49]. The models can also be categorized as: (i) long-term models that use steady-state solutions or long time-steps in a transient solution; (ii) short-term models that use hourly time-step variants, or less, in a transient solution. Alternatively, the models can be categorized according to the type of the ‘source’ heat (linear or nonlinear, finite or infinite), depending on the nature of the problem under investigation [32]. The typical analytical models are (i) the infinite line-source model (ILSM), (ii) the finite line-source model (FLSM) and (iii) the cylinder or cylindrical-source model (CSM). The numerical models are in general developed to represent the borehole geometry in more detail than analytical models. Additionally, the thermal properties of the borehole filling, pipe, fluid, ground, as well as varying heat transfer rates and influences of flowing groundwater can be considered. They can be developed in 1D, 2D, and 3D space in order to solve an energy balance equation and are based on the following three methods: finite differences (FDM), finite elements (FEM), and finite volumes (FVM).
The model input data in most cases come from the excitation function variations in time (heat power injection/extraction input). In some models additional input data are needed, for example a ground surface and ambient temperature functions. The known parameters values are usually the borehole effective length and radius, the U-pipe configuration and dimensions, and the volumetric heat capacity of ground. The model output is the temperature response of the borehole (average fluid temperature). The initial conditions include the undisturbed ground temperature. Regarding parameter calculation, this can be performed by either slope determination or parameter estimation techniques. Both techniques are based on fitting the test data to the model procedure through various available optimization algorithms.
Mechanical effects arising from thermal loads on SGE applications are relevant to all EGS (see above) and should be considered, in addition to structural design loads [46] (Note that BHE are of purely thermal character). As prime function of EGS (e.g., EP, energy wall or tunnel) is load resistance, one should be aware of potential for thermal loads impact on the structure performance with regard to serviceability or safety. Thus, the possible deformation of the building supported by EGS, the stress/strain variation in EGS due to cyclic heating and cooling and any possible changes in load resistance of the EGS ought to be studied [50].
There are currently more than 10 software tools, which are available and capable to model problems of energy gains/losses of SGES. They can be analytical, or response model based, considering the average fluid velocity of the borehole fluid. Some of them use computational fluid dynamics (CFD) modules and can consider the fluid flow of the heating/cooling fluid in more detail. Few of them can also tackle groundwater flow through porous media analysis, while some can be used for thermo-hydro-mechanical analyses. Depending on the software’s abilities, different configurations of GSHPs can be considered such as vertically or horizontally oriented pipes, U-tubes, double U-tubes, coaxial tubes, but also more complex geometries and EGS. The advantage of many software tools is that they can solve the short-time behavior of SGES (e.g., BHE) systems more accurately than analytical models; in particular, those using CFD modules.
The sequel of the paper is as follows. In Section 2, the main types of models for the design of SGES are overviewed. The boundary conditions, the effects of spatial dimensions, the techniques for the determination of SGES (borehole) parameters, a short-term vs. long-term analysis and the effects of groundwater in the system are discussed in Section 3, Section 4, Section 5, Section 6 and Section 7 respectively. Section 8 presents in a detailed manner related thermo-mechanical interactions and constitutive modeling. Then a comprehensive presentation of software tools and their use for solving problems of SGES is given in Section 9. Section 10 concludes with a discussion.

2. Types of Models for Thermal Analysis

2.1. Analytical Models

Generally, models can be divided into analytical (models with response factors included) and numerical models. The main differences between them are computational time, modeling accuracy, ability to define the exact location of the SGES and dynamic heating/cooling plant modeling capabilities.
Various assumptions can be implemented into analytical models [51] for simplification and further use. They are most feasible for evaluation of TRT results to determine soil and pile thermal characteristics. For instance, line or cylindrical source models are used to determine net heat rate per borehole length from the ground from the ground thermal diffusivity. Eskilson [52] introduced so-called ‘g-functions’, which are used to represent the dynamic thermal behavior of a specific borefield configuration. Furthermore, some of the analytical models can be used for transient simulations.

2.1.1. Infinite Line-Source Model

The ILSM was developed by Kelvin [53] and further elaborated by Whitehead [54]. Later, Ingersoll and Plass [47] and Carslaw and Jaeger [48] found an analytical solution in a form that is used today for TRT data analyses of vertical GHEs. Ever since, improved approaches, as regards accuracy, have been introduced in literature (e.g., [55]). This model represents the borehole as an infinite line heat source (sink) put into an infinite isotropic underground region, with the radial dimension of the borehole neglected. Heat transport is done through heat conduction. Performing evaluation formulas, the borehole thermal resistance and the ground thermal conductivity can be calculated or estimated. Then the solution of heat conduction is as follows (e.g., [56]):
T b T 0 = Q ˙ 4 π k Ei ( r b 2 4 α t ) = Q ˙ 4 π k r b 2 4 α t e u u d u
where Tb is the temperature [K] at rb [m], the borehole radius, T0 is the undistributed ground temperature or the initial borehole temperature, Q ˙ is the constant heat injection rate (W m−1), Ei(.) is the exponential integral, t is the time (s), k is the borehole thermal conductivity (W m−1 K−1), and α = k/(ρcp) is the ground thermal diffusivity (m2 s−1), ρ being the density (kg m−3) and cp the specific heat capacity (J kg−1 K−1). Equation (1) can be simplified into an analytic approximation without integrals (see [56]) that also include the borehole thermal resistivity Rb (K m W−1), the mass of the borehole m (kg), and its length L (m).
Note that the classical line-source theory does not consider the end effects of the heat source due to the assumption of infinite length. Ingersoll and Plass [47] recommend that the ILSM must be used for applications with Fourier number larger than 20. In such a situation, a distorted solution (due to the line source assumption) at smaller time steps is avoided [57].

2.1.2. Finite Line-Source Model

As the ILSM does not take into account the finite length of the borehole, the model error significantly rises in long-time simulations (of more than 10 years). The FLSM considers a heat flow rate on the vertical axis with a constant temperature gradient in the semi-infinite region, and a varying ground surface temperature [58]. The use of the FLSM for the duration of 50 h in TRT is not recommended. Eskilson [52] and Claesson and Eskilson [59] first developed an analytical so-called g-function expression, the solution being determined using a line source with finite length, assuming a constant temperature evaluated at the middle of the borehole [60]. The general solution is a function of radial coordinate r [m] and axial coordinate z [m] and is given by
T b T 0 = q ˙ 4 π k 0 L ( erfc ( r 2 + ( L 2 u ) 2 / 4 α t ) r 2 + ( L 2 u ) 2 erfc ( r 2 + ( L 2 + u ) 2 / 4 α t ) r 2 + ( L 2 + u ) 2 ) d u
where erfc(.) is the complementary error function.
There exist variations of the model described above. For example, Zeng et al. [61] proposed a modification, as they realized that the temperature in the middle of the borehole overestimates the reference temperature. Despite this being true, the specific modification is difficult to use for practical applications and does not offer a significant improvement in the accuracy of the solution. Lamarche and Beauchamp [60] adapted Zeng et al. [61] solution using the mean temperature for the evaluation of the analytical g-function in a faster way than Claesson and Eskilson [59].
Other thermal analysis models, specifically intended for composite media can be found, for example, in [62,63,64,65]. Hellstrom [63] described an LSM for a composite region as an approximation of some more general method. The instantaneous LSM for composite media was shown to be appropriate for modeling various GHE configurations and can be used for the analysis of transient heat conduction inside and outside boreholes [62]. A new approach of composite-media LSM assumed both an infinite-line and an infinite-composite cylinder to obtain temperature distributions around two BHEs and two pile-GHEs [64].

2.1.3. Cylindrical-Source Model

Another analytical model known as CSM was presented by Carslaw and Jaeger [48] (see also [66]). In such a model, the borehole GHE acts as an infinite hollow cylindrical heat injection/extraction source. The CSM, of which the ILSM is a simplified variation, may be used to approximate the GHE as an infinite cylinder with a constant heat flux. The general solution [63,67] (as a function of radial coordinate r [m]) is given by
T b T 0 = Q ˙ π 2 k r b 0 ( e α u 2 t 1 ) J 0 ( r u ) Y 1 ( r b u ) Y 0 ( r u ) J 1 ( r b u ) u 2 ( J 1 2 ( r b u ) + Y 1 2 ( r b u ) ) d u
where J0, J1, Y0, Y1 are Bessel functions of the first and second kind.
Under certain circumstances Equation (3) can be approximated to a form containing Rb (see, for instance, [68]). Modifications of the CSM for different evaluation procedures came in [69,70,71,72], for GSHP systems, in [63] for ground storage systems and in [73] for TRT.
In what concerns the spiral pipe GHE configuration, CSM has been the most widely used. When applied with a spiral configuration, it does not take into account the heat capacity of the spiral coil, but considers heat flux directly on the cylinder surface (borehole wall).
Man et al. [67], Man et al. [74] and Cui et al. [75] introduced new variations of CSM, called, depending on the situation, hollow CSM, solid CSM, or for energy piles, finite ring-coil source model, finite spiral source model, infinite ring-coil source model, and infinite spiral source model. In the CSM models the spiral coil is not considered as a spiral, but it is simplified as a continuous cylindrical heat source. In the ‘coil source’ models the discontinuity of the heat source and the spiral pitches, presented as separated rings, are taken into consideration. In the ‘spiral source’ models the coil is presented as a spiral line heat source, the asymmetry on the vertical axis is considered, requiring a 3D temperature distribution. Park et al. [76] validated the modified CSM and the ring-coil source model through experimental results. They found that for a certain spiral coil pitch (200 mm), both the modified CSM and the ring-coil model can be in good agreement with the TRT results for an acceptable overestimation of the heat exchange capacity of the piles, something that did not occur for pitches of increased size.

2.2. Numerical Models

Although analytical methods allow for an accurate prediction of the thermal stress and strain response of EGS (see Section 8), the methods exhibit shortcomings in their ability to simulate effects of cyclic heating and cooling, of transient pore water pressure generation and dissipation, as well as of radial stress changes [4]. Major drawbacks of current analytical tools consist in their inability to take account of the short-term thermal energy storage and thermal resistance of EGS, the effect of ground water flow and real surface boundary conditions [77], which might all affect significantly the mechanical response of EGS (see Chapter 4) whose geometry differs significantly from BHE.
On the other hand, numerical models have the ability to achieve very high detail of modeling accountancy for (i) phase change in soil due to freezing, (ii) variable boundary temperature conditions, (iii) pile internal components, (iv) thermal capacitance of pile inner materials, and (v) different configuration of piles (e.g., single/double/multi U-tube), etc. [78]. In principle heat exchangers with fluid, borehole walls and grout need special consideration, whereas the soil part (together with EP, if present) is considered separately, where both models are coupled via heat flux at the borehole boundary.
Numerical simulation offers also the possibility to assess quantity-wise the difference between 2D, 3D, and quasi 3D models. However, the cost for a detailed modeling consist of high computation time and modeling complexity, both barely applicable for everyday engineering design needs.
Analytical models consider soil and pile (when present) regions as a homogenous medium, as opposed to most of numerical models that can account for a multi-layered structure of ground. Some models, both analytical and numerical, can also account for water advection in the soil region.
It is worth noting that analytically and response factors-based software programs can instantly plot the results of a (say) 25-year period simulation, while numerical software programs can take days to finalize a similar simulation of large pile fields (see, e.g., [79]). Most of the numerical model-based software tools are designed in such a way that they are more precise and capable of calculations of (sub-)hourly time-steps for exact piles location. They can be coupled with whole building simulation software programs, enabling the user to model highly complex dynamic heating/cooling regimes. Now, in their original form, most of the analytical models can only simulate constant heat flux along the pile depth, as opposed to the numerical models that can handle a varying heat flux. Furthermore, software with whole building modeling capabilities exhibit a longer learning curve as opposed to other self-standing software packages designed for ground source heat exchanger simulations only [80]. For more, see Section 9.
Indicative and in no way exhaustive examples of 1D, 2D, and 3D numerical models for the thermal analysis of SGES are presented below.
An example of a 1D numerical model was presented by Shonder and Beck [81] in the framework of GPM (Geothermal Properties Measurements), a method developed for the Oak Ridge National Laboratory. By considering the two pipes of a U-loop as a single cylinder, the GPM model is governed by the CSM. The thermal resistances present in the real problem can be accounted for by adding to the model a thin film with resistance but no heat capacity, as well as a layer of grout that can have different thermal properties than the surrounding soil. Moreover, the model accounts for time-varying heat inputs, as it represents borehole transient heat conduction.
Based on the original Eskilson [52] idea so called state-space models were introduced by various FDM/FEM reductions to obtain a lower-order (1D) model [80,82]. Their authors claim them as very suitable for optimal control of a system.
A 2D model (in polar coordinates) was presented in [83]. The influences of the ambient air temperature, the underground temperature and moisture distribution over the borehole length are added into the model equations by superposition. The model is realized in the TRNSYS (trans-solar system simulation) software environment (see Section 9) and is a method of weighted residuals with aspects of FDM and FEM.
Yavuzturk et al. [84] and Yavuzturk [85] described a 2D radial-axial FVM model to match the TRT early-time data that are influenced by borehole effects. The authors developed an automated grid generation algorithm, in order to represent different borehole diameters, different U-tube sizes and different size of shank spacing. The model takes as inputs the specific geometric information and time series of heat injection rates.
Austin et al. [86] also proposed a 2D numerical model for borehole GHEs. A 2D (in polar coordinates) FVM is utilized for TRT interpretation with parameter estimation techniques. To avoid complex meshing, the actual cross section of the GHE is modeled by simplifying the U-tube pipes to pie-sector shaped elements. The heat flux over the pipes is assumed to be constant.
A different model in 2D in which the soil is divided into horizontal layers of uniform temperature was developed by Bojic et al. [87], whereby the convection heat transfer from the air to the soil and the solar irradiation were calculated. The heat flux between neighboring soil layer was also considered.
Bi et al. [88] constructed a 2D cylindrical coordinate FDM model of the vertical double spiral coil GHE to compare with the experimental TRT results.
A fully 3D FDM in a simple Cartesian coordinate system was described in [89]. The round pipes were replaced by square pipes of equivalent areas. Then, the influence of different layers in the ground, concrete foundations and insulation could be evaluated. Heat transfer in the pipes was dominated by convection in the axial direction, but it was also coupled with the temperature field in the ground through the boundary condition on the pipe surface.
Another 3D numerical model was developed by using the FEM code LAGAMINE in order to simulate an in-situ TRT [90,91]. It includes the feeder pipe influence. The ground was simulated with a 4-node 3D FEM to a depth of 220 m and a surface area of 0.11 km2 (80520 nodes). The heat loss through the building’s foundations was simulated by imposing a constant temperature over time.
More recently, Amanzholov et al. [92] developed a numerical technique, using a 3D FEM numerical model, which considers heterogeneities of the subsurface layers, underground water flows as well as conductive and convective heat transfer in porous media. Note that the model was elaborated in the COMSOL Multiphysics modeling software environment (see Section 9).
Other 3D models include the FEM of Breger et al. [93] in their study of thermal analysis of energy storage in the ground using U-tubes and boreholes and of Ozudogru et al. [94] in their study of vertical GHE, the FVM of Li and Zheng [95] for the development of vertical GHE modeling, and the work Gao et al. [96] for the assessment of the thermal performance of vertical energy piles.
There have also been some attempts to use artificial neural network models [97] or machine learning algorithms [98] for modeling the geothermal performance of SGES in 3D. Although they can accurately model existing complex SGES, they seem to be not very useful in long-term applications as their training requires long-term data.
Finally, of course, many numerical studies on thermal response modeling have been performed through the use of software tools, such as PILESIM 2 in [99], TOUGH(REACT) in [100], TRNSYS in [101], and IDA-ICE in [102]. For more, see Section 9.

2.3. On Models Comparison

There does not really exist an extended literature with regard to the comparison between the various analytical and numerical models as each method is suitable for problems of a specific type, with their pros and cons. An attempt for comparing models can be found, for example, in Aresti et al. [32], in the framework of a more general review study though.
Such attempts include comparisons between: (i) the ILSM and the CSM against experimental data [60,103,104,105]; (ii) design programs using the ILSM and the CSM [106]; (iii) a new analytical model and the ILSM and a composite model proposed [107,108]; (iv) the modified CSM, the ring-coil source model [76]; (v) the ILSM, a FEM based on a multipole method and a FVM [65,109,110]; (vi) analytical models, a FEM and a simulation model [93].
Further comparisons fall beyond the scope of the current study. Indeed, they constitute an important future task for researchers.

3. Boundary Conditions

The choice of boundary conditions majorly affects the predicted heat flows within the SGES and influences its outcomes. They can be defined via temperature conditions or via heat flows on boundaries. For instance, the simplest infinite LSM (1D axisymmetric) needs only initial ground temperature and heat flow at the energy geo-structure—the heat exchanger—as boundary conditions, with heat flow in the axial direction set to zero. Multidimensional models involve spatial temperature conditions as well as heat flows in various directions [51].
It is essential for the accurate estimation of the cyclic heat accumulation and heat extraction potential of SGES to consider proper boundary conditions. Particularly, the seasonal thermal storage and the additional interference between SGES ground surface boundary and the building floor temperature conditions must be considered, to avoid an unbalanced operation of SGES [102]. Such impacts might be less significant in case of deep BHEs, because only the depth near the ground surface is penetrated by the ambient surface temperature variations, and the ground surface boundary is usually modeled as adiabatic. Consideration of seasonal ground surface temperature variation becomes essential in case of shallow EP, whose whole length might lay within this climate-affected subsurface layer [78]. In general, seasonal temperature fluctuation affects ground surface temperature to the depth of 10–20 m. This surface temperature fluctuation can be modeled either as a time-varying temperature (i.e., a Dirichlet condition), or as a time-varying heat flux (i.e., a Neumann condition) [5,12]. Additionally, alternative boundary conditions must be considered in case of SGES interaction with other heat sources or heat sinks that may exist in the ground in urban areas between buildings, sewers, tunnels, etc.
Boundary conditions, such as ground surface temperature and far-field temperature are supposed to be constant in analytical models, while they can vary in numerical models. Also, as previously mentioned, most analytical models in their original form can only simulate constant heat flux along depth, while numerical models can handle varying heat flux [78]. These characteristics become very important when comparing modeling of EGS with various ground depth penetrations, e.g., EP vs. BHE. Similarly, special attention on boundary conditions is needed when additional solar assisted thermal storage is integrated within SGES [111].
De Rosa et al. [112] compared a steady-state model and a novel dynamic model in a real SGES (both implemented in TRNSYS). The inlet water temperature and the water mass flow rate, coming directly as experimental measurements (experimental setup in Valencia), were taken as boundary conditions in the BHE, which is often the case when comparing a mathematical model with an experiment.
Thoren [113] analyzed two case studies, namely the IKEA building complex in Uppsala, Sweden, and the Marine Corps Logistic base in Albany, Georgia, USA. In both cases heating and cooling load measurements on 2 boreholes were used as input for the simulations, done in TRNSYS.
Different considerations of boundary conditions were presented in the case of two common design approaches, namely the ASHRAE analytical method [71] and the GLHE-Pro commercial tool (a design software for ground loop vertical BHE) (see Section 9), based on g-functions, valid also for short-time calculations that were taken into account [114]. The two methods were used to design a BHE field for SGES in two case studies, namely a small-scale residential building (equipped with fan coils) and a medium-scale commercial building (equipped with radiant panels). ASHRAE requires peak loads only for the design heating/cooling month, while GLHE-Pro allows entering peak loads for each month together with the peak load duration. Partial load operation is considered in ASHRAE, while only full-load operation curves are input in GLHE-Pro. Heat pump inlet and outlet temperatures in design conditions are required for the ASHRAE approach, while the input for GLHE-Pro are the seasonal minimum and maximum temperatures (boundary conditions).
In principle, the boundary condition on the surface of the ground is not adiabatic in the case of EP, while the building floor temperature must be considered. Considering proper boundary conditions (Figure 2), a detailed model can be prepared with incorporated heat transfer between the fluid in the downward/upward pipes and the EP, between the EP and the surrounding soil, through the floor structure, between the EP and other heat sources/sinks and from ground surface to depth. Fadejev and Kurnitski [102] reported a 3D borehole model, validated with the actual borehole measurement data. BHE and EP solutions were in detail modeled in a whole building simulation software IDA-ICE (see Section 9). Different ground surface boundary conditions of EPs and a field of boreholes resulted in a more efficient performance of EPs by 23%, for the same field configuration. This result emphasizes the importance of accounting for the heat transfer through the floor structure when analyzing an EP case. Thus, software should support varying boundary conditions, which are mainly appropriate for software tools based on FEM and FDM, but they become problematic in the case of tools based on analytical functions or response factors, where adiabatic ground surface boundary condition must be used.

4. Spatial Dimensions Effect

Many of the analytical solutions assume an infinite length for the heat exchanger to simplify the calculations. However, with shorter heat exchangers (e.g., EP), the end-effects must be considered, particularly in the long term [5].
When using analytical models, thermal interactions between EP or BHE in group can be generally summed up by superposition. In standard TRNSYS modules, the temperature of each point in the storage region of the duct ground heat storage (DST) model is calculated by superposition of a global solution (via FD), a steady-flux solution (g-functions), and a local solution (via FD) (e.g., [112]).
On the other hand, solving the heat balance equations for an EP/BHE group can be conducted numerically too. Varying spatial thermal characteristics are enabled in such a case. Depending on the case, this kind of full numerical analysis might lead to excessive computational demands. Al-Khoury et al. [115] presented a special FEM technique for double U-pipe BHEs and surrounding soil mass with the focus placed on describing the capability of a BHE model to simulate 3D heat transfer processes in multiple heat pipe exchangers, embedded in multi-layer soil mass, using very low number of FE. The heat exchanger was model as 1D, whereas soil as 3D. The model was validated against measured data on a SGES consisting of seven double U-pipe BHEs inserted vertically in the ground. Only input and output fluid temperatures were measured. Numerical tests showed that for 16 BHEs of 100 m length, the system with soil can be modeled with 2210 elements for soil and 16, 2-noded, space-time elements for each individual BHE. Such a number of FEs is very low when compared to typical 3D FE simulations of such systems with millions of FEs. The procedure was implemented in the FEFLOW software [116] (see Section 9).
Thoren [113] tested and compared a new TRNSYS module (246) with experimental data (Uppsala case study). Module 246 is based on a load aggregation scheme and uses g-functions to take care of the borehole field geometry. The load aggregation scheme is based on three blocks, one large, one medium and one small block, regarding the load and time period. The number of individual heat extraction/injection rates to be added in one small load block, the number of small blocks that will be added in one medium block and the number of medium blocks to be added in one large block are set in the parameter properties of the components. On the other hand, in the Albany case study [113], the comparison of serial versus parallel coupled BHE fields showed relatively small differences in performance when simulated with type 557b for a specific case.
Xuedan et al. [117] simulated with TRNSYS the yearly performance of a library building located in Harbin, China, equipped with a GSHP system. To obtain 30-year operation conditions of the reference building, a TRNSYS model—based on hourly-load simulation results—was developed for the analysis of the underground thermal balance of the entire system. Simulations showed that the annual average storage temperature, the annual highest and lowest inlet temperatures, and the factors that influence the efficiency of the system, changed with the number of boreholes. Also, the optimization of the length of heat exchangers was studied. It turned out that more boreholes gave a better performance of the system. Temperature in the stored region increased by 8 °C (when cooling) in 30 years for 120 boreholes, whereas by 3 °C when 240 boreholes were used. Therefore, the overall performance depends significantly on the proper inclusion of all parts of the system into the model.
He [118] applied a 3D numerical model to investigate the characteristics of heat transfer of a BHE and the surrounding ground at both short- and long-time scales. Using the same numerical method, by also implementing a 2D model and comparing the results with those of the 3D model, the most significant 3D effects were identified and quantified. The models were based on an in-house model, called GEMS3D (General elliptical multi-block solver in 3D). Both the 3D and 2D models were validated against experimental data. The predicted outlet fluid temperatures showed a good agreement to the measured outlet fluid temperature for the simulation period. The delayed responses that were demonstrated in the experimental data were captured only by the 3D model, due to the explicit modeling of pipe fluid transport.

5. Methods and Techniques for Determination of SGES Parameters

Here TRTs, being a SGES modeling tool extensively used, serves as a base study toward understanding which techniques can be applied for the determination of SGES (borehole, here) parameters. The ideas discussed below can easily be extended to other SGES such as energy piles, and so on.
The philosophy of performing a TRT requires first putting the borehole from a steady state to a transient condition and then analyzing its thermal response. In theory, a system response is defined as a response of the system output when an impulse response function (Dirac function) is applied (as an excitation) on the system input. In practice, however, the above becomes an impossible scenario because it requires the application of an infinite thermal power to the borehole within a near zero-time interval. So, in all practical cases other than Dirac, excitation functions are used. Then, the same excitation functions must be applied to the mathematical TRT model toward the borehole thermal properties determination.
The thermal excitation may be conducted by various methods such as injection of heat at a constant rate during the test propagation, cyclic pulse injection of heat at a constant rate followed by a pause with zero heat injection or cyclic pulse injection and then extraction of heat at a constant rate.
A separate indication of the ground thermal conductivity can be obtained during the monitoring of the thermal recovery. The advantage here is that heating rate fluctuations than can occur during the TRT can be overcome, depending on the electricity source used for the heaters and on whether the surface pipework is sufficiently insulated to prevent the ambient air temperature from influencing the results.
The step-pulse excitation approach was firstly applied in TRT by Mogensen [34] and Mogensen [119]. Nowadays this is the most commonly used technique during in-situ TRTs. It requires keeping an injection heat rate at a constant level trough the all test duration time by regulation.
The evaluation of TRTs with several heat injection input step pulses is desirable and useful for several reasons. For example, the evaluation of power injection, extraction and recovery pulses can give information about groundwater flow and ground heat capacity. Also, the application of this test design can be necessary if the evaluation of a single heat pulse TRT is not valid for any reason, e.g., due to non-constant mass flow in the GHE. Applying the step-pulse method can allow for the repetition of the TRT in only small waiting periods [85,120]. The TRT evaluation includes the original, i.e., invalid, test, the recovery phase, and the repeated TRT. For an indicative example of such a test design the reader is referred to [121], where a repetition of an invalid test after a waiting period with recovery is shown.
For determination of borehole parameters, the so-called slope determination technique relies on the solution of the LSM problem. The temperature field equation, as a function of time and radius around a line source with constant heat injection rate (having sufficient thermal insulation of the installation and pipes to avoid thermal losses during TRT) can be used as an approximation of the heat injection from a GHE. It follows that the thermal conductivity can be determined from the gradient of the straight line resulting from plotting the fluid temperature against the logarithm of time. Then the values of other parameters such as borehole thermal resistivity can be determined using analytical formulas in which mean or average values are used (see Section 2). For a more interval-independent evaluation technique one needs to use a fitting function for the data, with the thermal conductivity and the borehole thermal resistivity remaining as the two variable parameters. This technique allows one to apply time variable power input to the model and to consider an ambient temperature change and underground water flow influence. The number of the estimated parameters depends on the model and may vary from 2 to 5 or more [122]. During the estimation procedure the simulation of the model runs repeatedly with different parameters set. Searching the best fit of the outputs, different optimization techniques can be applied such as: (i) the simplex method [123], (ii) direct search [124], (iii) gradient methods [125], (iv) evolutionary and stochastic search algorithms, and so on. Several examples of evaluation procedures, using parameter estimation techniques are listed below.
In Roth et al. [126] the commercial software Origin6 (https://www.originlab.com/) was used for the TRT data analysis. The nonlinear curve fitting to user input functions was performed by a Levenberg-Marquardt iteration algorithm. At each iteration, the variance–covariance matrix is computed using the previous iteration matrix. This matrix depends on the dataset assignments, the number of parameters and, of course, the fitting function; it becomes unusable once any of these properties are altered and the fitting session must end.
In Hemmingway and Long’s work [127], a technique (GPM) was applied with some success on the short duration TRT data from two mini-piles. The GPM model was developed by Shonder and Beck [128] to perform a method based on parameter-estimation, in combination with a 1D numerical model developed by Shonder and Beck [81]. The tool solves the equations directly and is hence theoretically valid at short-time periods when simpler analytical models are unusable. This procedure is well suited with borehole TRT data.
Wagner and Clauser [122] and Hemmingway and Long [127] presented a parameter estimation approach to overcome the fact that thermal capacity, usually assumed constant, varies within ±20%.
The evolutionary and stochastic search algorithms are optimization procedures, used in parameter estimation techniques. They allow for accurately finding a best fit, if the parameter space is highly nonlinear, having multiple local maxima, and/or discontinuities. An indicative example of applying such algorithms is the method presented by Popov et al. [129] suggesting the use of the input/output black box identification technique for TRT data analysis. The authors used a nonlinear autoregressive exogenous model structure and stochastic search algorithms for parameters’ estimation. Also, artificial intelligence, a genetic algorithm, and particle swarm optimization algorithm were all employed for avoiding local maxima problems. All analyses were performed in the MATLAB environment.
Finally, regarding estimation procedure evaluation, the so-called sequential forward evaluation [121] starts from the test data start time point or from another known point and goes stepwise forward in time. Thus, the last (in the time series) point of the resulting curve of the thermal conductivity estimates (see Figure 3a) is the resulting estimated value. The forward-convergence curve requires the use of the minimum time criterion at the beginning of the test. In an analogous manner, the sequential backward evaluation is carried out going stepwise backward (see Figure 3b).

6. Short Term vs. Long Term Analysis

The thermal storage characteristics, thermal capacity of the system or the thermal capacities of the circulating fluid, the pipe, the grout and finally the surrounding ground, affect mostly the time domain analysis. In the case of BHEs, as the length of a BHE is much larger compared to the diameter, LSM are usually used, with the thermal capacities of pipe, grout and fluid being neglected. On the other hand, the aspect ratio (length/diameter) in the case of EPs is quite small and, thus, the thermal capacitance of EP structure impacts significantly the short-term thermal behavior of the EP. Consequently, TRT of an EP should be carried out much longer compared to a BHE in order to overcome the issue of energy pile thermal capacitance [78]. Moreover, it is recommended that thermal storage is applied to maintain stable long-term operation of BHEs or EPs, which becomes very important, due to otherwise possible pre- or sub-dimensioning of the system [130].
Two different operating conditions were considered by de Rosa et al. [112], who compared experimental data with modeling results. In a step-test, seven hours of operation following by 12 h of recovery in cooling mode, a TRNSYS model gave a too low outlet temperature after longer time, while results of a B2G model, also implemented in TRNSYS, were in very good agreement with measurements (less than 0.1 °C difference). No time shift of temperature increase of the outlet water was implemented in the TRNSYS model (based on a steady-state model), while the B2G model had also an advection term and captured the time delay. A second type of operating conditions simulated a one-day heating performance. In that case B2G captured the peaks much better than the standard TRNSYS model. The observed differences in the response of the two models may not be relevant from a long-term point of view, but the results indicated some shortcomings of the TRNSYS model in the short-term analysis.
Back to the Uppsala case study [113], the thermal analysis result based on the TRNSYS DST model (types 557a, 557b) highlighted the underestimated heat transfer early on due to a poor consideration of the thermal capacity inside the borehole. At the same time, for the borehole model that considered borehole thermal capacity (only in grout), TRNSYS (module 246) overestimated the short-term heat transfer rate. On the other hand, numerical results for a long-term behavior were in good agreement with experimental data for all 3 model types/modules used here (557a, 557b and 246, see TRNSYS Manual Ref).
In a study conducted by He [118], due to the explicit modeling of the fluid transport in 3D model analysis, the time delay during short periods was captured well, which was not the case in 2D models. Long-term (months) response were captured well in both 2D and 3D models. The major disadvantage of He [118] very detailed 3D model was high computational time. The same experimental data were used for long-term analysis as well in the framework of a comparison performed by the TRNSYS, HVACSIM+ and EnergyPlus software tools (see Section 9) [106]. Although relatively large differences were found for months with high cooling and heating demands, the predicted values followed well the experimental data. The results from the EnergyPlus simulation software fitted best the experimental data.
Javed [131] developed an analytical solution valid also for short-term periods. It was compared to the results of a numerical analysis, showing a deviation of less than 0.01 °C. Long term (multi-year) simulation of SGES was also conducted using the same model based on analytical solutions. A difference of less than 0.3 °C was observed between the simulated and the experimental results.
A model, which combined the short time step g-function with the long-time step g-function was implemented in GLHE-Pro [132], as well as in EnergyPlus [133]. The biggest disadvantage of this model (short-time step g-model) was that the fluid inside the pipes was not explicitly modeled but was treated by a heat flux boundary condition at the pipe wall. Consequently, the dynamics of the fluid transport along the pipe loop could not be considered, and the thermal mass of the fluid was not measured.
Finally, a 20-year long-term simulation of a one-story commercial hall building [102] showed that the consideration of seasonal thermal storage could be useful. The validation of a borehole model showed that the model performed in the whole-building simulation software IDA-ICE (see Section 9) could accurately simulate a dynamic performance, and it could be highly suitable for a coupling with a dynamic plant model.

7. The Effect of Groundwater Flow

In a lot of cases groundwater may be present within the geological regime in which a SGES is installed. This may significantly affect the performance of the heat exchanging process since an additional mode of heat transfer, due to convection, is introduced as well. In general, groundwater flow is positive with regard to the thermal performance of SGES due to a moderating effect on fluid temperatures in both heating and cooling modes [134]. However, a precise considering is rather complex and depends heavily on the type of the model used.
Indicatively, the effect of groundwater flow on the thermal performance of a BHE for GSHP systems was studied in the form of small-scale experiments and numerical simulations [135]. The authors concluded that it was advantageous to use a U-tube BHE with an oval cross-section in regions of potentially fast groundwater flow.
One of the serious shortcomings of line/cylinder source models, analytical g-functions, and other various superposition models (e.g., superposition borehole model—SBM, DST, etc.) is their conduction-based heat transport approach, while convection type heat transport is not considered. The latter is important in the case of groundwater flow existence. The velocity of groundwater can vary significantly. If soil exhibits low permeability and, thus, low groundwater flow velocity, the process of convection due to groundwater might be neglectable. In that case, additional heat transfer due to the presence of water in soil pores can be considered by using a higher thermal conductivity value [136].
In any case, groundwater flow requires coupling of the heat conduction equation and the heat advection equation [137]. The authors used various seepage velocities to show that groundwater flow influenced the average BHE temperature, and in the water-baring layer the average temperature was less as opposed to the dry regions. It was also noticeable that the temperature of the affected ground layer reached a steady state much sooner than in other regions. Additionally, when multiple boreholes were present in the direction of flow, it was observed that there was interference. This interference influenced the downstream borehole that could reach a higher steady-state temperature as indicated in Figure 4. Figure 4 was simulated by the authors using COMSOL Multiphysics with the transient heat equation for incompressible fluid and coordinate scaling. The geometry was divided into regions. Therefore, the groundwater velocity was applied to one of the regions using seepage velocity.
Sutton et al. [138] and Diao et al. [139] coupled an ILSM with a moving line heat source theory to overcome this drawback of analytical and semi-analytical models. Furthermore, the PILESIM software (see Section 9) that is based on the DST model, can enable the consideration of groundwater flow by simple approximations [140], whose accuracy has not been checked yet [5]. Similarly, Zhang et al. [141] adopted infinite and finite cylindrical source models [67] to account for groundwater advection.
Wang et al. [142] improved their model for higher groundwater velocities by comparing results with the results of the ANSYS CFX FEM software (see Section 9). The authors analyzed EPs performance accounting also for a groundwater advection. Their study revealed that groundwater advection can improve the heat exchange performance of EPs by five times (or higher), compared to the no groundwater seepage case.
Several FE models are capable of accounting for water advection heat transfer [143,144,145,146]. Their use requires long computation time, but enables the consideration of non-homogenous ground, complex geometry, uneven boundary conditions and thermal loading, and provides good insight of the results allowing for detailed parametric studies.
The Uppsala case study [113] allowed an analysis of the accuracy of different types/modules (246, 557a, 557b) of TRNSYS, comparing temperature difference between the fluid, groundwater and the borehole wall. The simulation results indicated that the type 557b, where the borehole resistance is known from experimental data and is a pre-set as an input is the most accurate of the types for groundwater filled boreholes. On the other hand, calculating borehole resistance based on the borehole thermal properties and geometry showed that type 557a seldom represents the reality in groundwater filled boreholes.
Gehlin and Hellstrom [147] examined three different models for the estimation of the heat transfer effect of groundwater flow and compared the results with the case of no groundwater flow. The three model simulations gave different temperature field patterns around the borehole resulting in lower borehole temperatures compared to the case of no groundwater flow. The authors also showed that even a relatively narrow fracture close to a borehole can lead to higher effective thermal conductivity.
Wang et al. [148] investigated the performance of a BHE under groundwater flow positioned in Baoding, China. The authors showed that the groundwater enhanced the thermal performance of the BHE, with the enhancement depending mainly on the thickness of the groundwater layer and its characteristics. In this specific case the BHE heat injection was enhanced on average by 9.8%, while the BHE heat extraction by 12.9%, compared with the case of no groundwater flow, for a total thickness of the water enclosing layer as a percentage of the borehole depth at 10.6%.

8. Thermo-Mechanical Interactions and Constitutive Modeling

There exist three general approaches for thermomechanical analysis of EP, namely empirical rules, load transfer models and full numerical models (FEM, etc.). Regarding numerical modeling, while there exist a number of published thermal response modeling case studies, as presented in Section 2, publications reporting the thermomechanical response of EGS are largely absent. Still, the thermomechanical analysis of EP has attracted certain attention [4,5,26,145,149,150] while detailed thermomechanical analysis of other EGS remains rare [50,151].
The empirical rules are based on assumptions of maximum changes in axial load and maximum deformation caused by heating/cooling cycles. They generally lead into overly conservative design [4]. On the other hand, load transfer models, which represent pile-soil interaction by load transfer rules, are much more accurate and widely used.
Knellwolf et al. [26] developed the Thermo-Pile software (see Section 9) [152], which uses a load transfer approach. The software was coupled with thermal response analysis and verified by back-analysis of two well documented experimental tests—real scale energy piles at Swiss Federal Institute of Technology in Lausanne [153] and Lambeth College in London [154]. The latter case study was used also for fully numerical thermo-hydro-mechanical analysis by Adinolfi et al. [150]. In that case study soil was modeled as bi-phase (solid-air, solid-water), isotropic and linear elastic-perfectly plastic material with Drucker-Prager yield surface. The pile was modeled as an isotropic thermo-elastic non-porous material. The numerical model was solved by FEM software COMSOL Multiphysics, with mechanical boundary conditions as pinned constraints at the bottom and rollers on the side of the pile, distributed load on the top of the pile and thermal boundary conditions as fixed temperature on the top of external boundary and adiabatic conditions on all other external boundaries. The thermal load was a function of time. Additionally, hydraulic boundary conditions were set regarding the real groundwater table. A reasonably good matching of numerical and experimental results was observed [150].
Similarly, the COMSOL Multiphysics package was used also for 3D numerical modeling of five driven steel piles in a pilot test in Hämeenlinna of southern Finland [145]. The test proved the importance of a precise determination of the mechanical behavior of energy piles, as high temperature fluctuations were recognized over very short distances at the depth below the tube curve, as well as in the first meter of pile length.
Other fully numerical applications exist in literature, for parametric studies of certain problems, using various software. For instance, the numerical simulation of a SGES in layered soil foundations using the Code_Bright software highlighted long-term settlements of a building induced by thermal cycling [151] (see Section 9). Also, the thermomechanical analysis of the response of embedded retaining walls with the Abaqus Standard FEM software (see Section 9) showed that the wall response was dominated by climatic temperature changes, while any changes due to heat exchange in EGS were negligible [50].
As seen in Section 2, while there exist a number of published thermal response modeling case studies (see also, [80,134]), publications reporting the thermomechanical response (see also Section 9) of EGS are largely absent.
In the perspective of the use of ground as a thermal reservoir, specifically to SGES, the study of soil thermo-hydro-mechanical behavior (THM) becomes particularly relevant. Due to the range of the temperature change in the ground resulting within a SGES, excluding the case of freezing, no phase changes nor chemical reactions are expected to occur. It is therefore considered that the main physical processes, which occur in the ground during SGES operation are thermal processes mutually coupled with hydraulic and mechanical ones.
THM behavior might have consequences both in the structural conditions and in the SGES energy efficiency, thus when considered, it can deepen the knowledge of ground behavior under such conditions and lead to a proper design.
The interactions between thermal and mechanical behaviors (in the temperature range of SGES) is essentially unidirectional; i.e., temperature loading induces mechanical effects in soil, while the influence of mechanical actions on the temperature field can be neglected. In turn, thermal and hydraulic effects are mutually coupled, i.e., the thermal loads can induce changes in pore pressures and in the water flow regime, while the hydraulic conditions can also affect the thermal field (as the pore fluids conduct and transport heat). Lastly, changes in effective stress, induced by pore pressure variations, can cause the mutual interaction of mechanical and hydraulic effects.
Numerous experimental results show the influence of temperature on the behavior of soils. Briefly, the most relevant features are as follows.
(i)
Response to heating-cooling cycles under isotropic conditions (constant effective stress), where all test results show a large influence of the soil stress history (by means of the over-consolidation ratio, OCR) on its volumetric behavior. Under normally consolidated conditions, heating induces thermoplastic contractive volumetric behavior, while heavily over-consolidated clays tend to exhibit a thermo-elastic response (dilation during heating and contraction during cooling). However, as noted by some authors (e.g., [155,156]) inelastic deformations are also observed for high OCR values in the cooling cycles. Moreover, in some tests the whole thermal cycle indicates the irreversibility of strain [157,158] (Figure 5a) Additionally, after the first thermal cycle a trend to elastic behavior was observed [159].
(ii)
Changes in pre-consolidation pressure induced by temperature, as well, where consistently, test results show a reduction of the pre-consolidation pressure p’c (interception of the yield surface with the mean effective stress) with an increase in temperature (in Figure 5b are shown the results obtained by Eriksson [160] on intact samples of Lulea clay).
(iii)
Changes in friction angle (critical state angle) with temperature, where the reduction of p’c with temperature is definite, but the trend of strength evolution is yet to be confirmed. Some researchers have concluded that heating caused a decrease in strength, while others have reported the opposite.
In summary, a suitable approach for soil thermo-mechanical behavior should encompass the elastic and irreversible volumetric strains of soil significantly affected by OCR values, from thermal compaction at low OCR values to thermal expansion followed by contraction at higher OCR, together with a change of the yield surface size (shrinkage during heating) and a potential change in soil critical state angle. The complexity increases if soil viscous behavior is considered. The effect of the cooling rate on the shape of the contraction curve during cooling was noted by Sultan [161]. More recently, some research studies proposed the consideration of this relevant feature of soil behavior in order to properly reproduce non-isothermal behavior (e.g., [162,163,164]).
To include all the above-mentioned features of thermal behavior requires the use of complex constitutive models. It is worth noting however, that the use of simpler stress-strain-temperature relations may be sufficient and, in many cases, the only available option, for that matter. Different classes of constitutive models for thermo-mechanical soil behavior, with different complexity levels are referred below, with more focus given on SGES applications. As regards the numerical applications of complex soil models in these geothermal systems, few studies exist in the literature.

8.1. Thermo-Elastic Models

The action of temperature on the ground constituents (solid minerals, water, and air) is generally assumed as linear elastic, under SGES temperatures exploitation as in the case of concrete or steel. Thereby a linear relation between the temperature rate and the volumetric strain rate dictated by the volumetric coefficient of thermal expansion is often used. Under restrained or partially restrained conditions temperature loading induces stresses in the soil. In thermo-elastic models, these stress–strain relations are assumed to be totally reversible.

8.1.1. Analytical Solutions

There exist in literature analytical solutions for assessing the isotropic thermal-only elastic (reversible) volumetric strain rate of a soil element in unrestrained (free expansion) and saturated conditions, such as [165,166].
ε ˙ v T = β T ˙
where β is the coefficient of volumetric expansion that can be either a scalar value for isotropic materials, or a tensor for anisotropic thermal expansion, and T ˙ the temperature rate.
These solutions give the volumetric thermal coefficient, which depends on the water and soil particles’ thermal expansion, the bulk modulus of the soil skeleton and the Biot’s modulus of the drainage conditions. Under the same conditions, the pore pressure rate can also be obtained analytically. The trend is that an increase in temperature results in the soil skeleton restraining the water expansion (as the water thermal expansion is higher than that of the soil grains), the pore pressure increasing, and consequently the effective stress decreasing. Reciprocally, a decrease in temperature results in a decrease in pore pressure (and an increase in the effective stress).
β = [ K + ( 1 n ) K f ] β g + K f ( n β w ζ ˙ T ˙ ) K + K f ;             u ˙ = K f K [ n ( β w β g ) ζ ˙ T ˙ ] K + K f T ˙
where K is the volumetric elastic tangent modulus, generally assumed temperature independent in the restricted range of temperatures, Kf the Biot modulus, βg the volumetric thermal expansion of the solid particles, βw the volumetric thermal expansion of water, ζ ˙ the rate of water flow in the soil voids, n the porosity, and u ˙ the pore pressure rate.

8.1.2. Elastic Stress–Strain-Temperature Relations

The elastic models that consider the mechanical effect of temperature are actually stress-strain elastic constitutive relationships in which a non-isothermal reversible strain rate tensor is added to the isothermal one. This tensor is responsible for the reversible thermal expansion and contraction of the material and is defined as
ε ˙ T = β 3 T ˙ I
where ε ˙ T is the strain-rate tensor. This results in the stress rate tensor σ ˙ :
σ ˙ = D : ( ε ˙ e ε ˙ T )
where σ ˙ is the stress rate tensor, D the tensor of elastic moduli, and ε ˙ e the elastic strain rate tensor. The mean effective stress rate p   ˙ is thus obtained for isotropic elasticity by
p ˙ = K ( ε ˙ v e ε ˙ T ) = K ( ε ˙ v e + β T ˙ )
where ε ˙ v e is the volumetric strain.
Different types of elastic models extended to non-isothermal conditions can be found in literature. Extensions for isotropic elasticity under non-uniform temperature fields can be categorized as: linear thermoelasticity, thermohypoelasticity and thermohyperelasticity.
Laloui and François [157] used non-linear thermo-elasticity (thermohypoelastic) for the elastic component of the ACMEG-T model; this was considered also by many other authors who used critical-state based models.
Hyperelastic models considering the effect of temperature (thermohyperelastic) can also be found in applications for soils. This type of models was used for example by Hueckel and Borsetto [167] and more recently by Cui et al. [155] who proposed the integration of the following expressions for the evaluation of the thermoelastic potential Ψ:
p = p 0 exp ( ε v e β 3 Δ T κ ( 1 + e 0 ) ) = Ψ ε v e ;             q = q 0 + 3 G ε s e = Ψ ε s e
where p 0 is a reference isotropic effective stress, κ the slope of the elastic compression line, q the deviatoric stress, q 0 the reference value of the deviatoric stress, G the elastic shear modulus, e0 the initial void ratio, ε v e the elastic volumetric strain, ε s e the elastic shear strains strain, and ΔT the temperature difference.
Given the limited temperature range involved in the exploitation of SGES, as well as its amplitude, temperature independent thermal expansion coefficients are generally considered. There are however alternative propositions in literature, such as the one of Laloui et al. [158] (for more generalized applications including heat storage, geothermal structures, injection and production activities, petroleum drilling, zones around buried high-voltage cables, high-level nuclear waste disposal) according to which a thermal expansion coefficient of the soil skeleton βs depending on temperature is given as
β s = ( β s 0 + ζ T ) ξ
where βs0 is the isotropic thermal expansion coefficient at a reference temperature and ξ is the ratio between the preconsolidation pressure pcr0 and the effective mean pressure p at ambient temperature T. The proposed expression was used to explain observed high volumetric expansion deformation for high over-consolidated soils (see above), however there is a lack of experimental values to support their reversible nature.

8.2. Thermo-Elastoplastic Models

In thermo-plasticity, irreversible deformations can be induced by temperature. As mentioned above, there is systematic experimental evidence showing the occurrence of plastic strains caused by thermal loading, namely a significant contraction for normally consolidated clays when heated, with a significant part of this deformation being irreversible upon cooling. Also, in some cases, for highly over-consolidated soils, during heating important dilations are observed, not totally recovered during cooling.
For the numerical reproduction of this thermal hardening behavior, largely influenced by the soil stress history, adequate constitutive relations are required. A brief compilation of thermo-plastic -based formulation models is presented. These constitutive stress-strain-temperature relationships have different degrees of complexity, requiring different numbers of parameters. The applications of complex constitutive models to analyze the performance of SGES systems is still relatively limited.

8.2.1. Analytical Solutions

In the same manner as for the thermo-elastic behavior, analytical solutions can be obtained for the free expansion of a soil element under thermal loading, to assess its thermal expansion volumetric coefficient and the rate of the induced pore pressure.

8.2.2. Thermo-Elastoplastic Behavior with Critical-State-Based Structure (Extension of Modified Cam-Clay Model)

A large part of the constitutive thermo-elastoplastic models proposed in literature have a critical-state-based structure. Here are referred the simpler models that involve relatively small changes in the modified Cam-Clay model, in order to obtain thermal permanent deformations. These models take into account the dependence of the yield surface on temperature, which according to Hueckel and Borsetto [167] decreases as temperature increases.
Constitutive equations consisting of an extension of the critical-state model to non-isothermal conditions were firstly proposed in [167], based on the observations of tests over a wide range of temperatures (15 °C to 115 °C). The incremental relations of the critical-state model consisting on an elastic law, a plastic flow rule, a hardening law, and a yield condition, were generalized to explicitly depend on temperature. The basic concept of these model extension was the elastic domain dependence on temperature; it was assumed to shrink during heating (thermal softening) and to expand during cooling. At constant stress, thermal softening may be compensated by thermoplastic strain hardening.
The elastic law was generalized to thermal conditions by introducing a reversible thermal isotropic strain. For plastic behavior, the hypothesis of thermoplastic non-associativity was considered and the loading and unloading criteria were defined.
In the light of a series of experimental tests in three different clay soils, the thermoplastic model was analyzed and tested in [168]. For example, an evolution law of the pre-consolidation pressure with temperature was proposed, as follows.
p c = p c 0 exp ( 1 λ K T { e 1 ( 1 a 0 Δ T ) [ e g ( 1 + e 0 ) ε v T p ] } ) + 2 ( a 1 Δ T + a 2 Δ T | Δ T | )
where
K T = [ K 1 + e 0 + ( α 1 + α s Δ T ) Δ T ] ( 1 + e 0 )
and λ is the compressibility of the normal compression, a0, a1 and a2 are constant coefficients (usually negative, corresponding to the reduction of the semi-axis of the yield surface due to temperature alone), e0 and e1 values of the void ratio at the initial state at a hypothetical state corresponding to pc = pc0, T = T0 respectively, eg the void ratio at the maximum pre-consolidation pressure, ε v T p the thermoplastic volumetric strain, and a1 and a3 coefficients defined in Hueckel and Borsetto [168]).
A reasonably prediction of the clay behavior in such processes as thermoplastic consolidation and other aspects, like thermal pressurization and thermal ductilization of clay in triaxial compression, was attained [168].
On this conceptual basis, which combines the critical state theory with thermoplasticity directed to soils under thermomechanical conditions, several numerical models were subsequently proposed.
Robinet et al. [169] proposed a similar model considering two separate yield mechanisms. Abuel-Naga [170] proposed a model with extensions to thermoplastic behavior under isotropic conditions and subsequently to triaxial stress space. The model required merely a few constants more than those of the modified Cam-clay model. It was characterized by a non-associative temperature-dependent flow rule. The yield and potential surfaces evolve according to a combined thermal distortional and rotational kinematic rule additionally to the isotropic hardening rule. This extension to the triaxial stress space in the strain-hardening zone of the deviatoric stress plane (q—p) was based on assuming a temperature dependency of the shape of the yield surface attributed to the thermally induced soil fabric changes. The mathematical yield surface expression introduced in [171] was adopted, enabling a yield surface shape flexibility using only one more parameter than the conventional Cam-Clay model parameters.
Vieira and Maranha [166] also proposed a critical-state soil model with thermal hardening, with one yield surface that was calibrated with laboratory results. In this model the thermal effects were also introduced at two levels: reversible thermal volumetric strains (thermo-elasticity) inside the yield surface and thermal hardening (evolution of the yield surface with temperature). The model included also a change of the yield function shape on the super-critical region and the dependence on the Lode’s angle.
The pre-consolidation pressure (size of the yield surface) temperature-induced change is modeled using only one constant. An increase in temperature leads to a reduction in the size of the yield surface, with the opposite happening when temperature decreases. The evolution of the rate of the pre-consolidation pressure is obtained by
p ˙ c = p c ( υ λ κ ε ˙ v p d T T ˙ )
where dT is the constant that controls the shrinkage of the yield surface, v the specific volume, and ε ˙ v p the plastic volumetric strain rate.
Two types of isotropic hardening can thus occur: one caused by volumetric plastic strains and one by temperature variation. The model was applied to a floating pile thermoactive pile in a normally consolidated clay revealing the importance of considering thermoplasticity in SGES systems.
While introducing relevant features of soil thermal behavior, such is the case of thermal contraction of normally consolidated soils due to the shrinkage of the yield surface resulting from temperature increase, that these constitutive relationships have important limitations: namely their incapability of reproducing volumetric plastic strains for higher OCR values (under isotropic strains), cyclic thermal behavior and absence of inelastic strains inside the yield surface.
More complex formulations involving partial saturation were proposed in [172], and destructuring in [173]. However, these aspects will not be dealt with here.

8.2.3. Thermo-Elastoplastic Models with Multi-Mechanisms (Isotropic and Deviatoric Thermoplastic Strains)

The model proposed by Laloui et al. [158] includes a combination of two irreversible processes: one thermo-mechanical isotropic mechanism and one deviatoric mechanism. A new version of this model, termed ACMEG-T (Advanced Constitutive Model for Environmental Geomechanics—Thermal effect), was extended to bounding surface theory [157].
The model is based on the model of Hujeux [174], derived from the theory of multi-mechanisms plasticity [175]. Each mechanism is activated if the stress state reaches the yield function of a specific mechanism. Each dissipative process is described through an evolution law activated by a yield function, a dissipative potential, and a plastic multiplier. The thermo-plastic strain rate tensor is obtained by
ε ˙ T p = k = 1 2 λ k g k σ
where gk are the plastic potentials corresponding to each mechanism and λk the respective plastic multipliers.
The isotropic and the deviatoric yield functions define a closed domain in the effective stress space, where the soil behavior in reversible. Conversely, the mechanisms are activated if the stress state reaches the corresponding yield function, fiso and fdev for the isotropic and deviatoric mechanisms respectively. In that case strains are produced.

Isotropic Thermo-Plastic Mechanism

The yield limit fiso of the isotropic thermoplastic mechanism represented in the p′-T plane is expressed by
f iso = p p c r iso
where riso is a parameter related to the degree of plastification (mobilized hardening).

Deviatoric Thermo-Plastic Mechanism

The deviatoric yield limit, an extension of the original Cam-Clay model, also takes into account the effect of temperature as follows.
f dev = q M p ( 1 b ln p d p c ) r dev
where b is a material parameter defining the shape of the deviatoric yield limit, d the ratio between the pre-consolidation pressure and the critical pressure, M the slope of the critical state line in the p-q space and is temperature dependent, and rdev the parameter related to the degree of plastification (deviatoric surface).
The isotropic and the deviatoric yield limits are coupled through the hardening variable, ε v p , to establish a relation between the plastic multipliers. The most recent version of the model has 16 constants. A schematic representation of the yield surfaces is shown in Figure 6a, while Figure 6b presents numerical simulations of heating-cooling cycles at different OCR values.

8.3. Other Thermoplastic Models

In relation to nuclear waste management programs Cui et al. [155] proposed an elastoplastic model for non-isothermal conditions, with particular attention given to the volume-change behavior and the OCR effects, based on experimental data. The proposed model enables the prediction of thermo-plastic strains, compressive strains, as an inverse relation between pre-consolidation pressure and temperature. Proposed is a new plastic mechanism, based on a yield locus called TY, associated with the LY yield curve, where both curves constitute the limit of the elastic zone in a T–p′ plane. This mechanism allows the generation of thermal plastic strains for over-consolidated soils under heating.
Some numerical simulations carried out with this model are shown in Figure 7. As it can be noted, compressive thermal plastic strains are obtained for different OCR values. However, this model does not enable the occurrence of thermo-plastic expansion strains for over-consolidated soils.

8.4. Thermo-Viscoelastoplastic Models

Rate dependency is one of the most relevant features of soil behavior, and the possibility that rate-dependent behavior influences the thermo-mechanical response of soils was used by some authors to overcome limitations of the thermo-elastoplastic models. Experimentally, the effects of the loading rate on the shape of the contraction curve were investigated by Sultan [161] on boom clay. The tests results showed a significant effect of the cooling rate on the shape of the curve, with thermal expansion decreasing as the unloading (cooling) rate increased. Cui et al. [155] drew attention to the likelihood of some experimental results in literature, showing an expansion during cooling, which should be considered with care, particularly in cases where fast cooling was performed.
Also, some thermo-viscoplastic models for soils can be found in literature. One of the first was proposed in [178]. In this Perzyna-type viscoplastic model, also based on the multi-mechanisms’ theory, the following expressions were proposed for evaluating the irreversible strains, ε ˙ v T v p and ε ˙ d T p , (volumetric and distortional) of each mechanism k:
( ε ˙ v T v p ) k = f k Ξ f 0 n μ μ . ( Ψ v ) k ;             Ξ = m o r c ;             ( ε ˙ d T p ) k = f k Ξ f 0 n μ μ . ( Ψ d ) k
where μ is a measure of viscosity which may vary with temperature, f k Ξ is the yield function of each mechanism, m parameter related to deviatoric yield surface, c parameter related to cyclic yield surface, f0 a reference stress, Ψ v isotropic plastic flow of each mechanism, and Ψ d deviatoric plastic flow of each mechanism.
Some of the results obtained by this model under different stress paths are shown in Figure 8.
A semi-empirical elastic-thermoviscoplastic model for clay was proposed recently by Kurz et al. [162]. The plastic component includes viscous strains defined by a creep rate coefficient, which varies with plasticity index and temperature. The model formulation is based on the identification of a single material parameter, ψ a creep rate coefficient that depends on the mineralogy of the clay and temperature obtained from simple vertical compression tests. This relation is based on the similarities observed between isothermal tests at different strain rates and tests at different temperatures at the same strain rate [179] (Figure 9a).
The temperature dependent viscosity function adopted in the model is based on the proposal of Fox and Edil [180] based on the assumption that viscosity decreases as temperature increases, meaning that a clay at higher temperature will experience higher creep rates and more creep straining than one at lower temperature T.
Notwithstanding having a simple structure and a lack of some physical support, the model was able to reproduce some features of soil thermomechanical behavior as can be seen in Figure 9b, namely the expansion (dilatant viscoplastic strains) for highly consolidated soils. However, the trend of the curves is not well reproduced for low OCR.
A thermo-elasto-viscoplastic model, directed to the use of soft rocks as deep geological disposals for high level radioactive waste, was proposed in [181]. The model, modified from a previous version, includes the influence of intermediate stress. The model is based on the subloading concept [182]. The performance of the modified model was confirmed with drained triaxial compression tests and creep tests on soft sedimentary rocks and artificial soft rocks under different temperatures. Nine parameters are contained in this modified model, which is equivalent to the one proposed by Zhang et al. [183], with the only difference being that, in order to describe the thermal behavior of geomaterials, the thermal expansion coefficient, considered constant, is added to the model. The model does not include explicitly a shrinkage of the yield surface with temperature.
The subloading concept due to Hashiguchi [182] was also used to reproduce the influence of non-isothermal conditions on the stress–strain-time behavior of soils. A version of the viscoplastic subloading soil model proposed in Maranha et al. [184], restricted to isotropic stress and strain conditions, was extended to non-isothermal conditions in [184], introducing temperature dependence of the size of the yield surface and of the viscosity. This model was able to capture well the large expansion volumetric strains observed in highly over-consolidated clays in the initial stages on heating tests [185]. The strains were considered irreversible as the corresponding coefficient of thermal expansion is variable and much larger than those of the constituent minerals of the soil’s solid phase. In order to be able to reproduce these tests, it was necessary to assume that the previous unloading (to obtain the over-consolidated state) was fast enough so that creep was still occurring at the beginning of the heating stage. Using this model, three (3) heating only tests from Cekerevac and Laloui [185] could be precisely reproduced (Figure 10a). Enhanced versions of this model, introducing a temperature dependence motion to the center of homothety were presented in [164]. This modification was necessary for the reproduction of the large over-consolidated irreversible thermal expansions without assuming that creep strains were still occurring during the heating test. Using this model, a very good adjustment of the three heating and cooling tests described in [186], considering the variability between the three soil samples and the same material model constants were used for all tests (Figure 10b).

9. Software Tools

There exist different kind of software tools with implemented various models to predict the response of SGES with an acceptable accuracy. Most of these are for thermal studies. The ones suitable for THM applications are indicated in bold face. The basic general characteristics of the most common tools, including type of a model and referencing applications, are shown in the Table 1.
More detailed information of software listed in Table 1 is presented below, with attention given to physical models’ software. Detailed information about physical models of specific software programs are very difficult, often impossible (even a sensitive issue, as commercialism is involved), to find and compare. Therefore, the discussion given below is mostly informative.
EED computes vertical orientations (boreholes) only. A single instruction multiple data (SIMD) procedure is used to compute borehole fluid temperatures. The relevant g-functions used are dependent on the spacing between the boreholes and their depth, as well as on the tilt angle in the case of graded boreholes. The key ground parameters (such as thermal conductivity and specific heat), but also the properties of the heat carrier fluids and the pipe materials are provided by databases. The hourly heating and cooling or the monthly average loads are used as the input data, while an extra pulse for peak heating and cooling loads of several hours is possibly considered at the end of every month. The borehole thermal resistance is computed using the borehole geometry, the pipe material and geometry, and the grouting material.
EnergyPlus uses the model developed by Marcotte and Pasquier [191], a discretized LSM with boreholes discretized into segments. Each segment’s temperature response on all other segments determines the response factor for the selected geometry. The surface effects are estimated by the creation of “imaginary” boreholes that are mirrored about the ground surface.
The short time-step response is computed using the model of Xu and Spitler [192]. The model maps the 2D borehole geometry onto a 1D radial geometry, preserving the thermal mass, including the fluid, of the system. The multipole method of Claesson and Helström [193] is used to correct the pipe and grout conductivities such that the (correct) borehole resistance is preserved. The temperature response of the 1D domain is calculated through an FVM. The borehole wall temperature is then used to compute the short time-step g-function. EnergyPlus uses a load aggregation scheme too.
In EWS the ground can be vertically divided in at most 10 layers. The heat equation in cylindrical coordinates is solved at every layer, allowing the calculation of the common case of a ground consisting of several layers with different properties. The simulation of the ground temperatures near the boreholes (1.5–3 m) are performed by the Crank–Nicholson method, where the average fluid temperature of each layer acts as an inner boundary condition. An explicit time step procedure is used for the simulation of the unsteady fluid. Thus, the start-up borehole behavior can be calculated. The outer boundary conditions are obtained through the dimensionless thermal response factor (see g-functions). Note that one can choose between the methods of Eskilson [52] or of Carslaw and Jaeger [48].
In EWS it is possible to simulate a variable mass flow rate. Hence, the user has two options in relation to thermal borehole resistance: (i) either the thermal borehole resistance is recalculated for each calculation step (through the heat transfer coefficient) or (ii) it is fixed for the whole simulation. The issue of the non-constant heat extraction rate and the regeneration of the ground is treated through the superposition of several constant heat extraction rates that start at different times. The approach selected allows for the use of different time steps, that is to say the shortest time step is used for the unsteady calculation of the fluid, while a larger time step is used for the calculation in the simulation area by Crank–Nicholson. Note that a time step of one week suffices for the g-functions calculation of the ground outside the simulation area. The use of different time steps is needed to account for the temperature disturbances always coming from the boreholes (hence, the smallest time step is needed close to the boreholes), while only averaged heat extractions or inputs are observed far away from the boreholes.
FEFLOW simulates heat transfer, mass transfer and groundwater flow in fractured or porous media both in saturated and unsaturated conditions. FEFLOW can deal with: shallow and deep geothermal installations, open and closed loop systems, BHEs, aquifer thermal energy storage, GHE arrays, interaction with heating and cooling installations, advection-conduction/dispersion heat transport, free, forced and mixed convection, thermohaline convection, coupled density dependent simulation for varying temperature, linear or nonlinear temperature-density relationships, predefined or user defined temperature-viscosity relationships.
In FEFLOW a vertical closed loop can be modeled in two different approaches [196,197,198]: (i) following a built-in module, inserting a simplified 1D element at the center node of the BHE and coupling it with the rest of the model domain; (ii) discretizing, through a fully discretized 3D model (FD3DM), all borehole elements and assigning on a nodal/element basis flow and thermal material properties. The discretization approach requires an increased computational time and an increased amount of resources needed. This drawback is balanced by a more detailed temperature distribution, in and around the borehole, and an increased precision, especially near steady-state conditions [196].
FEHM simulates flows in large and geologically complex basins, but also complex coupled subsurface processes. Its features allow accurate representation of wellbores. FEHM is used to simulate groundwater flow in shallow or deep and fractured or un-fractured porous media. Subsurface physics ranges from single fluid/single-phase fluid flow, for simulations of basin scale groundwater aquifers, to complex multifluid/multi-phase fluid flow that includes phase change with boiling and condensing. While the FE method is used to obtain more accurate stress calculations, one of numerical methods used in FEHM is also the control volume (CV) method for fluid flow and heat transfer equations; this allows enforcing energy/mass conservation. A finite difference (FD) scheme is also available. FEHM uses analytic derivatives in the Newton–Raphson iterations.
GLD supports the following heat transfer theories: the ASHRAE standard, the European standard and hourly simulation function. It includes two models. The first one is based on the CSM and allows for fast length and/or temperature calculations. The second one is based on the LSM and can generate detailed monthly and/or hourly temperature profiles in time, for monthly loads/peak data and/or hourly loads data respectively. This second model can also model the impact of balanced or unbalanced loads.
In the GLD Borehole Module the first model uses the vertical borehole equations based on the heat transfer solution from Carslaw and Jaeger [48] and Ingersoll [66]. Additionally, the model adopts the adjustments of Kavanaugh and Deerman [201] on the methods of Ingersoll to account for hourly heat variations and a U-tube arrangement. The model employs the borehole resistance calculation techniques by Paul and Remund [202]. GLD’s second model can calculate the borehole wall temperature with respect to time for a constant heat flow rate extracted from the borehole.
The GLD Horizontal Module effectively combines the CSM of Carslaw and Jaeger [48] and the multiple pipe method by Parker et al. [200]. As in the case of the Borehole Module, adjustments suggested by Kavanaugh and Deerman [201] on the equations are adopted.
To determine the length of pipe necessary for different surface water systems, Kavanaugh and Rafferty [71] conducted experiments for different-size pipes in coiled and “slinky” configurations in both heating and cooling modes. GLD uses a polynomial fit of the above-mentioned experimental data when calculating surface water systems.
GLHE-Pro uses the multipole method [193] with 10th-order multipole solution, to compute the local borehole thermal resistance (fluid to borehole wall) and the internal thermal resistance (fluid to fluid) needed for calculating short-circuiting effects arising from the use of analytical expressions by Hellström [63]. The convective resistance is held constant.
Monthly peak and average entering fluid temperatures from the borehole(s) or horizontal Ground-Loop Heat Exchanger (GLHE) to the HP can be determined from hourly simulations based on average monthly or hourly loads. Boreholes can also be placed within irregular spacing in the horizontal plane, and they can be inclined from vertical. For the last design, the Free Placement Finite Line Source (FPFLS) model is used. For the standard vertical borehole model, the so-called Standing Column Well (SCW) model [204] acts as a configuration option. The SCW model uses an hourly numerical simulation based on mass and energy balances, rather than using g-functions. The modeling of any short circuiting (see above) between the annulus and the inside of the dip tube is analytical. Also, any well draw-down is neglected, as it is ‘small’ in most standing-column-well locations.
GshpCalc, formerly known as GchpCalc (version 4.2 and earlier), is a free software developed and provided by the University of Alabama. It includes vertical GHEs as well as groundwater heat pumps and surface water heat pumps. The software is based on the method suggested by Kavanaugh [69] with the use of an ASHRAE standard. The software supports single, double, and triple U-tube with series connected or parallel. It can also accommodate the design of hybrid GSHP with cooling tower and the water heating with heat pumps. It can use as an input heating and cooling load data from the TideLoad10v1 (or later versions) software. A comparison between five GHE design software programs [206], indicated that the GchpCalc provided the most accurate results [205]. It does however have its limitations with the most notable one regarding the maximum entering water temperature [207].
IDA–ICE is a 3D model combining vertical or leaning boreholes of equal length. The model uses the superposition of the 1D vertical field for the undisturbed ground temperature (including geothermal temperature) and the cylindrical 2D fields around each borehole to calculate the interactions between holes, the temperatures and the pressure drop in the brine liquid circuit. There are the following restrictions in the model: (i) the ground is one layer; (ii) all holes are of the same length; (iii) the GHEs are of the U-tube type; (iv) the borehole resistance is constant. Numerical errors can be eliminated allowing one to see how the equations truly behave with a time resolution of seconds, by the right choice of tolerance parameters.
In OpenGeoSys, the modeling of thermo-hydro-mechanical-chemical (THMC) processes of a vertical closed loop is possible through the approaches of Al-Khoury et al. [115], Diersch et al. [197], Diersch et al. [198] and Diersch [194]. They all suggested the separate inversion of the matrix system for the BHE domain, followed by the integration of its influence into the soil domain using a Schur complement [183]. The non-linearity of the governing equations cannot be eliminated without an iteration step. Piccard iterations are used [209], whereby convergence is reached after one iteration in nearly all simulations. Due to a sudden temperature change at the beginning of the simulations, a small-time step size in terms of minutes is proposed, otherwise many Piccard iterations will be needed for a time step. The heat transport process together with BHEs is simulated through a dual-continuum approach that is adopted to treat the soil and BHEs parts separately. Prism elements are used to discretize the 3D domain for the soil part, while chosen 1D line elements along the edge of the prism elements form the second domain, representing the BHE.
TOUGH3 deals with multiphase, multicomponent, and multidimensional systems, solving mass and energy balance equations for fluid and heat flow. It considers the transport of latent and sensible heat in conjunction to the movement of gaseous, aqueous, and non-aqueous phases, and the transition of components between the available phases. In each phase, advective fluid flow occurs under gravity pressure and viscous forces following the multiphase extension of Darcy’s law, which considers capillary pressure effects and relative permeability. In each phase, diffusive mass transport can occur too. Also included are Klinkenberg effects in the gas phase, and vapor pressure lowering due to capillary and phase adsorption effects. Heat flow occurs by conduction, convection, and radiative heat transfer, considering local thermal equilibria of all phases. TOUGH3 can simulate the injection/production of fluids/heat and can consider wellbore flow effects. For fractured media, implemented are methods of double-porosity, dual-permeability and multiple interacting continua (MINC) [212,213,214].
For systems of regular grid blocks the integral FD method is completely equivalent to conventional FD. In order to avoid impractical time-step limitations in flow problems with phase (dis-)appearances, the implicit time-stepping with the 100% upstream weighting of flux terms at interfaces are necessary; unconditional stability is hence achieved [211]. Newton-Raphson iterations for each time step are used to solve the resulting strongly coupled, nonlinear simultaneous algebraic equations.
TRNSYS is used to simulate photovoltaic, solar domestic hot water systems and thermal performance of buildings. One of the main advantages of TRNSYS is the ability of simultaneously transient calculation methods, so that the entire system can be simulated by connecting different components. In TRNSYS several types/modules for modeling boreholes as a component are included: 203 FLSM, 203 ILSM, 203 CSM, 243 CSM, 243 ILS, 243 CSM, 244 ILSM, 244 CSM, 246, 451a, 557a (no capacity), 557b (with capacity), 557c (capacity), 557d (no capacity).
Models of type 557 are most often used. They include the Duct Heat Storage (DST) model. In DST the ground temperature is a result of the superposition of three solutions: a global temperature, a local solution, and a steady-flux solution. An explicit FD method using 2D axial-symmetric formulations solves the global and local problems. The steady-flux solution for the storage volume is obtained analytically with pre-calculated g-functions.
A load aggregation scheme is used as a method to make long-term simulations more efficient; it is described in detail by Bernier [215]. The main principle behind the method is weighting the impact of the past load history on the current heat transfer. The distant history is aggregated in big (time) blocks and the recent history in smaller (time) blocks. The aggregated ground loads are then updated at each time step.
As already mentioned in previous Sections, the Uppsala case study [113], gave an evidence, that types 557a and 557b of DST model are considerable quicker (about 6 times) than 246 for long term simulations with the same time step in all cases. All 3 modules 246, 557a, and 557b seem to be suitable for installations with a balanced heating and cooling load over the year as confirmed in the Uppsala case study. The difference between average numerical and experimental fluid temperature in these cases is less than 1.1 °C.
Thermo-pile calculates temperature variations’ influences on stresses and strains with pile foundation. The basic assumptions are as follows: (i) the pile displacement calculation is done using a 1D FD scheme (with only the axial displacements considered); (ii) the pile properties, namely the diameter, the Young modulus and the coefficient of thermal expansion, are constant along the pile and with respect to temperature; (iii) the soil and soil-pile interaction properties are not affected by temperature; (iv) the following relationships are know: between friction/shaft displacement, between head stress/head displacement, between base stress/base displacement.
Moreover, the thermomechanical response of the pile is assumed to be thermoelastic and time independent. The load-transfer method is used to model the soil response. This is done using load-transfer curves to represent the soil response, with these curves linking the displacements of the pile to the mobilized bearing capacities.
The limitations of the Thermo-Pile are as follows: (i) safety factors are not included in the calculation process; (ii) only one pile with circular section is considered; (iii) tensile mechanical solicitation and bending moments are not considered; (iv) negative friction is not considered; (v) geothermal energy is not considered.
More software programs can be found in literature, with some of these falling back or made redundant in recent years, as bigger companies and open source software have been taking over. Such programs include ECA, GAEA, GS2000, HVACSIM+, PILESIM2, and Right-Loop. ECA is a commercial software, developed by the Elite Software Development Inc. The methodology used is described in the ASHRAE design and data manual for closed-loop ground heat pump systems. Vertical and horizontal pipes can be modeled with built-in heat pump performance data and weather data, but the user can input data manually. GAEA is a basic UI software developed according to models by Gnielinski [217], Grober et al. [218], Krischer [219] and the analytical model of Albers [220]. The GAEA model is described in [222] and can calculate only basic horizontal GHEs. GAEA has failed on updates for today’s standards and no updated or future versions are indicated. GS2000 was developed by Caneda Research, Canada, and was free to use. It was based on the LSM and CSM and could coop with both vertical and horizontal GHEs. Limitations of the software were the pre-determined limits, with no option of manual input by the user.
A comparison between EED, GLHE-Pro, and GS2000 for different scenarios, evaluated in TRNSYS, was presented in [230]. The author discussed on the different borehole lengths suggested by the different software and concluded that although EED, GLHE-Pro, and GS2000 provided undersized boreholes compared to TRNSYS simulations. Possible reasons for that are: (i) poor interpretation of the load data by the software or (ii) assumptions made about thermal mass of the borehole and groundwater movement.
Another software worth mentioning is the PILESIM2, developed initially by the Laboratory of Energy Systems (LASEN) [9] and the later version by Pahud [140], which is based on the simulation tools of TRNSYS and adapted to TRNSED format. The simulation models used in the software are the TRNVDST [226,231] and TRNSBM [225]. Although not “modern”, the software was used in relatively recent papers (e.g., [224]).
Furthermore, there are software developed by manufacturing companies or installation companies, either for their own line of production or for other external purposed; such example is the GeoLink Design Studio software developed by Water Furnace Intl, USA. Other software reported in older papers [207] are WFEA (Water Furnace International, Fort Wayne, IN USA), GL-Source (Kansas Electricity Utility, Topeka, KS, USA), GEOCALC (HVACR Programs, Ferris State Univ., Big Rapids, MI, USA) and GLGS (Oklahoma University).
It should also be noted that there exist web-based software tools, such as the notable case of LoopLink, developed by Geo-Connections Inc. Such options look very promising as they take advantage of the clouds solutions and are continuously updated.
There are cases, especially for research purposes, that a short-term analysis or detailed investigation and analysis on the GHEs is required. It could be for a new geometry, for a 2D or 3D study, or any other aspect that could influence the GHE. These cases require a CFD approach using software that are capable and suitable for such investigations. There are several commercially available software programs, such as COMSOL Multiphysics (COMSOL Inc.), FLUENT (ANSYS), STAR-CCM+, Autodesk CFD (Autodesk), Abaqus FEA (Dassault Systèmes), Altair AcuSolve (Altair Engineering) and SimScale (SimScale). Also, there are open source software, such as OpenFOAM (The OpenFOAM Foundation), SU2 (Stanford University Unstructured Project), Elmer (CSC—IT Centre for Science), and Advance Simulation Library (ASL by Avtech Scientific).
Another approach would be to solve specific equations, such as the LSM or the CSM, using mathematical based software. Minimum skills of programming would be required for the user of software that have friendly UI. Such software programs include MATLAB (MathWorks) or the open source alternative GNU Octave. Most of the CFD software mentioned above are also capable of performing such simulations.

10. Conclusions and Discussion

The current paper has aimed to give an overview of the modeling aspects of SGES and eventually contribute to the design of systematic guidelines. To this end, the main analytical and numerical methods used in literature for the investigation of the thermal behavior of SGES have been presented. It must be said that each method may be suitable for specific problems and exhibits its pros and cons. A comparison between the various analytical and numerical models may constitute an important future research task.
Auxiliary to the chosen model, are various factors that can affect the model design and the method used. Boundary conditions may come in the form of temperature conditions or heat flows on boundaries depending on the dimensions of the model chosen.
In their turn, spatial dimensions may be crucial. For instance, most analytical solutions assume an infinite length for the GHE to simplify the calculations, while for shorter GHEs, such as EP, the end-effects must be considered in the long term.
Regarding SGES parameters, it may be necessary to use some processes for their determination to overcome practical difficulties. Such processes include borehole and model excitation, parameter estimation techniques, evolutionary and stochastic search algorithms, sequential forward-evaluation, and sequential backward-evaluation.
Short-term and long-term analyses may be both necessary to perform when dealing with different regimes of SGES. For example, due to scale differences, a TRT should be carried out much longer for an EP than for a BHE.
The presence or not of groundwater may be crucial in the modeling of a SGES and may require the coupling of the heat conduction equation with the heat advection equation. In general, groundwater flow is beneficial to the thermal performance of SGES, but a precise considering is rather complex and depends heavily on the type of the model used.
The thermo-mechanical behavior of soils may assume relevance in SGES for the case of TAS systems, which exhibit the double functioning of GHEs and structural elements. Regarding thermo-mechanical interactions and constitutive modeling, a brief overview of the main features of soil stress-strain behavior under non-isothermal conditions, mostly directed to the use of SGES, has been presented. A remarkable influence of the soil stress history on its thermo-mechanical behavior and a significant irreversibility under thermal actions are observed. Namely, for the case of normally consolidated soils under constant stress conditions, an increase in the water pressure can be induced and consequently an effective stress drop can occur. Moreover, the rate dependent effects seem to play an important role. In order to reproduce all these features’ behavior, increasingly complex models from thermal elasticity to subloading or bounding surface thermo-viscoplasiticity have been proposed in the literature enabling the reproduction, with increasing accuracy, of the details of the loading mechanical and thermal sequence.
Complex constitutive models face the problem of the large number of constants required for calibration, and the scarcity of available tests for obtaining these constants. Additionally, issues related to the difficulties in obtaining these parameters, which interact in non-trivial ways, require proper algorithms, such as genetic algorithms.
Despite the difficulties in obtaining robust numerical models based on complex constitutive relations for soils under thermal actions (some of which have been mentioned), this is the way to accurately reproduce soil THM behavior and to have a proper knowledge of the behavior of SGE systems, even if an “adequate” design can be obtained with simpler models.
Finally, very useful tools for SGES modeling practices are the available commercial or open source developed software programs, web-based or not. These include software specifically built for geothermal applications, other for more general renewable energy sources applications, and other for more general purposes. It may even be possible or necessary to tackle certain SGES problems through a CFD approach using software that are capable and suitable for such investigations.
The end-result is that the nature of the SGES application and the required precision will guide the investigator/researcher/engineer as to which model/method/software they should choose.
This could be advantageous for future development of incorporating all types of SGES (see Introduction and Figure 1) in one software package. Such a package could include a wide range of mathematical models with the assistance of machine learning techniques to select the most appropriate model per study case. Any existing official guidelines and public restrictions and regulations could also be adopted within the software package. For effectiveness, the software could be offered as a regularly updated web-based platform to be used by engineers.

Author Contributions

Conceptualization: P.C., A.V.; Coordination: P.C.; Writing—original draft preparation: All authors; Writing-review and editing: P.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The authors gratefully acknowledge Transport and Urban Development COST Action TU1405—European Network for Shallow Geothermal Energy Applications in Buildings and Infrastructures (GABI; www.foundationgeotherm.org).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

ACMEG-TAdvanced Constitutive Model for Environmental Geomechanics—Thermal effect
BHEBorehole Heat Exchanger
CFDComputational Fluid Dynamics
CSMCylindrical-Source Model
DSTDuct ground heat Storage model
EGSEnergy Geo-structures
EPEnergy Piles
FDFinite Difference
FDMFinite Difference Model
FEMFinite Element Model
FLSMFinite Line-Source Model
FSFactor of Safety
FVMFinite Volume Model
GHEGround Heat Exchanger
GPMGeothermal Properties Measurements
GSHPGround-Source Heat Pumps
HPHeat Pump
ILSMInfinite Line-Source Model
LSMLine-Source Model
OCROver-Consolidation Ratio
SBMSuperposition Borehole Model
SCWStanding-Column Well
SGEShallow geothermal energy
SGESShallow Geothermal Energy Systems
TASThermo-Active Structure
THMThermo-Hydro-Mechanical
TRTThermal Response Test

Nomenclature

Ei(.)Exponential integral
GElastic shear modulus
J0,1Bessel functions of the first kind
KVolumetric elastic tangent modulus
KfBiot modulus
LBorehole length [m]
MSlope of the critical state line
Q ˙ Heat injection rate [W m−1]
RbThermal resistivity [K m W−1]
TTemperature [K]
T0Undistributed ground temperature or the initial borehole temperature [K]
TbTemperature at the borehole radius [K]
Y0,1Bessel functions of the second kind
cpSpecific heat capacity [J kg−1 K−1]
dRatio between the pre-consolidation pressure and the critical pressure
dTModel constant that determines the way the yield surface evolves with temperature
eVoid ratio
egVoid ratio at maximum pre-consolidation pressure
e0Initial void ratio
e0Initial void ratio
erfc(.)Complementary error function
fisoYield limit of the isotropic thermoplastic mechanism
gkPlastic potentials
kBorehole thermal conductivity [W m−1 K−1]
mMass of the borehole [kg]
nPorosity
pMean effective stress
p0Mean effective stress (reference value)
pcTPre-consolidation pressure (temperature dependent)
qDeviatoric stress
q0Deviatoric stress (reference value)
rbBorehole radius [m]
risoParameter related to the degree of plastification (isotropic yield limit)
rdevParameter related to the degree of plastification (deviatoric yield limit)
tTime [s]
uPorous water pressure
αGround thermal diffusivity [m2 s−1]
βVolumetric thermal expansion
βgVolumetric thermal expansion of the solid particles
βwVolumetric thermal expansion of water
εvVolumetric strain
ε s e Elastic shear strain
ε v e Elastic volumetric strain
ε e T p Thermoplastic volumetric strain
ε ˙ v e Elastic volumetric strain rate
ε ˙ v p Plastic volumetric strain rate
ε ˙ v T Volumetric thermal strain rate
ε ˙ v T Volumetric thermal strain rate
ε ˙ v T v p Volumetric thermal viscoplastic strain rate
ζ ˙ Rate of water flow in the soil voids
κSlope of elastic compression line
λSlope of plastic compression line
λkPlastic multipliers
μViscosity
ξRatio between the pre-consolidation pressure pc0 and the effective mean pressure
ρDensity [kg m−3]
vSpecific volume
ΨThermal elastic potential
DTensor of elastic moduli
ε ˙ e Elastic strain rate tensor
ε ˙ T Strain rate tensor
ε ˙ T p Thermoplastic strain rate tensor
ε ˙ d T p Thermoplastic deviatoric strain rate tensor
σ ˙ Stress rate tensor
ΨdDeviatoric plastic flow
ΨvIsotropic plastic flow

References

  1. Ilzarbe, J.M.B.; Malpartida, J.G.; Chandro, E.R. Recent Patents On Geothermal Power Extraction Devices. Recent Patents Eng. 2013, 7, 2–24. [Google Scholar] [CrossRef]
  2. Rees, S. Advances in Ground Source Heat Pump Systems; Woodhead Publishing: Sawston, UK, 2016. [Google Scholar]
  3. Shafagh, I.; Rees, S.J.; Mardaras, I.U.; Janó, M.C.; Carbayo, M.P. A Model of a Diaphragm Wall Ground Heat Exchanger. Energies 2020, 13, 300. [Google Scholar] [CrossRef] [Green Version]
  4. Bourne-Webb, P.J.; Pereira, J.-M.; Bowers, G.A.; Mimouni, T.; Loveridge, F.; Burlon, S.; Olgun, C.G.; McCartney, J.S.; Sutman, M. Design tools for thermoactive geotechnical systems. DFI J. J. Deep. Found. Inst. 2014, 8, 121–129. [Google Scholar] [CrossRef]
  5. Bourne-Webb, P.J.; Burlon, S.; Javed, S.; Kürten, S.; Loveridge, F. Analysis and design methods for energy geostructures. Renew. Sustain. Energy Rev. 2016, 65, 402–419. [Google Scholar] [CrossRef]
  6. Sani, A.K.; Singh, R.M.; Amis, T.; Cavarretta, I. A review on the performance of geothermal energy pile foundation, its design process and applications. Renew. Sustain. Energy Rev. 2019, 106, 54–78. [Google Scholar] [CrossRef]
  7. Ravera, E.; Sutman, M.; Laloui, L. Load Transfer Method for Energy Piles in a Group with Pile–Soil–Slab–Pile Interaction. J. Geotech. Geoenviron. Eng. 2020, 146, 04020042. [Google Scholar] [CrossRef]
  8. Parriaux, A.; Tacher, L. Utilisation de la Chaleur du Sol par des Ouvrages de Foundation et de Soutenement en Beton. Guide pour la Conception, la Realization et la Maintenance; Societe suisse des inenieurs et des architects, SIA: Zurich, Switzerland, 2005; (French and German). [Google Scholar]
  9. Pahud, D.; Fromentin, A. PILESIM: A simulation tool for pile and borehole heat exchanger systems. Bull. d’Hydrogéol. 1999, 17, 323–330. [Google Scholar]
  10. Brandl, H. Energy piles for heating and cooling of buildings. In Proceedings of the 7th International Conference Exhibition on Piling and Deep Foundations, Vienna, Austria, 15–17 June 1998. EU European Standards. [Google Scholar]
  11. Brandl, H.; Adam, D.; Markiewicz, R. Ground-sourced energy wells for heating and cooling of buildings. Acta Geotech. Slov. 2006, 3, 4–27. [Google Scholar]
  12. Brandl, H. Energy foundations and other thermo-active ground structures. Géotechnique 2006, 56, 81–122. [Google Scholar] [CrossRef]
  13. Katzenbach, R.; Olgun, C.G.; Loveridge, F.; Sutman, M.; Bowers, G.A.; McCartney, J.S.; Laloui, L.; Mimouni, T.; Dupray, F.; Spitler, J.D.; et al. New technologies and applications: Materials and equipment in near surface geothermal systems. DFI J. J. Deep. Found. Inst. 2014, 8, 93–107. [Google Scholar] [CrossRef]
  14. Adam, D.; Markiewicz, R. Energy from earth-coupled structures, foundations, tunnels and sewers. Géotechnique 2009, 59, 229–236. [Google Scholar] [CrossRef]
  15. Ghoreishi-Madiseh, S.A.; Kuyuk, A.F.; De Brito, M.A.R.; Baidya, D.; Torabigoodarzi, Z.; Safari, A. Application of Borehole Thermal Energy Storage in Waste Heat Recovery from Diesel Generators in Remote Cold Climate Locations. Energies 2019, 12, 656. [Google Scholar] [CrossRef] [Green Version]
  16. EN 15450. Heating Systems in Buildings—Design of Heat Pump Heating Systems; CEN/TC 228. 2007. Available online: https://shop.bsigroup.com/ProductDetail/?pid=000000000030144934 (accessed on 20 May 2020).
  17. VDI 4640. Thermal Use of the Underground, Part 2: Ground Source Heat Pump Systems; Beuth Verlag: Berlin, Germany, 2001. [Google Scholar]
  18. ÖWAV RB 207. Thermal use of the Groundwater and the Underground—Heating and Cooling; Austria. 2009. [Google Scholar]
  19. UNI. Sistemi Geotermici a Pompa di Calore—Requisiti per il Dimensionamento e la Progettazione; UNI: Milano, Italy, 2012. [Google Scholar]
  20. S. I. S.-F. CFMS. Recommandations pour la Conception, le Dimensionnement et la Mise en Oeuvre des Geostructures Paris, France Thermiques—Version 1. 2017. Available online: http://www.cfms-sols.org/sites/default/files/recommandations/Recommandations_Geostructures_Thermiques_CFMS_SYNTEC_Version_1_Janvier_2017.pdf (accessed on 20 May 2020).
  21. Kavanaugh, S.; Rafferty, K. Geothermal heating and cooling: Design of ground-source heat pump systems. ASHRAE Trans. 2014. [Google Scholar]
  22. Geotechnik, Deutsche Gesellschaft für Geowissen DGG. Shallow Geothermal Systems: Recommendations on Design, Construction, Operation and Monitoring; Wiley: New Jersey, NJ, USA, 2016. [Google Scholar]
  23. Ground Source Heat Pump Association GSHPA. Thermal Pile Design Installation and Materials Standards; Ground Source Heat Pump Association: Milton Keynes, UK, 2012; Issue 1.0. [Google Scholar]
  24. Geoexchange, B. Professional Guidelines for Geoexchange Systems in British Columbia, Part 2: Design; Geoexchange BC: Vancouver, BC, Canada, 2007. [Google Scholar]
  25. Peron, H.; Knellwolf, C.; Laloui, L. A Method for the Geotechnical Design of Heat Exchanger Piles. Geo-Frontiers 2011, 470–479. [Google Scholar] [CrossRef]
  26. Knellwolf, C.; Peron, H.; Laloui, L. Geotechnical Analysis of Heat Exchanger Piles. J. Geotech. Geoenviron. Eng. 2011, 137, 890–902. [Google Scholar] [CrossRef]
  27. Amatya, B.; Soga, K.; Bourne-Webb, P.J.; Amis, T.; Laloui, L. Thermo-mechanical behaviour of energy piles. Géotechnique 2012, 62, 503–519. [Google Scholar] [CrossRef]
  28. Xiao, J.; Luo, Z.; Martin, J.R.; Gong, W.; Wang, L. Probabilistic geotechnical analysis of energy piles in granular soils. Eng. Geol. 2016, 209, 119–127. [Google Scholar] [CrossRef]
  29. GABI Guideline. Geothermal energy Applications in Buildings and Infrastructure. 2020. Available online: https://www.foundationgeotherm.org/ (accessed on 20 May 2020).
  30. Raymond, J.; Therrien, R.; Gosselin, L.; Lefebvre, R. Numerical analysis of thermal response tests with a groundwater flow and heat transfer model. Renew. Energy 2011, 36, 315–324. [Google Scholar] [CrossRef]
  31. Bae, S.M.; Nam, Y.; Choi, J.M.; Lee, K.H.; Choi, J.S. Analysis on Thermal Performance of Ground Heat Exchanger According to Design Type Based on Thermal Response Test. Energies 2019, 12, 651. [Google Scholar] [CrossRef] [Green Version]
  32. Aresti, L.; Christodoulides, P.; Florides, G. A review of the design aspects of ground heat exchangers. Renew. Sustain. Energy Rev. 2018, 92, 757–773. [Google Scholar] [CrossRef]
  33. Buettner, K. Evaluation of soil heat conductivity with cylindrical test bodies. Trans. Am. Geophys. Union 1955, 36, 831–837. [Google Scholar] [CrossRef]
  34. Mogensen, P. Fluid to duct wall heat transfer in duct system heat storages. Geology 1983, 16, 652–657. [Google Scholar]
  35. Sanner, B.; Mands, E.; Sauer, M.K.; Grundmann, E. Thermal response test: A routine method to determine thermal ground properties for GSHP design. In Proceedings of the 9th International Energy Agency Heat Pump Conference, Zurich, Switzerland, 20–22 May 2008; pp. 04–35. [Google Scholar]
  36. Acuña, J. Improvements of U-Pipe Borehole Heat ExchI; KTH School of Industrial Engineering and Management: Stockholm, Sweden, 2010. [Google Scholar]
  37. Sharma, W.S.; Seki, A.; Angel, S.M.; Garvis, D. Field testing of an optical fiber temperature sensor in a geothermal. Geothermics 1990, 19, 285–294. [Google Scholar] [CrossRef]
  38. Jinguji, M.; Takehara, T.; Yamaguchi, T.; Kunimatsu, S. In-situ thermal conductivity exploration using penetrometer. J. Geother. Res. Soc. Jpn. 2002, 24, 349–356. [Google Scholar]
  39. Fujii, H.; Okubo, H.; Itoi, R. Thermal response tests using optical fiber thermometers. In Proceedings of the GRC 2006 Annual Meeting: Geothermal Resources-Securing Our Energy Future, San Diego, CA, USA, 10–13 September 2006; pp. 545–551. [Google Scholar]
  40. Márquez, M.I.V.; Raymond, J.; Blessent, D.; Philippe, M.; Simon, N.; Bour, O.; Lamarche, L. Distributed Thermal Response Tests Using a Heating Cable and Fiber Optic Temperature Sensing. Energies 2018, 11, 3059. [Google Scholar] [CrossRef] [Green Version]
  41. Fujii, H.; Okubo, H.; Nishi, K.; Itoi, R.; Ohyama, K.; Shibata, K. An improved thermal response test for U-tube ground heat exchanger based on optical fiber thermometers. Geothermics 2009, 38, 399–406. [Google Scholar] [CrossRef]
  42. Franco, A.; Conti, P. Clearing a Path for Ground Heat Exchange Systems: A Review on Thermal Response Test (TRT) Methods and a Geotechnical Routine Test for Estimating Soil Thermal Properties. Energies 2020, 13, 2965. [Google Scholar] [CrossRef]
  43. Acuña, J.; Palm, B. Distributed thermal response tests on pipe-in-pipe borehole heat exchangers. Appl. Energy 2013, 109, 312–320. [Google Scholar] [CrossRef]
  44. Dornstädter, J.; Heidinger, P.; Heinemann-Glutsch, B. Erfahrungen aus der praxis mit dem enhanced geothermal response test (EGRT). Karlsruhe, Germany. In Proceedings of the Der Geothermiekongress, Karlsruhe, Germany, 1–13 November 2008; pp. 271–279. [Google Scholar]
  45. Heske, C.; Kohlsch, O.; Dornstädter, J.; Heidinger, P. Der Enhanced-Geothermal-Response-Test als Auslegungsgrundlage und Optimierungstool. Bbr Fachmagazin Brunnen Leitungsbau 2011, 62, 36–43. (In German) [Google Scholar]
  46. Vieira, A.; Maranha, J.R.; Christodoulides, P.; Alberdi-Pagola, M.; Loveridge, F.; Nguyen, F.; Florides, G.; Radioti, G.; Cecinato, F.; Prodan, I.; et al. Characterisation of Ground Thermal and Thermo-Mechanical Behaviour for Shallow Geothermal Energy Applications. Energies 2017, 10, 2044. [Google Scholar] [CrossRef] [Green Version]
  47. Ingersoll, L.; Plass, H. Theory of the ground pipe source for the heat pump. ASHRAE Trans. 1948, 54, 339–348. [Google Scholar]
  48. Carslaw, H.; Jaeger, J. Conduction of Heat in Solids; Oxford University Press: Oxford, UK, 1959. [Google Scholar]
  49. Florides, G.; Kalogirou, S. Ground heat exchangers—A review of systems, models and applications. Renew. Energy 2007, 32, 2461–2478. [Google Scholar] [CrossRef]
  50. Bourne-Webb, P.; Freitas, T.M.B.; Gonçalves, R.D.C.; Bourne-Webb, P.J. Thermal and mechanical aspects of the response of embedded retaining walls used as shallow geothermal heat exchangers. Energy Build. 2016, 125, 130–141. [Google Scholar] [CrossRef]
  51. Al-Khoury, R. Computational Modeling of Shallow Geothermal Systems. In Computational Modeling of Shallow Geothermal Systems; CRC Press: Boca Raton, FL, USA, 2011; Volume 4. [Google Scholar] [CrossRef]
  52. Eskilson, P. Thermal Analysis of Heat Extraction Boreholes; University of Lund: Lund, Sweden, 1987. [Google Scholar]
  53. Kelvin, T. Mathematical and Physical Papers; Cambridge University Press: London, UK, 1982. [Google Scholar]
  54. Whitehead, S. Determining temperature distribution: A contribution to the evaluation of the flow of heat in isotropic media. Electrician 1927, 225–226. [Google Scholar]
  55. Hart, D.; Couvillion, R. Earth-Coupled Heat Transfer; National Water Well Association: Westerville, OH, USA, 1986. [Google Scholar]
  56. Christodoulides, P.; Florides, G.; Pouloupatis, P. A practical method for computing the thermal properties of a Ground Heat Exchanger. Renew. Energy 2016, 94, 81–89. [Google Scholar] [CrossRef]
  57. Javed, S.; Per, F.; Johan, C. Vertical ground heat exchangers: A review of heat flow models. In Proceedings of the Effstock 2009, Stockholm, Sweden, 14–17 June 2009. [CD-ROM]. [Google Scholar]
  58. Bandos, T.V.; Montero, A.; Fernández, E.; Gónzález-Santander, J.L.; Isidro, J.M.; Pérez, J.; De Córdoba, P.F.; Uchueguía, J. Finite line-source model for borehole heat exchangers: Effect of vertical temperature variations. Geothermics 2009, 38, 263–270. [Google Scholar] [CrossRef]
  59. Claesson, J.; Eskilson, P. Conductive heat extraction to a deep borehole: Thermal analyses and dimensioning rules. Energy 1988, 13, 509–527. [Google Scholar] [CrossRef]
  60. Lamarche, L.; Beauchamp, B. A new contribution to the finite line-source model for geothermal boreholes. Energy Build. 2007, 39, 188–198. [Google Scholar] [CrossRef]
  61. Zeng, H.Y.; Diao, N.R.; Fang, Z.H. A finite line-source model for boreholes in geothermal heat exchangers. Heat Transfer-Asian Res. 2002, 31, 558–567. [Google Scholar] [CrossRef]
  62. Jaeger, J. XVIII. Some problems involving line sources in conduction of heat. London, Edinburgh, Dublin Philos. Mag. J. Sci. 1944, 35, 169–179. [Google Scholar] [CrossRef]
  63. Hellström, G. Ground Heat Storage: Thermal Analyses of Duct Storage Systems. Ph.D. Thesis, Lund University, Lund, Sweeden, 1991. [Google Scholar]
  64. Li, M.; Lai, A.C.; Lai, A.C.K. New temperature response functions (G functions) for pile and borehole ground heat exchangers based on composite-medium line-source theory. Energy 2012, 38, 255–263. [Google Scholar] [CrossRef]
  65. Yang, Y.; Li, M. Short-time performance of composite-medium line-source model for predicting responses of ground heat exchangers with single U-shaped tube. Int. J. Therm. Sci. 2014, 82, 130–137. [Google Scholar] [CrossRef]
  66. Ingersoll, L.R.; Zobel, O.J.; Ingersoll, A.C. Heat Conduction with Engineering, Geological, and Other Applications. Phys. Today 1955, 8, 17. [Google Scholar] [CrossRef]
  67. Man, Y.; Yang, H.; Diao, N.; Liu, J.; Fang, Z. A new model and analytical solutions for borehole and pile ground heat exchangers. Int. J. Heat Mass Transf. 2010, 53, 2593–2601. [Google Scholar] [CrossRef]
  68. Monzo, P. Comparison of Different Line Source Model Approaches for Analysis of Thermal Response Test in a U-pipe Borehole heat Exchanger. Master’s Thesis, KTH School of Industrial Engineering and Management, Stockholm, Sweden, 2011. [Google Scholar]
  69. Kavanaugh, S. Simulation and Experimental Verification of Vertical Ground-Coupled Heat Pump Systems. Ph.D. Thesis, Oklahoma State University, Stillwater, OK, USA, 1985. [Google Scholar]
  70. Deerman, J.; Kavanaugh, S. Simulation of vertical U-tube ground-coupled heat pump systems using the cylindrical heat source solution. ASHRAE Trans. 1991, 97, 287–295. [Google Scholar]
  71. Kavanaugh, S.; Rafferty, K. Ground-Source Heat Pumps Design of Geothermal Systems for Commercial and Industrial Buildings; American Society of Heating: Peachtree Corners, GA, USA, 1997. [Google Scholar]
  72. Bernier, M. Ground-coupled heat pump system simulation/discussion. ASHRAE Trans. 2001, 107, 605. [Google Scholar]
  73. Gehlin, S. Thermal Response Test: Method Development and Evaluation. Ph.D. Thesis, Luleå University of Technology, Luleå, Sweden, 2002. [Google Scholar]
  74. Man, Y.; Yang, H.; Diao, N.; Cui, P.; Lu, L.; Fang, Z. Development of spiral heat source model for novel pile ground heat exchangers. HVAC&R Res. 2011, 17, 1075–1088. [Google Scholar]
  75. Cui, P.; Li, X.; Man, Y.; Fang, Z. Heat transfer analysis of pile geothermal heat exchangers with spiral coils. Appl. Energy 2011, 88, 4113–4119. [Google Scholar] [CrossRef]
  76. Park, S.; Lee, D.; Choi, H.-J.; Jung, K.; Choi, H. Relative constructability and thermal performance of cast-in-place concrete energy pile: Coil-type GHEX (ground heat exchanger). Energy 2015, 81, 56–66. [Google Scholar] [CrossRef]
  77. Loveridge, F.; Powrie, W. Temperature response functions (G-functions) for single pile heat exchangers. Energy 2013, 57, 554–564. [Google Scholar] [CrossRef]
  78. Fadejev, J.; Simson, R.; Kurnitski, J.; Haghighat, F. A review on energy piles design, sizing and modelling. Energy 2017, 122, 390–407. [Google Scholar] [CrossRef]
  79. Fadejev, J.; Kurnitski, J. Energy pile and heat pump modeling in whole building simulation model. In Proceedings of the 2nd IBPSA-England Conference on Building Simulation and Optimization, London, UK, 18 February 2014. [Google Scholar]
  80. Atam, E.; Helsen, L. Ground-coupled heat pumps: Part 1—Literature review and research challenges in modeling and optimal control. Renew. Sustain. Energy Rev. 2016, 54, 1653–1667. [Google Scholar] [CrossRef]
  81. Shonder, J.; Beck, J. Determining Effective Soil Formation Thermal Properties from Field Data Using a Parameter Estimation Technique; Oak Ridge National Lab: Oak Ridge, TN, USA, 1998. [Google Scholar]
  82. Verhelst, C.; Helsen, L. Low-order state space models for borehole heat exchangers. HVAC&R Res. 2011, 17, 928–947. [Google Scholar]
  83. Mihalakakou, G.; Santamouris, M.; Asimakopoulos, D. Modelling the thermal performance of earth-to-air heat exchangers. Sol. Energy 1994, 53, 301–305. [Google Scholar] [CrossRef]
  84. Yavuzturk, C.; Spitler, J.; Rees, S. A Transient two-dimensional finite volume model for the simulation of vertical U-tube ground heat exchangers. ASHRAE Trans. 1999, 105, 465–474. [Google Scholar]
  85. Yavuzturk, C. Modeling of Vertical Ground Loop Heat Exchangers for Ground Source Heat Pump Systems. Ph.D. Thesis, Oklahoma State University, Stillwater, OK, USA, 1999. [Google Scholar]
  86. Austin, W., III; Yavuzturk, C.; Spitler, J. Toward optimum sizing of heat exchangers for ground-source heat pump systems-development of an in-situ system and analysis procedure for measuring ground thermal. ASHRAE Trans. 2000, 106, 831–842. [Google Scholar]
  87. Bojić, M.; Trifunovic, N.; Papadakis, G.; Kyritsis, S. Numerical simulation, technical and economic evaluation of air-to-earth heat exchanger coupled to a building. Energy 1997, 22, 1151–1158. [Google Scholar] [CrossRef]
  88. Bi, Y.; Chen, L.; Wu, C. Ground heat exchanger temperature distribution analysis and experimental verification. Appl. Therm. Eng. 2002, 22, 183–189. [Google Scholar] [CrossRef]
  89. Gauthier, C.; Lacroix, M.; Bernier, H. Numerical simulation of soil heat exchanger-storage systems for greenhouses. Sol. Energy 1997, 60, 333–346. [Google Scholar] [CrossRef]
  90. Charlier, R.; Radu, J.-P.; Collin, F. Numerical modelling of coupled transient phenomena. Revue Française Génie Civil 2001, 5, 719–741. [Google Scholar] [CrossRef]
  91. Collin, F.; Li, X.; Radu, J.; Charlier, R. Thermo-hydro-mechanical coupling in clay barriers. Eng. Geol. 2002, 64, 179–193. [Google Scholar] [CrossRef] [Green Version]
  92. Amanzholov, T.; Akhmetov, B.; Georgiev, A.; Kaltayev, A.; Popov, R.; Dzhonova-Atanasova, D.; Tungatarova, M. Numerical modelling as a supplementary tool for Thermal Response Test. Bulg. Chem. Commun. 2016, 48, 109–114. [Google Scholar]
  93. Breger, D.S.; Hubbell, J.E.; El Hasnaoui, H.; Sunderland, J. Thermal energy storage in the ground: Comparative analysis of heat transfer modeling using U-tubes and boreholes. Sol. Energy 1996, 56, 493–503. [Google Scholar] [CrossRef]
  94. Ozudogru, T.Y.; Olgun, C.; Senol, A. 3D numerical modeling of vertical geothermal heat exchangers. Geothermics 2014, 51, 312–324. [Google Scholar] [CrossRef]
  95. Li, Z.; Zheng, M. Development of a numerical model for the simulation of vertical U-tube ground heat exchangers. Appl. Therm. Eng. 2009, 29, 920–924. [Google Scholar] [CrossRef]
  96. Gao, J.; Zhang, X.; Liu, J.; Li, K.; Yang, J. Numerical and experimental assessment of thermal performance of vertical energy piles: An application. Appl. Energy 2008, 85, 901–910. [Google Scholar] [CrossRef]
  97. Gang, W.; Wang, J. Predictive ANN models of ground heat exchanger for the control of hybrid ground source heat pump systems. Appl. Energy 2013, 112, 1146–1153. [Google Scholar] [CrossRef]
  98. Makasis, N.; Narsilio, G.; Bidarmaghz, A. A machine learning approach to energy pile design. Comput. Geotech. 2018, 97, 189–203. [Google Scholar] [CrossRef]
  99. Pahud, D.; Hubbuch, M. Measured thermal performances of the energy pile system of the dock midfield at Zürich airport. In Proceedings of the European Geothermal Congress 2007, Unterhaching, Germany, 30 May–1 June 2007. [Google Scholar]
  100. Kim, S.-K.; Bae, G.-O.; Lee, K.-K.; Song, Y. Field-scale evaluation of the design of borehole heat exchangers for the use of shallow geothermal energy. Energy 2010, 35, 491–500. [Google Scholar] [CrossRef]
  101. Emmi, G.; Zarrella, A.; De Carli, M.; Moretto, S.; Galgaro, A.; Cultrera, M.; Di Tuccio, M.; Bernardi, A. Ground source heat pump systems in historical buildings: Two Italian case studies. Energy Procedia 2017, 133, 183–194. [Google Scholar] [CrossRef]
  102. Fadejev, J.; Kurnitski, J. Geothermal energy piles and boreholes design with heat pump in a whole building simulation software. Energy Build. 2015, 106, 23–34. [Google Scholar] [CrossRef]
  103. Jun, L.; Zhang, X.; Jun, G.; Jie, Y. Evaluation of heat exchange rate of GHE in geothermal heat pump systems. Renew. Energy 2009, 34, 2898–2904. [Google Scholar] [CrossRef]
  104. Yang, H.; Cui, P.; Fang, Z. Vertical-borehole ground-coupled heat pumps: A review of models and systems. Appl. Energy 2010, 87, 16–27. [Google Scholar] [CrossRef]
  105. Yu, X.; Zhang, Y.; Deng, N.; Wang, J.; Zhang, N.; Wang, J. Thermal response test and numerical analysis based on two models for ground-source heat pump system. Energy Build. 2013, 66, 657–666. [Google Scholar] [CrossRef]
  106. Spitler, J.; Cullin, J.; Bernier, M.; Kummert, M.; Cui, P.; Liu, X.; Lee, E.; Fisher, D. Preliminary intermodel comparison of ground heat exchanger simulation models. In Proceedings of the 11th International Conference on Thermal Energy Storage, Stockholm, Sweden, 14–17 June 2009. [Google Scholar]
  107. Beier, R.A.; Smith, M.D. Minimum duration of in-situ tests on vertical boreholes. ASHRAE Trans. 2003, 109, 475. [Google Scholar]
  108. Beier, R.A.; Acuña, J.; Mogensen, P.; Palm, B. Transient heat transfer in a coaxial borehole heat exchanger. Geothermics 2014, 51, 470–482. [Google Scholar] [CrossRef]
  109. Lamarche, L.; Stanislaw, K.; Beauchamp, B. A review of methods to evaluate borehole thermal resistances in geothermal heat-pump systems. Geothermics 2010, 39, 187–200. [Google Scholar] [CrossRef]
  110. Claesson, J.; Javed, S. Explicit Multipole Formulas for Calculating Thermal Resistance of Single U-Tube Ground Heat Exchangers. Energies 2018, 11, 214. [Google Scholar] [CrossRef] [Green Version]
  111. Reda, F. Long term performance of different SAGSHP solutions for residential energy supply in Finland. Appl. Energy 2015, 144, 31–50. [Google Scholar] [CrossRef]
  112. De Rosa, M.; Ruiz-Calvo, F.; Corberan, J.M.; Montagud, C.; A Tagliafico, L. Borehole modelling: A comparison between a steady-state model and a novel dynamic model in a real ON/OFF GSHP operation. J. Phys. Conf. Ser. 2014, 547, 012008. [Google Scholar] [CrossRef] [Green Version]
  113. Thorén, Å. Practical Evaluation of Borehole Heat Exchanger Models in TRNSYS. Master’s Thesis, KTH, Stockholm, Sweden, 2016. [Google Scholar]
  114. Staiti, M.; Angelotti, A. Design of Borehole Heat Exchangers for Ground Source Heat Pumps: A Comparison between Two Methods. Energy Procedia 2015, 78, 1147–1152. [Google Scholar] [CrossRef] [Green Version]
  115. Al-Khoury, R.; Kölbel, T.; Schramedei, R. Efficient numerical modeling of borehole heat exchangers. Comput. Geosci. 2010, 36, 1301–1315. [Google Scholar] [CrossRef]
  116. Rühaak, W.; Renz, A.; Schätzl, P.; Diersch, H.-J. Numerical Modeling of Geothermal Applications. In Proceedings of the World Geothermal Congress 2010, Bali, Indonesia, 25–30 April 2010. [Google Scholar]
  117. Xuedan, Z.; Yiqiang, J.; Yang, Y. Improved method and case study of ground-coupled heat pump system design. Energy Procedia 2012, 14, 849–854. [Google Scholar] [CrossRef] [Green Version]
  118. He, M. Numerical Modelling of Geothermal Borehole Heat Exchanger Systems. Ph.D. Thesis, De Montfort University, Leicester, UK, 2012. [Google Scholar]
  119. Mogensen, P. Fullskaleförsök med berg som värmekälla för värmepump i Järfälla: Mätning och utvärdering; Statens råd för byggnadsforskning: Stockholm, Sweden, 1985. [Google Scholar]
  120. Witte, H.; van Gelder, A. Geothermal response tests using controlled multi-power level heating and cooling pulses (MPL-HCP): Quantifying ground water effects on heat transport around a borehole heat exchanger. In Proceedings of the Ecostock, the Tenth International Conference on Thermal Energy Storage, Galloway Township, NJ, USA, 31 May–2 June 2006. [Google Scholar]
  121. IEA ECES. Annex 21 Thermal Response Test (TRT); Final report; IEA Technology Collaboration Programmes: Paris, France, 2013. [Google Scholar]
  122. Wagner, R.; Clauser, C. Evaluating thermal response tests using parameter estimation for thermal conductivity and thermal capacity. J. Geophys. Eng. 2005, 2, 349–356. [Google Scholar] [CrossRef]
  123. Nelder, J.A.; Mead, R. A Simplex Method for Function Minimization. Comput. J. 1965, 7, 308–313. [Google Scholar] [CrossRef]
  124. Hooke, R.; Jeeves, T.A. Direct Search’’ Solution of Numerical and Statistical Problems. J. ACM 1961, 8, 212–229. [Google Scholar] [CrossRef]
  125. Wetter, M. GenOpt (R). In Generic Optimization Program, User Manual, Version 2.0.0; Lawrence Berkeley National Laboratory: Berkeley, CA, USA, 2003. [Google Scholar]
  126. Roth, P.; Georgiev, A.; Busso, A.; Barraza, E. First in situ determination of ground and borehole thermal properties in Latin America. Renew. Energy 2004, 29, 1947–1963. [Google Scholar] [CrossRef]
  127. Hemmingway, P.; Long, M. Energy piles: Site investigation and analysis. Proc. Inst. Civ. Eng. Geotech. Eng. 2013, 166, 561–575. [Google Scholar] [CrossRef]
  128. Shonder, J.; Beck, J. Field Test of a New Method for Determining Soil Formation Thermal Conductivity and Borehole Resistance; Oak Ridge National Lab.: Oak Ridge, TN, USA, 2000. [Google Scholar]
  129. Popov, R.; Georgiev, A.; Dzhonova-Atanasova, D. Stochastic search methods used for parameter estimation of thermal properties. In Proceedings of the World Geothermal Congress, Melbourne, Australia, 19–24 April 2015. [Google Scholar]
  130. Allaerts, K.; Coomans, M.; Salenbien, R. Hybrid ground-source heat pump system with active air source regeneration. Energy Convers. Manag. 2015, 90, 230–237. [Google Scholar] [CrossRef]
  131. Javed, S. Thermal Modelling and Evaluation of Borehole Heat Transfer. Ph.D. Thesis, Chalmers University of Technology, Göteborg, Sweden, 2012. [Google Scholar]
  132. Spitler, J.D. GLHEPRO-A design tool for commercial building ground loop heat exchangers. In Proceedings of the Fourth International Heat Pumps in Cold Climates Conference, Aylmer, QB, Canada, 17–18 August 2000. [Google Scholar]
  133. Crawley, D.B.; Lawrie, L.K.; Winkelmann, F.C.; Buhl, W.; Huang, Y.; Pedersen, C.O.; Strand, R.K.; Liesen, R.J.; Fisher, D.E.; Witte, M.J.; et al. EnergyPlus: Creating a new-generation building energy simulation program. Energy Build. 2001, 33, 319–331. [Google Scholar] [CrossRef]
  134. De Moel, M.; Bach, P.M.; Bouazza, A.; Singh, R.M.; Sun, J.O. Technological advances and applications of geothermal energy pile foundations and their feasibility in Australia. Renew. Sustain. Energy Rev. 2010, 14, 2683–2696. [Google Scholar] [CrossRef] [Green Version]
  135. Eldin, A.M.A.M.S.; Radwan, A.; Sakata, Y.; Katsura, T.; Nagano, K. The Effect of Groundwater Flow on the Thermal Performance of a Novel Borehole Heat Exchanger for Ground Source Heat Pump Systems: Small Scale Experiments and Numerical Simulation. Energies 2020, 13, 1418. [Google Scholar] [CrossRef] [Green Version]
  136. Lazzari, S.; Priarone, A.; Zanchini, E.; Antonella, P. Long-term performance of BHE (borehole heat exchanger) fields with negligible groundwater movement. Energy 2010, 35, 4966–4974. [Google Scholar] [CrossRef]
  137. Aresti, L.; Christodoulides, P.; Florides, G.A. Computational modelling of a ground heat exchanger with groundwater flow. Bulg. Chem. Commun. 2016, 48, 55–63. [Google Scholar]
  138. Sutton, M.G.; Nutter, D.W.; Couvillion, R.J. A Ground Resistance for Vertical Bore Heat Exchangers with Groundwater Flow. J. Energy Resour. Technol. 2003, 125, 183–189. [Google Scholar] [CrossRef]
  139. Diao, N.; Li, Q.; Fang, Z. Heat transfer in ground heat exchangers with groundwater advection. Int. J. Therm. Sci. 2004, 43, 1203–1211. [Google Scholar] [CrossRef]
  140. Pahud, D. PILESIM2, Simulation Tool For Heating/Cooling Systems with Heat Exchanger Piles or Borehole Heat Exchangers; Lugano, SUPSI Instory: Lugano, Switzerland, 2007. [Google Scholar]
  141. Zhang, W.; Yang, H.; Lu, L.; Fang, Z. The analysis on solid cylindrical heat source model of foundation pile ground heat exchangers with groundwater flow. Energy 2013, 55, 417–425. [Google Scholar] [CrossRef]
  142. Wang, D.; Lu, L.; Zhang, W.; Cui, P. Numerical and analytical analysis of groundwater influence on the pile geothermal heat exchanger with cast-in spiral coils. Appl. Energy 2015, 160, 705–714. [Google Scholar] [CrossRef]
  143. Muraya, N. Numerical Modeling of the Transient Thermal Interference of Vertical U-Tube Heat Exchangers. Ph.D. Thesis, Texas A&M University, College Station, TX, USA, 1994. [Google Scholar]
  144. Kohl, T.; Hopkirk, R. “FRACure”—A simulation code for forced fluid flow and transport in fractured, porous rock. Geothermics 1995, 24, 333–343. [Google Scholar] [CrossRef]
  145. Gashti, E.H.N.; Uotinen, V.-M.; Kujala, K. Numerical modelling of thermal regimes in steel energy pile foundations: A case study. Energy Build. 2014, 69, 165–174. [Google Scholar] [CrossRef]
  146. Dupray, F.; Laloui, L.; Kazangba, A. Numerical analysis of seasonal heat storage in an energy pile foundation. Comput. Geotech. 2014, 55, 67–77. [Google Scholar] [CrossRef]
  147. Gehlin, S.; Hellström, G. Influence on thermal response test by groundwater flow in vertical fractures in hard rock. Renew. Energy 2003, 28, 2221–2238. [Google Scholar] [CrossRef]
  148. Wang, H.; Qi, C.; Du, H.; Gu, J. Thermal performance of borehole heat exchanger under groundwater flow: A case study from Baoding. Energy Build. 2009, 41, 1368–1373. [Google Scholar] [CrossRef]
  149. Laloui, L.; di Donna, A. Understanding the thermo-mechanical behaviour of energy piles. Civil Eng. 2011, 164, 184–191. [Google Scholar]
  150. Adinolfi, M.; Maiorano, R.M.S.; Mauro, A.; Massarotti, N.; Aversa, S. On the influence of thermal cycles on the yearly performance of an energy pile. Geomech. Energy Environ. 2018, 16, 32–44. [Google Scholar] [CrossRef]
  151. Zymnis, D.M.; Whittle, A.J. Numerical Simulation of a Shallow Geothermal Heating/Cooling System. In From Soil Behavior Fundamentals to Innovations in Geotechnical Engineering; ASCE Press: Reston, VA, USA, 2014; pp. 2767–2776. [Google Scholar] [CrossRef] [Green Version]
  152. Mimouni, T.; Laloui, L. Thermo-pile: A numerical tool for the design of energy piles. In Energy Geostructures: Innovation in Underground Engineering; Willey: Hoboken, NJ, USA, 2013; Chapter 13; pp. 265–279. [Google Scholar]
  153. Laloui, L.; Nuth, M.; Vulliet, L. Experimental and numerical investigations of the behaviour of a heat exchanger pile. Int. J. Numer. Anal. Methods Geomech. 2006, 30, 763–781. [Google Scholar] [CrossRef]
  154. Bourne-Webb, P.J.; Amatya, B.; Soga, K.; Amis, T.; Davidson, C.; Payne, P. Energy pile test at Lambeth College, London: Geotechnical and thermodynamic aspects of pile response to heat cycles. Géotechnique 2009, 59, 237–248. [Google Scholar] [CrossRef]
  155. Cui, Y.; Sultan, N.; Delage, P. A thermomechanical model for saturated clays. Can. Geotech. J. 2000, 37, 607–620. [Google Scholar] [CrossRef]
  156. Maranha, J.R.; Pereira, C.; Vieira, A. Thermo-Viscoplastic Subloading Soil Model for Isotropic Stress and Strain Conditions. In Proceedings of the International Field Exploration and Development Conference 2019; Springer Science and Business Media LLC: Berlin, Germany, 2017; pp. 479–485. [Google Scholar]
  157. Laloui, L.; François, B. ACMEG-T: Soil Thermoplasticity Model. J. Eng. Mech. 2009, 135, 932–944. [Google Scholar] [CrossRef]
  158. Laloui, L.; Cekerevac, C.; François, B. Constitutive modelling of the thermo-plastic behaviour of soils. Revue Européenne Génie Civil 2005, 9, 635–650. [Google Scholar] [CrossRef]
  159. Di Donna, A.; Laloui, L. Response of soil subjected to thermal cyclic loading: Experimental and constitutive study. Eng. Geol. 2015, 190, 65–76. [Google Scholar] [CrossRef]
  160. Eriksson, L. Temperature effects on consolidation properties of sulphyde clays. In Proceedings of the 12th International Conference on Soil Mechanics and Foundation Engineering, Rio de Janeiro, Brazil, 13–18 August 1989; pp. 2087–2090. [Google Scholar]
  161. Sultan, N. Étude du Comportement Thermo-Méchanique de Lárgile de Boom: Expériences et Modélisation. Ph.D. Thesis, Ecole Nationale des Ponts et Chaussés, Paris, France, 1997. [Google Scholar]
  162. Kurz, D.; Sharma, J.; Alfaro, M.; Graham, J. Semi-empirical elastic–thermoviscoplastic model for clay. Can. Geotech. J. 2016, 53, 1583–1599. [Google Scholar] [CrossRef]
  163. Zhang, Z.-C.; Cheng, X. A thermo-mechanical coupled constitutive model for clay based on extended granular solid hydrodynamics. Comput. Geotech. 2016, 80, 373–382. [Google Scholar] [CrossRef]
  164. Maranha, J.R.; Pereira, C.; Vieira, A. Improved subloading thermo-viscoplastic model for soil under strictly isotropic conditions. Geomech. Energy Environ. 2018, 14, 38–47. [Google Scholar] [CrossRef]
  165. Mitchell, J.K. Fundamentals of Soil Behavior, 2nd ed.; Wiley International: New York, NY, USA, 1993. [Google Scholar]
  166. Vieira, A.; Maranha, J.R. Thermoplastic Analysis of a Thermoactive Pile in a Normally Consolidated Clay. Int. J. Geomech. 2017, 17, 04016030. [Google Scholar] [CrossRef]
  167. Hueckel, T.; Borsetto, M. Thermoplasticity of Saturated Soils and Shales: Constitutive Equations. J. Geotech. Eng. 1990, 116, 1765–1777. [Google Scholar] [CrossRef]
  168. Hueckel, T.; Baldi, G. Thermoplasticity of Saturated Clays: Experimental Constitutive Study. J. Geotech. Eng. 1990, 116, 1778–1796. [Google Scholar] [CrossRef]
  169. Robinet, J.-C.; Rahbaoui, A.; Plas, F.; Lebon, P. A constitutive thermomechanical model for saturated clays. Eng. Geol. 1996, 41, 145–169. [Google Scholar] [CrossRef]
  170. Abuel-Naga, H.M.; Bergado, D.T.; Bouazza, A.; Pender, M. Thermomechanical model for saturated clays. Géotechnique 2009, 59, 273–278. [Google Scholar] [CrossRef]
  171. Masad, E.; Muhunthan, B.; Chameau, J.L. Stress–strain model for clays with anisotropic void ratio distribution. Int. J. Numer. Anal. Methods Geomech. 1998, 22, 393–416. [Google Scholar] [CrossRef]
  172. Thomas, H.; He, Y.; Sansom, M.; Li, C. On the development of a model of the thermo-mechanical-hydraulic behaviour of unsaturated soils. Eng. Geol. 1996, 41, 197–218. [Google Scholar] [CrossRef]
  173. Hamidi, A.; Tourchi, S.; Kardooni, F. A critical state based thermo-elasto-plastic constitutive model for structured clays. J. Rock Mech. Geotech. Eng. 2017, 9, 1094–1103. [Google Scholar] [CrossRef]
  174. Hujeux, J.C. Calcul Numérique de Problémes de Consolidation Elastoplastique; École Centrale de Paris: Paris, France, 1979. [Google Scholar]
  175. Koiter, W.T. General theorems for elastic-plastic solids. In Progress of Solid Mechanics; North-Holland Publishing Company: Amsterdam, The Netherlands, 1960. [Google Scholar]
  176. Baldi, G.; Hueckel, T.; Pellegrini, R. Thermal volume changes of the mineral–water system in low-porosity clay soils. Can. Geotech. J. 1988, 25, 807–825. [Google Scholar] [CrossRef]
  177. Baldi, G.; Hueckel, T.; Peano, A.; Pellegrini, R. Developments in Modelling of Thermo-Hydro-Geomechanical Behaviour of Boom Clay and Clay-Based Buffer Materials; EUR 13365/1 and EUR 13365/2; Commission of the European Communities, Nuclear Science and Technology: Luxembourg, 1991. [Google Scholar]
  178. Modaressi, H.; Laloui, L. A THERMO-VISCOPLASTIC CONSTITUTIVE MODEL FOR CLAYS. Int. J. Numer. Anal. Methods Geomech. 1997, 21, 313–335. [Google Scholar] [CrossRef]
  179. Sällfors, G. Preconsolidation Pressure of Soft, High-Plastic Clays. Ph.D. Thesis, Chalmers University of Technology, Göteborg, Sweden, 1975. [Google Scholar]
  180. Fox, P.J.; Edil, T.B. Effects of stress and temperature on secondary compression of peat. Can. Geotech. J. 1996, 33, 405–415. [Google Scholar] [CrossRef]
  181. Xiong, Y.; Zhang, S.; Ye, G.; Zhang, F. Modification of thermo-elasto-viscoplastic model for soft rock and its application to THM analysis of heating tests. Soils Found. 2014, 54, 176–196. [Google Scholar] [CrossRef] [Green Version]
  182. Hashiguchi, K. Elastoplasticity Theory; Springer Science and Business Media LLC: Berlin, Germany, 2014. [Google Scholar]
  183. Zhang, F.; Yashima, A.; Nakai, T.; Ye, G.; Aung, H. An elasto-viscoplastic model for soft sedimentary rock based on tij concept and subloading yield surface. Soils Found. 2005, 45, 65–73. [Google Scholar]
  184. Maranha, J.R.; Pereira, C.; Vieira, A. A viscoplastic subloading soil model for rate-dependent cyclic anisotropic structured behaviour. Int. J. Numer. Anal. Methods Geomech. 2016, 40, 1531–1555. [Google Scholar] [CrossRef]
  185. Cekerevac, C.; Laloui, L. Experimental study of thermal effects on the mechanical behaviour of a clay. Int. J. Numer. Anal. Methods Geomech. 2004, 28, 209–228. [Google Scholar] [CrossRef]
  186. Cekerevac, C. Thermal Effects on the Mechanical Behaviour of Saturated Clays an Experimental and Constitutive Study; EPFL: Lausanne, Switzerland, 2003. [Google Scholar]
  187. Hellstrom, G.; Sanner, B. Software for dimensioning of deep boreholes for heat exctraction. In Proceedings of the Calorstock, Helsinki, Finland, 22–25 August 1994. [Google Scholar]
  188. Perego, R.; Pera, S.; Galgaro, A.; Santa, G.D.; Cultrera, M.; de Carli, M.; Emmi, G.; Bertermann, D.; Muller, J.; Mendrinos, D.; et al. Economic, geological and technical potential mapping test for GSHP systems in Europe. In Proceedings of the European Geothermal Congress, The Hague, The Netherlands, 11–14 June 2019. [Google Scholar]
  189. Kühl, L.; Renn, S.; Schneider, F.; Sanner, B.; Mands, E. Project geo:base—Energetical and ecological optimization of strategies for operation and control of complex energy supply systems based on shallow geothermal energy in commercial and non-residential buildings. In Proceedings of the European Geothermal Congress, The Hague, The Netherlands, 11–14 June 2019. [Google Scholar]
  190. Pinto, A.; Rodrigues, F.; Mota, A. Geothermal contribution on southern Europe climate for energy efficiency of university buildings. Winter season. Energy Procedia 2017, 134, 181–191. [Google Scholar] [CrossRef]
  191. Marcotte, D.; Pasquier, P. On the estimation of thermal resistance in borehole thermal conductivity test. Renew. Energy 2008, 33, 2407–2415. [Google Scholar] [CrossRef]
  192. Xu, X.; Spitler, J. Modeling of vertical ground loop heat exchangers with variable convectve resistance and thermal mass of the fluid. In Proceedings of the 10th International Conference on Thermal Energy Storage, Stockton, CA, USA, 31 May–2 June 2006. [Google Scholar]
  193. Claesson, J.; Hellström, G. Multipole method to calculate borehole thermal resistances in a borehole heat exchanger. HVAC&R Res. 2011, 17, 895–911. [Google Scholar]
  194. Diersch, H.-J.G. FEFLOW—Finite Element Modelling of Flow, Mass and Heat Transport in Porous and Fractured Media; Springer: Berlin/Heidelberg, Germany, 2014. [Google Scholar]
  195. Trefry, M.G.; Muffels, C. FEFLOW: A Finite-Element Ground Water Flow and Transport Modeling Tool. Ground Water 2007, 45, 525–528. [Google Scholar] [CrossRef]
  196. Diersch, H.; Bauer, D.; Heidemann, W. Finite element formulation for borehole heat exchangers in modeling geothermal heating systems by FEFLOW. WASY Softw. FEFLOW White Pap. 2010, 5, 5–96. [Google Scholar]
  197. Diersch, H.; Bauer, D.; Heidemann, W.; Rühak, W.; Schätzl, P. Finite element modeling of borehole heat exchanger systems: Part 1. Fundamentals. Comput. Geosci. 2011, 37, 1122–1135. [Google Scholar] [CrossRef]
  198. Diersch, H.; Bauer, D.; Heidemann, W.; Rühak, W.; Schätzl, P. Finite element modelling of borehole heat exchanger systems: Part 2. Numerical simulation. Comput. Geosci. 2011, 37, 1136–1147. [Google Scholar] [CrossRef]
  199. Zyvoloski, G. FEHM: A control Volume Finite Element Code for Simulating Subsurface Multi-Phase Multi-Fluid Heat and Mass Transfer; LAUR-07-3359; Los Alamos National Laboratory: Los Alamos, NM, USA, 2007. [Google Scholar]
  200. Parker, J.; Bose, J.; McQuiston, F. The ASHRAE Design/Data Manual for Ground-Coupled Heat Pumps; ASHRAE Transactions: Peachtree Corners, GA, USA, 1985. [Google Scholar]
  201. Kavanaugh, S.; Deerman, J. Simulation of vertical U tube ground coupled heat pump system. ASHRAE Trans. 1991, 97, 287–295. [Google Scholar]
  202. Paul, N.; Remund, C. The Effect of Grout Thermal Conductivity on Vertical Geothermal Heat Exchanger Design and Performance; Final Report No. TR-108529; Electric Power Research Institute: Washington, DC, USA, 1997. [Google Scholar]
  203. Cullin, J.R.; Spitler, J.D.; Montagud, C.; Ruiz-Calvo, F.; Rees, S.J.; Naicker, S.S.; Konečný, P.; Southard, L.E. Validation of vertical ground heat exchanger design methodologies. Sci. Technol. Built Environ. 2015, 21, 137–149. [Google Scholar] [CrossRef]
  204. Deng, Z.; Rees, S.J.; Spitler, J.D. A Model for Annual Simulation of Standing Column Well Ground Heat Exchangers. HVAC&R Res. 2005, 11, 637–655. [Google Scholar] [CrossRef]
  205. Hellström, G.; Sanner, B. PC-programs and modeling for borehole heat exchanger design. In Proceedings of the International Geothermal Days, Bad Urach, Germany, 17–22 September 2001. [Google Scholar]
  206. Shonder, J.; Hughes, P. Increasing Confidence in Geothermal Heat Pump Design Methods; Richard Stockton College: New Jersey, NJ, USA, 1998. [Google Scholar]
  207. Sanner, B.; Phetteplace, G.; Hellström, G. Introduction to computer models for geothermal heat pumps. In Proceedings of the International Geothermal Days, Klamath Falls, OR, USA , 10–16 October 1999. [Google Scholar]
  208. Shao, H.; Hein, P.; Sachse, A.; Kolditz, O. Geoenergy Modeling II; Springer Science and Business Media LLC: Berlin, Germany, 2016. [Google Scholar]
  209. Arnold, V. Ordinary Differential Equations; The MIT Press: Cambridge, MA, USA, 1978. [Google Scholar]
  210. Jung, Y.; Pau, G.H.; Finsterle, S.; Pollyea, R. A new efficient version of the TOUGH suite of multiphase flow and transport simulators. Comput. Geosci. 2017, 108, 2–7. [Google Scholar] [CrossRef]
  211. Peaceman, D. Fundamentals of Numerical Reservoir Simulation; Elsevier: Amsterdam, The Netherlands, 1977. [Google Scholar]
  212. Pruess, K. A Practical Method for Modeling Fluid and Heat Flow in Fractured Porous Media. Soc. Pet. Eng. J. 1985, 25, 14–26. [Google Scholar] [CrossRef] [Green Version]
  213. Pruess, K.; Narasimhan, T.N. On fluid reserves and the production of superheated steam from fractured, vapor-dominated geothermal reservoirs. J. Geophys. Res. Space Phys. 1982, 87, 9329–9339. [Google Scholar] [CrossRef]
  214. Pruess, K. GMINC: A Mesh Generator for Flow Simulations in Fractured Reservoirs; LBL-15227 ON: DE83014911; Lawrence Berkeley Lab: Berkley, CA, USA, 1983. [Google Scholar]
  215. Bernier, M.A.; Pinel, P.; Labib, R.; Paillot, R. A Multiple Load Aggregation Algorithm for Annual Hourly Simulations of GCHP Systems. HVAC&R Res. 2004, 10, 471–487. [Google Scholar] [CrossRef]
  216. Claesson, J.; Javed, S. A load-aggregation method to calculate extraction temperatures of borehole heat exchangers. ASHRAE Trans. 2012, 118, 530–539. [Google Scholar]
  217. Gnielinski, V. Neue Gleichungen für den Wärme- und den Stoffübergang in turbulent durchströmten Rohren und Kanälen. Forsch. Ingenieurwesen 1975, 41, 8–16. [Google Scholar] [CrossRef]
  218. Gröber, H.; Erk, S.; Grigull, U. Die Grundgesetzeder Wärmeübertragung; Springer: Berlin/Heidelberg, Germany, 1988. [Google Scholar]
  219. Krischer, O. Das temperaturfeld in der um-gebung von rohrleitungen, die in die erde verlegtsind. Gesundheits-Ingenieur 1936, 59, 537–539. [Google Scholar]
  220. Albers, K. Verhaltens-und elektrophysiologische Untersuchungen zur operanten Konditionierung von Antennenbewegungen der Honigbiene. Ph.D. Thesis, Technical University of Berlin, Berlin, Germany, 1991. [Google Scholar]
  221. Benkert, S.; Heidt, F.; Schöle, D. Calculation tool for Earth heat exchangers GAEA. In Proceedings of the Building Simulation 1997—The Fifth International IBPSA Conference, Prague, Czech Republic, 8–10 September 1997. [Google Scholar]
  222. Benkert, S.; Heid, F. Chapter 380—Designing Earth heat exchangers—Validation of the software GAEA. World Renew. Energy Congress 2000, VI, 1818–1821. [Google Scholar]
  223. Morrison, A. GS2000 software. In Proceedings of the Heat Pumps in Cold Climates—3rd International Conference, Wolfville, NS, Canada, 11–12 August 1997. [Google Scholar]
  224. Morrone, B.; Coppola, G.; Raucci, V. Energy and economic savings using geothermal heat pumps in different climates. Energy Convers. Manag. 2014, 88, 189–198. [Google Scholar] [CrossRef]
  225. Pahud, D.; Fromentin, A.; Hadorn, J. The Superposition Borehole Model for TRNSYS (TRNSBM). User Manuel, Internal Report; LASEN-EPFL: Lausanne, Switzerland, 1996. [Google Scholar]
  226. Pahud, D.; Fromentin, A.; Hadorn, J. The Duct Ground Heat Storage Model (DST) for TRNSYS Used for the Simulation of Heat Exchanger Piles; DGC-LASEN: Lausanne, Switzerland, 1996. [Google Scholar]
  227. Pola, M.; Fabbri, P.; Piccinini, L.; Zampieri, D.; Information, R. Conceptual and numerical models of a tectonically-controlled geothermal system: A case study of the Euganean Geothermal System, Northern Italy. Central Eur. Geol. 2015, 58, 129–151. [Google Scholar] [CrossRef]
  228. Laloui, L.; Moreni, M.; Vulliet, L. Comportement d’un pieu bi-fonction, fondation et échangeur de chaleur. Can. Geotech. J. 2003, 40, 388–402. [Google Scholar] [CrossRef]
  229. Laloui, L.; Nuth, M.; Vulliet, L. Experimental and Numerical Investigations of the Behavior of a Heat Exchanger Pile. In Elsevier Geo-Engineering Book Series; Elsevier: Amsterdam, The Netherlands, 2005; Volume 3, pp. 1065–1084. [Google Scholar] [CrossRef]
  230. Le Feuvre, P. An Investigation into Ground Source Heat Pump Technology, Its UK Market and Best Practise in System Design. Master’s Thesis, University of Stratchclyde, Glasgow, UK, 2007. [Google Scholar]
  231. Pahud, D.; Hellström, G.; Mazzarella, L. Duct Ground Heat Storage Model for TRNSYS (TRNVDST); Laboratory of Energy Systems: Lausanne, Switzerland, 1997. [Google Scholar]
Figure 1. Examples of SGES: (a) Energy pile, (b) Vertical GHE, (c) Surface pond/lake, (d) Horizontal GHE, (e) Tunnel lining, (f) Diaphragm wall, (g) Embankment.
Figure 1. Examples of SGES: (a) Energy pile, (b) Vertical GHE, (c) Surface pond/lake, (d) Horizontal GHE, (e) Tunnel lining, (f) Diaphragm wall, (g) Embankment.
Energies 13 04273 g001
Figure 2. EP boundary conditions and heat transfer as a case of SGES.
Figure 2. EP boundary conditions and heat transfer as a case of SGES.
Energies 13 04273 g002
Figure 3. (a) Sequential Forward-Evaluation, (b) Sequential Backward-Evaluation [121].
Figure 3. (a) Sequential Forward-Evaluation, (b) Sequential Backward-Evaluation [121].
Energies 13 04273 g003
Figure 4. Heat transfer interference between two boreholes. The direction of seepage water flow is upwards increasing the temperature of the second (top) borehole at steady state (modified from [137]).
Figure 4. Heat transfer interference between two boreholes. The direction of seepage water flow is upwards increasing the temperature of the second (top) borehole at steady state (modified from [137]).
Energies 13 04273 g004
Figure 5. (a) Temperature vs. volumetric strain [157,158]; (b) Temperature vs. pre-consolidation pressure [160].
Figure 5. (a) Temperature vs. volumetric strain [157,158]; (b) Temperature vs. pre-consolidation pressure [160].
Energies 13 04273 g005
Figure 6. (a) Coupled thermoplastic yield limits, (b) Numerical simulations and experimental results of a heating-cooling cycle at different degrees of consolidation [157].
Figure 6. (a) Coupled thermoplastic yield limits, (b) Numerical simulations and experimental results of a heating-cooling cycle at different degrees of consolidation [157].
Energies 13 04273 g006
Figure 7. Comparison between predicted and observed results: (a) on Pontida silty clay [176], (b) on Boom clay [177].
Figure 7. Comparison between predicted and observed results: (a) on Pontida silty clay [176], (b) on Boom clay [177].
Energies 13 04273 g007
Figure 8. Comparison of experimental results (from [178]) and numerical predictions with the proposed model: (a) Thermal loading/unloading under constant isotropic stress; (b) drained triaxial results and numerical simulations at 20 °C for two OCR values [178].
Figure 8. Comparison of experimental results (from [178]) and numerical predictions with the proposed model: (a) Thermal loading/unloading under constant isotropic stress; (b) drained triaxial results and numerical simulations at 20 °C for two OCR values [178].
Energies 13 04273 g008
Figure 9. (a) Differences in compression behavior with changes in strain rate [179], (b) Temperature versus volume strain during heating at constant isotropic pressure for various over-consolidation ratios (OCRs) from the ETVP model with the creep rate coefficient defined by an exponential relationship [162].
Figure 9. (a) Differences in compression behavior with changes in strain rate [179], (b) Temperature versus volume strain during heating at constant isotropic pressure for various over-consolidation ratios (OCRs) from the ETVP model with the creep rate coefficient defined by an exponential relationship [162].
Energies 13 04273 g009
Figure 10. Volumetric deformations of Kaolin clay under isotropic drained conditions laboratory results (discontinuous lines) and numerical simulations (a) heating for three different OCR values [184], (b) results for OCR = 6 (heating) and for OCR = 1 and 12 (heating and cooling cycle) [164].
Figure 10. Volumetric deformations of Kaolin clay under isotropic drained conditions laboratory results (discontinuous lines) and numerical simulations (a) heating for three different OCR values [184], (b) results for OCR = 6 (heating) and for OCR = 1 and 12 (heating and cooling cycle) [164].
Energies 13 04273 g010
Table 1. Comparison of the most common existing software tools for SGES modeling.
Table 1. Comparison of the most common existing software tools for SGES modeling.
SoftwareOpen Source/CommercialGeneral Notes on Method/Configuration/Time/Water/THMReferences
EED
(Earth Energy Designer)
https://www.buildingphysics.com
commercial
(Sweden, Germany, USA)
BHE thermal resistance calculation with multipoles;
stored g-functions obtained from numerical simulations (soil);
no CFD module included;
long-term models with monthly or hourly loads can be obtained from base or peak monthly loads;
surface water systems also included
[187,188,189]
(case studies)
EnergyPlus
https://energyplus.net
open source
(U.S. Department of Energy)
variable time-step model, uses both long time-step;
g-functions and short time-step;
g-functions for finite volume 1D radial soil;
whole building software with geothermal heat pump components;
monthly and hourly loads possible in short time-step g-functions;
no CFD module included
[190,191,192,193]
EWS
http://www.hetag.ch/
commercial
(Huber Energietechnik AG, Zurich)
finite difference model (EWS—Erdwärmesonden modelu); Eskilson’s g-functions used;
CFD module included;
single boreholes or fields of boreholes, hourly time step, up to 200 years, up to 600 boreholes in Pro version; double U-pipes, coaxial pipes, triangle, L-shape, square-shape, U-shape;
hydraulic linking of the BHE with the ventilation or the hydraulic cooling system possible;
influence of groundwater for one or two aquifers
FEFLOW
https://www.mikepoweredbydhi.com/products/feflow
commercial
(MIKE Powered by DHI Customer Success team, Denmark)
finite element analysis;
chemical kinetics for multi-component reaction systems;
able to calculate deep boreholes (700 m, for example);
optimized 3D unstructured CFD ground water model for pit dewatering plan and pore pressure estimation at geologically complex mine site
[194,195,196,197,198]
FEHM
(Finite Element Heat and Mass)
https://wikivividly.com/wiki/FEHM
open source
(Los Alamos National Laboratory); also, partly commercial
(SoilVision SVOFFICE™5/WR)
finite volume, finite difference or finite element method for heat and fluid flow (CFD);
finite elements for stress calculations;
not for standard BHEs but for groundwater consideration, from which the obtained heat can be extracted;
simulates groundwater and contaminant flow in deep and shallow porous media;
studies THM effects
[199]
GLD
(Ground Loop Design)
https://www.groundloopdesign.com/
commercial
(Thermal Dynamics
Inc., USA)
CFD module included;
infinite cylindrical source model (analytical) with few improvements, or finite line source model with g-functions (soil);
vertical, horizontal and pond systems, coaxial heat exchangers included;
surface water module for pond systems included
[200,201,202]
GLHE-Pro
https://hvac.okstate.edu/sites/default/files/pubs/glhepro/GLHEPRO_5.0_Manual.pdfspitler
commercial
(International Ground Source Heat Pump Association, Oklahoma State University)
multipole method (borehole);
finite difference g-function or free placement finite line source or standing column well model (soil);
no CFD module included;
max. 100 years; L-shaped, U-shaped, open rectangular and rectangular fields; max. 400 boreholes; vertical, horizontal and horizontal slinky systems
[132,193,203,204]
(validation)
GshpCalc
(name has been changed from Gchp—version 4.2 and earlier)
http://geokiss.com/
free to use, formerly commercial
(Geokiss, University of Alabama)
cylindrical source solution that uses four cyclic load pulses, multi-zones;
cooperated with the free software TideLoad to compute loads;
displays both imperial and metric units; only vertical GHEs;
surface water and groundwater included;
U-tubes (one, two or three) with series and parallel connection;
studies THM effects
[69,205,206,207]
IDA–ICE
(Indoor Climate and Energy)
http://www.equaonline.com
commercial
(EQUA Simulation AB, Stockholm)
finite difference 3D model;
no CFD module included;
multi-zone simulation application for studying thermal indoor climate as well as the energy consumption of the entire building;
no surface water considered;
pressure drop in the closed liquid circle calculated
[102]
(validation and case study)
OpenGeoSys
https://www.opengeosys.org/books/geoenergy-modeling-ii/
open sourcefinite elements, analytical infinite line source, numerical
infinite line source or numerical BHE models for soil simulations;
CFD module included;
U-tube, double U-tube, coaxial centered and coaxial annular configurations;
Groundwater flow simulated in porous and fractured media;
studies THM effects
[208,209]
(validation and case study)
TOUGH3
http://tough.lbl.gov/assets//files/Tough3/TOUGH3_Users_Guide_v2.pdf
commercial
(Berkeley Laboratory)
integral finite difference 1D, 2D or 3D method in space; irregular or regular discretization;
first order backward finite difference in time;
not for standard BHE, but for geothermal reservoir engineering;
developed for strongly heat driven flow;
water and air transport included
[210,211,212,213,214]
(case study)
TRNSYS
(Transient System Simulation Tool)
http://www.trnsys.com/
commercial
(Solar Energy Laboratory, University of Wisconsin, Madison)
finite difference model (EWS);
g-functions (two modules include DST model) for soil
whole building software with geothermal heat pump components;
global temperature + local solution + steady-flux solution (g-functions) = DST model;
no CFD module included
[130,203,215,216]
(case study)
Other models (older or redundant)
ECA
(Earth Coupled Analysis)
https://www.elitesoft.com/web/hvacr/ecaw.html
commercial
(Elite Software Development Inc.)
methodology described in Design/Data Manual for Closed-Loop Ground-Coupled Heat Pump Systems published ASHRAE;
vertical and horizontal pipes; built-in heat pump performance data;
no CFD module included
GAEA
(Graphical Design of Earth Heat Exchangers)
http://nesa1.uni-siegen.de/index.htm?/softlab/gaea.htm
commercial
(free full demo version limited to 10 days)
(University of Siegen, Germany)
developed according to models of Gnielisnki [217], Grober et al. [218], Krischer [219] and analytical model of Albers [220];
no vertical GHEs, and only basic horizontal GHEs;
no CFD module included
[221,222]
(validation)
GS2000
http://www.canetaenergy.com/research/software-development/
free to use
(Caneta Research, Canada)
line/cylindrical source models; vertical and horizontal configurations;
no CFD module included
[223]
HVACSIM+
https://nvlpubs.nist.gov/nistpubs/Legacy/IR/nistir7514.pdf
(NIST – National Institute of Standards and Technology, USA)capable to simulate entire HVAC building system
LoopLink
(Pro and RCL)
https://geoconnectionsinc.com/
commercial
(Geo-Connections Inc. SD, USA)
web-based; hourly simulation frequency
PILESIM2commercial
(Laboratory of Energy Systems—LASEN)
based on TRNSYS and models of TRNVDST and TRNSBM;
simulation time: 25 years
[9,99,224,225,226]
(case study)
Right-Loop
http://www.wrightsoft.com/icp/products/right-suite_universal/right-loop.aspx
commercial
(Wrightsoft Corporation, USA)
requires additional software, or purchased package
(right-j)
Ground and Geo-Structure Software
Code_Bright
https://deca.upc.edu/en/projects/code_bright
free to use (requires GiD pre-/post-prossessor)
(Universitat Politecnica de Catalunya)
FEM program capable of performing coupled THM analysis in geological media;
studies THM effects
[151]
HYDROTHERM
https://volcanoes.usgs.gov/software/hydrotherm/index.html
free to use
(U.S. Geological Survey—USGS)
can accommodate high fluid pressure and high temperatures;
multi-phase ground-water flow
[227]
Thermo-Pile
https://www.epfl.ch/labs/lms/wp-content/uploads/2018/08/Thermo-Pile_Documentation_Theory_V1-1.pdf
commercial
(École Polytechnique Ffédérale de Lausanne)
1D finite differences in radial direction; size of the piles; no water transfer taken into account; mesh analysis;
studies THM effects
[152,228,229]
(validation and case study)

Share and Cite

MDPI and ACS Style

Christodoulides, P.; Vieira, A.; Lenart, S.; Maranha, J.; Vidmar, G.; Popov, R.; Georgiev, A.; Aresti, L.; Florides, G. Reviewing the Modeling Aspects and Practices of Shallow Geothermal Energy Systems. Energies 2020, 13, 4273. https://doi.org/10.3390/en13164273

AMA Style

Christodoulides P, Vieira A, Lenart S, Maranha J, Vidmar G, Popov R, Georgiev A, Aresti L, Florides G. Reviewing the Modeling Aspects and Practices of Shallow Geothermal Energy Systems. Energies. 2020; 13(16):4273. https://doi.org/10.3390/en13164273

Chicago/Turabian Style

Christodoulides, Paul, Ana Vieira, Stanislav Lenart, João Maranha, Gregor Vidmar, Rumen Popov, Aleksandar Georgiev, Lazaros Aresti, and Georgios Florides. 2020. "Reviewing the Modeling Aspects and Practices of Shallow Geothermal Energy Systems" Energies 13, no. 16: 4273. https://doi.org/10.3390/en13164273

APA Style

Christodoulides, P., Vieira, A., Lenart, S., Maranha, J., Vidmar, G., Popov, R., Georgiev, A., Aresti, L., & Florides, G. (2020). Reviewing the Modeling Aspects and Practices of Shallow Geothermal Energy Systems. Energies, 13(16), 4273. https://doi.org/10.3390/en13164273

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop