Introduction

Setup

You will need these packages. Must use the development version of hydromad, requires a working compiler (RTools on Windows).

library(aqp)
library(soilDB)
library(sharpshootR)

# special installation instructions:
# https://github.com/josephguillaume/hydromad#installation
library(hydromad)

Simple Example

Simulate monthly water balances for 50cm of sandy loam vs. 50cm of silt loam soil material in a xeric/thermic climate.

data("ROSETTA.centroids")

idx <- which(ROSETTA.centroids$texture %in% c('sandy loam', 'silt loam'))
vars <- c('texture', 'pwp', 'fc', 'sat', 'awc')

knitr::kable(
  ROSETTA.centroids[idx, vars],
  row.names = FALSE
)
texture pwp fc sat awc
sandy loam 0.0624426 0.167354 0.3870001 0.1049114
silt loam 0.0858854 0.293965 0.4390000 0.2080796 
# near Sonora, CA
# thermic / xeric climate
# units are mm
PET <- c(0, 0, 5,80, 90, 120, 130, 140, 110, 90, 20, 5)
PPT <- c(0, 150, 200, 120, 20, 0, 0, 0, 10, 20, 30, 60)

# 50cm of soil material -> convert to mm
thick <- 50 * 10

# total water storage for sandy loam texture soil material
AWC.sl <- round(0.1049114 * thick)

# total water storage for silt loam texture soil material
AWC.sil <- round(0.2080796 * thick)

# simulate water balances:
# * run for 5 cycles
# * start in January
# * start with "completely dry" soil conditions
# * keep only the last cycle
wb.sl <- monthlyWB(AWC.sl, PPT, PET, S_init = 0, starting_month = 1, rep = 5, keep_last = TRUE, distribute = TRUE)
wb.sil <- monthlyWB(AWC.sil, PPT, PET, S_init = 0, starting_month = 1, rep = 5, keep_last = TRUE, distribute = TRUE)

All values are mm. U: surplus, S: soil water storage, ET: actual ET, D: deficit.

kable(wb.sil, digits = 2, row.names = FALSE)
PPT PET U S ET D month mo .sim
0 0 0.00 95.94 0.00 0.00 1 Jan 5
150 0 141.94 104.00 0.00 0.00 2 Feb 5
200 5 195.00 104.00 5.00 0.00 3 Mar 5
120 80 40.00 104.00 80.00 0.00 4 Apr 5
20 90 0.00 54.72 69.28 -20.72 5 May 5
0 120 0.00 16.06 38.66 -81.34 6 Jun 5
0 130 0.00 4.22 11.83 -118.17 7 Jul 5
0 140 0.00 1.00 3.23 -136.77 8 Aug 5
10 110 0.00 6.02 4.98 -105.02 9 Sep 5
20 90 0.00 15.01 11.01 -78.99 10 Oct 5
30 20 0.00 39.36 5.65 -14.35 11 Nov 5
60 5 0.00 95.94 3.42 -1.58 12 Dec 5 
# compare
par(mar=c(4,4,4,1), bg = 'white', mfcol = c(1, 2))
plotWB(WB = wb.sl, legend.cex = 0.75)
title('50cm of sandy loam texture soil material', line = 3)

plotWB(WB = wb.sil, legend.cex = 0.75)
title('50cm of silt loam texture soil material', line = 3)

# compare
par(mar=c(4,4,4,1), bg = 'white', mfcol = c(1, 2))
plotWB_lines(WB = wb.sl)
title('50cm of sandy loam texture soil material', line = 3)

plotWB_lines(WB = wb.sil)
title('50cm of silt loam texture soil material', line = 3)

Distributing Monthly Values

As of sharpshootR 2.2, it is possible to “distribute” monthly data into k bins within each month.

Run a single water balance cycle, starting with a completely dry soil. Test the effect of distributing monthly values into 10 bins per month.

wb.sl1 <- monthlyWB(AWC.sil, PPT, PET, S_init = 0, starting_month = 1, rep = 1, distribute = FALSE)

wb.sl2 <- monthlyWB(AWC.sil, PPT, PET, S_init = 0, starting_month = 1, rep = 1, , distribute = TRUE, k = 10, method = 'equal')
par(mar=c(4,4,4,1), bg = 'white', mfcol = c(1, 2))
plotWB(WB = wb.sl1, legend.cex = 0.75)
title('50cm of silt loam texture soil material\nMonthy', line = 1.5)

plotWB(WB = wb.sl2, legend.cex = 0.75)
title('50cm of silt loam texture soil material\nEqual Distribution, 10 bins', line = 1.5)

This time, test the effects of using ‘equal’ vs. ‘random’ distribution of monthly PPT and PET into 10 bins per month.

wb.sl1 <- monthlyWB(AWC.sil, PPT, PET, S_init = 0, starting_month = 1, rep = 1, distribute = TRUE, k = 10, method = 'equal')

wb.sl2 <- monthlyWB(AWC.sil, PPT, PET, S_init = 0, starting_month = 1, rep = 1, distribute = TRUE, k = 10, method = 'random')
par(mar=c(4,4,4,1), bg = 'white', mfcol = c(1, 2))
plotWB(WB = wb.sl1, legend.cex = 0.75)
title('50cm of silt loam texture soil material\nEqual Distribution, 10 bins', line = 1.5)

plotWB(WB = wb.sl2, legend.cex = 0.75)
title('50cm of silt loam texture soil material\nRandom Distribution, 10 bins', line = 1.5)

Real Data

Series-level summaries of climate provided by fetchOSD().

AWC estimated from SSURGO data, via SDA.

s <- 'amador'
x <- fetchOSD(s, extended = TRUE)

# get representative, profile-total AWC from SSURGO
sql <- sprintf("
SELECT chorizon.cokey AS cokey, 
SUM(awc_r * (hzdepb_r - hzdept_r)) AS ws 
FROM 
legend JOIN mapunit ON legend.lkey = mapunit.lkey
JOIN component ON mapunit.mukey = component.mukey
JOIN chorizon ON component.cokey = chorizon.cokey 
WHERE compname = '%s'
AND areasymbol != 'US'
GROUP BY chorizon.cokey;", s
)

# get via SDA
res <- SDA_query(sql)

# median AWC in mm
# over all components correlated to named series 
AWC <- round(median(res$ws, na.rm = TRUE) * 10)


# monthly climate data from series summary
PPT <- x$climate.monthly$q50[x$climate.monthly$variable == 'Precipitation (mm)']
PET <- x$climate.monthly$q50[x$climate.monthly$variable == 'Potential ET (mm)']

Water year.

# tighter margins
par(mar=c(4,4,2,1), bg = 'white')

# water year
# last iteration
x.wb <- monthlyWB(AWC, PPT, PET, S_init = 0, starting_month = 9, rep = 3, keep_last = TRUE, distribute = TRUE)
plotWB(WB = x.wb)

# convert total ETa into inches
sum(x.wb$ET) * 0.03937
## [1] 9.246288

Calendar year.

# tighter margins
par(mar=c(4,4,2,1), bg = 'white')

# calendar year
# last iteration
x.wb <- monthlyWB(AWC, PPT, PET, S_init = 0, starting_month = 1, rep = 3, keep_last = TRUE, distribute = TRUE)
plotWB(WB = x.wb)

Bar vs. Line Graph

s <- 'auburn'
x <- fetchOSD(s, extended = TRUE)

# get representative, profile-total AWC from SSURGO
sql <- sprintf("SELECT chorizon.cokey AS cokey, 
SUM(awc_r * (hzdepb_r - hzdept_r)) AS ws 
FROM 
legend JOIN mapunit ON legend.lkey = mapunit.lkey
JOIN component ON mapunit.mukey = component.mukey
JOIN chorizon ON component.cokey = chorizon.cokey 
WHERE compname = '%s'
AND areasymbol != 'US'
GROUP BY chorizon.cokey;", s)

# get via SDA
res <- SDA_query(sql)

# median AWC in mm
# over all components correlated to named series 
AWC <- round(median(res$ws, na.rm = TRUE) * 10)

# monthly climate data from series summary
PPT <- x$climate.monthly$q50[x$climate.monthly$variable == 'Precipitation (mm)']
PET <- x$climate.monthly$q50[x$climate.monthly$variable == 'Potential ET (mm)']

# 3 warm-up cycles
# keep last iteration
# calendar year
x.wb <- monthlyWB(AWC, PPT, PET, S_init = 0, starting_month = 1, rep = 3, keep_last = TRUE, distribute = TRUE)

# tighter margins
par(mar=c(4,4,3,1), bg = 'white', mfcol = c(1, 2))

plotWB(x.wb)
title(sprintf('Monthly Water Balance: %s Series', toupper(s)), line = 2)

plotWB_lines(x.wb)
title(sprintf('Monthly Water Balance: %s Series', toupper(s)), line = 2)

Comparisons

Abstract into a function, and use to compare multiple soils / climatic parameters. A more efficient approach would load all of the data then iterate over chunks.

# x: soil series name
compareMonthlyWB <- function(x, ylim) {
  
  # get climate data
  osd <- fetchOSD(x, extended = TRUE)
  
  # get representative, profile-total AWC from SSURGO
  sql <- sprintf("
SELECT chorizon.cokey AS cokey, 
SUM(awc_r * (hzdepb_r - hzdept_r)) AS ws 
FROM 
legend JOIN mapunit ON legend.lkey = mapunit.lkey
JOIN component ON mapunit.mukey = component.mukey
JOIN chorizon ON component.cokey = chorizon.cokey 
WHERE compname = '%s'
AND areasymbol != 'US'
GROUP BY chorizon.cokey;", x
  )
  
  # get via SDA
  res <- suppressMessages(SDA_query(sql))
  
  # median AWC in mm
  # over all components correlated to named series 
  AWC <- round(median(res$ws, na.rm = TRUE) * 10)
  
  
  # monthly climate data from series summary
  PPT <- osd$climate.monthly$q50[osd$climate.monthly$variable == 'Precipitation (mm)']
  PET <- osd$climate.monthly$q50[osd$climate.monthly$variable == 'Potential ET (mm)']
  
  # do water balance simulation
  x.wb <- monthlyWB(AWC, PPT, PET, S_init = 0, starting_month = 1, rep = 3, keep_last = TRUE, distribute = TRUE)
  
  # compose WB figure
  plotWB(WB = x.wb, legend.cex = 0.9, ylim = ylim)
  
  # annotate with series name / family level classification
  s <- site(osd$SPC)
  txt <- sprintf("%s\n%s", s$id, s$family)
  title(txt, line = 3, cex.main = 1)
}
# custom margins / multi-figure
par(mar=c(4,4,5,1), bg = 'white', mfrow = c(2, 2))

# iterate over these soil series names
soils <- c('amador', 'lucy', 'drummer', 'pierre')
for(i in soils) {
  # adjust ylim manually
  compareMonthlyWB(i, ylim = c(-250, 180))
}

# custom margins / multi-figure
par(mar=c(4,4,5,1), bg = 'white', mfrow = c(2, 2))

# iterate over these soil series names
soils <- c('leon', 'tristan', 'miami', 'freznik')
for(i in soils) {
  # adjust ylim manually
  compareMonthlyWB(i, ylim = c(-180, 190))
}

This document is based on sharpshootR version 2.2 and soilDB version 2.8.1.