Setup

library(aqp)
library(soilDB)

library(dendextend)
library(cluster)
library(ape)

library(latticeExtra)
library(grid)
library(tactile)
library(viridis)

Theoretical Water Retention Curves

Average soil hydraulic parameters by USDA soil texture class, extracted from the ROSETTA documentation. Approximate water-retention data have been computed for an idealized volume of homogenous soil material via van Genuchten model.

texture theta_r theta_s alpha npar theta_r_sd theta_s_sd alpha_sd npar_sd sat fc pwp awc
clay 0.098 0.459 -1.825 0.098 0.107 0.079 0.68 0.07 0.459 0.332 0.189 0.143
clay loam 0.079 0.442 -1.801 0.151 0.076 0.079 0.69 0.12 0.442 0.255 0.116 0.139
loam 0.061 0.399 -1.954 0.168 0.073 0.098 0.73 0.13 0.399 0.235 0.091 0.144
loamy sand 0.049 0.390 -1.459 0.242 0.042 0.070 0.47 0.16 0.390 0.103 0.052 0.051
sand 0.053 0.375 -1.453 0.502 0.029 0.055 0.25 0.18 0.375 0.054 0.053 0.001
sandy clay 0.117 0.385 -1.476 0.082 0.114 0.046 0.57 0.06 0.385 0.278 0.190 0.087
sandy clay loam 0.063 0.384 -1.676 0.124 0.078 0.061 0.71 0.12 0.384 0.228 0.111 0.117
sandy loam 0.039 0.387 -1.574 0.161 0.054 0.085 0.56 0.11 0.387 0.167 0.062 0.105
silt 0.050 0.489 -2.182 0.225 0.041 0.078 0.30 0.13 0.489 0.283 0.069 0.214
silty clay 0.111 0.481 -1.790 0.121 0.119 0.080 0.64 0.10 0.481 0.320 0.174 0.146
silty clay loam 0.090 0.482 -2.076 0.182 0.082 0.086 0.59 0.13 0.482 0.304 0.121 0.183
silt loam 0.065 0.439 -2.296 0.221 0.073 0.093 0.57 0.14 0.439 0.294 0.086 0.208

Develop theoretical water retention curves via van Genuchten model and texture class centroids.

# load average soil hydraulic parameters
data("ROSETTA.centroids")

# re-level texture class by approximate AWC
ll <- ROSETTA.centroids$texture[order(ROSETTA.centroids$awc)]
ROSETTA.centroids$texture <- factor(ROSETTA.centroids$texture, levels=ll)

# iterate over horizons and generate VG model curve
res <- lapply(1:nrow(ROSETTA.centroids), function(i) {
  # fit model using parameters from centroids
  # model bounds are given in kPA of suction
  vg <- KSSL_VG_model(VG_params = ROSETTA.centroids[i, ], phi_min = 10^-3, phi_max=10^6)
  
  # extract curve and add texture ID
  m <- vg$VG_curve
  m$texture <- ROSETTA.centroids$texture[i]
  
  return(m)
})

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

# check: OK
head(res)
##           phi     theta texture
## 1 0.001000000 0.4589988    clay
## 2 0.001232847 0.4589984    clay
## 3 0.001519911 0.4589980    clay
## 4 0.001873817 0.4589974    clay
## 5 0.002310130 0.4589966    clay
## 6 0.002848036 0.4589955    clay
tps <- tactile.theme(superpose.line=list(col=viridis::viridis(12), lwd=2))

xyplot(
  phi ~ theta, groups=texture, data=res, 
  type = 'l', 
  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 = tps, 
  strip = strip.custom(bg = grey(0.85)), 
  auto.key = list(columns = 4, lines = TRUE, points = FALSE, cex = 0.8),
  panel = function(...) {
    .y <- log(c(0.1, 33, 1500), base = 10)
    panel.grid(-1, -1)
    panel.abline(h = .y, lty = 3, col = 1)
    panel.text(x = c(0.025), y = .y, label = c('SAT', 'FC', 'PWP'), pos = 1, font = 2, cex = 1)
    panel.superpose(...)
  }
)

Plot water retention curves and annotate with permanent wilting point (approximately 1500kPa suction), field capacity (approximately 33kPa suction), and saturation (close to 0 kPa suction). Note that matric potential (typically negative pressure by convention) are displayed as suction (positive values).

trellis.par.set(superpose.symbol=list(col=c('black', 'black', 'black'), fill=c('firebrick', 'orange', 'royalblue'), cex=1, pch=22))

sk <- simpleKey(text=c('permanent wilting point', 'field capacity', 'saturation'), points = TRUE, columns=3)

xyplot(
  phi ~ theta | texture, data=res, 
  type=c('l', 'g'), 
  scales=list(alternating=3, x=list(tick.number=6), y=list(log=10, tick.number=6)), 
  yscale.components=yscale.components.logpower, 
  ylab=expression(Matric~~Potential~~(-kPa)), 
  xlab=expression(Volumetric~Water~Content~~(cm^3/cm^3)), 
  par.settings = tactile.theme(plot.line=list(col='black', lwd=2)), 
  strip=strip.custom(bg=grey(0.85)), 
  as.table=TRUE,
  key=sk,
  main='Idealized Water Retention\nUSDA-ARS ROSETTA Model Centroids', 
  # sub='van Genuchten model',
  subscripts = TRUE,
  panel= function(x, y, subscripts=subscripts, ...) {
    # idealized water retention curve
    panel.xyplot(x, y, ...)
    
    # get texture (i.e. the current panel)
    tx <- res$texture[subscripts][1]
    # look-up centroid data for this texture
    centroid.data <- ROSETTA.centroids[match(tx, ROSETTA.centroids$texture), ]
    
    # add PWP, note that y-scale is log_10-transformed
    panel.points(x=centroid.data$pwp, y=log(1500, base = 10), col='black', fill='firebrick', cex=0.85, pch=22)
    # add PWP, note that y-scale is log_10-transformed
    panel.points(x=centroid.data$fc, y=log(33, base = 10), col='black', fill='orange', cex=0.85, pch=22)
    # add SAT, note that y-scale is log_10-transformed
    panel.points(x=centroid.data$sat, y=log(0.1, base = 10), col='black', fill='royalblue', cex=0.85, pch=22)
    
    # annotate with approx AWC
    grid.text(sprintf('AWC: %s', round(centroid.data$awc, 3)), x = unit(0.95, 'npc'), y = unit(0.95, 'npc'), gp=gpar(cex=0.75), hjust=1)
    
  }
)