Let $T$ be the temperature of the solid earth. In preparation, the Crank-Nicolson method is noted with $n$ as the current time step:
\[ \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} \]
At this point we are only concerned with the vertical heat conduction; the other terms are calculated explicitly:
\[ \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} \]
Here $\lambda$ is the thermal conductivity, $\rho$ is the mass density and $c$ is the specific heat capacity. These quantities are assumed to be known and constant. Furthermore, the vertical expansion of the grid boxes is neglected due to the small thickness, so one makes use of the shallow atmosphere approximation according to Eq. (33.25) use. The share
\[ \begin{align} \frac{\lambda}{2\rho c}\frac{\partial}{\partial z}\frac{\partial T^{(n)}}{\partial z} \end{align} \]
is absorbed into the explicit tendency, i.e. defined
\[ \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} \]
around. The vertical discretization requires a layer index $1 \leq i \leq N_S$ with $N_S$ as the number of soil layers:
\[ \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} \]
The deleted sizes refer to the layer boundaries, i.e
\[ \begin{align} \lambda'_i \coloneqq \frac{\lambda'_i + \lambda'_{i+1}}{2} \end{align} \]
and
\[ \begin{align} \Delta z'_i \coloneqq z_i - z_{i+1}. \end{align} \]
At $i = 1$ there is no overlying soil layer, the interaction with the air is in the explicit tendency:
\[ \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} \]
At $i = N_S$, $T_{N_S + 1}$ is the lower boundary condition assumed to be known and constant:
\[ \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} \]
The term with $T_{N_S + 1}$ can be absorbed into the explicit tendency, which results in
\[ \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} \]
leads.
The vector $\mathbf{x}$ of the unknowns is defined by
\[ \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} \]
for this a linear system of equations applies
\[ \begin{align} A\cdot\mathbf{x} = \mathbf{r} \end{align} \]
with a matrix $A$ and a right-hand side $\mathbf{r}$. For this right side applies
\[ \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} \]
For the matrix $A$ one obtains
\[ \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} \]
with vectors $\mathbf{c}, \mathbf{e} \in \mathbb{R}^{N_S - 1}$, $\mathbf{d} \in \mathbb{R}^{N_S}$. For this you get
\[ \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} \]