Closed
Description
There may be a problem with hydromad::bucket.sim: hydromad/hydromad#188
Until that time, consider using an R version of the algorithm, borrowed / modified from: https://github.com/josephguillaume/hydromad/blob/master/R/bucket.R
TODO
- review code, related literature, specifically behavior at
S_prev = 0,P[t]= small, andE[t] >> P[t]. - contact hydromad authors (bucket model, AET values confusing when S=0 hydromad/hydromad#188)
- contact authors of Bai et al., 2009
- test against manually performed WB
hydromad::bucket.sim() implements model "S2" as defined in Bai et al., 2009. The possible "bug" occurs in step 2:
Until resolved, all related functions in sharpshootR are going to be a lot slower with the (current) R-based implementation.
Consider implementation of the S1-S4 and M1-M4 version once all of this is sorted out.
library(sharpshootR)
## monthly water balance
# example data
d <- structure(list(P = c(65, 59, 57, 28, 13, 3, 0, 1, 4, 20, 33,
53), E = c(14, 22, 38, 54, 92, 125, 154, 140, 106, 66, 29, 14
)), class = "data.frame", row.names = c(NA, -12L))
# original implementation
m <- hydromad::hydromad(d, sma = "bucket", routing = NULL)
m <- update(m, Sb = 47, fc = 1, S_0 = 1, a.ss = 0.01, M = 0, etmult = 1, a.ei = 0)
res <- predict(m, return_state = TRUE)
# combine original PPT,PET with results
res <- data.frame(d, res)
# cleanup names
names(res) <- c('PPT', 'PET', 'U', 'S', 'ET')
# with "fix" that ensures AET cannot be > S[t]
res2 <- sharpshootR:::.leakyBucket(d, Sb = 47, fc = 1, S_0 = 1, a.ss = 0.01, M = 0, etmult = 1, a.ei = 0)
# combine original PPT,PET with results
res2 <- data.frame(d, res2)
# cleanup names
names(res2) <- c('PPT', 'PET', 'U', 'S_new', 'ET_new')
# note: U does not change
# compare differences
x <- data.frame(
PPT = res$PPT,
PET = res$PET,
U = res$U,
S = res$S,
S_new = res2$S_new,
dS = ifelse(res$S != res2$S_new, '*', ''),
ET = res$ET,
ET_new = res2$ET_new,
dET = ifelse(res$ET != res2$ET_new, '*', '')
)
knitr::kable(x, digits = 1)| PPT| PET| U| S| S_new|dS | ET| ET_new|dET |
|---:|---:|----:|----:|-----:|:--|----:|------:|:---|
| 65| 14| 51.0| 47.0| 47.0| | 14.0| 14.0| |
| 59| 22| 37.0| 47.0| 47.0| | 22.0| 22.0| |
| 57| 38| 19.0| 47.0| 47.0| | 38.0| 38.0| |
| 28| 54| 0.0| 21.0| 28.0|* | 54.0| 47.0|* |
| 13| 92| 0.0| 0.0| 0.0| | 66.6| 41.0|* |
| 3| 125| 0.0| 0.0| 0.0| | 8.0| 3.0|* |
| 0| 154| 0.0| 0.0| 0.0| | 0.0| 0.0| |
| 1| 140| 0.0| 0.0| 0.0| | 3.0| 1.0|* |
| 4| 106| 0.0| 0.0| 0.0| | 9.0| 4.0|* |
| 20| 66| 0.0| 0.0| 0.0| | 28.1| 20.0|* |
| 33| 29| 0.0| 12.6| 12.6| | 20.4| 20.4| |
| 53| 14| 4.6| 47.0| 47.0| | 14.0| 14.0| |
Metadata
Assignees
Labels
No labels