Seasonal thermal energy storage in shallow geothermal systems: thermal equilibrium stage

This paper is dedicated to the study of seasonal heat storage in shallow geothermal installations in unsaturated soils for which hydrothermal properties such as degree of saturation and thermal conductivity vary with time throughout the profile. In the model, a semi-analytical model which estimates time-spatial thermal conductivity is coupled with a 2D cylindrical heat transfer modeling using finite difference method. The variation of temperature was obtained after 3 heating and cooling cycles for the different types of loads with maximum thermal load of qmax = 15 W.m−1 with variable angular frequency (8 months of heating and 4 months of cooling).and constant angular frequency (6 months of heating and 6 months of cooling) to estimate the necessary number of cycles to reach the thermal equilibrium stage. The results show that we approach a thermal equilibrium stage where the same variation of temperature can be observed in soils after several heating and cooling cycles. Based on these simulations, the necessary number of cycles can be related to the total applied energy on the system and the minimum number of cycles is for a system with the total applied energy of 1.9qmax .


INTRODUCTION
The thermal properties of the ground are generally recognized to be the most important parameters that affect heat transfer in soil.Most existing heat transfer calculation methods assume that the soil is homogeneous with constant thermal properties to avoid the mathematical complexity associated with modeling heat transfer in a non-homogenous medium.Unfortunately, data for soil thermal properties are often very difficult to obtain.Indeed, soil thermal properties are influenced by a myriad of factors such as soil type, soil density and soil water content.
The entire process of heat extraction/deposition is a transient one, due to the weather-dependent ground surface boundary conditions and heating/cooling load.The soil water content in close vicinity to the Borehole Heat Exchanger (BHE) can be influenced by numerous factors, such as: soil structure, temperature gradient, water gradient and irrigation (Krarti et al., 1995).In particular, the temperature gradient in the soil surrounding the BHE plays a significant role in the combined heat and water flow.essential to both the design and the operation of ground heat pump systems (Leong et al., 1998).Due to the very complex character of the ground, the actual design of the BHE should be based on a detailed mathematical model of simultaneous heat and water flow in soils.The parameters influencing the soil thermal properties are essential for testing the BHE length and carrying out a parametric analysis on the long-term effects of BHE operation on the surrounding ground (Ruševljan et al., 2009).The results of laboratory testing of soil samples (soil texture, soil dry density, water content) should also be available prior to the evaluation of soil thermal properties.
A few papers dealing with the influence of soil types and moisture conditions on ground heat pump performance have been published in the last decade (Leong et al., 1998 andSong et al., 2006).Hailey et al., (1990) analyzed the behavior of the thermal conductivity of soils adjacent to the BHE.
It was found that soil water content was a dominant factor responsible for seasonal thermal conductivity variations.
Another important point, less investigated in the literature, is the effect of several cooling and heating cycles on the performance of the BHE.If the surrounding ground may reach a thermal equilibrium state after several cooling and heating cycles, the same as the hydraulic equilibrium stage observed in the swelling clay soils after several wetting and drying cycles or the mechanical equilibrium stage observed in the granular soils after several loading and unloading cycles.
In this study, a finite difference model is used to simulate the temperature variations outside the shallow borehole in the ground while hydrothermal properties (such as degree of saturation and the texture of soil) of soil varied with time and depth.Due to the fact that the simulation is performed at shallow levels, a surface temperature variation is considered.
A time-spatial boundary condition is applied at the far depth of the borehole which estimates variation of temperature with time and depth.In the last part, the model is applied to case study of a Borehole Heat Exchanger (BHE) in Alsace region in France to investigate the effect of several cooling and heating cycles on the BHE operation.

NUMERICAL MODEL
Borehole Heat Exchanger (BHE) is considered as a geothermal system storage at shallow levels (less than 100 m).
Figure 1 shows a schematic geothermal system installed in a nonhomogenous soil while the groundwater level varies with time and depth.
Hydrothermal properties of the soil change with depth and time.Thus, a 2D axisymmetric model is performed with finite difference method.The borehole is subjected to a heat flux q.We assume that such heat flux q is equal to the heat flux q' per unit length at the borehole radius.
The heat transfer problem is governed by the following partial differential equation in cylindrical coordinates: where α is the thermal diffusivity, T is the temperature, t is the time and z is the depth.
The following equations are defined where i is the counter of r; j is the counter of z; and k is the counter of t. (2) (3) (4) (5) (6) By replacing above equations 2-5 in Equation 1, the temperature T at time k + 1 computed in each iteration until the difference of temperature between two iteration becomes less than 10 -4 .Normally for the convergence we should have αdt/(dr 2 + dz 2 ) < 0.5.
For the boundary conditions, we consider that at t = 0, T = T g0 .No heat flux was also considered at large depth (z=L):  Results from Abu-Hamdeh (2003) and Ochsner et al. (2001) show that there is no effect of initial dry density on the variation of thermal diffusivity.Therefore, we propose a linear relationship for saturated and dry thermal diffusivities (α sat and α dry ) with respect to sand content x s : sat=exs+f (11a) α αdry=gxs+h (11b) where e, f, g and h are constant parameters of values 4.5 × 10 −7 m 2 s −1 , 3.5 × 10 −7 m 2 s −1 , 2 × 10 −8 m 2 s −1 and 1.8 × 10 −7 m 2 s −1 , respectively.

Field description
Table 1 summarizes the geological and geotechnical description of the studied field in Hangenbieten such as the soil type, the sand content and the initial dry density (BRGM, 2000).W sat is the saturated water content calculated by the initial dry density value (γ d ).The λ sat and α sat of the different layers are calculated based on the equations 9, 10 and 11 by considering their sand content x s and initial dry density γ d .

Numerical simulation of a Borehole Heat Exchanger (BHE) installed in Hangenbieten
In this part, the storage capacity of a BHE (Figure 1    period of 3 years at surface (z = 0) was adapted from Nowamooz et al. (2015).At large depth L, there is no heat flux and soil temperature is considered as the same as the average temperature value of the ground (Tg0) at initial time t = 0 (equal to 10°).Two types of thermal load q as presented in figure 2  Figure 3 shows the variation of temperature with radius at the depth of 4.65 m for the maximum thermal loads applied after 4 months of heating and after 10 months (2 months of cooling).It can be observed that the thermal load with qmax=25 W/m produces the negative temperature in the soil during its heating phase.Because of the negative temperatures and to prevent the freezing soil phenomenon, the maximum thermal load of 15 W/m was selected for the following calculations.for the q max = 15 w/m.We can observe that there is a decrease of temperature after the first two cycles and then this decrease will be less important before reaching a thermal equilibrium stage where the soil presents the same temperatures after several heating and cooling cycles.

THERMAL EQUILIBRIUM STAGE
The maximum applied thermal load as well as the number cycles to reach the thermal equilibrium stage will be dependent on the climate of the studied regions.In this case, these findings are useful for the selection of the thermal load applied to the − 63 − shallow Borehole Heat Exchanger (BHE) installed in the different regions.

CONCLUSION
This paper is dedicated to the study of seasonal heat storage in shallow Borehole Heat Exchanger (BHE) for insitu soils in Alsace region in France.The thermal properties of each layer were considered as new input values in our proposed 2D finite difference model.Generally, we can observe that a thermal equilibrium stage can be reached after several (3 to 5) heating and cooling cycles.This equilibrium stage is important to evaluate the efficiency of shallow geothermal system in long period after several heating and cooling cycles.
Figure 1.Schematic diagram of a shallow geothermal system installed in a non-homogenous soil with variable groundwater level.
In this section, we studied the effect of two different thermal loads on the storage capacity of a field in Hangenbieten near to Strasbourg located in Alsace in France.Estimation of soil thermal diffusivity with regard to degree of saturation, texture, and dry densityNowamooz et al. (2015) proposed a set of relationships for the fitting functions of thermal conductivity as well as thermal diffusivity adapted to different soils in the literature etal., 2001, Abu-Hamdeh, 2003and Lu et al., 2007): the degree of saturation (-), x s is the sand content (-), γ d is the intial dry density (KN/m -3 ), Kλ is the normalized thermal conductivity, λ dry and λ sat are the thermal conductivity of dry and saturated soils (Wm -1 K -1 ), a and b are fitted parameters for saturated thermal conductivity equal to 0.53 Wm -1 K -1 and 0.1 Wm 2 K -1 .(kN) - , c and d are fitted parameters for unsaturated thermal conductivity equal to 0.087 Wm -1 K -1 and 0.019 Wm 2 K -1 .(kN) - , κ is dependent to the soil texture related to constant parameters e and f equal to 4.4 and 0.4.b)thermal diffusivityWe proposed the following equation for soil thermal diffusivity with regard to the degree of saturation by considering the saturated thermal diffusivity (α sat ) and the dry thermal diffusivity (α dry ): a constant parameter taken equal to 4.
) installed in Hangenbieten is presented.The height of the borehole (L) was taken 20 m.The radius of the borehole (rb) is 5 cm.The water table level (H) is 8 m.At the surface of soil, the temperature variation during the − 62 −

Figure 2 .
Figure 2. Two different thermal loads used in this study.

Figure 3 .
Figure 3. Variation of temperature with radius at the depth of 4.65 m for both thermal loads.

Figure 4 .
Figure 4. Equilibrium stage after several heating and cooling cycles a) with depth at the radius of 0.34m and b) with radius at the depth of 4.65m.
were applied to this borehole.For these loads, a period of 8 months of heating and a period of months of cooling were considered corresponding to the same existing period of heating and cooling in Alsace region.The maximum thermal loads during the heating phase (qmax) are 15 and 25 (W/m).The minimum values for the cooling phase (qmin =-1/3qmax) are respectively -5 and -8.33 W/m.

Figure 4 -
Figure 4-a and 4-b show the variation of temperature with depth (at the radius of 0.34 m) and with radius (at the depth of 4.65 m) after three successive heating and cooling cycles

Table 1 :
Profile of the multilayered subsurface in Hangenbieten (Alsace region).