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/cm3

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

  • optional, volumetric water content at 1500 kPa (15 bar): 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)

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.075 0.042 0.035 0.367 -1.244 0.273 2.705 5 1
Aluf 50 H2 117 251 61.0 19.0 20.0 1.55 0.254 0.117 0.041 0.381 -1.906 0.133 1.229 5 1
Hitilo 30 H1 0 137 94.9 0.6 4.5 1.55 0.077 0.040 0.032 0.381 -1.228 0.258 2.736 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'
)