library(aqp)
library(soilDB)
library(dendextend)
library(cluster)
library(ape)
library(latticeExtra)
library(grid)
library(tactile)
library(viridis)
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)
}
)