Firstly, we may declare the thermal conductivity of the soil layers as constant value (2nd assumption). Indeed, if heat transfer is plane and steady, and if is not large (in our case ), one can make a reasonably accurate approximation using a constant average value of (Lienhard &Lienhard 2003, p. 51).
Then, we may declare the density and the specific heat of the soil layers as constants and (3rd assumption). It is necessary to note that such assumption is somewhat groundless, especially for the cases of strictly inhomogeneous soils. However, only after this we may introduce a constant diffusion coefficient () of the soil (Bird et al. 2002, p. 268).
Finally, we may declare that the rate of internal energy conversion is negligibly small (4th assumption). In fact, such assumption is declaration of absence of the heat generation or consumption within the soil. Again we note that this assumption can be groundless for the inhomogeneous soils with stone inclusions, fluid- or air-filled interstices with internal convective flows. Moreover, the heat transfer in such porous and composite media is very difficult to analyze (Bird et al. 2002, p. 281-283).
for the steady boundary conditions and ; is the thickness of the soil upon the rock background. Solution (5) is easy to derive analytically (Haberman 1983, p. 13-14), so we will use it for checking our numerical model by approximation at .
Heat transfer model parameters
The simplified problem (4) is stated by following values: m, m2s-1.
Boundary conditions are: , (1st case), and (2nd case).
Initial values are stated by equation .
Model discretization is stated by number of soil layers , their thicknesses m, and the timestep (in seconds) which we can modify (240s, or 550s).
Model geometry is shown at the figure 1. In our model zero-level () is located at the surface of the rock background because the soil thickness is rather unstable parameter. Indeed, thickness of the real soil cover is a function and for the small areas only. Therefore, we will use more "stable" rock surface to count out -values of the soil layers.
Model dynamics (i.e. heat conduction process) is described by (4) which is transformed in a form of difference equation
here, denotes ; , ; , .
The work equation (6) is derived from (4) by FTCS scheme, when forward differentiation was used for and centered differentiation was used for (Boyce & DiPrima 2001, p. 419f).
Figure 1 - Model geometry
1. Let us transform the FTCS scheme (6) into an explicit form:
For the bottom soil layer () we have
because of .
For the layer near the soil surface () we have
because of .
2. To create a Matlab script for solving equation (4) in the explicit FTCS form (7), we can use both initial script and examples of (Mathews & Fink 1999, p. 526-536). Work model code is in