A Simple Monthly Water Balance Simulation

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.

knitr::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.265973

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)
  
  return(x.wb)
}

# 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')
res <- lapply(soils, compareMonthlyWB, ylim = c(-250, 180))

Experimental: summarize water balance results.

res <- do.call('rbind', lapply(res, monthlyWB_summary))
res$series <- soils

knitr::kable(res, digits = 2)

dry	moist	wet	dry_con	moist_con	wet_con	total_deficit	total_surplus	total_AET	annual_AET_PET_ratio	series
61	183	122	61	91	91	-618.64	225.14	235.36	0.28	amador
0	183	183	0	183	122	-121.66	434.66	833.34	0.87	lucy
0	152	213	0	152	152	-65.11	335.11	633.89	0.91	drummer
0	304	61	0	244	61	-223.56	25.56	407.44	0.65	pierre

# 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')
res <- lapply(soils, compareMonthlyWB, ylim = c(-180, 190))