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.