In this exercise you will:

0. Load libraries.

library(dataRetrieval)
library(hydroTSM)

I. Flow duration curves

A. Load USGS daily streamflow data for Los Penasquitos

See Exercise 5 (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 “Qin” (instead of “x” as in Exercise 5).

B. Create flow duration curves for different time periods.

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

  1. 1965-1989 (pre urban)

  2. 1990-2005 mid-urban

  3. 2006-2019 full urban

1. Create subsets of the full dataset

I’m showing the code for 1965-1989…you will also do the code for 1990-2005 and 2005-2015 on your own. If you call the time series for e.g. 1990-2005 “Q.1990.2005” it will work with the code presented later…

Q.1965.1989 = Qin[(Qin$Date >= as.Date("1964-10-01")) & (Qin$Date <= as.Date("1989-09-30")),]
3. Find the range of all dates, and plot discharge for the three intervals.

Note: my html rendition of this plot is missing 2006-2010 and not rending properly…your plot should show a continuous time series.

begin.end = range(Qin$Date)  # Get the date range for the whole dataset
plot(as.Date(Q.1965.1989$Date),Q.1965.1989$Q.ft.s,xlim=as.Date(begin.end),ylab="Q ft3/s",type="l",xlab="")
lines(Q.1990.2005$Date,Q.1990.2005$Q.ft.s,col="red")
lines(Q.2006.2019$Date,Q.2006.2019$Q.ft.s,col="blue")
legend("topright",c("1965-1989","1990-2005","2006-2019"),col=c("black","red","blue"),lty=c(1,1,1))

The legend cuts off some of the later time series. Figure out how to change the ylimits so that the legend doesn’t cover up the time series.

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 fdc to plot flow duration curves for three time periods.

fdc.1914.1935 = fdc(Q.1965.1989$Q.ft.s,new=TRUE,ylab="Q ft3/s")
fdc.1945.1965 = fdc(Q.1990.2005$Q.ft.s,new=FALSE,col="red")
fdc.1995.2019 = fdc(Q.2006.2019$Q.ft.s,new=FALSE,col="blue")
legend("topright",c("1965-1989","1990-2005","2006-2019"),col=c("black","red","blue"),lty=c(1,1,1))

Make the graph prettier:

# First determine what your y-limits (ylim) should be:
print(range(Qin$Q.ft.s))
## [1]    0 2270
fdc.1965.1989 = fdc(Q.1965.1989$Q.ft.s,new=TRUE,ylab="Q ft3/s",las=1, cex.lab=1, cex.axis=0.8,yat=c(0.01,0.1,1,1,10,100,1000,3000),ylim=c(0.01,3000))
fdc.1990.2005 = fdc(Q.1990.2005$Q.ft.s,new=FALSE,col="red", thr.shw=FALSE,xat=NA,yat=NA,ylim=c(0.01,3000),cex.axis=0.8)
fdc.2006.2019 = fdc(Q.2006.2019$Q.ft.s,new=FALSE,col="blue", thr.shw=FALSE, xat=NA,yat=NA,ylim=c(0.01,3000),cex.axis=0.8)
legend("topright",c("1965-1989","1990-2005","2006-2019"),col=c("black","red","blue"),lty=c(1,1,1))

You’ll interpret these FDCs in the homework.

II. Stormflow

A. Load your daily precipitation series. See Exercise 5.

  • Note that when you saved your rainfall time series, it may have changed the format of the date, for example, to “2010-10-30” instead of “10/30/2010”, so you may need to change the format in the as.Date line…

B. Select a single year to find rainfall events (doing all years would be too much).

Based on the discharge plot, it looks like the 1998 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:

rainfall.1998wy = daily.rainfall[daily.rainfall$Date %in% seq(as.Date("1997-10-01"),as.Date("1998-09-30"),by="day"),]

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

par(new=T)
plot(Qin$Date,Qin$Q.ft.s,type="l",col="red",yaxt="n",ylim=c(0,1500),ylab="",xlim=as.Date(c("1997-10-01","1998-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. This uses the “cumsum” function and the “diff” function (use ? to see what they are).

rain.over.0 = rainfall.1998wy[rainfall.1998wy$Precip.inches>0,]
rainindex = c(0,cumsum(diff(rain.over.0$Date)>1)) 
  # Rain index labels all days with rainfall > 0 with an event number, starting with 0.
# 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.inches,by=list(rainindex),FUN="sum")
ndays.in.event = aggregate(rain.over.0$Precip.inches,by=list(rainindex),FUN="length")
  
#  Extract discharge for just 2011 water year:
Q1998 = Qin[(Qin$Date>=as.Date("1997-10-01")) & Qin$Date<=as.Date("1998-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?
Q1998.rain.over.0 = Q1998[rainfall.1998wy$Precip.inches>0,]
Q.mean.by.event.cfs = aggregate(Q1998.rain.over.0$Q.ft.s,by=list(rainindex),FUN="mean")

E. Convert mean discharge in ft3/s to mm for each storm event. I’m not showing the code for this…

To get the total mm of runoff for the storm, remember that the mean is the mean daily runoff, and you’ll want to get the total over every N day event.

1. Note that you need to find the drainage area of the watershed (see Exercise 2).

From exercise 4:

2. Now convert the units from ft/s to mm per event.
3. Now convert precipitation into mm per event (you can put it in a new variable in rain.total.by.event). I’m not showing the code…

F. Plot storm event runoff vs rainfall. I used plot symbol pch=19.

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

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