In dynamic-R/RTM: Reaction Transport Modelling: Course Tutorials and Exercises

knitr::opts_chunk$set(echo = TRUE)

library(ReacTran)
library(marelac)

Solutions to exercises

Using R as a calculator

Examples of how to use R as a powerful calculator. Note the use of spaces around operators such as + or -, or around parentheses. They make the code much more legible.

(4/6*8-1)^(2/3)
( 4/6*8 - 1)^(2/3)
log(20)
log2(4096)
2*pi*3
sqrt( 2.3^2 + 5.4^2 - 2*2.3*5.4 * cos (pi/8) )

Vectors and sequences

Mean of a vector

mean(c(9,17))

Sediment depth profiles

depth <- seq(from = 0.05, to = 9.95, by = 0.1)
depth

In the following code, note the use of ==. This operator returns a logical value (TRUE or FALSE) depending on whether or not the value on the left is equal to the value on the right. This logical value is then used within [] to select which elements of the vector are to be considered. (Similar effects are achieved using the logical operators <=, >=, <, >).

porosity <- 0.7 + (1-0.7) * exp(-1*depth)
V <- c(porosity[depth == 0.05], porosity[depth == 9.95])
V
mean(porosity)
mean(porosity[depth <= 1])

When plotting a graph, note that by default the first input variable is plotted on the x axis and the second variable is plotted on the y axis. On the third line we change the orientation of the axis by setting the scale in the "opposite" direction, which results in a more natural display of a depth profile.

par(mfrow=c(1,2))
plot(depth, porosity)
plot(porosity, depth, ylim=c(10,0), xlab="porosity (vol/vol)", ylab="depth (m)")

Estuarine morphology

In the following code, we avoid using constants such as 200 in the line calculating the positions of the middle of the boxes. Instead, we assign this value to a variable (N) and use this variable when necessary and appropriate. Also note that it is not necessary --- and in fact rather confusing, as shown on the "confusing line" in the code below --- to use parentheses to indicate priority of operators such as ^ and *. R adheres to priorities that we are used to from algebra: first ^, then * or /, then + or -. Parentheses should, therefore, only be used if you want to modify these conventional priority rules. Otherwise, it's better to avoid them to improve code legibility.

L  <- 100000 # metres
N  <- 200
dx <- L/N
x  <- seq(from = dx/2, length.out = N, by = dx)
Ar <- 4000
As <- 76000
p  <- 5
ks <- 50000
Area <- Ar + (As-Ar) * ((x^p)/(x^p+ks^p)) # very confusing, too many unnecessary ()
Area <- Ar + (As-Ar) * x^p/(x^p+ks^p)     # much easier to read, equal result
Volume <- Area*dx
sum(Volume)  # m3

Plotting observed data

Oxygen <- c(210, 250, 260, 289, 280, 260, 270, 260)
Hour   <- 1:length(Oxygen)
plot(Hour, Oxygen, type = "l", main = "Oxygen concentration at Jetty", ylab = "mmol/m3")

R-functions

R-function to estimate saturated oxygen concentrations

SatOx <- function(TC = 20, S = 35){ # TC in deg-Celsius
  TK <- TC+273.15                   # TK in Kelvin
  A <- -173.9894 + 25559.07/TK + 146.4813 * log(TK/100) - 22.204*TK/100 + S *
       (-0.037362 + 0.016504*TK/100 - 0.0020564 * TK/100 * TK/100)
  return(exp(A))
}
SatOx()
SatOx(TC = 0:30)

Molecular diffusion coefficient

require(marelac)
DC <- diffcoeff(S=20, t=10, species = c("O2", "CO2"))  # m2/sec
DC*1e4*3600*24                                         # cm2/d

t.seq  <- 1:30                                    # temperature sequence
DC.O2  <- diffcoeff(S=20, t=t.seq)$O2             # m2/sec
DC.CO2 <- diffcoeff(S=20, t=t.seq)[["CO2"]]       # m2/sec (note the alternative to $CO2)
DC.O2.fresh  <- diffcoeff(S=0, t=t.seq)[["O2"]]   # m2/sec

m2_sTOcm2_d <- 1e4*3600*24
yrange <- c(0.5, 2.5)                             # A suitable range for y-axis
plot (t.seq, DC.O2*m2_sTOcm2_d, type = "l", xlab = "Temperature, dgC", 
      ylim = yrange, ylab = "cm2/d", main = "diffusion coefficients")
lines(t.seq, DC.O2.fresh*m2_sTOcm2_d, col = "blue")
points(t.seq, DC.CO2*m2_sTOcm2_d, col = "red")
legend("topleft", lty = 1, col = c("black", "blue", "red"), lwd=c(1,1,NA),
       pch=c(NA,NA,1), legend = c("O2, S=20","O2, S=0","CO2, S=20"), bty="n")

R-function sphere

N  <- 100    # number of boxes
L  <- 100/2  # radius, micrometer
dr <- L/N
r  <- seq(from = dr/2, by = dr, length.out = N)
Sphere <- function(r) { return(4*pi*(r/1000)^2) }  # surface area, mm^2
plot(r, Sphere(r), xlab="radius (mm)", ylab = "area (mm2)")

Porosity profile and estuarine morphology as a function

# porosity as a function of depth
Porfun <- function(depth) {
  return( 0.7 + (1-0.7)*exp(-1*depth) )
}
depth <- seq(from = 0, to = 10, length.out = 100)
plot(Porfun(depth), depth, ylim = c(10,0), xlab = "porosity (vol/vol)", ylab = "depth (cm)")

# estuarine function, returns the cross-section area and volume as a function of distance
Estfun <- function(x, Ar = 4000, As = 76000, p = 5, ks = 50000, dx = 1000){
  Area <- Ar + (As-Ar) * x^p/(x^p+ks^p)
  Volume <- Area*dx
  return(list(Area = Area, Volume = Volume))
}
dx_new <- 2000
x_new  <- seq(from = dx_new/2, length.out = 50, by = dx_new)
# evaluate for different x and dx input values, the others are kept at default
Estfun(x = x_new, dx = dx_new)

Estuarine morphology using ReacTran

Grid <- setup.grid.1D(N = 100, L = 100000)

EstArea <- function(x, Ar = 4000, As = 76000, p = 5, ks = 50000, dx = 1000){
  return( Ar + (As-Ar) * x^p/(x^p+ks^p) )
}

Area <- setup.prop.1D(grid = Grid, func = EstArea)
# example how to pass different input values to the function within setup.prop.1D
Area2 <- setup.prop.1D(grid = Grid, func = EstArea, Ar=10000, As=50000)
par(mfrow=c(1,2))
plot(Area, grid = Grid, ylab = "Area (m2)")
plot(Area2, grid = Grid, ylab = "Area (m2)")

Solving differential equations in R

Lotka-Volterra model

LVmodel <- function(t, state, pars) {
  with (as.list(c(state, pars)), {

    dX.dt <- a*X*(1-X/k) - b*X*Y
    dY.dt <- g*b*X*Y - e*Y

    return (list(c(dX.dt, dY.dt), 
                 sum = X+Y))
  })
}

state.ini <- c(X = 670, Y = 610)
pars      <- c(a = 0.04, k = 1000, b = 5e-5, g = 0.8, e = 0.008)
times <- 1:100
out   <- ode(y = state.ini, func = LVmodel, times = times, parms = pars)
# change initial conditions
state2.ini <- c(X = 100, Y = 540) 
out2       <- ode(y = state2.ini, func = LVmodel, times = times, parms = pars)
# plot results in one graph
plot(out, out2, mfrow=c(1,3), col=c("black","red"), lty=c(1,1))

Now we experiment a little and explore how the predator-prey dynamics look like at longer time scales and for different values of the parameter b.

times <- 1:1500
pars1 <- pars
pars1["b"] <- 1e-5
out1  <- ode(y = state.ini, func = LVmodel, times = times, parms = pars1)
pars2 <- pars
pars2["b"] <- 1.3e-5
out2  <- ode(y = state.ini, func = LVmodel, times = times, parms = pars2)
pars3 <- pars
pars3["b"] <- 5e-5
out3  <- ode(y = state.ini, func = LVmodel, times = times, parms = pars3)
pars4 <- pars
pars4["b"] <- 10e-5
out4  <- ode(y = state.ini, func = LVmodel, times = times, parms = pars4)
plot(out1,out2,out3,out4, mfrow=c(1,3), col=c(1,2,3,4), lty=1)
legend("right",title="Parameter b", col=1:4, lty=1,
       legend=c(pars1["b"],pars2["b"],pars3["b"],pars4["b"]) )

Lorenz model

Lorenz <- function(t, state, parms) {
  with (as.list(c(state, parms)), {
    dX.dt <- -8/3*X + Y*Z
    dY.dt <- -10*(Y-Z)
    dZ.dt <- -X*Y + 28*Y - Z
    return (list(c(dX.dt, dY.dt, dZ.dt), 
                 sum = X+Y+Z))
  })
}
pars      <-  NULL
state.ini <- c(X = 1, Y = 1, Z = 1)
time.seq  <- seq(from = 0, to = 100, by = 0.005)
# solve it
out <- ode(y = state.ini, times = time.seq, func = Lorenz, parms=pars )
head(out, n=4)
plot(out, xlab = "time", lwd = 2)

dynamic-R/RTM documentation built on Feb. 28, 2025, 1:23 p.m.