- Email: [email protected]
A new modeling approach for improved ground temperature profile determination
Renewable Energy 85 (2016) 436e444
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
A new modeling approach for improved ground temperature profile determination Messaoud Badache a, *, Parham Eslami-Nejad a, Mohamed Ouzzane a, Zine Aidoun a, Louis Lamarche b a b
CanmetENERGY Natural Resources Canada, 1615 Lionel Boulet Blvd., P.O.Box 4800, Varennes, Qu ebec J3X1S6, Canada Department of Mechanical Engineering, Ecole de technologie sup erieure, 1100 Rue Notre-Dame Ouest, Montr eal, Qu ebec H3C1K3, Canada
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 December 2014 Received in revised form 5 June 2015 Accepted 9 June 2015 Available online xxx
The knowledge of the ground temperature profile with respect to time and depth is very important in many technological fields like geothermal heat pumps, solar energy systems and geotechnical applications. Many researches were performed in the past in order to evaluate this profile. The most common ones, known as energy balance models, use the energy balance equation as a boundary condition. Unfortunately the performance of these models strongly depends on an accurate estimation of several input factors. The objective of this paper is to develop an improved model for the prediction of the ground temperature profile in which the energy balance equation at the ground surface is supplemented by an empirical correlation for the annual average ground surface temperature calculation. This model is less sensitive to uncertainties of input factors. Furthermore, unlike the previous models, a periodic variation of the sky temperature is introduced instead of a previously assumed constant value. The model is validated against measured data in a site located in Varennes (Montreal-Canada) and two further sites, Fort Collins (Colorado) and Temple (Texas) in the U.S.A. © 2015 Crown Copyright and ELSEVIER Ltd. Published by Elsevier Ltd. All rights reserved.
Keywords: Ground temperature profile Energy balance model Undisturbed ground temperature
1. Introduction The assessment of thermal performance in many engineering applications requires the estimation of the ground temperature as a function of time, location, and depth. Examples of these applications are buried ground heat exchangers connected to ground source heat pump systems or energy storage systems [1,2]. Earthto-air heat exchangers for heating and cooling of buildings and agricultural greenhouses [3e5] are two other applications. Ground temperature is also important in the design of airport, road pavements, pipelines, as well as for buried high-voltage power cables and nuclear waste disposal facilities. Many factors such as slope orientation, terrain, solar radiation, wind, rain, etc., can affect the ground thermal behavior to a more or less important extent [6] and a reliable assessment approach is highly needed. Since few measured data are available on ground temperature,
* Corresponding author. E-mail addresses: [email protected] (M. Badache), Parham. [email protected] (P. Eslami-Nejad), [email protected] (M. Ouzzane), [email protected] (Z. Aidoun), [email protected] (L. Lamarche).
several models have been developed for its estimation as a function of time and depth. Determination of this temperature is generally based on the solution of the transient one dimensional heat conduction problem in the ground (Equation (1)), and depends on the type of boundary condition at the ground surface. A full discussion of this subject can be found in Ref. [7]. The boundary condition can take two forms depending on the available data [8]. If the ground surface temperature is known, its temporal variation is the boundary condition. If this information is not available, the energy balance equation can be used as boundary condition at the ground surface. The present work has been performed in this latter context. It is worth noting that, both kinds of boundary conditions lead to solutions of the same general form (Equation (2)) but require different inputs in order to determine their respective parameters. For a given surface temperature, the most commonly used model is that presented by Refs. [9,10]. More details about this equation including all parameters and variables are given in Ref. [9]. This model seems to be very simple to use. However it requires the knowledge of the ground temperature at or near the surface which is relatively scarce and without it this model cannot be used. For the surface energy balance models, meteorological data are used in the boundary condition accounting for all heat gains and
http://dx.doi.org/10.1016/j.renene.2015.06.020 0960-1481/© 2015 Crown Copyright and ELSEVIER Ltd. Published by Elsevier Ltd. All rights reserved.
M. Badache et al. / Renewable Energy 85 (2016) 436e444
Nomenclature
d u f
annual damping depth, (m) emissivity angular frequency, (radian/day) phase angle, (radian)
F
heat flux, (W/m2)
ε
cp G h k rh t T V Z
specific heat at constant pressure, (J/kg,K) Intensity of solar radiation, (W/m2) surface convective heat transfer coefficient, (W/m2,K) thermal conductivity, (W/m,K) relative humidity of the air time, (day) temperature, (K) wind velocity, (m/s) depth of the ground, (m)
Greek symbols thermal diffusivity, (m2/s) absorption coefficient StefaneBoltzmann constant, (W/m2,K4)
a as s
losses over the ground surface [11]. were among the earliest who adopted an energy balance model to investigate the daily and annual variation of the ground temperature in Kuwait. They developed an expression of the ground temperature based on the periodic variation of solar radiation and air temperature, making use of the sol-air temperature definition. Latent heat exchange due to evaporation has not been taken into account in their energy balance equation [12]. reproduced the analysis of [11] by introducing the sol-air-evaporation temperature definition. A parametric study was carried out to investigate the effect of various parameters (air relative humidity, ground absorptivity, wind speed, and evaporation fraction) on the heat transfer in the ground [13,14] presented two other models, very similar to that of [11]. They added another parameter to the energy balance equation which is the latent heat exchange due to evaporation. A drawback of the aforementioned energy balance models is that their performances depend strongly on the accuracy of the weather data and the estimation of the ground input factors. Many of these input factors are determined by empirical correlations (e.g., convective heat transfer coefficients), some are experimental constants and others are effective mean values (e.g., thermal conductivity and diffusivity of the ground). In addition, while the long wave radiation (Fsky ) depends on soil radiative properties, air relative humidity and effective sky temperature it was assumed in the aforementioned models [12e14] to be constant and equal to 63 W/m2 for all locations in the world. While this value fits conditions on which it is selected and similar conditions elsewhere, it becomes a rough approximation when different conditions have to be considered. The objective of this paper is to develop an improved model that predicts the ground temperature as a function of depth and time, based on the meteorological data. This model is less affected by the inaccuracy of input factors estimation. The energy balance equation is used as a boundary condition at the ground surface in order to determine the amplitude and the phase angle of the temperature at this location while the annual average ground surface temperature is calculated using an empirical correlation. Unlike previous energy balance models, besides assuming a periodic variation of solar radiation and ambient air temperature, this model assumes a periodic variation of the sky temperature.
2. Model development The equation describing the transient, one dimensional heat conduction in the ground is:
437
Subscripts Air air amb ambient dp dew point s ground surface sky sky conv convection evap evaporation
dT dT 2 ¼a 2 dt dz
(1)
The ground is assumed to be a semi-infinite medium with constant physical properties. Considering the following initial and boundary conditions: Tðz; 0Þ ¼ T s at t ¼ 0 Ts ¼ Tð0; tÞ ¼ T s As cosð ut fs Þ at z ¼ 0 Tð∞; tÞ ¼ T s at z ¼ ∞ The solution of this problem is given by Equation (2) [8,9].
z z Tðz; tÞ ¼ T s As :exp d cos ut fs d
(2)
This form is used in several software: TRNSYS (2005) [15], DOE2 (1982) [16] and RETScreen (2005) [17].where: T(z,t) is the ground temperature at time t (day or hour) and depth z (m). T s is the annual average ground surface temperature, equivalent of the undisturbed ground temperature. As is the annual amplitude of the ground surface temperature. d is the damping depth (m) of annual fluctuation of the ground temperature which can be calculated from:
d¼
rffiffiffiffiffiffi 2a u
(3)
where a is the thermal diffusivity and u ¼ 2 p/365 (radian/days) is the angular frequency.
fs is the phase angle (radian). At the ground surface (z ¼ 0) the energy balance is applied as a boundary condition (Equation (4)) in order to calculate the parameters T s ,As,4s. The details of the energy balance are given in Fig. 1. This figure represents the main heat flux contributions at the surface of the ground including the conduction heat flux into the ground, the convective heat flux transferred between the surface and the ambient air, the short-wave global solar radiation (radiant flux) absorbed by the ground surface, the long wave radiant flux exchanged with the surroundings (sky), and the latent heat flux of evaporation.
438
M. Badache et al. / Renewable Energy 85 (2016) 436e444
humidity. It can be evaluated from the following equation [13,18]:
Fevap ¼ c:f :hconv ½ðaTs þ bÞ rh ðaTamb þ bÞ
(8)
where a ¼ 103 (Pa/K), b ¼ 609 (Pa) and c ¼ 0.0168 (K/Pa) rh is the relative humidity of the air, while f is a fraction of evaporation rate, varying between 0 and 1, and depending mainly on the ground cover and moisture. For bare soil, the fraction f can be estimated as follow: For saturated soil, f ¼ 1, for moist soil, f ¼ 0.6e0.8. For dry soil, f ¼ 0.4e0.5, for arid soil, f ¼ 0.1e0.2. By substituting Equations (5)e(8) into Equation (4), the boundary condition is written as follow:
hr Tamb þ as G þ hrad Tsky c:hconv f :b:ð1 rh Þ ¼ he Ts k
vT vz z¼0 (9)
where hr ¼ hconv(1 þ c . a.f.rh) and he ¼ hconv(1 þ c .a .f)þ hrad Fig. 1. Main heat flux contributions at the surface of the ground.
2.1. Fourier series approximation
The sign convention used in the energy balance assumed that heat transfer to the ground surface is positive, while heat transfer away is negative.
Fnet ¼ k
vT ¼ Fconv þ as G Fsky Fevap vz z¼0
(4) Tðz; tÞ ¼ T s þ Real
where Fnet represents the net heat flux transferred by conduction across the ground surface. The equations described below provide the necessary relationships representing each of the terms in Equation (4). Fconv is the sensible heat flux exchanged between the air and the ground surface, expressed as:
Fconv ¼ hconv ðTamb Ts Þ
(5)
where hconv is the convective heat transfer coefficient. asG is the short-wave global solar radiation (radiant flux) absorbed by the ground surface. Fsky is the long-wave radiation (radiant flux) emitted by the ground surface to the sky and calculated by Equation (6). Since air temperature and solar radiation, show a cyclic behavior throughout the day and the year, it is reasonable to assume that the sky temperature also follows such a cyclic pattern. Therefore, a sinusoidal function can be used to correlate the sky temperature with the meteorological data. This representation of the sky temperature gives a more accurate evaluation of Fsky in comparison to that given by a constant value.
Fsky ¼ hrad Ts Tsky
with
For the sake of simplicity, the solution of Equation (1) is expressed using the complex form of Fourier series (Equation (10)). This expression can be applied to the surface as well as to the various depths. All assumptions made above have been considered herein.
hrad ¼ εs s Ts þ Tsky
N X
1þi where g ¼ 0 and d0 ¼ d
2 Ts2 þ Tsky
The term hrad represents the constant part of the linearized term, of short wave solar radiation, εs is the emissivity of the ground surface and s is the StefaneBoltzmann constant (5.670 108 W/ m2 K4). Fevap is the evaporation heat exchange flux that depends on several factors including ground surface and ambient air temperatures, wind speed, surface cover, soil moisture content and air
rffiffiffiffiffiffi 2a nu
(11)
equal to, Arctan (an/bn), and an and bn are Fourier coefficients for the wave temperature given as:
an ¼
365 2 X TðtÞsinðnutÞ 365 j¼1
(12)
bn ¼
365 2 X TðtÞcosðnutÞ 365 j¼1
(13)
d0 is the nth harmonic damping depth and ‘’Real’’ denotes the real part of the complex number. The ground surface temperature Ts (0,t) in the complex form can be written as: Ts ð0; tÞ ¼ T s þ Real
(7)
(10)
where n is the number of the harmonic, t is the time (number of the day), As is the amplitude of the nth harmonic and is equal qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi to a2n þ b2n , fn (radians) is the nth harmonic phase angle and is
(6)
As einutzgfn
n¼1
N X
As eifs einut
(14)
n¼1
The complex form of Fourier series is also used to express the periodic behavior of the solar radiative flux, sky and air temperature as follows:
Tamb ðtÞ ¼ T amb þ Real
N X n¼1
Aa eifa einut
(15)
M. Badache et al. / Renewable Energy 85 (2016) 436e444
GðtÞ ¼ G þ Real
N X
AG e
ifG inut
e
(16)
n¼1
Tsky ðtÞ ¼ T sky þ Real
N X
Asky eifsky einut
(17)
n¼1
By substituting Equations 10e17 into Equation (9) and equating the time dependent parts, the amplitude of the ground surface temperature, As, and the phase angle difference between the air and the ground surface temperature, fs- fa in Equation (2) are given by:
h A þ a A eiðfG fa Þ þ h A eiðfsky fa Þ s G r a rad sky As ¼ ðhe þ k: g Þ "
hr Aa þ as AG eiðfG fa Þ þ hrad Asky eiðfsky fa Þ fs fa ¼ Arg ðhe þ k:g Þ
(18) #
(19) In Equations (18) And (19) the symbols kk and Arg are respectively used to indicate the modulus and the argument of a complex number. By equating the time independent parts, the annual mean ground surface temperature can be calculated by the following equation:
Ts ¼
i 1h hr T amb þ hrad T sky þ as G c:hconv :f :bð1 rh Þ he
(20)
The Equations (18)e(20) provide the necessary parameters of Equation (2) to calculate the ground temperature as a function of depth and time with respect to the energy balance equation applied at the ground surface. These equations are different from those of [13,14] because they include more realistic effects through additional terms which are fsky, Asky, hrad, and Tsky. 2.2. Model sensitivity to factors estimation Many input factors such as hconv, as, εs, f and rh are involved to determine the parameters T s ,As, 4s in Equations (18)e(20). While rh and f are relatively simple to estimate, there is a great number of available correlations in the literature for hconv, as, εs. For the absorption factor, as, different correlations relating this factor to the surface water content are used [19,20]. The emissivity of the soil surface, 3s, generally depends on the water content and vegetation of the surface, the commonly used correlations are [21e25] reviewed a large number of convective heat transfer coefficient correlations in linear, power law and boundary layer form. Generally hconv correlations depend on the wind velocity and several specific conditions under which they have been correlated. As there is a large number of correlations available it is difficult to assess the impact of all correlations factors on the present model predictive performance. As a result, the following analysis stresses the impact of convective heat transfer coefficient values using three representative correlations on the ground temperature profile. The ground characteristics of a site located in Varennes are taken as the base case example for calculations. Equation (2) is used to predict the ground temperature for two different days (21st march and 21st September) at different depths. Correlations of Givoni and Mostrel [26], Kasuma [27] and Mc Adams [28] have been selected for the convective heat transfer coefficient calculations. Fig. 2 represents the ground temperature profile, calculated using the aforementioned correlations. The difference between the calculated values
439
using Kasuma and Mc Adams correlations are relatively insignificant. However, there is a significant difference between the values calculated using Givoni and the other two correlations (Kasuma and Mc Adams). It shows also that Givoni gives the highest ground temperature prediction for the shallow (1 me8 m), and the deep zone (below 8 m). This result can be explained by the fact that the convective heat transfer coefficient may be underestimated by Givoni and Mostrel correlation (3 W/m2.K) or overestimated by Kasuma correlation (47 W/m2 K). An underestimation of the convection heat transfer coefficient overestimates the ground temperature. The opposite is true when overestimating the convective heat transfer coefficient (Fig. 2). In summary, the above analysis clearly shows that, an inaccurate estimation of the convective heat transfer coefficient, due to inappropriate choice of the correlation factors results in inaccurate model parameters (T s , As, fs) and ground temperature prediction. Therefore it is important for designers and modelers to select reliable convective heat transfer coefficients correlations for T s, As and fs computations. The same remark holds true for as and εs factors. 2.3. Model enhancement features As stated in Equation (2), the present model is a function of three main parametersT s , As and fs which are strongly dependent on several input variables. The importance of accurate estimation of input variable has been stressed upon in the previous section. It is the purpose of this section to improve the model in order to develop an accurate and easy to use correlation for ground temperature profile. Equation (2) implies that the undisturbed ground temperature is constant, while the amplitude and phase angle are damped with a factor of e(z/d). Therefore the potential inaccuracy in As and fs decreases very rapidly in the first few meters from the down from the ground surface, while inaccuracy in T s remains unchanged with depth. This means that, calculating accurate values of T s is crucial for a correct model temperature prediction. It is obvious that accurate values of input factors result in accurate calculation of T s . Unfortunately, as shown above it is not always straightforward because of variety of available correlations. A solution to this problem is to use a form for T s that is independent of the aforementioned input factors. Recently a new correlation (Equation (21)) has been developed by Ref. [29] for calculation of T s . It is only a function of the mean annual air ambient temperature
Fig. 2. Ground temperature prediction using three different correlations for the calculation of the convective heat transfer coefficient.
440
M. Badache et al. / Renewable Energy 85 (2016) 436e444
(available factor in all meteorological databases). This correlation can be used in Equation (2) instead of Equation (20). Therefore, the improved ground temperature model proposed in this paper is based on Equations (18), (19) and (21).
T s ¼ 17:898 þ 0:951T amb
(21)
where T amb is the annual average ambient air temperature in degree Kelvin. 3. Results and discussion In this section, first the developed model is applied to predict the daily ground temperature of the Varennes' (Montreal region) site for various depths and time. The model predictions are compared with experimental results and the profile of Varennes ground temperature is then analysed. Further, the model is validated for two other sites in the U.S.A: Fort Collins (Colorado) and Temple (Texas). The following is a description of the procedure used for ground temperature prediction of Varennes using this model. The data of a typical meteorological year (TMY) [30] of Montreal (global solar radiation (G), ambient dry bulb (Tamb), dew point temperatures (Tdp), and wind velocity (Vw)) are used for the determination of the model inputs, including the mean values (T amb , G, T sky ), the phase angles (4a,4G,4sky), and the amplitudes (Aa AG Asky ). Fourier analysis (FA) is used to extract these parameters. The calculations are limited to the base case of the first harmonic (n ¼ 1). The first harmonic accounts for 95% of the variance in the worst cases. Table 1 gives the calculated mean values, the phase angles, and amplitudes of the parameters G, Tamb and Tsky for the site of Varennes. The sky temperature is calculated by Equation (22) [31].
h i0:25 2 Tsky ¼ Tamb 0:711þ0:0056Tdp þ0:000073Tdp þ0:013 cosð15tÞ (22)
was fractured black shale from 0 m to 19.5 m, slate from 19.5 m to 23.2 m and black shale from 23.2 m to 31 m. Twenty four 24 calibrated thermocouples (T type) were attached on the surface of both legs of the U tube at different depths. The thermocouples were used to measure the ground temperature profile from February to November 2014. Because of the important variation of the temperature near the ground surface, more thermocouples were located in this region (Fig. 3) while beyond 11 m deep they were fewer and equally spaced. Temperature data were recorded every 5 min, using data acquisition system (DAQ) and then compiled into the daily-average format. Daily-average formats are calculated by integrating 5 min data over the period (24 h) using an in-house Fortran code. The acquisition system was fitted with two extensions of twenty entries each and communications ports including Ethernet and USB. The system used for the temperature measurement is composed of five independent sources of error: the DAQ, the temperature sensor itself, the cold junction compensation (CJC) sensor, thermocouples position, and cabling. The uncertainty on the temperature measurements is ±0.42%. It is noted that the error due to thermocouples position and cabling is ignored so much it is so small (3 105). The following observations were made: The thermocouples at the same level in both borehole legs give the same temperature measurements. Therefore only one leg measurements are considered. The daily mean temperature curves (at depths less than 3 m) against the day number vary approximately in a non-sinusoidal form. This can be due to many parameters such as short term weather variations, seasonal variations, ground water content, and non-homogeneous ground thermal properties. The amplitude of the ground temperature oscillation reduces from around 15 C at 1 m to 0.17 C at 31 m. The shift in phase angle gradually increases with depth. Temperature measurements show that at depths greater than about 8 m below the surface, the ground temperature remains relatively constant (9.5 C).
where Tsky and Tamb are in degree Kelvin and Tdp is in degrees Celsius. Parameter t is the hour from midnight. Calculation of As and 4s parameters (from Equations (18) and (19)) involves the determination of f, εs, as, d0 , hconv, and hrad. Parameter d0 is a function of the thermal properties of the ground and the frequency of the thermal oscillation. Effective ground thermal properties obtained from a thermal response test (TRT) report are k ¼ 2.65 W/m$K and a ¼ 0.0948 m2/day. The TRT test was carried out by a contractor to determine these properties based on ASHRAE standards with less than 5% inaccuracy. Recent guidelines for the test procedure can be found in Ref. [32]. It is assumed that the coefficients hconv, f, rh, εs, as do not vary during the year and their annual average values are used. The coefficients f, εs, and as are assumed to be 0.7, 0.9, and 0.9, respectively. Therefore the annual value of d0 calculated from Equation (3) is 3.31 m. In the case of Montreal (cold climate) [33], suggest the use of the McAdams correlation, in order to estimate hconv. The mean annual wind velocity is given to be 4.5 m/s. Since hrad is a function of ground surface temperature (Equation (7)), an iterative method is used to calculate the daily values of this factor. 3.1. Experimental set up in Varennes site A U-tube borehole (31 m deep) was drilled at CanmeENERGYVarennes. This set-up was fitted with spacers which were placed in the center of the U-tube to keep the tube legs apart. Super Grout was used for the filling material. The soil stratigraphy encountered
Fig. 3. Position of the thermocouples in the borehole.
M. Badache et al. / Renewable Energy 85 (2016) 436e444
441
Table 1 Mean annual value, phase angle, and amplitude of the meteorological parameters G, Tamb and Tsky for site of Varennes. Parameter
Global irradiance (G)
Dry bulb temperature (Tamb)
Sky temperature (Tsky)
Amplitude Phase angle Mean annual value
160.57 0.28 103.12
15.59 0.35 6.17
21.33 0.37 15.25
3.2. Comparison with experimental measurements Figs. 4 and 5 show the comparison between the measured and predicted ground temperatures at eight different depths during the period from February 4th to November 25th, 2014. Generally, the agreement between the predicted and the measured values is good. However, the predicted ground temperature values appear to be overestimated in the spring time (day 50e150) as clearly shown for depths 1 and 2 m below the soil surface. This temperature difference vanishes gradually with depth. It can be seen also that the measured ground temperature variation does not quite take a sinusoidal form. This may be due to the nonhomogeneous thermal properties of the ground. Thermal properties change with time and depth due to the volume increase of the ground water content during the thaw period (spring time). In the present study however, properties are assumed to be constant over time and location, as measured in November 2010. During the thaw period the thermal diffusivity and therefore the damping depth increases due to the increase of ground water content. At any given depth, the relative value of the ground temperature amplitude is proportional to the damping depth as shown in equation (23). Consequently, an increase of damping depth brings the surface and deep ground temperatures closer together, thus resulting in a slight drop in ground temperature for this period of the year (day 50e150). ðzz0 Þ= d
Az=A ¼ exp s
(23)
To provide a more quantitative basis for this comparison we used a statistical analysis. The Root Mean Square Error (RMSE) and the Normalized Root Mean Square Error (NRMSE) statistics are calculated. RMSE was generally 2.3 K in the worst cases, i.e. at depth of 1 m and less than 1.5 K for depths below 2 m from the soil surface. NRMSE is less than 0.92%. Despite the fact that the data were measured over approximately ten months and the RMSE is expected to be lower for an entire year, the results are very satisfactory. In other words, the present model predicts quite accurately the ground temperature profile over time and depth.
3.3. Analysis of the ground temperature profile in Varennes Figs. 4 and 5 above clearly show the damping of the annual ground temperature in which the amplitude reduces from around 15 C at 1 m to 0.17 C at 26 m. Equation (2) shows that the amplitude decreases exponentially with distance from the surface at a characteristic length known as the damping depth, d. The damping depth is the depth at which the amplitude of the temperature is 1/e times the amplitude of the temperature at the surface. The mean value of the annual damping depth of Varennes is estimated to be 2.35 m from the measured ground temperature. According to this value the amplitude of the ground temperature is almost completely (95%) damped at 3d (i.e. at 7 m). This is confirmed in Fig. 5 (b), (c) and (d), where the ground temperatures are essentially constant throughout the year for depths below 6 m.
Fig. 4. Measured and predicted ground temperature at depths 1 m, 2 m, 3 m, and 4 m at Varennes for the period from February 4th to November 25th, 2014.
442
M. Badache et al. / Renewable Energy 85 (2016) 436e444
Fig. 5. Measured and predicted ground temperature at depths 6 m, 8.5 m, 16 m, and 26 m at Varennes for the period from February 4th to November 25th, 2014.
These figures show also that below the depth of 1 m, the maximum or minimum temperature occurs later than it does at 1 m, meaning that the time lag increases with the depth. For example the ground temperature is at its minimum on days 38e40 (7th to 9th of February) for the depth 1 m while it is the lowest 2.5 months later (end of April) for the depth of 4 m. Finally, it is interesting to analyse the model predictions on a daily basis. Therefore the profiles of the ground temperature for two days (April 10th and September 2nd) are presented in Fig. 6. This figure shows that due to the temperature fluctuations on the ground surface, the temperature profiles are different. For each day, two different regions are observed; near the ground surface (up to 8.5 m) and deep down to 31 m. The ground temperature profile
Fig. 6. Measured and calculated ground temperature at different depthat Varennes for the April 10th and September 2nd, 2014.
reaches the same value for both days, which corresponds to the undisturbed ground temperature. This confirms our previous observation in Fig. 5c and d. The undisturbed ground temperature value, predicted to be 10.5 C, is 1 C higher than the experimental value, which is about 9.5 C.
3.4. Comparison with additional experimental data In order to extend the application of the new model to different climate conditions, it was applied for the sake of example to predict monthly average ground temperature of two sites in the U.S.: Fort Collins (Colorado) and Temple (Texas). Temple has a warm humid temperate climate with hot summers while Fort Collins has a humid continental climate. The ground temperature data for these two locations are given in the reference article of Kusuda et al. (1965). The mean, the amplitude and phase angle for solar radiation, sky and ambient air temperatures are determined using FA method from TMY data (Table 2). The FA coefficients are calculated for the first harmonic (n ¼ 1). The thermal properties of the ground at Fort Collins and Temple are not available but they can be estimated from the measured temperature data, Equations (3) and (23). Equation (23) is used to estimate the damping depth between 0.0254 m and 1.83 m Equation (3) is used to calculate the corresponding thermal diffusivity. Table 3 shows the parameters used to predict the ground temperature of the two sites. The data (hconv, rh, f, k, a, as, εs) of Fort Collins are obtained from Ref. [13] while those of Temple are from Ref. [30]. Equations (18) (19) and (21) are used to calculate the monthly ground temperature for Fort Collins at 0.1 m and 1.83 m, and for Temple at 0.1 m and 1.23 m from the surface. Figs. 7 and 8 provide the measured and predicted monthly ground temperature for Fort Collins at 0.1 m and 1.83 m, and for Temple at 0.1 m and 1.23 m from the surface. Despite all potential uncertainties that may be accounted for assuming some factors
M. Badache et al. / Renewable Energy 85 (2016) 436e444
443
Table 2 Calculated mean annual value, phase angle, and amplitude of the meteorological parameters G, Tamb and Tsky for Fort Collins, Colorado and Temple, Texas. Factor
T amb (ºC)
fa (rad)
Aa (Rad)
G (W/m2)
fG (rad)
AG (W/m2)
T sky (ºC)
fsky (rad)
Asky (ºC)
Fort collins (Colorado) Temple (Texas)
6.14 18.8
0.23 0.23
14.5 11.15
157 86.5
0.13 0.04
132 202.1
14.9 2.65
0.23 0.23
16.1 16.4
Table 3 Parameters used to predict Fort collins and temple ground temperature. Factor
hconv (W/m2.K)
rh
f
ɑs
ks (W/m.K)
ɑ (m2/day)
d (m)
3s
Fort collins (Colorado) Temple (Texas)
22.0 23.6
0.55 0.70
0.6 0.3
0.8 0.8
1 1
0.0362 0.0279
2.48 2.10
0.9 0.9
used as a boundary condition at the ground surface in order to determine the amplitude and the phase angle of the temperature at this location. The mean ground surface temperature is then calculated using an empirical correlation. Unlike the previous energy balance models, this model is less sensitive to input factors inaccuracy. Furthermore, besides assuming a periodic variation of the solar radiation and ambient air temperature, it also assumes a periodic variation of the sky temperature. Validation against measured data of a site located in Varennes (Montreal-Canada) and two sites in the U.S.A; Fort Collins (Colorado) and Temple (Texas) was performed. Generally the agreement between measurements and the model predictions is good. Therefore, the developed model is proposed as useful tool for the ground temperature profile estimation in the context of various environmental and energy applications. Acknowledgments Fig. 7. Predicted and measured values of ground temperature of Fort Collins (Colorado) at 0.1 m and 1.83 m from the surface.
values, the model still predicts adequately the ground temperature variation for both sites. 4. Conclusion A new model was developed to improve the prediction of the temperature variation in the ground as a function of depth and time based on the meteorological data. The energy balance equation is
Fig. 8. Predicted and measured values of ground temperature of Temple (Texas) at 0.1 m and 1.22 m from the surface.
Financial support for this work was provided by the Canadian Federal Government's Program on EcoENERGY Innovation Initiative (Natural Resources Canada) (EEBI 022). References [1] G. Florides, S. Kalogirou, Ground heat exchangers-A review of systems, models and applications, Renew. Energy 32 (2007) 2461e2478. [2] H. Esen, M. Inalli, M. Esen, K. Pihtili, Energy and exergy analysis of a groundcoupled heat pump system with two horizontal ground heat exchangers, Build. Environ. 42 (10) (2007) 3606e3615. [3] G. Mihalakakou, J.O. Lewis, M. Santamouris, The influence of different ground covers on the heating potential of earth-to-air heat exchangers, Renew. Energy 7 (1996) 33e46. [4] M. Esen, T. Yuksel, Experimental evaluation of using various renewable energy sources for heating a greenhouse, Energy Build. 65 (2013) 340e351. [5] H. Esen, M. Inalli, M. Esen, K. Pihtili, Numerical and experimental analysis of a horizontal ground-coupled heat pump system, Build. Environ. 42 (3) (2007) 1126e1134. [6] G.P. Williams, L.W. Gold, Ground Temperatures, National Research Council Canada, Institute for Research in Construction, 1976. Available from: http://irc. nrc-cnrc.gc.ca/pubs/cbd/cbd180e.html. [7] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Physics Today, 1962. [8] B. Givoni, L. Katz, Earth temperatures and underground buildings, Energ. Build. 8 (1985) 15e25, http://dx.doi.org/10.1016/0378-7788(85)90011-8. [9] T. Kasuda, P.R. Archenbach, Earth temperature and thermal diffusivity at selected stations in the United States, ASHRAE Trans. 71 (1965). Part1. Available from: http://oai.dtic.mil/oai/oai? verb¼getRecord&metadataPrefix¼html&identifier¼AD0472916. [10] W.R. Van Wijk, D.A. De Vries, Periodic Temperature Variations in an Homogeneous Soil. Physics of Plant Environment, North-Holland Publ, Amsterdam, 1966, pp. p.102ep.143. [11] A.K. Khatry, M.S. Sodha, M.A. Malik, Periodic variation of ground temperature with depth, Sol. Energy 20 (1978) 425e427. [12] M.M. Salah El-Din, On the heat flow into the ground, Renew. Energy 18 (1999) 473e490. [13] M. Krarti, C. Lopez-Alonzo, D.E. Claridge, J.F. Kreider, Analytical model to predict annual soil surface temperature variation, J. Sol Energy Eng. 117 (2) (1995) 91e99.
444
M. Badache et al. / Renewable Energy 85 (2016) 436e444
[14] G. Mihalakakou, M. Santamouris, J.O. Lewis, D.N. Asimakopoulos, On the application of the energy balance equation to predict ground temperature profiles, Sol. Energy 60 (3) (1997) 181e190. [15] TRNSYS. Version 16.a. Transient System Simulation Program Solar Energy Laboratory, University of Wisconsin, Madison, USA, 2005. [16] DOE-2, Technical Information Center. Engineers Manual Version 2.1A, Department of Energy, USA, 1982. [17] RETscreen 4, Clean Energy Project Analysis: RESTscreen Engineering & Cases Textbook, Ground-Source Heat Pump Project Analysis Chapter, CanmetENERGY at Varennes, Natural Resources Canada, 2005. [18] G. Mihalakakou, On estimating soil surface temperature profiles, Energy Build. 34 (3) (2005) 251e259. [19] C.H. Van Bavel, D.I. Hillel, Calculating potential and actual evaporation from a bare soil surface by simulation of concurrent flow of water and heat, Agric. Meteorol. 17 (1976) 453e476. [20] J. Kondo, N. Saigusa, T. Sato, A model and experimental study of evaporation from bare soil surfaces, J. Appl. Meteorol. 31 (3) (1992) 304e312. [21] S. Ortega-Farias, R. Antonioletti, A. Olioso, Net radiation model evaluation at an hourly time step for Mediterranean conditions, Agronomie 20 (2000) 157e164. [22] M. Fuchs, C.B. Tanner, Infrared thermometry of vegetation, Agron. J. 58 (1976) 597e601. [23] K. Noborio, K.J. McInnes, J.L. Heilman, Two-dimensional model for water, heat,
[24] [25] [26] [27] [28] [29]
[30] [31] [32] [33]
and solute transport in furrow irrigated soil: II. Field evaluation, Soil Sci. Soc. Am. 60 (1996) 1010e1021. W. Brutsaert, On a derivable formula for long-wave radiation from clear skies, Water Resour. Res. 11 (1976) 742e744. J.A. Palyvos, A survey of wind convection coefficient correlations for building envelope energy system's modeling, Appl. Therm. Eng. 28 (2008) 801e808. B. Givoni, M. Mostrel, Windscreens in radiant cooling, Passive Sol. J. 1 (1982) 229e238. G. Kusuma, Estimation of the exchange coefficient of heat during low wind convective conditions, Bound. Layer. Meteorol. 111 (2004) 247e273. W.H. McAdams, Heat Transmission, McGraw-Hill, New York, 1954. M. Ouzzane, P. Eslami-Nejad, M. Badache, Z. Aidoun, New correlations for the prediction of the undisturbed ground temperature, Geothermics 53 (2015) 379e384. National Solar Radiation Data Base, web address: http://rredc.nrel.gov/solar/ old_data/nsrdb/1991-2005/tmy3/by_USAFN.html. J.A. Duffie, W.A. Beckman, Solar Engineering of Thermal Processes, third ed., John Wiley & Sons, Inc, USA, 2006. ISBN-13 978-0-471-69867-8. ASHRAE, ASHRAE Hnadbook: HVAC Applications. Chapter 32, ASHRAE Inc., Atlanta, GA, 2007. M. Ouzzane, P. Eslami-Nejad, Z. Aidoun, L. Lamarche, Analysis of the convective heat exchange effect on the undisturbed ground temperature, Sol. Energy 108 (2014) 340e347.
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
A new modeling approach for improved ground temperature profile determination Messaoud Badache a, *, Parham Eslami-Nejad a, Mohamed Ouzzane a, Zine Aidoun a, Louis Lamarche b a b
CanmetENERGY Natural Resources Canada, 1615 Lionel Boulet Blvd., P.O.Box 4800, Varennes, Qu ebec J3X1S6, Canada Department of Mechanical Engineering, Ecole de technologie sup erieure, 1100 Rue Notre-Dame Ouest, Montr eal, Qu ebec H3C1K3, Canada
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 December 2014 Received in revised form 5 June 2015 Accepted 9 June 2015 Available online xxx
The knowledge of the ground temperature profile with respect to time and depth is very important in many technological fields like geothermal heat pumps, solar energy systems and geotechnical applications. Many researches were performed in the past in order to evaluate this profile. The most common ones, known as energy balance models, use the energy balance equation as a boundary condition. Unfortunately the performance of these models strongly depends on an accurate estimation of several input factors. The objective of this paper is to develop an improved model for the prediction of the ground temperature profile in which the energy balance equation at the ground surface is supplemented by an empirical correlation for the annual average ground surface temperature calculation. This model is less sensitive to uncertainties of input factors. Furthermore, unlike the previous models, a periodic variation of the sky temperature is introduced instead of a previously assumed constant value. The model is validated against measured data in a site located in Varennes (Montreal-Canada) and two further sites, Fort Collins (Colorado) and Temple (Texas) in the U.S.A. © 2015 Crown Copyright and ELSEVIER Ltd. Published by Elsevier Ltd. All rights reserved.
Keywords: Ground temperature profile Energy balance model Undisturbed ground temperature
1. Introduction The assessment of thermal performance in many engineering applications requires the estimation of the ground temperature as a function of time, location, and depth. Examples of these applications are buried ground heat exchangers connected to ground source heat pump systems or energy storage systems [1,2]. Earthto-air heat exchangers for heating and cooling of buildings and agricultural greenhouses [3e5] are two other applications. Ground temperature is also important in the design of airport, road pavements, pipelines, as well as for buried high-voltage power cables and nuclear waste disposal facilities. Many factors such as slope orientation, terrain, solar radiation, wind, rain, etc., can affect the ground thermal behavior to a more or less important extent [6] and a reliable assessment approach is highly needed. Since few measured data are available on ground temperature,
* Corresponding author. E-mail addresses: [email protected] (M. Badache), Parham. [email protected] (P. Eslami-Nejad), [email protected] (M. Ouzzane), [email protected] (Z. Aidoun), [email protected] (L. Lamarche).
several models have been developed for its estimation as a function of time and depth. Determination of this temperature is generally based on the solution of the transient one dimensional heat conduction problem in the ground (Equation (1)), and depends on the type of boundary condition at the ground surface. A full discussion of this subject can be found in Ref. [7]. The boundary condition can take two forms depending on the available data [8]. If the ground surface temperature is known, its temporal variation is the boundary condition. If this information is not available, the energy balance equation can be used as boundary condition at the ground surface. The present work has been performed in this latter context. It is worth noting that, both kinds of boundary conditions lead to solutions of the same general form (Equation (2)) but require different inputs in order to determine their respective parameters. For a given surface temperature, the most commonly used model is that presented by Refs. [9,10]. More details about this equation including all parameters and variables are given in Ref. [9]. This model seems to be very simple to use. However it requires the knowledge of the ground temperature at or near the surface which is relatively scarce and without it this model cannot be used. For the surface energy balance models, meteorological data are used in the boundary condition accounting for all heat gains and
http://dx.doi.org/10.1016/j.renene.2015.06.020 0960-1481/© 2015 Crown Copyright and ELSEVIER Ltd. Published by Elsevier Ltd. All rights reserved.
M. Badache et al. / Renewable Energy 85 (2016) 436e444
Nomenclature
d u f
annual damping depth, (m) emissivity angular frequency, (radian/day) phase angle, (radian)
F
heat flux, (W/m2)
ε
cp G h k rh t T V Z
specific heat at constant pressure, (J/kg,K) Intensity of solar radiation, (W/m2) surface convective heat transfer coefficient, (W/m2,K) thermal conductivity, (W/m,K) relative humidity of the air time, (day) temperature, (K) wind velocity, (m/s) depth of the ground, (m)
Greek symbols thermal diffusivity, (m2/s) absorption coefficient StefaneBoltzmann constant, (W/m2,K4)
a as s
losses over the ground surface [11]. were among the earliest who adopted an energy balance model to investigate the daily and annual variation of the ground temperature in Kuwait. They developed an expression of the ground temperature based on the periodic variation of solar radiation and air temperature, making use of the sol-air temperature definition. Latent heat exchange due to evaporation has not been taken into account in their energy balance equation [12]. reproduced the analysis of [11] by introducing the sol-air-evaporation temperature definition. A parametric study was carried out to investigate the effect of various parameters (air relative humidity, ground absorptivity, wind speed, and evaporation fraction) on the heat transfer in the ground [13,14] presented two other models, very similar to that of [11]. They added another parameter to the energy balance equation which is the latent heat exchange due to evaporation. A drawback of the aforementioned energy balance models is that their performances depend strongly on the accuracy of the weather data and the estimation of the ground input factors. Many of these input factors are determined by empirical correlations (e.g., convective heat transfer coefficients), some are experimental constants and others are effective mean values (e.g., thermal conductivity and diffusivity of the ground). In addition, while the long wave radiation (Fsky ) depends on soil radiative properties, air relative humidity and effective sky temperature it was assumed in the aforementioned models [12e14] to be constant and equal to 63 W/m2 for all locations in the world. While this value fits conditions on which it is selected and similar conditions elsewhere, it becomes a rough approximation when different conditions have to be considered. The objective of this paper is to develop an improved model that predicts the ground temperature as a function of depth and time, based on the meteorological data. This model is less affected by the inaccuracy of input factors estimation. The energy balance equation is used as a boundary condition at the ground surface in order to determine the amplitude and the phase angle of the temperature at this location while the annual average ground surface temperature is calculated using an empirical correlation. Unlike previous energy balance models, besides assuming a periodic variation of solar radiation and ambient air temperature, this model assumes a periodic variation of the sky temperature.
2. Model development The equation describing the transient, one dimensional heat conduction in the ground is:
437
Subscripts Air air amb ambient dp dew point s ground surface sky sky conv convection evap evaporation
dT dT 2 ¼a 2 dt dz
(1)
The ground is assumed to be a semi-infinite medium with constant physical properties. Considering the following initial and boundary conditions: Tðz; 0Þ ¼ T s at t ¼ 0 Ts ¼ Tð0; tÞ ¼ T s As cosð ut fs Þ at z ¼ 0 Tð∞; tÞ ¼ T s at z ¼ ∞ The solution of this problem is given by Equation (2) [8,9].
z z Tðz; tÞ ¼ T s As :exp d cos ut fs d
(2)
This form is used in several software: TRNSYS (2005) [15], DOE2 (1982) [16] and RETScreen (2005) [17].where: T(z,t) is the ground temperature at time t (day or hour) and depth z (m). T s is the annual average ground surface temperature, equivalent of the undisturbed ground temperature. As is the annual amplitude of the ground surface temperature. d is the damping depth (m) of annual fluctuation of the ground temperature which can be calculated from:
d¼
rffiffiffiffiffiffi 2a u
(3)
where a is the thermal diffusivity and u ¼ 2 p/365 (radian/days) is the angular frequency.
fs is the phase angle (radian). At the ground surface (z ¼ 0) the energy balance is applied as a boundary condition (Equation (4)) in order to calculate the parameters T s ,As,4s. The details of the energy balance are given in Fig. 1. This figure represents the main heat flux contributions at the surface of the ground including the conduction heat flux into the ground, the convective heat flux transferred between the surface and the ambient air, the short-wave global solar radiation (radiant flux) absorbed by the ground surface, the long wave radiant flux exchanged with the surroundings (sky), and the latent heat flux of evaporation.
438
M. Badache et al. / Renewable Energy 85 (2016) 436e444
humidity. It can be evaluated from the following equation [13,18]:
Fevap ¼ c:f :hconv ½ðaTs þ bÞ rh ðaTamb þ bÞ
(8)
where a ¼ 103 (Pa/K), b ¼ 609 (Pa) and c ¼ 0.0168 (K/Pa) rh is the relative humidity of the air, while f is a fraction of evaporation rate, varying between 0 and 1, and depending mainly on the ground cover and moisture. For bare soil, the fraction f can be estimated as follow: For saturated soil, f ¼ 1, for moist soil, f ¼ 0.6e0.8. For dry soil, f ¼ 0.4e0.5, for arid soil, f ¼ 0.1e0.2. By substituting Equations (5)e(8) into Equation (4), the boundary condition is written as follow:
hr Tamb þ as G þ hrad Tsky c:hconv f :b:ð1 rh Þ ¼ he Ts k
vT vz z¼0 (9)
where hr ¼ hconv(1 þ c . a.f.rh) and he ¼ hconv(1 þ c .a .f)þ hrad Fig. 1. Main heat flux contributions at the surface of the ground.
2.1. Fourier series approximation
The sign convention used in the energy balance assumed that heat transfer to the ground surface is positive, while heat transfer away is negative.
Fnet ¼ k
vT ¼ Fconv þ as G Fsky Fevap vz z¼0
(4) Tðz; tÞ ¼ T s þ Real
where Fnet represents the net heat flux transferred by conduction across the ground surface. The equations described below provide the necessary relationships representing each of the terms in Equation (4). Fconv is the sensible heat flux exchanged between the air and the ground surface, expressed as:
Fconv ¼ hconv ðTamb Ts Þ
(5)
where hconv is the convective heat transfer coefficient. asG is the short-wave global solar radiation (radiant flux) absorbed by the ground surface. Fsky is the long-wave radiation (radiant flux) emitted by the ground surface to the sky and calculated by Equation (6). Since air temperature and solar radiation, show a cyclic behavior throughout the day and the year, it is reasonable to assume that the sky temperature also follows such a cyclic pattern. Therefore, a sinusoidal function can be used to correlate the sky temperature with the meteorological data. This representation of the sky temperature gives a more accurate evaluation of Fsky in comparison to that given by a constant value.
Fsky ¼ hrad Ts Tsky
with
For the sake of simplicity, the solution of Equation (1) is expressed using the complex form of Fourier series (Equation (10)). This expression can be applied to the surface as well as to the various depths. All assumptions made above have been considered herein.
hrad ¼ εs s Ts þ Tsky
N X
1þi where g ¼ 0 and d0 ¼ d
2 Ts2 þ Tsky
The term hrad represents the constant part of the linearized term, of short wave solar radiation, εs is the emissivity of the ground surface and s is the StefaneBoltzmann constant (5.670 108 W/ m2 K4). Fevap is the evaporation heat exchange flux that depends on several factors including ground surface and ambient air temperatures, wind speed, surface cover, soil moisture content and air
rffiffiffiffiffiffi 2a nu
(11)
equal to, Arctan (an/bn), and an and bn are Fourier coefficients for the wave temperature given as:
an ¼
365 2 X TðtÞsinðnutÞ 365 j¼1
(12)
bn ¼
365 2 X TðtÞcosðnutÞ 365 j¼1
(13)
d0 is the nth harmonic damping depth and ‘’Real’’ denotes the real part of the complex number. The ground surface temperature Ts (0,t) in the complex form can be written as: Ts ð0; tÞ ¼ T s þ Real
(7)
(10)
where n is the number of the harmonic, t is the time (number of the day), As is the amplitude of the nth harmonic and is equal qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi to a2n þ b2n , fn (radians) is the nth harmonic phase angle and is
(6)
As einutzgfn
n¼1
N X
As eifs einut
(14)
n¼1
The complex form of Fourier series is also used to express the periodic behavior of the solar radiative flux, sky and air temperature as follows:
Tamb ðtÞ ¼ T amb þ Real
N X n¼1
Aa eifa einut
(15)
M. Badache et al. / Renewable Energy 85 (2016) 436e444
GðtÞ ¼ G þ Real
N X
AG e
ifG inut
e
(16)
n¼1
Tsky ðtÞ ¼ T sky þ Real
N X
Asky eifsky einut
(17)
n¼1
By substituting Equations 10e17 into Equation (9) and equating the time dependent parts, the amplitude of the ground surface temperature, As, and the phase angle difference between the air and the ground surface temperature, fs- fa in Equation (2) are given by:
h A þ a A eiðfG fa Þ þ h A eiðfsky fa Þ s G r a rad sky As ¼ ðhe þ k: g Þ "
hr Aa þ as AG eiðfG fa Þ þ hrad Asky eiðfsky fa Þ fs fa ¼ Arg ðhe þ k:g Þ
(18) #
(19) In Equations (18) And (19) the symbols kk and Arg are respectively used to indicate the modulus and the argument of a complex number. By equating the time independent parts, the annual mean ground surface temperature can be calculated by the following equation:
Ts ¼
i 1h hr T amb þ hrad T sky þ as G c:hconv :f :bð1 rh Þ he
(20)
The Equations (18)e(20) provide the necessary parameters of Equation (2) to calculate the ground temperature as a function of depth and time with respect to the energy balance equation applied at the ground surface. These equations are different from those of [13,14] because they include more realistic effects through additional terms which are fsky, Asky, hrad, and Tsky. 2.2. Model sensitivity to factors estimation Many input factors such as hconv, as, εs, f and rh are involved to determine the parameters T s ,As, 4s in Equations (18)e(20). While rh and f are relatively simple to estimate, there is a great number of available correlations in the literature for hconv, as, εs. For the absorption factor, as, different correlations relating this factor to the surface water content are used [19,20]. The emissivity of the soil surface, 3s, generally depends on the water content and vegetation of the surface, the commonly used correlations are [21e25] reviewed a large number of convective heat transfer coefficient correlations in linear, power law and boundary layer form. Generally hconv correlations depend on the wind velocity and several specific conditions under which they have been correlated. As there is a large number of correlations available it is difficult to assess the impact of all correlations factors on the present model predictive performance. As a result, the following analysis stresses the impact of convective heat transfer coefficient values using three representative correlations on the ground temperature profile. The ground characteristics of a site located in Varennes are taken as the base case example for calculations. Equation (2) is used to predict the ground temperature for two different days (21st march and 21st September) at different depths. Correlations of Givoni and Mostrel [26], Kasuma [27] and Mc Adams [28] have been selected for the convective heat transfer coefficient calculations. Fig. 2 represents the ground temperature profile, calculated using the aforementioned correlations. The difference between the calculated values
439
using Kasuma and Mc Adams correlations are relatively insignificant. However, there is a significant difference between the values calculated using Givoni and the other two correlations (Kasuma and Mc Adams). It shows also that Givoni gives the highest ground temperature prediction for the shallow (1 me8 m), and the deep zone (below 8 m). This result can be explained by the fact that the convective heat transfer coefficient may be underestimated by Givoni and Mostrel correlation (3 W/m2.K) or overestimated by Kasuma correlation (47 W/m2 K). An underestimation of the convection heat transfer coefficient overestimates the ground temperature. The opposite is true when overestimating the convective heat transfer coefficient (Fig. 2). In summary, the above analysis clearly shows that, an inaccurate estimation of the convective heat transfer coefficient, due to inappropriate choice of the correlation factors results in inaccurate model parameters (T s , As, fs) and ground temperature prediction. Therefore it is important for designers and modelers to select reliable convective heat transfer coefficients correlations for T s, As and fs computations. The same remark holds true for as and εs factors. 2.3. Model enhancement features As stated in Equation (2), the present model is a function of three main parametersT s , As and fs which are strongly dependent on several input variables. The importance of accurate estimation of input variable has been stressed upon in the previous section. It is the purpose of this section to improve the model in order to develop an accurate and easy to use correlation for ground temperature profile. Equation (2) implies that the undisturbed ground temperature is constant, while the amplitude and phase angle are damped with a factor of e(z/d). Therefore the potential inaccuracy in As and fs decreases very rapidly in the first few meters from the down from the ground surface, while inaccuracy in T s remains unchanged with depth. This means that, calculating accurate values of T s is crucial for a correct model temperature prediction. It is obvious that accurate values of input factors result in accurate calculation of T s . Unfortunately, as shown above it is not always straightforward because of variety of available correlations. A solution to this problem is to use a form for T s that is independent of the aforementioned input factors. Recently a new correlation (Equation (21)) has been developed by Ref. [29] for calculation of T s . It is only a function of the mean annual air ambient temperature
Fig. 2. Ground temperature prediction using three different correlations for the calculation of the convective heat transfer coefficient.
440
M. Badache et al. / Renewable Energy 85 (2016) 436e444
(available factor in all meteorological databases). This correlation can be used in Equation (2) instead of Equation (20). Therefore, the improved ground temperature model proposed in this paper is based on Equations (18), (19) and (21).
T s ¼ 17:898 þ 0:951T amb
(21)
where T amb is the annual average ambient air temperature in degree Kelvin. 3. Results and discussion In this section, first the developed model is applied to predict the daily ground temperature of the Varennes' (Montreal region) site for various depths and time. The model predictions are compared with experimental results and the profile of Varennes ground temperature is then analysed. Further, the model is validated for two other sites in the U.S.A: Fort Collins (Colorado) and Temple (Texas). The following is a description of the procedure used for ground temperature prediction of Varennes using this model. The data of a typical meteorological year (TMY) [30] of Montreal (global solar radiation (G), ambient dry bulb (Tamb), dew point temperatures (Tdp), and wind velocity (Vw)) are used for the determination of the model inputs, including the mean values (T amb , G, T sky ), the phase angles (4a,4G,4sky), and the amplitudes (Aa AG Asky ). Fourier analysis (FA) is used to extract these parameters. The calculations are limited to the base case of the first harmonic (n ¼ 1). The first harmonic accounts for 95% of the variance in the worst cases. Table 1 gives the calculated mean values, the phase angles, and amplitudes of the parameters G, Tamb and Tsky for the site of Varennes. The sky temperature is calculated by Equation (22) [31].
h i0:25 2 Tsky ¼ Tamb 0:711þ0:0056Tdp þ0:000073Tdp þ0:013 cosð15tÞ (22)
was fractured black shale from 0 m to 19.5 m, slate from 19.5 m to 23.2 m and black shale from 23.2 m to 31 m. Twenty four 24 calibrated thermocouples (T type) were attached on the surface of both legs of the U tube at different depths. The thermocouples were used to measure the ground temperature profile from February to November 2014. Because of the important variation of the temperature near the ground surface, more thermocouples were located in this region (Fig. 3) while beyond 11 m deep they were fewer and equally spaced. Temperature data were recorded every 5 min, using data acquisition system (DAQ) and then compiled into the daily-average format. Daily-average formats are calculated by integrating 5 min data over the period (24 h) using an in-house Fortran code. The acquisition system was fitted with two extensions of twenty entries each and communications ports including Ethernet and USB. The system used for the temperature measurement is composed of five independent sources of error: the DAQ, the temperature sensor itself, the cold junction compensation (CJC) sensor, thermocouples position, and cabling. The uncertainty on the temperature measurements is ±0.42%. It is noted that the error due to thermocouples position and cabling is ignored so much it is so small (3 105). The following observations were made: The thermocouples at the same level in both borehole legs give the same temperature measurements. Therefore only one leg measurements are considered. The daily mean temperature curves (at depths less than 3 m) against the day number vary approximately in a non-sinusoidal form. This can be due to many parameters such as short term weather variations, seasonal variations, ground water content, and non-homogeneous ground thermal properties. The amplitude of the ground temperature oscillation reduces from around 15 C at 1 m to 0.17 C at 31 m. The shift in phase angle gradually increases with depth. Temperature measurements show that at depths greater than about 8 m below the surface, the ground temperature remains relatively constant (9.5 C).
where Tsky and Tamb are in degree Kelvin and Tdp is in degrees Celsius. Parameter t is the hour from midnight. Calculation of As and 4s parameters (from Equations (18) and (19)) involves the determination of f, εs, as, d0 , hconv, and hrad. Parameter d0 is a function of the thermal properties of the ground and the frequency of the thermal oscillation. Effective ground thermal properties obtained from a thermal response test (TRT) report are k ¼ 2.65 W/m$K and a ¼ 0.0948 m2/day. The TRT test was carried out by a contractor to determine these properties based on ASHRAE standards with less than 5% inaccuracy. Recent guidelines for the test procedure can be found in Ref. [32]. It is assumed that the coefficients hconv, f, rh, εs, as do not vary during the year and their annual average values are used. The coefficients f, εs, and as are assumed to be 0.7, 0.9, and 0.9, respectively. Therefore the annual value of d0 calculated from Equation (3) is 3.31 m. In the case of Montreal (cold climate) [33], suggest the use of the McAdams correlation, in order to estimate hconv. The mean annual wind velocity is given to be 4.5 m/s. Since hrad is a function of ground surface temperature (Equation (7)), an iterative method is used to calculate the daily values of this factor. 3.1. Experimental set up in Varennes site A U-tube borehole (31 m deep) was drilled at CanmeENERGYVarennes. This set-up was fitted with spacers which were placed in the center of the U-tube to keep the tube legs apart. Super Grout was used for the filling material. The soil stratigraphy encountered
Fig. 3. Position of the thermocouples in the borehole.
M. Badache et al. / Renewable Energy 85 (2016) 436e444
441
Table 1 Mean annual value, phase angle, and amplitude of the meteorological parameters G, Tamb and Tsky for site of Varennes. Parameter
Global irradiance (G)
Dry bulb temperature (Tamb)
Sky temperature (Tsky)
Amplitude Phase angle Mean annual value
160.57 0.28 103.12
15.59 0.35 6.17
21.33 0.37 15.25
3.2. Comparison with experimental measurements Figs. 4 and 5 show the comparison between the measured and predicted ground temperatures at eight different depths during the period from February 4th to November 25th, 2014. Generally, the agreement between the predicted and the measured values is good. However, the predicted ground temperature values appear to be overestimated in the spring time (day 50e150) as clearly shown for depths 1 and 2 m below the soil surface. This temperature difference vanishes gradually with depth. It can be seen also that the measured ground temperature variation does not quite take a sinusoidal form. This may be due to the nonhomogeneous thermal properties of the ground. Thermal properties change with time and depth due to the volume increase of the ground water content during the thaw period (spring time). In the present study however, properties are assumed to be constant over time and location, as measured in November 2010. During the thaw period the thermal diffusivity and therefore the damping depth increases due to the increase of ground water content. At any given depth, the relative value of the ground temperature amplitude is proportional to the damping depth as shown in equation (23). Consequently, an increase of damping depth brings the surface and deep ground temperatures closer together, thus resulting in a slight drop in ground temperature for this period of the year (day 50e150). ðzz0 Þ= d
Az=A ¼ exp s
(23)
To provide a more quantitative basis for this comparison we used a statistical analysis. The Root Mean Square Error (RMSE) and the Normalized Root Mean Square Error (NRMSE) statistics are calculated. RMSE was generally 2.3 K in the worst cases, i.e. at depth of 1 m and less than 1.5 K for depths below 2 m from the soil surface. NRMSE is less than 0.92%. Despite the fact that the data were measured over approximately ten months and the RMSE is expected to be lower for an entire year, the results are very satisfactory. In other words, the present model predicts quite accurately the ground temperature profile over time and depth.
3.3. Analysis of the ground temperature profile in Varennes Figs. 4 and 5 above clearly show the damping of the annual ground temperature in which the amplitude reduces from around 15 C at 1 m to 0.17 C at 26 m. Equation (2) shows that the amplitude decreases exponentially with distance from the surface at a characteristic length known as the damping depth, d. The damping depth is the depth at which the amplitude of the temperature is 1/e times the amplitude of the temperature at the surface. The mean value of the annual damping depth of Varennes is estimated to be 2.35 m from the measured ground temperature. According to this value the amplitude of the ground temperature is almost completely (95%) damped at 3d (i.e. at 7 m). This is confirmed in Fig. 5 (b), (c) and (d), where the ground temperatures are essentially constant throughout the year for depths below 6 m.
Fig. 4. Measured and predicted ground temperature at depths 1 m, 2 m, 3 m, and 4 m at Varennes for the period from February 4th to November 25th, 2014.
442
M. Badache et al. / Renewable Energy 85 (2016) 436e444
Fig. 5. Measured and predicted ground temperature at depths 6 m, 8.5 m, 16 m, and 26 m at Varennes for the period from February 4th to November 25th, 2014.
These figures show also that below the depth of 1 m, the maximum or minimum temperature occurs later than it does at 1 m, meaning that the time lag increases with the depth. For example the ground temperature is at its minimum on days 38e40 (7th to 9th of February) for the depth 1 m while it is the lowest 2.5 months later (end of April) for the depth of 4 m. Finally, it is interesting to analyse the model predictions on a daily basis. Therefore the profiles of the ground temperature for two days (April 10th and September 2nd) are presented in Fig. 6. This figure shows that due to the temperature fluctuations on the ground surface, the temperature profiles are different. For each day, two different regions are observed; near the ground surface (up to 8.5 m) and deep down to 31 m. The ground temperature profile
Fig. 6. Measured and calculated ground temperature at different depthat Varennes for the April 10th and September 2nd, 2014.
reaches the same value for both days, which corresponds to the undisturbed ground temperature. This confirms our previous observation in Fig. 5c and d. The undisturbed ground temperature value, predicted to be 10.5 C, is 1 C higher than the experimental value, which is about 9.5 C.
3.4. Comparison with additional experimental data In order to extend the application of the new model to different climate conditions, it was applied for the sake of example to predict monthly average ground temperature of two sites in the U.S.: Fort Collins (Colorado) and Temple (Texas). Temple has a warm humid temperate climate with hot summers while Fort Collins has a humid continental climate. The ground temperature data for these two locations are given in the reference article of Kusuda et al. (1965). The mean, the amplitude and phase angle for solar radiation, sky and ambient air temperatures are determined using FA method from TMY data (Table 2). The FA coefficients are calculated for the first harmonic (n ¼ 1). The thermal properties of the ground at Fort Collins and Temple are not available but they can be estimated from the measured temperature data, Equations (3) and (23). Equation (23) is used to estimate the damping depth between 0.0254 m and 1.83 m Equation (3) is used to calculate the corresponding thermal diffusivity. Table 3 shows the parameters used to predict the ground temperature of the two sites. The data (hconv, rh, f, k, a, as, εs) of Fort Collins are obtained from Ref. [13] while those of Temple are from Ref. [30]. Equations (18) (19) and (21) are used to calculate the monthly ground temperature for Fort Collins at 0.1 m and 1.83 m, and for Temple at 0.1 m and 1.23 m from the surface. Figs. 7 and 8 provide the measured and predicted monthly ground temperature for Fort Collins at 0.1 m and 1.83 m, and for Temple at 0.1 m and 1.23 m from the surface. Despite all potential uncertainties that may be accounted for assuming some factors
M. Badache et al. / Renewable Energy 85 (2016) 436e444
443
Table 2 Calculated mean annual value, phase angle, and amplitude of the meteorological parameters G, Tamb and Tsky for Fort Collins, Colorado and Temple, Texas. Factor
T amb (ºC)
fa (rad)
Aa (Rad)
G (W/m2)
fG (rad)
AG (W/m2)
T sky (ºC)
fsky (rad)
Asky (ºC)
Fort collins (Colorado) Temple (Texas)
6.14 18.8
0.23 0.23
14.5 11.15
157 86.5
0.13 0.04
132 202.1
14.9 2.65
0.23 0.23
16.1 16.4
Table 3 Parameters used to predict Fort collins and temple ground temperature. Factor
hconv (W/m2.K)
rh
f
ɑs
ks (W/m.K)
ɑ (m2/day)
d (m)
3s
Fort collins (Colorado) Temple (Texas)
22.0 23.6
0.55 0.70
0.6 0.3
0.8 0.8
1 1
0.0362 0.0279
2.48 2.10
0.9 0.9
used as a boundary condition at the ground surface in order to determine the amplitude and the phase angle of the temperature at this location. The mean ground surface temperature is then calculated using an empirical correlation. Unlike the previous energy balance models, this model is less sensitive to input factors inaccuracy. Furthermore, besides assuming a periodic variation of the solar radiation and ambient air temperature, it also assumes a periodic variation of the sky temperature. Validation against measured data of a site located in Varennes (Montreal-Canada) and two sites in the U.S.A; Fort Collins (Colorado) and Temple (Texas) was performed. Generally the agreement between measurements and the model predictions is good. Therefore, the developed model is proposed as useful tool for the ground temperature profile estimation in the context of various environmental and energy applications. Acknowledgments Fig. 7. Predicted and measured values of ground temperature of Fort Collins (Colorado) at 0.1 m and 1.83 m from the surface.
values, the model still predicts adequately the ground temperature variation for both sites. 4. Conclusion A new model was developed to improve the prediction of the temperature variation in the ground as a function of depth and time based on the meteorological data. The energy balance equation is
Fig. 8. Predicted and measured values of ground temperature of Temple (Texas) at 0.1 m and 1.22 m from the surface.
Financial support for this work was provided by the Canadian Federal Government's Program on EcoENERGY Innovation Initiative (Natural Resources Canada) (EEBI 022). References [1] G. Florides, S. Kalogirou, Ground heat exchangers-A review of systems, models and applications, Renew. Energy 32 (2007) 2461e2478. [2] H. Esen, M. Inalli, M. Esen, K. Pihtili, Energy and exergy analysis of a groundcoupled heat pump system with two horizontal ground heat exchangers, Build. Environ. 42 (10) (2007) 3606e3615. [3] G. Mihalakakou, J.O. Lewis, M. Santamouris, The influence of different ground covers on the heating potential of earth-to-air heat exchangers, Renew. Energy 7 (1996) 33e46. [4] M. Esen, T. Yuksel, Experimental evaluation of using various renewable energy sources for heating a greenhouse, Energy Build. 65 (2013) 340e351. [5] H. Esen, M. Inalli, M. Esen, K. Pihtili, Numerical and experimental analysis of a horizontal ground-coupled heat pump system, Build. Environ. 42 (3) (2007) 1126e1134. [6] G.P. Williams, L.W. Gold, Ground Temperatures, National Research Council Canada, Institute for Research in Construction, 1976. Available from: http://irc. nrc-cnrc.gc.ca/pubs/cbd/cbd180e.html. [7] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Physics Today, 1962. [8] B. Givoni, L. Katz, Earth temperatures and underground buildings, Energ. Build. 8 (1985) 15e25, http://dx.doi.org/10.1016/0378-7788(85)90011-8. [9] T. Kasuda, P.R. Archenbach, Earth temperature and thermal diffusivity at selected stations in the United States, ASHRAE Trans. 71 (1965). Part1. Available from: http://oai.dtic.mil/oai/oai? verb¼getRecord&metadataPrefix¼html&identifier¼AD0472916. [10] W.R. Van Wijk, D.A. De Vries, Periodic Temperature Variations in an Homogeneous Soil. Physics of Plant Environment, North-Holland Publ, Amsterdam, 1966, pp. p.102ep.143. [11] A.K. Khatry, M.S. Sodha, M.A. Malik, Periodic variation of ground temperature with depth, Sol. Energy 20 (1978) 425e427. [12] M.M. Salah El-Din, On the heat flow into the ground, Renew. Energy 18 (1999) 473e490. [13] M. Krarti, C. Lopez-Alonzo, D.E. Claridge, J.F. Kreider, Analytical model to predict annual soil surface temperature variation, J. Sol Energy Eng. 117 (2) (1995) 91e99.
444
M. Badache et al. / Renewable Energy 85 (2016) 436e444
[14] G. Mihalakakou, M. Santamouris, J.O. Lewis, D.N. Asimakopoulos, On the application of the energy balance equation to predict ground temperature profiles, Sol. Energy 60 (3) (1997) 181e190. [15] TRNSYS. Version 16.a. Transient System Simulation Program Solar Energy Laboratory, University of Wisconsin, Madison, USA, 2005. [16] DOE-2, Technical Information Center. Engineers Manual Version 2.1A, Department of Energy, USA, 1982. [17] RETscreen 4, Clean Energy Project Analysis: RESTscreen Engineering & Cases Textbook, Ground-Source Heat Pump Project Analysis Chapter, CanmetENERGY at Varennes, Natural Resources Canada, 2005. [18] G. Mihalakakou, On estimating soil surface temperature profiles, Energy Build. 34 (3) (2005) 251e259. [19] C.H. Van Bavel, D.I. Hillel, Calculating potential and actual evaporation from a bare soil surface by simulation of concurrent flow of water and heat, Agric. Meteorol. 17 (1976) 453e476. [20] J. Kondo, N. Saigusa, T. Sato, A model and experimental study of evaporation from bare soil surfaces, J. Appl. Meteorol. 31 (3) (1992) 304e312. [21] S. Ortega-Farias, R. Antonioletti, A. Olioso, Net radiation model evaluation at an hourly time step for Mediterranean conditions, Agronomie 20 (2000) 157e164. [22] M. Fuchs, C.B. Tanner, Infrared thermometry of vegetation, Agron. J. 58 (1976) 597e601. [23] K. Noborio, K.J. McInnes, J.L. Heilman, Two-dimensional model for water, heat,
[24] [25] [26] [27] [28] [29]
[30] [31] [32] [33]
and solute transport in furrow irrigated soil: II. Field evaluation, Soil Sci. Soc. Am. 60 (1996) 1010e1021. W. Brutsaert, On a derivable formula for long-wave radiation from clear skies, Water Resour. Res. 11 (1976) 742e744. J.A. Palyvos, A survey of wind convection coefficient correlations for building envelope energy system's modeling, Appl. Therm. Eng. 28 (2008) 801e808. B. Givoni, M. Mostrel, Windscreens in radiant cooling, Passive Sol. J. 1 (1982) 229e238. G. Kusuma, Estimation of the exchange coefficient of heat during low wind convective conditions, Bound. Layer. Meteorol. 111 (2004) 247e273. W.H. McAdams, Heat Transmission, McGraw-Hill, New York, 1954. M. Ouzzane, P. Eslami-Nejad, M. Badache, Z. Aidoun, New correlations for the prediction of the undisturbed ground temperature, Geothermics 53 (2015) 379e384. National Solar Radiation Data Base, web address: http://rredc.nrel.gov/solar/ old_data/nsrdb/1991-2005/tmy3/by_USAFN.html. J.A. Duffie, W.A. Beckman, Solar Engineering of Thermal Processes, third ed., John Wiley & Sons, Inc, USA, 2006. ISBN-13 978-0-471-69867-8. ASHRAE, ASHRAE Hnadbook: HVAC Applications. Chapter 32, ASHRAE Inc., Atlanta, GA, 2007. M. Ouzzane, P. Eslami-Nejad, Z. Aidoun, L. Lamarche, Analysis of the convective heat exchange effect on the undisturbed ground temperature, Sol. Energy 108 (2014) 340e347.