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.4.

Theoretical Water Retention Curves

D.E. Beaudette

2024-07-30

Setup

Theoretical Water Retention Curves

Soil Texture Class