WEST — Modellierung des Wasser-, Energie und Stofftransports in Böden

Inhaltsverzeichnis

> index / notes : west

status: in Arbeit

Auftakt

Informationen für Studierende

Dies sind Notizen, die im Modul Modellierung des Wasser-, Energie und Stofftransports in Böden, welches jedes Wintersemester angeboten wird, benutzt werden. In diesem Modul, halten wir sechs Vorlesungen im zweiten Teil des Semesters. Der erste Teil dieses Moduls wird von der AG Bodenwissenschaften am Institut für Geoökologie übernommen.

Ein grober Zeitplan für unsere sechs gemeinsamen Wochen sieht wie folgt aus:

Woche 1 Numerische Methoden: FDM, FVM
Woche 2 Bodenwasserfluss: Richards(on) Gleichung
Woche 3 Stofftransport: Advektions–Dispersions-Gleichung
Woche 4 Reaktiver Stofftransport
Woche 5 Wurzelwasseraufnahme
Woche 6 Pflanzenhydraulik, System Boden–Pflanze–Atmosphäre

Einleitung

I kannt Da ned song wia lang
I hock und schau und nix versteh
Wia's Wasser geht
Weich und gleich und immer wieder nei

— Dreiviertelblut, Auf und davo

Motivation

Satellitenbilder von Trockengebieten zeigen oft interessante räumliche Muster von Vegetation: Kreise, Punkte, Bänder, Labyrinthe und sogenannte Feenkreise. Der Ursprung dieser Muster ist auf die lokale Wasserverfügbarkeit im Gebiet zurückzuführen, wobei die Pflanzen jedoch den lokalen Wasserhaushalt maßgeblich beeinflussen. Dieser Prozess wird "Selbst-organisation der Vegetation" genannt, siehe z.B., Rietkerk et al. (2002)°. Die mathematische Beschreibung dieses Prozesses der Selbst-organisation wird Turing'sche Instabilität° genannt, welche auch viele andere biologische Prozesse erklärt.

Bodenwasser reguliert neben der Ausbreitung der Vegetation auch die Vitalität der Pflanzen. So ist Baumsterben oft auf wasserstressbedingte Kavitation zurückzuführen, die dann auftritt wenn der Boden sehr trocken fällt, siehe Ruffault et al. (2022)°.

Diese enge Beziehung zwischen Verfügbarkeit von Bodenwasser und Pflanzenvitalität motiviert die modellbasierte Untersuchung des Bodenwasserhaushaltes in der Ökohydrologie. Dies beinhaltet auch die Modellierung des Bodenwasserflusses und den Transportprozessen in Böden—eine Schnittstelle zur Bodenphysik.

Julia

Das Werkzeug spielt keine Rolle, wir müssen jedoch eins wählen. Hier verwenden wir Julia°, eine Programmiersprache, die sich auf wissenschaftliches Rechnen spezialisiert. Wer Vorkenntnisse in den Programmiersprachen Matlab oder Python hat, wird Julia ohne Probleme lesen können. Eine detaillierte Auseinandersetzung mit der Sprache ist nicht erforderlich, da in diesem Kurs nicht programmiert wird.

Alle in diesem Kurs benutzten Implementierungen sind archiviert unter https://gitlab.com/iozgen/west-vorlesungsbeispiele.

Empfohlene Literatur

Es gibt mehrere gute Lehrbücher, wir empfehlen:

Zum Abschluss sei Ökten (2019)° hier als Lehrbuch für alle empfohlen, die an Julia interessiert sind und sich selbständig einarbeiten möchten.

Numerik partieller Differenzialgleichungen

Grundlagen

Differenzialgleichungen liegen den meisten prozess-basierten Modellen der Umweltnaturwissenschaften zugrunde. So beschreiben wir den Wasserfluss und den Stofftransport anhand partieller Differenzialgleichungen. Diese Gleichungen sind oft so komplex, dass eine analytische Lösung nur in sehr speziellen Fällen möglich ist. In diesem Abschnitt wollen wir uns deshalb Verfahren erarbeiten, welche die Lösung von partiellen Differenzialgleichungen numerisch approximieren können.

Wir sind Umweltnaturwissenschaftler:innen und die folgende Darstellung soll grundlegende Konzepte zweier numerischer Verfahren für die Umweltnaturwissenschaften erläutern. Die mathematische Stringenz ist an manchen Stellen vernachlässigt worden, da es sonst den Rahmen dieses Kurses sprengen würde. Wir empfehlen (LeVeque, 2007)° und (LeVeque, 2002)° für Studierende, die Numerik weiter vertiefen möchten.

Bei der Lösung von Differenzialgleichungen suchen wir nach einer Funktion \(u\), welche die Gleichung erfüllt. Dies ist anders als in algebraischen Gleichungen, wo wir nach Werten suchen, welche die algebraische Gleichungen erfüllen. Die Numerik versucht Differenzialgleichungen durch algebraische Gleichungen zu ersetzen und somit die Funktion \(u\) mit diskreten Werten zu beschreiben.

In diesem Kurs betrachten wir gitter-basierte Verfahren, in denen das Lösungsgebiet in diskrete Punkte oder Zellen aufgeteilt wird. An jedem diskreten Punkt oder in jeder diskreter Zelle wird die Funktion \(u\) approximiert. Unsere kontinuierliche Funktion wird somit durch diskrete Punkte im Raum dargestellt. Dies wird räumliche Diskretisierung genannt und ist der erste Schritt zu einer algebraischen Gleichung. Falls eine zeitliche Ableitung vorliegt, muss diese nun in einem zweiten Schritt auch numerisch beschrieben werden. Wir behandeln die Zeit hierbei jedoch als kontinuierliche Größe. Diesen Schritt nennen wir zeitliche Integration. Diese beiden Schritte ergeben das numerische Verfahren zur Lösung der Differenzialgleichung.

Allgemein fordern wir von unseren numerischen Verfahren folgende Eigenschaften: (i) Stabilität, (ii) Konsistenz und (iii) Konvergenz. Diese Eigenschaften werden über die numerischen Fehler eines Verfahrens definiert. Wir unterscheiden zwischen einem globalen Fehler—Differenz zwischen der exakten und numerischen Lösung im gesamten Lösungsgebiet—und einem lokalen Fehler—Differenz zwischen der exakten und numerischen Lösung an einem Punkt.

Ein numerisches Verfahren ist stabil, wenn der globale Fehler \(E\)

\begin{equation*} \lVert E^h \rVert \leq C \lVert \tau^h \rVert \end{equation*}

erfüllt, wobei \(h\) die Diskretisierungslänge, \(C\) eine Konstante und \(\tau\) der lokale Fehler unseres Verfahrens ist.

Das Verfahren ist konsistent mit der zu approximierenden Differentialgleichung, wenn die Bedingung

\begin{equation*} \lVert \tau^h \rVert \rightarrow 0 \, \text{ wenn } \, h \rightarrow 0 \end{equation*}

erfüllt ist.

Das Verfahren konvergiert, wenn es stabil und konsistent ist:

\begin{equation*} \lVert E^h \rVert \rightarrow 0 \, \text{ wenn } \, h \rightarrow 0 \end{equation*}

Finite Differenzen Methode

Wir halten uns an die Einleitung von LeVeque (2007)° und setzen uns die Approximation von \(u'(\bar{t})\), also der ersten Ableitung der Funktion \(u(t)\) an einem Punkt \(\bar{t}\) zum Ziel. Diese Approximation soll nur von der Funktion \(u\) selbst in der Nähe von \(\bar{t}\) abhängen.

Ein erster Ansatz wäre

\begin{equation*} D_+ u(\bar{t}) \equiv \frac{u(\bar{t} + h) - u(\bar{t})}{h}, \end{equation*}

wobei \(h\) eine sehr kleine Distanz ist. Wir sehen dass im Falle \(h \rightarrow 0\), diese Gleichung die Standarddefinition der Ableitung wird, die wir in Analysis gelernt haben. Das gleiche kann natürlich auch in die andere Richtung formuliert werden als

\begin{equation*} D_- u(\bar{t}) \equiv \frac{u(\bar{t}) - u(\bar{t} - h)}{h}, \end{equation*}

oder als ungerichtete Differenz der Form

\begin{equation*} D_0 u(\bar{t}) \equiv \frac{u(\bar{t} + h) - u(\bar{t} - h)}{2h} = \frac{1}{2} \left(D_+ u(\bar{t}) + D_- u(\bar{t})\right). \end{equation*}

Zweite Ableitungen \(u''(\bar{t})\) können mit

\begin{equation*} D^2 u(\bar{t}) \equiv \frac{1}{h^2} \left( u(\bar{t} - h) - 2 u(\bar{t}) + u(\bar{t} + h) \right) \end{equation*}

approximiert werden.

Approximationen dieser Art, in der Differenziale durch Differenzen ersetzt werden, nennen wir Finite Differenzen Methode. Die Herleitung dieser Approximationen kann über eine Taylorreihenentwicklung um den Punkt \(\bar{t}\) erfolgen. Das Restglied der Taylorreihenentwicklung gibt dann Information über die Höhe des lokalen Fehlers der Approximation.

Taylorreihenentwicklung zur Herleitung von \(D_+ u\)

Als Beispiel sei hier die sogenannte Vorwärtsdifferenz \(D_+ u\) hergeleitet. Eine Taylorreihenentwicklung° um \(\bar{t}\) ergibt

\begin{equation*} u(\bar{t} + h) = u(\bar{t}) + \frac{d_t u(\bar{t})}{1!} h + \frac{d_t^ 2u(\bar{t})}{2!} h^2 + \dots + \frac{d_t^n u(\bar{t})}{n!} h^n + R_n(t), \end{equation*}

mit \(R_n\) als das n-te Restglied. Wir schneiden ab der zweiten Ableitung ab und erhalten

\begin{equation*} u(\bar{t} + h) = u(\bar{t}) + \frac{d_t u(\bar{t})}{1!} h + R_1(t). \end{equation*}

\(D_+ u(\bar{t})\) ergibt sich nun durch Umformung:

\begin{equation*} D_+ u(\bar{t}) \approx \frac{u(\bar{t} + h) - u(\bar{t})}{h} \end{equation*}

Weitere Herleitungen von \(D_-\) und \(D_0\) erfolgen sinngemäß.

Beispiel: Lösung einer gewöhnlichen Differenzialgleichung

Wir möchten die gewöhnliche Differenzialgleichung

\begin{equation*} d_t u(t) = \cos(t) \end{equation*}

mit der Anfangsbedingung \(u(t=0) = 0\) numerisch approximieren. Damit wäre die exakte allgemeine Lösung welche wir suchen \(u(t) = \sin (t) + C\), und durch die Anfangsbedingung determinierte spezielle Lösung \(u(t) = \sin (t)\).

Approximation von \(u(t)\) an einem Punkt \(t_i\) durch die Finite Differenzen Methode lautet

\begin{equation*} d_t u_i \approx D_+ u_i = \frac{u_{i+1} - u_i}{\Delta t}. \end{equation*}

Somit ergibt sich die diskrete Gleichung

\begin{equation*} \frac{u_{i+1} - u_i}{\Delta t} = \cos (t_i), \end{equation*}

welche umgestellt werden kann zu

\begin{equation*} u_{i+1} = u_i + \Delta t \cos (t_i). \end{equation*}

In Julia implementieren wir dies wie folgt:

t[1] = t0
u[1] = u0

for i in 2:n
    t[i] = t[i-1] + Δt
    u[i] = u[i-1] + Δt * cos(t[i])
end

Das Ergebnis—diskretisiert mit 100 Punkten—ist in Abb. 1 zu sehen. Wir sehen, dass die Approximation von der exakten Lösung abweicht. Eine Diskretisierung mit einem kleineren \(\Delta t\), also mit mehr Punkten, würde eine bessere Approximation ergeben.

pde_100.png

Abbildung 1: 01_ode_fdm.jl: Finite Differenzen Methode Lösung für \(d_t u = \cos (t)\).

Weitere Ressourcen

Beispielhafte Implementierungen der Finite Differenzen Methode für kanonische partielle Differentialgleichungen sind in Barba et al. (2024)° gut zusammengefasst.

Finite Volumen Methode

Für viele strömungsmechanische Anwendungen hat sich die Finite Volumen Methode als Diskretisierungsmethode durchgesetzt. Im Gegensatz zur Finite Differenzen Methode approximiert die Finite Volumen Methode die Integralform der Gleichungen, nähert also eine schwache Lösung an.

control-volume_c.png

Abbildung 2: Ein dreidimensionales Kontrolvolumen. f1 bis f6 repräsentieren Flüsse in das bzw. aus dem Kontrolvolumen.

Eine anschauliche Herleitung der Methode basiert auf dem Konzept eines Kontrolvolumen, wie in Abb. 2 beispielhaft aufgezeichnet. Wir betrachten nun eine Größe \(e\) innerhalb dieses Kontrolvolumens. Die zeitliche Veränderung des räumlich aggregierten \(e\) kann durch eine Bilanz der Flüsse über den Rand dieses Kontrolvolumens beschrieben werden. In dem in Abb. 2 skizzierten Fall könnten wir also

\begin{equation*} d_t \int_\Omega e d\Omega = \sum_{i=1}^6 f_i n_i A_i \end{equation*}

schreiben, wobei die Flüsse \(f_i\) von der Größe \(e\) abhängen können und \(n\) der Normalenvektor ist, definiert als positiv wenn er aus dem Kontrolvolumen herauszeigt. Das Integral aggregiert \(e\) im ganzen Kontrolvolumen \(\Omega\). \(A_i\) ist die Fläche der Kanten. Im allgemeineren Fall wird die Summe zu einem Integral über den Rand:

\begin{equation*} d_t \int_\Omega e d \Omega = \oint_{\Gamma} f n d\Gamma \end{equation*}

Diese Art von Gleichungen nennen wir auch Erhaltungssätze, da sie den Erhalt einer Größe \(e\) innerhalb eines Kontrolvolumens beschreiben.

Im einfachsten Fall ist \(e\) über das Kontrolvolumen konstant. Das Integral über \(\Omega\) kann dann direkt evaluiert werden als

\begin{equation*} \Omega d_t e = \oint_{\Gamma} f n d\Gamma. \end{equation*}

Die Grundidee der Finite Volumen Methode ist nun das Gebiet in eine finite Anzahl von Kontrolvolumen \(j\) zu zerlegen. Dies ist im Gegensatz zu den einfachen Punkten, welche wir in der Finite Differenzen Methode zur Diskretisierung hatten. Unter Annahme dass \(f\) konstant über die Kante \(i\) ist, wird das Integral über den Rand nun wieder zur Summe über die Zellkanten \(i\), also

\begin{equation} \label{orgb688fbe} d_t e_j = \frac{1}{\Omega_j} \sum_i^N f_i n_i A_i. \end{equation}

Über die Flüsse über die Kanten der einzelnen Zellen fließt \(e\) nun von einer Zelle in die nächste und in jeder Zelle kann eine lokale Bilanz von \(e\) nach Gl. \eqref{orgb688fbe} aufgestellt werden.

Die zeitliche Integration kann nun z.B. über eine Vorwärtsdifferenz gelöst werden, was den Ausdruck

\begin{equation*} d_t e_j \approx D_+ e_j = \frac{e^{n+1}_j - e^n_j}{\Delta t} = \frac{1}{\Omega_j} \sum_i^N f_i n_i A_i \end{equation*}

ergibt. Hier steht \(n\) für den Index eines zeitlichen Schrittes. Schlussendlich erhalten wir nach Umstellung die komplett algebraische Form

\begin{equation*} e^{n+1}_j = e^n_j + \frac{\Delta t}{\Omega_j} \sum_i^N f_i^\eta n_i A_i, \end{equation*}

mit \(\eta \in \{n, n+1\}\) als ein Zeitschrittindex der gewählt werden muss. Wenn \(\eta = n\), nennen wir die zeitliche Integration explizit. Falls \(\eta = n+1\), nennen wir sie implizit.

\(f_i^a\) nennen wir einen numerischen Fluss, der an der Kante \(i\) der Zelle \(j\) berechnet werden muss. Hier bieten sich oft mehrere Möglichkeiten an, welche die Stabilität und die Konvergenz des Verfahrens beeinflussen.

Beispiel: Lösung einer partiellen Differenzialgleichung

Wir möchten die partielle Differenzialgleichung

\begin{equation*} \partial_t u = - \partial_x (v u) \end{equation*}

für \(v > 0\) und \(x \in [0, 1]\) mit periodischen Randbedingungen und der Anfangsbedingung

\begin{equation*} u(t = 0, x) = \begin{cases} 1, & x < 0.2\\ 0, & \text{sonst} \end{cases} \end{equation*}

numerisch approximieren. Wir zerlegen gedanklich das Gebiet in eine finite Anzahl von Zellen der Länge \(\Omega_j\) und integrieren die obige Gleichung über \(\Omega_j\). Wir erhalten

\begin{equation*} \int_\Omega \partial_t u d\Omega = - \int_\Omega \partial_x (v u) d\Omega. \end{equation*}

Wir benutzen den Satz von Green–Gauß° um das Flächenintegral in ein Kurvenintegral über den Rand zu transformieren:

\begin{equation*} \int_\Omega \partial_t u d\Omega = - \oint_\Gamma v u n d\Gamma \end{equation*}

Wir sehen dass dies der Gl. \eqref{orgb688fbe} entspricht, wobei \(f = vu\). Wir setzen nun Voraus, dass \(u\) im Kontrollvolumen \(j\) durch ein konstantes Mittel \(u_j\) repräsentiert werden kann und erhalten

\begin{equation*} \Omega_j d_t u_j = - \oint_\Gamma v u n d\Gamma. \end{equation*}

Des weiteren ersetzen wir nun das Kurvenintegral durch eine diskrete Summe um

\begin{equation*} \Omega_j d_t u_j = - \sum_i^N v_i u_i n_i A_i \end{equation*}

zu bekommen. Im eindimensionalen Fall haben wir nun als "Zelle" eine Linie und zwei "Kanten"—eigentlich Punkte an beiden Enden einer Linie. Wir schreiben die obige Gleichung noch einmal als

\begin{equation*} \Delta x d_t u_j = v_{j-1/2} u_{j-1/2} - v_{j+1/2} u_{j+1/2}. \end{equation*}

Wie oben besprochen, benutzen wir eine Vorwärtsdifferenz für die zeitliche Integration und schreiben ein explizites Finite Volumen Verfahren als

\begin{equation*} u_j^{n+1} = u_j^{n} + \frac{\Delta t}{\Delta x} \left( v_{j-1/2} u_{j-1/2} - v_{j+1/2} u_{j+1/2}\right)^n. \end{equation*}

Was benutzen wir nun als \(f_{j\pm1/2} = v_{j\pm1/2} u_{j\pm1/2}\)? Eine recht gängige Methode ist das sogenannte upwinding, in dem wir den Fluss basierend auf dem Wert stromaufwärts berechnen. Da in unserem Fall \(v > 0\) gilt, berechnen wir \(f_{j+1/2} = f(u_j)\).

In Julia, liest sich die Methodik wie folgt:

t = 0.0 # Zeit
while t < tEnd
    global un = copy(u)
    for i in 2:n
        global un[i] = u[i] + (Δt / Δx) * (v * u[i-1] - v * u[i])
    end
    # Randbedingungen
    un[1] = u[1] + (Δt / Δx) * (v * u[n] - v * u[1])
    # Fertig.  Nur noch das Original ersetzen und den Zeitschritt
    # aktualisieren
    global u = copy(un)
    if t + Δt > tEnd
        global Δt = tEnd - t
    end
    global t = t + Δt
end

Das Ergebnis—diskretisiert mit 201 Zellen—ist in Abb. 3 zu sehen. Wir sehen, dass die Approximation von der exakten Lösung abweicht. Eine Diskretisierung mit einem kleineren \(\Delta x\), also mit mehr Zellen, würde eine bessere Approximation ergeben.

pde_201c.png

Abbildung 3: 02_pde_fvm.jl: Finite Volumen Methode Lösung von \(\partial_t u = - \partial_x (v u)\).

Linienmethode

Die Linienmethode° ist ein Verfahren zum Lösen von partiellen Differenzialgleichungen, bei dem alle bis auf eine Dimension diskretisiert werden. Dies reduziert die partielle Differenzialgleichung zu einer gewöhnlichen Differenzialgleichung.

Modellierung des Bodenwasserflusses

Richards(on) Gleichung

Die Richards(on) Gleichung beschreibt den Bodenwasserfluss unter teilgesättigten Bedingungen. Die Gleichung kann in verschiedenen Formen geschrieben werden, jede mit eigenen Vor- und Nachteilen (Caviedes-Voullième et al., 2013)°.

In sogenannter Potentialform liest die eindimensionale Richards(on) Gleichung

\begin{equation*} \rho_l C(\psi) \partial_t \psi = \nabla \left( K(\psi) \nabla(\psi + gz) \right), \end{equation*}

wobei \(\psi~[\mathrm{J/kg}]\) das Matrixpotential, \(\rho_l~[\mathrm{kg/m^3}]\) die Dichte des Wassers, \(K~[\mathrm{kgs/m^3}]\) die hydraulische Leitfähigkeit (im mehrdimensionalen Fall ein Tensor), \(g~[\mathrm{m/s^2}]\) ist die Erdbeschleunigung und \(z~[\mathrm{m}]\) die geodätische Höhe ist. \(C~[\mathrm{kg/J}]\) ist die sogenannte hydraulische Kapazität, welche definiert ist als

\begin{equation*} C(\psi) = \partial_\psi \theta, \end{equation*}

und \(\theta\) die Bodenfeuchte ist. Die Potentialform hat für numerische Verfahren einige Vorteile, ist jedoch nicht per se massenkonservativ.

In Wassergehaltform ist die eindimensionale Richards(on) Gleichung

\begin{equation*} \partial_t \theta = \nabla \left( D(\theta) \nabla \theta \right) + \nabla K(\theta), \end{equation*}

wo \(D~[\mathrm{Js/m^3}]\) die Diffusivität ist. Diese Form der Richards(on) Gleichung ist massenkonservativ, kann jedoch keine Gradienten im gesättigten Boden abbilden. Heterogenitäten in bodenhydraulischen Kenngrößen würden in dieser Form zu Diskontinuitäten in der Lösung führen, was für numerische Verfahren ein Nachteil sein kann.

Als letztes besprechen wir nun die gemischte Form der Richards(on) Gleichung. In eindimensionaler Form liest sich diese als

\begin{equation} \label{org0b90e2f} \partial_t \theta = \nabla \left( K(\psi) \nabla (\psi + z) \right). \end{equation}

Als Kombination von beiden Formen, erlaubt die gemische Form die massenkonservative Modellierung des gesättigten Bodenwasserflusses. Die Gleichung muss jedoch für zwei Unbekannte — \(\psi\) und \(\theta\) — gelöst werden.

Ein Finite Volumen Modell für die Richards(on) Gleichung

Wir integrieren Gl. \eqref{org0b90e2f} über ein eindimensionales Kontrolvolumen der Länge \(\Delta z\):

\begin{equation*} \int_\Omega \partial_t \theta d\Omega = \int_\Omega \partial_z \left( K(\psi) \partial_z (\psi + z) \right) d\Omega \end{equation*}

Nach den üblichen Umformungen erhalten wir

\begin{equation*} \partial_t \theta = \frac{1}{\Delta z} \oint_\Gamma K(\psi) D_0 (\psi + z) n_\Gamma d\Gamma, \end{equation*}

wobei die räumliche Ableitung innerhalb des Integrals durch eine zentrale Differenz \(D_0 (\psi + z)\) ersetzt wurde.

Bodenmodelle

Mualem–van Genuchten Bodenmodell

Die Beziehung zwischen \(\psi\) und \(\theta\), sowie die Abhängigkeit von \(K\) von diesen beiden Variablen, wird durch Bodenmodelle beschrieben. Das wahrscheinlich am häufigsten benutzte Bodenmodell ist das Mualem–van Genuchten (MvG) Bodenmodell (van Genuchten, 1980)°.

Dieses Modell beschreibt \(\theta\) als Funktion von \(\psi\) als

\begin{equation*} \begin{aligned} \theta(\psi) = \begin{cases} \frac{\theta_s - \theta_r}{(1 + (\alpha \vert \psi \vert)^n)^m} + \theta_r & \psi \leq 0,\\ \theta_s & \psi > 0, \end{cases}\\ m = 1 - \frac{1}{n}, \end{aligned} \end{equation*}

wobei \(\alpha\) und \(n\) die beiden Modellparameter sind. \(\alpha~[\mathrm{1/m}]\) beschreibt den Lufteintritt in den Boden und \(n\) steht in Beziehung mit der Porendurchmesserverteilung im Boden. Diese Beziehung wird auch Wasserretentionskurve genannt.

Die hydraulische Leitfähigkeit \(K\) ist

\begin{equation*} \begin{aligned} K = \begin{cases} K_s S_e^L \left( 1 - (1 - S_e^{\frac{1}{m}})^m \right)^2 & \psi \leq 0, \\ K_s & \psi > 0, \end{cases}\\ S_e = \frac{\theta - \theta_r}{\theta_s - \theta_r}, \end{aligned} \end{equation*}

wobei \(L\) ein freier Parameter ist, der oft als \(L=0.5\) gewählt wird.

Peters–Durner–Iden Bodenmodell

Das Braunschweiger Peters–Durner–Iden (PDI) Bodenmodell (Peters, 2013)° hat gegenüber dem MvG den Vorteil, dass es das gesamte Intervall der möglichen Bodenfeuchte abbilden kann. Das MvG hat Schwierigkeiten, bodenhydraulische Eigenschaften unter sehr trockenen Bedingungen abzubilden. PDI liefert keine Prozessbeschreibung für verschiedenen Komponenten, sondern benutzt bestehende Bodenmodelle.

PDI beschreibt die Wasserretention als eine Kombination aus Wasserrückhalt in vollständig gefüllten Kapillaren und an der Oberfläche von Feststoffen. Die Retentionskurve lautet

\begin{equation*} \theta(\psi) = \theta_c(\psi) + \theta_{nc}(\psi) = (\theta_s - \theta_r) S_c (\psi) + \theta_r S_{nc} (\psi), \end{equation*}

mit \(\theta_c\) als die Bodenfeuchte in den gefüllten Kapillaren und \(\theta_{nc}\) als die Bodenfeuchte die an der Oberfläche von Feststoffen adsorbiert wird. Der kapillare Speicher \(S_c\) ist

\begin{equation*} S_c(\psi) = \frac{U(\psi) - U_0}{1 - U_0}, \end{equation*}

wobei \(U_0\) die Sättigung bei Ofentrockenheit. Wenn das MvG Bodenmodell benutzt wird, lautet \(U\)

\begin{equation*} U(\psi) = \left(1 + (\alpha \psi)^n \right)^{\frac{1}{m}}. \end{equation*}

Die gesamte hydraulische Leitfähigkeit des Bodens wird als Folge von drei Strömungsprozessen betrachtet: (i) Kapillarströmung, (ii) Filmströmung, (iii) isothermische Dampfströmung.

Modellierung des Stofftransports in Böden

Advektions–Dispersions-Gleichung

Die Advektions–Dispersions-Gleichung beschreibt konservativen Stofftransport in Böden. In eindimensionaler Form schreibt sie sich als

\begin{equation*} \partial_t C_l + \partial_z v_f C_l - \partial_z \left( D \partial_z C_l \right) = 0, \end{equation*}

wobei \(t~[\mathrm{s}]\) die Zeit und \(z\) die räumliche Koordinate in der vertikalen Richtung sind. \(C_l~[\mathrm{mg/m^3}]\) ist die gelöste Konzentration des transportierten Stoffes und \(D~[\mathrm{m^2/s}]\) ist der Dispersionskoeffizient.

Reaktiver Transport

Reaktiver Transport bezeichnet Stofftransport und die dadurch ausgelösten chemischen Reaktionen in Böden. Wir betrachten zwei Prozesse des reaktiven Transports: (i) Sorption und (ii) Abbau.

Sorption

Sorption ist die Anreicherung eines Stoffes innerhalb einer Phase (Absorption) oder einer Grenzfläche (Adsorption). Der sorbierende Stoff wird Sorbent genannt, das Gesamtsystem aus Sorbat und dem sorbiertem Stoff wird Sorbat genannt. Sorption ist ein Gleichgewichtsprozess zwischen Aufnahme und Abgabe eines Stoffes.

sorption_sketch.png

Abbildung 4: Konzeptskizze von Sorption an der Grenzfläche des Bodens

Im Boden wird vor allem Adsorption an der Grenzfläche zwischen Feststoff und Fluid—also Bodenpartikel und Wasser—beobachtet. Die spezifische Oberfläche° ist hierbei der wichtigste Parameter für Adsorption. Je größer die spezifische Oberfläche, desto mehr Adsorptionsplätze stehen zur Verfügung. Es gilt, dass die spezifische Oberfläche mit abnehmendem Korndurchmesser steigt—siehe Tab. 1. Eine Konzeptskizze für Sorption im Boden ist in Abb. 4 gegeben.

Tabelle 1 Typische Werte für spezifische Oberflächen je Bodenart nach Hillel (1998)°
Material spezifische Oberfläche (m2/g)
Sand > 0.1
Schluff 0.1–1
Ton 5–500
organische Substanz 800–1200

Sorptionsisotherme

Sorptionsisotherme sind Gleichgewichtszustände der Sorption eines bestimmten Stoffes bei konstanter Temperatur. Als solches, Sorptionsisotherme sind sowohl vom sorbiertem Stoff als auch von der Temperatur abhängig. Für eine detaillierte Übersicht, siehe Limousin et al. (2007)°.

Im Allgemeinen kann eine Sorptionsisotherme als eine Funktion der Form

\begin{equation*} S = f(C_l) \end{equation*}

geschrieben werden, wobei \(S~[\mathrm{mg/kg}]\) die Konzentration des sorbierten Stoffes in Relation zur Masse der Festsubstanz des Bodens und \(C_l~[\mathrm{mg/m^3}]\) die gelöste Konzentration des Stoffes in Relation zum Wasservolumen im Boden ist.

Giles et al. (1974) definieren vier Hauptformen von Isothermen, die in der Abb. 5 zusammengefasst sind.

giles_1974_from_limousin_2007.png

Abbildung 5: Die vier Hauptformen von Isothermen, reproduziert von Giles et al. (1974)

Die "C" Isotherme—Abb. 5 (a)—kann durch ein lineares Modell (Henry Modell) in der Form

\begin{equation*} S = K_d C_l \end{equation*}

abgebildet werden. Hierbei ist \(K_d~[\mathrm{m^3 kg}]\) der Verteilungskoeffizient.

"L" und "H" Isotherme—Abb. 5 (b) und (c), resp.—können mit Freundlich Modellen abgebildet werden. Dies sind auch die am häufigsten beobachteten Isotherme. Das einfache Freundlich Modell ist

\begin{equation*} S = K_f C_l^n, \end{equation*}

wobei \(K_f~[\mathrm{m^3 kg}]\) der Freundlich-Parameter und \(n\) ein einheitsloser Modellparameter ist.

Ein weiteres Modell welches "L" und "H" Isotherme—Abb. 5 (b) und (c), resp.—beschreibt ist das einfache Langmuir Modell der Form

\begin{equation*} S = S_{\mathrm{max}} \frac{K_l C_l}{1 + K_l C_l}, \end{equation*}

wobei \(S_{\mathrm{max}}~[\mathrm{mg/kg}]\) die maximale Konzentration und \(K_l\) die Langmuir-Konstante sind. \(K_l\) beschreibt die Affinität eines Stoffes zum Sorbat.

Die "S" Isotherme resultiert aus der Interaktion zweier entgegengesetzter Mechanismen. Ein Beispiel hierzu sind unpolare organische Verbindungen (Limousin et al., 2007)°. Solche Verbindungen haben eine geringe Affinität zu Tonen. Sobald aber eine Tonoberfläche mit diesen Verbindungen bedeckt ist, werden andere organische Moleküle leichter adsorbiert.

Advektions–Dispersions-Gleichung mit Sorption

Die Advektions–Dispersions-Gleichung (§Advektions–Dispersions-Gleichung) kann um einen Sorptionsterm erweitert werden. Damit ergibt sich

\begin{equation} \label{org5a711a9} \partial_t \theta C_l + \partial \rho_b S + \partial_z v_f \theta C_l - \partial_z \left( D \partial_z \theta C_l \right) + R = 0, \end{equation}

wobei \(\theta\) der Volumenanteil der Fluidphase und \(\rho_b~[\mathrm{kg/m^3}]\) die Lagerungsdichte des Bodens sind. \(S\) wird anhand von Isotherme Modellen als Funktion von \(C_l\) berechnet. \(R~[\mathrm{mg/m^3s}]\) ist ein Reaktionsterm.

Retardationsfaktor

Der Retardationsfaktor \(r\) kann unter Annahme eines konstanten \(\rho_b\) mit der Kettenregel hergeleitet werden. Unter diesen Annahmen, kann \(\rho_b\) aus dem zweiten Term der linken Seite in Gl. \eqref{org5a711a9} aus der zeitlichen Ableitung herausgezogen werden, was

\begin{equation*} \partial_t \rho_b S(C_l) = \rho_b \partial_t S(C_l) \end{equation*}

ergibt. Nun kann die Kettenregel auf \(S(C_l)\) angewandt werden um

\begin{equation*} \rho_b \partial_t S(C_l) = \rho_b \partial_{C_l} S \, \partial_t C_l \end{equation*}

zu erhalten. Unter der Annahme eines konstanten \(\theta\), können die ersten beiden Terme der linken Seite in Gl. \eqref{org5a711a9} noch weiter vereinfacht werden. Wir teilen alles durch \(\theta\) und erhalten

\begin{equation*} \partial_t \theta C_l + \partial \rho_b S = \partial_t C_l + \frac{\rho_b}{\theta} \partial_{C_l} S \, \partial_t C_l. \end{equation*}

Hier können wir die \(\partial_t C_l\) Terme ausklammern und kriegen

\begin{equation*} \partial_t C_l + \frac{\rho_b}{\theta} \partial_{C_l} S \, \partial_t C_l = \left( 1 + \frac{\rho_b}{\theta} \partial_{C_l} S \right) \partial_t C_l = r \partial_t C_l, \end{equation*}

wobei \(r\) der sogenannte Retardationsfaktor ist. Wir erhalten statt Gl. \eqref{org5a711a9}

\begin{equation*} \left\{ \begin{aligned} r \partial_t C_l + \partial_z v_f C_l - \partial_z \left( D \partial_z C_l \right) + \frac{R}{\theta} = 0,\\ r = 1 + \frac{\rho_b}{\theta} \partial_{C_l} S. \end{aligned} \right. \end{equation*}

Im Falle einer linearen Sorptionsisotherme, der Retardationsfaktor ist

\begin{equation*} r = 1 + \frac{\rho_b}{\theta} K_d \end{equation*}

und im Falle des Freundlich Modells

\begin{equation*} r = 1 + \frac{\rho_b}{\theta} K_f n C^{n-1}. \end{equation*}

Nichtlineare Sorption wie das Freundlich Modell führt zu selbstverschärfenden Fronten—siehe Abb. 6.

sorption.png

Abbildung 6: Beispielrechnung mit nichtlinearer Sorption nach Freundlich mit \(K_f = 0.15~\mathrm{m^3 kg}\), \(n = 3\), \(\rho_b = 1.2~\mathrm{kg/m^3}\) und \(\theta = 0.4\). Transportparameter sind \(v_f = 10^{-3}~\mathrm{m/s}\) und \(D = 10^{-5}~\mathrm{m^2/s}\).

Rechenbeispiel

In einem sandigen Boden mit Lagerungsdichte \(\rho_b = 1.5~\mathrm{g/cm^3}\) und mittlerem Wassergehalt \(\theta = 0.1\) und Grundwasserflurabstand von \(100~\mathrm{cm}\) wurde vom 1.11. bis 31.4. (0.5 a) eine Grundwasserneubildung von \(10~\mathrm{cm}\) gemessen. Der Henry-Koeffizient des Systems ist \(K_d = 1~\mathrm{cm^3/g}\).

  1. Wie groß ist der Retardationsfaktor \(r\)?
  2. Wie groß ist die Abstandsgeschwindigket (\(v^* = v/r\))?
  3. Wie lange braucht ein infiltrierender Stoff bis zum Grundwasser?

Abbau

Ein weiterer wichtiger Prozess des reaktiven Transports ist der Abbau von Substanzen im Boden. So spielt der Abbau von Stickstoff und Kohlenstoff im Boden eine wichtige Rolle im globalen Kreislauf dieser Elemente—siehe, z.B. Porporato and Yin (2022)°.

Der Abbau von einem Stoff wird in Gl. \eqref{org5a711a9} im Term \(R\) berücksichtigt. Reaktionsterme verschiedener Ordnung können benutzt werden, um den Abbau eines Stoffes zu beschreiben. Ein Abbauterm der 0-ten Ordnung ist

\begin{equation*} R = - \gamma, \end{equation*}

wobei \(\gamma~[\mathrm{kg/(cm^{3}d)}]\) der Abbauparameter ist—siehe Abb. 7. Ein Abbauterm 1-ster Ordnung wäre

\begin{equation*} R = - \mu C_l, \end{equation*}

wobei \(\mu~[\mathrm{1/d}]\) der Abbauparameter ist—siehe Abb. 8.

abbau_0.png

Abbildung 7: Abbau 0-ter Ordnung

abbau_1.png

Abbildung 8: Abbau 1-ster Ordnung

Modellierung des Transports stabiler Wasserisotope in Böden

Stabile Wasserisotope werden in der Ökohydrologie gerne als natürliche Tracer eingesetzt, die Auskunft über die Herkunft und das Alter des Wassers geben können. Als Beispiel für isotopen-gestützte Ökohydrologie verweisen wir auf unsere Schwester-Arbeitsgruppe isodrones° an der Technischen Universität Braunschweig.

Die Transportmodellierung stabiler Wasserisotope beinhaltet die Darstellung von zwei Prozessen (Kuppel et al., 2018)°: (i) Durchmischung und (ii) Fraktionierung.

Wurzelwasseraufnahme

Einleitung

Die Wurzeln verankern die Pflanzen im Boden und nehmen Wasser und Nährstoffe aus dem Boden auf und leiten diese Stoffe dann nach oben zum Stamm (Nobel, 2020)°. An der äußersten Spitze einer Wurzel befindet sich die Wurzelhaube. Die Zellwände der Wurzelhaube sind—um die Reibung mit Bodenpartikeln verringern—oft schleimig. Neben der Wurzelhaube befindet sich eine Region, in der schnelle Zellteilung stattfindet, was die Wurzelhaube mechanisch durch den Boden schiebt.

flach-wurzler.jpg

Abbildung 9: Flach wurzelnder Baum im Prinz-Albrecht-Park, Braunschweig (gemäßigtes Klima)

Die Architektur der Wurzel, also ihre Länge, Verästelung, Ausbreitung und Neigung ist von der Wasserverfügbarkeit des Standortes geprägt. Unter trockenen Randbedingungen investiert die Pflanze eher in das Wurzelsystem als in das Sprossystem (Hsiao & Yu, 2000)°. Je moderater das Umfeld der Pflanze, desto mehr Biomasse wird in das Sprossystem allokiert—siehe, z.B. Abb. 9. Deshalb korreliert zum Beispiel die Wurzeltiefe auch deutlich mit der Ökozone und schwach mit dem Blattflächenindex (LAI)—siehe Abb. 10.

lai-depth.png

Abbildung 10: Blattflächenindex und Wurzeltiefe in Abhängigkeit von der Ökozone nach dem Datensatz von Jackson et al. (1997)°.

Wurzelwasseraufnahme (WWA) erfolgt über die Wurzelhaare und feinen Seitenwurzeln. Ältere Wurzeln dienen oft nur der Wasserleitung. WWA wird vom Potentialunterschied zwischen Wurzel und Boden angetrieben. Das Wasserpotential in den Wurzeln ist in der Regel negativer als im Boden. WWA ist proportional zur aktiven Wurzelfläche, hängt deshalb auch von der vertikalen Verteilung der Wurzeln, also der Wurzelarchitektur, im Boden ab. Jedoch ist eine hohe Wurzelbiomasse in einer bestimmten Bodenschicht kein eindeutiger Beleg, dass aus dieser Bodenschicht tatsächlich viel Wasser aufgenommen wird (Lobet et al., 2014)°.

Trotz zahlreicher Untersuchungen auf multiplen Skalen—siehe, z.B. Maurel et al. (2008)°; Steudle (2001)°; Jarvis (2011)°—ist der Prozess der WWA immer noch nicht vollständig verstanden, siehe auch Carminati et al. (2013)°.

Modellierung der Wurzelwasseraufnahme

Für den Transport von Wasser in Böden ist die WWA ein Senkenterm in der Richards(on) Gleichung, den wir mit \(S\) notieren als

\begin{equation*} \partial_t \theta = \nabla \left( K(\psi) \nabla (\psi + z) \right) - S(\psi). \end{equation*}

Im folgenden betrachten wir uns die Berechnung dieses Senkenterms \(S\) genauer. Die Berechnung von \(S\) erfolgt meistens in fünf Schritten:

  1. Berechnung der potenziellen Evapotranspirationsrate für einen bestimmten Bestand (\(ET_c~[\mathrm{cm/d}]\)) aus der Grasreferenzverdunstung (\(ET_0~[\mathrm{cm/d}]\))
  2. Partitionierung von \(ET_c\) in potenzielle Transpirationsrate und potenzielle Evaporationsrate
  3. Verteilung der potenziellen Transpirationsrate zur potenziellen WWA auf die Bodentiefe
  4. Reduktion der potenziellen WWA auf eine aktuelle WWA durch eine Stressfunktion
  5. Kompensation für den Fall, dass in bestimmten Tiefen wenig Wasser, aber in anderen Tiefen genügend Wasser vorhanden ist

Berechnung der Grasreferenzverdunstung nach FAO

Die Formel für Grasreferenezverdunstung der Food and Agriculture Organization of the United States (FAO) ist eine Vereinfachung der Penman–Monteith Verdunstungsgleichung unter Annahme einer Referenzfläche. Die Penman–Monteith Gleichung ist

\begin{equation*} \lambda ET = \frac{\Delta (R_n - G) + \rho_a c_p \frac{e_s - e_a}{r_a}}{\Delta + \gamma \left( 1 + \frac{r_s}{r_a} \right)}, \end{equation*}

wobei \(\lambda~[\mathrm{J/kg}]\) die latente Wärme, \(\Delta~[\mathrm{Pa/°K}]\) die Steigung der Dampfdruckkurve, \(R_n~[\mathrm{W/m^2}]\) die Nettostrahlung, \(G~[\mathrm{W/m^2}]\) der Bodenwärmefluss, \(\rho_a~[\mathrm{kg/m^3}]\) die Luftdichte, \(c_p~[\mathrm{J/kg°K}]\), \(r_a~[\mathrm{s/m}]\) der aerodynamische Widerstand, \(r_s~[\mathrm{s/m}]\) der Bestandeswiderstand, \(e_s~[\mathrm{Pa}]\) der Sättigungsdampfdruck, \(e_a~[\mathrm{Pa}]\) der aktuelle Dampfdruck und \(\gamma~[\mathrm{Pa/°K}]\) die psychrometrische Konstante ist. Siehe Calanca et al. (2011)° für eine detailliertere Betrachtung.

Die FAO° definiert nun eine Referenzfläche wie folgt:

A hypothetical reference crop with an assumed crop height of \(0.12~\mathrm{m}\), a fixed surface resistance of \(70~\mathrm{s m^{-1}}\) and an albedo of \(0.23\).

Daraus ergibt sich die Grasreferenzverdunstung als

\begin{equation*} ET_0 = \frac{0.408 \Delta (R_n - G) + \gamma \frac{900}{T + 273} u_2 (e_s - e_a)}{\Delta + \gamma (1 + 0.34 u_2)}, \end{equation*}

mit \(u_2\) als die Windgeschwindigkeit gemessen auf \(2~\mathrm{m}\) Höhe.

Um von der Referenzverdunstung \(ET_0\) auf die aktuelle Standverdunstung \(ET_c\) zu kommen, with \(ET_0\) mit einem Bestandskoeffizienten \(k_c\) multipliziert, der die Unterschiede zwischen der Referenzfläche und der betrachteten Fläche berücksichtigt. Die \(k_c\) Werte können Tabellen in der Literatur entnommen werden.

Partitionierung der Bestandsverdunstung in Transpiration und Evaporation

Nachdem \(ET_c\) berechnet wurde, muss es in seine Komponenten Bodenverdunstung \(E_c\) und Transpiration \(T_c\) aufgeteilt werden. Dies geschieht auf Basis des sogenannten Blattflächenindex (LAI), der Blattfläche pro Bodenoberfläche. LAI kann Artbezogen oder Biombezogen angegeben werden—siehe Tab. 2.

Tabelle 2 Beispielhafte LAI Werte für einige Arten und Biome (wiki)°.
Art LAI Biom LAI
Douglasie 10–13 Tropischer Regenwald 6–16
Fichte 5–10 Sommergrüner Laubwald 3–12
Buche 6–8 Nadelwald 5
Eiche 5–7 Savanne 1–5
Waldkiefer 3–4 Tundra 0.5–2.5
Europäische Lärche 2–4 Buchenwald (Winter) 0.2

Die Aufteilung der Evapotranspiration kann nach dem Lambert–Beer Gesetz erfolgen:

\begin{equation*} \left\{ \begin{aligned} E_c = ET_c \exp(- k \mathrm{LAI})\\ T_c = ET_c (1 - \exp(- k \mathrm{LAI})) \end{aligned} \right. \end{equation*}

Verteilung der potenziellen Transpirationsrate über die Bodentiefe

Die räumliche Verteilung der potenziellen Transpirationsrate ist empirisch und kann von der Wurzelarchitektur hergeleitet werden. Oft sind hier jedoch grobe Annahmen erforderlich. Die räumlich verteilte Wurzelwasseraufnahme \(S_c\), die sich aus der Verteilung der Transpirationsrate ergibt, berechnet sich als

\begin{equation*} \left\{ \begin{aligned} S_c(z,t) = T_c(t) b(z), \\ \int_l^0 b(z) dz = 1. \end{aligned} \right. \end{equation*}

Hier ist \(b\) eine normierte Aufnahmeverteilung, welche sicherstellt, das die Summe der aufgeteilten Transpirationsraten die gesamte Transpirationsrate ergibt.

Reduktion der potenziellen Wurzelwasseraufnahme durch eine Stressfunktion

Wenn der Boden trocken fällt sinkt das Bodenwasserpotential und es fällt der Pflanze schwieriger, das Wasser aus dem Boden zu heben. Über längere Zeit verursacht dies in der Pflanze Wasserstress und die in den oberen Abschnitten errechnete Transpirationsrate kann nicht erfüllt werden. Modelliert werden kann dieser Prozess zum Beispiel nach Feddes et al. (1978)° durch eine Reduktion der Transpirationsrate anhand eines Wasserstresskoeffizienten \(\alpha\):

\begin{equation*} S_u(z,t) = \alpha(z,t) S_c(z,t) \end{equation*}

Die Form von \(\alpha\) ist in Abb. 11 gegeben. Mit steigendem Wasserpotential (feuchterem Boden) steigt auch der Wert von \(\alpha\), bis eine optimale Bodenfeuchte erreicht ist (\(h_3\) bis \(h_2\)), bei der die potentielle Transpirationsrate erreicht wird. Wird der Boden weiterhin gesättigt, sinkt \(\alpha\) jedoch aufgrund von Vernässung°, was zu Luftmangel und Reduktionserscheinungen führt. In der Landwirtschaft sehen wir deshalb oft Drainagen, die der Vernässung der Kulturpflanzen vorbeugen sollen.

feddes.jpg

Abbildung 11: Form des \(\alpha\) Koeffizienten. Typische Werte für die Wasserpotentialhöhen sind \(h_1 = 10\), \(h_2 = 25\), \(h_3 = 300-1000\) und \(h_4 = 8000~\mathrm{cm}\). Reproduziert von Feddes et al. (1978)°.

Kompensation

Wenn in bestimmten Bodentiefen wenig Wasser, in anderen Tiefen jedoch genügend Wasser vorhanden ist, kann die Reduktion der Transpirationsrate in den trockenen Tiefen durch die feuchteren Tiefen kompensiert werden.

Die unkompensierte Transpirationsrate errechnet sich als

\begin{equation*} T_u = \int_l^0 S_u dz. \end{equation*}

Wir definieren einen Stressindex \(\omega\) als

\begin{equation*} \omega = \frac{T_u}{T_c} \end{equation*}

und berechnen nun die kompensierten (Evapo)transpirationsraten als

\begin{equation*} \begin{aligned} \hat{T}_c = T_u \min \left( \frac{1}{\omega}, \frac{1}{\omega_c} \right), \\ \hat{S}_c = S_u \min \left( \frac{1}{\omega}, \frac{1}{\omega_c} \right), \end{aligned} \end{equation*}

wobei \(\omega_c\) ein kritischer Stressindex ist.

kompensation.png

Abbildung 12: Effekt der Kompensation: Auf der linken Seite ist die räumliche Verteilung von \(\alpha\) in einer Bodensäule abgebildet. Auf der rechten Seite sind zwei Vergleichsrechnungen gemacht worden: (i) unkompensierte und (ii) kompensierte WWA. Die unkompensierte WWA holt eine reduzierte Menge Wasser aus den nicht gestressten Tiefen des Bodens. In der kompensierten Rechnung ist die WWA-Rate in den nicht gestressten Tiefen erhöht und kompensiert dadurch die reduzierte WWA aus den gestressten Tiefen. Rechnung und Plot: Andre Peters.

Hydraulische Umverteilung

Die hydraulische Umverteilung ist ein passiver Mechanismus, bei dem Wasser über unterirdische Netze von feuchten zu trockenen Böden transportiert wird. Zwar kann hydraulische Umverteilung zur Kompensation beitragen, ist jedoch nicht mit der oben diskutierten Kompensation zu verwechseln. Es handelt sich hier um die Umverteilung von Wasser in trockenere Bodenschichten und nicht um eine erhöhte WWA in feuchteren Schichten.

Der Mechanismus der hydraulischen Umverteilung basiert auf passivem Transport entlang des Gradienten im Wasserpotential. Wasser strömt hauptsächlich vom Boden in die Wurzel, da das Wasserpotential in der Wurzel niedriger ist als im Boden.

Unter verschiedenen Umständen kann es dazu kommen, das in den oberen Bodenschichten das Bodenwasserpotential niedriger ist als das Wurzelwasserpotential (Neumann et al., 2012)°. Zum Beispiel kann es in der Nacht, wenn die Stomata sich schliessen, steigt das Wasserpotential in der gesamten Pflanze, da keine Transpiration stattfindet. Wenn nun die obere Bodenschicht während des Tages besonders trocken gefallen ist, kann das Bodenwasserpotential das Pflanzenwasserpotential unterschreiten. In diesem Falle wird das Wasser aus den tieferen und feuchteren Bodenschichten in die Wurzel gezogen, tritt jedoch in den oberen trockenen Bodenschichten wieder aus. Der Transport ist deutlich schneller und kann größere Distanzen überwinden als kapillarer Aufstieg. Eine schnelle Umverteilung in tiefere Schichten mit demselben Mechanismus ist natürlich auch möglich.

Die hydraulische Umverteilung hat sehr wahrscheinlich großen Einfluss auf das Ökosystem—die genauen Auswirkungen werden jedoch immer noch untersucht. Einige Arbeitshypothesen nach Neumann et al. (2012)° zur Auswirkung sind: (i) höhere Transpirationsraten während Trockenphasen wenn Wasser nach oben umverteilt wird, wo sich der Großteil der Wurzeln befindet; (ii) Unterstützung von Sämlingen; (iii) erhöhte Aufnahme von Nährstoffen da die Bodenschicht befeuchtet wird, was auch die Mikroben im Boden aktiviert; (iv) verlängerte Lebensdauer von feinen Wurzeln; (v) Speicherung von Wasser im Boden durch Umverteilung nach unten.

Hydraulische Umverteilung kann durch eine Kopplung mit der Bodenhydraulik mit einem Pflanzenhydraulischen Modell simuliert werden (Bittelli et al., 2015)°.

System Boden–Pflanze–Atmosphäre

Soil–plant–atmosphere continuum (SPAC)

Das System Boden–Pflanze–Atmosphäre wird im englischen soil–plant–atmosphere continuum (SPAC) genannt. Der Begriff geht auf Philip (1966)° zurück, welcher ihn zum ersten mal einführte. Die grundsätzliche Theorie der SPAC Modellierung ist, dass das Wasserpotential in dem System Boden–Pflanze–Atmosphäre kontinuierlich abnimmt, was den Wasserfluss durch dieses System antreibt.

spac-pot.png

Abbildung 13: Wasserpotentiale im System Boden–Pflanze–Atmosphäre. Der Wasserfluss durch dieses System wir von der Differenz des Wasserpotentials zwischen den verschiedenen Kompartmenten angetrieben.

Das Wasserpotential ist die Energie die im Wasser gespeichert ist. Das Wasserpotential ist die Summe verschiedener Komponenten:

\begin{equation*} \psi = G + \Pi + P + \Psi \end{equation*}

Ausblick

Wer nun Interesse an der Modellierung von Wasser-, Energie und Stofftransport im System Boden–Pflanze–Atmosphäre bekommen hat ist herzlich zu unserem Mastermodul Ecohydrological Modelling eingeladen.

Danksagung

Dank an Dr. Mikael Gillefalk, der mir geholfen hat, diese Notizen für den Kurs WEST an der Technischen Universität Braunschweig zu entwickeln. Des weiteren auch Dank an Dr. Andre Peters, dessen Vorlesungsfolien ich benutzen durfte. Diese Notizen wurden im Rahmen des vom Land Niedersachsen finanzierten Juniorprofessurprogramm geschrieben.


Impressum

Autor: ilhan özgen xian

Created: 2024-02-05 Mon 09:18