Fitting van Genuchten Model using R

Fitting van Genuchten Model using R

The following functions can be used to efficiently fit van Genuchten (1981) water retention curve to paired water content and water potential data. The non-linear fitting routine requires the package minpack.lm.

# Function returns effective saturation, given water potential and the parameters alpha and n
vG = function(P,alpha,n){
 m =1-1/n
 vg= (1+(alpha*P)^n)^(-m)
 return(vg)
}
# Function returns water content, given water potential and the parameters alpha, n, theta_s, and theta_r
predict.vG = function(P,theta.r,theta.s,alpha,n){
 theta = vG(P,alpha,n)*(theta.s-theta.r)+theta.r
 return(theta)
}
# Fitting Routine ----------------------------------------------------
# The function returns best-fit parameters of the van Genuchten model, 
# given paired water content and water potential data. If any of the 
# of the four parameters are already known they can be supplied in the
# function call. The out put is a data frame with all the four parameters.
# If the parameters are fitted, the corresponding statistics are also provided
#
# Example usage: assuming n is known a priori
#
#  fit.data = fit.vG(theta.data,P.data, n=2.3)

fit.vG = function(theta,P,theta.s=-999,theta.r=-999,alpha=-999,n=-999){
 # REQUIRES
 # - minpack.lm | Library
 # - vG | Function
 
 # Initialize starting guess and bounds of free parameters
 s.par = c(theta.s=0.5,theta.r=0.01,alpha=0.1,n=2)
 l.par = c(theta.s=0.1,theta.r=0.,alpha=0.000001,n=1.05)
 u.par = c(theta.s=0.9,theta.r=0.3,alpha=10,n=8) 
 
 # If parameter is supplied, assume it is known and remove from free parameters
 if(theta.s!=-999){
 s.par = s.par[-1]
 l.par = l.par[-1]
 u.par = u.par[-1]
 }
 if(theta.r!=-999){
 s.par = s.par[-2]
 l.par = l.par[-2]
 u.par = u.par[-2]
 }
 if(alpha!=-999){
 s.par = s.par[-3]
 l.par = l.par[-3]
 u.par = u.par[-3]
 }
 if(n!=-999){
 s.par = s.par[-4]
 l.par = l.par[-4]
 u.par = u.par[-4]
 }
 
 # Setup Objective Function
 fit = nlsLM( theta ~ vG(P,alpha,n)*(theta.s-theta.r)+theta.r, start=s.par,lower=l.par, upper=u.par,trace=FALSE, control = nls.lm.control(maxiter=100))
 
 # Summarize fitted results. Output a list of:
 # 1. all parameters (so that predict function can use them all regardless of whether fitted or not)
 # 2. Statistics of fitted parameters
 param = coef(summary(fit))
 param = data.frame(t(param))
 
 if(theta.s!=-999){
 param$theta.s = c(theta.s,NA,NA,NA)
 }
 if(theta.r!=-999){
 param$theta.r = c(theta.r,NA,NA,NA)
 }
 if(alpha!=-999){
 param$alpha = c(alpha,NA,NA,NA)
 }
 if(n!=-999){
 param$n = c(n,NA,NA,NA)
 }
 # 3. Number of iterations and SSE
 
 
 
 
 return(param)
}

Comments are closed, but trackbacks and pingbacks are open.