This paper presents a simple and fast methodology to consider changing ther- mal properties in, for example, the design of shallow geothermal systems, by using thermal response factors. The simulation period is split into multiple sections according to the changes in thermal properties. By transforming the temperature response from one section to the next and subsequent su- perposition, any changes in the thermal conductivity properties of the media through which heat travels can be taken into account. The method is verified by numerical simulation and its efficiency is demonstrated in an application example. Results show that even though the computational effort increases exponentially with the variation of thermal parameters in time, the compu- tational time is significantly shorter than comparable numerical simulations. Keywords: Thermal response factor, g-function, thermal properties, variable conductivity, time dependent properties 1. Introduction Thermal response factors, also known as g-functions [1], are one of the main methods for designing shallow geothermal systems such as borehole heat exchangers, energy piles or horizontal geothermal collectors. The g- functions describe the dimensionless, time-dependent thermal resistance of ∗Corresponding author Email address: dueber@gut.rwth-aachen.de (Stephan D¨ uber) Preprint submitted to Elsevier November 20, 2023 the ground and allow the calculation of temperature changes due to heat conduction [2]: ∆T(t) =q 2πλ·g(t). (1) Here λis the thermal conductivity of the ground, qthe thermal load, and gis the g-function of the heat exchanger, depending on its geometry and a dimensionless time. Starting from the infinite line source [2], more and more analytical approaches have been developed in recent years to calculate the g-functions. For borehole heat exchangers and energy piles, the heat source radius [3], finite length and surface effects [4, 5], groundwater flow [6], layered soil types [7] and more realistic interaction between multiple boreholes [8] can now be considered. Approaches for horizontal collectors can also take into account different geometries such as the slinky heat exchanger [9] and the finite length of the pipes [10]. All of the approaches described above assume though that the thermal properties of the ground are constant over time. This assumption is generally justified in many cases, especially for deeper systems. However, it may be too simplistic for e.g., horizontal geothermal collectors or even energy piles, since the thermal properties may vary because of ground moisture content fluctuations due to seasonal changes or climate change [11]. In these cases, time-demanding, more complex, numerical models are necessary to obtain an accurate response. This paper fills this gap and presents a simple method to account for changing thermal properties using the thermal response factor method. The methodology is presented in detail in section 2, followed by a verification and application example in section 3, and finishing with conclusions in section 4. A Python implementation of the proposed method is available online. 2. Methodology The proposed method works with any type of g-function, but for illustra- tion purposes, we use here a single buried pipe, for which the g-function is calculated using the horizontal finite line source (HFLS) [10]: g(t) =Z∞ 1√ 4αt e−r2s2−e−(r2+4z2)s2 Hs2 Hserf(Hs)−1√π 1−e−H2s2 ds. (2) 2 Here αdenotes the thermal diffusivity of the ground, tthe time, Hthe length of the line source, zits depth and rthe horizontal distance where the temperature is evaluated. The temperature change ∆ T(t) due to a time- varying thermal load q(t) is calculated using the g-function method and the superposition principle [12] as: ∆T(tk) =1 2πλkX i=1∆q(ti)·g(tk−i+1) (3) where ∆ q(ti) =q(ti)−q(ti−1) is the load increment, λis the thermal conductivity of the ground and gthe g-function according to Equation 2. Neither Equation 3 nor the formulation of the g-function account for time dependent thermal properties of the ground. It means that, applying Equation 3 for a constant load q(t) for tstart≤t≤tcresults in a ∆ T(t) as shown in the top part of Figure 1. α1α = const. α=α21α<α21 α>α21α2 Figure 1: Temperature change ∆ Taccording to a constant load qbetween tstartandtcfor a typically assumed constant thermal parameter αof the heat transfer medium (top) and for the case when αchanges for t > t c(bottom). However, in reality, the thermal properties of the ground can change. Hence, if the thermal diffusivity αof the ground changes instantaneously at 3 t=tc, this change will only affect the temperature decay for t > t c, as shown conceptually in the bottom part of Figure 1. As the temperatures for t > t chave been calculated using Equation 3 for α1, we suggest that there is no need to recalculate them for α2. In fact, we propose that the change of αonly causes a compression or stretching of the already calculated temperature change along the time axis. Hence, we can write for t > t c: ∆T(α2, t) = ∆ T(α1, tα1 α2−tc(α1 α2−1)). (4) The factorα1 α2accounts for a time stretching/compression of the temper- ature curve while tc(α1 α2−1) corrects for the thereby introduced offset. Figure 2 shows conceptually what happens if the thermal property reduces from α1toα2at time t=tc(see solid line). The temperature will grow faster from t=tcuntil the time when q= 0, from which it will start a normal decay. The resulting curve is the superposition of the ∆ Tcalculated using Equation 3 for α1, from t=tstarttot=tc, plus the same temperature response transformed for t > t c, plus one calculated for α2from t=tc onward. α1 α<α21 Figure 2: Temperature change for a constant load with changing αattcas the sum of ∆T1and ∆ T2. The method described can be used for any type of thermal load profile and for any variation in thermal conductivity or heat capacity. For an efficient implementation, we use the Fast Fourier Transform (FFT) as presented by [13], which replaces the summation over time in Equation 3 with a single multiplication in the Fourier domain: 4 ∆T=F−1 F∆q 2πλ ·F(g) (5) whereFis the direct and F−1the inverse FFT. To ensure that Equation 4 can be applied for any change in α, the simulation period must be extended. An increase in αresults in a compression of the temperature response along the time axis, as can be seen in the bottom part of Figure 1. To ensure that the transformed temperature response covers the entire simulation period, the maximum time tmaxmust be increased to t′ maxas: t′ max= max( {f(αj) :j= 1, . . . , n }) (6) with: f(αj) =1 αjn−1X i=jαitc,i+1−nX i=j+1αitc,i+αntmax (7) where nstands for the number of changes in αandtc,ifor the times at which αchanges. Finally, the implementation follows the flowchart presented in Figure 3. 3. Verification and application To verify the described approach, we compare it to a simulation using a 2D numerical finite volume model. The numerical model is spatially dis- cretised with 600 ×600 cells each 10 ×10 mm in size. The temperatures on the cell walls are calculated using the central difference scheme while the explicit Euler scheme is used for the numerical integration. All parameters and boundary conditions are listed in Table 1. A large value was chosen for the length of the line source to minimise the influence of the finite length, which is not considered in the 2D numerical model. The radius in Table 1 denotes the horizontal distance from the line source where the temperature is evaluated. Using two thermal property changes over the calculation period, the results in Figure 4 show perfect agreement between the numerical simu- lation and our approach. The numerical simulation took 181 s on a standard personal computer, while the computational time for the analytical approach was just 1.5 s. 5 p > n? j > n?Extend time array (Eq. 6) j = 0 Set q = 0 outside current section Calc ΔT for whole time domain ΔT += ΔT section j section j p = j p += 1 Endj += 1Transform ΔT (Eq. 4) ΔT += ΔTsection p transformed, section p ΔT = ΔT transformedSpilt simulation time in n sections according to changes in α no noyes yesFigure 3: Flowchart for the algorithm to account for arbitrary variations of α. To show how the method performs for a realistic scenario, we consider a single pipe of a horizontal geothermal collector buried 1.2 m below the surface exposed to the hourly load profile shown in Figure 5. The load profile is derived from measurements of an operating geothermal system in Germany, all other parameters are given in Table 2. For the ground it is assumed that the thermal properties change sea- 6 Table 1: Parameters used for verification simulation. Parameter Symbol Value Units HFLS length H 500 m HFLS depth z 0.85 m Radius r 0.2 m Thermal load q 30 W m−1 Thermal conductivity λ 1.0 for t < t c,1 W m−1K−1 2.0 for tc,1≤t < t c,2 1.5 for t≥tc,2 Volumetric heat capacity ρc 1 000 000 J m−3K−1 Simulation time tmax 35 h Time step ∆ t 12.5 s Time of change 1 tc,1 41 675 s Time of change 2 tc,2 83 337.5 s λ=1 λ=2 λ=1.5 Figure 4: Comparison of the proposed method with and numerical finite volume simula- tion. sonally due to a change in the water content of the ground. The ground water content over central-western Europe can be approximated with a si- nusoidal curve [14]. Here we assume for the thermal conductivity λa value of 1.0 W m−1K−1for the dry and 2.2 W m−1K−1for the fully saturated case. The volumetric heat capacity ρcis considered to be 1.5 MJ m−3K−1for the dry and 2.2 MJ m−3K−1for the fully saturated case, resulting in the pro- file shown in Figure 6. To understand the importance of ground property variations, we approximate the profile of the thermal properties with daily, 7 /uni00000013 /uni00000015 /uni00000017 /uni00000019 /uni0000001b /uni00000014/uni00000013 /uni00000014/uni00000015 /uni00000037/uni0000004c/uni00000050/uni00000048/uni00000003/uni0000003e/uni00000050/uni00000057/uni0000004b/uni00000040/uni00000015/uni00000013 /uni00000014/uni00000013 /uni00000013/uni00000014/uni00000013/uni00000015/uni00000013/uni00000054/uni00000003/uni0000003e/uni0000003a/uni00000012/uni00000050/uni00000040Figure 5: Hourly load profile used for the application example. Table 2: Parameters used for the application example. Parameter Symbol Value Units HFLS length H 20 m HFLS depth z 1.2 m Radius r 0.02 m Thermal conductivity λ see Figure 6 Volumetric heat capacity ρc see Figure 6 Simulation time tmax 8760 h Time step ∆ t 1200 s monthly and three-monthly averages, resulting in 4, 12 and 365 values for α in the simulation. As we are only interested in the effect of the changes in the thermal properties, we neglect the course of the undisturbed ground tem- perature and restrict the analysis to the temperature change at the outside of the pipe introduced by the heat exchanger. The top part of Figure 7 shows the results for the simulation with aver- aged thermal properties ( λ= 1.6 W m−1K−1,ρc= 1.85 MJ m−3K−1) for the whole simulation period as reference. The lower parts show the temperature difference introduced by considering the different time resolution for changes in thermal properties. For the scenarios with 4 and 12 different values of α, the changes are clearly visible compared to the smooth curve for 365 val- ues of α. The temperature difference introduced by the change in thermal properties is more than 2 °C, which is almost half of the temperature change 8 Figure 6: Sinusoidal approximation of annual profile of the thermal properties. Daily, monthly and 3 monthly averaged values used as input for the simulation. caused by the heat exchanger when considering averaged properties (Figure 7, top). Finally, as Table 3 shows, the computational times of the model increase exponentially with the number of changes in α. This is because for each change, the response has to be transformed (Equation 4) for all subsequent values of α. However, even with daily updated values the total computation time is still only 44.8 s, clearly superior to numerical simulations. Table 3: Computational times depending on the number of values for α. nαcomp. time [s] 1 0.017 4 0.065 12 0.31 365 44.8 4. Conclusions The use of thermal response factors is widespread in the simulation of shallow geothermal systems. While the thermal properties of the ground are rightly assumed to be constant for the simulation of borehole heat exchangers, typically 100 m or deeper, they can vary closer to the surface for a number of reasons over the life-time operation. This change, caused, for example, by varying groundwater levels or moisture content, can be relevant for systems such as horizontal geothermal collectors or shorter energy piles. 9 /uni00000018 /uni00000013/uni00000018/uni00000037/uni0000000b/uni00000051/uni00000020/uni00000014/uni0000000c /uni00000013/uni00000015/uni00000037/uni00000048/uni00000050/uni00000053/uni00000048/uni00000055/uni00000044/uni00000057/uni00000058/uni00000055/uni00000048/uni00000003/uni0000003e/uni00000083/uni00000026/uni00000040/uni00000037/uni0000000b/uni00000051/uni00000020/uni00000017/uni0000000c/uni00000037/uni0000000b/uni00000051/uni00000020/uni00000014/uni0000000c /uni00000013/uni00000015/uni00000037/uni0000000b/uni00000051/uni00000020/uni00000014/uni00000015/uni0000000c/uni00000037/uni0000000b/uni00000051/uni00000020/uni00000014/uni0000000c /uni00000013 /uni00000015 /uni00000017 /uni00000019 /uni0000001b /uni00000014/uni00000013 /uni00000014/uni00000015 /uni00000037/uni0000004c/uni00000050/uni00000048/uni00000003/uni0000003e/uni00000050/uni00000057/uni0000004b/uni00000040/uni00000013/uni00000015/uni00000037/uni0000000b/uni00000051/uni00000020/uni00000016/uni00000019/uni00000018/uni0000000c/uni00000037/uni0000000b/uni00000051/uni00000020/uni00000014/uni0000000c Figure 7: Results of the application example with nαvalues for αfor the simulation period. The proposed approach is a simple, yet efficient, method to account for changing thermal properties using the response factor method. By dividing the simulation domain into several sections according to the changes in ther- mal properties and applying the superposition principle, arbitrary variations in thermal properties can be considered. The computational effort increases exponentially with the number of changes, but it is still considered superior to numerical simulations in many cases. 10