Modeling the Temperature Field in the Ground with an Installed Slinky-Coil Heat Exchanger

: In order to ﬁnd the temperature ﬁeld in the ground with a heat exchanger, it is necessary to determine temperature responses of the ground caused by heat sources and the inﬂuence of the environment. To determine the latter, a new model of heat transfer in the ground under natural conditions was developed. The heat ﬂux of the evaporation of moisture from the ground was described by the relationship taking into account the annual amount of rainfall. The analytical solution for the equations of this model is presented. Under the conditions for which the calculations were performed, the following data were obtained: the average ground surface temperature T sm = 10.67 ◦ C, the ground surface temperature amplitude A s = 13.88 K, and the phase angle P s = 0.202 rad. This method makes it possible to easily determine the undisturbed ground temperature at any depth and at any time. This solution was used to ﬁnd the temperature ﬁeld in the ground with an installed slinky-coil heat exchanger that consisted of 63 coils. The results of calculations according to the presented model were compared with the results of measurements from the literature. The 3D model for the ground with an installed heat exchanger enables the analysis of the inﬂuence of miscellaneous parameters of the process of extracting or supplying heat from/to the ground on its temperature ﬁeld. This presents a new model of heat transfer in the ground under natural conditions. The heat ﬂux for the of was described with a linearized relationship taking account the annual of rainfall. The analytical solution of the model equations enables users to easily determine the parameters of Equation (2): T sm , A s , P s , and hence the undisturbed ground temperature proﬁle. The proﬁle was used to ﬁnd the temperature ﬁeld in the ground with an installed slinky-coil exchanger. When determining the ground temperature ﬁeld, the intermediate results of calculations of ground temperature responses ϑ and undisturbed ground temperature T 0 were used to compare their time courses.


Introduction
Among the many known types of heat exchangers, ground heat exchangers play an important role in installations with renewable energy sources. In recent years, there has been a significant development in terms of new designs and new ways of modeling and designing these devices. This progress is the result of numerous theoretical and experimental research studies on heat transfer between the working fluid flowing through the exchanger and the ground where the exchanger pipes are installed. Various ways of arranging the pipes of the ground heat exchangers, both horizontally and with vertical boreholes, have been developed. Recently, there have been attempts to use nanofluids as working fluids [1], like in the case of heat transfer between fluids [2]. Another modification involves the introduction of phase change materials with the ability to accumulate heat [3].
The ground heat exchangers are parts of the heat pump installations. The advantages of exchangers with horizontal pipes, compared to exchangers with vertical pipes, include lower investment costs and the use of the environment to supplement thermal deficiencies occurring in the ground (or to get rid of the excess heat).
Although many devices for monitoring the operation of ground heat exchangers have been installed in recent years around the world, few of them have been used to study the thermal behavior of the ground. Most of the measurement data relate to the obtained COP values. A small amount of data of temperature of the ground with an installed heat exchanger leads to an incomplete understanding of how the ground responds to extracting/supplying heat from/to a ground heat pump system.
Hepburn et al. [4] conducted studies on the ground temperature near a ground heat exchanger coupled with a heat pump. Their observations showed that most of the ground temperature changes near the heat exchanger pipes were caused by the variable temperature of the ground surface, and not by the heat extraction.
In the work by Wu et al. [5], measurements of the temperature of the ground near a slinky-coil heat exchanger coupled with a heat pump for space heating were presented. The horizontal spiral coils of the exchanger pipes were arranged in four rows, each one was 1 m wide and 80 m long, they were placed at a depth of 1.14 m.
Pauli et al. [6] monitored the ground temperature with installed horizontal ground heat exchangers of various configurations. The research concerned the collection of heat from the ground. The conclusions from the research were qualitative: (a) the ground temperatures rarely dropped below zero during most of the heating season, (b) the lowest temperatures occurred near the ground surface and not at the level of the heat exchanger installation, and (c) only at the end of the heating season did the ground temperature drop below the ambient temperature.
Neuberger and Adamovsky [7,8] conducted measurements of the ground temperature near the slinky-coil and straight pipe exchangers. The research was conducted both during the heating season and for the thermal regeneration of the ground. It was found that, in contrast to the exchanger with straight pipes, the ground temperature near the slinky-coil exchanger dropped below 0 • C during the heating season.
Jeon [9] compared the operation of the slinky-coil and straight pipe exchangers using the concept of a scale factor that was defined as the ratio of the mean thermal energy transferred for the spiral-coil-type GHE to that for the straight-line-type GHE.
The results of ground temperature measurements near the ground heat exchangers are often not confronted with appropriate mathematical models of heat transfer processes. Therefore, an attempt was made to create such a model.
The superposition principle [10] is the basis for modeling heat transfer in the ground with an installed heat exchanger. The complicated heat transfer process in the pipe system is composed of two independent elementary heat transfer processes. The first one is the heat transfer between the source (the surface of the pipes) and the semi-infinite medium (ground) with an isothermal boundary. The second process is the transfer of heat in this medium with a changing temperature on its surface. Therefore, to calculate the temperature of the ground, the temperature responses of the ground due to heat sources ϑ and the temperature responses of the varying temperature of the ground surface ϑ 0 must be taken into account: In the absence of heat sources, i.e., for ϑ = 0, the ground temperature is T 0 = T init + ϑ 0 , and therefore: In designing and modeling of the ground heat exchangers, it is necessary to know the ground temperature field T 0 under natural conditions, i.e., in the ground without a heat exchanger installed. The temperature under these conditions varies with time and location and is called the undisturbed ground temperature (the same name is also used for the ground temperature at great depths; however, in this work, the undisturbed ground temperature should be understood as the temperature previously described).
The solution of the heat conduction equation for a semi-infinite body with periodic changes of its surface temperature is known (Carslaw and Jaeger [10]): The results of extensive studies on the ground temperature at various depths under different climatic conditions in the USA were presented by Kusuda and Achenbach [11]. Based on the measurement results, the coefficients of Equation (2) were determined and the results were in good agreement with the model. This paper presents a new model of heat transfer in the ground under natural conditions. The heat flux for the evaporation of moisture from the ground was described with a linearized relationship taking into account the annual amount of rainfall. The analytical solution of the model equations enables users to easily determine the parameters of Equation (2): T sm , A s , P s , and hence the undisturbed ground temperature profile. The profile was used to find the temperature field in the ground with an installed slinky-coil exchanger. When determining the ground temperature field, the intermediate results of calculations of ground temperature responses ϑ and undisturbed ground temperature T 0 were used to compare their time courses.
The aim of this work was to present and implement a model describing the temperature field for the ground in which a slinky-coil exchanger is placed. The simulation concerned the climatic conditions for East-Central Europe. Moreover, the results of calculations according to the presented model were compared with the results of measurements from the literature.
The model for the ground with an installed heat exchanger enables the analysis of the influence of various parameters of the process of extracting or supplying heat from/to the ground on its temperature field. The knowledge of the temperature field near the ground heat exchanger is important for the intensive extraction or supply of heat from/to the ground. Both excessive cooling and excessive ground overheating are unfavorable for the further operation of the exchanger.

Heat Fluxes on the Ground Surface
The following heat fluxes occur on the ground surface: convective flux, radiant flux, moisture evaporation flux, and conductive flux.

Convective Flux
The convective flux can be described by the heat transfer equation: The average daily values of air temperature and ground surface temperature change with time according to the following dependencies: The heat transfer coefficient between the ground surface and the surroundings, h, depends on the wind speed at the ground surface and the type of ground surface covering. Among the many empirical relationships from the literature, the following formula is often used to determine the h coefficient: where ρ a and c pa are the density and heat capacity of air, and r a is the aerodynamic resistance to heat transfer. For example, when the ground is covered with 0.12 m high plants, the following dependency can be used [12]: where u is the wind velocity 2 m above the ground.

Radiant Flux
The average daily solar radiation (short-wave) flux absorbed by the ground changes in the annual cycle according to the relationship: The long-wave radiation flux refers to the transport of heat from the ground surface to the sky (Earth radiation). Since the surface temperature of the ground T s is higher than the temperature of the sky T sky , the ground radiates thermal energy. Due to the low temperature of the radiation source, it is the long-wave type radiation. The average daily net flux LW of the long-wave radiation is: The relationship (9) was transformed into the linear form with respect to the temperature of the ground surface [13]. Since: where T m is the mean from the interval (T sky , T s ). Thus, the average daily long-wave radiation flux is: where: Figure 1 shows the relationship between the expressions on the right-hand sides of Equations (9) and (11). For the T s − T sky interval between 10 and 20 K and the value of T s ≈ 10 • C, the relationship was approximated by the linear equation (coefficient of determination R 2 = 0.99975, standard deviation σ = 0.20 × 10 8 K 4 ). On the basis of the coefficient of this equation, the value of T m in Formula (10) was determined: T m = 275.1 K and the value of the constant C LW was found C LW = 4.72 W/(m 2 K).

Radiant Flux
The average daily solar radiation (short-wave) flux absorbed by the ground in the annual cycle according to the relationship: The long-wave radiation flux refers to the transport of heat from the ground to the sky (Earth radiation). Since the surface temperature of the ground Ts is high the temperature of the sky Tsky, the ground radiates thermal energy. Due to temperature of the radiation source, it is the long-wave type radiation. The avera net flux LW of the long-wave radiation is: The relationship (9) was transformed into the linear form with respect temperature of the ground surface [13]. Since: Therefore: where Tm is the mean from the interval (Tsky, Ts). Thus, the average daily lon radiation flux is: where: Figure 1 shows the relationship between the expressions on the right-hand Equations (9) and (11). For the Ts − Tsky interval between 10 and 20 K and the valu 10 °C, the relationship was approximated by the linear equation (coeffic determination R 2 = 0.99975, standard deviation σ = 0.20 × 10 8 K 4 ). On the basi coefficient of this equation, the value of Tm in Formula (10) was determined: Tm = and the value of the constant CLW was found CLW = 4.72 W/(m 2 K). The temperature of the sky Tsky is related to the air temperature by the relat resulting from the Stefan-Boltzmann equation [14,15]:  The temperature of the sky T sky is related to the air temperature by the relationship resulting from the Stefan-Boltzmann equation [14,15]: where T sky and T a are expressed in [K] and ε sky is the effective emissivity of the sky. In order to determine ε sky , the results of the measurements presented in the work by Bryś et al. were used [16], which showed that the average annual net long-wave radiation flux is LW = 45.7 W/m 2 . This value applies to the ground surface not covered with plant under the conditions of the measurements (Wrocław, Poland). For calculations, the average annual air temperature and the average annual temperature of the ground surface under similar climatic conditions were adopted: T am = 8.95 • C [17] and T sm ≈ 10 • C. The unit emissivity of the ground surface was assumed and the annual average temperature of the sky was calculated from Formula (9): T sky,m = 273.5 K, and the effective emissivity of the sky-according to Formula (13): ε sky = 0.893.

Evaporative Flux
Heat is also transferred between the ground surface and its surroundings as a result of the evaporation of water contained in the ground. Convective transfer coefficients are related to each other by the Chilton-Colburn heat and mass transfer analogy [18]. The resulting Penman-Monteith equation [19,20] determines the potential of the evaporated water stream: where ∆ = dp sat /dT, γ 0 (= 59.5 Pa/K)-psychrometric constant, r c -crop canopy resistance, p sat,s and p-partial pressures of water vapor on the ground surface and in the atmospheric air, respectively. Equation (14) relates to the potentially possible rate of water evaporation from the ground. The actual evaporation rate E may be lower than E 0 due to the lack of moisture in the ground. To estimate the actual rate of evaporation, the inflow of water to the ground, mainly determined by the annual rainfall P, should be taken into account. Therefore, the relationship between the aridity index Φ = E 0 /P and the E/P ratio is considered, i.e., the Budyko curve [21][22][23][24]. For Φ < 1 the evaporation rate is limited by E 0 (i.e., by the energy supply), while for Φ > 1 the evaporation rate is limited by the available amount of water P. Pike [23] proposed the following empirical relationship: It should be noted that most of the experimental studies on which the Budyko curve is based concerned the range 0.5 < Φ < 2.0 and showed a significant dispersion. Figure 2 shows the approximations of the Budyko curve according to the research of various authors. The relationship (15) was approximated by a straight line for the range of Φ from 0.5 to 2, and then transformed to the form: where C 1 = 0.287 and C 2 = 0.357. From this dependence, it can be seen what share in the actual flux of the evaporated water has the evaporation potential E 0 and the annual amount of precipitation P.
In order to express the heat flux for water evaporation in a linear form with respect to the temperature of the ground surface, the dependence between the saturated water vapor pressure and temperature was approximated with a linear equation [25]  where a P = 103 Pa/K, b P = 609 Pa. Taking into account dependencies (14), (16), and (17) and that EV = L evap ·E, one can obtain the expression for the evaporative heat flux: where: Energies 2021, 14, x FOR PEER REVIEW 6 of 19 Figure 2. Linearization of the relationship between the E/P ratio and the aridity index Ф.
In order to express the heat flux for water evaporation in a linear form with respect to the temperature of the ground surface, the dependence between the saturated water vapor pressure and temperature was approximated with a linear equation [25] where aP = 103 Pa/K, bP = 609 Pa. Taking into account dependencies (14), (16), and (17) and that EV = Levap•E, one can obtain the expression for the evaporative heat flux: where: Expression (18) consists of three components. The first is the result of the mass transfer driving force, i.e., the difference in partial pressures of water vapor on the ground surface and in the bulk of the gas phase. The second component depends on the heat supplied to the system externally (from radiation). The third component takes into account the effect of the amount of water supplied to the system in the form of rainfall. The significance of the latter effect was demonstrated by the authors in [26] by formulating an empirical equation. Figure 3 shows the effect of the average annual solar radiation flux Sm and the annual precipitation P on the time courses of the evaporative heat flux EV. The increase in the intensity of solar radiation, as well as the increase in the amount of precipitation, causes the evaporative heat flux to increase. The maximum value of the evaporative flux EV occurs in the period of the year close to the period with the maximum temperature of the ground surface. Expression (18) consists of three components. The first is the result of the mass transfer driving force, i.e., the difference in partial pressures of water vapor on the ground surface and in the bulk of the gas phase. The second component depends on the heat supplied to the system externally (from radiation). The third component takes into account the effect of the amount of water supplied to the system in the form of rainfall. The significance of the latter effect was demonstrated by the authors in [26] by formulating an empirical equation. Figure 3 shows the effect of the average annual solar radiation flux S m and the annual precipitation P on the time courses of the evaporative heat flux EV. The increase in the intensity of solar radiation, as well as the increase in the amount of precipitation, causes the evaporative heat flux to increase. The maximum value of the evaporative flux EV occurs in the period of the year close to the period with the maximum temperature of the ground surface.

Thermal Balance on Ground Surface-The Analytical Solution
The heat flux q cond resulting from the algebraic summation of the above-mentioned fluxes is transported (as a conductivity flux) deep into the ground or in the opposite direction depending on the sign. With the use of the average daily values, the heat balance on the ground surface can be written as follows [13,[25][26][27][28][29]: On the other hand, the balance equation for the average annual fluxes has the form: Under natural conditions, the ground is in a cyclic steady state, and therefore the average annual conductivity flux is zero:

Thermal Balance on Ground Surface-The Analytical Solution
The heat flux qcond resulting from the algebraic summation of the above-mentioned fluxes is transported (as a conductivity flux) deep into the ground or in the opposite direction depending on the sign. With the use of the average daily values, the heat balance on the ground surface can be written as follows [13,[25][26][27][28][29]: On the other hand, the balance equation for the average annual fluxes has the form: Under natural conditions, the ground is in a cyclic steady state, and therefore the average annual conductivity flux is zero: On the basis of the relationships for the individual heat fluxes, presented in Sections 2.1-2.3, the average annual fluxes Hm, (LW)m, and (EV)m can be expressed as a function of the average annual temperature of the ground surface Tsm, and hence after the substitution to Equation (21), one can determine this temperature: where: The values of the amplitude As and the phase shift Ps of the ground surface temperature can be determined on the basis of the balance of average daily fluxes (20). For the conductivity flux, it should be substituted in accordance with the Fourier equation:  (21), one can determine this temperature: where: The values of the amplitude A s and the phase shift P s of the ground surface temperature can be determined on the basis of the balance of average daily fluxes (20). For the conductivity flux, it should be substituted in accordance with the Fourier equation: The phase shift and amplitude are determined by the formulas: A s = p 3 1 + p 3 (p 1 cos P s + p 2 sin P s ), where: The derivation of dependencies (28) and (29) was discussed in [13,29].

Examples of Time Courses of Heat Fluxes on the Ground Surface and Temperatures
The values presented in Table 1 were taken for calculations. From dependencies (23), (28), and (29) the following parameters of Equation (2) were obtained: T sm = 10.67 • C, A s = 13.88 K, P s = 0.202 rad. The coefficient r a was calculated according to Formula (7), for which the value from Table 1 was used as the wind speed. The time courses of temperatures are shown in Figure 4, while the courses of individual heat fluxes are presented in Figure 5.   For most of the year, the temperature of the ground surface is higher than the air temperature. Moreover, the temperature of the ground surface has a greater amplitude than the air temperature. The temperature of the sky (daily average) is lower than the air temperature and the temperature of the ground surface. The ground surface temperature gradient is negative in the summer months, and positive in the winter months.
The solar radiation flux is always directed towards the ground surface. On the other hand, the long-wave and evaporative (daily average) radiation fluxes are directed in the opposite direction. The average value of the long-wave radiation flux is about 50 W/m 2 . The convective and conductive fluxes change their direction twice a year; the convective

The Heat Flux Extracted from the Ground
In order to simulate the heat transfer in the ground with an installed heat exchanger, the time dependence of the heat flux taken from the ground should be known. A simplified model was adopted in which it was assumed that:

•
Heat is extracted from the ground only for space heating purposes.

•
Outside the heating season, no heat is supplied to the ground.

•
The average daily value of the heat flux extracted from the ground qGHE is assigned to each day of the year.

•
The maximum value of the heat flux extracted from the ground is qmax.

•
The number of days in the year hd for which heat is extracted from the ground for heating purposes is known; for the remaining days, qGHE = 0.
With such assumptions, the heat flux extracted from the ground, qGHE, varies with time according to the relationship [29]: where:

The Heat Flux Extracted from the Ground
In order to simulate the heat transfer in the ground with an installed heat exchanger, the time dependence of the heat flux taken from the ground should be known. A simplified model was adopted in which it was assumed that:

•
Heat is extracted from the ground only for space heating purposes.

•
Outside the heating season, no heat is supplied to the ground.

•
The average daily value of the heat flux extracted from the ground q GHE is assigned to each day of the year.

•
The maximum value of the heat flux extracted from the ground is q max .

•
The number of days in the year h d for which heat is extracted from the ground for heating purposes is known; for the remaining days, q GHE = 0.
With such assumptions, the heat flux extracted from the ground, q GHE , varies with time according to the relationship [29]: where: Figure 6 shows the annual course of the heat flux extracted from the ground for different values of h d , in accordance with dependencies (33)-(35). Integrating the dependence of q GHE on t, in the range corresponding to the heating season, allows one to determine the amount of heat Q extracted from 1 m 2 of the ground during the year.

Application of the Ring Heat Source Model
It was assumed that the thermal power of each ring as a heat source is the same. This power is determined by the relationship: where A j is the ground area per ring. The ground temperature response at time t at the point P with coordinates x, y, and z, caused by the action of N ring ring heat sources, was determined according to the equation [30,31]: where I 0 is the modified Bessel function of the first kind of order zero, τ is the integration variable, E real,j and E virt,j are related to real and virtual heat sources as follows: where r j is the distance of the point P from the centers of the individual rings with the coordinates x 0j , y 0j :

Application of the Ring Heat Source Model
It was assumed that the thermal power of each ring as a heat power is determined by the relationship: where Aj is the ground area per ring. The ground temperature r point P with coordinates x, y, and z, caused by the action of Nrin

Calculations and Simulations of Heat Transfer in a Slinky-Coil Ground Heat Exchanger
The system of 63 rings arranged in seven rows was considered in the calculations and simulations. Figure 7 shows a diagram of the arrangement of the rings in the analyzed system.
Based on dependence (37), the ground temperature responses were determined for different values of the location variables and for three full calendar years of the heat exchanger operation. The 1st of F was assumed as the start of the simulation. Then, the values of the undisturbed ground temperature were generated depending on the position and time according to Equation (2). As parameters of this equation, the values determined in Section 2.5 were used. Finally, the ground temperature was determined by summing the values of ϑ and T 0 , in accordance with relationship (1). The geometric parameters and other data necessary for the calculations are presented in Table 2. Based on dependence (37), the ground temperature r different values of the location variables and for three f exchanger operation. The 1st of F was assumed as the sta values of the undisturbed ground temperature were gener and time according to Equation (2). As parameters of this e in Section 2.5 were used. Finally, the ground temperature the values of  and T0, in accordance with relationship and other data necessary for the calculations are presented Table 2. Numerical values used in calculations of the ground tem slinky-coil exchanger.

Quantity
Value Quantit   Figure 8a refers to the temperature response profiles along the x-axis where the centers of the rings are located. The simulation was performed for the 772nd day of the process (11th ofFebruary). The profiles refer to different distances from the ground surface z. For the depth z = 1.5 m, where the rings are placed, the fluctuations in temperature responses are very distinct. They are caused by the changing distance from the heat sources (pipe walls) when moving along the x axis. The values are always negative, which results from the heat extraction from the ground. Edge effects are visible in the form of decreasing absolute values as one moves away from the center of the system (x = 0). It is worth noting that changes in temperature responses are not symmetrical with respect to the installation level of the exchanger; for z = 0.5 m, the absolute values are smaller than for z = 2.5 m. It results from the short distance from the ground surface treated as a semiinfinite body.

Overview of the Designated Ground Temperature Fields
The undisturbed ground temperature at a given depth depends only on time, and therefore the lines in Figure 8b are horizontal. The undisturbed ground temperature  and 42nd days of the calendar year. For the day under consideration, the undisturbed ground temperature rises with distance from the ground surface.
The ground temperature profiles along the x-axis are shown in Figure 8c. The temperatures shown in this figure are the sum of the values of  and T0 on the vertical axes of Figure 8a,b. For the given depth z, the temperature changes as a function of location (x) have the same shape as the changes of  shown in Figure 8a.  Figure 9a shows the time course of the ground temperature responses resulting from the heat extraction from the ground. For the 3-year exploitation period of the exchanger, the ground temperatures were determined at points lying at different depths and The undisturbed ground temperature at a given depth depends only on time, and therefore the lines in Figure 8b are horizontal. The undisturbed ground temperature courses are the same in subsequent years, and therefore do not differ for the 772nd, 407th, and 42nd days of the calendar year. For the day under consideration, the undisturbed ground temperature rises with distance from the ground surface.
The ground temperature profiles along the x-axis are shown in Figure 8c. The temperatures shown in this figure are the sum of the values of ϑ and T 0 on the vertical axes of Figure 8a,b. For the given depth z, the temperature changes as a function of location (x) have the same shape as the changes of ϑ shown in Figure 8a. in Figure 9b. In this case, the curves for individual annual periods are identical. The amplitudes of individual curves decrease with depth, which is significant for ground temperatures under natural conditions. The time courses of ground temperature are shown in Figure 9c. They result from summing the courses of  and T0. The highest temperatures refer to the smallest depths (z = 1 m), and the lowest ones refer to the depth at which the exchanger is installed (z = 1.5 m).
(a) (b) (c) Figure 10. Ground temperature profiles along the z axis: (a) temperature response caused by the heat source, (b) undisturbed ground temperature, (c) ground temperature. Figure 10a shows the ground temperature response profiles resulting from the heat extraction from the ground for different process times. The profiles were determined along the depth coordinate z, i.e., for x = 0 and y = 0. The range of values for t covers the  Figure 9a shows the time course of the ground temperature responses resulting from the heat extraction from the ground. For the 3-year exploitation period of the exchanger, the ground temperatures were determined at points lying at different depths and corresponding to the values of the x = 0; y = 0 coordinates. The falling part of the curves corresponds to the extraction of heat from the ground, and the rising part refers to thermal regeneration.
The first year of operation is characterized by a slightly different course than the others. During the period of heat extraction from the ground, the highest absolute value occurs near the depth z = 1.5 m, at which the heat exchanger is installed.
The time course of the undisturbed ground temperature at different depths is shown in Figure 9b. In this case, the curves for individual annual periods are identical. The amplitudes of individual curves decrease with depth, which is significant for ground temperatures under natural conditions.
The time courses of ground temperature are shown in Figure 9c. They result from summing the courses of ϑ and T 0 . The highest temperatures refer to the smallest depths (z = 1 m), and the lowest ones refer to the depth at which the exchanger is installed (z = 1.5 m). Figure 10a shows the ground temperature response profiles resulting from the heat extraction from the ground for different process times. The profiles were determined along the depth coordinate z, i.e., for x = 0 and y = 0. The range of values for t covers the heating season, and the value for t = 772 days (11 February) corresponds to the greatest demand for space heating. For z = 0 it is ϑ = 0, which corresponds to the isothermal ground surface conditions. During the heating season, the highest absolute values occur at the level of the exchanger installation (z = 1.5 m).
The undisturbed ground temperature profile is shown in Figure 10b. The slope of the curves for z = 0 determines in which direction the heat is transported across the ground's surface under natural condition; for 680 < t < 800 it is the direction from the ground to the surroundings, and for 800 < t < 880-from the surroundings to the ground. Figure 10c shows the ground temperature profiles with the exchanger installed. This figure shows the minimum ground temperature of approx. −2.5 • C. What is more, the obtained courses confirm the observations presented in the work by Pauli et al. [6] concerning the ground temperature at the level of the exchanger installation-despite the heat consumption, this temperature is generally not lower than the ground surface temperature; moreover, its value does not fall significantly below 0 • C.
The analysis shows that the ground temperature is constantly changing. These changes are caused by two processes:

•
Heat consumption by the working medium flowing through the ground heat exchanger.

•
Heat transfer between the ground surface and its surroundings.
In the considered case, the first of the processes occurs periodically ( Figure 6). Ground temperature changes caused by each of the processes separately can be determined by the analytical relationships. The global change in the ground temperature results from a simple summation (Equation (1)), according to the superposition principle. The conducted simulations make it possible to assess the contribution of individual processes to the total change of the ground temperature. The amplitudes of temperature changes shown in Figure 9b are greater than the amplitudes of temperature changes shown in Figure 9a. Thus, for the conditions for which the simulations were conducted, it can be seen that the heat transfer between the ground surface and the surroundings has the dominant influence on the ground temperature. A similar conclusion can be drawn from the comparison of the temperature profiles shown in Figure 10a,b.

Ground Temperature Calculation Methodology
The calculation scheme of the ground temperature with the installed slinky-coil exchanger is presented below. The ground temperature T at the point x, y, z is determined after time t.

I.
The undisturbed ground temperature T 0 should be calculated according to Formula (2) as follows: • Using the data from Table 1, calculate r 1 , r 2 , r 3 , and then T sm (according to Equations (23)-(26)); • Calculate successively p 1 , p 2 , p 3 , and then P s and A s (according to Equations (28)-(32)); • Substitute the calculated values of T sm , A s , and P s into Formula (2).

II.
In order to determine the temperature response ϑ, one must know values presented in Table 2. The following steps should be taken:

III.
Calculate the ground temperature T from Formula (1).

Comparison with Experimental Results
The computational values of temperature responses are the basis for determining the temperature field for the ground with an installed heat exchanger. The values determined from Equation (37) were compared with the experimental data presented by Wu et al. [5]. The research concerned a slinky-coil exchanger, placed at a depth of 1.14 m, which worked continuously for 52 days with a heat extraction rate equal to 4220 W. The rings were placed in several rows, significantly spaced from each other (4 m).
The calculations of ground temperature responses were conducted for the point corresponding to the location of the temperature measurement sensor (at the wall of one of the rings). The calculations took into account the influence of seven adjacent rings on both sides of the ring under consideration (N ring = 15).
In order to eliminate the influence of daily changes of ground temperature, only the results of measurements conducted at a depth below 0.5 m were used for comparison. The results of measurements of ground temperature at various levels above the exchanger T and at a large distance from the exchanger T 0 are presented in Table 3. Table 3. The results of the experiments presented in [5], and the values ϑ used in Figure 11.

III.
Calculate the ground temperature T from Formula (1).

Comparison with Experimental Results
The computational values of temperature responses are the basis for d temperature field for the ground with an installed heat exchanger. The valu from Equation (37) were compared with the experimental data presented b The research concerned a slinky-coil exchanger, placed at a depth of 1.14 m, continuously for 52 days with a heat extraction rate equal to 4220 W. The rin in several rows, significantly spaced from each other (4 m).
The calculations of ground temperature responses were conducted corresponding to the location of the temperature measurement sensor (at t of the rings). The calculations took into account the influence of seven adj both sides of the ring under consideration (Nring = 15).
In order to eliminate the influence of daily changes of ground temper results of measurements conducted at a depth below 0.5 m were used for co results of measurements of ground temperature at various levels above th and at a large distance from the exchanger T0 are presented in Table 3.  Figure 11. Comparison of experimental and computational results. Figure 11 shows a comparison of the values from Table 3 and those Equation (37). The compatibility of calculated and experimental values is g

Conclusions
As a result of the theoretical analysis, calculations and simulations of h the ground under natural conditions and under conditions with an  Figure 11. Comparison of experimental and computational results. Figure 11 shows a comparison of the values from Table 3 and those obtained from Equation (37). The compatibility of calculated and experimental values is good.

Conclusions
As a result of the theoretical analysis, calculations and simulations of heat transfer in the ground under natural conditions and under conditions with an installed heat exchanger, it was found that:

•
Parameters of Equation (2) for determining the undisturbed ground temperature-T sm , A s , and P s -can be calculated from dependencies (23), (28). and (29) The ring heat source model was used to find the 3D ground temperature field near a slinky-coil exchanger. In calculations, dependence (37) was used to determine the temperature response caused by the heat sources and Equation (2)-to find the undisturbed ground temperature. The parameters of Equation (2) were specified on the basis of dependencies (23), (28) and (29).

•
The influence of spatial variables x and z and time t on the ground temperature was investigated in various combinations. With the use of the superposition principle, it was also shown graphically how the ground temperature is influenced by heat sources and variable ground surface temperature.

•
The results of calculations lead to conclusions consistent with those presented in the work by Pauli et al. [6]. Quantitative comparisons were also made by comparing the results of calculations with the results of experiments presented in the study by Wu et al. [5].

•
The model used to calculate the temperature field can be used to study the influence of the process parameters on the temperature field near the exchanger pipes. This enables finding such process parameters that prevent the occurrence of dangerously low ground temperature values. Informed Consent Statement: Not applicable.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
A amplitude of temperature, K A j ground surface area per 1 ring, m 2 A s annual amplitude of temperature fluctuations on the ground surface, K C 1 , C 2 constants C LW constants defined by Equation (12)