#-------------------------------------------------------
#  Tiefenverlauf von Wassergehalt, Matrixpotential 
#                und Gravitationspotential
#  Annahme: hydraulisches Gleichgewicht (q(z) = 0)
#-------------------------------------------------------
# Aufgaben: a) Variieren Sie die Lage des Bezugsniveaus  
#           b) Verändern Sie die Richtung des z-Pfeils
#           c) Variieren Sie die Parameter des MVG-Modells
#
# --- Parametrisierung der Wasserspannungskurve ---
#
theta_VG <- function(alpha,n,theta_s,theta_r,psi) {
    m <- 1-1/n
    theta_VG <- theta_r + (theta_s-theta_r)/(1+(alpha*psi)^n)^m
}
#
# --- Mualemmodell für die relative hydraulische Leitfähigkeit ---
#
K_VGrel <- function(alpha,n,theta_s,theta_r,psi) {
    m <- 1-1/n
    K_VGrel <- (1-(alpha*psi)^(n-1)*(1+(alpha*psi)^n)^(-m))^2/
                 (1+(alpha*psi)^n)^(m/2)
}
#
# --- Modellparameter ---       
#
theta_s <- 0.55   # Sättigungswassergehalt
theta_r <- 0.15   # Residualwassergehalt
K_s <- 100        # Wasserleitfähigkeit in cm/d
alpha <- 0.1      # emp. Formparameter
n <- 2.5          # emp. Formparameter
#
gwl <- 110        # Tiefe des Grundwasserspiegels in cm
bzn <-  0         # Tiefe des Bezugsniveaus
iz  <- -1         # Indikator für Richtung der z-Achse: 
                  # +1 nach oben, -1 nach unten
par(mfrow=c(2,3))
# 
# Profile von Matrixpotential, Gravitationspotential, Gesamtpotential,
# Wassergehalt und hydraulischer Leitfähigkeit
#
z   <- -iz*seq(0,gwl,0.5) # Tiefe
gwl <- -iz*gwl
bzn <- -iz*bzn                                             
psi <- -iz*(z - gwl)
grp <- -iz*(bzn - z)
H   <-  psi + grp
theta <- theta_VG(alpha,n,theta_s,theta_r,abs(psi))
K_rel <- K_VGrel(alpha,n,theta_s,theta_r,abs(psi))
K <- K_s*K_rel
#
# Graphische Umsetzung
#
if(iz==1) {z.u <- min(z); z.o <- max(z)} else {z.u <- max(z); z.o <- min(z)}
plot(psi,z,     ylim=c(z.u,z.o),type="l",lwd=3,col="blue",main="Matrixpotential")
plot(grp,z,     ylim=c(z.u,z.o),type="l",lwd=3,col="orange",main="Gravitationspotential")
plot(H,z,       ylim=c(z.u,z.o),type="l",lwd=3,col="black",main="Gesamtpotential")
plot(theta,z,   ylim=c(z.u,z.o),type="l",lwd=3,col="blue",main="Vol. Wassergehalt")
plot(K,z,       ylim=c(z.u,z.o),type="l",lwd=3,col="blue",main="Leitfähigkeit")
plot(log10(K),z,ylim=c(z.u,z.o),type="l",lwd=3,col="blue",main="Log10(K)")
stop("++++ hier ++++")
# Skriptende