slumping model

June 8, 2014
By

(This article was first published on 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.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
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")
    }
}

center

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

To 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...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.