# 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)
}```