Introduction

This document outlines how to access a new, online, interface to the ROSETTA model for predicting hydraulic parameters from soil properties. The ROSETTA API was developed by Dr. Todd Skaggs (USDA-ARS) and links to the work of Zhang and Schaap, (2017). The ROSETTA function in soilDB is one possible approach to using the API from R.

The ROSETTA model relies on a minimum of 3 soil properties, with increasing (expected) accuracy as additional properties are included:

required, sand, silt, clay: USDA soil texture separates (percentages) that sum to 100%
optional, bulk density at 33 kPa (1/3 bar): mass per volume after accounting for >2mm fragments, units of gram/cm³
optional, volumetric water content at 33 kPa (1/3 bar): roughly “field capacity” for most soils, units of cm³/cm³
optional, volumetric water content at 1500 kPa (15 bar): roughly “permanent wilting point” for most plants, units of cm³/cm³

Soil properties must be described, in order, via vars argument. The API does not use the names but column ordering must follow: sand, silt, clay, bulk density, volumetric water content at 33kPa (1/3 bar), and volumetric water content at 1500kPa (15 bar).

Consider using the interactive version, with copy/paste functionality.

Model Versions

Version 1: Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology 251(3-4): 163-176. doi: 10.1016/S0022-1694(01)00466-8.
Version 2: Schaap, M.G., A. Nemes, and M.T. van Genuchten. 2004. Comparison of Models for Indirect Estimation of Water Retention and Available Water in Surface Soils. Vadose Zone Journal 3(4): 1455-1463. doi: https://doi.org/10.2136/vzj2004.1455.
Version 3: Zhang, Y., and M.G. Schaap. 2017. Weighted recalibration of the Rosetta pedotransfer model with improved estimates of hydraulic parameter distributions and summary statistics (Rosetta3). Journal of Hydrology 547: 39-53. doi: 10.1016/j.jhydrol.2017.01.004.

Setup

Install all of the packages that this tutorial depends on via:

# run these commands in the R console
install.packages('httr', dependencies =  TRUE)
install.packages('jsonlite', dependencies =  TRUE)
install.packages('latticeExtra', dependencies = TRUE)
install.packages('tactile', dependencies =  TRUE)
install.packages('soilDB', dependencies = TRUE)
install.packages('aqp', dependencies = TRUE)

A Very Basic Example

The following demonstrates how to push soil property data from SDA through the ROSETTA API.

# required libraries
library(soilDB)

# SQL to be submitted to SDA
q <- "SELECT
-- contextual data
compname, comppct_r,
-- horizon morphology
hzname, hzdept_r, hzdepb_r, 
-- parameters used by ROSETTA, may contain NULL
sandtotal_r, silttotal_r, claytotal_r, dbthirdbar_r, 
wthirdbar_r/100 AS wthirdbar_decimal, wfifteenbar_r/100 AS wfifteenbar_decimal
-- tables of interest
FROM 
mapunit AS mu
-- implicit filtering 
INNER JOIN component AS co ON mu.mukey = co.mukey
INNER JOIN chorizon AS ch ON co.cokey = ch.cokey
-- single map unit
WHERE mu.mukey = '373596'
-- order for visual inspection
ORDER BY co.cokey, ch.hzdept_r ASC;"

# submit query to SDA
x <- SDA_query(q)

# attempting to use all possible soil properties
vars <- c('sandtotal_r', 'silttotal_r', 'claytotal_r', 'dbthirdbar_r', 'wthirdbar_decimal', 'wfifteenbar_decimal')

# call ROSETTA API 
r <- ROSETTA(x, vars = vars)

compname	comppct_r	hzname	hzdept_r	hzdepb_r	sandtotal_r	silttotal_r	claytotal_r	dbthirdbar_r	wthirdbar_decimal	wfifteenbar_decimal	theta_r	theta_s	alpha	npar	ksat	.rosetta.model	.rosetta.version
Aluf	50	H1	0	117	94.4	0.6	5.0	1.60	0.117	0.049	0.027	0.371	-1.281	0.160	2.285	5	1
Aluf	50	H2	117	251	61.0	19.0	20.0	1.55	0.224	0.134	0.055	0.385	-1.484	0.119	1.465	5	1
Hitilo	30	H1	0	137	94.9	0.6	4.5	1.55	0.108	0.040	0.025	0.384	-1.255	0.163	2.428	5	1
Hitilo	30	H2	137	178	55.4	17.6	27.0	1.58	0.252	0.181	0.080	0.388	-1.387	0.111	1.324	5	1
Hitilo	30	H3	178	203	60.2	18.3	21.5	1.50	0.221	0.135	0.060	0.400	-1.435	0.124	1.629	5	1

Theoretical Water Retention Curves

Develop theoretical water retention curves (van Genuchten model) from soil texture class centroids.

# required libraries
library(aqp)
library(latticeExtra)
library(tactile)

# make a table of soil textures
# sand sand/silt/clay values at geometric centroids
tex <- SoilTextureLevels(simplify = FALSE)
x <- texcl_to_ssc(tex)
x <- cbind(x, tex)

sand	silt	clay	tex
92	5	3	cos
92	5	3	s
92	5	3	fs
92	5	3	vfs
82	12	6	lcos
82	12	6	ls

Submit example data to the ROSETTA API and generate theoretical water retention curves.

r <- ROSETTA(x, vars = c('sand', 'silt', 'clay'))

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

# flatten to data.frame
res <- do.call('rbind', res)

phi	theta	texture
0.001	0.379	cos
0.001	0.379	cos
0.002	0.379	cos
0.002	0.379	cos
0.002	0.379	cos
0.003	0.379	cos

Plot water retention curves. Note that matric potential (typically negative pressure by convention) are displayed as suction (positive values).

tps <- tactile.theme(
  plot.line=list(col='royalblue', lwd=2)
  )

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 = tps, 
  strip = strip.custom(bg=grey(0.85)), 
  as.table = TRUE,
  layout = c(7, 3),
  main='Idealized Water Retention'
)

Automatic Model Selection

# data from Texas
q <- "SELECT TOP 10 musym, co.cokey, compname, comppct_r,
    hzname, hzdept_r, hzdepb_r, sandtotal_r, silttotal_r, claytotal_r, dbthirdbar_r,
    wthirdbar_r/100 AS wthirdbar_decimal, wfifteenbar_r/100 AS wfifteenbar_decimal
    FROM legend
    INNER JOIN mapunit mu ON mu.lkey = legend.lkey
    INNER JOIN component co ON mu.mukey = co.mukey
    INNER JOIN chorizon ch ON co.cokey = ch.cokey
    WHERE legend.areasymbol LIKE 'TX%'
    ORDER BY musym, co.cokey, ch.hzdept_r ASC;"

x <- SDA_query(q)

# simulate missing data
x$dbthirdbar_r[1] <- NA
x$wthirdbar_decimal[2] <- NA
x$wfifteenbar_decimal[3] <- NA
x$sandtotal_r[9] <- NA
x[10, ] <- NA

# attempting to use all possible soil properties
vars <- c('sandtotal_r', 'silttotal_r', 'claytotal_r', 'dbthirdbar_r', 'wthirdbar_decimal', 'wfifteenbar_decimal')

# automatic model selection by API
# based on available data
# m = "0"
r <- ROSETTA(x, vars = vars)

musym	cokey	compname	comppct_r	hzname	hzdept_r	hzdepb_r	sandtotal_r	silttotal_r	claytotal_r	dbthirdbar_r	wthirdbar_decimal	wfifteenbar_decimal	theta_r	theta_s	alpha	npar	ksat	.rosetta.model	.rosetta.version
1	22242478	Aluf	47	H1	0	127	94.4	0.6	5.0	NA	0.117	0.049	0.058	0.371	-1.536	0.490	2.779	2	1
1	22242478	Aluf	47	H2	127	203	91.0	1.0	8.0	1.55	NA	0.062	0.058	0.379	-1.541	0.402	2.488	3	1
1	22242479	Hitilo	22	H1	0	117	94.1	1.4	4.5	1.55	0.106	NA	0.055	0.381	-1.504	0.408	2.494	4	1
1	22242479	Hitilo	22	H2	117	137	55.4	17.6	27.0	1.58	0.252	0.181	0.080	0.388	-1.387	0.111	1.324	5	1
1	22242479	Hitilo	22	H3	137	157	59.3	13.7	27.0	1.58	0.261	0.156	0.056	0.390	-1.649	0.108	1.300	5	1
1	22242479	Hitilo	22	H4	157	203	62.8	19.2	18.0	1.50	0.211	0.122	0.052	0.394	-1.446	0.126	1.635	5	1
1	22243860	Aransas	100	H1	0	152	23.3	29.2	47.5	1.45	0.369	0.285	0.099	0.456	-1.457	0.072	0.974	5	1
1	22245205	Muskogee	20	H1	0	38	11.7	69.8	18.5	1.38	0.286	0.134	0.050	0.414	-1.983	0.144	1.519	5	1
1	22245205	Muskogee	20	H2	38	64	NA	63.2	30.0	1.35	0.308	0.165	NA	NA	NA	NA	NA	-1	1
NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	-1	1

Differences Between Versions

Water Retention Curves

# make a table of soil textures
# sand sand/silt/clay values at geometric centroids
tex <- SoilTextureLevels(simplify = FALSE)
x <- texcl_to_ssc(tex)
x <- cbind(x, tex)

# request results from all three versions
r1 <- ROSETTA(x, vars = c('sand', 'silt', 'clay'), v = '1')
r2 <- ROSETTA(x, vars = c('sand', 'silt', 'clay'), v = '2')
r3 <- ROSETTA(x, vars = c('sand', 'silt', 'clay'), v = '3')

# stack
r <- rbind(r1, r2, r3)

# iterate over results and generate VG model curve
res <- lapply(1:nrow(r), function(i) {
  
  # model bounds are given in kPA of suction
  vg <- KSSL_VG_model(VG_params = r[i, ], phi_min = 10^-3, phi_max=10^6)
  
  # extract curve and add texture ID
  m <- vg$VG_curve
  m$texture <- r$tex[i]
  
  # save model version
  m$.rosetta.version <- r$.rosetta.version[i]
  
  return(m)
})

# flatten to data.frame
res <- do.call('rbind', res)

# set factor levels / labels for model version
res$version <- factor(res$.rosetta.version, levels = c('1', '2', '3'), labels = c('ROSETTA 1', 'ROSETTA 2', 'ROSETTA 3'))

# plot style adjustments
tps <- tactile.theme(
  plot.line=list(col='royalblue', lwd = 2) 
)

# panels are soil texture classes
xyplot(
  phi ~ theta | texture, data = res, groups = version,
  type = c('l', 'g'), 
  auto.key = list(lines = TRUE, points = FALSE, columns = 3),
  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 = tps, 
  strip = strip.custom(bg=grey(0.85)), 
  as.table = TRUE,
  layout = c(7, 3),
  main = 'Idealized Water Retention'
)

Approximate AWC

# iterate over results
# generate VG model curve
# compute a simplistic AWC
awc <- lapply(1:nrow(r), function(i) {
  
  # model bounds are given in kPA of suction
  vg <- KSSL_VG_model(VG_params = r[i, ], phi_min = 10^-3, phi_max=10^6)
  
  # extract VWC at specific matric potentials (kPa)
  d <- data.frame(
    texture = r$tex[i], 
    sat = vg$VG_function(0),
    fc = vg$VG_function(33),
    pwp = vg$VG_function(1500)
  )
  
  # simplistic determination of AWC using 33 kPa -> 1500 kPa interval
  d$awc <- with(d, fc - pwp)
  
  # save model version and texture class
  d$.rosetta.version <- r$.rosetta.version[i]
  d$tex <- r$tex[i]
  
  return(d)
})

awc <- do.call('rbind', awc)

dotplot(
  factor(.rosetta.version) ~ awc | tex, data = awc,
  ylab = 'Rosetta Model Version', 
  scales = list(alternating = 1, x = list(tick.number = 5, cex = 0.66, rot = 90)), 
  xlab = expression(AWC~~(cm^3/cm^3)), 
  par.settings = tps, 
  strip = strip.custom(bg=grey(0.85)), 
  as.table = TRUE,
  layout = c(7, 3),
  main='Simplistic AWC (FC - PWP)',
  panel = function(...) {
    panel.abline(v = seq(0, 0.25, by = 0.05), col = trellis.par.get()$dot.line$col)
    panel.dotplot(...)
  }
)

Effect of Bulk Density

Water Retention Curves

# make a table of soil textures
# sand sand/silt/clay values at geometric centroids
tex <- SoilTextureLevels(simplify = FALSE)
x <- texcl_to_ssc(tex)
x <- cbind(x, tex)
x$Db1 <- 1.1
x$Db2 <- 1.25

# request results from all three versions
r1 <- ROSETTA(x, vars = c('sand', 'silt', 'clay', 'Db1'))
r2 <- ROSETTA(x, vars = c('sand', 'silt', 'clay', 'Db2'))

# stack
r <- make.groups(Db1 = r1, Db2 = r2)

# set factor levels / labels
r$which <- factor(r$which, levels = c('Db1', 'Db2'), labels = c('Db = 1.1 g/cc', 'Db = 1.25 g/cc'))


# iterate over results and generate VG model curve
res <- lapply(1:nrow(r), function(i) {
  
  # model bounds are given in kPA of suction
  vg <- KSSL_VG_model(VG_params = r[i, ], phi_min = 10^-3, phi_max=10^6)
  
  # extract curve and add texture ID
  m <- vg$VG_curve
  m$texture <- r$tex[i]
  
  # save Db set version
  m$which <- r$which[i]
  
  return(m)
})

# flatten to data.frame
res <- do.call('rbind', res)


# plot style adjustments
tps <- tactile.theme(
  plot.line=list(col='royalblue', lwd = 2) 
)

# panels are soil texture classes
xyplot(
  phi ~ theta | texture, data = res, groups = which,
  type = c('l', 'g'), 
  auto.key = list(lines = TRUE, points = FALSE, columns = 2),
  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 = tps, 
  strip = strip.custom(bg=grey(0.85)), 
  as.table = TRUE,
  layout = c(7, 3),
  main = 'Idealized Water Retention'
)

Approximate AWC

# iterate over results
# generate VG model curve
# compute a simplistic AWC
awc <- lapply(1:nrow(r), function(i) {
  
  # model bounds are given in kPA of suction
  vg <- KSSL_VG_model(VG_params = r[i, ], phi_min = 10^-3, phi_max=10^6)
  
  # extract VWC at specific matric potentials (kPa)
  d <- data.frame(
    texture = r$tex[i], 
    sat = vg$VG_function(0),
    fc = vg$VG_function(33),
    pwp = vg$VG_function(1500)
  )
  
  # simplistic determination of AWC using 33 kPa -> 1500 kPa interval
  d$awc <- with(d, fc - pwp)
  
  # save Db set and texture class
  d$which <- r$which[i]
  d$tex <- r$tex[i]
  
  return(d)
})

awc <- do.call('rbind', awc)

dotplot(
  which ~ awc | tex, data = awc,
  ylab = '', 
  scales = list(alternating = 1, x = list(tick.number = 5, cex = 0.66, rot = 90)), 
  xlab = expression(AWC~~(cm^3/cm^3)), 
  par.settings = tps, 
  strip = strip.custom(bg=grey(0.85)), 
  as.table = TRUE,
  layout = c(7, 3),
  main='Simplistic AWC (FC - PWP)',
  panel = function(...) {
    panel.abline(v = seq(0, 0.25, by = 0.05), col = trellis.par.get()$dot.line$col)
    panel.dotplot(...)
  }
)

Bulk Density and Ksat

TODO

This document is based on soilDB version 2.7.6.

ROSETTA Model API

Dylan Beaudette, Richard Reid, Todd Skaggs

2022-11-29

Introduction

Model Versions

Setup

A Very Basic Example

Theoretical Water Retention Curves

Automatic Model Selection

Differences Between Versions

Water Retention Curves

Approximate AWC

Effect of Bulk Density

Water Retention Curves

Approximate AWC

Bulk Density and Ksat

ROSETTA Model API

Dylan Beaudette, Richard Reid, Todd Skaggs

2022-11-29

Introduction

Model Versions

Related Documentation

Setup

A Very Basic Example

Theoretical Water Retention Curves

Automatic Model Selection

Differences Between Versions

Water Retention Curves

Approximate AWC

Effect of Bulk Density

Water Retention Curves

Approximate AWC

Bulk Density and Ksat