Reviewing the Modeling Aspects and Practices of Shallow Geothermal Energy Systems

: Shallow geothermal energy systems (SGES) may take di ﬀ erent 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 ﬁrst presenting the main analytical and numerical models and methods related to the thermal behavior of SGES. Then, the most important supplementary factors a ﬀ ecting such modeling are discussed. These include: (i) the boundary conditions, in the form of temperature variation or heat ﬂow, that majorly a ﬀ ect the predicted thermal behavior of SGES; (ii) the spatial dimensions that may be crucial when relaxing the inﬁnite length assumption for short heat exchangers such as energy piles (EP); (iii) the determination of SGES parameters that may need employing speciﬁc techniques to overcome practical di ﬃ culties; (iv) a short-term vs. long-term analysis depending on the thermal storage characteristics of GHE of di ﬀ erent sizes; (v) the inﬂuence of groundwater that can have a moderating e ﬀ ect on ﬂuid 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 uniﬁed software package incorporating all related features of all SGES may be a future aim.


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

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]): where T b is the temperature [K] at r b [m], the borehole radius, T 0 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/(ρc p ) is the ground thermal diffusivity (m 2 s −1 ), ρ being the density (kg m −3 ) and c p 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 R b (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].

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

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 where J 0 , J 1 , Y 0 , Y 1 are Bessel functions of the first and second kind. Under certain circumstances Equation (3) can be approximated to a form containing R b (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.

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 Energies 2020, 13, 4273 7 of 45 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 Energies 2020, 13, 4273 8 of 45 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 km 2 (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.

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.

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, Energies 2020, 13, 4273 9 of 45 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.

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, 2noded, 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

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 Energies 2020, 13, 4273 11 of 45 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.

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 Energies 2020, 13, 4273 13 of 45 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).

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.

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.

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. Energies 2020, 13, x FOR PEER REVIEW 15 of 47 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 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 Energies 2020, 13, 4273 16 of 45 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%.

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]. 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]. 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  [157,158]; (b) Temperature vs. pre-consolidation pressure [160].
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 Energies 2020, 13, 4273 18 of 45 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.

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.

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].
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).
where K is the volumetric elastic tangent modulus, generally assumed temperature independent in the restricted range of temperatures, K f 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.

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 Energies 2020, 13, 4273 19 of 45 isothermal one. This tensor is responsible for the reversible thermal expansion and contraction of the material and is defined as where . ε T is the strain-rate tensor. This results in the stress rate tensor . σ: . σ = D : .
σ 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 where . ε e v 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 Ψ: 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, e 0 the initial void ratio, ε e v the elastic volumetric strain, ε e s 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 where β s0 is the isotropic thermal expansion coefficient at a reference temperature and ξ is the ratio between the preconsolidation pressure p cr0 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.

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, Energies 2020, 13, 4273 20 of 45 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.

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.

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. where and λ is the compressibility of the normal compression, a 0 , a 1 and a 2 are constant coefficients (usually negative, corresponding to the reduction of the semi-axis of the yield surface due to temperature alone), e 0 and e 1 values of the void ratio at the initial state at a hypothetical state corresponding to p c = p c0 , T = T 0 respectively, e g the void ratio at the maximum pre-consolidation pressure, ε v Tp the thermoplastic volumetric strain, and a 1 and a 3 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 where d T 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.

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 Energies where g k 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, f iso and f dev for the isotropic and deviatoric mechanisms respectively. In that case strains are produced.
Isotropic Thermo-Plastic Mechanism The yield limit f iso of the isotropic thermoplastic mechanism represented in the p -T plane is expressed by where r iso 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.
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 r dev 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.
Energies 2020, 13, x FOR PEER REVIEW 22 of 47 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 where 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.
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 pq space and is temperature dependent, and the parameter related to the degree of plastification (deviatoric surface).
The isotropic and the deviatoric yield limits are coupled through the hardening variable, , 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

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.

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.
Energies 2020, 13, x FOR PEER REVIEW 23 of 47 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.

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, and , (volumetric and distortional) of each mechanism k: where μ is a measure of viscosity which may vary with temperature, is the yield function of each mechanism, parameter related to deviatoric yield surface, parameter related to cyclic yield surface, a reference stress, isotropic plastic flow of each mechanism, and deviatoric plastic flow of each mechanism.
Some of the results obtained by this model under different stress paths are shown in Figure 8.

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 Tvp and . ε d Tp , (volumetric and distortional) of each mechanism 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, f 0 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) (b) 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]. 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). 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) (b) 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].
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). 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).

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.

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. 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 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] [132,193,203,204] (validation)  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) 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.

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.