====== Brooks-Corey Funktion ======

Die Brooks Corey Funktion stellt eine Beziehung zwischen Bodenfeuchte oder Sättigung einerseits und der Saugspannung andererseits her:

$$ \frac{\theta-\theta_r}{\theta_s-\theta_r} = {\left[\frac{\psi_0}{\psi}\right]}^\lambda$$

mit den Parametern $\theta$ Bodenfeuchte, $\theta_s$ Sättigungs-Bodenfeuchte, $\theta_r$ Rest-Bodenfeuchte, $\psi_0$ Lufteintrittspunkt, $\psi$ Saugspannung, Kapillarspannung.

Die ungesättigte Leitfähigkeit wird nach Brooks-Corey durch folgende Beziehung dargestellt:

$$ \frac{K({\theta})}{K_s} = {\left[{\frac{\theta-\theta_r}{\theta_s-\theta_r}}\right]}^{3+2*\lambda}$$

mit den Parametern $K_\theta$ ungesättigte Bodenfeuchte, $K_s$ gesättigte Bodenfeuchte, $\lambda$ Porengrößenindex.

<code S+ | Brooks-Corey.R>

# Brooks-Corey (1966)
brookscorey <- function(theta,tred,tmax,ksat,pi){
# Standard values
ths <- seq(0.05,0.4,0.005)
trs <- 0.035
tms <- 0.45                # maximale Feuchte
kss <- 1E-5                # hydraul. Leitfähigkeit in m/s
ips <- 0.2                 # pore size index
# Define values
if (is.null(theta)) {
        theta <- ths       # Default Bodenfeuchte
  }
if (is.null(tred))  {
        tred <- trs        # Default Restfeuchte
  }
if (is.null(tmax))  {
        tmax <- tms        # Default max. Feuchte
  }
if (is.null(ksat))  {
        ksat <- kss        # Default Durchlaessigkeit
  }
if (is.null(pi))    {
        pi <- ips          # Default Porenindex
  }
# Calculate unsaturated hydraulic conductivity
Se <- (theta-tred)/(tmax-tred)  
if ((Se < 0) || (Se > 1)) {
       result <- "Error in the definition of moisture."
       }
else{
  kunsat <- Se^(2+3*pi)
  if (is.null(kunsat)) {
    
                result <- "Error in calculation of unsaturated conductivity."
            }
  else{
    ratio  <- Se^(2+3*pi)
    kunsat <- ratio * ksat
    result <- list(Se,ratio,kunsat)}
  }
return(result)    
}

#Test
unsat <- brookscorey(0.25,0.05,0.4,1E-4,0.1)

# Values
ths <- seq(0.05,0.4,0.005)
trs <- 0.045
tms <- 0.45                # maximale Feuchte
kss <- 1E-5                # hydraul. Leitfähigkeit in m/s
ips <- 0.2                 # pore size index

unsat <- brookscorey(ths,trs,tms,kss,ips)
Se <- unsat[[1]]
ku <- unsat[[3]]

plot(Se,ku)

</code>