# Simulating surface runoff of water down a slope.
# Owen Jones, July 2014
#
# We suppose that the system is in a stable state.
# The slope is divided into an m*n grid.
# The rate at which rain falls on each cell is "rain";
# the rate at which water is absorbed by a cell is random,
# with distribution given by function "sleep", with mean 1/rho.
# If the rate at which water enters a cell (from rain or runoff)
# exceeds the infiltration rate, then the excess becomes runoff.
# Runoff flows either to the cell immediately below (probability 1 - 2*delta),
# or to the cell below and left, or below and right (each with probability delta).
#
# x.mat[i,j] is the nett runoff from cell (i,j),
# where i is the position along the slope, and j the position down the slope.
# We calculate x.mat one row at a time, working our way down the slope.
# The runoff is plotted so that heavier runoff is darker
set.seed(100)
rho <- 0.6
seep <- function(k = 1) rexp(k, rho)
rain <- 1
delta <- 0.3
m <- 300
n <- 150
x.mat <- matrix(nrow = n, ncol = m)
x <- pmax(0, rain - seep(m)) # the first row has no runoff from above
x.mat[1,] <- x # save first row
for (i in 2:n) {
# i is the row
# x is the runoff from the row above, y is the runoff for the current row
# first calculate runoff from above
y <- rep(0, m)
# left hand cell of x can only runoff down or down and right
if (runif(1) < delta) { # x[1] runs down and right
y[2] <- y[2] + x[1]
} else { # x[1] runs straight down
y[1] <- y[1] + x[1]
}
# cells x[2] through x[m-1] can run down, down and left, or down and right
for (j in 2:(m-1)) {
u <- runif(1)
if (u < delta) { # x[j] runs down and left
y[j-1] <- y[j-1] + x[j]
} else if (u < 2*delta) { # x[j] runs down and right
y[j+1] <- y[j+1] + x[j]
} else { # x[j] runs straight down
y[j] <- y[j] + x[j]
}
}
# right hand cell of x acn only runoff down or down and left
if (runif(1) < delta) { # x[m] runs down and left
y[m-1] <- y[m-1] + x[m]
} else { # x[m] runs straight down
y[m] <- y[m] + x[m]
}
# next add rainfall and infiltration, and remove excess infiltration
y <- pmax(0, y + rain - seep(m))
x.mat[i,] <- y # save row i
x <- y
}
# output
image(1:m, 1:n, t(x.mat[n:1,]), breaks=c(0:11, max(x.mat)),
col = hsv(h=4/6, alpha=0:11/11), xlab="", ylab="", axes=FALSE)