EVSC 3600 Homework 4
Charles Hanley
03/15/2015
Problem 1
# Input data
year <- seq(1940, 1949)
discharge <- c(4060,1111,2100,3780,5320,1887,1223,1089,2338,3752)
ln_discharge <- log(discharge)
discharge_rank <- 11 - rank(discharge) #Subtract from 11 to invert ranks
# Calculate exceedence probabilities
exceedance_probability <- rep(0, length(discharge))
for (i in 1:length(exceedance_probability)){
exceedance_probability[i] <- sum(discharge > discharge[i]) / length(discharge)
}
# View requested chart
cbind(year, discharge, ln_discharge, discharge_rank, exceedance_probability)## year discharge ln_discharge discharge_rank exceedance_probability
## [1,] 1940 4060 8.308938 2 0.1
## [2,] 1941 1111 7.013016 9 0.8
## [3,] 1942 2100 7.649693 6 0.5
## [4,] 1943 3780 8.237479 3 0.2
## [5,] 1944 5320 8.579229 1 0.0
## [6,] 1945 1887 7.542744 7 0.6
## [7,] 1946 1223 7.109062 8 0.7
## [8,] 1947 1089 6.993015 10 0.9
## [9,] 1948 2338 7.757051 5 0.4
## [10,] 1949 3752 8.230044 4 0.3# Return interval in years for peak flow of 1948:
1/exceedance_probability[which(year==1948)]## [1] 2.5Problem 2
# Input data
vol.init <- 4*10^6 #m^3
outflow <- 12.5 #m^3/s
inflow <- c(12,12,12.6,15.9,22.3,35.8,28.2,24.1,21.3,18.2,17,16.2,15.5,14.9,
14.5,14.1,13.7,13.5,13.3,13.1,13,12.9,12.85,12.8,12.82,12.81,12.8,
12.75,12.7,12.72,12.71,12.69,12.68) #m^3/s
netflow <- inflow - outflow
time <- seq(from=0, to=8, by=0.25) #days
# One timestep equals how many seconds?
timestep <- 60*60*6 #seconds/step
netflow.per_step <- netflow * timestep #m^3/step
# Calculate volume at each timestep iteratively.
volume <- rep(0, 33)
for (i in 1:length(volume)){
if (i==1){volume[i] <- vol.init
} else volume[i] <- volume[i-1] + netflow.per_step[i-1]
}
# View table (time in days, flows in m^3/s, volume in m^3)
cbind(time, inflow, outflow, netflow, volume)## time inflow outflow netflow volume
## [1,] 0.00 12.00 12.5 -0.50 4000000
## [2,] 0.25 12.00 12.5 -0.50 3989200
## [3,] 0.50 12.60 12.5 0.10 3978400
## [4,] 0.75 15.90 12.5 3.40 3980560
## [5,] 1.00 22.30 12.5 9.80 4054000
## [6,] 1.25 35.80 12.5 23.30 4265680
## [7,] 1.50 28.20 12.5 15.70 4768960
## [8,] 1.75 24.10 12.5 11.60 5108080
## [9,] 2.00 21.30 12.5 8.80 5358640
## [10,] 2.25 18.20 12.5 5.70 5548720
## [11,] 2.50 17.00 12.5 4.50 5671840
## [12,] 2.75 16.20 12.5 3.70 5769040
## [13,] 3.00 15.50 12.5 3.00 5848960
## [14,] 3.25 14.90 12.5 2.40 5913760
## [15,] 3.50 14.50 12.5 2.00 5965600
## [16,] 3.75 14.10 12.5 1.60 6008800
## [17,] 4.00 13.70 12.5 1.20 6043360
## [18,] 4.25 13.50 12.5 1.00 6069280
## [19,] 4.50 13.30 12.5 0.80 6090880
## [20,] 4.75 13.10 12.5 0.60 6108160
## [21,] 5.00 13.00 12.5 0.50 6121120
## [22,] 5.25 12.90 12.5 0.40 6131920
## [23,] 5.50 12.85 12.5 0.35 6140560
## [24,] 5.75 12.80 12.5 0.30 6148120
## [25,] 6.00 12.82 12.5 0.32 6154600
## [26,] 6.25 12.81 12.5 0.31 6161512
## [27,] 6.50 12.80 12.5 0.30 6168208
## [28,] 6.75 12.75 12.5 0.25 6174688
## [29,] 7.00 12.70 12.5 0.20 6180088
## [30,] 7.25 12.72 12.5 0.22 6184408
## [31,] 7.50 12.71 12.5 0.21 6189160
## [32,] 7.75 12.69 12.5 0.19 6193696
## [33,] 8.00 12.68 12.5 0.18 6197800# Plot reservoir volume
library(ggplot2)
ggplot(data=data.frame(time=time, volume=volume/1000000), aes(x=time, y=volume)) +
geom_line(color="blue") +
geom_point(color="blue") +
scale_x_continuous("Time (days)") +
scale_y_continuous("Volume (million cubic meters)", limits=c(0,6.5)) +
ggtitle("Reservoir Volume through Flood") +
theme_bw()
Problem 3
# Input data
U <- 2 #m/s
L <- 25000 #m
x <- 0.2
delta_t <- 6*60*60 #s
outflow.init <- 12 #m^3/s
muskingum <- function(inflow, U, L, x, delta_t, outflow.init){
# Muskingum method of flood routing. Requires inflow hydrograph, returns
# outflow hydrograph. Assumes kt = L/U.
kt <- L/U
c0 <- (-kt*x + 0.5*delta_t)/(kt - kt*x + 0.5*delta_t)
c1 <- (kt*x + 0.5*delta_t)/(kt - kt*x + 0.5*delta_t)
c2 <- (kt - kt*x - 0.5*delta_t)/(kt - kt*x + 0.5*delta_t)
outflow <- rep(0, length(inflow))
outflow[1] <- outflow.init
for (i in 2:length(inflow)){
outflow[i] <- c0*inflow[i] + c1*inflow[i-1] + c2*outflow[i-1]
}
return(outflow)
}
# A:
hydrograph.25m <- muskingum(inflow, U, L, x, delta_t, outflow.init)
ggplot(data=data.frame(flow=hydrograph.25m, time=time), aes(x=time, y=flow)) +
geom_line() +
geom_point() +
scale_x_continuous("Time (days)") +
scale_y_continuous("Discharge (cubic meters per second)") +
ggtitle("Flood Hydrograph 25km Downstream") +
theme_bw()
# B:
# The peak outflow (in m^3/s) 25m downstream is:
max(hydrograph.25m)## [1] 33.07377# C:
hydrograph.50m <- muskingum(hydrograph.25m, U, L, x, delta_t, outflow.init)
ggplot(data=data.frame(flow=hydrograph.50m, time=time), aes(x=time, y=flow)) +
geom_line() +
geom_point() +
scale_x_continuous("Time (days)") +
scale_y_continuous("Discharge (cubic meters per second)") +
ggtitle("Flood Hydrograph 50km Downstream") +
theme_bw()
# D:
# The peak outflow (in m^3/s) 50m downstream is:
max(hydrograph.50m)## [1] 30.51434# Plot all three hydrographs:
three_hydrographs <- as.data.frame(cbind(c(inflow, hydrograph.25m, hydrograph.50m),
rep(time, 3), c(rep("0 m", 33),
rep("25 m", 33),
rep("50 m", 33))))
names(three_hydrographs) <- c("flow", "time", "distance")
three_hydrographs$flow <- as.numeric(as.character(three_hydrographs$flow))
three_hydrographs$time <- as.numeric(as.character(three_hydrographs$time))
three_hydrographs$distance <- as.character(three_hydrographs$distance)
ggplot(data=three_hydrographs, aes(x=time, y=flow, color=distance)) +
geom_line() +
geom_point() +
theme_bw() +
scale_y_continuous("Discharge (cubic meters per second)") +
scale_x_continuous("Time (days)") +
ggtitle("Flood Hydrographs at Three Reaches")
This document was made with RMarkdown.