2.2 Analytical Solutions of the Heat Conduction Equation

Analytical solutions of the non-linear, parabolic partial differential Eq. (2.3) are hard to find and making a numerical approach is very appealing. Nevertheless, analytical solutions may contain a significant amount of physical insight and can be used to verify more complex numerical models of soil heat flow.

In order to use Fourier theory, we assume that the heat capacity $c_{h}$ and the thermal conductivity $k_{h}$ are independent of temperature and time. Equation (2.3) can then be Fourier transformed with respect to time to give an ordinary differential equation for the $n$th harmonic, namely

$${c_{h}}(z) {a_{n}}(z) i \omega n - \frac{d}{d z} \left( {k_{h}}(z) \frac{d}{d z} {a_{n}}(z) \right)$$ (2.5)

where $\omega = \frac{2 \pi}{P}$ is the fundamental frequency and $a_n$ is the Fourier coefficient of the $n$th harmonic. $P$ is the period of the surface temperature signal and is usually chosen equal to one day. The soil temperature is given by the Fourier series

$$T(z,t) = \sum_{n=0}^{\infty} {a_{n}(z) e^{i w n t}}$$ (2.6)

The physical values are given by the real parts. The boundary conditions for the $a_{n}$'s are determined by the Fourier transform of the boundary conditions imposed on $T(z,t)$.

2.2.1 Constant Thermal Properties

Assuming depth, time and temperature independent thermal properties throughout the soil, and imposing the boundary conditions $T(0,t) = T_{0}(t)$ and $T(z \to \infty,t) = T_{\infty}$, the physically relevant solution of Eq. (2.5) immediately follows and yields

$$T(z,t) = T_{\infty} + \sum_{n=1}^{\infty} e^{- \frac{z}{d_{n... ...c_{n} \sin \left( \omega n t - \frac{z}{d_{n}} \right) \right)$$ (2.7)

where $b_n$ and $c_n$ are Fourier coefficients determined by

$$b_{n} = \frac{2}{P} \int_{0}^{P} {T_{0}(t) \cos(\omega n t) dt}$$ (2.8)

and

$$c_{n} = \frac{2}{P} \int_{0}^{P} {T_{0}(t) \sin(\omega n t) dt}$$ (2.9)

The damping depth $d_n$ is determined by

$$d_{n} = \sqrt{\frac{2 D_{h,app}}{\omega n}}$$ (2.10)

with the thermal diffusivity $D_{h,app} = \frac{k_{h,app}}{c_h}$. Thus, the diurnal temperature wave of a typical soil has a damping depth of approximately $5 cm$ to $20 cm$. For the annual temperature variations these values are roughly $\sqrt{365} \approx 19$ times larger. Rapid temperature variations at the soil surface, for instance due to temporary cloud cover, have a smaller damping depth.

2.2.2 Depth Dependent Thermal Properties

In natural environments, thermal properties are rarely homogeneous throughout the soil. Soil type, soil composition and moisture content may vary considerably with depth and location. Several analytical (Shao et al., 1998; Massman, 1993; Wiltshire, 1982; Wiltshire, 1983) and approximate solutions (Nerpin and Chudnovskii, 1984; Thakur and Momoh, 1983) of the heat conduction equation with depth dependent coefficients have been studied. All analytical solutions either depend on additional assumptions about the periodicity of the solution or assume a parametric model of the thermal properties as a function of depth. According to a study by Massman (1993), the most successful approximate solution is the one proposed by Nerpin and Chudnovskii (1984). It is based on the WKB^2.1 approximation which describes the propagation of waves in inhomogeneous media.

For the sake of brevity, only a very simplified situation will be considered here. We assume that the functional form of the thermal properties is known to be linear in a certain soil layer

$$c_h(z) = c_0 (1 + \alpha z)$$ (2.11)

$$k_h(z) = k_0 (1 + \alpha z)$$ (2.12)

where $\alpha$ is the slope of the $c_h$ and $k_h$ gradient in the soil. The local coordinate system is chosen so that $c_h=c_0$ and $k_h=k_0$ at the top of the layer ($z=0$). Using the substitution

$$\xi = \sqrt{\frac{\omega n c_0}{\alpha^2 k_0}} (1+\alpha z)$$ (2.13)

in Eq. (2.5) leads to a differential equation which is solved by Bessel functions (Abramowitz and Stegun, 1972). The solution for the $n$th complex amplitude is

$$a_n(z) = \beta_n J_0(\xi e^{\frac{3}{4} \pi i}) + \gamma_n N_0(\xi e^{\frac{1}{4} \pi i})$$ (2.14)

where $J_0$ and $N_0$ are Bessel functions of the first and second kind respectively, and $\beta_n$ and $\gamma_n$ are coefficients which are determined by the boundary condition.