Sei $T$ die Temperatur der festen Erde. Vorbereitend wird das Crank-Nicolson-Verfahren mit $n$ als gegenwärtigem Zeitschritt notiert:
\[ \begin{align} T^{(n + 1)} = T^{(n)} + \Delta t\frac{1}{2}\left(\frac{\partial T^{(n)}}{\partial t} + \frac{\partial T^{(n + 1)}}{\partial t}\right) \end{align} \]
An dieser Stelle geht es dabei nur um die vertikale Wärmeleitung, die anderen Terme werden explizit berechnet:
\[ \begin{align} T^{(n + 1)} = T^{(n)} + \Delta t\newdot{T}^{(n)}_{\text{expl}} + \Delta t\frac{\lambda}{2\rho c}\frac{\partial}{\partial z}\left(\frac{\partial T^{(n)}}{\partial z} + \frac{\partial T^{(n + 1)}}{\partial z}\right) \end{align} \]
Hierbei sind $\lambda$ die Wärmeleitfähigkeit, $\rho$ die Massendichte und $c$ die spezifische Wärmekapazität. Diese Größen werden als bekannt und konstant vorausgesetzt. Weiterhin wird das vertikale Aufweiten der Gitterboxen aufgrund der geringen Dicke vernachlässigt, man macht also von der Shallow-Atmosphere-Approximation nach Glg. (33.25) Gebrauch. Den Anteil
\[ \begin{align} \frac{\lambda}{2\rho c}\frac{\partial}{\partial z}\frac{\partial T^{(n)}}{\partial z} \end{align} \]
absorbiert man in die explizite Tendenz, definiert also
\[ \begin{align} \newdot{T}^{(n)}_{\text{expl}} \to \newdot{T}^{(n)}_{\text{expl}} - \frac{1}{2\rho c}\frac{\partial}{\partial z}\frac{\lambda\partial T^{(n)}}{\partial z} \end{align} \]
um. Die vertikale Diskretisierung erfordert einen Schichtindex $1 \leq i \leq N_S$ mit $N_S$ als Anzahl der Soil-Schichten:
\[ \begin{align} T_i^{(n + 1)} = T_i^{(n)} + \Delta t\newdot{T_i}^{(n)}_{\text{expl}} + \frac{\Delta t}{2\rho_ic_i\Delta z_i}\left(\lambda'_{i-1}\frac{T^{(n + 1)}_{i-1} - T^{(n + 1)}_i}{\Delta z'_{i-1}} - \lambda'_i\frac{T^{(n + 1)}_i - T^{(n + 1)}_{i+1}}{\Delta z'_i}\right)\tag{38.5}\label{eq:heat_conduc_deriv_1} \end{align} \]
Die gestrichenen Größen beziehen sich auf die Schichtgrenzen, also
\[ \begin{align} \lambda'_i \coloneqq \frac{\lambda'_i + \lambda'_{i+1}}{2} \end{align} \]
und
\[ \begin{align} \Delta z'_i \coloneqq z_i - z_{i+1}. \end{align} \]
Bei $i = 1$ gibt es keine darüberliegende Soil-Schicht, die Wechselwrikung mit der Luft steckt in der expliziten Tendenz:
\[ \begin{align} T_1^{(n + 1)} = T_1^{(n)} + \Delta t\newdot{T_1}^{(n)}_{\text{expl}} + \frac{\Delta t}{2\rho_1c_1\Delta z_1}\lambda'_1\frac{T^{(n + 1)}_1 - T^{(n + 1)}_2}{\Delta z'_{1}} \end{align} \]
Bei $i = N_S$ ist $T_{N_S + 1}$ die als bekannt und konstant angenommene untere Ranbedingung:
\[ \begin{align} T_{N_S}^{(n + 1)} = T_{N_S}^{(n)} + \Delta t\newdot{T_{N_S}}^{(n)}_{\text{expl}} + \frac{\Delta t}{2\rho_{N_S}c_{N_S}\Delta z_{N_S}}\left(\lambda'_{N_S-1}\frac{T^{(n + 1)}_{N_S-1} - T^{(n + 1)}_{N_S}}{\Delta z'_{N_S-1}} - \lambda'_{N_S}\frac{T^{(n + 1)}_{N_S} - T_{N_S + 1}}{\Delta z'_{N_S}}\right) \end{align} \]
Der Term mit $T_{N_S + 1}$ kann in die explizite Tendenz abosrbiert werden, was auf
\[ \begin{align} T_{N_S}^{(n + 1)} = T_{N_S}^{(n)} + \Delta t\newdot{T_{N_S}}^{(n)}_{\text{expl}} + \frac{\Delta t}{2\rho_{N_S}c_{N_S}\Delta z_{N_S}}\left(\lambda'_{N_S-1}\frac{T^{(n + 1)}_{N_S-1} - T^{(n + 1)}_{N_S}}{\Delta z'_{N_S-1}} - \lambda'_{N_S}\frac{T^{(n + 1)}_{N_S}}{\Delta z'_{N_S}}\right) \end{align} \]
führt.
Man definiert den Vektor $\mathbf{x}$ der Unbekannten durch
\[ \begin{align} \mathbf{x} = \left(\begin{array}{c} T^{(n + 1)}_{1}\\ T^{(n + 1)}_{2}\\ \vdots\\ T^{(n + 1)}_{N_S} \end{array}\right), \end{align} \]
für diesen gilt ein lineares Gleichungssystem
\[ \begin{align} A\cdot\mathbf{x} = \mathbf{r} \end{align} \]
mit einer Matrix $A$ und einer rechten Seite $\mathbf{r}$. Für diese rechte Seite gilt
\[ \begin{align} \mathbf{r} = \left(\begin{array}{c} T^{(n)}_{1} + \Delta t\newdot{T}^{(n)}_{1, \text{expl}}\\ T^{(n)}_{2} + \Delta t\newdot{T}^{(n)}_{2, \text{expl}}\\ T^{(n)}_{3} + \Delta t\newdot{T}^{(n)}_{3, \text{expl}}\\ \vdots\\ T^{(n)}_{N_S - 1} + \Delta t\newdot{T}^{(n)}_{N_L - 1, \text{expl}}\\ T^{(n)}_{N_S} + \Delta t\newdot{T}^{(n)}_{N_L, \text{expl}} \end{array}\right). \end{align} \]
Für die Matrix $A$ erhält man
\[ \begin{align} A &= \left(\begin{array}{cccc} d_1 & e_1 & \dots & 0 \\ c_1 & d_2 & e_2 \hspace{2 cm}\dots & 0 \\ \vdots & \hspace{2 cm}\ddots & \ddots & \vdots \\ 0 & \dots & c_{N_L - 1} & d_{N_L} \end{array}\right) \end{align} \]
mit Vektoren $\mathbf{c}, \mathbf{e} \in \mathbb{R}^{N_S - 1}$, $\mathbf{d} \in \mathbb{R}^{N_S}$. Für diese erhält man
\[ \begin{align} c_i &= -\frac{\Delta t}{2\rho_{i+1}c_{i+1}\Delta z_{i+1}}\frac{\lambda'_i}{\Delta z'_i},\\ d_1 &= 1 + \frac{\Delta t}{2\rho_1c_1\Delta z_1}\frac{\lambda'_1}{\Delta z'_1},\\ d_i &= 1 + \frac{\Delta t}{2\rho_ic_i\Delta z_i}\left(\frac{\lambda'_{i-1}}{\Delta z'_{i-1}} + \frac{\lambda'_i}{\Delta z'_i},\right),\\ e_i &= -\frac{\Delta t}{2\rho_ic_i\Delta z_i}\frac{\lambda'_i}{\Delta z'_i}. \end{align} \]