In dynamic-R/RTM: Reaction Transport Modelling: Course Tutorials and Exercises
knitr::opts_chunk$set(echo = TRUE)
Choice of data and model parameters
dataset <- 3
d <- c(dataset==1, dataset==2, dataset==3)
All data sets and default model parameters
Data set 1
# experimental data, set1
t.exp <- c(0, 22, 94, 166, 262, 358, 430, 526)
CH4.exp <- c(22.18, 21.74, 20.08, 18.46, 14.08, 11.31, 10.20, 8.91)
CO2.exp <- c(2.95, 4.51, 10.10, 12.10, 13.03, 14.04, 17.35, 17.95)
xCO2.exp <- c(NA, 12.90, 17.65, 21.60, 24.83, 26.26, 27.21, 28.10)*1e-2
DATA1 <- data.frame(time = t.exp, CH4 = CH4.exp, CO2 = CO2.exp, xCO2 = xCO2.exp)
# default model parameters
pars1 <- c(
CH4.ini = 22.18, # all values in umol/bottle
CO2.ini = 2.95,
Cdead.ini = 1.15*8.69-2.95,
Clive.ini = 1,
rMO = 0.0274, # CH4 removal rate (umol/bottle/h),
kMin = 3*3/530, # mineralization of dead biomass (/h)
eff = 0.45, # assimilation efficiency
xCH4.src = 0.5, #0.41,
xCdead.src = 0.2, #0.23,
xCO2.ini = 0.01,
xClive.ini = 0.01
)
Data set 2
# experimental data, set2
t.exp <- c(0, 22, 94, 166, 262, 358, 430, 526)
CH4.exp <- c(23.80, 23.00, 20.28, 19.10, 14.89, 11.98, 10.73, 9.38)
CO2.exp <- c(2.00, 4.87, 9.67, 12.20, 14.10, 16.05, 18.87, 21.71)
xCO2.exp <- c(NA, 15.95, 16.37, 19.70, 23.30, 23.86, 23.56, 23.22)*1e-2
DATA2 <- data.frame(time = t.exp, CH4 = CH4.exp, CO2 = CO2.exp, xCO2 = xCO2.exp)
# default model parameters
pars2 <- c(
CH4.ini = 23.80, # all values in umol/bottle
CO2.ini = 2.0,
Cdead.ini = 6.2,
Clive.ini = 1,
rMO = 0.031, # CH4 removal rate (umol/bottle/h),
kMin = 3*3/530, # mineralization of dead biomass (/h)
eff = 0.2, # assimilation efficiency
xCH4.src = 0.26, #0.41,
xCdead.src = 0.27, #0.23,
xCO2.ini = 0.01,
xClive.ini = 0.01
)
Data set 3
# experimental data, set3
t.exp <- c(0, 22, 94, 166, 262, 358, 430, 526)
CH4.exp <- c(23.66, 23.48, 21.22, 19.77, 16.41, 14.58, 13.43, 10.87)
CO2.exp <- c(2.96, 3.74, 11.86, 15.38, 17.68, 20.69, 23.54, 24.17)
xCO2.exp <- c(NA, 9.46, 10.40, 13.81, 16.32, 17.62, 18.24, 18.91)*1e-2
DATA3 <- data.frame(time = t.exp, CH4 = CH4.exp, CO2 = CO2.exp, xCO2 = xCO2.exp)
# default model parameters
pars3 <- c(
CH4.ini = 23.65, # all values in umol/bottle
CO2.ini = 2.96,
Cdead.ini = 11.63-2.96,
Clive.ini = 1, # arbitrary
rMO = 0.0245, # CH4 removal rate (umol/bottle/h),
kMin = 3*3/530, # mineralization of dead biomass (/h)
eff = 0, # assimilation efficiency
xCH4.src = 0.278,
xCdead.src = 0.13,
xCO2.ini = 0.01,
xClive.ini = 0.01
)
Model function
C13model <-function(t, state, parms) {
with (as.list(c(state, parms)),{
xCH4 <- CH4_C13/CH4
xCO2 <- CO2_C13/CO2
xCdead <- Cdead_C13/Cdead
xClive <- Clive_C13/Clive
rMin <- kMin*Cdead # mineralization of dead Corg
# total C balances
dCH4.dt <- -rMO
dCO2.dt <- (1-eff)*rMO + rMin
dCdead.dt <- -rMin
dClive.dt <- eff*rMO
# 13C balances
dCH4_C13.dt <- -rMO * xCH4.src
dCO2_C13.dt <- (1-eff)*rMO * xCH4.src + rMin * xCdead.src
dCdead_C13.dt <- -rMin * xCdead.src
dClive_C13.dt <- eff*rMO * xCH4.src
list(
c(dCH4.dt, dCH4_C13.dt, dCO2.dt, dCO2_C13.dt,
dCdead.dt, dCdead_C13.dt, dClive.dt, dClive_C13.dt),
xCH4 = xCH4, xCO2 = xCO2, xCdead = xCdead, xClive = xClive
)
})
}
# function for calculating initial state from model parameters
state.ini <- function(parms){
with(as.list(parms), {
c(CH4 = CH4.ini, CH4_C13 = CH4.ini*xCH4.src,
CO2 = CO2.ini, CO2_C13 = CO2.ini*xCO2.ini,
Cdead = Cdead.ini, Cdead_C13 = Cdead.ini*xCdead.src,
Clive = Clive.ini, Clive_C13 = Clive.ini*xClive.ini
)})
}
Selected data set: L+r paste(which(d))
if (d[1]){ DATA <- DATA1; pars <- pars1 }
if (d[2]){ DATA <- DATA2; pars <- pars2 }
if (d[3]){ DATA <- DATA3; pars <- pars3 }
Model parameters for sensitivity analysis
p0 <- p1 <- p2 <- p3 <- pars
p1["Cdead.ini"] <- 0
p1["xCH4.src"] <- 0.4
p2["eff"] <- 0.0
p2["xCH4.src"] <- 0.38
p3["xCH4.src"] <- 0.4
p3["xCdead.src"] <- 0.26
p0 <- p1 <- p2 <- p3 <- pars
p1["Cdead.ini"] <- 0
p1["xCH4.src"] <- 0.28
p2["eff"] <- 0.0
p3["xCH4.src"] <- 0.5
p0 <- p1 <- p2 <- p3 <- pars
p1["Cdead.ini"] <- 0
p1["xCH4.src"] <- 0.23
p2["eff"] <- 0.7
p2["xCH4.src"] <- 0.5
p3["xCH4.src"] <- 0.5
Solve the model
require(deSolve)
outtimes <- seq(from = 0, to = 530, length.out = 100)
out0 <- ode(y = state.ini(p0), parms = p0, func = C13model, times = outtimes)
out1 <- ode(y = state.ini(p1), parms = p1, func = C13model, times = outtimes)
out2 <- ode(y = state.ini(p2), parms = p2, func = C13model, times = outtimes)
out3 <- ode(y = state.ini(p3), parms = p3, func = C13model, times = outtimes)
Full model results
# display complete output (for debugging)
plot(out0, out1, out2, out3,
which = c("CH4", "CO2", "Cdead", "Clive",
"CH4_C13", "CO2_C13", "Cdead_C13", "Clive_C13",
"xCH4", "xCO2", "xCdead", "xClive"),
xlab="time (h)", lty=1, lwd=c(1,2,1,2), mfrow=c(3,4))
Model results with experimental data
# display exp. data and model predictions
plot(out0, out1, out2, out3,
obs = DATA,
obspar = list(pch=16, col=which(d)),
mfrow=c(1,3), lty=c(1,2,1,2), lwd=c(1,2,1,2))
legend("bottomright", legend=0:3, lty=c(1,2,1,2), lwd=c(1,2,1,2), col=1:4)
All experimental data and model results
out1 <- ode(y = state.ini(pars1), parms = pars1, func = C13model, times = outtimes)
out2 <- ode(y = state.ini(pars2), parms = pars2, func = C13model, times = outtimes)
out3 <- ode(y = state.ini(pars3), parms = pars3, func = C13model, times = outtimes)
pars1_50 <- pars1; pars1_50["xCH4.src"] <- 0.5
pars2_50 <- pars2; pars2_50["xCH4.src"] <- 0.5
pars3_50 <- pars3; pars3_50["xCH4.src"] <- 0.5
out1_50 <- ode(y = state.ini(pars1_50), parms = pars1_50, func = C13model, times = outtimes)
out2_50 <- ode(y = state.ini(pars2_50), parms = pars2_50, func = C13model, times = outtimes)
out3_50 <- ode(y = state.ini(pars3_50), parms = pars3_50, func = C13model, times = outtimes)
# display exp. data and model predictions
plot(out1, out2, out3, out1_50, out2_50, out3_50,
obs = list(DATA1, DATA2, DATA3),
obspar = list(pch=16, col=1:3),
mfrow=c(1,3), lty=c(1,1,1,2,2,2), col=c(1,2,3,1,2,3), lwd=1)
legend("bottomright", legend=1:3, lty=1, lwd=2, col=1:3)
Reference
- Rob J. M. van Spanning, Qingtian Guan, Chrats Melkonian, James Gallant, Lubos Polerecky, Jean-François Flot, Bernd W. Brandt, Martin Braster, Paul Iturbe Espinoza, Joost W. Aerts, Marion M. Meima-Franke, Sander R. Piersma, Catalin M. Bunduc, Roy Ummels, Arnab Pain, Emily J. Fleming, Nicole N. van der Wel, Vasile D. Gherman, Serban M. Sarbu, Paul L. E. Bodelier, and Wilbert Bitter (2022) Methanotrophy by a Mycobacterium species that dominates a cave microbial ecosystem. Nature Microbiology 7, 2089–2100. https://doi.org/10.1038/s41564-022-01252-3
dynamic-R/RTM documentation built on Feb. 28, 2025, 1:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
Choice of data and model parameters
dataset <- 3 d <- c(dataset==1, dataset==2, dataset==3)
All data sets and default model parameters
Data set 1
# experimental data, set1 t.exp <- c(0, 22, 94, 166, 262, 358, 430, 526) CH4.exp <- c(22.18, 21.74, 20.08, 18.46, 14.08, 11.31, 10.20, 8.91) CO2.exp <- c(2.95, 4.51, 10.10, 12.10, 13.03, 14.04, 17.35, 17.95) xCO2.exp <- c(NA, 12.90, 17.65, 21.60, 24.83, 26.26, 27.21, 28.10)*1e-2 DATA1 <- data.frame(time = t.exp, CH4 = CH4.exp, CO2 = CO2.exp, xCO2 = xCO2.exp) # default model parameters pars1 <- c( CH4.ini = 22.18, # all values in umol/bottle CO2.ini = 2.95, Cdead.ini = 1.15*8.69-2.95, Clive.ini = 1, rMO = 0.0274, # CH4 removal rate (umol/bottle/h), kMin = 3*3/530, # mineralization of dead biomass (/h) eff = 0.45, # assimilation efficiency xCH4.src = 0.5, #0.41, xCdead.src = 0.2, #0.23, xCO2.ini = 0.01, xClive.ini = 0.01 )
Data set 2
# experimental data, set2 t.exp <- c(0, 22, 94, 166, 262, 358, 430, 526) CH4.exp <- c(23.80, 23.00, 20.28, 19.10, 14.89, 11.98, 10.73, 9.38) CO2.exp <- c(2.00, 4.87, 9.67, 12.20, 14.10, 16.05, 18.87, 21.71) xCO2.exp <- c(NA, 15.95, 16.37, 19.70, 23.30, 23.86, 23.56, 23.22)*1e-2 DATA2 <- data.frame(time = t.exp, CH4 = CH4.exp, CO2 = CO2.exp, xCO2 = xCO2.exp) # default model parameters pars2 <- c( CH4.ini = 23.80, # all values in umol/bottle CO2.ini = 2.0, Cdead.ini = 6.2, Clive.ini = 1, rMO = 0.031, # CH4 removal rate (umol/bottle/h), kMin = 3*3/530, # mineralization of dead biomass (/h) eff = 0.2, # assimilation efficiency xCH4.src = 0.26, #0.41, xCdead.src = 0.27, #0.23, xCO2.ini = 0.01, xClive.ini = 0.01 )
Data set 3
# experimental data, set3 t.exp <- c(0, 22, 94, 166, 262, 358, 430, 526) CH4.exp <- c(23.66, 23.48, 21.22, 19.77, 16.41, 14.58, 13.43, 10.87) CO2.exp <- c(2.96, 3.74, 11.86, 15.38, 17.68, 20.69, 23.54, 24.17) xCO2.exp <- c(NA, 9.46, 10.40, 13.81, 16.32, 17.62, 18.24, 18.91)*1e-2 DATA3 <- data.frame(time = t.exp, CH4 = CH4.exp, CO2 = CO2.exp, xCO2 = xCO2.exp) # default model parameters pars3 <- c( CH4.ini = 23.65, # all values in umol/bottle CO2.ini = 2.96, Cdead.ini = 11.63-2.96, Clive.ini = 1, # arbitrary rMO = 0.0245, # CH4 removal rate (umol/bottle/h), kMin = 3*3/530, # mineralization of dead biomass (/h) eff = 0, # assimilation efficiency xCH4.src = 0.278, xCdead.src = 0.13, xCO2.ini = 0.01, xClive.ini = 0.01 )
Model function
C13model <-function(t, state, parms) { with (as.list(c(state, parms)),{ xCH4 <- CH4_C13/CH4 xCO2 <- CO2_C13/CO2 xCdead <- Cdead_C13/Cdead xClive <- Clive_C13/Clive rMin <- kMin*Cdead # mineralization of dead Corg # total C balances dCH4.dt <- -rMO dCO2.dt <- (1-eff)*rMO + rMin dCdead.dt <- -rMin dClive.dt <- eff*rMO # 13C balances dCH4_C13.dt <- -rMO * xCH4.src dCO2_C13.dt <- (1-eff)*rMO * xCH4.src + rMin * xCdead.src dCdead_C13.dt <- -rMin * xCdead.src dClive_C13.dt <- eff*rMO * xCH4.src list( c(dCH4.dt, dCH4_C13.dt, dCO2.dt, dCO2_C13.dt, dCdead.dt, dCdead_C13.dt, dClive.dt, dClive_C13.dt), xCH4 = xCH4, xCO2 = xCO2, xCdead = xCdead, xClive = xClive ) }) } # function for calculating initial state from model parameters state.ini <- function(parms){ with(as.list(parms), { c(CH4 = CH4.ini, CH4_C13 = CH4.ini*xCH4.src, CO2 = CO2.ini, CO2_C13 = CO2.ini*xCO2.ini, Cdead = Cdead.ini, Cdead_C13 = Cdead.ini*xCdead.src, Clive = Clive.ini, Clive_C13 = Clive.ini*xClive.ini )}) }
Selected data set: L+r paste(which(d))
if (d[1]){ DATA <- DATA1; pars <- pars1 } if (d[2]){ DATA <- DATA2; pars <- pars2 } if (d[3]){ DATA <- DATA3; pars <- pars3 }
Model parameters for sensitivity analysis
p0 <- p1 <- p2 <- p3 <- pars p1["Cdead.ini"] <- 0 p1["xCH4.src"] <- 0.4 p2["eff"] <- 0.0 p2["xCH4.src"] <- 0.38 p3["xCH4.src"] <- 0.4 p3["xCdead.src"] <- 0.26
p0 <- p1 <- p2 <- p3 <- pars p1["Cdead.ini"] <- 0 p1["xCH4.src"] <- 0.28 p2["eff"] <- 0.0 p3["xCH4.src"] <- 0.5
p0 <- p1 <- p2 <- p3 <- pars p1["Cdead.ini"] <- 0 p1["xCH4.src"] <- 0.23 p2["eff"] <- 0.7 p2["xCH4.src"] <- 0.5 p3["xCH4.src"] <- 0.5
Solve the model
require(deSolve) outtimes <- seq(from = 0, to = 530, length.out = 100) out0 <- ode(y = state.ini(p0), parms = p0, func = C13model, times = outtimes) out1 <- ode(y = state.ini(p1), parms = p1, func = C13model, times = outtimes) out2 <- ode(y = state.ini(p2), parms = p2, func = C13model, times = outtimes) out3 <- ode(y = state.ini(p3), parms = p3, func = C13model, times = outtimes)
Full model results
# display complete output (for debugging) plot(out0, out1, out2, out3, which = c("CH4", "CO2", "Cdead", "Clive", "CH4_C13", "CO2_C13", "Cdead_C13", "Clive_C13", "xCH4", "xCO2", "xCdead", "xClive"), xlab="time (h)", lty=1, lwd=c(1,2,1,2), mfrow=c(3,4))
Model results with experimental data
# display exp. data and model predictions plot(out0, out1, out2, out3, obs = DATA, obspar = list(pch=16, col=which(d)), mfrow=c(1,3), lty=c(1,2,1,2), lwd=c(1,2,1,2)) legend("bottomright", legend=0:3, lty=c(1,2,1,2), lwd=c(1,2,1,2), col=1:4)
All experimental data and model results
out1 <- ode(y = state.ini(pars1), parms = pars1, func = C13model, times = outtimes) out2 <- ode(y = state.ini(pars2), parms = pars2, func = C13model, times = outtimes) out3 <- ode(y = state.ini(pars3), parms = pars3, func = C13model, times = outtimes) pars1_50 <- pars1; pars1_50["xCH4.src"] <- 0.5 pars2_50 <- pars2; pars2_50["xCH4.src"] <- 0.5 pars3_50 <- pars3; pars3_50["xCH4.src"] <- 0.5 out1_50 <- ode(y = state.ini(pars1_50), parms = pars1_50, func = C13model, times = outtimes) out2_50 <- ode(y = state.ini(pars2_50), parms = pars2_50, func = C13model, times = outtimes) out3_50 <- ode(y = state.ini(pars3_50), parms = pars3_50, func = C13model, times = outtimes)
# display exp. data and model predictions plot(out1, out2, out3, out1_50, out2_50, out3_50, obs = list(DATA1, DATA2, DATA3), obspar = list(pch=16, col=1:3), mfrow=c(1,3), lty=c(1,1,1,2,2,2), col=c(1,2,3,1,2,3), lwd=1) legend("bottomright", legend=1:3, lty=1, lwd=2, col=1:3)
Reference
- Rob J. M. van Spanning, Qingtian Guan, Chrats Melkonian, James Gallant, Lubos Polerecky, Jean-François Flot, Bernd W. Brandt, Martin Braster, Paul Iturbe Espinoza, Joost W. Aerts, Marion M. Meima-Franke, Sander R. Piersma, Catalin M. Bunduc, Roy Ummels, Arnab Pain, Emily J. Fleming, Nicole N. van der Wel, Vasile D. Gherman, Serban M. Sarbu, Paul L. E. Bodelier, and Wilbert Bitter (2022) Methanotrophy by a Mycobacterium species that dominates a cave microbial ecosystem. Nature Microbiology 7, 2089–2100. https://doi.org/10.1038/s41564-022-01252-3
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.