Heat Transfer Model to Predict Temperature Distribution 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. 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.


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.

Heat Conduction in the Ground
Heat transfer in the ground occurs mainly as a result of heat conduction.Heat conduction equation is given by: where: T-temperature of the ground ( 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: where: q cond -conductive heat flux (W/m 2 ); H-convective heat flux (W/m 2 ); S-solar radiation heat flux (W/m 2 ); LW-long-wave radiation heat flux (W/m 2 ); EV-evaporative heat flux (W/m 2 ).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: 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: 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: where: ω-frequency of temperature fluctuations.For phenomena occurring in the annual cycle, the frequency is equal to ω = 2π/365 days −1 .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].

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.

Convective Heat Flux
Convective heat flux is determined by the heat transfer equation: ( ) where: h-convective heat transfer coefficient (W/(m 2 K)); Ta-temperature of the air (°C); Tstemperature 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: where: u-wind velocity (m/s).Daily average air temperature changes in the annual cycle as follows: ( ) where: Tam-annual average of the air temperature (°C); Aa-amplitude of air temperature (K); Paphase angle (rad).
An example temporal course of daily average air temperature is shown in Figure 2 (blue line).

Convective Heat Flux
Convective heat flux is determined by the heat transfer equation: where: h-convective heat transfer coefficient (W/(m 2 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: h = 5.7 + 3.8 u for u < 4.88 m/s h = 7.2 u 0.78 for u ≥ 4.88 m/s (7) where: u-wind velocity (m/s).Daily average air temperature changes in the annual cycle as follows: 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(ωt − P s ) 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.The daily average temperature of the surface of the ground is periodically variable, like the air temperature: ( ) where: Tsm-annual average of the ground surface temperature (°C): Aa-amplitude of the ground surface temperature (K); Ps-phase angle (rad).The parameters of Equation ( 9) Tsm, Ps, As are different from the corresponding parameters Tam, Pa and Aa for the air.

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: ( ) where Sm is the annual solar radiation flux absorbed by the ground (W/m 2 ), Asol is the amplitude of the daily average of this radiation flux (W/m 2 ), whereas Psol is the phase angle (rad).

Long-Wave Radiation Heat Flux
The ground radiates heat energy since the temperature of the ground surface Ts is higher than the sky temperature Tsky.The daily average net flux of long-wave radiation LW equals: where σ = 5.67 × 10 −8 W/(m 2 K 4 ) is the Stefan-Boltzmann constant, ε is the emissivity of the surface of the ground and CLW is a constant equal to 4.83 W/(m 2 K).The second part of Equation ( 11) is a linearized form, convenient for modelling.The approximation concerns small absolute differences of Ts and Tsky temperatures.
The sky temperature Tsky 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]:

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: where S m is the annual solar radiation flux absorbed by the ground (W/m 2 ), A sol is the amplitude of the daily average of this radiation flux (W/m 2 ), whereas P sol is the phase angle (rad).

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: where σ = 5.67 × 10 −8 W/(m 2 K 4 ) 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 2 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]: 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: 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).

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 p T + b p (a p = 103 Pa/K, b p = 609 Pa [13]), the evaporative heat flux amounts to: 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 (-).

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: where: In order to determine the P s value, the nonlinear algebraic equation should be solved: Amplitude A s can be determined from the formula: wherein the constants p 1 -p 4 are:

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: 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: 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 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  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  • 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 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.

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.

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 2 , 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 2 (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).
Energies 2018, 11, x FOR PEER REVIEW 10 of 17 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 2 , 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 2 (maximum values).However, there is no discrepancy in the values of conductive heat flux for dry and moist soil, qcond changes slightly (marked as brown symbols).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.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 qcond 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.

• 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 2 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 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 2 , wherein for moist soil is equal to 31.77kWh/m 2 (calculated for h = 13.3W/m 2 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 2 the amount of transferred heat increases by 46% compared to A sol = 101 W/m 2 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 2 K), the emissivity in the 0.5-1.0range does not affect the amount of heat transferred inside the ground.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 2 K), the emissivity in the 0.5-1.0range does not affect the amount of heat transferred inside the ground.
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 2 K), the emissivity in the 0.5-1.0range does not affect the amount of heat transferred inside the ground.

• Direction of Heat Transfer
The directions of the fluxes qcond and H (Figure 10) are variable over time.If Ta > Ts, 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 qcond < 0 and conductive heat flux is directed from the ground surface to its interior (Figure 10b,d).

• 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: 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.8K, P s = 0.166 rad (calculated for a = 0.6•10 −6 m 2 /s).Under these conditions, convective heat flux is negative during 3  4 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.
changes in the direction of thermal fluxes qcond and H.For annual average climatic conditions in Cracow parameters of Equation ( 4) were determined: Tsm = 10.9 °C, As = 13.8K, Ps = 0.166 rad (calculated for a = 0.6•10 −6 m 2 /s).Under these conditions, convective heat flux is negative during ¾ of the year (February 3rd to November 8th) due to the relationship Ta < Ts-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.   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.

• 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 Energies 2019, 12, 25 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.
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 October 15th and December 15th respectively.The calculations were carried out for a = 0.6•10 −6 m 2 /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.

•
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 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 −6 m 2 /s.As can be seen from Figure 12 the change in the value of the emissivity of the ground surface does not affect the qcond 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.
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.

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 qcond 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 qcond 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 (qcond).Amplitude of the daily average solar radiation flux has a large influence on conductive flux, but this value (Asol) 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: Tsm, As, Ps 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.

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.

Energies 2018 ,
11, x FOR PEER REVIEW 4 of 17 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.

Figure 1 .
Figure 1.Heat fluxes on the surface of the ground.

Figure 1 .
Figure 1.Heat fluxes on the surface of the ground.

Figure 2 .
Figure 2. Ambient air and sky temperature in Cracow.

Figure 2 .
Figure 2. Ambient air and sky temperature in Cracow.

Figure 3 .
Figure 3.Comparison of experimental and computational results in Lemont, USA.(a) temperature profiles in the ground; (b) the consistency of the model with the measurement results.•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 Figure4the 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 (Figure4b).The following values of parameters of Equation

Figure 3 .
Figure 3.Comparison of experimental and computational results in Lemont, USA.(a) temperature profiles in the ground; (b) the consistency of the model with the measurement results.• 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.
) have been obtained: T sm = 21.3 • C, A s = 11.9K, P s = 0.232 rad, L = 1.91 m.The sum of squares equals to SS = 38.08 and 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 Figure4a.

Figure 3 .
Figure 3.Comparison of experimental and computational results in Lemont, USA.(a) temperature profiles in the ground; (b) the consistency of the model with the measurement results.

Figure 4 .
Figure 4. Comparison of experimental and computational results in Zarqa, Jordan.(a) temperature profiles in the ground; (b) the consistency of the model with the measurement results.

Figure 5 .
Figure 5. Temporal courses of heat fluxes on the surface of the ground for different values of f.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 2 (during the whole year).Changes in the values of the discussed thermal

Figure 5 .
Figure 5. Temporal courses of heat fluxes on the surface of the ground for different values of f.

Figure 5 .
Figure 5. Temporal courses of heat fluxes on the surface of the ground for different values of f.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 2 (during the whole year).Changes in the values of the discussed thermal fluxes cause minimal changes in the conductive heat flux values.

Figure 6 .
Figure 6.Temporal courses of heat fluxes on the surface of the ground for different values of u.

Energies 2018 , 17 Figure 6 .
Figure 6.Temporal courses of heat fluxes on the surface of the ground for different values of u.

Figure 7 .
Figure 7. Temporal courses of heat fluxes on the surface of the ground for different values of ε.

Figure 7 .
Figure 7. Temporal courses of heat fluxes on the surface of the ground for different values of ε.

17 Figure 8 .
Figure 8. Temporal courses of conductive heat flux for different Sm and f.

Figure 8 .
Figure 8. Temporal courses of conductive heat flux for different S m and f.

Figure 9 .
Figure 9.The impact of h, f and ε on the total amount of heat transferred between the ground and the environment.

Figure 9 .
Figure 9.The impact of h, f and ε on the total amount of heat transferred between the ground and the environment.

Figure 11 .
Figure 11.Temporal courses of fluxes qcond and H.

Figure 11 .
Figure 11.Temporal courses of fluxes q cond and H.

Figure 12 .
Figure 12.Impact of u, f and ε on the temperature profiles in the ground.

Figure 13 .
Figure 13.Temporal changes of the ground temperature at a depth of 0, 1 and 5 m for different values of u.

Figure 13 .
Figure 13.Temporal changes of the ground temperature at a depth of 0, 1 and 5 m for different values of u.

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