Numerical simulation and experimental investigation of the thermal diffusivity of frozen soil under different moisture contents conditions

Soil thermal diffusivity is considered as the most important thermal characteristic of the
soil which indicates the gradient of its warming due to a unit change in its temperature. Several methods are available to determine soil thermal diffusivity from observed temperature variations. Most of these methods are based on solutions of the one-dimensional conduction heat equation with constant diffusivity and thus apply to uniform soils only.
In the absence of local heat sources or sinks, the equation describing conductive heat transfer in a one-dimensional isotropic medium is:
                                                      (1)
where  is temperature,  is time,  is the soil depth, and () is the thermal diffusivity of the soil , equal to  and  being the volumetric heat capacity . Several methods have been developed for estimating the soil diffusivity using equation (1). Horton et al. (1983) have tested six methods and concluded that Harmonic Equation and Numerical Method provide the most accurate results among all. The finite difference is considered as the most applicable method for numerical solution of the heat conduction equation in soils. For approximation of partial derivatives using finite differences, different algorithms may be used. In the present research, the Crank-Nicolson method which has a high degree of accuracy was employed. Using the above method, equation (1) can be discretised as:
                                (2)
with            
and              .
where  and  indicate the depth node and time step, respectively.
In the present work, for numerical solution of the equation, time intervals of 1 second (Δt = 1s) and spatial intervals of 1cm (Δz = 1cm) have been employed. With application of the Crank-Nicolson method, a set of simultaneous equations will be produced for each time interval. This set of simultaneous equations can be solved using different methods. In the present research the Tri-Diagonal Matrix Algorithm (TDMA) method has been employed. When the initial and boundary conditions are known and soil temperatures at different points have been measured, the soil thermal diffusivity can be determined using a trial and error technique. The approach is based on solving equation (2) iteratively by changing , and determining the -value based on which the calculated values of temperature best match observations. In the present work the criterion used for choosing  is minimizing the Root Mean Square Error (RMSE) of the calculated (C_i) against observed (Mi) temperature:
 
The above procedure was conducted to evaluate the thermal diffusivity (α) of a silty soil with mass moisture contents of 5, 10, 15 and 20 per cent. A chamber with dimensions of 500×500×800 mm was made and its walls and bottom were carefully insulated using layers of plasto-foam sheets with a thickness of 100mm, to minimize the exchange of heat with the surrounding environment. The chamber was filled with the soil of given texture and moisture. To eliminate evaporation from the top surface of the soil, it was covered by a plastic sheet. Temperature was measured using seven thermometers installed at the depths of 50, 110, 170, 250, 350, and 500 mm, as well as at the top surface of the soils. The sensors were connected to a computer, where soil temperatures were recorded at 1 min intervals. A frost condition in soil was simulated with the use of a cooling system located at the top of the soil that was able to produce temperatures as low as -20 ^oC.
Thermal diffusivity was estimated for two different thermal conditions in the soil profile: one with temperatures lower than -2 ^oC throughout the soil at all times, and the other having temperatures lower than zero degrees centigrade at some and higher than zero at other depths.
According to the results, the model used in this study led to a low RMSE (between 0.41 to 0.71 ^oC) and reasonable predictions of soil temperature for the first case (i.e. temperatures lower than -2 ^oC at all times and depths). The results showed that the values of α increased with increasing moisture content up to a critical point and then decreased. The maximum value of α occurred at 15 percent moisture content.
The model failed to estimate the soil temperature profile within an acceptable range of error in the second case. The range of RMSE values between the simulated and measured temperature in this case was found between 1.58 to 2.76 ^oC. This failure was attributed to the fact that the assumptions made in solving the heat conduction equation, namely homogeneity of soil and lack of sources and sinks of heat within the soil, were not fulfilled for the second case.