Theoretical Water Retention Curves

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

Soil Texture Class

Align (order) soil texture classes based on (average) estimated water retention parameters.

data("ROSETTA.centroids")

# local copy
x <- ROSETTA.centroids

# cluster based on parameters
m <- ROSETTA.centroids[, c('sat', 'fc', 'pwp')]
row.names(m) <- x$texture
h <- as.hclust(agnes(daisy(m)))

# attemp to align with total AWC
h <- dendextend::rotate(h, order(x$awc))
h2 <- as.phylo(h)

# check: OK!
par(mar = c(0, 0, 0, 0))
plot(h2, label.offset = 0.005)

h$labels[h$order]

##  [1] "silt"            "silty clay loam" "silt loam"       "clay loam"       "silty clay"     
##  [6] "clay"            "sandy clay"      "loam"            "sandy clay loam" "sandy loam"     
## [11] "loamy sand"      "sand"

Explore several sorting strategies.

# clustering approach
x$texture_sort1 <- factor(x$texture, levels = rev(h$labels[h$order]))

# just AWC
x$texture_sort2 <- factor(x$texture, levels = x$texture[order(x$awc)])

Develop figures.

p.1 <- segplot(texture ~ pwp + fc, data=x, 
        horizontal = TRUE, 
        level = awc, col.regions=viridis,
        main='Available Water Holding Capacity\nUSDA-ARS ROSETTA Model Centroids', 
        sub='Sorted According to Field Book v3.0',
        xlab=expression(Volumetric~Water~Content~~(cm^3/cm^3)),
        par.settings = tactile.theme()
)

p.2 <- segplot(texture_sort1 ~ pwp + fc, data=x, 
        horizontal = TRUE, level = awc, col.regions=viridis,
        main='Available Water Holding Capacity\nUSDA-ARS ROSETTA Model Centroids',
        sub='Sorted According to {PWP, FC, SAT} Values',
        xlab=expression(Volumetric~Water~Content~~(cm^3/cm^3)),
        par.settings = tactile.theme()
)

p.3 <- segplot(texture_sort2 ~ pwp + fc, data=x, horizontal = TRUE, 
        level = awc, col.regions=viridis,
        main='Available Water Holding Capacity\nUSDA-ARS ROSETTA Model Centroids',
        sub='Sorted According to AWC',
        xlab=expression(Volumetric~Water~Content~~(cm^3/cm^3)),
        par.settings = tactile.theme()
        )

---

This document is based on aqp version 2.0.