A Numerical Model to Estimate the Soil Thermal Conductivity Using Field Experimental Data

Abstract

Soil thermal conductivity is an important parameter for understanding soil heat transfer. It is difficult to measure in situ with available instruments. This work aims to propose a numerical model to estimate the thermal conductivity from the experimental measurements of soil heat flux and soil temperature. The new numerical model is based on the Fourier Law adding a constant empirical parameter to minimize the uncertainties contained in the data from field experiments. Numerically, the soil thermal conductivity is obtained by experimental linear data fitting by the Least Squares Method (LSM). This method avoids numerical indetermination when the soil temperature gradient or soil heat flux is very close to zero. The new model is tested against the different numerical methodology to estimate the soil heat flux and validated with field experimental data. The results indicate that the proposed model represents the experimental data satisfactorily. In addition, we show the influence of the different methodologies on evaluating the dependence of the thermal conductivity on the soil water content.

1. Introduction

Soil thermal conductivity is an important thermal property that governs the transfer of the soil heat flux. The complexity of soil makes its thermal conductivity dependent on physical or chemical properties, such as soil water content, porosity, density, texture, and mineral composition [1]. Different empirical models have been suggested for the estimation of soil thermal conductivity, taking into account these properties [2,3,4,5,6,7,8]. However, these models are mainly obtained using data from controlled experiments in the laboratory by measurement techniques including steady-state methods such as the guarded hot plate method [9], transient methods including the single line heat source probe [10,11] and the dual-probe heat-pulse method [12,13,14]. A challenge has been to obtain estimates of the soil thermal conductivity in situ (or in field experiments), in which soil properties are not well characterized as a result of technical limitations [15,16] and/or climatic conditions.

Theoretically, in situ experimental measurements of soil heat flux and soil temperature can be used to estimate the soil thermal conductivity by the Fourier Law of heat conduction, defined by [17]:

$$G_z= -\lambda \frac{\partial T}{\partial z}$$

(1)

where $G_z$ (${\mathrm{Wm}}^{-2}$) is the soil heat flux at a given depth $z$ ($\mathrm{m}$), $T$ ($\mathrm{K}$) is the soil temperature and $\lambda$ (${\mathrm{Wm}}^{-1}{\mathrm{K}}^{-1}$) is the thermal conductivity of the soil layer. However, because the data are discrete (measured in time intervals), the soil temperature gradient, $\frac{\partial T}{\partial z}$, can be obtained by the variation rate between the soil temperature and the depth, $\frac{\Delta T}{\Delta z}=\frac{T_2-T_1}{z_2-z_1}$. This methodology is called the Gradient Method [4,5,16,17,18,19] Therefore, if we know the measurement of the soil heat flux at a given depth $z$ ($G_z$) and the soil temperature at depths relatively close to $G_z$, Equation (1) can be written as:

$$G_z= -\lambda \left[\frac{T_2-T_1}{z_2-z_1}\right]$$

(2)

with $z_1<z<z_2$, and where $T_1=T\left(z_1\right)$ and $T_2=T\left(z_2\right)$ represent the temperatures at depths $z_1$ and $z_2$, respectively.

The experimental estimate of soil heat flux in $z$ is often performed by the heat flux plate method. Soil heat flux plates are small, rigid, disc-shaped sensors of known and constant thermal properties that are placed horizontally in the soil near the surface. The heat flux through a calibrated plate is used to estimate $G_z$ in the surrounding soil at the plate depth [20]. That is, Equation (2) is used to estimate the soil heat flux in a small soil layer if the $\lambda$ of the plate is known.

Heat flux plate manufacturers describe important factors to be considered in an experiment to estimate $G_z$, especially in a field experiment in which several procedures are not controlled, such as the exact installation depth and the thermal contact resistance at the plate-soil interface [21,22,23]. Soil heat flux sensors are often installed at depths ranging from $0.05 \mathrm{m}$ to $0.20 \mathrm{m}$, considered to be mechanically stable to ensure good measurement conditions besides not suffering significant time delays that can reduce the accuracy of the measurements [24,25,26]. Ochsner et al. [15] concluded that the heat flux plate underestimates soil heat flux because of the low plate thermal conductivity, thermal contact resistance, and latent heat transfer effects. However, this is the main method to obtain the experimental soil heat flux and widely used in studies of surface energy balance closure [27,28,29].

Using experimental measures of soil heat flux and soil temperature, several authors estimated the soil thermal conductivity by isolating $\lambda$ in Equation (2). However, this procedure can lead to possible numerical divergences, as shown by Li et al. [19,30] that eliminated data when soil temperature gradient or soil heat flux is very close to zero to avoid this problem.

In this work, we propose a new numerical model to estimate thermal conductivity from experimental measurements of soil heat flux and soil temperature. The model uses the Fourier Law adding a constant empirical parameter ($\epsilon$) to minimize the uncertainties contained in the data from field experiments. Numerically, the solution model is a linear data fitting of the experimental data using the Least Squares Method (LSM). This methodology avoids numerical indetermination and possible data loss. The soil data collected in a field experiment on native vegetation of the Brazilian Pampa biome will be used to calibrate and validate the proposed model.

This work is divided into three sections: Materials and Methods; Results and Discussions; and Conclusions. The Materials and Methods section presents the description of the experimental site, the proposed model, the justification for the proposition of the model and the method used to estimate $\lambda$, in addition to the analysis sequence. The Results and Discussions section presents the results of the proposed model in comparison to the different models used in here. This section also shows the influence of the different methodologies on evaluating the dependence of the thermal conductivity on the soil water content. Finally, the Conclusions present the pertinent reflections on the results obtained.

2. Materials and Methods

2.1. Experimental Site

The experimental measurements were obtained in an area of native vegetation designated for cattle raising, in the city of Pedras Altas, Rio Grande do Sul state, Brazil, in the geographical location ${31}^{°}{43}^{\prime }{56}^{″}\mathrm{S}$, ${53}^{°}{43}^{\prime }{36}^{″}\mathrm{W}$, $395 \mathrm{m}$, located in the Pampa biome region (hereafter PAS). The soil in the study area is characterized as sandy loam with $59.30\%$ sand; $0.81\%$ clay; $39.89\%$ silt. The average values of the soil physical properties at $0.05 \mathrm{m}$ and $0.15 \mathrm{m}$ depths is: field capacity ${\theta }_{FC}=0.28 {\mathrm{m}}^3{\mathrm{m}}^{-3}$; permanent wilting point ${\theta }_{WP}=0.03{\mathrm{m}}^3{\mathrm{m}}^{-3}$; soil porosity ${\theta }_s=0.42 {\mathrm{m}}^3{\mathrm{m}}^{-3}$; and, soil bulk density ${\rho }_s=440 {\mathrm{kg}\mathrm{m}}^{-3}$. More details about soil composition and experimental site are described in [31].

The typical climate in the region is classified as humid subtropical Cfa, according to Köppen [32]. The mean annual temperature is of 17 °C, presenting negative temperatures at morning in the winter (~−3 °C) and almost 40 °C in the afternoon in the summer. The precipitation regime is well distributed throughout the year.

Soil temperatures were measured with temperature sensors (T108, Campbell Scientific Inc., Logan, Utah, USA) positioned at $0.05 \mathrm{m}$ and $0.15 \mathrm{m}$ below the surface, and the soil heat flux was measured with a flux plate (HFP01, Hukseflux Thermal Sensor B.V., Delft, The Netherlands) positioned at $0.10 \mathrm{m}$ below the surface of the soil, so that it was halfway between the soil temperature sensors. The soil water content ($\theta$) was measured with a water content reflectometer (CS616, Campbell Scientific Inc., Logan, UT, USA) installed at $0.10 \mathrm{m}$ below the surface. Data were collected during the year 2015 at a frequency of 1 min and then averaged every 30 min. Due to technical problems, partial data collection occurred during the period from January 01 to 12 (only the soil temperatures and water content were collected), and there was no data collection from September 17 to November 05. If the measurement system failed for one variable, all the other variables were removed from the data set.

2.2. Proposed Model

The relationship between $G_z$ and $\frac{\Delta T}{\Delta z}$, represented by Equation (2), is described by a linear function of the variables, with $\lambda$ being the proportionality constant. However, experimentally, we do not obtain a relation defined by a linear function described by $y=ax$; instead, we obtain a relation defined by a function described by $y=ax+b$, as shown in Figure 1, independent of the soil water content. This figure illustrates the relationship between the $G_z$ measured at $z=-0.10 \mathrm{m}$ and the variation rate between the soil temperature and the depth $\left[\frac{T_2-T_1}{z_2-z_1}\right],$ with $z_1=-0.05 \mathrm{m}$ and $z_2=-0.15 \mathrm{m}$, does not satisfy Equation (2); since the graph $G_z\times \frac{\Delta T}{\Delta z}$ does not pass through the origin, there will be an underestimate or overestimate of $G_z$ with respect to $\left[\frac{T_2-T_1}{z_2-z_1}\right]$ for values negative or positive of $G_z$, respectively. We assume that this behavior may be related to the finite difference method and the experimental measurements of T and/or $G_z$. Therefore, to estimate $\lambda$ from the experimental data on the soil heat flux ($G_{exp}$) and the soil temperature at different depths, we propose to include an empirical parameter ($\epsilon$) in Equation (2). The new model for $G_z$ can be expressed as:

$$G_z =-\lambda \left[\frac{T_2-T_1}{z_2-z_1}\right]+\epsilon$$

(3)

where $\epsilon$ is a constant empirical parameter to be determined, representing a correction factor of an energy flux that is not accounted for by the field experimental measurements. The addition of $\epsilon$ becomes the model able to fit the observed behavior in the experimental data (Figure 1).

In analytical solutions of the Fourier Law, although not explicitly, this correction is also present, for example in [17,33,34]. In these cases, the estimation of this correction is quite complex, requiring a greater number of parameters obtained, in general, from adjustments to experimental data of soil heat flux and soil temperature.

Data Fitting: Least Squares Approximation

Figure 2 shows the daily cycle for the soil heat flux ($G_z$) at $z=-0.10 \mathrm{m}$ and the temperatures at $z_1=-0.05 \mathrm{m}$ and $z_2=-0.15 \mathrm{m}$ for four days in different seasons (Summer DOY 21; Autumn DOY 111; Winter DOY 202; and Spring, DOY 325). Precipitation of $0.1 \mathrm{mm}$ occurred in DOY 325 at 01:00 AM (local time). These days were chosen to emphasize the variability in the amplitude of $G_z$ and soil temperature values and the hours when the soil temperature gradient is zero ($T_{z_1}\approx T_{z_2}$), in different seasons. It may be noted that at hours when the soil temperatures approach each other at different depths, the $\lambda$ obtained from Equation (2), by the division of $G_z$ by $\left[\frac{T_2-T_1}{z_2-z_1}\right]$, can result in nonfinite values or lead to numerical divergences resulting in a divergent solution.

The daily average of $G_{exp}$ and $T_{soil}$ for the different depths can show variability throughout the year, as shown in Figure 3. The existence of days that have low-temperature variability at different depths or values of $G_{exp}$ near zero can also cause numerical divergences, especially when isolating $\lambda$ in Equation (2) [19,30]. In Figure 3, the daily soil water content ($\theta$) is also shown, presenting a high variability throughout the year, with small values in summer and autumn seasons.

To avoid possible divisions by zero using Equation (2), different works in the literature estimate the soil thermal conductivity using, for example, daily averages [4,18] or specific times [18]; some even eliminate data in which the temperatures are approaching zero [19]; and, others eliminate soil heat flux data near signal changes (before and after) because these values are relatively small and consequently increasing errors [19,30]. In contrast to these techniques, in this work, we propose the use of the Linear Least Squares Method (LLSM) to obtain the coefficient of soil thermal conductivity and the additional term ($\epsilon$) presented in Equation (3). This technique allows the estimation of $\lambda$ without inverting Equation (2). This procedure avoids the exclusion of data, allowing a more realistic estimation of the soil thermal conductivity.

2.3. Numerical Models Calibration and Validation

Data from $G_{exp}$ and $T_{soil}$ were selected for different depths on 3 random days (an entire day of data) of each month of the year, totaling 33 days (since during October, data were not collected), and were used to estimate $\lambda$ by different methods. The following methods were used:

• Method 1 ($M1$): Isolating $\lambda$ in Equation (2) and calculating the average value of $\lambda$;
• Method 2 ($M2$): Isolating $\lambda$ in Equation (2) and calculating the average value of $\lambda$ for the observed data whose temperature variation ($\Delta \mathrm{T}$) is equal to or greater than $\left|0.1\right|$;
• Method 3 ($M3$): Isolating $\lambda$ in Equation (2) and calculating the average value of $\lambda$ for the 13-h data (local time). This time was chosen because it presents, in general, a greater variation of the temperature in relation to the depth than the other times;
• Method 4 ($M4$): Least squares method applied to Equation (2);
• Method 5 ($M5$): Least squares method applied to Equation (3).

The method validation was performed using the experimental dataset, excluding the days used to estimate $\lambda$. Analyses were performed considering the whole period (day and night) and separating the periods into day and night. For this, the daytime period was considered as the interval between 06:00 and 18:00 and the night, the interval between 18:00 and 06:00.

The $\lambda$ values obtained by each calibration were entered into the respective equations of each method together with the experimental values of the soil temperature, thereby obtaining estimated $G_z$ values that were used in the different methods to estimate $\lambda$ ($G_{M1}$, $G_{M2}$,…). The method validation was performed using the values of $G_z$ with the values obtained experimentally ($G_{exp}$) through the statistical indices RMSE and $R^2$. These indices were calculated by the following equations:

$$RMSE= \sqrt{\frac{{\displaystyle \sum }{\left(Y_{exp}-Y_{est}\right)}^2}{n}}$$

(4)

$$R^2=1-\frac{{\displaystyle \sum }{\left(Y_{exp}-Y_{est}\right)}^2}{{\displaystyle \sum }{\left(Y_{exp}-\overline{Y}\right)}^2}$$

(5)

where $Y_{exp}$ represents the observed values, $Y_{est}$ represents the estimated values, $\overline{Y}$ represents the mean of the observations and $n$ represents the number of observations.

3. Results and Discussions

Table 1. Values of $\lambda$ and $\epsilon$, obtained by the different methods and respective statistical analyses. Obs.: The Daily column represents the statistics for the entire period (Day and Night).

	Calibration		Validation
	$\mathit{\lambda }$ $\left(\mathit{W}{\mathit{m}}^{-1}{\mathit{K}}^{-1}\right)$	$\mathit{\epsilon }$ $\left(\mathit{W}{\mathit{m}}^{-2}\right)$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$$\left(\mathit{W}{\mathit{m}}^{-2}\right)$			${\mathit{R}}^2$
			Daily	Day	Night	Daily	Day	Night
Method 1 ($M1)$	1.331	-	5.56	7.31	2.92	0.82	0.71	0.77
Method 2 ($M2)$	1.148	-	3.86	5.08	1.97	0.92	0.86	0.89
Method 3 ($M3)$	1.014	-	3.24	3.97	2.28	0.94	0.91	0.86
Method 4 ($M4)$	1.011	-	3.23	3.95	2.30	0.94	0.91	0.85
Method 5 ($M5)$	0.993	-2.216	2.43	3.06	1.58	0.97	0.95	0.93

shows the soil thermal conductivity and the empirical parameter calibrated by the different methods with the experimental data of the Pedras Altas-RS site. The results of $\lambda$ obtained using Equation (2) by the different methods ($M1$, $M2$, $M3$, and $M4$) show variability of almost 25%, with a reduction of the RMSE of approximately 40% between $M1$ and $M4$ and an increase of $R^2$ of approximately 12% for the daily period.

From all analyzed methods and periods, $M1$ presented the worst statistical indices. This result is perhaps due to the fact all data were used in the simulation, which can generate numerical divergences when $\left|\Delta T\right|\approx 0$, consequently, on average, can overestimate $\lambda$. For $M2$ were eliminated data when $\left|\Delta T\right| \le 0.1$, representing $5\%$ of total data and resulting in a decrease of 30% in RMSE compared to $M1$ (Table 1). However, in periods when $G$ is around zero, this methodology can cause unrealistic $\lambda$ values. The $M3$, estimating $\lambda$ value for a reference hour (13:00—local time), had a decrease of RMSE around 15%. However, also in this method, some unrealistic $\lambda$ values can be obtained. An example is a day when precipitation occurs, or the sky is very cloudiness, increasing hardly the solar radiation at the surface. If in this hour $G$ is close to zero or the soil temperature in both depths are similar, the same problems described for $M1$ and $M2$ can happen.

The methods $M3$ and $M4$ presented results in which both the $\lambda$ and the statistical indices were very close, suggesting that using the hour that represents the largest difference between the temperatures can be considered a good alternative for the estimation of $\lambda$. However, using $M3,$ we do not have information about the diurnal cycle of $\lambda$, as mentioned by [30]. Here we do not analyze the diurnal cycle of $\lambda$ for anyone method.

The $M5$ presented the best statistical indices for all analyzed periods, proving to be the best method for $\lambda$ estimation. In the daily period, a decrease of $3.13 {\mathrm{Wm}}^{-2}$ in RMSE between $M1$ and $M5$ is observed. The greatest errors were found in the day periods for all methods since it is during this period that the greatest temperature variations occur (Table 1). The results of $\lambda$ obtained by Equation (2) and Equation (3) by the LSM ($M4$ and $M5,$ respectively) show very close values. However, the RMSE of $M4$ is higher than that of $M5$, indicating that $M5$ improved the estimate of $G_z$. Additionally, the parameter $\epsilon$ represents almost 20% of the average annual value of $G_{exp}$ that was $12.8 {\mathrm{Wm}}^{-2}$.

Figure 4 shows the dispersion diagrams of $G_z$ using the $\lambda$ values in Table 1 with the validation data. In general, the $M1$ and $M2$ models underestimate the observations when they are negative and overestimate them when they are positive. The $M3$ and $M4$ models show an overestimation of the values for both situations. Finally, the $M5$ model shows a behavior very similar to the observed data. Furthermore, the applied least square method (LSM) is adequated for this estimation.

Figure 5 shows the daily average cycle of the $G_{exp}$ values obtained by the methods $M4$ and $M5$, with the values of $\lambda$ from Table 1 used for these methods. The $M5$ represent the daily average cycle better than $M4$ for all hours, especially at night. The greatest difference between the experimental data and $M4$ and $M5$, during the nighttime, occurred at 05:30 and was $2.4 {\mathrm{Wm}}^{-2}$ and $0.3 {\mathrm{Wm}}^{-2}$, respectively; and during the daytime occurred at 14:30 and was $3.1 {\mathrm{Wm}}^{-2}$ and $0.9 {\mathrm{Wm}}^{-2}$, respectively. The diurnal peak of $G_z$ was better represented by $M5$, which had a difference of $0.4 {\mathrm{Wm}}^{-2}$ at 12:30, while $M4$ showed a difference of $2.5 {\mathrm{Wm}}^{-2}$ at the same hour.

The best results obtained by $M5$ prove that the addition of the $\epsilon$ parameter includes corrections that interfere in the estimation of the soil heat flux, i.e., $G_z$ is not simply proportional to the soil temperature gradient. Thus, the inclusion of ε also influences the soil thermal conductivity estimation, although the model still considers ε constant with depth.

Further analysis of the influence of the soil water content on the estimation of $\lambda$ is shown in Figure 6. All experimental data were used to present the relationship between the soil thermal conductivity ($\lambda$) and the soil water content ($\theta$). In this case, we used only the LSM to solve Equations (2) and (3) for each day, i.e., $M4$ and $M5,$ obtaining $\lambda$ daily values. Then, these values were grouped by intervals of $0.01$ of the soil water content daily average. With the grouped data, the $\lambda$ mean value and the respective standard deviation of each cluster were calculated, as shown in Figure 6a. Note that the use of different equation results in a different behavior between the soil thermal conductivity and the soil water content. In general, $M4$ presented higher mean values and standard deviations than $M5$.

Using another approach, the experimental half-hour data of the days presenting similar soil water content daily averages ($0.01$ intervals) were grouped. With the grouped data, a single value of $\lambda$ was estimated for each soil water content range, using $M4$ and $M5$, as shown in Figure 6b. For low and high values of soil water content, there are few experimental data; however, the models present the similar values of $\lambda$ in the intervals of low soil water content and differ considerably for intervals of high values soil water content. For the rest of the soil water content values, both equations showed similar behaviors.

The results present in Figure 6a,b show that for high values of soil moisture, both methods generate a greater variability on the $\lambda$ dataset, concluding that other variables commanding $\lambda$ values can exist when the soil water content is high, i.e., close to the field capacity.

The $\epsilon$ values of Equation (3) using this latter methodology are shown in Figure 6c. A nonlinear variation of $\epsilon$ in relation to the soil water content is observed, with values between $-3 {\mathrm{Wm}}^{-2}$ and $-1.5 {\mathrm{Wm}}^{-2}$ and an average value of $-2.1 {\mathrm{Wm}}^{-2}$.

Figure 6d shows the relationships between soil thermal conductivity and soil water content—$\lambda \left(\theta \right)$, obtained using $M4$, $M5$ (both estimating $\lambda$ for soil water content intervals), Johansen Model [35], Lu-Ren Model [36] and Tong Model [6] (see Appendix A). For Tong Model, the empirical parameters were calibrated using $\lambda \left(\theta \right)$, for $M4$ and $M5$. The values of these parameters are presented in

Table 2. Original and Calibrated parameters of Tong Model [6], and statistical index (R² and RMSE).

	Original—Tong Model	$\mathit{M}4$	$\mathit{M}5$
$a$	$1.88$	$1.14$	$1.25$
$b$	$1.67$	$1.07$	$1.01$
$c$	$3.90$	$14.06$	$10.02$
R2	$0.65$	$0.66$	$0.88$
RMSE	$0.0953$	$0.0941$	$0.0545$

together with the original values for sandy loam soil type. The results show that the variation of $\lambda \left(\theta \right)$ depends on the method applied, consequently influencing the estimation parameters describing empirical equations for $\lambda \left(\theta \right)$. The Tong Model calibrated using $M5$ presented the best statistic index (Table 2).

4. Conclusions

In this paper, we proposed a numerical model to estimate the thermal conductivity from experimental measurements of the soil heat flux and soil temperature at different depths. The proposal is based on the Fourier Law adding a constant empirical parameter ($\epsilon$). The $\epsilon$ parameter minimizes the uncertainties contained in the data from field experiments, leaving the model able to fit the behavior observed in the experimental data.

Different numerical techniques of solving Equation (2) were used and were compared to the method used in the proposed model (Equation (3)). This study shows that isolate $\lambda$ in Equation (2) ($M1$, $M2$, $M3$) can affect the accuracy of the soil thermal conductivity estimates. The use of the LSM to estimate the thermal conductivity ($M4$ and $M5$) presented better results than the other methods, which allows us to conclude that LSM is an adequate technique for this estimation. However, the $M5$ presented the best results, thus confirming the hypothesis that the addition of an empirical parameter is necessary and influences the accuracy of the soil heat flux estimate.

Currently, many micrometeorological experimental fields are measuring atmospheric and soil variables to improve the weather and climate models. For the soil, the data are obtained similarly to the experimental design used in this work. Therefore, we believe the results obtained here will contribute to this issue since the methodology proposed improved the estimation on soil thermal properties, consequently the soil heat flux, an essential component for studies on surface energy balance closure.