distribute.firstflush <- function(firstflush, #a firstflush depth (mm/d) hrrunoff){ #24 hourly runoff depth values (mm/d) # this function distributes daily firstflush depth values (as provided to the # daily tankmodel.R model) into the first hourly runoff records for that day. ## extract hourly records in hrrunoff with >0 runoff and order them ## order them from first hour to last hour rn <- hrrunoff[which(hrrunoff>0)[order(which(hrrunoff>0),decreasing = FALSE)]] ##create a vector of the cumulative sum of rn cumrn <- cumsum(rn) ##set up vector for resulting distribution of firstflush ff <- vector(mode = "numeric",length = 24) if(firstflush > sum(rn)) stop(cat("i = ",i,"dover", dover[i],"rn", rn,"firstflush depth exceeds rainfall depth"), call. = FALSE) if(firstflush <= cumrn[1]){ ##if it is less than the last hour's rain, all of dover can go in the last hour ff1 <- firstflush} else { ##otherwise firstflush is distributed across all the hours with a cumulative value < firstflush ff1 <- c(rn[cumrn < firstflush], #and the remainder (dover - last cumulative sum value) is distributed to the next earliest hour firstflush - cumrn[cumrn < firstflush][sum(cumrn < firstflush)]) } ff[which(hrrunoff>0)[order(which(hrrunoff>0), decreasing = FALSE)]][1:length(ff1)] <- ff1 # ff } distribute.overflow <- function(dover, #a daily overflow depth (mm/d) hrrunoff){ #a vector of 24 hourly runoff depth values (mm/h) #this function distributes daily overflow depth values (as produced by the #daily tankmodel.R model) into the last hourly runoff records for that day. ## extract hourly records in hrrunoff with >0 runoff and order them ## order them from last hour to first hour rn <- hrrunoff[which(hrrunoff>0)[order(which(hrrunoff>0),decreasing = TRUE)]] ##create a vector of the cumulative sum of rn cumrn <- cumsum(rn) ##set up matrix for resulting distributions of dover across each value in dover over <- rep(0,24) if(round(sum(rn),2) < round(dover,2)) stop("overflow depth exceeds rainfall depth", call. = FALSE) if(dover <= cumrn[1]){ ##if it is less than the last hour's rain, all of dover can go in the last hour over1 <- dover} else { ##otherwise dover is distributed across all the hours with a cumulative value < dover #this line avoids for rounding errors in determination of dover depth if(dover - cumrn[cumrn < dover][sum(cumrn < dover)] < 0.01){ over1 <- rn[cumrn < dover]}else { #and the remainder (dover - last cumulative sum value) is distributed to the next earliest hour over1 <- c(rn[cumrn < dover], dover - cumrn[cumrn < dover][sum(cumrn < dover)]) } } over[which(hrrunoff>0)[order(which(hrrunoff>0), decreasing = TRUE)]][1:length(over1)] <- over1[over1 > 0] # over } compile.inflow <- function(hrunoff, #data.frame of hourly runoff in mm (called runoff.mm), #with cols datetime,date,year,month, (and final col daten) formatted using chron) carea, #catchment area for hrunoff in sq m doverflow = NULL, #data frame with first field daten (formated with chron) #and 2nd column from daily overflow data from, #tank model $budget$overflow (in L/d) firstflush = NA, #if !is.na, = firstflush in L from tankmodel.R carea.tank = 0){ #catchment area for tank producing #overflow in sq m if(sum(unique(hrunoff$date) != as.numeric(doverflow$date))>0) stop("Dates for the two periods do not match", call. = FALSE) if(is.null(doverflow) & (is.na(firstflush) | firstflush == 0)){ inflow <- data.frame(datetime = hrunoff$datetime, date = hrunoff$date, year = hrunoff$year, month = hrunoff$month, inflow.L = hrunoff$runoff.mm*carea) runoffandover <- NULL } else { drunoff <- aggregate(hrunoff$runoff.mm, by = list(daten = hrunoff$daten), sum) doverflow <- as.vector(doverflow[,2])/carea.tank if(carea > 0){real.hro <- hrunoff$runoff.mm} else {real.hro <- rep(0, length(hrunoff$runoff.mm))} runoffandover <- data.frame(daten = hrunoff$daten, runoff.mm = real.hro, over = rep(0,length(hrunoff$runoff.mm))) #runoffandover$runoff.mm[runoffandover$runoff.mm < 1] <- 0 #names(runoffandover)[3:(dim(doverflow)[2]+2)] <- letters(1:nd) if(is.na(firstflush) | firstflush == 0){ if(sum(doverflow >0) == 0) { runoffandover[,3] <- 0 } else { indices <- which(doverflow >0) dates <- drunoff$daten[doverflow > 0] for(i in 1:length(indices)){ index <- indices[i] runoffandover[runoffandover$daten == dates[i], 3] <- distribute.overflow(doverflow[index], hrunoff$runoff.mm[hrunoff$daten == dates[i]]) } } }else{ indices <- which(drunoff[,2] >0) dates <- drunoff$daten[drunoff[,2] > 0] ff <- as.vector(apply(cbind(drunoff[,2],rep(firstflush/carea.tank,dim(drunoff)[1])), 1, min)) ff[drunoff[,2] < ff] <- drunoff[drunoff[,2] < ff, 2] for(i in 1:length(indices)){ index <- indices[i] if(doverflow[index] == 0){ # c(firstflush, rep(0,23)) runoffandover[runoffandover$daten == dates[i],3] <- distribute.firstflush(ff[index], hrunoff$runoff.mm[hrunoff$daten == dates[i]]) }else{ runoffandover[runoffandover$daten == dates[i],3] <- distribute.overflow(as.numeric(doverflow[index]), hrunoff$runoff.mm[hrunoff$daten == dates[i]]) + # c(firstflush, rep(0,23)) distribute.firstflush(ff[index], hrunoff$runoff.mm[hrunoff$daten == dates[i]]) } } } inflow <- data.frame(datetime = hrunoff$datetime, date = hrunoff$date, year = hrunoff$year, month = hrunoff$month, inflow.L = runoffandover[,3]*carea.tank + real.hro*carea) } list(inflow = inflow, runoffandover = runoffandover) } gardenmodel <- function(inflow, #data.frame with hourly inflow vol in L and #first three columns datetime, date, year, month #all class chron (compiled #using function compile.inflow()) and one column of inflow data hrunoff, #as for compile inflow - only used to calculate urban N load et,#vector the same length as inflow$datetime with #matching evapotranspiration values in mm/h (sets to zero if isveg is true) Af, Pf, Hf, #filter area (sq m) perimeter (m) and depth (m) Ap, Hp, #pond area (sq m) and depth (m) isveg = TRUE, lined.bottom = TRUE, # logical Ksu = 0 if TRUE, 0.0005 if FALSE (dm/h) lined.side = TRUE, # logical - can only be true if lined.bottom is TRUE medium = "loamy sand", #filter medium, must be sandy loam, sand, or gravel (scoria) #decides Ksf = 0.5,2.5, 36(dm/h) respectively (assumes halving of K for # but for tim's scenario, Ksf of loamy sand = 1.5 #loamy sand and sand over time) Ho, #distance from base to invert of outlet pipe (m) #if no outlet pipe set Ho to Hf + Hp Vstart, # Vol of water in system at t1 in L carea, carea.tank, # as in compile.inflow() prop.pave = 1#proportion of carea that is paving (as opposed to roof) ){ if(lined.bottom == 1 & lined.side == 0){ stop("If the sides are unlined, the bottom should be too", call. = FALSE) } if(isveg & medium == "gravel (scoria)"){ stop("Sorry: the calculator does not allow vegetated gravel systems", .call = FALSE) } Ksu <- ifelse(lined.bottom,0,0.001) #dm/h #cw 0.1 mm/h in Coffey report Ho <- Ho*10; Hf <- Hf*10; Hp <- Hp*10; Pf <- Pf*10#dm if(Ho < Hf & medium != "gravel (scoria)") Ho <- Ho + (Hf - Ho)*0.5 #revision 21 May to simulate field capacity of 20% for sand and loamy sand (i.e. *0.5 of 0.4 porosity) Qin <- inflow[,5] if(isveg){et <- et*Af*50} else {et <- rep(0,dim(inflow)[1])} #L/h #*50 revision 19 May 2008 #to simulate ET losses from known experimental raingardens (probably to do with #leaf area being much greater than Af, and ET draining interstitial soil after Vf = 0. Ksf <- c(1,2.5,36)[match(medium,c("loamy sand", "sand", "gravel (scoria)"))] #dm/h ################need to change back to 0.5 for loamy sand #set for Tim's scenarios - need to change back to 1.5 for loamy sand Porosity <- c(0.4,0.4,0.6)[match(medium,c("loamy sand", "sand", "gravel (scoria)"))] pondV <- Hp*(Af + Ap)*100/2 #total volume of pond in L filterV <- Porosity*Hf*Af*100 #total volume of filter in L suboV <- Porosity*Ho*Af*100 #volume of filter below outlet pipe in L #volume of water in each section at t0 Vp0 <- max(Vstart - filterV, 0)#L Vf0 <- min(filterV, Vstart) #L #depth of water in each section at t0 Hwf0 <- Hf*Vf0/filterV #dm Hwp0 <- ifelse(pondV == 0, 0, Hp*Vp0^(2/3)/pondV^(2/3)) #dm #Exfiltration from filter medium to u/lying soil: (= 0 if filter is lined or empty) #Area available for infiltration depends on Hwf if lined.side = FALSE Af0 <- ifelse(lined.side, Af*100, Af*100 + Hwf0*Pf*1.5) #sq dm ################assuming full Hwf available for infiltration Qexf1 <- ifelse(Vstart == 0 | lined.bottom, 0, Af0*Ksu) #L/h #cw changed from Vstart== 0 #Potential infiltration rate at t0 (any water falling on dry pond will be at Hf) Qinfi.pot <- ifelse(Ho == Hf + Hp & Vf0 == filterV, Qexf1, 100*Af*Ksf*max(Hf,Hwf0 + Hwp0)/Hf) #L/h #if the system has no underpipe, if(Ho == Hf + Hp){ #the volume that can flow from the pond into the filter is either limited by Qin, Qinf t1losses <- min(Vp0 + Qin[1], Qinfi.pot, #or by the available vol in the filter (+ exfiltration and et losses) filterV - Vf0 + Qexf1 + et[1]) #and Qo (flow out of the pipe in t1) is (obviously) zero Qo <- 0} else { #but if there is a pipe, space in filter is no longer a limiting factor t1losses <- min(Qin[1], Qinfi.pot) #and if t1losses > 0, outflow through pipe = t1losses - any void that was filled below pipe - Qexf - et Qo <- max(0, t1losses - max(0, suboV - Vf0) - Qexf1 - et[1]) } #set antecendent dry weather period as zero (assumes that the model starts immediately after the last rain) adwp <- ifelse(Qin[1] == 0, 1, 0) #N removal estimates (and adwp.trigger and dry.penalty are adjustments for adwp) adwp.trigger <- ifelse(isveg, 4*24, 5*24) ####assumes hourly timestep dry.penalty <- ifelse(isveg, 0.03, 0.02) if(isveg){ N1 <- c(0.59, 0.72, NA)[match(medium,c("loamy sand", "sand", "gravel (scoria)"))] } else { N1 <- c(-1.01, 0.32, 0.39)[match(medium,c("loamy sand", "sand", "gravel (scoria)"))] } if(medium == "loamy sand"){ Nred <- ifelse(adwp <= adwp.trigger, N1, N1 - dry.penalty*ceiling((adwp - adwp.trigger)/24)) }else{ Nred <- N1} Nout1 <- 2.2*Qo*(1-Nred)/1000000 ###kg of N in Qo Qover <- max(0, Vp0 + Qin[1] - t1losses - pondV) #L/h #Volume in pond at end of t1 Vp1 <- max(0, Vp0 + Qin[1] - t1losses - Qover) #L/h #Volume in filter at end of t1 #if everything is correct to now, there shouldn't be a need to allow for the #possibility of this being greater than filterV, should there? Vf1 <- max(0, Vf0 + t1losses - Qo - Qexf1 - et[1]) budget <- list(flowin = Qin, out = Qo, Nout = Nout1, over = Qover, store.pond = Vp1, store.filter = Vf1, Qexf = Qexf1) for(i in 2:dim(inflow)[1]){ if(i == 2) {Vfprev <- Vf1; Vpprev <- Vp1} else { Vfprev <- budget$store.filter[i-1] Vpprev <- budget$store.pond[i-1] } adwp <- ifelse(Qin[i] == 0, adwp + 1, 0) Hwfprev <- Hf*Vfprev/filterV #dm Hwpprev <- ifelse(pondV == 0, 0, Hp*Vpprev^(2/3)/pondV^(2/3)) #dm Afprev <- ifelse(lined.side, Af*100, Af*100 + Hwfprev*Pf*1.5) # sq dm budget$Qexf <- c(budget$Qexf, ifelse(Vfprev == 0 | lined.bottom, 0, Afprev*Ksu)) #L/h Qinfi.pot <- ifelse(Ho == Hf + Hp & Vfprev == filterV, budget$Qexf[i] + et[i], 100*Af*Ksf*max(Hf,Hwfprev + Hwpprev)/Hf) #L/h if(Ho == Hf + Hp){ tilosses <- min(Vpprev + Qin[i], Qinfi.pot, filterV - Vfprev + budget$Qexf[i] + et[i]) budget$out <- c(budget$out, 0)} else { tilosses <- min(Qin[i], Qinfi.pot) budget$out <- c(budget$out, max(0,tilosses - max(0, suboV - Vfprev) - budget$Qexf[i] - et[i])) } if(medium == "sandy loam"){ Nred <- ifelse(adwp <= adwp.trigger, N1, N1 - dry.penalty*ceiling((adwp - adwp.trigger)/24)) }else{ Nred <- N1} budget$Nout <- c(budget$Nout, 2.2*budget$out[i]*(1-Nred)/1000000) ###kg of N in budget$out[i] budget$over <- c(budget$over, max(0, Vpprev + Qin[i] - tilosses - pondV)) budget$store.pond <- c(budget$store.pond, max(0,Vpprev + Qin[i] - tilosses - budget$over[i])) #L/h budget$store.filter <- c(budget$store.filter, max(0, Vfprev + tilosses - budget$out[i] - budget$Qexf[i] - et[i])) } acceptable.out <- (carea + carea.tank)*0.1 #was 0.068 L/h/m2 in first draft; changed to 0.1 L/h/m2 (16 May 08) #We assume that an acceptable outflow from a bioretention system equals the #the infiltration rate of the surrounding soil mulitplied by the system's catchment area #As a comforting reality check, this is close to equivalent to the #depth of flow observed in a reference creek (Olinda Creek at Mt Evelyn) #at the peak of an event that is sufficient to commence overland flow: #about 6 m3/s for a 31.8 km2 catchment = 0.068 mm/h rofh <- budget$out > acceptable.out | budget$over > 0 drof <- aggregate(rofh, by = list(date = inflow$date), FUN = sum) roverh <- budget$over > 0 drover <- aggregate(roverh, by = list(date = inflow$date), FUN = sum) roexfh <- budget$out > acceptable.out droexf <- aggregate(roexfh, by = list(date = inflow$date), FUN = sum) days.runoff <- round(sum(drof$x > 0)*365/length(unique(inflow$date)),0) days.over <- round(sum(drover$x > 0)*365/length(unique(inflow$date)),0) days.oexf <- round(sum(droexf$x > 0)*365/length(unique(inflow$date)),0) # days.runoff <- round(apply(budget$overflow > 0,2,sum)*365/dim(dailyrain)[1],0) nat.rof <- 12 #see tankmodel.R urban.rof <- 121 #ditto urban.Nload <- 2.2*sum(hrunoff$runoff.mm)*(carea + carea.tank)*365*24/(10^6*length(hrunoff$runoff.mm)) garden.Nload <- (2.2*sum(budget$over)/10^6 + sum(budget$Nout))*365*24/(length(hrunoff$runoff.mm)) #assumes hourly timestep #runoff index calculated for entire catchment area, including that of the tank rof.index <- (1 - max(0,(days.runoff - nat.rof)/(urban.rof - nat.rof))) * (carea + carea.tank)/100 N.index <- (1- garden.Nload/urban.Nload)*(carea + carea.tank)/100 #weight paved areas 37.5:62.5 N and rof, and roof areas 30:50 (missing 20 is water savings) EBIt <- (prop.pave*(0.625*rof.index + 0.375*N.index) + (1 - prop.pave)*(0.5*rof.index + 0.3*N.index)) list(budget = budget, EBstats = data.frame(Nkg.reduction = round(urban.Nload - garden.Nload,3), days.runoff = days.runoff, days.over = days.over, days.oexf = days.oexf, stream.index = rof.index, N.index = N.index, EBI = round(EBIt, 2))) }