1 Introduction

This is a simple demonstration of how to use (estimated) van Genuchten Model parameters from the KSSL snapshot (added 2016-11-17) as provided by fetchKSSL() to generate water retention curves. These parameters are based on a model fit using the Rosetta software. See the Rosetta Texture Class Averages for additional documentation.

A note from our KSSL staff:

The Rosetta computer program must have at a minimum the particle size distribution. Then it will use the following properties, if provided, bulk density, water content at 33 kPa, and water content at 1500 kPa.

As a practice, the KSSL uses the Rosetta computer program for a horizon only if all of these properties are present. For PSD, we only use the routine method, like pipette. For bulk density, we only use the clod at 33 kPa. For water content at 33 kPa, we only use the clod. For water content at 1500 kPa, we only use the <2 mm sieved sample.

1.1 Setup

With a recent version of R (>= 3.5), it is possible to get all of the packages that this tutorial depends on via:

# run these commands in the R console
install.packages('latticeExtra', dep=TRUE)
install.packages('reshape2', dep=TRUE)
install.packages('aqp', dep=TRUE)
install.packages('soilDB', dep=TRUE)

You will also need the latest version of soilDB and aqp:

remotes::install_github("ncss-tech/aqp", dependencies=FALSE, upgrade_dependencies=FALSE)
remotes::install_github("ncss-tech/soilDB", dependencies=FALSE, upgrade_dependencies=FALSE)

1.2 A Very Basic Example

Here is a very quick example of how to build a water retention curve from KSSL data. Lets get the van Genuchten parameters for the top horizon of pedon ID S08NV003003. See ?KSSL_VG_model for details on how this function works.

# load required libraries
library(aqp)
library(soilDB)
library(latticeExtra)
library(reshape2)
library(grid)

# you can get the data like this:
# s <- fetchKSSL(pedon_id = 'S08NV003003')

# or, we can manualy enter it in like this
vg.params <- data.frame(theta_r=0.0337, theta_s=0.4864, alpha=-1.5814, npar=0.1227)
vg.model <- KSSL_VG_model(vg.params)

p.model <- xyplot(phi ~ theta, data=vg.model$VG_curve, type=c('l', 'g'), scales=list(alternating=3, x=list(tick.number=10), y=list(log=10, tick.number=10)), yscale.components=yscale.components.logpower, ylab=expression(Matric~~Potential~~(kPa)), xlab=expression("Volumetric Water Content" (cm^{3}/cm^{3})), par.settings = list(plot.line=list(col='RoyalBlue', lwd=2)))

update(p.model, main='Estimated Water Retention Curve\nS08NV003003\n0-9cm', sub='van Genuchten Model Parameters fit by USDA-ARS Rosetta')

2 Generate Curves for the Cecil Soil Series

Continuing from above, this time plotting water retention curves for all horizons associated with the Cecil soil series. Horizon designations are generalized for grouping of similar data.

# get KSSL data
s <- fetchKSSL(series = 'cecil')

# subset VG parameters
s.vg.hz <- horizons(s)[, c('labsampnum', 'hzn_desgn', 'hzn_top', 'hzn_bot', 'lab_texture_class', 'theta_r', 'theta_s', 'alpha', 'npar')]
# remove rows with missing value
s.vg.hz <- na.omit(s.vg.hz)

# simple generalization of horizons
s.vg.hz$hz <- generalize.hz(s.vg.hz$hzn_desgn, new=c('A', 'Bt', 'C'), pat=c('A', 'Bt', '^C'))

# melt into long format, for viz. of VG parameters
s.long <- melt(s.vg.hz, id.vars = 'hz', measure.vars = c('theta_r', 'theta_s', 'alpha', 'npar'))
# bw plot of VG parameters by generalized horizon label
bwplot(hz ~ value | variable, data=s.long, scales = list(x=list(relation='free', alternating=3)), strip=strip.custom(bg=grey(0.85)), as.table=TRUE, xlab='', main='van Genuchten Model Parameters fit by USDA-ARS Rosetta', sub='Pedons Correlated to the Cecil Soil Series')

# iterate over horizons and generate VG model curve
res <- lapply(1:nrow(s.vg.hz), function(i) {
  m <- KSSL_VG_model(VG_params = s.vg.hz[i, ], phi_min = 10^-3, phi_max=10^6)$VG_curve
  # copy generalized hz label
  m$hz <- s.vg.hz$hz[i]
  # copy ID
  m$.id <- s.vg.hz$labsampnum[i]
  return(m)
})

# copy over lab sample number as ID
res <- do.call('rbind', res)

# plot each curve, panel by genhz
p.model <- xyplot(phi ~ theta | hz, groups=.id, data=res, type=c('l', 'g'), scales=list(alternating=3, x=list(tick.number=10), y=list(log=10, tick.number=10)), yscale.components=yscale.components.logpower, ylab=expression(Suction~~(kPa)), xlab=expression(Volumetric~Water~Content~~(cm^3/cm^3)), par.settings = list(superpose.line=list(col='RoyalBlue', lwd=1)), strip=strip.custom(bg=grey(0.85)), as.table=TRUE)

update(p.model, main='Estimated Water Retention Curves', sub='van Genuchten Model Parameters fit by USDA-ARS Rosetta')

3 Additional Examples

See the SCAN/SNOTEL tutorial for more ideas.

This document is based on aqp version 1.43 and soilDB version 2.7.2.