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_S - 1, \text{expl}}\\ T^{(n)}_{N_S} + \Delta t\newdot{T}^{(n)}_{N_S, \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_S - 1} & d_{N_S} \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} \]
Der Wassergehalt ist die gemittelte Massendichte des Wassers $\rho_s$ im Boden. Diese Größe ist leider etwas komplizierter als man zunächst denken würde, da der Boden nicht einfach das Wasser durchleitet (analog zur Wärmeleitung oder allgemein Diffusion), sondern das Wasser durch thermodynamische Effekte an sich bindet. Weiterhin spielen auch die Pflanzen hier eine Rolle, die Fähigkeiten haben, auch aus relativ trockenen Böden Wasser aufzunehmen. Um diese Prozesse beschreiben zu können, sind einige weitere Größen notwendig, die einfachste davon ist die Porosität $\Theta_s$. Dies ist einfach der maximale Wassergehalt dividiert durch die Dichte von Wasser $\rho_\text{H2O}$. Es gilt also
\[ \begin{align} 0\leq \rho_s \leq \Theta_s\cdot\rho_\text{H2O}. \end{align} \]
Für die Beschreibung der Thermodynamik des Wassergehaltes im Boden werden thermodynamische Potentiale $\psi_i$ verwendet, s. Abschn. 5.2.1. Dies notiert man für Pflanzen und Soil separat:
\[ \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} \]
wobei $p$ für Pflanzen steht und $s$ für Soil. Hierbei sind
$\psi_p$ das in Glg. (5.251) beschriebene hydraulische Potential ohne den Schwereterm, der hier abgetrennt wurde,
$\psi_o$ das in Abschn. 5.4.2 beschriebene osmotische Potential und
$\psi_g$ das Schwerepotential.
Das vollständige hydraulische Potential lässt sich hieraus durch $\psi_h = \psi_p + \psi_g$ ableiten. Für Pflanzen bezeichnet man $\psi_\text{me}\coloneqq\psi_p+\psi_o$ auch als Mesophyll-Potential. Für Soil bezeichnet man $\psi_\text{ma}\coloneqq\psi_p+\psi_o$ auch als Matrixpotential.
Für das Matrixpotential von Soil wurde empirisch folgender Ansatz als adäquat für viele verschiedene Arten von Böden erkannt:
\[ \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} \]
Hierbei sind $\Theta_s$ die oben definierte Porosität, $\psi_s<0$ das Matrixpotential bei Sättigung ($\rho_s = \Theta_s\rho_\text{H2O}$), und $b$ der sogenannte b-Parameter.
Um Potentialdifferenzen mit der Massenflussdichte zu verknüpfen, benötigt man die hydraulische Leitfähigkeit $\kappa$, für diese gilt
\[ \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} \]
wobei $\kappa_s$ die hydraulische Leitfähigkeit bei Sättigung ist. Setzt man hier Glg. (38.22) ein, erhält man
\[ \begin{align} \kappa = \kappa_s\left(\frac{\psi_s}{\psi_\text{ma}}\right)^{2 + 3/b}. \end{align} \]
Hieraus erhält man mit dem Darcy-Gesetz Glg. (5.253) eine Evolutionsgleichung für den Wassergehalt im Boden:
\[ \begin{align} \frac{\partial\rho_s}{\partial t} = \frac{\partial}{\partial z}\left(\kappa\frac{\partial\psi}{\partial z}\right). \end{align} \]
Für $\kappa$ notiert man nun
\[ \begin{align} \kappa =: \frac{\rho_\text{H2O}}{g}\kappa'\Leftrightarrow \kappa'\coloneqq\frac{g}{\rho_\text{H2O}}\kappa. \end{align} \]
$\kappa'$ hat die SI-Einheit $\frac{\text{m}}{\text{s}}$. Hieraus folgt
\[ \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} \]
Der Strich wird in der Notation von nun an vernachlässigt, also
\[ \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} \]
Mit der Kettenregel kann man dies umformen zu
\[ \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} \]
hierbei ist
\[ \begin{align} D \coloneqq\frac{\rho_\text{H2O}\kappa}{g}\frac{\partial\psi}{\partial\rho_s} \end{align} \]
die Diffusivität. Für jedes der oben aufgelisteten Potentiale kann eine eigene Flussdichte bzw. Diffusivität hergeleitet werden. Für den Anteil des Matrixpotentials kann man mit Glg. (38.22) die Diffusivität wie folgt formulieren:
\[ \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} \]
Mit Glg. (38.23) folgt weiter
\[ \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} \]
Für den Anteil des Schwerepotentials gilt, wiederum mit dem Darcy-Gesetz Glg. (5.253)
\[ \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} \]
Setzt man hier Glg. (38.23) ein, erhält man
\[ \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} \]
Zusammenfassend gilt
\[ \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} \]
Hier fehlen jedoch noch einige Terme. Zunächst einmal sind hier die Effekte von Pflanzen noch nicht berücksichtigt. Diese werden durch einen zusätzlichen Quellterm $Q_E\leq 0$ repräsentiert, wobei $E$ für Evapotranspiration steht. Weiterhin sind Niederschlag und Runoff zu berücksichtigen. Runoff ist hierbei das Wasser, das in Bäche oder Flüsse abfließt, bevor es versickern kann. Da dies Oberflächenterme sind, werden diese hier jedoch nicht mitnotiert. Somit gilt innerhalb des Bodens
\[ \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} \]
Nun diskretisiert man diese Gleichung örtlich, zunächst für eine der inneren Schichten $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} \]
Der Index $s$ wurde dabei durch den Schichtindex $i$ ersetzt, und $\newtilde{z}_i\coloneqq\frac{z_{i-1} + z_i}{2}$ ist die vertikale Position der Begrenzungsfläche zwischen Schicht $i-1$ und $i$, analog sind $\newtilde{\kappa}_i\coloneqq\frac{\kappa_{i-1} + \kappa_i}{2}$ die an die Schichtgrenze $i$ gemittelte hydraulische Leitfähigkeit und $\newtilde{D}_i\coloneqq\frac{D_{i-1} + D_i}{2}$ die an die Schichtgrenze $i$ gemittelte Diffusivität. In der ersten Schicht gilt
\[ \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} \]
hierbei ist $P_l$ die Massenflussdichte des flüssigen Niederschlags, der nicht von Pflanzen aufgefangen wurde, $L$ die des von Pflanzen heruntertropfenden Wassers, $S$ die des Schmelzwassers und $R$ die des Runoffs. In der Schicht $i = N_S$ gilt
\[ \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} \]
dabei wurde der diffusive Fluss an der Untergrenze der Schicht vernachlässigt. Dies entspricht der Annahme, dass unterhalb der untersten Schicht der Wassergehalt vertikal homogen ist.
Für die aus der Evapotranspiration folgende Quelldichte setzt man $Q_{E,i} = w_iQ_E$ an, wobei $w_i$ ein Gewicht ist und $Q_E$ die Gesamt-Evapotranspiration. Für die Summe der Gewichte gilt
\[ \begin{align} \sum_{i=1}^{N_S}w_i = 1. \end{align} \]
Üblicherweise setzt man $w_i = 0$ für Schichten unterhalb einer gewissen Tiefe.
Nun diskretisiert man die Glg.en (38.37) - (38.39) zusätzlich zeitlich:
\[ \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} \]
Hierbei ist $n$ der Index des Zeitschritts. Dies kann man umformen zu
\[ \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} \]
Notiert man für die expliziten Tendenzen
\[ \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} \]
kann man dies als
\[ \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} \]
notieren. Schreibt man die Mittelung an die Schichtgrenzen explizit, führt dies auf
\[ \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} \]
Dabei wurde $\newtilde{\kappa}_{N_S+1} = \kappa_{N_S}$ angesetzt. Um die $\kappa^{(n+1)}_{i}$ durch die $\rho_{i}^{(n+1)}$ auszudrücken, verwendet man Glg. (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} \]
Mit $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}$ kann man dies als
\[ \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} \]
formulieren. Absorbiert man die Terme $\propto-\left(2b + 2\right)\kappa_i^{(n)}$ in die expliziten Tendenzen, also
\[ \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} \]
führt dies auf
\[ \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} \]
Analog zu Abschn. 38.1 definiert man den Vektor $\mathbf{x}$ der Unbekannten durch
\[ \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} \]
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} \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} \]
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_S - 1} & d_{N_S} \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\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} \]