**Dan Kelley Blog/R**, and kindly contributed to R-bloggers)

# Introduction

I got interested in layered sedimentation from viewing a video and decided it would be interesting to code this into R. More on that in due course, but my first step was to code a syatem with one sediment “type”.

# Procedure

The following code drops sediment particles at x=1, and lets them roll downhill until they reach the bottom or a ledge. It draws numbers at the sedimented particles’ final positions. Since the numbers start at 1, the values are like inverse ages.

| ```
m <- 51 # number of particles
n <- 10 # grid width
debug <- FALSE # put TRUE for debugging
info <- function(...) if (debug) cat(...)
pch <- 20
cex <- 4/log2(n)
type <- "t"
set.seed(1)
rollDownhill <- function(X, Z) {
info("rollDownhill(", X, ",", Z, ")\n", sep = "")
if (Z == 1)
return(list(x = X, z = Z))
## Particles roll down-slope until they hit the bottom... ... or a ledge
## comprising two particles.
XX <- X
ZZ <- Z
while (0 == S[XX + 1, ZZ - 1]) {
# move down and to right
info(" XX:", XX, " ZZ:", ZZ, "\n")
XX <- XX + 1
ZZ <- which(0 == S[XX, ])[1]
if (ZZ == 1)
break
if (XX == n)
break
}
return(list(x = XX, z = ZZ))
}
S <- matrix(0, nrow = n, ncol = n) # 'S' means 'space'
par(mar = c(3, 3, 1, 1), mgp = c(2, 0.7, 0))
plot(1:n, 1:n, type = "n", xlab = "", ylab = "")
xDrop <- 1 # location of drop; everything goes here or to right
for (i in 1:m) {
# 'p' means partcle
while (0 == length(zDrop <- which(0 == S[xDrop, ])[1])) {
info("in while line 72\n")
xDrop <- xDrop + 1
if (xDrop == n) {
message("RHS")
break
}
}
info("particle:", i, " ")
p <- rollDownhill(xDrop, zDrop)
S[p$x, p$z] <- 1
if (type == "p") {
points(p$x, p$z, col = "gray", pch = pch, cex = cex)
} else {
text(p$x, p$z, i, col = "gray")
}
}
``` |

# Discussion and conclusions

Reading the numbers on the graph as inverse age, one can see an interesting age structure.

Viewed along diagonals, ages increase by 1 time unit with every lateral step away from the source.

Viewed along Z levels, though, the time step is more interesting. You can see this at a glance, by first-differencing the values along z=1, and then at z=2, etc.

I suppose that if something came along and sliced the sediment mound along z levels, we’d see this more interesting pattern of time variation in the lateral.

I wonder if these patterns (or the code) are of interest to geologists?

# Resources

- Source code: 2014-06-08-slumping-model.R

**leave a comment**for the author, please follow the link and comment on his blog:

**Dan Kelley Blog/R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...