Appendix L — Numerical solution heat flow equation

The discretized form of the heat flow equation as described in Chapter 9, is:

\[ C^{j+1/2}_i\left(T^{j+1}_i-T^j_i\right)=\frac{\Delta t^j}{\Delta z_i} \left[\lambda^{j+1/2}_{i-1/2}\frac{T^{j+1}_{i-1}-T^{j+1}_i}{\Delta z_u} - \lambda^{j+1/2}_{i+1/2}\frac{T^{j+1}_i-T^{j+1}_{i+1}}{\Delta z_l}\right] \tag{L.1}\]

where for notational convenience the subscript heat of thermal conductivity \(\lambda\) and soil heat capacity \(C\) is omitted. Equation L.1 can be rewritten as:

\[ \begin{align} &-\frac{\Delta t^j}{\Delta z_i \Delta z_u} \lambda^{j+1/2}_{i-1/2} T^{j+1}_{i-1}+ \left[C^{j+1/2}_i+\frac{\Delta t^j}{\Delta z_i \Delta z_u} \lambda^{j+1/2}_{i-1/2}+\frac{\Delta t^j}{\Delta z_i \Delta z_l} \lambda^{j+1/2}_{i+1/2}\right]T^{j+1}_i\\ &- \frac{\Delta t^j}{\Delta z_i \Delta z_l} \lambda^{j+1/2}_{i+1/2} T^{j+1}_{i+1} = C^{j+1/2}_i T^j_i \end{align} \tag{L.2}\]

Application of Equation L.2 to each node results in a tri-diagonal matrix:

\[ \left[ \begin{array}{ccccccc} \beta_1 & \gamma_1 & \cdot & \cdot & \cdot & \cdot & \cdot\\ \alpha_2 & \beta_2 & \gamma_2 & \cdot & \cdot & \cdot & \cdot\\ \cdot & \alpha_3 & \beta_3 & \gamma_3 & \cdot & \cdot & \cdot\\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot\\ \cdot & \cdot & \cdot & \cdot & \alpha_{n-1} & \beta_{n-1} & \gamma_{n-1}\\ \cdot & \cdot & \cdot & \cdot & \cdot & \alpha_n & \beta_n\\ \end{array} \right] \left[ \begin{array}{c} T^{j+1}_1\\ T^{j+1}_2\\ T^{j+1}_3\\ \cdot\\ T^{j+1}_{n-1}\\ T^{j+1}_n\\ \end{array} \right]= \left[ \begin{array}{c} f_1\\ f_2\\ f_3\\ \cdot\\ f_{n-1}\\ f_n\\ \end{array} \right] \tag{L.3}\]

where \(n\) is the number of nodal points. Next the coefficients \(\alpha_i\), \(\beta_i\), \(\gamma_i\), and \(f_i\) are explained for the intermediate nodes and for the top and bottom node.

Intermediate nodes

From Equation L.2 and Equation L.3 we may derive the coefficients:

\[ \alpha_i = -\frac{\Delta t^j}{\Delta z_i \Delta z_u} \lambda^{j+1/2}_{i-1/2} \tag{L.4}\]

\[ \beta_i = C^{j+1/2}_i+\frac{\Delta t^j}{\Delta z_i \Delta z_u} \lambda^{j+1/2}_{i-1/2}+\frac{\Delta t^j}{\Delta z_i \Delta z_l} \lambda^{j+1/2}_{i+1/2} \tag{L.5}\]

\[ \gamma_i = \frac{\Delta t^j}{\Delta z_i \Delta z_l} \lambda^{j+1/2}_{i+1/2} \tag{L.6}\]

\[ f_i = C^{j+1/2}_i T^j_i \tag{L.7}\]

Top node

The temperature at soil surface is set equal to the daily average air temperature, \(T_{avg}\). Therefore, in case of the top node, Equation L.1 transforms to:

\[ C^{j+1/2}_i\left(T^{j+1}_i-T^j_i\right)=\frac{\Delta t^j}{\Delta z_i} \left[\lambda^{j+1/2}_{1/2}\frac{T_{avg}-T^{j+1}_i}{\Delta z_u} - \lambda^{j+1/2}_{1/2}\frac{T^{j+1}_1-T^{j+1}_{2}}{\Delta z_l}\right] \tag{L.8}\]

which can be written as:

\[ \begin{align} &\left[C^{j+1/2}_1+\frac{\Delta t^j}{\Delta z_1 \Delta z_u} \lambda^{j+1/2}_{1/2}+\frac{\Delta t^j}{\Delta z_1 \Delta z_l} \lambda^{j+1/2}_{1+0.5}\right]T^{j+1}_1 - \frac{\Delta t^j}{\Delta z_1 \Delta z_l} \lambda^{j+1/2}_{1+0.5} T^{j+1}_2=\\ &C^{j+1/2}_1 T^j_1 + \frac{\Delta t^j}{\Delta z_1 \Delta z_u} \lambda^{j+1/2}_{1/2}T_{avg} \end{align} \tag{L.9}\]

Combination of Equation L.9 and Equation L.3 gives the following coefficients:

\[ \beta_1 = C^{j+1/2}_1+\frac{\Delta t^j}{\Delta z_1 \Delta z_u} \lambda^{j+1/2}_{1/2}+\frac{\Delta t^j}{\Delta z_1 \Delta z_l} \lambda^{j+1/2}_{1+0.5} \tag{L.10}\]

\[ \gamma_1 = \frac{\Delta t^j}{\Delta z_1 \Delta z_l} \lambda^{j+1/2}_{1+0.5} \tag{L.11}\]

\[ f_1 = C^{j+1/2}_1 T^j_1 + \frac{\Delta t^j}{\Delta z_1 \Delta z_u} \lambda^{j+1/2}_{1/2}T_{avg} \tag{L.12}\]

In case of a snow cover at the soil surface, the average air temperature can no longer be used as a top boundary condition, because of the insulating properties of snow. The corrected top boundary temperature, \(T_{top}\) (\(^\circ\)C), is calculated as follows:

\[ T_{top} = \frac{T_{soil\left(1\right)}aT_{avg}}{1+a} \tag{L.13}\]

where:

\[ a = 0.5\frac{\lambda_{snow}}{\lambda \left(1\right)}\frac{\Delta z\left(1\right)}{dz_{snow}} \tag{L.14}\]

\[ \lambda_{snow} = 2.86 \space \dot \space 10^{-6} \space \text{x} \space 864 \space \text{x} \space R^2_{snow} \tag{L.15}\]

\[ dz_{snow} = \frac{S_{snow}}{0.17} \tag{L.16}\]

where \(T_{soil}(1)\) is the current soil temperature of the first soil layer (\(^\circ\)C), \(T_{avg}\) is the average air temperature (\(^\circ\)C), \(\Delta z(1)\) is the thickness of the first soil layer (cm), \(dz_{snow}\) is the snow layer thickness (cm; see Chapter 10), \(S_{snow}\) is the amount of snow (equivalent liters of water), \(R_{snow}\) is the density of snow (kg m^-3) (Granberg et al. 1999), and \(\lambda_{snow}\) is the thermal conductivity of snow (J cm^-1 d^-1 \(^\circ\)C^-1). The expression for \(\lambda_{snow}\) was taken as one of the possibilities listed in Fukusako (1990) (see also Chapter 10). In case \(\lambda(1)\) < 10^-10, \(\lambda(1)\) is set equal to 100.

Bottom node

SWAP adopts a heat flow rate \(q_{heat,bot}\) (J cm^-2 d^-1) at the bottom of the soil profile. At the bottom node, the general heat flow equation, Equation L.1, transforms to:

\[ C^{j+1/2}_n\left(T^{j+1}_n-T^j_n\right)=\frac{\Delta t^j}{\Delta z_n} \left[\lambda^{j+1/2}_{n-1/2}\frac{T^{j+1}_{n-1}-T^{j+1}_n}{\Delta z_u} - q_{heat,bot} \right] \tag{L.17}\]

which can be written as:

\[ -\frac{\Delta t^j}{\Delta z_n \Delta z_u} \lambda^{j+1/2}_{n-1/2} T^{j+1}_{n-1}+ \left[C^{j+1/2}_n+\frac{\Delta t^j}{\Delta z_n \Delta z_u} \lambda^{j+1/2}_{n-1/2}\right]T^{j+1}_n = C^{j+1/2}_n T^j_n + \frac{\Delta t^j}{\Delta z_n} q_{heat,bot} \tag{L.18}\]

Combination of Equation L.18 and Equation L.3 gives the following coefficients:

\[ \alpha_n = -\frac{\Delta t^j}{\Delta z_n \Delta z_u} \lambda^{j+1/2}_{n-1/2} \tag{L.19}\]

\[ \beta_n = C^{j+1/2}_n+\frac{\Delta t^j}{\Delta z_n \Delta z_u} \lambda^{j+1/2}_{n-1/2} \tag{L.20}\]

\[ f_n = C^{j+1/2}_n T^j_n - \frac{\Delta t^j}{\Delta z_n} q_{heat,bot} \tag{L.21}\]

In case of prescribed temperature \(T_{bot}\) at the soil profile bottom, Equation L.1, transforms to:

\[ C^{j+1/2}_n\left(T^{j+1}_n-T^j_n\right)=\frac{\Delta t^j}{\Delta z_n} \left[\lambda^{j+1/2}_{n-1/2}\frac{T^{j+1}_{n-1}-T^{j+1}_n}{\Delta z_u} - \lambda^{j+1/2}_{n+1/2}\frac{T^{j+1}_n-T_{bot}}{\Delta z_l} \right] \tag{L.22}\]

which can be written as:

\[ \begin{align} &-\frac{\Delta t^j}{\Delta z_n \Delta z_u} \lambda^{j+1/2}_{n-1/2} T^{j+1}_{n-1}+ \left[C^{j+1/2}_n+\frac{\Delta t^j}{\Delta z_n \Delta z_u} \lambda^{j+1/2}_{n-1/2}+\frac{\Delta t^j}{\Delta z_n \Delta z_l} \lambda^{j+1/2}_{n+1/2}\right]T^{j+1}_n = \\ &C^{j+1/2}_n T^j_n + \frac{\Delta t^j}{\Delta z_n \Delta z_l} \lambda^{j+1/2}_{n+1/2}T_{bot} \end{align} \tag{L.23}\]

Combination of Equation L.23 and Equation L.3 gives the following coefficients:

\[ \alpha_n = -\frac{\Delta t^j}{\Delta z_n \Delta z_u} \lambda^{j+1/2}_{n-1/2} \tag{L.24}\]

\[ \beta_n = C^{j+1/2}_n + \frac{\Delta t^j}{\Delta z_n \Delta z_u} \lambda^{j+1/2}_{n-1/2} + \frac{\Delta t^j}{\Delta z_n \Delta z_l} \lambda^{j+1/2}_{n+1/2} \tag{L.25}\]

\[ f_n = C^{j+1/2}_n T^j_n + \frac{\Delta t^j}{\Delta z_n \Delta z_l} \lambda^{j+1/2}_{n+1/2}T_{bot} \tag{L.26}\]

Fukusako, S. 1990. “Thermophysical Properties of Ice, Snow, and Sea Ice.” Journal Article. International Journal of Thermophysics 11: 353–72. https://doi.org/10.1007/BF01133567.

Granberg, G., H. Grip, M. O. Löfvenius, I. Sundh, B. H. Svensson, and M. Nilsson. 1999. “A Simple Model for Simulation of Water Content, Soil Frost, and Soil Temperatures in Boreal Mixed Mires.” Journal Article. Water Resources Research 35 (12): 3771–82. https://doi.org/10.1029/1999wr900216.