In this exercise you will:

0. Load libraries.

library(dataRetrieval)
library(hydroTSM)

I. Flow duration curves

A. Load USGS streamflow data for the lower San Diego River at Mast Blvd (siteNumber 11022480)

See Exercise 4, II.A.1-4 (I’m not showing the code). Don’t forget to rename the Q and QA columns.
I called the discharge timeseries downloaded from the USGS “x”, as in Exercise 4.

B. Summarize the data availability by water year from the beginning to the end of the discharge time series.

# x is the discharge timeseries from the USGS.
breaks <- seq(from=as.Date("1912-10-01"), to=as.Date("2015-10-01"), by="year")  
years.breaks = as.numeric(format(breaks,"%Y"))
labels.wy = years.breaks[2:length(breaks)]
x$WaterYear <- cut(x$Date, breaks,labels=labels.wy)
data.avail.by.wy = aggregate(x$Q.ft.s,by=list(x$WaterYear),FUN=length)
data.avail.by.wy$WYear = as.numeric(as.character(data.avail.by.wy$Group)) # Get the water year as a number

C. Create flow duration curves for different time periods.

We will divide your time series into three groups for analysis:

  1. 1914-1935 (pre El Capitan Dam)

  2. 1945-1965 (post El Capitan and San Vicente dams but pre-urbanization)

  3. 1995-2015 (post dams and post-urbanization)

1. Find water years with data for all days in the water year.
goodyears = data.avail.by.wy[data.avail.by.wy$x>=365,"WYear"]
2. Intersect those water years of good data with the year ranges you want to subset by:
wyears.1914.1935.good = intersect(seq(1914,1935),goodyears)
wyears.1945.1965.good = intersect(seq(1945,1965),goodyears)
wyears.1995.2015.good = intersect(seq(1995,2015),goodyears)

Qsub.1914.1935 = x[x$WaterYear %in% wyears.1914.1935.good,]
Qsub.1945.1965 = x[x$WaterYear %in% wyears.1945.1965.good,]
Qsub.1995.2015 = x[x$WaterYear %in% wyears.1995.2015.good,]

#  Do head() and tail() on each to make sure it got the ranges correct.
3. Find the range of all dates, and plot discharge for the three intervals.
begin.end = range(x$Date)
plot(Qsub.1914.1935$Date,Qsub.1914.1935$Q.ft.s,xlim=begin.end,ylab="Q ft3/s",type="l")
lines(Qsub.1945.1965$Date,Qsub.1945.1965$Q.ft.s,col="red")
lines(Qsub.1995.2015$Date,Qsub.1995.2015$Q.ft.s,col="blue")
legend("topright",c("1914-1935","1945-1965","1995-2015"),col=c("black","red","blue"),lty=c(1,1,1))

D. Flow duration curves.

1. Develop flow-duration curves for each period using the “fdc” function, and plot together on the same graph.

Look at the help file for fdc:

?fdc

Then use it to plot flow duration curves for three time periods.

fdc.1914.1935 = fdc(Qsub.1914.1935$Q.ft.s,new=TRUE,ylab="Q ft3/s")
fdc.1945.1965 = fdc(Qsub.1945.1965$Q.ft.s,new=FALSE,col="red")
fdc.1995.2015 = fdc(Qsub.1995.2015$Q.ft.s,new=FALSE,col="blue")
legend("topright",c("1914-1935","1945-1965","1995-2015"),col=c("black","red","blue"),lty=c(1,1,1))

What is your interpretation of these FDCs?

II. Stormflow

A. Load your DAYMET precipitation series

outdir = << Y directory, "rainfall/DAYMET/" >>
daymet.lowerSDR = read.table(paste0(outdir,"DAYMET_lowerSDR.txt"))
daymet.lowerSDR$Date = as.Date(daymet.lowerSDR$Date,format="%Y-%m-%d")

B. Select a single year to find rainfall events.

plot(daymet.lowerSDR$Date,daymet.lowerSDR$Precip.mm,type="l")

Looks like the 2011 water year was pretty wet and had several events. Zoom in on that year, and plot both rainfall and discharge on a single plot. This is how hydrologists often plot rainfall and runoff:

daymet.2011wy = daymet.lowerSDR[daymet.lowerSDR$Date %in% seq(as.Date("2010-10-01"),as.Date("2011-09-30"),by="day"),]

par(mar=c(5,5,2,5))  # Sets the margins of the plot.
barplot(daymet.2011wy$Precip.mm,daymet.2011wy$Date,space=c(0,1),width=0.5,ylim=rev(c(0,4*max(daymet.2011wy$Precip.mm))),ylab="Rainfall, mm",col="blue",border="blue",xaxt="n")

par(new=T)
plot(x$Date,x$Q.ft.s,type="l",col="red",yaxt="n",ylim=c(0,5000),ylab="",xlim=as.Date(c("2010-10-01","2011-09-30")),xlab="")
axis(side=4)  #  Add and axis to the right side for discharge.
mtext("Streamflow, ft3/s",side=4,line=3)  # Adds text to the second y-axis.

C. Find individual rainfall events.

rain.over.0 = daymet.2011wy[daymet.2011wy$Precip.mm>0,]
rainindex = c(0,cumsum(diff(rain.over.0$Date)>1)) # The 1 specifies that the difference in date between days with rainfall must be greater than 1 to count as a separate event.
#  How many events were identified?

D. Aggregate rainfall and runoff by rainfall event.

rain.total.by.event = aggregate(rain.over.0$Precip.mm,by=list(rainindex),FUN="sum")
ndays.in.event = aggregate(rain.over.0$Precip.mm,by=list(rainindex),FUN="length")
  
#  Extract discharge for just 2011 water year:
x2011 = x[(x$Date>=as.Date("2010-10-01")) & x$Date<=as.Date("2011-09-30"),]

#  And for days where rainfall was greater than zero
  # What assumption are we making about the timing of runoff, in order to compare the daily rainfall and runoff series?
x2011.rain.over.0 = x2011[daymet.2011wy$Precip.mm>0,]
x.mean.by.event.cfs = aggregate(x2011.rain.over.0$Q.ft.s,by=list(rainindex),FUN="mean")

E. Convert mean discharge in ft3/s to mm for each storm event.

#  1ft3    1 m3             1           1 mi2          1000mm     3600x24 s     ndays.in.event
#  ___  x  ____         x  ________  x  _____      x   ______  x  _______    x  _____
#   s     35.3157 ft3       A mi2       1609^2 m2        m        day           event
1. Note that you need to find the drainage area of the watershed.

From exercise 4:

site.data.sdr.mast = readNWISsite("11022480")
drain.area.mi2.sdr.mast = site.data.sdr.mast$drain_area_va
site.data.svicente = readNWISsite("11022100")
drain.area.mi2.svicente = site.data.svicente$drain_area_va
site.data.elcap = readNWISsite("11020600")
drain.area.mi2.elcap = site.data.elcap$drain_area_va
drainarea.sdr.mast.noup.mi2 = drain.area.mi2.sdr.mast-drain.area.mi2.svicente-drain.area.mi2.elcap
2. Now convert the units:
x.mm.by.event = (x.mean.by.event.cfs$x)*(1000*3600*24)/(35.3157*1609.34*1609.34*drainarea.sdr.mast.noup.mi2)*ndays.in.event$x

F. Plot storm event runoff vs rainfall:

plot(rain.total.by.event$x,x.mm.by.event,xlab="Rainfall, mm",ylab="Runoff, mm",pch=19,cex=1.4)

I like to use big dots…“cex”" adjusts the dot size.

You’ll intepret this rainfall-runoff relationship in the homework.