Level pool routing refers to one of the more straightforward methods for calculating reservoir outflow given an input hydrograph (time vs inflow) along with information about basin discharge relative to elevation. Here the term reservoir is used in the technical engineering context and does not exclude the use of this method for natural lakes. I believe this method is applicable to any waterbody where you can assume a onetoone relationship between discharge and elevation. Said another way, we assume that discharge does not have a hysteresis component and does not depend on the direction (rising vs falling limb) of the input hydrograph.
For more a more rigorous mathematical treatment of this method see here, here, and here. The following demonstration will lay out the computational details and describe code to implement the method on three test datasets from the above links.
Computation
The level pool routing computation procedure involves:

Determining the volume of the reservoir at a range of depths (if not known aheadoftime) given its area.

Developing a relationship between outflow and reservoirchangeinstorage.

The reference examples use the classical approach of linear interpolation. Here I use nonlinear generalized additive modelling to more flexibly represent the shape of this relationship.
Then for each timestep we:

Calculate the sum of inflows for the current and previous timesteps.

Calculate changeinstoragewithtime as a function of outflow during the previous timestep.

Use our fitted relationship to calculate outflows for the current timestep as a function of reservoir changeinstorage.

Subtract these outflows from the running reservoir storage term.
#' level_pool_routing
#' @param lt data.frame with time and inflow columns
#' @param qh data.frame with elevation and discharge columns.
#' Storage column optional.
#' @param area numeric reservoir area
#' @param delta_t numeric time step interval in seconds.
#' @param initial_outflow numeric
#' @param initial_storage numeric
#' @param linear.fit logical operator specifying a linear
#' relationship between outflow and reservoirchangeinstorage
#' @importFrom mgcv gam
level_pool_routing < function(lt, qh, area, delta_t,
initial_outflow, initial_storage,
linear.fit){
lagpad < function(x, k) {
c(rep(NA, k), x)[1 : length(x)]
}
lt$ii < apply(cbind(lagpad(lt$inflow, 1), lt$inflow), 1, sum)
if(is.null(qh$storage)){
qh$storage < area * qh$elevation
}
qh$stq < ((2 * qh$storage) / (delta_t)) + qh$discharge
lt$sjtminq < NA
lt$sj1tplusq < NA
lt$outflow < NA
lt[1, c("sj1tplusq")] < c(NA)
lt[1, c("sjtminq")] < ((2 * initial_storage) / delta_t) 
initial_outflow
lt[1, "outflow"] < initial_outflow
if(linear.fit == TRUE){
fit < lm(discharge ~ stq, data = qh)
}else{
fit < mgcv::gam(discharge ~ s(stq, k = 3), data = qh)
plot(qh$stq, qh$discharge, xlab = "Changeinstoragewithtime",
ylab = "Discharge")
lines(qh$stq, predict(fit))
}
for(i in seq_len(nrow(lt))[1]){
lt[i, "sj1tplusq"] < lt[i1, "sjtminq"] + lt[i, "ii"]
lt[i, "outflow"] < predict(fit,
data.frame(stq = lt[i, "sj1tplusq"]))
lt[i, "sjtminq"] < lt[i, "sj1tplusq"] 
(lt[i, "outflow"] * 2)
}
lt
}
Example 1
Linear relationship between discharge and changeinstorage
The data for this example comes from this levelpool routing walkthrough. We are given an inflow hydrograph in 6 hour increments so we will specify a delta_t
timestep of 6 * 3600
seconds. The problem statement specifies an intial storage of 0 and an initial outflow of 20. There is no need to specify reservoir area because we are given storage as a function of discharge.
input_hydro < data.frame(
time = seq(0, 162, by = 6),
inflow = c(0, 50, 130, 250, 350, 540, 735, 1215,
1800, 1400, 1050, 900, 740, 620, 510,
420, 320, 270, 200, 150, 100, 72, 45,
25, 10, 0, 0, 0))
reservoir_char < data.frame(
elevation = c(130:134, 136:139),
discharge = c(20, 34, 57, 96, 162, 463, 781,
1318, 2226),
storage = c(1, 1.69, 2.85, 4.8, 8.12, 23.1,
39.1, 65.9, 111))
reservoir_char$storage < reservoir_char$storage * 1000000
delta_t < 6 * 3600
res_linear < level_pool_routing(input_hydro, reservoir_char,
area = NA, delta_t = delta_t, initial_outflow = 20,
initial_storage = 0, linear.fit = FALSE)
plot(res_linear$time, res_linear$inflow,
xlab = "Time (h)", ylab = "Flow (m3/s)")
lines(res_linear$time, res_linear$outflow)
legend("topleft", legend = c("Inflow", "Outflow"), lty = c(NA, 1),
pch = c(21, NA))
Example 2
Curvilinear relationship between discharge and changeinstorage
The data for this example comes from this levelpool routing walkthrough pdf. I scraped the data tables from the pdf using the pdftools package. We are given storage as a function of discharge so we have no need for information on the area of the reservoir. Also, we are given an inflow hydrograph in 2 hour increments so we will specify a delta_t
timestep of 2 * 3600
seconds. Unlike the previous example, our reservoir has a nonzero initial storage.
library(pdftools)
txt < strsplit(pdf_text(
"https://www.caee.utexas.edu/prof/maidment/CE374KSpr12/Docs/Hmwk5Soln.pdf"), "\n")[[1]]
parse_table < function(tbl, tbl_names){
tbl < lapply(tbl, function(x) read.table(text = x,
stringsAsFactors = FALSE))
tbl < lapply(tbl, function(x) gsub(",", "", x))
inds < lapply(tbl,
function(x) ifelse(max(grep(")", x)) == Inf,
1, max(grep(")", x))))
inds < lapply(seq_along(tbl),
function(x) c(1, (inds[[x]] + 1):length(tbl[[x]])))
tbl < lapply(seq_along(tbl),
function(x) as.numeric(tbl[[x]][inds[[x]]]))
tbl < data.frame(t(do.call("rbind", tbl)))
tbl < tbl[2:nrow(tbl),]
names(tbl) < tbl_names
tbl
}
tbl1 < suppressWarnings(rbind(
parse_table(txt[6:7], tbl_names = c("time", "inflow")),
data.frame(time = seq(20, 36, by = 2), inflow = 0)))
tbl2 < suppressWarnings(
parse_table(txt[9:10], tbl_names = c("storage", "discharge")))
tbl2$storage < tbl2$storage * 1000000
res_curv < level_pool_routing(lt = tbl1, qh = tbl2, area = NA,
delta_t = 7200, initial_outflow = 57,
initial_storage = 75000000, linear.fit = FALSE)
plot(res_curv$time, res_curv$inflow, xlab = "Time (h)",
ylab = "Flow (m3/s)")
lines(res_curv$time, res_curv$outflow)
legend("topleft", legend = c("Inflow", "Outflow"), lty = c(NA, 1),
pch = c(21, NA))
Example 3
Oscillating relationship between discharge and changeinstorage
The data for this example comes from this levelpool routing walkthrough ppt. Here we are given discharge as a function of reservoir elevation but we are not given the corresponding storage. As a result, we must specify a reservoir area. We are given an inflow hydrograph in 10 minute increments so we will specify a delta_t
timestep of 10 * 60
seconds. In this case, our reservoir a zero initial storage and initial outflow.
lt < data.frame(time = seq(0, 210, by = 10),
inflow = c(seq(0, 360, by = 60), seq(320, 0, by = 40),
rep(0, 6)))
qh < data.frame(elevation = seq(0, 10, by = 0.5),
discharge = c(0, 3, 8, 17, 30, 43, 60, 78, 97, 117, 137,
156, 173, 190, 205, 218, 231, 242, 253, 264, 275))
res_osc < level_pool_routing(lt, qh, area = 43560, delta_t = 600,
initial_outflow = 0, initial_storage = 0,
linear.fit = FALSE)
plot(res_osc$time, res_osc$inflow, xlab = "Time (h)",
ylab = "Flow (cfs)")
lines(res_osc$time, res_osc$outflow)
legend("topleft", legend = c("Inflow", "Outflow"), lty = c(NA, 1),
pch = c(21, NA))