# 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)