Let $T$ be the temperature of the solid earth. In preparation, the Crank-Nicolson method is written down 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, one is 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 owing to their small thickness, so one makes use of the shallow-atmosphere approximation according to Eq. (33.25). The part
\[ \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. one redefines
\[ \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} \]
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 primed quantities 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 contained 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} \]
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 which a linear system of equations holds
\[ \begin{align} A\cdot\mathbf{x} = \mathbf{r} \end{align} \]
with a matrix $A$ and a right-hand side $\mathbf{r}$. For this right-hand side, one has
\[ \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_S - 1, \text{expl}}\\ T^{(n)}_{N_S} + \Delta t\newdot{T}^{(n)}_{N_S, \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_S - 1} & d_{N_S} \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 these, one obtains
\[ \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} \]
The water content is the mean mass density of water $\rho_s$ in the soil. Unfortunately, this quantity is somewhat more involved than one might initially think, since the soil does not simply conduct the water (analogously to heat conduction or diffusion in general), but rather binds the water through thermodynamic effects. Furthermore, plants also play a role here, as they are able to draw water even from relatively dry soils. In order to describe these processes, a few additional quantities are required, the simplest of which is the porosity $\Theta_s$. This is simply the maximum water content divided by the density of water $\rho_\text{H2O}$. Thus it holds that
\[ \begin{align} 0\leq \rho_s \leq \Theta_s\cdot\rho_\text{H2O}. \end{align} \]
To describe the thermodynamics of the soil water content, thermodynamic potentials $\psi_i$ are used, see Sect. 5.2.1. These are written down separately for plants and soil:
\[ \begin{align} \psi^{(p)} &= \psi^{(p)}_p + \psi^{(p)}_o + \psi^{(p)}_g,\\ \psi^{(s)} &= \psi^{(s)}_p + \psi^{(s)}_o + \psi^{(s)}_g, \end{align} \]
where $p$ denotes plants and $s$ denotes soil. Here,
$\psi_p$ is the hydraulic potential described in Eq. (5.251), but without the gravity term, which has been separated off here,
$\psi_o$ is the osmotic potential described in Sect. 5.4.2, and
$\psi_g$ is the gravitational potential.
The complete hydraulic potential can be recovered from these via $\psi_h = \psi_p + \psi_g$. For plants, $\psi_\text{me}\coloneqq\psi_p+\psi_o$ is also called the mesophyll potential. For soil, $\psi_\text{ma}\coloneqq\psi_p+\psi_o$ is also called the matrix potential.
For the matrix potential of the soil, the following empirical ansatz has been found to be adequate for many different types of soil:
\[ \begin{align} \psi_\text{ma} = \psi_s\left(\frac{\Theta_s\rho_\text{H2O}}{\rho_s}\right)^b.\tag{38.22}\label{eq:ansatz_matrixpotential} \end{align} \]
Here, $\Theta_s$ is the porosity defined above, $\psi_s<0$ is the matrix potential at saturation ($\rho_s = \Theta_s\rho_\text{H2O}$), and $b$ is the so-called $b$-parameter.
To relate potential differences to the mass flux density, one needs the hydraulic conductivity $\kappa$, for which
\[ \begin{align} \kappa = \kappa_s\left(\frac{\rho_s}{\Theta_s\rho_\text{H2O}}\right)^{2b + 3},\tag{38.23}\label{eq:kappa_from_kappa_s} \end{align} \]
holds, where $\kappa_s$ is the hydraulic conductivity at saturation. Inserting Eq. (38.22) here yields
\[ \begin{align} \kappa = \kappa_s\left(\frac{\psi_s}{\psi_\text{ma}}\right)^{2 + 3/b}. \end{align} \]
From this, together with Darcy's law, Eq. (5.253), one obtains an evolution equation for the water content in the soil:
\[ \begin{align} \frac{\partial\rho_s}{\partial t} = \frac{\partial}{\partial z}\left(\kappa\frac{\partial\psi}{\partial z}\right). \end{align} \]
For $\kappa$, one now writes
\[ \begin{align} \kappa =: \frac{\rho_\text{H2O}}{g}\kappa'\Leftrightarrow \kappa'\coloneqq\frac{g}{\rho_\text{H2O}}\kappa. \end{align} \]
$\kappa'$ has the SI unit $\frac{\text{m}}{\text{s}}$. From this it follows that
\[ \begin{align} \frac{\partial\rho_s}{\partial t} = \frac{\partial}{\partial z}\left(\frac{\rho_\text{H2O}\kappa'}{g}\frac{\partial\psi}{\partial z}\right). \end{align} \]
From now on, the prime is omitted from the notation, so that
\[ \begin{align} \frac{\partial\rho_s}{\partial t} = \frac{\partial}{\partial z}\left(\frac{\rho_\text{H2O}\kappa}{g}\frac{\partial\psi}{\partial z}\right). \end{align} \]
Using the chain rule, this can be recast as
\[ \begin{align} \frac{\partial\rho_s}{\partial t} = \frac{\partial}{\partial z}\left(\frac{\rho_\text{H2O}\kappa}{g}\frac{\partial\psi}{\partial\rho_s}\frac{\partial\rho_s}{\partial z}\right) =: \frac{\partial}{\partial z}\left(D\frac{\partial\rho_s}{\partial z}\right), \end{align} \]
where
\[ \begin{align} D \coloneqq\frac{\rho_\text{H2O}\kappa}{g}\frac{\partial\psi}{\partial\rho_s} \end{align} \]
is the diffusivity. For each of the potentials listed above, a separate flux density, and hence a separate diffusivity, can be derived. For the contribution of the matrix potential, Eq. (38.22) allows the diffusivity to be written as
\[ \begin{align} D = -\frac{\kappa\rho_\text{H2O}}{g}\psi_sb\left(\frac{\Theta_s\rho_\text{H2O}}{\rho_s}\right)^{b-1}\frac{\Theta_s\rho_\text{H2O}}{\rho_s^2} = -\frac{\kappa\rho_\text{H2O}\psi_sb}{g\rho_s}\left(\frac{\Theta_s\rho_\text{H2O}}{\rho_s}\right)^b \end{align} \]
With Eq. (38.23), it further follows that
\[ \begin{align} D = -\frac{\kappa_s\rho_\text{H2O}\psi_sb}{g\rho_s}\left(\frac{\rho_s}{\Theta_s\rho_\text{H2O}}\right)^{b + 3} = -\frac{\kappa_s\psi_sb}{g\Theta_s}\left(\frac{\rho_s}{\Theta_s\rho_\text{H2O}}\right)^{b + 2}. \end{align} \]
For the contribution of the gravitational potential, again using Darcy's law, Eq. (5.253), one obtains
\[ \begin{align} \frac{\partial\rho_s}{\partial t} = -\frac{\partial}{\partial z}\left(-\frac{\rho_\text{H2O}\kappa}{g}\frac{\partial\psi_g}{\partial z}\right) = \frac{\partial}{\partial z}\left(\frac{\rho_\text{H2O}\kappa}{g}\frac{\partial\left(gz\right)}{\partial z}\right) = \rho_\text{H2O}\frac{\partial}{\partial z}\left(\kappa\frac{\partial z}{\partial z}\right) = \rho_\text{H2O}\frac{\partial\kappa}{\partial z}. \end{align} \]
Inserting Eq. (38.23) here yields
\[ \begin{align} \frac{\partial\rho_s}{\partial t} = \rho_\text{H2O}\frac{\partial}{\partial z}\left[\kappa_s\left(\frac{\rho_s}{\Theta_s\rho_\text{H2O}}\right)^{2b + 3}\right] = \frac{\kappa_s}{\Theta_s}\left(2b + 3\right)\left(\frac{\rho_s}{\Theta_s\rho_\text{H2O}}\right)^{2b + 2}\frac{\partial\rho_s}{\partial z}. \end{align} \]
In summary,
\[ \begin{align} \frac{\partial\rho_s}{\partial t} = \frac{\partial}{\partial z}\left(D\frac{\partial\rho_s}{\partial z}\right) + \rho_\text{H2O}\frac{\partial\kappa}{\partial z}. \end{align} \]
However, some terms are still missing here. First of all, the effects of plants have not yet been taken into account. They are represented by an additional source term $Q_E\leq 0$, where $E$ stands for evapotranspiration. Furthermore, precipitation and runoff must be accounted for. Runoff is the water that drains into streams or rivers before it can seep into the ground. Since these are surface terms, however, they are not written out here. Thus, within the soil,
\[ \begin{align} \frac{\partial\rho_s}{\partial t} = \frac{\partial}{\partial z}\left(D\frac{\partial\rho_s}{\partial z}\right) + \rho_\text{H2O}\frac{\partial\kappa}{\partial z} + Q_E. \end{align} \]
This equation is now discretized in space, first for one of the inner layers $2\leq i\leq N_S - 1$:
\[ \begin{align} \frac{\partial\rho_i}{\partial t} = \frac{1}{\newtilde{z}_i - \newtilde{z}_{i+1}}\left(\newtilde{D}_i\frac{\rho_{i-1} - \rho_i}{z_{i-1} - z_i} - \newtilde{D}_{i+1}\frac{\rho_i - \rho_{i+1}}{z_i - z_{i+1}}\right) + \rho_\text{H2O}\frac{\newtilde{\kappa}_i - \newtilde{\kappa}_{i+1}}{\newtilde{z}_i - \newtilde{z}_{i+1}} + Q_{E,i}.\tag{38.37}\label{eq:soil_moisture_discrete_1} \end{align} \]
Here, the index $s$ has been replaced by the layer index $i$, and $\newtilde{z}_i\coloneqq\frac{z_{i-1} + z_i}{2}$ is the vertical position of the interface between layers $i-1$ and $i$; analogously, $\newtilde{\kappa}_i\coloneqq\frac{\kappa_{i-1} + \kappa_i}{2}$ and $\newtilde{D}_i\coloneqq\frac{D_{i-1} + D_i}{2}$ are the hydraulic conductivity and the diffusivity averaged onto the layer interface $i$, respectively. In the first layer,
\[ \begin{align} \frac{\partial\rho_1}{\partial t} = \frac{1}{\newtilde{z}_1 - \newtilde{z}_{2}}\left(P_l + L + S - R - \newtilde{D}_2\frac{\rho_1 - \rho_{2}}{z_1 - z_{2}}\right) - \rho_\text{H2O}\frac{\newtilde{\kappa}_{2}}{\newtilde{z}_1 - \newtilde{z}_{2}} + Q_{E,1},\tag{38.38}\label{eq:soil_moisture_discrete_2} \end{align} \]
where $P_l$ is the mass flux density of the liquid precipitation that has not been intercepted by plants, $L$ that of the water dripping down from plants, $S$ that of the meltwater, and $R$ that of the runoff. In the layer $i = N_S$,
\[ \begin{align} \frac{\partial\rho_{N_S}}{\partial t} = \frac{1}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}}\newtilde{D}_{N_S}\frac{\rho_{N_S-1} - \rho_{N_S}}{z_{N_S-1} - z_{N_S}} + \rho_\text{H2O}\frac{\newtilde{\kappa}_{N_S} - \newtilde{\kappa}_{N_S+1}}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}} + Q_{E,N_S},\tag{38.39}\label{eq:soil_moisture_discrete_3} \end{align} \]
where the diffusive flux across the lower boundary of the layer has been neglected. This corresponds to the assumption that the water content is vertically homogeneous below the lowest layer.
For the source density resulting from evapotranspiration, one makes the ansatz $Q_{E,i} = w_iQ_E$, where $w_i$ is a weight and $Q_E$ is the total evapotranspiration. The weights sum to unity,
\[ \begin{align} \sum_{i=1}^{N_S}w_i = 1. \end{align} \]
Usually, one sets $w_i = 0$ for layers below a certain depth.
Eqs. (38.37) - (38.39) are now additionally discretized in time:
\[ \begin{align} \frac{\rho_i^{(n+1)} - \rho_i^{(n)}}{\Delta t} &= \frac{1}{\newtilde{z}_i - \newtilde{z}_{i+1}}\left(\frac{\newtilde{D}_{i}^{(n)}}{2}\frac{\rho_{i-1}^{(n)} - \rho_i^{(n)}}{z_{i-1} - z_i} - \frac{\newtilde{D}_{i+1}^{(n)}}{2}\frac{\rho_i^{(n)} - \rho_{i+1}^{(n)}}{z_i - z_{i+1}} + \frac{\newtilde{D}_{i}^{(n)}}{2}\frac{\rho_{i-1}^{(n+1)} - \rho_i^{(n+1)}}{z_{i-1} - z_i} - \frac{\newtilde{D}_{i+1}^{(n)}}{2}\frac{\rho_i^{(n+1)} - \rho_{i+1}^{(n+1)}}{z_i - z_{i+1}}\right)\nonumber\\ &+ \frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n)}_i - \newtilde{\kappa}^{(n)}_{i+1}}{\newtilde{z}_i - \newtilde{z}_{i+1}} + \frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n+1)}_i - \newtilde{\kappa}^{(n+1)}_{i+1}}{\newtilde{z}_i - \newtilde{z}_{i+1}} + Q_{E,i}\\ \frac{\rho_1^{(n+1)} - \rho_1^{(n)}}{\Delta t} &= \frac{1}{\newtilde{z}_1 - \newtilde{z}_{2}}\left(P_l + L + S - R - \frac{\newtilde{D}_{2}^{(n)}}{2}\frac{\rho_1^{(n)} - \rho_{2}^{(n)}}{z_1 - z_{2}} - \frac{\newtilde{D}_{2}^{(n)}}{2}\frac{\rho_1^{(n+1)} - \rho_{2}^{(n+1)}}{z_1 - z_{2}}\right) - \frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n)}_{2}}{\newtilde{z}_1 - \newtilde{z}_{2}}\nonumber\\ &- \frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n+1)}_{2}}{\newtilde{z}_1 - \newtilde{z}_{2}} + Q_{E,1}\\ \frac{\rho_{N_S}^{(n+1)} - \rho_{N_S}^{(n)}}{\Delta t} &= \frac{1}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}}\left(\frac{\newtilde{D}_{N_S}^{(n)}}{2}\frac{\rho_{N_S-1}^{(n)} - \rho_{N_S}^{(n)}}{z_{N_S-1} - z_{N_S}} + \frac{\newtilde{D}_{N_S}^{(n)}}{2}\frac{\rho_{N_S-1}^{(n+1)} - \rho_{N_S}^{(n+1)}}{z_{N_S-1} - z_{N_S}}\right) + \frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n)}_{N_S} - \newtilde{\kappa}^{(n)}_{N_S+1}}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}}\nonumber\\ &+ \frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n+1)}_{N_S} - \newtilde{\kappa}^{(n+1)}_{N_S+1}}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}} + Q_{E,N_S} \end{align} \]
Here, $n$ is the index of the time step. This can be rearranged to
\[ \begin{align} \rho_i^{(n+1)} &= \rho_i^{(n)} + \frac{\Delta t}{\newtilde{z}_i - \newtilde{z}_{i+1}}\left(\frac{\newtilde{D}_{i}^{(n)}}{2}\frac{\rho_{i-1}^{(n)} - \rho_i^{(n)}}{z_{i-1} - z_i} - \frac{\newtilde{D}_{i+1}^{(n)}}{2}\frac{\rho_i^{(n)} - \rho_{i+1}^{(n)}}{z_i - z_{i+1}} + \frac{\newtilde{D}_{i}^{(n)}}{2}\frac{\rho_{i-1}^{(n+1)} - \rho_i^{(n+1)}}{z_{i-1} - z_i} - \frac{\newtilde{D}_{i+1}^{(n)}}{2}\frac{\rho_i^{(n+1)} - \rho_{i+1}^{(n+1)}}{z_i - z_{i+1}}\right)\nonumber\\ &+ \Delta t\frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n)}_i - \newtilde{\kappa}^{(n)}_{i+1}}{\newtilde{z}_i - \newtilde{z}_{i+1}} + \Delta t\frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n+1)}_i - \newtilde{\kappa}^{(n+1)}_{i+1}}{\newtilde{z}_i - \newtilde{z}_{i+1}} + \Delta tQ_{E,i},\\ \rho_1^{(n+1)} &= \rho_1^{(n)} + \frac{\Delta t}{\newtilde{z}_1 - \newtilde{z}_{2}}\left(P_l + L + S - R - \frac{\newtilde{D}_{2}^{(n)}}{2}\frac{\rho_1^{(n)} - \rho_{2}^{(n)}}{z_1 - z_{2}} - \frac{\newtilde{D}_{2}^{(n)}}{2}\frac{\rho_1^{(n+1)} - \rho_{2}^{(n+1)}}{z_1 - z_{2}}\right)\nonumber\\ &- \Delta t\frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n)}_{2}}{\newtilde{z}_1 - \newtilde{z}_{2}} - \Delta t\frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n+1)}_{2}}{\newtilde{z}_1 - \newtilde{z}_{2}} + \Delta tQ_{E,1},\\ \rho_{N_S}^{(n+1)} &= \rho_{N_S}^{(n)} + \frac{\Delta t}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}}\left(\frac{\newtilde{D}_{N_S}^{(n)}}{2}\frac{\rho_{N_S-1}^{(n)} - \rho_{N_S}^{(n)}}{z_{N_S-1} - z_{N_S}} + \frac{\newtilde{D}_{N_S}^{(n)}}{2}\frac{\rho_{N_S-1}^{(n+1)} - \rho_{N_S}^{(n+1)}}{z_{N_S-1} - z_{N_S}}\right)\nonumber\\ &+ \Delta t\frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n)}_{N_S} - \newtilde{\kappa}^{(n)}_{N_S+1}}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}} + \Delta t\frac{\rho_\text{H2O}}{2}\frac{\newtilde{\kappa}^{(n+1)}_{N_S} - \newtilde{\kappa}^{(n+1)}_{N_S+1}}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}} + \Delta tQ_{E,N_S}. \end{align} \]
Denoting the explicit tendencies by
\[ \begin{align} \newdot{\rho}^{(n)}_{i, \text{expl}} &\coloneqq \frac{1}{\newtilde{z}_i - \newtilde{z}_{i+1}}\left[\frac{\newtilde{D}_{i}^{(n)}}{2}\frac{\rho_{i-1}^{(n)} - \rho_i^{(n)}}{z_{i-1} - z_i} - \frac{\newtilde{D}_{i+1}^{(n)}}{2}\frac{\rho_i^{(n)} - \rho_{i+1}^{(n)}}{z_i - z_{i+1}} + \frac{\rho_\text{H2O}}{2}\left(\newtilde{\kappa}^{(n)}_i - \newtilde{\kappa}^{(n)}_{i+1}\right)\right] + Q_{E,i},\\ \newdot{\rho}^{(n)}_{1, \text{expl}} &\coloneqq \frac{1}{\newtilde{z}_1 - \newtilde{z}_{2}}\left(P_l + L + S - R - \frac{\newtilde{D}_{2}^{(n)}}{2}\frac{\rho_1^{(n)} - \rho_{2}^{(n)}}{z_1 - z_{2}} - \frac{\rho_\text{H2O}}{2}\newtilde{\kappa}^{(n)}_{2}\right) + Q_{E,1},\\ \newdot{\rho}^{(n)}_{N_S, \text{expl}} &\coloneqq \frac{1}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}}\left[\frac{\newtilde{D}_{N_S}^{(n)}}{2}\frac{\rho_{N_S-1}^{(n)} - \rho_{N_S}^{(n)}}{z_{N_S-1} - z_{N_S}} + \frac{\rho_\text{H2O}}{2}\left(\newtilde{\kappa}^{(n)}_{N_S} - \newtilde{\kappa}^{(n)}_{N_S+1}\right)\right] + Q_{E,N_S}, \end{align} \]
one can write this as
\[ \begin{align} \rho_i^{(n+1)} &= \rho_i^{(n)} + \Delta t\newdot{\rho}^{(n)}_{i, \text{expl}} + \frac{\Delta t}{2\left(\newtilde{z}_i - \newtilde{z}_{i+1}\right)}\left[\newtilde{D}_{i}^{(n)}\frac{\rho_{i-1}^{(n+1)} - \rho_i^{(n+1)}}{z_{i-1} - z_i} - \newtilde{D}_{i+1}^{(n)}\frac{\rho_i^{(n+1)} - \rho_{i+1}^{(n+1)}}{z_i - z_{i+1}} + \rho_\text{H2O}\left(\newtilde{\kappa}^{(n+1)}_i - \newtilde{\kappa}^{(n+1)}_{i+1}\right)\right],\\ \rho_1^{(n+1)} &= \rho_1^{(n)} + \Delta t\newdot{\rho}^{(n)}_{1, \text{expl}} + \frac{\Delta t}{2\left(\newtilde{z}_1 - \newtilde{z}_{2}\right)}\left(-\newtilde{D}_{2}^{(n)}\frac{\rho_1^{(n+1)} - \rho_{2}^{(n+1)}}{z_1 - z_{2}} - \rho_\text{H2O}\newtilde{\kappa}^{(n+1)}_{2}\right),\\ \rho_{N_S}^{(n+1)} &= \rho_{N_S}^{(n)} + \Delta t\newdot{\rho}^{(n)}_{N_S, \text{expl}} + \frac{\Delta t}{2\left(\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}\right)}\left[\newtilde{D}_{N_S}^{(n)}\frac{\rho_{N_S-1}^{(n+1)} - \rho_{N_S}^{(n+1)}}{z_{N_S-1} - z_{N_S}} + \rho_\text{H2O}\left(\newtilde{\kappa}^{(n+1)}_{N_S} - \newtilde{\kappa}^{(n+1)}_{N_S+1}\right)\right] \end{align} \]
Writing out the averaging at the layer boundaries explicitly leads to
\[ \begin{align} \rho_i^{(n+1)} &= \rho_i^{(n)} + \Delta t\newdot{\rho}^{(n)}_{i, \text{expl}} + \frac{\Delta t}{2\left(\newtilde{z}_i - \newtilde{z}_{i+1}\right)}\left[\newtilde{D}_{i}^{(n)}\frac{\rho_{i-1}^{(n+1)} - \rho_i^{(n+1)}}{z_{i-1} - z_i} - \newtilde{D}_{i+1}^{(n)}\frac{\rho_i^{(n+1)} - \rho_{i+1}^{(n+1)}}{z_i - z_{i+1}} + \frac{\rho_\text{H2O}}{2}\left(\kappa^{(n+1)}_{i-1} - \kappa^{(n+1)}_{i+1}\right)\right],\\ \rho_1^{(n+1)} &= \rho_1^{(n)} + \Delta t\newdot{\rho}^{(n)}_{1, \text{expl}} + \frac{\Delta t}{2\left(\newtilde{z}_1 - \newtilde{z}_{2}\right)}\left[-\newtilde{D}_{2}^{(n)}\frac{\rho_1^{(n+1)} - \rho_{2}^{(n+1)}}{z_1 - z_{2}} - \frac{\rho_\text{H2O}}{2}\left(\kappa^{(n+1)}_{1} + \kappa^{(n+1)}_{2}\right)\right],\\ \rho_{N_S}^{(n+1)} &= \rho_{N_S}^{(n)} + \Delta t\newdot{\rho}^{(n)}_{N_S, \text{expl}} + \frac{\Delta t}{2\left(\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}\right)}\left[\newtilde{D}_{N_S}^{(n)}\frac{\rho_{N_S-1}^{(n+1)} - \rho_{N_S}^{(n+1)}}{z_{N_S-1} - z_{N_S}} + \frac{\rho_\text{H2O}}{2}\left(\kappa^{(n+1)}_{N_S-1} - \kappa^{(n+1)}_{N_S}\right)\right]. \end{align} \]
Here, $\newtilde{\kappa}_{N_S+1} = \kappa_{N_S}$ has been assumed. To express the $\kappa^{(n+1)}_{i}$ in terms of the $\rho_{i}^{(n+1)}$, one uses Eq. (38.23):
\[ \begin{align} \kappa_i^{(n+1)} &= \kappa_i^{(n)} + \kappa_s\left(2b + 3\right)\left(\frac{\rho_i^{(n)}}{\Theta_s\rho_\text{H2O}}\right)^{2b + 2}\frac{\rho^{(n+1)}_i - \rho^{(n)}_i}{\Theta_s\rho_\text{H2O}} + O\left(\Delta\rho^2_i\right)\nonumber\\ &= \kappa_i^{(n)} + \kappa_s\left(2b + 3\right)\left(\frac{\rho_i^{(n)}}{\Theta_s\rho_\text{H2O}}\right)^{2b + 2}\frac{\rho^{(n+1)}_i}{\Theta_s\rho_\text{H2O}} - \kappa_s\left(2b + 3\right)\left(\frac{\rho_i^{(n)}}{\Theta_s\rho_\text{H2O}}\right)^{2b + 2}\frac{\rho^{(n)}_i}{\Theta_s\rho_\text{H2O}} + O\left(\Delta\rho^2_i\right)\nonumber\\ &= \kappa_i^{(n)} + \kappa_s\left(2b + 3\right)\left(\frac{\rho_i^{(n)}}{\Theta_s\rho_\text{H2O}}\right)^{2b + 2}\frac{\rho^{(n+1)}_i}{\Theta_s\rho_\text{H2O}} - \kappa_s\left(2b + 3\right)\left(\frac{\rho_i^{(n)}}{\Theta_s\rho_\text{H2O}}\right)^{2b + 3} + O\left(\Delta\rho^2_i\right)\nonumber\\ &= \kappa_i^{(n)} + \kappa_s\left(2b + 3\right)\left(\frac{\rho_i^{(n)}}{\Theta_s\rho_\text{H2O}}\right)^{2b + 2}\frac{\rho^{(n+1)}_i}{\Theta_s\rho_\text{H2O}} - \left(2b + 3\right)\kappa_i^{(n)} + O\left(\Delta\rho^2_i\right) \end{align} \]
With $f_i\coloneqq\kappa_s\left(2b + 3\right)\left(\frac{\rho_i^{(n)}}{\Theta_s\rho_\text{H2O}}\right)^{2b + 2}\frac{1}{\Theta_s}$, one can write this as
\[ \begin{align} \kappa_i^{(n+1)}&= \kappa_i^{(n)} + f_i\frac{\rho^{(n+1)}_i}{\rho_\text{H2O}} - \left(2b + 3\right)\kappa_i^{(n)} + O\left(\Delta\rho^2_i\right) = -\left(2b + 2\right)\kappa_i^{(n)} + f_i\frac{\rho^{(n+1)}_i}{\rho_\text{H2O}} + O\left(\Delta\rho^2_i\right) \end{align} \]
Absorbing the terms $\propto-\left(2b + 2\right)\kappa_i^{(n)}$ into the explicit tendencies, i.e.
\[ \begin{align} \newdot{\rho}^{(n)}_{i, \text{expl}} &\coloneqq \frac{1}{\newtilde{z}_i - \newtilde{z}_{i+1}}\left[\frac{\newtilde{D}_{i}^{(n)}}{2}\frac{\rho_{i-1}^{(n)} - \rho_i^{(n)}}{z_{i-1} - z_i} - \frac{\newtilde{D}_{i+1}^{(n)}}{2}\frac{\rho_i^{(n)} - \rho_{i+1}^{(n)}}{z_i - z_{i+1}} - \left(2b + 1\right)\frac{\rho_\text{H2O}}{4}\left(\kappa^{(n)}_{i-1} - \kappa^{(n)}_{i+1}\right)\right] + Q_{E,i},\\ \newdot{\rho}^{(n)}_{1, \text{expl}} &\coloneqq \frac{1}{\newtilde{z}_1 - \newtilde{z}_{2}}\left[P_l + L + S - R - \frac{\newtilde{D}_{2}^{(n)}}{2}\frac{\rho_1^{(n)} - \rho_{2}^{(n)}}{z_1 - z_{2}} + \left(2b + 1\right)\frac{\rho_\text{H2O}}{4}\left(\kappa^{(n)}_{1} + \kappa^{(n)}_{2}\right)\right] + Q_{E,1},\\ \newdot{\rho}^{(n)}_{N_S, \text{expl}} &\coloneqq \frac{1}{\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}}\left[\frac{\newtilde{D}_{N_S}^{(n)}}{2}\frac{\rho_{N_S-1}^{(n)} - \rho_{N_S}^{(n)}}{z_{N_S-1} - z_{N_S}} - \left(2b + 1\right)\frac{\rho_\text{H2O}}{4}\left(\kappa^{(n)}_{N_S-1} - \kappa^{(n)}_{N_S}\right)\right] + Q_{E,N_S}, \end{align} \]
leads to
\[ \begin{align} \rho_i^{(n+1)} &= \rho_i^{(n)} + \Delta t\newdot{\rho}^{(n)}_{i, \text{expl}} + \frac{\Delta t}{2\left(\newtilde{z}_i - \newtilde{z}_{i+1}\right)}\left[\newtilde{D}_{i}^{(n)}\frac{\rho_{i-1}^{(n+1)} - \rho_i^{(n+1)}}{z_{i-1} - z_i} - \newtilde{D}_{i+1}^{(n)}\frac{\rho_i^{(n+1)} - \rho_{i+1}^{(n+1)}}{z_i - z_{i+1}} + \frac{1}{2}\left(f_{i-1}\rho^{(n+1)}_{i-1} - f_{i+1}\rho^{(n+1)}_{i+1}\right)\right],\\ \rho_1^{(n+1)} &= \rho_1^{(n)} + \Delta t\newdot{\rho}^{(n)}_{1, \text{expl}} + \frac{\Delta t}{2\left(\newtilde{z}_1 - \newtilde{z}_{2}\right)}\left[-\newtilde{D}_{2}^{(n)}\frac{\rho_1^{(n+1)} - \rho_{2}^{(n+1)}}{z_1 - z_{2}} - \frac{1}{2}\left(f_{1}\rho^{(n+1)}_{1} + f_{2}\rho^{(n+1)}_{2}\right)\right],\\ \rho_{N_S}^{(n+1)} &= \rho_{N_S}^{(n)} + \Delta t\newdot{\rho}^{(n)}_{N_S, \text{expl}} + \frac{\Delta t}{2\left(\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}\right)}\left[\newtilde{D}_{N_S}^{(n)}\frac{\rho_{N_S-1}^{(n+1)} - \rho_{N_S}^{(n+1)}}{z_{N_S-1} - z_{N_S}} + \frac{1}{2}\left(f_{N_S-1}\rho^{(n+1)}_{N_S-1} - f_{N_S}\rho^{(n+1)}_{N_S}\right)\right]. \end{align} \]
Analogously to Sect. 38.1, one defines the vector $\mathbf{x}$ of the unknowns by
\[ \begin{align} \mathbf{x} = \left(\begin{array}{c} \rho^{(n + 1)}_{1}\\ \rho^{(n + 1)}_{2}\\ \vdots\\ \rho^{(n + 1)}_{N_S} \end{array}\right), \end{align} \]
for which a linear system of equations holds
\[ \begin{align} A\cdot\mathbf{x} = \mathbf{r} \end{align} \]
with a matrix $A$ and a right-hand side $\mathbf{r}$. For this right-hand side, one has
\[ \begin{align} \mathbf{r} = \left(\begin{array}{c} \rho^{(n)}_{1} + \Delta t\newdot{\rho}^{(n)}_{1, \text{expl}}\\ \rho^{(n)}_{2} + \Delta t\newdot{\rho}^{(n)}_{2, \text{expl}}\\ \rho^{(n)}_{3} + \Delta t\newdot{\rho}^{(n)}_{3, \text{expl}}\\ \vdots\\ \rho^{(n)}_{N_S - 1} + \Delta t\newdot{\rho}^{(n)}_{N_S - 1, \text{expl}}\\ \rho^{(n)}_{N_S} + \Delta t\newdot{\rho}^{(n)}_{N_S, \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_S - 1} & d_{N_S} \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 these, one obtains
\[ \begin{align} c_i &= \frac{\Delta t}{2\left(\newtilde{z}_{i+1} - \newtilde{z}_{i+2}\right)}\left(-\frac{\newtilde{D}_{i+1}^{(n)}}{z_{i} - z_{i+1}} - \frac{f_{i}}{2}\right),\\ d_1 &= 1 + \frac{\Delta t}{2\left(\newtilde{z}_1 - \newtilde{z}_{2}\right)}\left(\frac{\newtilde{D}_{2}^{(n)}}{z_{1} - z_{2}} + \frac{1}{2}f_1\right),\\ d_i &= 1 + \frac{\Delta t}{2\left(\newtilde{z}_i - \newtilde{z}_{i+1}\right)}\left(\frac{\newtilde{D}_{i}^{(n)}}{z_{i-1} - z_i} + \frac{\newtilde{D}_{i+1}^{(n)}}{z_{i} - z_{i+1}}\right),\\ d_{N_S} &= 1 + \frac{\Delta t}{2\left(\newtilde{z}_{N_S} - \newtilde{z}_{N_S+1}\right)}\left(\frac{\newtilde{D}_{N_S}^{(n)}}{z_{N_S-1} - z_{N_S}} + \frac{f_{N_S}}{2}\right),\\ e_i &= \frac{\Delta t}{2\left(\newtilde{z}_i - \newtilde{z}_{i+1}\right)}\left(-\frac{\newtilde{D}_{i+1}^{(n)}}{z_{i} - z_{i+1}} + \frac{f_{i+1}}{2}\right). \end{align} \]