A numerical model for ground temperature determination

The ground surface temperature and the temperature with respect to depth are one of the most important issues for geotechnical and environmental applications as well as for plants and other living organisms. In geothermal systems, temperature is directly related to the energy resources in the ground and it influences the efficiency of the ground source system. The ground temperature depends on a very large number of parameters, but it often needs to be evaluated with good accuracy. In the present work, models for the prediction of the ground temperature with a focus on the surface temperature at which all or selected important ground and environmental phenomena are taken into account have been analysed. It has been found that the simplest models and the most complex model may result in a similar temperature variation, yet at a very low depth and for specific cases only. A detailed analysis shows that taking into account different types of pavement or a greater depth requires more complex and advanced models.


Introduction
The ground temperature prediction is a very important issue for a large number of engineering and geothermal applications which use heat pump systems. However, this topic is also important for a number of various environmental applications [1], agricultural applications and thermal energy storage applications [16]. The soil as the heat source and energy storage/sink is widely used for passive heating and cooling applications [18].
The ground surface temperature is also an important issue for the health, epidemic diseases and athlete performance [6]. A high ground surface temperature could raise the human body temperature to a harmful level and generate a heat stroke, heat stress or heat injury.
From the energy's point of view, the ground temperature is the key parameter describing the potential of the thermal energy in the soil and it directly influences the ground base heat pump system efficiency COP (Coefficient of Performance) and SPF (Seasonal Performance Factor) [4]. A proper design of the ground source system is a relatively complex task and requires knowledge from many different fields. The measurement of the ground temperature at different depths is not easy and, for this reason, modelling is a very useful tool for providing the annular soil temperature variation.
In the literature, many different types of models have been proposed in order to model the ground temperature or ground surface temperature [7,8]. Typical modelling approaches consider empirical models [21,22], analytical and semi-analytical solutions [12] or numerical solutions [14,23]. For this purpose, less common approaches are also used: in [13], even neutral networks were implemented in order to estimate the soil surface temperature profiles and [20] used the Fourier method. This allows a ground temperature evaluation, which is helpful in the designing of a ground heat exchanger or in general systems for the ground lower heat source. Because of the availability of many modelling approaches, it is not easy to select the proper model for the energy analysis.
In most energy applications, it is important to know the thermal potential of the ground, but, what is more important, also the realistic local weather conditions. In the system analysis, it is also important to know the impact of the main parameters on the system's behaviour. In order to perform a computer simulation, it is necessary to have information about the local soil's geological structure and the local environmental conditions. Without the above parameters, the local analysis is not able to provide sufficient information for the selection of the system size.
The analysis presented in the literature [12] shows that the near-surface region exhibits a significant temperature fluctuation, dependent on the air temperature [4]. However, in reality, the ground surface temperature and, in consequence, also the ground temperature at different depths, depend on a combination of different weather elements e.g. solar radiation, wind speed, evaporation, precipitation etc. One of the important parameters independent of the weather as well as soil structure which may significantly influence the ground temperature profile is the ground pavement type. The purpose of this research is to investigate an accurate mathematical and numerical model in order to calculate the ground surface and ground temperature profiles, keeping the number of required data to the minimum. The proposed different unsteady models will be tested and the results will be compared for the cases of various ground surface covers and various surface cover structures. The proposed models' formulation includes a description of the heat transfer processes inside the ground and a description of the environmental processes between the soil surface and the surrounding area. The analysis will be performed for the actual environmental and geological conditions. The numerical results will be also compared with the experimental measurements of the ground at different depths.

Mathematical and numerical model
The three-dimensional, unsteady, mathematical model presented in this paper takes into account the transport phenomena inside the ground as well as above the top ground surface. The major model components are listed below and a graphical representation of the energy components including all the important transport phenomena is presented in Figure 1.  The symbols in Figure 1 represent the heat transfer fluxes related to the following phenomena: natural convection ( ), forced convection ( ), solar direct radiation ( ), solar indirect radiation ( , Earth thermal radiation ( ), water evaporation ( ), precipitation ( ), groundwater flow ( ), water phase change ( ) and the Earth natural heat flux ( ).
Assuming unsteadiness of the phenomena, the ground water flow and the phase changes in the soil as well as the variations of the soil properties with the depth, a model equation for the soil in Cartesian coordinates can be written as follows: (1) where are the ground and water densities, respectively, is the specific heat of the ground and water, are the Darcy velocity components of the groundwater flow in x y and z direction. The thermal conductivity of the ground at a selected depth is denoted as . In order to consider water phase change the apparent heat capacity approach was implemented and the effects of the latent heat was represented by an additional source terms s= containing equivalent (apparent) specific heat. It was assumed that the phase change takes place only over a small temperature range (0-0.5 o C). Therefore, the equivalent specific heat was used , where L latent heat of fusion, phase change range. Additionally water property will change according to the formula presented in [11]. The lithological profile was verified in the course of drilling at the University campus and the results are presented in Table 1. The thickness of the soil's lithological layers and the thermo-physical properties are also shown. Table 1. Lithological-stratigraphic profile of the soil with the thermal properties [11].

No. Top
Bottom Thickness Lithology Thermal conductivity Volumetric specific heat 0 0.0 m 0.05 m 0.05 m see Table 2  In general, the flow of the groundwater in the soil is in the order of 10 to 300 m/year [9]. In this study, it is assumed that the underground water flow occurs in the horizontal direction at depths between 4 and 6 m with the Darcy velocity equalling 40 m/year (low range of class 3 flows) [9]. The total heat flux on the top ground surface takes into account the environmental phenomena presented in Figure 1 and it can be written as follows: (2) The first two elements in eq. (2) are related to the direct and indirect solar radiation and the heat fluxes reaching the Earth's surface can be calculated from the following formulas: (3) (4) where is the local Earth's albedo varying in the present paper from 0.1 up to 0.6, depending on the type of surface (see Table 2 for details), is the direct solar radiation and is the diffusive part of the solar radiation in the atmosphere (indirect radiation).
The third component in eq. (2) is related to the Earth thermal radiation. It is assumed in this analysis that the Earth's surface is flat and can be treated as a grey body with the emissivity in the  Table 2 for details). This component is due to the differences between the ambient temperature (also called the sky temperature) and the temperature of the ground surface, and the Earth radiation thermal heat flux can be determined from the following formula: (5) where is the sky temperature. The sky temperature in the present mathematical model is calculated with the use of the following formula [19]: (6)  The next two components in eq.(2) are the forced convection and the natural convection . When the temperature difference between the surface temperature and the air temperature is positive, the natural convection occurs. If the air velocity module is non-zero, the forced convection is also taken into account. Assuming that, at the local scale, the ground can be approximated by a flat plate, the convective heat transport may be calculated as follows: ) (7) where: , are the heat transfer coefficients, in the case of the forced and the natural convection, respectively and, , are the air and the ground surface temperature. In the current model, a semi-empirical solution based on the characteristic numbers (Nusselt, Reynolds, Prandtl and Rayleigh number) was used.
where is the magnitude of the temporal local value of air velocity, is the characteristic dimension defined as the length of the area under consideration, , , and and are the kinematic viscosity, the thermal expansion, the thermal diffusivity of the air and the heat transfer coefficient, respectively. The Reynolds number and the Rayleigh number was in the range Re=0÷6.2·10 6 and Ra=0÷3.8·10 11 respectively. The value of the Nusselt number for the forced convection and the natural convection, is determined based on the following formula for the laminar and turbulent flow over the flat plate [3]: for Re<10 5 and Re>10 5 (9) for 10 4 < Ra < 10 7 and N for 10 7 > Ra The heat flux related to the evaporation and the precipitation depends on several factors (soil moisture, air temperature, velocity and humidity). Because, in the present work, different types of surfaces are considered, in some cases, fully resistant, for the clarity of the analysis, those components will not be taken into account in the present analysis. With the use of a mathematical formulation, three different mathematical models were constructed. The first, Model 1, which takes into account only the air temperature, assumes that the ground top surface has the air temperature. Model 2 assumes only the convective heat transfer and Model 3 takes into account all the considered components. To solve the mathematical model eq.(1)-(7), a numerical algorithm was developed in Fortran 90, based on the finite volume method. For the spatial derivative, three-level methods were implemented and for the convective terms, the central difference scheme was used. The geometry for the analysis contains the area with the dimensions of 20x20x30m, in the directions x, y and z, respectively. The applied mesh (50x50x100) was irregular, fine close to the surface, where the temperature changes were significant. The initial temperature was =9.6 o C (local mean annual temperature). This value was also used as the side walls boundary condition. The temporary values of air temperature , wind speed and solar radiation were measured for a local area: Kraków, Poland(performed in 2014). In order to get periodic solution at presented depths (with the time period of 365 days), the same weather information was repeated 10 times for each case and the time of each simulation was 10 years.

Numerical and experimental results
The numerical model was validated with the experimental measurements [5] of the ground temperature carried out in the region of Bialystok, Poland. The measurement was performed with two differently covered locations (lawn and car park). The first station was situated on a lawn at the distance of about 30 m from a 2-storey building. The second station was situated on a car park covered with bricks at the distance of about 20 m from a low building. For temperature measurements, Cu-Konstantan thermocouples of 0,5 mm diameter ,with the thermocouple meter having the resolution of 0,1 K, were used. The temperature readings were taken every 7 or 14 days at about 1 p.m. For the given soil's thermal properties and the measured air temperature, the ground temperature at two different depths (z=1,z=2m) for the time of one year (2001) with the mean annual local average temperature 8.5 o C (set-up as a initial temperature for verification case only) is presented in Figure 2 and compared with the measurements. Taking into account the fact that the exact environmental conditions were not known, the agreement between the measurements and the numerical prediction is relatively good. The significant source of error is also the unknown initial soil temperature.   In Figure 4, for a lawn surface, and in Figure 5, for an asphalt surface, the temperature of the ground calculated in the 10 th year of the calculation period at the depths z=0.5, 1.0, 2.0, 5.0 m is presented with the use of Model 1-Model 3. It can be seen that in the calculation with the simple models, i.e. Model 1 and Model 2, it is not possible to follow the temperature variation, calculated with the most advanced model, i.e. Model 3. This is independent of the pavement type. The temperature predicted with Model 3 is higher than the temperature obtained with the other models in summer time and it is lower in winter time. The temperature fluctuation amplitude predicted with Model 3 is much higher than that in the case of the other models. This can be well seen in Figure 3, where the calculated cumulative distributions of the temperature at the depth z=5cm for the lawn and asphalt pavement is shown. In general, with Model 3, higher as well as lower temperatures are observed more often. This is observed for the lawn and for the asphalt. In Figure 5, the numerical results for Model 3 and different pavement types at the depth z=0.5 and z=2m are presented together with the temperature histogram of the surface for (b) artificial turf and (d) lawn. It can be seen that the artificial turf temperature is similar to the asphalt temperature and for that pavement type, the observed temperatures during the whole year are the highest. For the lawn surface, the temperature during the year is significantly lower. For the white concrete observed, the temperature is the lowest. These phenomena cannot be correctly reproduced with Model 1 or Model 2.

Conclusion
In this paper, three mathematical models for the ground surface temperature and the ground temperature have been studied. The influence of the environmental model on the calculated results has been tested in order to develop an accurate and well convergent model. It has been found that the simplest models and the most complex model may result in a similar temperature variation, yet at a very low depth and for specific cases only. A detailed analysis shows that taking into account different types of pavement or greater depths requires more complex and advanced models. Simple models are not able to reproduce the correct temperature fluctuation amplitude or the temperature level. This means that, in order to have an accurate correlation between the environmental condition and the surface temperature, the air temperature information is not sufficient and has to be supplemented with other data (wind speed and solar radiation) The use of simple models results in incorrect maximum and minimum values and, even if the mean annual temperature is similar, it will still significantly influence the ground heat exchanger design and the theoretical COP, SPF estimation.