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 (any moisture basis): mass per volume after accounting for >2mm fragments, units of gram/cm3

  • optional, volumetric water content at 33 kPa: roughly “field capacity” for most soils, units of cm3/cm3

  • optional, volumetric water content at 1500 kPa: roughly “permanent wilting point” for most plants, units of cm3/cm3

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)

You will also need the latest version of soilDB:

install.packages('remotes', dep=TRUE)
remotes::install_github("ncss-tech/soilDB", dependencies = FALSE, upgrade_dependencies=FALSE)

A Very Basic Example

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

# required libraries
library(soilDB)
library(httr)
library(jsonlite)

# 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()
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 20392868 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 20392868 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 20392870 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 20392870 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 20392870 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 20392870 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 20393517 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 20394343 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 20394343 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()
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
  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(...)
  }
)