# Wrapper function for the Estimation of Soil Hydrologic Property Model Parameters Estimates model parameters of implemented soil hydraulic property functions. This function sets up the parameter estimation, given a set of arguments, and enables minimisation of (weighted) sum of squared residuals, assuming independent and identically distributed model residuals. More information on the options is given in the **Details** ```r shypEstFun( shpmodel = "01110", parL, retdata, condata, ivap = NULL, hclip = FALSE, weightmethod = "none", LikModel = "rss", ALG = "DE", set.itermax = 200, ALGoptions = NULL, lhs.query = FALSE ) ``` ## Arguments - `shpmodel`: |||| |:--|:--|:--| |Character specifying the soil hydraulic property model. Currently valid models as documented in shypFun and are:| |`01110` constrained unimodal van Genuchten-Mualem.| |`01210` constrained bimodal van Genuchten-Mualem.| |`01310` constrained trimodal van Genuchten-Mualem.| |`02110` unimodel Kosugi 2 parametric-Mualem model (Kosugi, 1996)| |`03110` unimodel Fredlund-Xing-Mualem model, with the contraint of m = 1-1/n (Fredlund D.G., and A. Xing, 1994)| - `parL`: `list` of 4 vectors with named vectors, the order in the list is sensitive. |||| |:--|:--|:--| |`p`|vector with length `l` of model specific initial parameters, has to coincide with the chosen soil hydraulic property model.| |`psel`|vector with length `l` identifying which parameters are to be estimated| |`plo`|vector of lower bounds (non-transformed parameter boundaries)| |`pup`|vector of upper bounds (non-transformed parameters boundaries)| - `retdata`: A dataframe or matrix with 2 columns. The first with log10 values of pressure head values in [cm] and the second with volumetric water contents in [cm cm-3]. - `condata`: A dataframe or matrix with 2 columns. The first with log10 values of pressure head values in [cm] and the second with hydraulic conductivity values log10[cm d-1]. - `ivap`: Specification if isothermal vapour conductivity after Saito et al. (2006) is accounted, defaults to `NULL` and no isothermal vapour conducitvity is considered. See Details. - `hclip`: Implemented for future development reasons and is not yet functional. Specification if the hydraulic conductivity model should be 'clipped', i.e. constrained to a maxium pore diamater as introduced by Iden et al. (2015), defaults to `FALSE`. - `weightmethod`: Specification of weight method. The implemented methods are |||| |:--|:--|:--| |`none`|no weights are considered, i.e. no measurement error assumed| |`range`|normalization of observations to the intervall [0,1]| |`fix1`|fixed scalar weight for THETA is 0.05^2 and weight for log10K is 1| |`est1`|Two scalar model weights (\$1/sigma^2) are treated as free parameters to be estimated by inversion, one for `THETA` and one for `log10K`| Alternatively, a list of vectors can be provided specifying the user given model weights (\$1/sigma^2). Either as skalar for each data class, or a vector with the same length as the number of data points given for each of the measurements in the respective data class. The length of the list has to coincide with the data groups. - `LikModel`: Specification of inverse modelling type. Has to be specified but implemented for future compatability) |||| |:--|:--|:--| |`rss`|Default for the optimisation algorithm `DEoptim`. `resFun` returns skalar sum of squared (weighted) residuals| |`-2loglik`|Specified if `ALG == -2*log-likelihood` value, which is minimised assuming Gaussian and i.i.d (weighted) residuals| - `ALG`: Select global optimisation algorithm or a Markov chain Monte Carlos (MCMC) sampler. |||| |:--|:--|:--| |`DE`|Default for the optimisation algortihm DEoptim. resFun returns a skalar sum of squared (weighted) residuals if `LikModel == "rss"`.| |`modMCMC`|Default for the DRAM (Delayed Rejection Adaption Metropolis) algrothim implemented in `modMCMC` of the FME package. resFun returns a skalar -2loglik and `LikModel = "-2logLik"` has to be specified.| - `set.itermax`: Integer specifying the maximum number of iterations `default = 200`. - `ALGoptions`: A list with named entries setting the algorithm options. Each list element name is required to be identical with the names as documented in the respective algortihm help DEoptim.control and modMCMC . `set.itermax` overrides the maximum iterations argument. - `lhs.query`: default `FALSE`, `TRUE` will produce a Latin Hypercube Sample for the initial population when using `DEoptim`. ## Returns `list` returns the result of the optimisation algrorithm or MCMC sampler and all settings. - **settings**: a `list` with output of the optimisation and summary of settings: |||| |:--|:--|:--| |`weigth`|the `list` with weights for the retention and conductivity data.| |`parL`|the `list` of initial and selected model parameters, and upper and lower bounds.| |`transL`|list of parameter transformation rules used| |`shpmodel`|the used soil hydraulic property model| |`ivap`|isothermal vapour conductivity model| |`hclip`|for future compatability| |`LikModel`|the adopted method to calculate the objective function value| |`data`|a `list` of 2 objects with a) retention data and b) conductivity data used for the parameter estimation.| - **out**: result of algorithm function `DEoptim` or `modMCMC` ## Details Several in-built methods for weighting the (multi-) objective function residuals are available, they may be specified, or estimated as nuisance parameters for the two data groups. More details see `weightFun`. Weights are the inverse of the squared standard deviation of the residuals (variance). Generally, soil hydraulic property model parameters are estimated as transformed parameters: log10 for alpha_i, Ks, and log10 for n_i-1, Kc, Knc For model codes in **ivap** please refer to KvapFun . Parallel computing for package `DEoptim` is not supported. And the optional arguments in `modMCMC` are not supported. ## Examples ```r ## Not run: data("shpdata1") retdata <- shpdata1$TS1$wrc condata <- shpdata1$TS1$hcc condata <- condata[!is.na(condata[,1]),] weightmethod <- "range" ivap <- NULL set.itermax <- 1 LikModel <- "rss" # ALTERNATIVE OPTION: LikModel = "-2logLik" ALG <- "DE" # ALTERNATIVE OPTION: ALG = "modMCMC" parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ks"=100,"tau"=.5), "psel" = c(1, 1, 1, 1, 1, 1), "plo"= c(0.001 , 0.2 , 0.001 , 1.1, 1, -2), "pup"= c(0.3 , 0.8 , .1, 11 , 1e4, 10)) out <- shypEstFun(shpmodel = "01110", parL = parL, retdata = retdata, condata = condata, ivap = ivap, hclip = FALSE, weightmethod = weightmethod, LikModel = LikModel, ALG = ALG, set.itermax = set.itermax, lhs.query = FALSE) \dontshow{ } \donttest{ data("shpdata1") retdata <- ret <- shpdata1$TS1$wrc condata <- con <- shpdata1$TS1$hcc condata <- condata[!is.na(condata[,1]),] --- # 1 SET VARIABLES -------------------- # VARIABLES FOR PLOTTING {pF <- seq(-3, 6.8, length = 201) h <- 10^pF ticksatmin <- -2 tcllen <- 0.4 ticksat <- seq(ticksatmin,5,1) pow <- ticksatmin:6 # VARIABLES FOR THE FITTING ALGORITHM weightmethod = "range" ivap = NULL set.itermax = 3e1 LikModel = "rss" # ALTERNATIVE OPTION: LikModel = "-2logLik" ALG = "DE" # ALTERNATIVE OPTION: ALG = "modMCMC" shpmodel.v <- c("01110", "01110FM") plot.query = FALSE no.shps <- length(shpmodel.v) # initialising lists out.L <- vector("list", no.shps) gof.L <- vector("list", no.shps) } # Run comparison for (i in 1:2) { shpmodel = shpmodel.v[i] # INITIAL PARAMETERS, BOUNDS, and SELECTED PARAMETERS FOR FITTING switch(shpmodel, "01110" = { # van Genuchten-Mualem Model parameters parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ks"=100,"tau"=.5), "psel" = c(1, 1, 1, 1, 1, 1), "plo"= c(0.001 , 0.2 , 0.001 , 1.1, 1, -2), "pup"= c(0.3 , 0.8 , .1, 11 , 1e4, 10) ) }, "01110FM" = { # van Genuchten-Mualem Model parameters + BRUNSWICK MODEL parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ksc"=100, "tau"=.5,"Ksnc"=1e-4,"a"=1.5,"h0"=6.8), "psel" = c( 1,1, 1 ,1 , 1,1,1, 0, 0), "plo"= c(0.001 , 0.1 , 0.001 , 1.1, 1,0,1e-6 , 1, 6.5), "pup"= c(0.35, 0.7 , .1, 11 , 1e4,10 ,1e0, 2, 6.9) ) }, stop("Enter a meaningful shpmodel") ) out <- shypEstFun(shpmodel = shpmodel, parL = parL, retdata = retdata, condata = condata, ivap = ivap, hclip = FALSE, weightmethod = weightmethod, LikModel = LikModel, ALG = ALG, set.itermax = set.itermax ,lhs.query = FALSE) out$model <- shpmodel.v[i] out.L[[i]] <- out # Calculate the soil hydraulic properties for the plot if(ALG == "DE"){ p <- out$out$optim$phattrans } if(ALG == "modMCMC"){ p <- out$out$bestpartrans } if(weightmethod =="est1"){ np <- length(p) p <- p[-c(np-1, np)] if(ALG =="modMCMC"){ parL$p[which(parL$psel==1)] <- p p <- parL$p } } if(plot.query==TRUE){ shyp.L<-shypFun(p,h,shpmodel=shpmodel.v[i],ivap.query=ivap) if(shpmodel == c("01110")){ wrc<-shyp.L$theta hcc<-log10(shyp.L$Kh) # PLOT THE WATER RETENTION CURVE par(mfrow=c(1,2),tcl=tcllen) plot(retdata,ylim=c(.0,.50),xlim=c(0,6.8),ylab="",xlab="", col="darkgrey",axes=FALSE,main="WaterRetentionCurve",cex=2) lines(log10(abs(h)),wrc,col="darkblue",lwd=2) legend("topright",c("observed","simulated"),pch=c(1,NA), lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue")) axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen) axis(2,tcl=tcllen) axis(4,labels=NA) axis(3,labels=NA) mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8) mtext("vol.watercontent[-]",2,cex=1.2,line=2.8) box() # PLOT THE MEASURED HYDRAULIC CONDUCTIVITY CURVE plot(condata,ylim=c(-8,2),xlim=c(0,6.8),ylab="",xlab="",col="darkgrey", axes=FALSE,main="HydraulicConductivityCurve",cex=2) lines(log10(abs(h)),hcc,col="darkblue",lwd=2) legend("topright",c("observed","simulated"),pch=c(1,NA), lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue")) axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen) axis(2) axis(4,labels=NA) axis(3,labels=NA) mtext("log10K[cm/d]",2,cex=1.2,line=2.8) mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8) box() par(mfrow=c(1,1)) } if(shpmodel == "01110FM"){ wrc<-shyp.L$theta wrccap<-shyp.L$thetacap wrcnc<-shyp.L$thetanc hcc<-log10(shyp.L$Kh) hcccap<-log10(shyp.L$Kcap) hccnc<-log10(shyp.L$Knc) hcvap<-log10(shyp.L$Kvap) par(mfrow=c(1,2),tcl=tcllen) plot(retdata,ylim=c(.0,.50),xlim=c(0,6.8),ylab="",xlab="", col="darkgrey",axes=FALSE,main="WaterRetentionCurve",cex=2) lines(log10(h),wrccap,col="brown",lwd=2) lines(log10(h),wrcnc,col="brown",lwd=2,lty=2) lines(log10(h),wrc,col="darkblue",lwd=2) legend("topright",c("observed","simulated"),pch=c(1,NA), lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue")) axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen) axis(2,tcl=tcllen) axis(4,labels=NA) axis(3,labels=NA) mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8) mtext("vol.watercontent[-]",2,cex=1.2,line=2.8) box() # PLOT THE HYDRAULIC CONDUCTIVITY CURVE plot(condata,ylim=c(-8,max(max(condata[,1]),max(hcc))) ,xlim=c(0,6.8),ylab="",xlab="",col="darkgrey", axes=FALSE,main="HydraulicConductivityCurve",cex=2) lines(log10(h),hcccap,col="brown",lwd=2) lines(log10(h),hccnc,col="brown",lwd=2,lty=2) lines(log10(h),hcc,col="darkblue",lwd=2) lines(log10(h),hcvap,col="darkblue",lwd=2) legend("topright",c("observed","simulated"),pch=c(1,NA), lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue")) axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen) axis(2) axis(4,labels=NA) axis(3,labels=NA) mtext("log10K[cm/d]",2,cex=1.2,line=2.8) mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8) box() par(mfrow=c(1,1)) } } phattrans.m <- out$out$optim$phattrans gof.L[[i]]<-gofFun(phattrans.m,shpmodel=shpmodel.v[i],retdata=retdata,condata=condata, out.L[[i]]$settings$weight,parL$psel,ivap.query=NULL,hclip.query=FALSE) } statstab3 <- cbind("th_rmse" = unlist(lapply(lapply(gof.L, `[[`, "th"), '[[', "rmse")), "log10Kh_rmse" = unlist(lapply(lapply(gof.L, `[[`, "log10Kh"), '[[', "rmse")) ) } ## End(Not run) ```

shypEstFun function