rain barrel

rain.Rmd

# Rain Barrel Simulation

I want to decide how big my rain barrel has to be. So I use NOAA's version of YYZ's weather

        mkdir 71624099999
        for f in `seq 2006 2022`; do
                wget https://ncei.noaa.gov/data/global-hourly/access/$f/71624099999.csv -O $f
        done
        rm 71624099999/2005.csv # is missing for some reason

This gives rain, wind, dry bulb temperature and dew point at 6 or 12 hourly intervals:

```{r}
library(tidyverse)
library(ggplot2)
library(lubridate)
library(Rcpp)
#Date.daily - date in daily time step,
#Date.monthly - date in monthly time step,
#J - julian days,
#i - month,
#ndays - days in month,
#Tmax - daily maximum temperature in degree Celcius,
#Tmin - daily minimum temperature in degree Celcius,
#RHmax - daily maximum relative humidity in percentage,
#RHmin - daily minimum relative humidity in percentage,
#uz - daily wind speed in meters per second,
#Rs - daily solar radiation in Megajoule per square meter. 
# solar radiation from GF1 sky cover. But how?

if (file.exists("71624099999.xz")) load("71624099999.xz") else {
rain <- Sys.glob("71624099999/*.csv") %>%
        map(read_csv) %>%
        map_dfr(~ select(., DATE, WND, DEW, TMP, AA1, GF1)) %>%
        separate(WND, c("wdeg", "wdegq", "wtype","wspeed","wspeedq"), convert=T) %>%
        separate(DEW, c("dew", "dewq"), sep=",", convert=T) %>%
        separate(TMP, c("tmp", "tmpq"), sep=",", convert=T) %>%
        separate(AA1, c("rainint","rainmm", "raintrace", "rainquality"), convert=T) %>%
        separate(GF1, c("oktat","oktato", "oktatq", "oktatl", "oktattlq", "cloudgenus", "cloudgenusq",
                        "lcbh", "lcbhq", "mcgh", "mcgq", "hcg", "hcgq"), convert=T)
save(rain, file="71624099999.xz")
}
rain <- rain %>% select(-oktato) %>% # all 99s
        filter(dew != 9999,
                        tmp != 9999, wspeed != 9999, wdeg != 999,
                        rainmm != 9999,
                        rainquality != 2,
                        # !is.na(oktat), oktat != 99,
                        dewq != 2,
                        year(DATE) >= 2006) %>%
        mutate(dew = dew/10, tmp = tmp/10, wspeed = wspeed/10,
               rainmm = rainmm/ 10)
```

### Linacre 1993 equation 12
```{r}
# z is elevation in metres
linacreEvap <- function(tmp, dew, u, A=45, z=100) {
        (0.0015 + 4e-4 * tmp + 1e-6 * z) *
                (480 * (tmp + 0.006*z) / (84 - A) - 40 + 2.3 * u * (tmp - dew))

}
```




## Irrigation needs
One rule of thumb is to water grass 2.5cm per week. However I do a better method based on [pan evaporation](https://www.agric.wa.gov.au/water-management/evaporation-based-irrigation-scheduling?page=0%2C1). Unfortunately, the weather station does not provide pan evaporation. One option was to learn pan evaporation rates from other weather stations. But the data has anomalies such as 1m evaporation per day without a clear separation between reasonable and unreasonable values. Instead I use separate correlations for [heat transfer](https://www.engineeringtoolbox.com/convective-heat-transfer-d_430.html) and [mass transfer](https://www.engineeringtoolbox.com/evaporation-water-surface-d_690.html) from [engineeringtoolbox.com](http://www.engineeringtoolbox.com). If the pan is initially at the dry bulb temperature, it cools down by evaporation. As the temperature drops, evaporation slows down and heat supplied by convection increases. Eventually the two equal, which determines the pan water temperature and evaporation rate.

```{r}
# Stull 1947 coefficients via NIST https://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185&Mask=4&Type=ANTOINE&Plot=on#ANTOINE
antoine <- function(celsius) {
        with(list(A = 4.6543, B = 1435.264, C = -64.848),
          10^(A - B/(celsius+C+273.15)))
}
# https://www.engineeringtoolbox.com/convective-heat-transfer-d_430.html (2b)
hcW <- function(v) 12.12 - 1.16*v + 11.6*sqrt(v)

# (1) https://www.engineeringtoolbox.com/evaporation-water-surface-d_690.html
theta <- function(v) 25 + 19*v
evapW <- function(v, twater, RH, P=1) {
        Psat <- antoine(twater)
        xs <- 0.62198 * Psat / (P - Psat)
        x <- 0.62198 * Psat * RH / (P - Psat*RH)
        theta(v) * (xs - x) / 3600 # kg/m2/s
}
hvap <- 2400 * 1000 # j/kg

twaterObj <- function(twater, celsius=20, v=3, RH=0.5) hcW(v)*(celsius - twater) - hvap*evapW(v, twater, RH)
getTwater <- function(celsius, ...) uniroot(function(x) twaterObj(x, celsius=celsius, ...), c(-30, celsius))

evaporationRate <- function(tmp, dew, v) {
        RH <- antoine(dew) / antoine(tmp)
        tw <- getTwater(tmp, v, RH)$root
        evapW(v, tw, RH) * 3600 # mm / hr
}
rain <- within(rain, {
        evap.mmph.me <- Vectorize(evaporationRate)(tmp, dew, wspeed)
        evap.mmph <- linacreEvap(tmp, dew, wspeed) / 24
          })
```
```{r}
ggplot(rain, aes(evap.mmph.me, evap.mmph)) + geom_point() + geom_smooth()
```

Total evaporation is 3x rain. Which I think is wrong. Linacre's version says evaporation is less than rain,
but very close. This is consistent with lakes and rivers being common. But Linacre's formula is for lake evaporation
not for pan evaporation.

```{r}
sum(c(6, diff(rain$DATE)) * rain$evap.mmph)
sum(rain$rainmm)
```


The roof is 36m2 for the outlet closest to the garden, 64 on the southeast, and the front is around 25m2 for each east and west. 

Next is the simulation of the rain barrel.

```{r}
cppFunction('List barrelSoilCPP(
        IntegerVector hours, // 5 vectors derived from the `rain` tibble
        IntegerVector yday,  // day of the year 1 to 366
        IntegerVector month, // month 1 (January) to 12
        NumericVector mm,
        NumericVector evapmmph, // mm per hour evaporation
        double consume24h,
        double roofm2,
        double cropm2, // crop parameters
        double cropevapfac,
        double soilmmMax,
        double barrelL, // barrel size
        NumericVector barrelDrinkL) {
  // barrelDrinkL[j] is the amount in the barrel
  // that cannot be used to water plants on day j
  // j goes from 1 to 366
  int n = hours.size();
  NumericVector soilmm(n), holdup(n), watering(n), distill(n);
  holdup[0] = 2;
  soilmm[0] = soilmmMax;
  for (int i=1; i < n; i++) {
          holdup[i] = holdup[i-1] + mm[i]*roofm2 - (hours[i] /24.0) * consume24h ;
          if (holdup[i] > barrelL) holdup[i] = barrelL; // overflow
          if (holdup[i] < 0) {
                  distill[i] = -holdup[i];
                  holdup[i] = 0;
          }

          soilmm[i] = soilmm[i-1] + mm[i] - evapmmph[i] * hours[i] * cropevapfac;
          if (soilmm[i] > soilmmMax) soilmm[i] = soilmmMax; // overflow
          if (soilmm[i] < 0) {
                double h0 = holdup[i];
                holdup[i] = holdup[i] + soilmm[i] * cropm2;
                // refill
                if (holdup[i] < barrelDrinkL[yday[i]]) {
                        if (h0 > barrelDrinkL[yday[i]]) h0 = barrelDrinkL[yday[i]];
                        watering[i] = h0 - holdup[i];
                        holdup[i] = h0;
                }
                soilmm[i] = 0;
          }
    }
    return List::create(
           _("holdup") = holdup,
           _("distill") = distill,
           _("watering") = watering,
           _("soilmm") = soilmm);
}')

barrelSoil <- function(raindf = rain,
                       consume24h=2,
                       roofm2=36,
                       cropm2=25,
                       crop.evap.fac=1.1,
                       soilmm.max=100,
                       barrelL=55,
                       barrelDrinkL=4,
                       cpp=TRUE) {
    barrelDrinkL <- spline(
           x = seq(1, 365, length.out = 1 +length(barrelDrinkL)),
           y = c(barrelDrinkL, barrelDrinkL[1]),
           method = "periodic",
           xout = 1:366)$y
    barrelSoilCPP(c(6, diff(raindf$DATE)),
                  yday(raindf$DATE),
                  month(raindf$DATE),
                  raindf$rainmm,
                  raindf$evap.mmph,
                  consume24h,
                  roofm2,
                  cropm2,
                  crop.evap.fac,
                  soilmm.max,
                  barrelL,
                  barrelDrinkL) %>%
            as_tibble
}
b1 <- barrelSoil()
```
```{r}
ggplot(cbind(b1, rain) %>% filter(distill < 1000),
       aes(DATE, cumsum(watering))) + geom_smooth() + geom_point()
```


```{r}
# quite slow possibly RCPP for the barrelSoil function?
bWaterParam <- expand_grid(
                       barrelL=c(8, 18, 55, 110, 200, 300, 500, 1000),
                       barrelDrinkL=c(4:30, seq(32, 100, by=5) )) %>%
        filter(barrelL > barrelDrinkL*2) %>%
        rowwise %>%
        mutate(sim = barrelSoil(barrelL = barrelL, barrelDrinkL = barrelDrinkL) %>% list)

bwp <- bWaterParam %>%
        mutate(watering = sim$watering %>% sum,
               wateringw = sim[month(rain$DATE) %in% c(5,6,7,8,9), "watering" ]  %>% sum,
                distill = sim$distill %>% sum,
               distillw = sim[month(rain$DATE) %in% c(5,6,7,8,9), "distill" ]  %>% sum,
               waterCharge = distillw * 0.15 + wateringw * 0.002 )
ggplot(bwp, aes(barrelL, waterCharge, col=barrelDrinkL )) + geom_point()
```

So I minimize waterCharge, by selecting barrelDrinkL

```{r}
library(nloptr)
bdlObj <- function(bdl, barrelL, ...) { 
        bs <- barrelSoil(barrelL=barrelL,
                         barrelDrinkL = bdl,
                         ...)[month(rain$DATE) %in% (5:9),]
        sum(bs$watering*0.002 + bs$distill*0.15)
}

# now do multiple variables
bdlOpt <- function(barrelL, ndf, ...) {
        o <- optimize(bdlObj, c(0, barrelL), barrelL = barrelL, ...)
        p <- cobyla(x0 = rep(o$minimum, ndf), fn = bdlObj,
               lower=rep(0, ndf), upper=rep(barrelL, ndf), barrelL = barrelL, ...)
        # and then compare o$objective with p$value
        p$policy12Gain <- o$objective - p$value
        p
}



bdlOpts <- expand_grid(barrelL=c(18, 36, 55),
                       roofm2=c(36, 100),
                       ndf = c(4, 8, 12)) %>%
        rowwise %>%
        mutate(out = list(bdlOpt(barrelL, ndf=ndf, roofm2 = roofm2)),
               waterCharge = out$value,
               policy12Gain = out$policy12Gain,
               monthly = list(tibble(barrelDrinkL = out$par,
                                     pari = seq_along(out$par))))
```
```{r}
ggplot(bdlOpts, aes(barrelL, waterCharge, col=roofm2, group=roofm2)) + geom_point() + geom_line() +
        geom_abline(slope=1, intercept=0)
```

```{r}
bdlOptsProf <- bdlOpts %>% unnest(monthly) %>% group_by(barrelL, roofm2, ndf, waterCharge) %>%
 transmute(xy = spline(seq(1, 366, length=length(barrelDrinkL)+1), c(barrelDrinkL, barrelDrinkL[1]), xout=1:365, method="periodic") %>% as_tibble %>% 
        list) %>% unnest(xy)

ggplot(bdlOptsProf, aes(x,y, col=barrelL, group=interaction(barrelL, roofm2, ndf))) +
        geom_line() + facet_wrap(~ ndf) +
        ggtitle("monthly barrelL is overfit")
```

```{r}
ggplot(bdlOpts, aes(barrelL, policy12Gain, group=interaction(roofm2, ndf), lty=factor(ndf), col=roofm2)) + geom_line() +
        ggtitle("barely any savings ($2 over 20 years) for a different drinking volume each month")
```