---
author: "Michael Koohafkan"
title: "Reproducing the R Journal Publication"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Publication Reproduction}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r include = FALSE}
library(knitr)
opts_chunk$set(echo = TRUE, eval = FALSE)
```

This document provides the code needed to reproduce all results and figures 
from 

> Michael C. Koohafkan and Bassam A. Younis (2015). 
> Open-Channel Computation with R. The R Journal, 7(2), 249-262. 
> URL https://journal.r-project.org/archive/2015-2/koohafkan-younis.pdf

The .Rnw file used to generate the full publication is hosted on the
[package development repository](https://github.com/mkoohafkan/rivr/tree/master/inst/pub).

```{r the-code}
## ---- load-libs -------------------------------------------------------
library(knitr)
opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = "fig/", 
  warning = FALSE, message = FALSE, error = FALSE, results = "asis", 
  out.width = '\\textwidth', cache = TRUE, cache.path = "cache/")
library(rivr)
library(ggplot2)
library(scales)
library(RColorBrewer)
library(dplyr)
library(xtable)

## ---- startparams -----------------------------------------------------
plotopts = list(theme_bw(), xlab(expression(Distance~from~control~section~~(ft))))
g = 32.2
Cm = 1.486
slope = 0.001
mannings = 0.045
flow = 250  
width = 100
sideslope = 0
# calculate control depth as 1ft above the normal depth
depth.m1 = round(1 + normal_depth(slope, mannings, flow, 2, Cm, 
  width, sideslope), 2) 
depth.m2 = round(1.1*critical_depth(flow, 1, g, width, sideslope), 2)
rivdist = 3000

## ---- restest ---------------------------------------------------------
# Test sensitivity of step size
stepsizes = c(500, 100, 50, 10, 1)
resolution.test.m1 = list()
resolution.test.m2 = list()
for(r in stepsizes){
  resolution.test.m1[[paste0('dx=', r)]] = compute_profile(slope, mannings, 
    flow, depth.m1, Cm, g, width, sideslope, stepdist = r, 
    totaldist = rivdist)
  resolution.test.m2[[paste0('dx=', r)]] = compute_profile(slope, mannings, 
    flow, depth.m2, Cm, g, width, sideslope, stepdist = r, 
    totaldist = rivdist)
}
resolution.plot = NULL
for(lbl in names(resolution.test.m1)){
  f1 = resolution.test.m1[[lbl]]
  f2 = resolution.test.m2[[lbl]]
  f1['type'] = 'M1'
  f2['type'] = 'M2'
  f = rbind(f1, f2)
  f['run'] = lbl
  f['res'] = as.numeric(substr(lbl, 4, 10L))
  resolution.plot = rbind(resolution.plot, f)
}
resolution.plot['res'] = factor(resolution.plot$res, levels = rev(sort(stepsizes)))  
resolution.plot['run'] = factor(resolution.plot$run, 
  levels = names(resolution.test.m1))
resolution.plot['type'] = factor(resolution.plot$type, levels = c('M1', 'M2'))
stepsize = 10

## ----resfig ----------------------------------------------------------
ggplot(resolution.plot, aes(x = x, y = y + z)) + plotopts + 
geom_line(aes(color = res, linetype = res), size = 1) + 
facet_wrap(~type) + ylab(expression(River~stage~~(ft))) + 
scale_linetype_manual("resolution (ft)", labels = levels(resolution.plot$res), 
  values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +  
scale_color_manual("resolution (ft)", values = brewer.pal(7, "YlGnBu")[3:7], 
  labels = levels(resolution.plot$res))

## ----roughness -------------------------------------------------------
roughness.test.m1 = list()
roughness.test.m2 = list()
roughnesses = 0.045*seq(0.5, 1.5, length = 5)
for(n in roughnesses){
  roughness.test.m1[[paste0('n=',n)]] = compute_profile(slope, n, flow, 
    depth.m1, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
  roughness.test.m2[[paste0('n=',n)]] = compute_profile(slope, n, flow, 
    depth.m2, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
}
roughness.plot = NULL
for(lbl in names(roughness.test.m1)){
  f1 = roughness.test.m1[[lbl]]
  f2 = roughness.test.m2[[lbl]]
  f1['type'] = 'M1'
  f2['type'] = 'M2'
  f = rbind(f1, f2)
  f['run'] = lbl
  f['mannings'] = as.numeric(substr(lbl, 4, 10L))
  roughness.plot = rbind(roughness.plot, f)
}
roughness.plot['run'] = factor(roughness.plot$run, 
  levels = names(roughness.test.m1))
roughness.plot['mannings'] = factor(roughness.plot$mannings, 
  levels = sort(roughnesses))
roughness.plot['type'] = factor(roughness.plot$type, levels=c('M1', 'M2'))

## ---- roughplot -------------------------------------------------------
ggplot(roughness.plot, aes(x = x, y = y + z)) + 
geom_line(aes(linetype = mannings, color = mannings), size = 1) +  plotopts + 
facet_wrap(~ type) + ylab(expression(River~stage~~(ft))) + 
scale_linetype_manual("Bed roughness", labels = levels(roughness.plot$mannings), 
  values = c("solid", "longdash", "dashed", "dotdash", "dotted")) + 
scale_color_manual("Bed roughness", labels = levels(roughness.plot$mannings), 
  values = brewer.pal(6, "Oranges")[2:6])

## ---- relroughplot ----------------------------------------------------
roughness.rel.m1 = list()
roughness.rel.m2 = list()
for(n in roughnesses){
  thisyn = normal_depth(slope, n, flow, 2, Cm, width, sideslope)
  roughness.rel.m1[[paste0('n=',n)]] = compute_profile(slope, n, flow, 
    1.25*thisyn, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
  roughness.rel.m2[[paste0('n=',n)]] = compute_profile(slope, n, flow, 
    0.75*thisyn, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
  roughness.rel.m1[[paste0('n=',n)]]['pd.yn'] = 
    (roughness.rel.m1[[paste0('n=',n)]]$y - thisyn)/thisyn
  roughness.rel.m2[[paste0('n=',n)]]['pd.yn'] = 
    (roughness.rel.m2[[paste0('n=',n)]]$y - thisyn)/thisyn
}
roughness.rel.plot = NULL
for(lbl in names(roughness.rel.m1)){
  f1 = roughness.rel.m1[[lbl]]
  f2 = roughness.rel.m2[[lbl]]
  f1['type'] = 'M1'
  f2['type'] = 'M2'
  f = rbind(f1, f2)
  f['run'] = lbl
  f['mannings'] = as.numeric(substr(lbl, 4, 10L))
  roughness.rel.plot = rbind(roughness.rel.plot, f)
}
roughness.rel.plot['mannings'] = factor(roughness.plot$mannings, 
  levels = sort(roughnesses))
roughness.rel.plot['run'] = factor(roughness.rel.plot$run, 
  levels = names(roughness.rel.m1))
roughness.rel.plot['type'] = factor(roughness.rel.plot$type, 
  levels = c('M1', 'M2'))
ggplot(roughness.rel.plot, aes(x = x, y = pd.yn)) + 
geom_line(aes(linetype = mannings, color = mannings), size = 1) + 
plotopts + facet_wrap(~type, scales = 'free_y') + 
scale_y_continuous(expression(Percent~difference~from~normal~depth), labels = percent) + 
scale_linetype_manual("Bed roughness", labels = levels(roughness.rel.plot$mannings), 
  values = c("solid", "longdash", "dashed", "dotdash", "dotted")) + 
scale_color_manual("Bed roughness", labels = levels(roughness.rel.plot$mannings), 
  values = brewer.pal(6, "Oranges")[2:6])


## ---- loadup-kwm ------------------------------------------------------
oldscipen = options('scipen')
options(scipen = 1000)

slope = 0.001
extent = 150000
mannings = 0.045
B = 100
SS = 0
Cm = 1.486
g = 32.2
iflow = 250

# keep Courant number at 0.7 to balance temporal and spatial resolution
idepth = normal_depth(slope, mannings, 250, 10, 1.485, 100, SS)  
iarea = channel_geom(idepth, B, SS)["A"]
cn = 0.7*iarea/iflow
# define upstream boundary condition assuming a timestep in seconds
bcfunc = function(x)
  ifelse(x < 9000, 250 + (750/pi)*(1 - cos(pi*x/(60*75))), 250)
bctime = 76000
xnodes = extent/c(25, 50, 125, 250, 500, 1000, 2000, 5000) + 1
myxres = extent/(xnodes - 1)
mytres = cn*myxres
plotopts = list(theme_bw())

# run the model
modtime = list()
modelresults = list()
for(i in seq(length(xnodes))){
  numnodes = xnodes[i]
  xstep = myxres[i]
  tstep = mytres[i]
  bc = bcfunc(seq(0, bctime, by=tstep))
  mp = c(1, as.integer(50000/xstep + 1), as.integer(100000/xstep + 1), numnodes)
  mt = as.integer(round(seq(1, length(bc), length.out=10)))[seq(7)]
  if(xstep == min(myxres)){
    mt = as.integer(round(seq(1, length(bc), length.out=125)))
    mtslice = mt[round(seq(1, length(mt), length.out=7))]
  }
  lbl = paste0('dx=',xstep)
  modtime[[lbl]] <- system.time(modelresults[[lbl]] <- 
    route_wave(slope, mannings, Cm, g, B, SS, iflow, bc, 
	timestep=tstep, spacestep=xstep, numnodes=numnodes, 
	  monitor.nodes=mp, monitor.times=mt, engine="Kinematic"))
  modelresults[[lbl]]['deltax'] = xstep
  modelresults[[lbl]]['deltat'] = tstep  
  modelresults[[lbl]]['computationtime'] = modtime[[lbl]][[3]]
  modelresults[[lbl]]['label'] = lbl
}
allresults = do.call(rbind.data.frame, modelresults)
row.names(allresults) = NULL
names(allresults)[which(names(allresults)=='flow')] = "Q"
names(allresults)[which(names(allresults)=='distance')] = "x"
names(allresults)[which(names(allresults)=='time')] = "t"

## ---- floodhydrograph -------------------------------------------------
floodhydro = data.frame(time = seq(15000), flow=bcfunc(seq(15000)))
ggplot(floodhydro, aes(x=time/60, y=flow)) + geom_line(color="darkblue", 
  size = 1) + plotopts +
xlab(expression(time~~(minutes))) + ylab(expression(flow~~(ft^3~~s^-1)))

## ---- routeplot-kinematic ---------------------------------------------
# still plot
surfaces = group_by(filter(allresults, monitor.type=='timestep', 
  deltax==min(deltax)))
ggplot(filter(surfaces, step %in% mtslice[2:5]), 
  aes(x=x, y=Q, linetype=factor(round(t/60)), color=factor(round(t/60)))) + 
geom_line(size = 1) + plotopts + 
scale_color_manual("time\n(minutes)", values = brewer.pal(7, "YlGnBu")[3:6]) + 
scale_linetype_manual("time\n(minutes)", values = c("solid", 
  "dashed", "dotdash", "dotted")) + 
xlab(expression(distance~~downstream~~(ft))) + 
ylab(expression(flow~~(ft^3~~s^-1)))
# animation
routeylims = c(min(surfaces$Q), max(surfaces$Q))
routexlims = c(min(surfaces$x), max(surfaces$x))
thisdt = unique(surfaces$deltat)
for(n in unique(surfaces$step))
  print(
    ggplot(filter(surfaces, step==n), aes(x=x, y=Q)) + 
    geom_line(color = "darkblue", size = 1) + plotopts + 
    scale_linetype_discrete("time\n(minutes)") + 
    scale_x_continuous(expression(distance~~downstream~~(ft)), 
      limits=routexlims) + 
	scale_y_continuous(expression(flow~~(ft^3~~s^-1)), limits=routeylims) + 
	ggtitle(paste(format(round((n - 1)*thisdt/60), width = 4), "minutes"))
  )

## ---- xsectionplot-kinematic ------------------------------------------
xsections = filter(allresults, monitor.type == "node", 
  x %in% c(0,50000))[c("x","t","Q", "deltax")]
xsections['deltax'] = factor(paste(xsections$deltax, "feet"), 
  levels=paste(sort(unique(xsections$deltax)), "feet"))
ggplot(xsections, aes(x=t/60, y=Q, linetype=factor(format(x, big.mark=",", 
  trim=TRUE), levels = format(sort(unique(x)), big.mark=",", trim=TRUE)), 
  color = deltax)) + 
geom_path(size = 1) + facet_grid(deltax~.) + plotopts +
scale_y_continuous(expression(flow~~(ft^3~~s^-1)), breaks=c(300, 500,700)) + 
scale_linetype_manual("downstream\ndistance (ft)", values = c("solid", 
  "dashed", "dotted")) + 
scale_x_continuous(expression(time~~(min)), limits=c(0,600)) + 
scale_color_manual(values = c("#7570B3", "#8DA0CB", "#66C2A5", "#A6D854",
  "#E5C494", "#E6AB02", "#FC8D62", "#E78AC3"), guide = FALSE)  

## ---- peakstable-kinematic --------------------------------------------
peaks = as.data.frame(summarize(group_by(filter(allresults, x==50000, 
  monitor.type=='node'), computationtime, deltax, deltat, monitor.type), 
  peak.flow=max(Q), time.to.peak=t[which(Q==peak.flow)]))
real.peak = max(floodhydro$flow)
peaks['flow.percent.error'] = (peaks$peak.flow - real.peak)/real.peak
tbld = as.data.frame(peaks[c("deltax", "deltat", "flow.percent.error", 
  "computationtime")])
tbld["deltax"] = round(tbld$deltax)
tbld["deltat"] = round(tbld$deltat, 2)
tbld["flow.percent.error"] = round(100*tbld$flow.percent.error, 2)
names(tbld) = c("$\\Delta x$ (ft)", "$\\Delta t$ (s)", "\\% error (peak flow)", 
  "cost (s)")
ptable = xtable(tbld[order(tbld[,1]),], digits = c(0,0,2,2,2))
print(ptable, include.rownames = FALSE, sanitize.colnames.function = identity,
  sanitize.text.function = function(x){x}, floating=FALSE, hline.after=NULL, 
  add.to.row=list(pos=list(-1,0, nrow(ptable)),
  command=c('\\toprule ', '\\midrule ', '\\bottomrule ')), 
  format.args = list(big.mark = ","))  

## ---- loadup-dwm ------------------------------------------------------
oldscipen = options('scipen')
options(scipen = 1000)

slope = 0.001
extent = 150000
mannings = 0.045
B = 100
SS = 0
g = 32.2
Cm = 1.486
iflow = 250

# define upstream boundary condition assuming a timestep in seconds
idepth = normal_depth(slope, mannings, 250, 10, Cm, 100, SS)  
iarea = channel_geom(idepth, B, SS)[["A"]]
# keep Courant number at 0.06 to balance temporal and spatial resolution
cn = 0.06*iarea/250
bcfunc = function(x)
  ifelse(x < 9000, 250 + (750/pi)*(1 - cos(pi*x/(60*75))), 250)
bctime = 76000
xnodes = extent/c(50, 125, 250, 500, 1000) + 1
myxres = extent/(xnodes - 1)
mytres = cn*myxres
plotopts = list(theme_bw())

# run the model
modtime = list()
modelresults = list()
for(i in seq(length(xnodes))){
  numnodes = xnodes[i]
  xstep = myxres[i]
  tstep = mytres[i]
  bc = bcfunc(seq(0, bctime, by=tstep))
  dc = rep(-1, length(bc))
  mp = c(1, as.integer(c(1000, 50000, 100000, 149000)/xstep + 1), numnodes)
  mt = as.integer(round(seq(1, length(bc), length.out=10)))[seq(7)]
  if(xstep == min(myxres)){
	mt = as.integer(round(seq(1, length(bc), length.out=125)))	
    mtslice = mt[round(seq(1, length(mt), length.out=7))]
  }
  lbl = paste0('dx=',xstep)
  modtime[[lbl]] <- system.time(modelresults[[lbl]] <- 
    route_wave(slope, mannings, Cm, g, B, SS, iflow, bc, dc, timestep=tstep, 
      spacestep=xstep, numnodes=numnodes, monitor.nodes=mp, monitor.times=mt, 
	  engine="Dynamic", boundary.type="QQ"))
  modelresults[[lbl]]['deltax'] = xstep
  modelresults[[lbl]]['deltat'] = tstep  
  modelresults[[lbl]]['computationtime'] = modtime[[lbl]][[3]]
  modelresults[[lbl]]['label'] = lbl
}
# combine results
allresults = do.call(rbind.data.frame, modelresults)
row.names(allresults) = NULL
names(allresults)[which(names(allresults)=='flow')] = "Q"
names(allresults)[which(names(allresults)=='distance')] = "x"
names(allresults)[which(names(allresults)=='time')] = "t"
allresults['t'] = (allresults$step - 1)*allresults$deltat
allresults['x'] = (allresults$node - 1)*allresults$deltax

## ---- routeplot-characteristic ----------------------------------------
# still plot
surfaces = filter(allresults, monitor.type=='timestep', 
  deltax==min(deltax))
ggplot(filter(surfaces, step %in% mtslice[2:5]), 
  aes(x=x, y=Q, linetype=factor(round(t/60)), color=factor(round(t/60)))) + 
geom_line(size = 1) + plotopts + 
scale_color_manual("time\n(minutes)", values = brewer.pal(7, "YlGnBu")[3:6]) + 
scale_linetype_manual("time\n(minutes)", values = c("solid", 
  "dashed", "dotdash", "dotted")) +
  xlab(expression(distance~~downstream~~(ft))) + 
  ylab(expression(flow~~(ft^3~~s^-1)))
# animation
routeylims = c(min(surfaces$Q), max(surfaces$Q))
routexlims = c(min(surfaces$x), max(surfaces$x))
thisdt = unique(surfaces$deltat)
for(n in unique(surfaces$step))
  print(
    ggplot(filter(surfaces, step==n), aes(x=x, y=Q)) + 
      geom_line(color = "darkblue", size = 1) + plotopts + 
      scale_linetype_discrete("time\n(minutes)") + 
      scale_x_continuous(expression(distance~~downstream~~(ft)), 
        limits=routexlims) + 
	  scale_y_continuous(expression(flow~~(ft^3~~s^-1)), limits = routeylims) + 
	  ggtitle(paste(format(round((n - 1)*thisdt/60), width = 4), "minutes"))
  )

## ---- xsectionplot-characteristic -------------------------------------
xsections = filter(allresults, monitor.type == "node", 
  x %in% c(0, 50000, 150000))[c("x","t","Q", "deltax")]
xsections['deltax'] = factor(paste(xsections$deltax, "feet"), 
  levels=paste(sort(unique(xsections$deltax)), "feet"))
data(waterolympics)
realdat = NULL
for(x in unique(xsections$deltax)){
  thisdat = waterolympics[waterolympics$t > 200*60,]
  thisdat['deltax'] = x
  realdat = rbind(realdat, thisdat)
}
realdat['deltax'] = factor(realdat$deltax)
ggplot(xsections, aes(x=t/60, y=Q, linetype=factor(format(x, big.mark=",", 
  trim=TRUE), levels = format(sort(unique(x)), big.mark=",", trim=TRUE)), 
  color = deltax)) + 
geom_path(size = 1) + geom_point(data = realdat[realdat$t > 200*60,], color = "black") + 
facet_grid(deltax~.) +scale_y_continuous(expression(flow~~(ft^3~~s^-1)), 
  breaks=c(300, 500,700)) + 
scale_linetype_manual("downstream\ndistance (ft)", values = c("solid", 
  "dashed", "dotted")) + plotopts +
scale_x_continuous(expression(time~~(min)), limits=c(0,1250)) +
scale_color_manual(values = c("#8DA0CB", "#66C2A5", "#A6D854", "#E5C494", 
  "#E6AB02"), guide = FALSE)  

## ---- peakstable-characteristic ---------------------------------------
peaks = as.data.frame(summarize(group_by(filter(allresults, x==50000, 
  monitor.type=='node'), computationtime, deltax, deltat, monitor.type), 
  peak.flow=max(Q), time.to.peak=t[which(Q==peak.flow)]))
real.peak.flow = max(waterolympics$Q)
real.time.to.peak = waterolympics$t[which(waterolympics$Q==real.peak.flow)][1]
peaks['flow.percent.error'] = 
  (peaks$peak.flow - real.peak.flow)/real.peak.flow
peaks['time.percent.error'] = 
  (peaks$time.to.peak - real.time.to.peak)/real.time.to.peak
tbld = as.data.frame(peaks[c("deltax", "deltat", "flow.percent.error", 
  "time.percent.error", "computationtime")])
tbld["flow.percent.error"] = round(100*tbld$flow.percent.error, 2)
tbld["time.percent.error"] = round(100*tbld$time.percent.error, 2)
names(tbld) = c("$\\Delta x$ (ft)", "$\\Delta t$ (s)", "\\% error (peak flow)", 
  "\\% error (time to peak)", "cost (s)")
ptable = xtable(tbld[order(tbld[,1]),], digits = c(0,0,2,2,2,2))
print(ptable, include.rownames=FALSE, sanitize.colnames.function=identity,
  sanitize.text.function = function(x){x}, floating=FALSE, hline.after=NULL, 
  add.to.row=list(pos=list(-1,0, nrow(ptable)),  
  command=c('\\toprule ', '\\midrule ', '\\bottomrule ')), 
  format.args = list(big.mark = ","))  

## ---- loadup-boundaries -----------------------------------------------
oldscipen = options('scipen')
options(scipen = 1000)

plotopts = list(theme_bw())
slope = 0.00008
extent = 5000
mannings = 0.013
B = 6.1
SS = 1.5
g = 9.81
Cm = 1

ic = 126
id = 5.79
ia = channel_geom(id, B, SS)[["A"]]
CN = 0.9

dx = 10
dt = round(dx*CN/(ic/ia + sqrt(id*g)), 2)
numnodes = extent/dx + 1

bctime = 2000
bc = rep(id, round(bctime/dt) + 1)
dc = rep(0, length(bc))
dt = round(bctime/(length(bc) - 1), 2)

CN = dt*(ic/ia + sqrt(id*g))/dx

mp = c(1, as.integer(c(1500, 2500, 3000, 5000)/dx + 1))
mt = as.integer(round(seq(0, length(bc)-1, by=25))) + 1L
mtslice = c(1, as.integer(c(500, 1000, 1500, 2000)/dt + 1))

dclose.lax = route_wave(slope, mannings, Cm, g, B, SS, ic, bc, dc, 
  timestep=dt, spacestep=dx, numnodes=numnodes, monitor.nodes=mp, 
  monitor.times=mt, engine="Dynamic", scheme="Lax", boundary.type="yQ")
dclose.mac = route_wave(slope, mannings, Cm, g, B, SS, ic, bc, dc, 
  timestep=dt, spacestep=dx, numnodes=numnodes, monitor.nodes=mp, 
  monitor.times=mt, engine="Dynamic", scheme="MacCormack", boundary.type="yQ")
dclose.lax['scheme'] = "Lax diffusive"
dclose.mac['scheme'] = "MacCormack predictor-corrector"
dclose = rbind(dclose.lax, dclose.mac)
row.names(dclose) = NULL
dclose["CN"] = dt*(dclose$flow/dclose$area + sqrt(dclose$depth*g))/dx

## ---- through-time ----------------------------------------------------
dclose.times = filter(dclose, monitor.type=="timestep")
ggplot(filter(dclose.times, step %in% mtslice), 
  aes(x=distance, y=depth, linetype = factor(round(time)), 
    color = factor(round(time)))) + geom_path(size = 1) + 
scale_y_continuous(expression(depth~~(ft))) + facet_wrap(~scheme) +
scale_x_continuous(expression(distance~~downstream~~(ft))) + plotopts +
scale_color_manual("time (s)", values = brewer.pal(7, "YlGnBu")[3:7]) + 
scale_linetype_manual("time (s)", values = c("solid", "longdash", "dashed", 
  "dotdash", "dotted"))
# animation
routeylims = c(min(dclose.times$depth), max(dclose.times$depth))
routexlims = c(min(dclose.times$distance), max(dclose.times$distance))
for(n in unique(dclose.times$step))
  print(
    ggplot(filter(dclose.times, step==n), aes(x=distance, y=depth)) + 
    geom_line(color = "darkblue", size = 1) + plotopts + 
    scale_linetype_discrete("time\n(minutes)") + 
    scale_x_continuous(expression(distance~~downstream~~(ft)), 
      limits=routexlims) + 
	scale_y_continuous(expression(depth~~(ft)), limits=routeylims) + 
	ggtitle(paste(format(round((n - 1)*dt), width = 4), "seconds")) + 
    facet_wrap(~scheme)
  )
```