Heat Transfer Model to Predict Temperature Distribution in the Ground

Abstract

Knowledge of the temperature of the ground in time and space as well as its thermal properties gives basic information about physical phenomena concerning the transfer and accumulation of heat in the ground. It can be also used for evaluation of the heating possibilities of heat pumps; to proper design the size of the ground exchangers and the depth, at which they should be installed. For this purpose, a mathematical model based on the heat balance equation on the ground surface was developed. The basis of the model is the Carslaw-Jaeger equation regarding the temperature profile in the ground. The model was verified using experimental results for two different locations (different climatic conditions, moderate and arid climate)—the standard deviation is equal 0.62 K and 0.92 K, respectively. In this work, the impact of several parameters on the ground temperature profiles and thermal fluxes was determined. It was found that among the examined parameters the amplitude of the daily average solar radiation flux strongly effects on the total amount of heat transferred between the ground and the environment during the year, wherein the other parameters have a negligible effect.

1. Introduction

The lower heat source for heat pumps is usually air or ground. The ground is an advantageous heat source, due to a much more stable temperature and a high specific heat. Ground temperature and its fluctuations (apart from other physical quantities of the ground) is a factor that determines the course of physical, chemical and biological processes in the ground. Knowledge of the temperature of the ground in time and space as well as its thermal properties gives basic information about physical phenomena concerning the transfer and accumulation of heat in the ground. This is particularly important in the design, modelling and exploitation of ground heat exchangers as well as underground buildings and other structures related to heat transport in the ground (for example, pipes transporting heat carriers).

Experimental determination of the temperature distribution in the ground should be carried out continuously for a sufficiently long period of time, so that the measured values are not random, which affects the high cost of measurements. However, the results of the ground temperature measurements at various depths and in different places around the world are published in the literature [1,2,3,4,5]. Popiel and Wojtkowiak [1] presented the results of the temperature distributions of the ground monitored during 10 years in the region of Poznan City. The ground temperature was measured at a depth from 0 to 6.9 m (car park) and from 0 to 17.3 m (lawn). The results of extensive studies of temporal temperature changes of the ground at various depths under different climatic conditions in the USA have been presented by Kusuda and Achenbach [4] and Neuberger and Adamovský [5] carried out the measurements during three heating periods in Prague. To be useful, the measurement results must relate to depth >1 m, because for shallow layers of the ground diurnal changes impose on the effect of seasonal temperature changes.

Thus, the development a mathematical model for heat transfer in the ground and the estimation the subsurface temperature distribution is a key issue. Therefore, the boundary condition on the surface of the ground should be defined carefully. The annual variation of the ground temperature at different depths can be determined using a sinusoidal temperature variation as a first kind boundary condition [4,6,7]. In numerical studies, Piechowsky [8] included the heat flux associated with mass transfer in order to take into account the effects of the soil moisture content and migration. Kupiec et al. [9] considered only the convective heat flux between air and ground. Jaszczur et al. [10] studied the impact of individual model elements on the temperature of the ground. It has been found that the simplest models and the most complex model result in a similar temperature variation over the simulation period, but only at a low depth. Also, Bortoloni et al. [11,12] analysed the effect of the first, second and third kind boundary conditions (Dirichlet, Neumann and Robin boundary condition, respectively) imposed at the ground surface in modelling heat transfer process in the ground. They improved the convective component by introducing the aerodynamic resistance of the canopy layer. Most of other works are based on the heat fluxes balance equation on the ground surface, i.e., short-wave and long-wave radiation fluxes, convective heat flux, evaporative heat flux and conductive flux [13,14,15,16,17], however some of these works assume that the long-wave radiation flux has a constant value or show the results of temperature measurements or heat transfer rates when a ground heat exchanger is installed in the soil. For example, Nam et al. [17] developed a numerical model that combines a heat transport model with groundwater flow using heat fluxes balance on the ground surface and a heat exchanger model with a U-tube shape.

Although the influence of the type of boundary condition in the literature has been thoroughly analysed, there are still no reports on the impact of various parameters on the distribution of temperature in the ground, or even on the amount of heat conducted through the ground under natural conditions (except for the water content in the soil and its thermal conductivity). This applies particularly to those parameters that are difficult to determine accurately, for example the heat transfer coefficient between the surface of the ground and the surroundings.

The aim of this work is to present a universal, simple mathematical model to predict temperature distribution in the ground and to analyse the size of thermal fluxes occurring on the surface of the ground (primarily in terms of variability of these fluxes over time). In the developed mathematical model, it was taken into account that the long-wave radiation flux is also variable over time. The impact of coefficients: amplitude of annual solar radiation flux, heat transfer coefficient, emissivity of the ground surface and evaporation rate coefficient on the ground temperature as well as on the thermal fluxes, especially on the total amount of heat transferred between the ground and the environment during the year resulted from the conductive heat flux in the ground, is determined.

The presented simulation results relate to average climatic conditions occurring in Cracow. In addition, based on data from the literature, the compatibility of the Carslaw-Jaeger equation with the results of ground temperature measurements at different depths in different climatic conditions (Lemont, USA and Zarqa, Jordan) was evaluated.

2. Heat and Mass Transfer Model

2.1. Heat Conduction in the Ground

Heat transfer in the ground occurs mainly as a result of heat conduction. Heat conduction equation is given by:

$$\frac{{\partial }^{2}T}{\partial x^2}=\frac{1}{a}\cdot \frac{\partial T}{\partial t}$$

(1)

where: T—temperature of the ground (°C); a—thermal diffusivity of the ground (m²/s); x—position coordinate (m); t—time (s).

In numerical calculations the ground is treated as a plate with a finite thickness. On the upper surface of the plate (the ground surface) periodically variable heat fluxes occur. These fluxes are variable both in the annual and daily cycles. Changes resulting from daily cycles occur at depths of less than 1 m. The boundary condition on the surface of the ground is given by:

$$x=0:\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{q}_{cond}=H-LW+S-EV$$

(2)

where: q_cond—conductive heat flux (W/m²); H—convective heat flux (W/m²); S—solar radiation heat flux (W/m²); LW—long-wave radiation heat flux (W/m²); EV—evaporative heat flux (W/m²).

Geothermal flux has, in general, a slight impact on the ground temperature profile and will not be taken into account in this work. Therefore, the second boundary condition is as follow:

$$x=x_{inf}:\hspace{1em}\hspace{1em}\hspace{1em}\frac{dT}{dx}=0$$

(3)

where x_inf is the depth at which the temperature of the ground is independent on position coordinate (m). Due to neglecting the geothermal flux the undisturbed ground temperature and the average temperature of the ground surface are the same.

The solution of the heat transfer Equation (1) with boundary conditions (for x = 0: T_s = f(t) and for x_inf → ∞: T → T_s) was provided by Carslaw and Jaeger [18]. The daily average temperature of the ground T is the following function of position coordinate x and time t:

$$T(x,t)=T_{sm}-A_s\cdot \exp \left(-\frac{x}{L}\right)\cdot \cos \left(\omega t-P_s-\frac{x}{L}\right)$$

(4)

The parameters of the above relationship are: T_sm—annual average temperature of the surface of the ground (°C); A_s—amplitude of daily average temperature of the ground surface (K); P_s—phase angle (rad); L—damping depth (m) defined as:

$$L=\sqrt{\frac{2a}{\omega }}$$

(5)

where: ω—frequency of temperature fluctuations. For phenomena occurring in the annual cycle, the frequency is equal to ω = 2π/365 days⁻¹.

Damping depth, defined by Equation (5), is a constant characterizing the decrease in amplitude with an increase in distance from the ground surface. The ground can be treated as a system consisting of a subsurface layer, in which there are interactions related to changing weather conditions, and a deeper layer in which these impacts do not occur. The thickness of subsurface layer depends on the thermal diffusivity of the ground. For low values of thermal diffusivity, the subsurface layer has a small thickness, but when the thermal diffusivity of the ground is high, the stabilization of the ground temperature occurs at larger depths [15].

2.2. Heat Balance on the Ground Surface

Based on a theoretical analysis of thermal fluxes occurring on the surface of the ground, the temperature of the ground surface can be linked to the air temperature. Both of these quantities are time-dependent and defined by analytical relationships containing the cosine function. The parameters of these relationships are: annual average temperatures, annual fluctuation amplitudes and phase angles. These parameters are easily accessible for the air. In order to determine the parameters for the ground surface with known parameters for air, the heat balance on the ground surface should be considered, and the dependence for determining individual heat fluxes should be formulated. Based on this heat balance, a mathematical model can be developed.

On the surface of the ground are the following, presented in Figure 1, heat fluxes: conductive, convective and thermal radiation fluxes. Moreover, on the surface of the ground there is moisture evaporation accompanied by the phase change heat transfer. The heat balance on the surface of the ground has the form as the boundary condition for the ground surface.

2.2.1. Convective Heat Flux

Convective heat flux is determined by the heat transfer equation:

$$H=h\left(T_a-T_s\right)$$

(6)

where: h—convective heat transfer coefficient (W/(m²K)); T_a—temperature of the air (°C); T_s—temperature of the ground surface (°C).

The heat transfer coefficient h is difficult to determine precisely. There are many different empirical relationships to determine this coefficient; they take into account that the coefficient h is a function of wind velocity. In this study the McAdams formula [19] is used:

$$\begin{array}{l}h=5.7+3.8\hspace{0.17em}u\hspace{1em}\hspace{1em}\mathrm{for}\hspace{1em}\hspace{1em}u<4.88\hspace{0.17em}\mathrm{m}/\mathrm{s}\\ h=7.2\hspace{0.17em}{u}^{0.78}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{1em}\hspace{1em}u\ge 4.88\hspace{0.17em}\mathrm{m}/\mathrm{s}\end{array}$$

(7)

where: u—wind velocity (m/s).

Daily average air temperature changes in the annual cycle as follows:

$$T_a=T_{am}-A_a\cos \left(\omega t-P_a\right)$$

(8)

where: T_am—annual average of the air temperature (°C); A_a—amplitude of air temperature (K); P_a—phase angle (rad).

An example temporal course of daily average air temperature is shown in Figure 2 (blue line).

The daily average temperature of the surface of the ground is periodically variable, like the air temperature:

$$T_s=T_{sm}-A_s\cos \left(\omega t-P_s\right)$$

(9)

where: T_sm—annual average of the ground surface temperature (°C): A_a—amplitude of the ground surface temperature (K); P_s—phase angle (rad). The parameters of Equation (9) T_sm, P_s, A_s are different from the corresponding parameters T_am, P_a and A_a for the air.

2.2.2. Short-Wave Radiation Heat Flux

The short-wave flux comes from solar radiation. For balance calculations it is convenient to use a daily average solar radiation flux absorber by the ground S. The solar radiation flux S changes in the annual cycle as follows:

$$S=S_m-A_{sol}\cos \left(\omega t-P_{sol}\right)$$

(10)

where S_m is the annual solar radiation flux absorbed by the ground (W/m²), A_sol is the amplitude of the daily average of this radiation flux (W/m²), whereas P_sol is the phase angle (rad).

2.2.3. Long-Wave Radiation Heat Flux

The ground radiates heat energy since the temperature of the ground surface T_s is higher than the sky temperature T_sky. The daily average net flux of long-wave radiation LW equals:

$$LW=\epsilon \sigma \left(T_s^4-T_{sky}^4\right)\cong \epsilon C_{LW}\left(T_s-T_{sky}\right)$$

(11)

where σ = 5.67 × 10⁻⁸ W/(m²K⁴) is the Stefan-Boltzmann constant, ε is the emissivity of the surface of the ground and C_LW is a constant equal to 4.83 W/(m²K). The second part of Equation (11) is a linearized form, convenient for modelling. The approximation concerns small absolute differences of T_s and T_sky temperatures.

The sky temperature T_sky depends on the air temperature and its humidity; moreover, it is variable during the day and night. The possibilities to determine it include the use of the empirical formula [15]:

$$T_{sky}=T_a\hspace{0.17em}{\left[0.711+0.0056 \hspace{0.17em}{T}_{dp}+0.000073\hspace{0.17em}{T}_{dp}^2+0.013\cdot \cos \left(15t\right)\right]}^{\hspace{0.17em}1/4}$$

(12)

where T_dp is a dew point temperature (°C), T_sky and T_a are expressed in Kelvins, while t is time measured since midnight expressed in hours. The temporary variability of the sky temperature is determined by the dependence:

$$T_{sky}=T_{sky,m}-A_{sky}\cos \left(\omega t-P_a\right)$$

(13)

where: T_sky,m—annual average temperature of the sky (°C); A_sky—amplitude of the sky temperature (K). The values of the phase angle of the temperature of air and sky are similar. An exemplary temporal course of daily average sky temperature is shown in Figure 2 (black line).

2.2.4. Evaporative Heat Flux

The stream of moisture is caused by a difference in the partial pressures of water vapour on the surface of the ground and in the bulk of the gaseous phase. This mass flux is connected with the heat flux transfer required for the evaporation of water—the heat being lost by the ground [15]. Assuming a linear form of the relationship of saturated vapour pressure and temperature: P_sat = a_pT + b_p (a_p = 103 Pa/K, b_p = 609 Pa [13]), the evaporative heat flux amounts to:

$$EV\cong C_{EV}fh\hspace{0.17em}\left[\left(a_p{T}_s+b_p\right)-RH\cdot \hspace{0.17em}\left(a_p{T}_a+b_p\right)\right]$$

(14)

where C_EV is treated as a constant and equals 0.0168 °K/Pa [13], f is evaporation rate coefficient (-) and RH is relative humidity of the air (-).

2.3. Determination of T_sm, P_s, A_s

The basis for calculation of thermal fluxes on the ground surface is the surface temperature determined by Equation (9). The parameters of this equation are computed based on equations resulting from the heat balance Equation (2) and the Carslaw–Jaeger Equation (4).

Yearly average values of fluxes H, S, EV and LW result in the value of the average annual temperature of the ground surface:

$$T_{sm}=\frac{\epsilon C_{LW}{T}_{sky,m}+h p_r{T}_{am}+S_m-C_{EV}fh\hspace{0.17em}{b}_p \left(1-RH\right)}{h{p}_e+\epsilon C_{LW}}$$

(15)

where:

$$p_e=1+C_{EV}f{a}_p$$

(16)

$$p_r=1+C_{EV}f{a}_p\cdot RH$$

(17)

In order to determine the P_s value, the nonlinear algebraic equation should be solved:

$$\left[\sin \left(P_s\right)+\left(p_1+1\right)\cdot \cos \left(P_s\right)\right]\cdot \left[p_2\cdot \cos \left(P_s-P_a\right)+\cos \left(P_s-P_{sol}\right)\right]-p_3=0$$

(18)

Amplitude A_s can be determined from the formula:

$$A_s=p_4\cdot \frac{p_2\cos \left(P_s-P_a\right)+\cos \left(P_s-P_{sol}\right)}{p_1+1}$$

(19)

wherein the constants p₁–p₄ are:

$$p_1=\frac{L}{k}\left(p_{e}h+\epsilon C_{LW}\right)$$

(20)

$$p_2=\frac{p_{r}h{A}_a+\epsilon C_{LW}{A}_{sky}}{A_{sol}}$$

(21)

$$p_3=\left[p_2\cdot \cos \left(P_a\right)+\cos \left(P_{sol}\right)\right]\cdot \left(1+p_1\right)$$

(22)

$$p_4=\frac{A_{sol}L}{k}$$

(23)

3. Experimental Verification of the Model

In order to verify the presented mathematical model, based on the Carslaw-Jaeger Equation (4), temperature profiles generated computationally were compared with the measurements presented by Kusuda and Anechbach [4] and Al-Hinti et al. [2].

In both cases nonlinear regression using the Solver application (Excel) was utilized to develop the results—the parameters of Equation (4). The sum of squares (SS) of differences between experimental and calculated temperatures according to Equation (4) was minimized:

$$SS={\displaystyle \sum _{i=1}^n{\left(T_{calculated}-T_{experimental}\right)}^2}$$

(24)

where n is the number of measurements used to develop results.

The standard deviation between experimental and computational temperatures (based on the values of parameters determined by nonlinear regression) is equal to:

$$\sigma =\sqrt{\frac{SS}{n-m}}$$

(25)

where m—the number of determined parameters.

• Lemont, USA

The measurements were carried out under moderate climate conditions (Lemont, Illinois, USA). Annual average of ambient air temperature in Lemont equals to 10.0 °C, and yearly average sunshine duration is 2508 h [20]. The temperature of the ground to a depth of 8.84 m was measured and results were presented in [4]. The results of measurements carried out for whole the year were utilized (n = 84).

The following values of parameters of Equation (4) have been obtained: T_sm = 11.2 °C, A_s =13.1 K, P_s = 0.664 rad, L = 2.44 m. In Figure 3, the comparison of experimental values and values calculated from Equation (4) for the values of the parameters presented above is shown. Temperature profiles for selected months of the year: January, April, July and October are presented. The symbols represent the experimental values read out from the drawings shown in [4]. The good compliance of experimental and computational values confirms that Equation (4) correctly describes the temperature distribution in the ground and its variability over time. The ground temperature stabilizes throughout the year at around 11 °C (undisturbed ground temperature) starting from a depth of 10 m below the surface.

The sum of squares equals to SS = 30.28. The standard deviation according to (25) is equal:

$$\sigma =\sqrt{\frac{30.28}{84-4}}=0.62\hspace{0.17em}\mathrm{K}$$

• Zarqa, Jordan

The measurements were conducted at an open field, near the city of Zarqa, Jordan. The climate is arid and characterized by low annual rainfall, and around 3300 h of sunshine per year. Yearly average ambient air temperature equals to 19.2 °C. Thus, climatic conditions are completely different from Lemont, USA.

The measurements were carried out over a one-year period extending from October 2014 to August 2015. The temperature of the ground at five selected depths (0.5, 1.0, 2.0, 5.0, and 10.0 m) as well as the ambient air and the ground surface temperature was measured.

In Figure 4 the comparison of experimental values and values calculated from Equation (4) for the values of the parameters presented above is presented. The results of measurements carried out at different times of the year: January 18th, April 28th, August 6th and October 10th were utilized (n = 49). Also, in this case, the ground temperature obtained by the presented model is in good agreement with the experimental results (Figure 4b). The following values of parameters of Equation (4) have been obtained: T_sm = 21.3 °C, A_s = 11.9 K, P_s = 0.232 rad, L = 1.91 m. The sum of squares equals to SS = 38.08 and the standard deviation according to (25) is equal σ = 0.91 K.

The calculated values of the ground temperature confirm the author’s observations [2] that underground temperature stabilizes throughout the year at around 21°C (calculated undisturbed ground temperature equals to 21.3 °C) starting from a depth of 8 m below the surface, as can be seen in Figure 4a.

4. Climatic Conditions in Cracow

Cracow, Lesser Poland Voivodeship (latitude of 52°08’ N and longitude of 19°92’ E) is located on the lower boundary of the moderate warm climatic zone of the Carpathians as a variation of valley climate (according to Hess, 1969) [21]. Based on data (1981–2010) [22] the annual average air temperature is 8.5 °C, wherein the minimum temperature is in January: from −10.6 to 3.2 °C (the temperature of the coldest and warmest month), and maximum temperatures occur in July: from 15.6 to 21.6 °C. The chart of annual change in ambient air temperature (marked as blue symbols) and sky temperature (marked as grey symbols) is presented in Figure 2.

The sculpture of the terrain undoubtedly influences the anemological conditions in Cracow. The predominant wind direction is western and then southwest. The average wind velocity is approx. 2.5 m/s. Most days with strong wind (>8 m/s) occur in the winter. The annual average of relative humidity of the air equals to 79%. The cloud cover is on average 68% [23].

The annual average total sunshine duration based on data from 1884–2014 is about 1555 h. Due to astronomical reasons, the highest insolation occurs in July (221.5 h), and the smallest (37.5) in December [24].

The monthly average climatic conditions in Cracow are presented in

Table 1. Annual and monthly average climatic conditions in Cracow, Poland.

Month	Jan.	Feb.	Mar.	Apr.	May	June	July	Aug.	Sept.	Oct.	Nov.	Dec.	Year
Temperature (°C)	−2.1	−0.8	3.1	8.7	14.0	16.8	18.8	18.1	13.5	8.7	3.1	−0.9	8.5
Precipitation (mm)	37.5	29.7	40.2	46.4	81.2	86.4	87.9	75.7	62.4	43.2	42.1	39.2	671.9
Wind velocity (m/s)	3.1	2.5	3.0	2.4	2.2	2.5	2.1	1.8	2.3	2.3	2.7	2.9	2.5
Relative humidity (%)	83	82	76	74	71	76	75	78	81	83	87	86	79
Sunshine duration (h)	47.9	64.8	108.4	150.9	203.7	207.0	221.5	203.8	151.6	105.9	52.0	37.5	1555.3

.

5. Results of Calculations

The presented mathematical model contains parameters that are difficult to determine:

• Heat transfer coefficient h, depending on the wind velocity,
• Emissivity of the surface of the ground ε, depending on the type of coverage of this surface,
• Evaporation rate coefficient f, depending on the humidity of the ground (this humidity is influenced by the amount of precipitation and ground permeability for water).

Parametric analysis of this model was carried out and the impact of h, ε and f on the temporal ground temperature courses was determined. The calculations were performed for moderate climate conditions (Cracow, Poland).

Knowing the temporal courses of the ground temperature on its surface, it is possible to determine the temporal courses of thermal fluxes. Knowledge of LW, EV and H fluxes (in addition to the knowledge of the solar radiation flux S, independent of the ground temperature) allows to determine the heat flux associated with the heat transfer to the subsurface q_cond.

• Analysis of Heat Fluxes on the Ground Surface

The evaporation rate coefficient f takes into account that the rate of evaporation of water from the ground surface is lower than the rate of evaporation from the water surface; this factor ranges from 0.1–0.2 for dry soils up to 0.4–0.5 for humid soils [13]. The impact of the evaporation rate coefficient on temporal courses of heat fluxes on the surface of the ground is presented in Figure 5. The value of evaporative heat flux varies significantly for dry and for moist soil, especially during the summer: for dry soil (f = 0.1) the maximum value of EV is equal to 30 W/m², wherein for moist soil (f = 0.4) this value is three times as large as for dry soil. However, differences in the values of H are smaller than EV, and vary from 65 to 110 W/m² (maximum values). However, there is no discrepancy in the values of conductive heat flux for dry and moist soil, q_cond changes slightly (marked as brown symbols).

In Figure 6 temporal courses of heat fluxes on the surface of the ground for different values of wind velocity u are shown. According to Equation (7) as the wind velocity increases, the value of the heat transfer coefficient also increases, which results in lower value of convective heat flux (applies to negative values). The maximum difference between the values of H compared to u = 2 m/s and for u = 4 m/s is equal to 10 W/m² (during the whole year). Changes in the values of the discussed thermal fluxes cause minimal changes in the conductive heat flux values.

Emissivity of the surface of the ground ε, depending on the type of coverage of this surface, adopts different values: for example, for asphalt ε ranges from 0.9 to 0.98, for the ground ε ranges from 0.92 to 0.96, for vegetation ε = 0.95 and for snow ε = 0.83 [25]. The impact of the emissivity of the ground surface on the temporal courses of the heat fluxes is presented in Figure 7. The changes mainly concern the LW flux, from Equation (11). Courses of convective and evaporative fluxes varies slightly. However, even for a large change in the value of emissivity of the ground surface, changes in the values of individual fluxes are not as significant as in the case of u or f. Furthermore, these changes do not affect the change in a value of the conductive heat flux.

• Analysis of Conductive Heat Flux

A heat flux transported to (or from) the interior of the ground q_cond results from the algebraic summation of all thermal fluxes on the surface of the ground (Equation (2)). In spring and summer, the conductive heat flux is directed from the surface to the subsurface, while in autumn and winter, heat from the deeper layers of the ground is transported towards the surface.

In spite of the range in the values of fluxes H and EV for different value of f (Figure 5) and u (Figure 6) value of conductive heat flux is changed slightly over time. However, as the wind speed increases, the value of the conductive heat flux varies slightly (in comparison with values of other fluxes): q_cond is equal to 6.5, 7.0, 9.0 W/m² for u = 1, 3, 5 m/s, respectively (calculated for f = 0.3). More heat is conducted deep into the ground for dry soil (Figure 8), because the evaporation rate coefficient is lower than for moist soil, which in turn leads to a reduction in the heat flux associated with evaporation of moisture (provided that k = constant). It results in the total amount of heat transferred between the ground and the environment during the year Q. As can be seen from Figure 6, Q for dry soil is equal to 34.56 kWh/m², wherein for moist soil is equal to 31.77 kWh/m² (calculated for h = 13.3 W/m²K). Furthermore, with the higher value of amplitude of the daily average of this radiation flux (A_sol) the value of Q increases: for A_sol = 188 W/m² the amount of transferred heat increases by 46% compared to A_sol = 101 W/m² for moist soil.

The impact of h, f and ε on the total amount of heat transferred between the ground and the environment during the year is shown in Figure 9. With an increase in the value of h, the amount of heat transferred through the ground surfaces decreases, regardless of the value of the ground surface emissivity or the evaporation rate coefficient (for k = constant). This amount of heat is 20% greater for dry soil than for moist soil. A significant change in the value of the emissivity of the ground surface does not cause a large difference in the amount of transferred heat. Furthermore, these differences decrease with the increase of the heat transfer coefficient and the evaporation rate coefficient. In the case of moist soil for wind velocity above 2.5 m/s (i.e., h > 15 W/m²K), the emissivity in the 0.5–1.0 range does not affect the amount of heat transferred inside the ground.

• Direction of Heat Transfer

The directions of the fluxes q_cond and H (Figure 10) are variable over time. If T_a > T_s, then H > 0 and convective heat flux is directed from the air to the ground surface (Figure 10a,b). If (dT/dx)_x=0 < 0, then q_cond < 0 and conductive heat flux is directed from the ground surface to its interior (Figure 10b,d). Temperature gradient for the ground surface can be determined by the differentiation of Equation (4). Hence:

$${\left(\frac{\partial T}{\partial x}\right)}_{x=0}=-\frac{A_s}{L}\left[\sin \left(\omega t-P_s\right)-\cos \left(\omega t-P_s\right)\right]$$

(26)

Because the ground is in a cyclic steady state, the yearly average value of the conductive flux is zero.

In Figure 11 the temporal courses of fluxes q_cond and H are presented. There are the following theoretically directions of heat fluxes throughout the year: both H and q_cond are positive or negative, and both of these fluxes are opposite signs. However, in reality, during the year there is no case for which both H > 0 and q_cond > 0 occurs (Figure 10a). Climatic conditions directly affect the date of changes in the direction of thermal fluxes q_cond and H. For annual average climatic conditions in Cracow parameters of Equation (4) were determined: T_sm = 10.9 °C, A_s = 13.8 K, P_s = 0.166 rad (calculated for a = 0.6·10⁻⁶ m²/s). Under these conditions, convective heat flux is negative during ¾ of the year (February 3rd to November 8th) due to the relationship T_a < T_s—the heat flux is directed from the ground surface to the air. But the ground cools down in shorter period (half a year): February 24th to August 25th, due to the fact that (dT/dx)_x=0 < 0.

• Temperature Distribution in the Ground

The change in the values of coefficients h, f and ε affects not only the values of individual thermal fluxes occurring on the surface of the ground, but also the temperature distribution in the ground. The ground temperature profile also depends on the thermal diffusivity of the ground and the associated damping depth as well as on the atmospheric conditions. On the basis of the heat balance on the ground surface, it is possible to determine the parameters of Equation (4), which will allow determination the temperature distribution in the ground with known thermal parameters.

Based on climatic data for Cracow-Balice (1981–2010) [21] the parameters of Equations (8), (9) and (12) were obtained: T_am = 8.5 °C, A_a = 10.6 K, P_a = 0.270 rad, T_skym = −0.3 °C A_sky = 11.6 K, S_m = 119 W/m², A_sol = 101 W/m², P_sol = −0.153 rad.

In Figure 12 the impact of u, f and ε on the temperature profiles in the ground are shown. Example temperature profiles for the 49th (green lines), 288th (red lines) and 349th (dashed lines) day of the year which is February 18th, October 15th and December 15th respectively. The calculations were carried out for a = 0.6·10⁻⁶ m²/s. As can be seen from Figure 12 the change in the value of the emissivity of the ground surface does not affect the q_cond value (Figure 7) but also the course of the ground temperature profile, which remains practically unchanged. In the subsurface layers of the ground (x < 3 m) the wind velocity also does not cause a significant change in the ground temperature. However, for larger depths, the ground temperature varies even by 1.2 °C, even at a depth of 16 m, which affects the value of the undisturbed temperature of the ground. The ground at any depth has a higher temperature for lower wind velocity as well as for lower value of evaporative rate coefficient, which applies to dry soils.

In Figure 13, temporal changes of ground temperature at a depth of 0 m, 1 m and 5 m for different values of wind velocity are presented. As can be seen, the wind velocity significantly affects the temporal course of temperature in the ground at any depth; for x = 5 m the difference in the ground temperature is up to 2 °C during the whole year. For the ground surface and a depth of 1 m, significant differences in the ground temperature occur only in the summer—during fall and winter the increase in the wind speed does not cause a significant difference in the temperature of the ground. Moreover, at a depth higher than 1 m, the ground temperature is positive throughout the year.

6. Conclusions

The presented, verified model has been used to simulate the temperature distribution in the ground under moderate climatic conditions. This model can be also used to analyse the influence of various climatic conditions on the efficiency of the subsurface as a heat source or sink for ground coupled heat pumps, to predict the ground temperature profile as function of time and depth at any place (as long as climatic conditions are known) as well as to analyse the size of thermal fluxes occurring on the surface of the ground.

Heat fluxes H, EV and LW are interrelated. The change of one of these fluxes affects the value of the other, which in consequence, when the S value is determined at a given time, causes the q_cond value to remain unchanged. Solar radiation flux (S) does not have this feature. Although for large S values, the values of opposite directed fluxes also increase, but not enough to keep q_cond unchanged.

In addition, even if the values of coefficients h, f and ε cannot be accurately defined, they have minimal effect on the conductive heat flux between the environment and the ground (q_cond). Amplitude of the daily average solar radiation flux has a large influence on conductive flux, but this value (A_sol) is generally quite well known.

The temperature profiles in the ground at different depths determined with the use of the presented simple mathematical model are consistent with the results of the measurements shown in the literature for different climatic conditions—moderate and arid climate. The parameters of the Carslaw-Jaeger equation: T_sm, A_s, P_s and L was determined using nonlinear regression. The good compliance of experimental and computational values confirms that the Carslaw-Jaeger equation correctly describes the temperature distribution in the ground and its variability over time.