In dynamic-R/RTM: Reaction Transport Modelling: Course Tutorials and Exercises
knitr::opts_chunk$set(echo = TRUE)
Answers
Rate expressions
$$\begin{array}{l}
AeroMin = r_{Min} \cdot \frac{O_2}{O_2+k_{O_2}} \cdot Corg \[5mm]
Nitri = r_{Nit} \cdot \frac{O_2}{O_2+k_{O_2}} \cdot NH_3 \[5mm]
Denitri = r_{Min} \cdot \frac{k_{O_2}}{O_2+k_{O_2}} \cdot \frac{NO_3^-}{NO_3^-+k_{NO_3}} \cdot Corg \[5mm]
Sred = r_{Min} \cdot \frac{k_{O_2}}{O_2+k_{O_2}} \cdot \frac{k_{NO_3}}{NO_3^-+k_{NO_3}} \cdot \frac{SO_4^{2-}}{SO_4^{2-}+k_{SO_4}} \cdot Corg \[5mm]
Sox = r_{Sox} \cdot \frac{O_2}{O_2+k_{O_2}} \cdot H_2S \[5mm]
Methanogen = r_{Min} \cdot \frac{k_{O_2}}{O_2+k_{O_2}} \cdot \frac{k_{NO_3}}{NO_3^-+k_{NO_3}} \cdot \frac{k_{SO_4}}{SO_4^{2-}+k_{SO_4}} \cdot Corg
\end{array}
$$
Units of the rates:
-
$AeroMin$, $Denitri$, $Sred$ and $Methanogen$ are in $mol~C~m^{-3}_S~d^{-1}$.
-
$Nitri$ is in $mol~N~m^{-3}_L~d^{-1}$.
-
$Sox$ is in $mol~S~m^{-3}_L~d^{-1}$.
Note that $AeroMin+Denitri+Sred+Methanogen$ is equal to $rMIN \cdot Corg$, as follows directly from the fact that the sum of the limitation and inhibition terms for a specific substrate is 1. Thus, the total mineralisation rate is first-order with respect to $Corg$ and independent of the terminal electron acceptor. The limitation and inhibition terms only dictate the contributions of the different oxidants to the overall mineralisation process.
Mass balance equations
The updated mass balance equations that include all processes mentioned are:\footnote{Note the factors $\phi$ or $1-\phi$, which cannot be canceled out from the spatial derivative because the porosity varies with depth.}
$$
\begin{array}{l}
\frac{\partial Corg}{\partial t} = \frac{1}{1-\phi}\cdot \frac{\partial}{\partial x}\left[(1-\phi)\cdot D_b\cdot \frac{ \partial Corg}{\partial x}\right]
-\frac{1}{1-\phi}\cdot \frac{ \partial }{\partial x}\left[(1-\phi)\cdot v\cdot Corg\right]
- AeroMin - Denitri - Sred - Methanogen \[4mm]
\frac{\partial O_2}{\partial t} = \frac{1}{\phi}\cdot \frac{\partial}{\partial x}\left[\phi \cdot D_{O2}\cdot \frac{ \partial O_2}{\partial x}\right]
-\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot O_2\right] - AeroMin \cdot f_{2L} - 2 \cdot Nitri - 2 \cdot Sox \[4mm]
\frac{\partial NH_3}{\partial t} = \frac{1}{\phi}\cdot \frac{\partial}{\partial x}\left[\phi \cdot D_{NH_3}\cdot \frac{ \partial NH_3}{\partial x}\right]
-\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot NH_3\right]
+ \frac{16}{106} \cdot (AeroMin + Denitri + Sred + Methanogen) \cdot f_{2L} - Nitri \[4mm]
\frac{\partial NO_3^-}{\partial t} = \frac{1}{\phi} \cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{NO_3^-}\cdot \frac{ \partial NO_3^-}{\partial x}\right]
-\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot NO_3^-\right]- 0.8 \cdot Denitri \cdot f_{2L} + Nitri \[4mm]
\frac{\partial SO_4^{2-}}{\partial t} = \frac{1}{\phi}\cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{SO_4^{2-}}\cdot \frac{ \partial SO_4^{2-}}{\partial x}\right]
-\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot SO_4^{2-}\right]
- 0.5 \cdot Sred\cdot f_{2L} + Sox \[4mm]
\frac{\partial H_2S}{\partial t} = \frac{1}{\phi}\cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{H_2S}\cdot \frac{ \partial H_2S}{\partial x}\right]
-\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot H_2S\right]
+ 0.5\cdot Sred\cdot f_{2L}- Sox \[4mm]
\frac{\partial CO_2}{\partial t} = \frac{1}{\phi} \cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{CO_2}\cdot \frac{ \partial CO_2}{\partial x}\right]
-\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot CO_2\right]
+ (OxicMin + Denitri + Sred + 0.5\cdot Methanogen)\cdot f_{2L} \[4mm]
\frac{\partial CH_4}{\partial t} = \frac{1}{\phi} \cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{CH_4}\cdot \frac{ \partial CH_4}{\partial x}\right]
-\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot CH_4\right]
+ 0.5\cdot Methanogen\cdot f_{2L}
\end{array}
$$
where $f_{2L}=(1-\phi)/\phi$.
R implementation
require(marelac)
require(ReacTran)
The model grid and associated properties
# units: time=days, space=meters, amount=moles, concentration=mol/m3
# spatial domain: total length of 0.10 m, 500 boxes
Length <- 0.1 # [m]
N <- 500
# grid with an exponentially increasing grid size, first cell 0.05 cm (=5e-4 m)
Grid <- setup.grid.1D(L = Length, N = N, dx.1 = 0.05e-4)
# function describing the variation of porosity (volume fraction of LIQUID) with depth
porFun.L <- function(x, por.SWI, por.deep, porcoef)
return( por.deep + (por.SWI-por.deep)*exp(-x*porcoef) )
# function describing the SOLID volume fraction (svf = 1-porosity)
porFun.S <- function(x, por.SWI, por.deep, porcoef)
return( 1-porFun.L(x, por.SWI, por.deep, porcoef) )
# calculate porosity and svf on the grid (mid-points and box interfaces, etc.)
porLiquid <- setup.prop.1D(func=porFun.L, grid=Grid, por.SWI=0.9, por.deep=0.7, porcoef=100)
porSolid <- setup.prop.1D(func=porFun.S, grid=Grid, por.SWI=0.9, por.deep=0.7, porcoef=100)
# molecular diffusion coefficients (m2/d) from the marelac package
diff.O2 <- diffcoeff(S=35, t=20)$O2 * 3600*24
diff.CO2 <- diffcoeff(S=35, t=20)$CO2 * 3600*24
diff.CH4 <- diffcoeff(S=35, t=20)$CH4 * 3600*24
diff.NH3 <- diffcoeff(S=35, t=20)$NH3 * 3600*24
diff.NO3 <- diffcoeff(S=35, t=20)$NO3 * 3600*24
diff.SO4 <- diffcoeff(S=35, t=20)$SO4 * 3600*24
diff.H2S <- diffcoeff(S=35, t=20)$H2S * 3600*24
# effective diffusion coefficients in the sediment (corrected for tortuosity)
porInt <- porLiquid$int # porosity at the box interfaces
diffO2 <- diff.O2 /(1-log(porInt^2))
diffCO2 <- diff.CO2/(1-log(porInt^2))
diffCH4 <- diff.CH4/(1-log(porInt^2))
diffNO3 <- diff.NO3/(1-log(porInt^2))
diffNH3 <- diff.NH3/(1-log(porInt^2))
diffSO4 <- diff.SO4/(1-log(porInt^2))
diffH2S <- diff.H2S/(1-log(porInt^2))
Other model parameters
pars <- c(
Dbio = 5e-4/365, # [m2/d] bioturbation mixing coefficient
v_adv = 5e-6, # [m/d] advection velocity = sediment accretion rate
rMin = 0.005, # [/d] POC mineralisation rate constant
depoPOC = 1e-3, # [mol/m2/d] POC deposition rate at SWI (sediment-water interface)
rNit = 10, # [/d] Nitrification rate constant
rSox = 10, # [/d] H2S reoxidation rate constant
kO2 = 0.001, # [mol/m3] half-saturation O2 concentration
kNO3 = 0.0001, # [mol/m3] half-saturation NO3 concentration
kSO4 = 0.1, # [mol/m3] half-saturation SO4 concentration
bwO2 = 0.300, # [mol/m3] O2 concentration at SWI
bwCO2 = 2, # [mol/m3] CO2 concentration at SWI
bwNO3 = 0.010, # [mol/m3] NO3 concentration at SWI
bwNH3 = 0.001, # [mol/m3] NH3 concentration at SWI
bwSO4 = 28, # [mol/m3] SO4 concentration at SWI
bwH2S = 0, # [mol/m3] H2S concentration at SWI
bwCH4 = 0 # [mol/m3] CH4 concentration at SWI
)
Definition and initialisation of state variables
names <- c("POC", "CO2", "O2", "NH3", "NO3", "SO4", "H2S", "CH4")
nspec <- length(names)
POC.ini <- rep(0, length = N) # initial conditions
CO2.ini <- rep(0, length = N)
O2.ini <- rep(0, length = N)
NH3.ini <- rep(0, length = N)
NO3.ini <- rep(0, length = N)
SO4.ini <- rep(0, length = N)
H2S.ini <- rep(0, length = N)
CH4.ini <- rep(0, length = N)
state.ini <- c(POC.ini, CO2.ini, O2.ini, NH3.ini, NO3.ini, SO4.ini, H2S.ini, CH4.ini)
Definition of the model function
Diamodel <- function (t, state, parms) # state is a LONG vector, at time t
{
with (as.list(parms),{
# unpack state variables
POC <- state[(0*N+1):(1*N)] # first N elements: POC
CO2 <- state[(1*N+1):(2*N)] # next N elements: CO2
O2 <- state[(2*N+1):(3*N)] # next N elements: O2
NH3 <- state[(3*N+1):(4*N)] # next N elements: NH3
NO3 <- state[(4*N+1):(5*N)] # next N elements: NO3
SO4 <- state[(5*N+1):(6*N)] # next N elements: SO4
H2S <- state[(6*N+1):(7*N)] # next N elements: H2S
CH4 <- state[(7*N+1):(8*N)] # next N elements: CH4
# === transport rates ===
# note: zero gradient by default at lower boundaries
# solid substances, VF = solid volume fraction = 1-porosity!
tran.POC <- tran.1D(C = POC, flux.up = depoPOC, # upper boundary: flux
dx = Grid, VF = porSolid, # grid and volume fraction (1-por)
D = Dbio, v = v_adv) # mixing (bioturbation) and advection
# dissolved substances, VF = liquid volume fraction = porosity!
tran.CO2 <- tran.1D(C = CO2, C.up = bwCO2, # upper boundary: concentration
dx = Grid, VF = porLiquid, # grid and volume fraction (por)
D = diffCO2, v = v_adv) # diffusive mixing and advection
tran.O2 <- tran.1D(C = O2, C.up = bwO2,
dx = Grid, VF = porLiquid,
D = diffO2, v = v_adv)
tran.NH3 <- tran.1D(C = NH3, C.up = bwNH3,
dx = Grid, VF = porLiquid,
D = diffNH3, v = v_adv)
tran.NO3 <- tran.1D(C = NO3, C.up = bwNO3,
dx = Grid, VF = porLiquid,
D = diffNO3, v = v_adv)
tran.SO4 <- tran.1D(C = SO4, C.up = bwSO4,
dx = Grid, VF = porLiquid,
D = diffSO4, v = v_adv)
tran.H2S <- tran.1D(C = H2S, C.up = bwH2S,
dx = Grid, VF = porLiquid,
D = diffH2S, v = v_adv)
tran.CH4 <- tran.1D(C = CH4, C.up = bwCH4,
dx = Grid, VF = porLiquid,
D = diffCH4, v = v_adv)
# === reaction rates ===
# [mol/m3 SOLID/d] (per volume of solid)
AeroMin <- rMin * O2/(O2+kO2) * POC
Denitri <- rMin * kO2/(O2+kO2) * NO3/(NO3+kNO3) * POC
SO4reduction <- rMin * kO2/(O2+kO2) * kNO3/(NO3+kNO3) * SO4/(SO4+kSO4) * POC
Methanogen <- rMin * kO2/(O2+kO2) * kNO3/(NO3+kNO3) * kSO4/(SO4+kSO4) * POC
Mineralisation <- AeroMin + Denitri + SO4reduction + Methanogen
# [mol/m3 LIQUID/d] (per volume of liquid)
Nitri <- rNit * O2/(O2+kO2) * NH3
Sox <- rSox * O2/(O2+kO2) * H2S
# === mass balances : dC/dt = transport + reactions ===
# solid substances [mol/m3 SOLID/d]
dPOC.dt <- ( tran.POC$dC # transport
- AeroMin - Denitri - SO4reduction - Methanogen) # reactions,
# dissolved substances [mol/m3 LIQUID/d]
poro <- porLiquid$mid
f2Liquid <- (1-poro)/poro # to convert from /solid to /liquid
dCO2.dt <- ( tran.CO2$dC
+ (AeroMin + Denitri + SO4reduction + 0.5*Methanogen)*f2Liquid )
dO2.dt <- ( tran.O2$dC - 2*Nitri - 2*Sox - AeroMin*f2Liquid )
dNH3.dt <- ( tran.NH3$dC - Nitri
+ 16/106*(AeroMin + Denitri + SO4reduction + Methanogen)*f2Liquid )
dNO3.dt <- ( tran.NO3$dC + Nitri - 4/5*Denitri*f2Liquid )
dSO4.dt <- ( tran.SO4$dC - 0.5*SO4reduction*f2Liquid + Sox)
dH2S.dt <- ( tran.H2S$dC + 0.5*SO4reduction*f2Liquid - Sox)
dCH4.dt <- ( tran.CH4$dC + 0.5*Methanogen*f2Liquid)
# === depth-integrated rates: [mol/m2 BULK/d]
TotalAero <- sum(AeroMin * Grid$dx * porSolid$mid)
TotalDenit <- sum(Denitri * Grid$dx * porSolid$mid)
TotalSO4red <- sum(SO4reduction * Grid$dx * porSolid$mid)
TotalCH4gen <- sum(Methanogen * Grid$dx * porSolid$mid)
TotalMin <- sum(Mineralisation* Grid$dx * porSolid$mid)
TotalNit <- sum(Nitri * Grid$dx * porLiquid$mid)
TotalSox <- sum(Sox * Grid$dx * porLiquid$mid)
return(list(c(dPOC.dt, dCO2.dt, dO2.dt, dNH3.dt, dNO3.dt,
dSO4.dt, dH2S.dt, dCH4.dt), # time-derivatives
AeroMin = AeroMin, # rate profiles
Denitri = Denitri,
SO4reduction = SO4reduction,
Methanogen = Methanogen,
Mineralisation = Mineralisation,
Nitri = Nitri,
Sox = Sox,
# part of mineralisation due to aero, denit, SO4red, CH4genesis
pAero = TotalAero/TotalMin,
pDenit = TotalDenit/TotalMin,
pSO4red = TotalSO4red/TotalMin,
pCH4gen = TotalCH4gen/TotalMin,
# for creating budgets - all in [mol/m2 BULK/d]
TotalAero = TotalAero,
TotalDenit = TotalDenit,
TotalSO4red = TotalSO4red,
TotalCH4gen = TotalCH4gen,
TotalMin = TotalMin,
TotalNit = TotalNit,
TotalSox = TotalSox,
CO2.SWI.Flux = tran.CO2$flux.up,
CO2.Deep.Flux = tran.CO2$flux.down,
O2.SWI.Flux = tran.O2$flux.up,
O2.Deep.Flux = tran.O2$flux.down,
NH3.SWI.Flux = tran.NH3$flux.up,
NH3.Deep.Flux = tran.NH3$flux.down,
NO3.SWI.Flux = tran.NO3$flux.up,
NO3.Deep.Flux = tran.NO3$flux.down,
SO4.SWI.Flux = tran.SO4$flux.up,
SO4.Deep.Flux = tran.SO4$flux.down,
H2S.SWI.Flux = tran.H2S$flux.up,
H2S.Deep.Flux = tran.H2S$flux.down,
CH4.SWI.Flux = tran.CH4$flux.up,
CH4.Deep.Flux = tran.CH4$flux.down,
POC.SWI.Flux = tran.POC$flux.up,
POC.Deep.Flux = tran.POC$flux.down))
})
}
Model application
Steady-state solution
Here, we find steady state solutions of the model for the three POC fluxes.
std <- steady.1D (y=state.ini, func=Diamodel, parms=pars,
nspec=nspec, dimens=N, names=names,
positive = TRUE) # to have only positive values!
p1 <- pars
p1["depoPOC"] <- 20e-3 # [mol/m2/d] POC deposition rate (flux at SWI)
std1 <- steady.1D (y=state.ini, func=Diamodel, parms=p1,
nspec=nspec, dimens=N, names=names, positive = TRUE)
p2 <- pars
p2["depoPOC"] <- 200e-3 # [mol/m2/d] POC deposition rate (flux at SWI)
std2 <- steady.1D (y=state.ini, func=Diamodel, parms=p2,
nspec=nspec, dimens=N, names=names, positive = TRUE)
Plotting
First, we plot depth profiles of the concentrations of state variables, all runs in one graph.
plot(std, std1, std2,
grid=Grid$x.mid, xyswap=TRUE, ylim=c(Length,0), lty=1, lwd=2, las=1,
ylab="depth (m)", xlab=c("mol/m3 Solid", rep("mol/m3 Liquid",nspec-1)))
Next, we plot depth profiles of the process rates, all runs in one graph.
plot(std, std1, std2,
which = c("AeroMin","Denitri","SO4reduction","Methanogen","Mineralisation","Nitri","Sox"),
grid=Grid$x.mid, xyswap=TRUE, ylab="depth (m)", ylim=c(Length,0), lty=1, lwd=2, las=1,
xlab=c(rep("molC/m3 Solid/d",5), "molN/m3 Liquid", "molS/m3 Liquid"))
Finally, we plot depth profiles of all mineralisation rates in one graph, separately for each model run. The graphs show that the organic matter is mineralised by aerobic mineralisation at the very top of the sediment (where $O_2$ is available), then by denitrification in a narrow band (where $NO_3^-$ is available but $O_2$ is depleted), then by sulphate reduction in the deeper sediment (where $SO_4^{2-}$ is available but $O_2$ and $NO_3^-$ are depleted).\footnote{Note that the latter two processes do not occur in the run with the lowest POC deposition flux.} Finally, methanogenesis occurs at depths where $O_2$, $NO_3^-$ and $SO_4^{2-}$ have been depleted and organic matter is still available (only visible in the run with the largest POC deposition flux). The plots thus illustrate the so-called biogeochemical cascade.
# all process rates in one graph
par(mfrow=c(1,3))
s <- std
matplot(y=Grid$x.mid, x=cbind(s$AeroMin, s$Denitri, s$SO4reduction, s$Methanogen),
type="l", ylim=c(Length,0), col=1:4, lwd=2, lty=1,
xlab="molC/m3 Solid/d", ylab="depth (m)")
s <- std1
matplot(y=Grid$x.mid, x=cbind(s$AeroMin, s$Denitri, s$SO4reduction, s$Methanogen),
type="l", ylim=c(Length,0), col=1:4, lwd=2, lty=1,
xlab="molC/m3 Solid/d", ylab="depth (m)", main="biogeochemical cascade")
s <- std2
matplot(y=Grid$x.mid, x=cbind(s$AeroMin, s$Denitri, s$SO4reduction, s$Methanogen),
type="l", ylim=c(Length,0), col=1:4, lwd=2, lty=1,
xlab="molC/m3 Solid/d", ylab="depth (m)")
legend("bottomright",legend=c("AeroMin", "Denitri", "SO4reduction", "Methanogen"),
lwd=2, col=1:4, lty=1, bty="n")
Contributions of mineralisation pathways
The following table shows that for the lowest POC deposition flux (std), organic matter mineralisation is primarily aerobic (almost 99% of total mineralisation). For the second largest POC deposition flux (std1), about 80% of organic matter is mineralised via sulphate reduction, while aerobic mineralisation and denitrification contribute by about 15% and 6%, respectively. Methanogenesis is significant (about 16%) only for the largest POC deposition flux (std2). Note that the depth-integrated total mineralisation (last line, values in $mol~C~m^{-2}~d^{-1}$) is equal to the POC deposition flux, showing that at the bottom of the modelled sediment column, all deposited POC was mineralised and converted to either $CO_2$ or $CH_4$.
toselect <- c("pAero", "pDenit", "pSO4red", "pCH4gen", "TotalMin")
PATH <- std[toselect]
PATH1 <- std1[toselect]
PATH2 <- std2[toselect]
CONTRIB <- data.frame(std=unlist(PATH), std1=unlist(PATH1), std2=unlist(PATH2))
knitr::kable(CONTRIB, digits = 4)
Oxygen budget
The following table ($O_2$ fluxes in $mol~O_2~m^{-2}~yr^{-1}$) shows that for the lowest POC deposition flux (std), about 2/3 of the sedimentary $O_2$ consumption occurs due to the aerobic mineralisation of the organic matter, while the remaining 1/3 is due to nitrification. In contrast, for the larger POC deposition fluxes, sedimentary $O_2$ consumption occurs primarily (50% or more) due to oxidation of the free sulphide produced by sulphate reduction, while the contributions of aerobic mineralisation and nitrification are lower.
toselect <- c("TotalAero", "TotalNit", "TotalSox", "O2.SWI.Flux", "O2.Deep.Flux")
fac <- 365 # to molO2/m2/yr
BUDGET <- std[toselect]
BUDGET1 <- std1[toselect]
BUDGET2 <- std2[toselect]
O2budget <- rbind(std=unlist(BUDGET), std1=unlist(BUDGET1), std2=unlist(BUDGET2))*fac
O2BUDG <-
data.frame(O2FluxIn = O2budget[,"O2.SWI.Flux"],
O2FluxOut = -O2budget[,"O2.Deep.Flux"],
O2consByAeroMin = -O2budget[,"TotalAero"],
O2consByNitri = -O2budget[,"TotalNit"]*2, # 2 moles per mole of NH3 nitrified
O2consBySox = -O2budget[,"TotalSox"]*2) # 2 moles per mole of H2S oxidised
knitr::kable(O2BUDG, digits = 2)
The output below demonstrates that the difference between the $O_2$ fluxes at the domain boundaries and the sum of all depth-integrated rates of $O_2$ consumption is equal to (nearly) zero, consistent with the mass conservation law.
rowSums(O2BUDG)
Sulphur budget
The following table ($S$ fluxes in $mol~S~m^{-2}~yr^{-1}$) shows that if sulphate reduction is significant, the total amount of $SO_4^{2-}$ consumed in the sediment by sulphate reduction is always greater than the amount of $SO_4^{2-}$ that enters the sediment from the overlying water. This is only possible if $SO_4^{2-}$ is additionally produced within the sediment, i.e., if the free sulphide produced by sulphate reduction is recycled back to sulphate. The process responsible for this recycling is the aerobic oxidation of sulphide, which, as noted in the previous section, is a major contributor to the sedimentary oxygen consumption.
toselect <- c("TotalSO4red", "TotalSox", "SO4.SWI.Flux", "SO4.Deep.Flux",
"H2S.SWI.Flux", "H2S.Deep.Flux")
fac <- 365 # to molS/m2/yr
BUDGET <- std[toselect]
BUDGET1 <- std1[toselect]
BUDGET2 <- std2[toselect]
SBUDG <- rbind(std=unlist(BUDGET), std1=unlist(BUDGET1), std2=unlist(BUDGET2))*fac
SBUDG[,"TotalSO4red"] <- SBUDG[,"TotalSO4red"]*0.5 # 0.5 moles of SO4 per mole of C
knitr::kable(SBUDG, digits = 4)
This conclusion is confirmed by the following output, which shows that the net flux of sulphate consumed in the sediment (which is equal to SO4.SWI.Flux - SO4.Deep.Flux) together with the gross flux of sulphide production (TotalSox) is equal to the gross (total) rate of sulphate reduction in the sediment (TotalSO4red).
Scheck <- cbind(SBUDG[,"SO4.SWI.Flux"]-SBUDG[,"SO4.Deep.Flux"] # net SO4 consumed
+ SBUDG[,"TotalSox"], # + H2S recycled
SBUDG[,"TotalSO4red"]) # gross SO4 consumed
Scheck
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)
Answers
Rate expressions
$$\begin{array}{l} AeroMin = r_{Min} \cdot \frac{O_2}{O_2+k_{O_2}} \cdot Corg \[5mm] Nitri = r_{Nit} \cdot \frac{O_2}{O_2+k_{O_2}} \cdot NH_3 \[5mm] Denitri = r_{Min} \cdot \frac{k_{O_2}}{O_2+k_{O_2}} \cdot \frac{NO_3^-}{NO_3^-+k_{NO_3}} \cdot Corg \[5mm] Sred = r_{Min} \cdot \frac{k_{O_2}}{O_2+k_{O_2}} \cdot \frac{k_{NO_3}}{NO_3^-+k_{NO_3}} \cdot \frac{SO_4^{2-}}{SO_4^{2-}+k_{SO_4}} \cdot Corg \[5mm] Sox = r_{Sox} \cdot \frac{O_2}{O_2+k_{O_2}} \cdot H_2S \[5mm] Methanogen = r_{Min} \cdot \frac{k_{O_2}}{O_2+k_{O_2}} \cdot \frac{k_{NO_3}}{NO_3^-+k_{NO_3}} \cdot \frac{k_{SO_4}}{SO_4^{2-}+k_{SO_4}} \cdot Corg \end{array} $$
Units of the rates:
-
$AeroMin$, $Denitri$, $Sred$ and $Methanogen$ are in $mol~C~m^{-3}_S~d^{-1}$.
-
$Nitri$ is in $mol~N~m^{-3}_L~d^{-1}$.
-
$Sox$ is in $mol~S~m^{-3}_L~d^{-1}$.
Note that $AeroMin+Denitri+Sred+Methanogen$ is equal to $rMIN \cdot Corg$, as follows directly from the fact that the sum of the limitation and inhibition terms for a specific substrate is 1. Thus, the total mineralisation rate is first-order with respect to $Corg$ and independent of the terminal electron acceptor. The limitation and inhibition terms only dictate the contributions of the different oxidants to the overall mineralisation process.
Mass balance equations
The updated mass balance equations that include all processes mentioned are:\footnote{Note the factors $\phi$ or $1-\phi$, which cannot be canceled out from the spatial derivative because the porosity varies with depth.} $$ \begin{array}{l} \frac{\partial Corg}{\partial t} = \frac{1}{1-\phi}\cdot \frac{\partial}{\partial x}\left[(1-\phi)\cdot D_b\cdot \frac{ \partial Corg}{\partial x}\right] -\frac{1}{1-\phi}\cdot \frac{ \partial }{\partial x}\left[(1-\phi)\cdot v\cdot Corg\right] - AeroMin - Denitri - Sred - Methanogen \[4mm] \frac{\partial O_2}{\partial t} = \frac{1}{\phi}\cdot \frac{\partial}{\partial x}\left[\phi \cdot D_{O2}\cdot \frac{ \partial O_2}{\partial x}\right] -\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot O_2\right] - AeroMin \cdot f_{2L} - 2 \cdot Nitri - 2 \cdot Sox \[4mm] \frac{\partial NH_3}{\partial t} = \frac{1}{\phi}\cdot \frac{\partial}{\partial x}\left[\phi \cdot D_{NH_3}\cdot \frac{ \partial NH_3}{\partial x}\right] -\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot NH_3\right] + \frac{16}{106} \cdot (AeroMin + Denitri + Sred + Methanogen) \cdot f_{2L} - Nitri \[4mm] \frac{\partial NO_3^-}{\partial t} = \frac{1}{\phi} \cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{NO_3^-}\cdot \frac{ \partial NO_3^-}{\partial x}\right] -\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot NO_3^-\right]- 0.8 \cdot Denitri \cdot f_{2L} + Nitri \[4mm] \frac{\partial SO_4^{2-}}{\partial t} = \frac{1}{\phi}\cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{SO_4^{2-}}\cdot \frac{ \partial SO_4^{2-}}{\partial x}\right] -\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot SO_4^{2-}\right] - 0.5 \cdot Sred\cdot f_{2L} + Sox \[4mm] \frac{\partial H_2S}{\partial t} = \frac{1}{\phi}\cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{H_2S}\cdot \frac{ \partial H_2S}{\partial x}\right] -\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot H_2S\right] + 0.5\cdot Sred\cdot f_{2L}- Sox \[4mm] \frac{\partial CO_2}{\partial t} = \frac{1}{\phi} \cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{CO_2}\cdot \frac{ \partial CO_2}{\partial x}\right] -\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot CO_2\right] + (OxicMin + Denitri + Sred + 0.5\cdot Methanogen)\cdot f_{2L} \[4mm] \frac{\partial CH_4}{\partial t} = \frac{1}{\phi} \cdot \frac{\partial}{\partial x}\left[\phi\cdot D_{CH_4}\cdot \frac{ \partial CH_4}{\partial x}\right] -\frac{1}{\phi}\cdot \frac{ \partial }{\partial x}\left[\phi\cdot v\cdot CH_4\right] + 0.5\cdot Methanogen\cdot f_{2L} \end{array} $$ where $f_{2L}=(1-\phi)/\phi$.
R implementation
require(marelac) require(ReacTran)
The model grid and associated properties
# units: time=days, space=meters, amount=moles, concentration=mol/m3 # spatial domain: total length of 0.10 m, 500 boxes Length <- 0.1 # [m] N <- 500 # grid with an exponentially increasing grid size, first cell 0.05 cm (=5e-4 m) Grid <- setup.grid.1D(L = Length, N = N, dx.1 = 0.05e-4) # function describing the variation of porosity (volume fraction of LIQUID) with depth porFun.L <- function(x, por.SWI, por.deep, porcoef) return( por.deep + (por.SWI-por.deep)*exp(-x*porcoef) ) # function describing the SOLID volume fraction (svf = 1-porosity) porFun.S <- function(x, por.SWI, por.deep, porcoef) return( 1-porFun.L(x, por.SWI, por.deep, porcoef) ) # calculate porosity and svf on the grid (mid-points and box interfaces, etc.) porLiquid <- setup.prop.1D(func=porFun.L, grid=Grid, por.SWI=0.9, por.deep=0.7, porcoef=100) porSolid <- setup.prop.1D(func=porFun.S, grid=Grid, por.SWI=0.9, por.deep=0.7, porcoef=100) # molecular diffusion coefficients (m2/d) from the marelac package diff.O2 <- diffcoeff(S=35, t=20)$O2 * 3600*24 diff.CO2 <- diffcoeff(S=35, t=20)$CO2 * 3600*24 diff.CH4 <- diffcoeff(S=35, t=20)$CH4 * 3600*24 diff.NH3 <- diffcoeff(S=35, t=20)$NH3 * 3600*24 diff.NO3 <- diffcoeff(S=35, t=20)$NO3 * 3600*24 diff.SO4 <- diffcoeff(S=35, t=20)$SO4 * 3600*24 diff.H2S <- diffcoeff(S=35, t=20)$H2S * 3600*24 # effective diffusion coefficients in the sediment (corrected for tortuosity) porInt <- porLiquid$int # porosity at the box interfaces diffO2 <- diff.O2 /(1-log(porInt^2)) diffCO2 <- diff.CO2/(1-log(porInt^2)) diffCH4 <- diff.CH4/(1-log(porInt^2)) diffNO3 <- diff.NO3/(1-log(porInt^2)) diffNH3 <- diff.NH3/(1-log(porInt^2)) diffSO4 <- diff.SO4/(1-log(porInt^2)) diffH2S <- diff.H2S/(1-log(porInt^2))
Other model parameters
pars <- c( Dbio = 5e-4/365, # [m2/d] bioturbation mixing coefficient v_adv = 5e-6, # [m/d] advection velocity = sediment accretion rate rMin = 0.005, # [/d] POC mineralisation rate constant depoPOC = 1e-3, # [mol/m2/d] POC deposition rate at SWI (sediment-water interface) rNit = 10, # [/d] Nitrification rate constant rSox = 10, # [/d] H2S reoxidation rate constant kO2 = 0.001, # [mol/m3] half-saturation O2 concentration kNO3 = 0.0001, # [mol/m3] half-saturation NO3 concentration kSO4 = 0.1, # [mol/m3] half-saturation SO4 concentration bwO2 = 0.300, # [mol/m3] O2 concentration at SWI bwCO2 = 2, # [mol/m3] CO2 concentration at SWI bwNO3 = 0.010, # [mol/m3] NO3 concentration at SWI bwNH3 = 0.001, # [mol/m3] NH3 concentration at SWI bwSO4 = 28, # [mol/m3] SO4 concentration at SWI bwH2S = 0, # [mol/m3] H2S concentration at SWI bwCH4 = 0 # [mol/m3] CH4 concentration at SWI )
Definition and initialisation of state variables
names <- c("POC", "CO2", "O2", "NH3", "NO3", "SO4", "H2S", "CH4") nspec <- length(names) POC.ini <- rep(0, length = N) # initial conditions CO2.ini <- rep(0, length = N) O2.ini <- rep(0, length = N) NH3.ini <- rep(0, length = N) NO3.ini <- rep(0, length = N) SO4.ini <- rep(0, length = N) H2S.ini <- rep(0, length = N) CH4.ini <- rep(0, length = N) state.ini <- c(POC.ini, CO2.ini, O2.ini, NH3.ini, NO3.ini, SO4.ini, H2S.ini, CH4.ini)
Definition of the model function
Diamodel <- function (t, state, parms) # state is a LONG vector, at time t { with (as.list(parms),{ # unpack state variables POC <- state[(0*N+1):(1*N)] # first N elements: POC CO2 <- state[(1*N+1):(2*N)] # next N elements: CO2 O2 <- state[(2*N+1):(3*N)] # next N elements: O2 NH3 <- state[(3*N+1):(4*N)] # next N elements: NH3 NO3 <- state[(4*N+1):(5*N)] # next N elements: NO3 SO4 <- state[(5*N+1):(6*N)] # next N elements: SO4 H2S <- state[(6*N+1):(7*N)] # next N elements: H2S CH4 <- state[(7*N+1):(8*N)] # next N elements: CH4 # === transport rates === # note: zero gradient by default at lower boundaries # solid substances, VF = solid volume fraction = 1-porosity! tran.POC <- tran.1D(C = POC, flux.up = depoPOC, # upper boundary: flux dx = Grid, VF = porSolid, # grid and volume fraction (1-por) D = Dbio, v = v_adv) # mixing (bioturbation) and advection # dissolved substances, VF = liquid volume fraction = porosity! tran.CO2 <- tran.1D(C = CO2, C.up = bwCO2, # upper boundary: concentration dx = Grid, VF = porLiquid, # grid and volume fraction (por) D = diffCO2, v = v_adv) # diffusive mixing and advection tran.O2 <- tran.1D(C = O2, C.up = bwO2, dx = Grid, VF = porLiquid, D = diffO2, v = v_adv) tran.NH3 <- tran.1D(C = NH3, C.up = bwNH3, dx = Grid, VF = porLiquid, D = diffNH3, v = v_adv) tran.NO3 <- tran.1D(C = NO3, C.up = bwNO3, dx = Grid, VF = porLiquid, D = diffNO3, v = v_adv) tran.SO4 <- tran.1D(C = SO4, C.up = bwSO4, dx = Grid, VF = porLiquid, D = diffSO4, v = v_adv) tran.H2S <- tran.1D(C = H2S, C.up = bwH2S, dx = Grid, VF = porLiquid, D = diffH2S, v = v_adv) tran.CH4 <- tran.1D(C = CH4, C.up = bwCH4, dx = Grid, VF = porLiquid, D = diffCH4, v = v_adv) # === reaction rates === # [mol/m3 SOLID/d] (per volume of solid) AeroMin <- rMin * O2/(O2+kO2) * POC Denitri <- rMin * kO2/(O2+kO2) * NO3/(NO3+kNO3) * POC SO4reduction <- rMin * kO2/(O2+kO2) * kNO3/(NO3+kNO3) * SO4/(SO4+kSO4) * POC Methanogen <- rMin * kO2/(O2+kO2) * kNO3/(NO3+kNO3) * kSO4/(SO4+kSO4) * POC Mineralisation <- AeroMin + Denitri + SO4reduction + Methanogen # [mol/m3 LIQUID/d] (per volume of liquid) Nitri <- rNit * O2/(O2+kO2) * NH3 Sox <- rSox * O2/(O2+kO2) * H2S # === mass balances : dC/dt = transport + reactions === # solid substances [mol/m3 SOLID/d] dPOC.dt <- ( tran.POC$dC # transport - AeroMin - Denitri - SO4reduction - Methanogen) # reactions, # dissolved substances [mol/m3 LIQUID/d] poro <- porLiquid$mid f2Liquid <- (1-poro)/poro # to convert from /solid to /liquid dCO2.dt <- ( tran.CO2$dC + (AeroMin + Denitri + SO4reduction + 0.5*Methanogen)*f2Liquid ) dO2.dt <- ( tran.O2$dC - 2*Nitri - 2*Sox - AeroMin*f2Liquid ) dNH3.dt <- ( tran.NH3$dC - Nitri + 16/106*(AeroMin + Denitri + SO4reduction + Methanogen)*f2Liquid ) dNO3.dt <- ( tran.NO3$dC + Nitri - 4/5*Denitri*f2Liquid ) dSO4.dt <- ( tran.SO4$dC - 0.5*SO4reduction*f2Liquid + Sox) dH2S.dt <- ( tran.H2S$dC + 0.5*SO4reduction*f2Liquid - Sox) dCH4.dt <- ( tran.CH4$dC + 0.5*Methanogen*f2Liquid) # === depth-integrated rates: [mol/m2 BULK/d] TotalAero <- sum(AeroMin * Grid$dx * porSolid$mid) TotalDenit <- sum(Denitri * Grid$dx * porSolid$mid) TotalSO4red <- sum(SO4reduction * Grid$dx * porSolid$mid) TotalCH4gen <- sum(Methanogen * Grid$dx * porSolid$mid) TotalMin <- sum(Mineralisation* Grid$dx * porSolid$mid) TotalNit <- sum(Nitri * Grid$dx * porLiquid$mid) TotalSox <- sum(Sox * Grid$dx * porLiquid$mid) return(list(c(dPOC.dt, dCO2.dt, dO2.dt, dNH3.dt, dNO3.dt, dSO4.dt, dH2S.dt, dCH4.dt), # time-derivatives AeroMin = AeroMin, # rate profiles Denitri = Denitri, SO4reduction = SO4reduction, Methanogen = Methanogen, Mineralisation = Mineralisation, Nitri = Nitri, Sox = Sox, # part of mineralisation due to aero, denit, SO4red, CH4genesis pAero = TotalAero/TotalMin, pDenit = TotalDenit/TotalMin, pSO4red = TotalSO4red/TotalMin, pCH4gen = TotalCH4gen/TotalMin, # for creating budgets - all in [mol/m2 BULK/d] TotalAero = TotalAero, TotalDenit = TotalDenit, TotalSO4red = TotalSO4red, TotalCH4gen = TotalCH4gen, TotalMin = TotalMin, TotalNit = TotalNit, TotalSox = TotalSox, CO2.SWI.Flux = tran.CO2$flux.up, CO2.Deep.Flux = tran.CO2$flux.down, O2.SWI.Flux = tran.O2$flux.up, O2.Deep.Flux = tran.O2$flux.down, NH3.SWI.Flux = tran.NH3$flux.up, NH3.Deep.Flux = tran.NH3$flux.down, NO3.SWI.Flux = tran.NO3$flux.up, NO3.Deep.Flux = tran.NO3$flux.down, SO4.SWI.Flux = tran.SO4$flux.up, SO4.Deep.Flux = tran.SO4$flux.down, H2S.SWI.Flux = tran.H2S$flux.up, H2S.Deep.Flux = tran.H2S$flux.down, CH4.SWI.Flux = tran.CH4$flux.up, CH4.Deep.Flux = tran.CH4$flux.down, POC.SWI.Flux = tran.POC$flux.up, POC.Deep.Flux = tran.POC$flux.down)) }) }
Model application
Steady-state solution
Here, we find steady state solutions of the model for the three POC fluxes.
std <- steady.1D (y=state.ini, func=Diamodel, parms=pars, nspec=nspec, dimens=N, names=names, positive = TRUE) # to have only positive values! p1 <- pars p1["depoPOC"] <- 20e-3 # [mol/m2/d] POC deposition rate (flux at SWI) std1 <- steady.1D (y=state.ini, func=Diamodel, parms=p1, nspec=nspec, dimens=N, names=names, positive = TRUE) p2 <- pars p2["depoPOC"] <- 200e-3 # [mol/m2/d] POC deposition rate (flux at SWI) std2 <- steady.1D (y=state.ini, func=Diamodel, parms=p2, nspec=nspec, dimens=N, names=names, positive = TRUE)
Plotting
First, we plot depth profiles of the concentrations of state variables, all runs in one graph.
plot(std, std1, std2, grid=Grid$x.mid, xyswap=TRUE, ylim=c(Length,0), lty=1, lwd=2, las=1, ylab="depth (m)", xlab=c("mol/m3 Solid", rep("mol/m3 Liquid",nspec-1)))
Next, we plot depth profiles of the process rates, all runs in one graph.
plot(std, std1, std2, which = c("AeroMin","Denitri","SO4reduction","Methanogen","Mineralisation","Nitri","Sox"), grid=Grid$x.mid, xyswap=TRUE, ylab="depth (m)", ylim=c(Length,0), lty=1, lwd=2, las=1, xlab=c(rep("molC/m3 Solid/d",5), "molN/m3 Liquid", "molS/m3 Liquid"))
Finally, we plot depth profiles of all mineralisation rates in one graph, separately for each model run. The graphs show that the organic matter is mineralised by aerobic mineralisation at the very top of the sediment (where $O_2$ is available), then by denitrification in a narrow band (where $NO_3^-$ is available but $O_2$ is depleted), then by sulphate reduction in the deeper sediment (where $SO_4^{2-}$ is available but $O_2$ and $NO_3^-$ are depleted).\footnote{Note that the latter two processes do not occur in the run with the lowest POC deposition flux.} Finally, methanogenesis occurs at depths where $O_2$, $NO_3^-$ and $SO_4^{2-}$ have been depleted and organic matter is still available (only visible in the run with the largest POC deposition flux). The plots thus illustrate the so-called biogeochemical cascade.
# all process rates in one graph par(mfrow=c(1,3)) s <- std matplot(y=Grid$x.mid, x=cbind(s$AeroMin, s$Denitri, s$SO4reduction, s$Methanogen), type="l", ylim=c(Length,0), col=1:4, lwd=2, lty=1, xlab="molC/m3 Solid/d", ylab="depth (m)") s <- std1 matplot(y=Grid$x.mid, x=cbind(s$AeroMin, s$Denitri, s$SO4reduction, s$Methanogen), type="l", ylim=c(Length,0), col=1:4, lwd=2, lty=1, xlab="molC/m3 Solid/d", ylab="depth (m)", main="biogeochemical cascade") s <- std2 matplot(y=Grid$x.mid, x=cbind(s$AeroMin, s$Denitri, s$SO4reduction, s$Methanogen), type="l", ylim=c(Length,0), col=1:4, lwd=2, lty=1, xlab="molC/m3 Solid/d", ylab="depth (m)") legend("bottomright",legend=c("AeroMin", "Denitri", "SO4reduction", "Methanogen"), lwd=2, col=1:4, lty=1, bty="n")
Contributions of mineralisation pathways
The following table shows that for the lowest POC deposition flux (std), organic matter mineralisation is primarily aerobic (almost 99% of total mineralisation). For the second largest POC deposition flux (std1), about 80% of organic matter is mineralised via sulphate reduction, while aerobic mineralisation and denitrification contribute by about 15% and 6%, respectively. Methanogenesis is significant (about 16%) only for the largest POC deposition flux (std2). Note that the depth-integrated total mineralisation (last line, values in $mol~C~m^{-2}~d^{-1}$) is equal to the POC deposition flux, showing that at the bottom of the modelled sediment column, all deposited POC was mineralised and converted to either $CO_2$ or $CH_4$.
toselect <- c("pAero", "pDenit", "pSO4red", "pCH4gen", "TotalMin") PATH <- std[toselect] PATH1 <- std1[toselect] PATH2 <- std2[toselect] CONTRIB <- data.frame(std=unlist(PATH), std1=unlist(PATH1), std2=unlist(PATH2)) knitr::kable(CONTRIB, digits = 4)
Oxygen budget
The following table ($O_2$ fluxes in $mol~O_2~m^{-2}~yr^{-1}$) shows that for the lowest POC deposition flux (std), about 2/3 of the sedimentary $O_2$ consumption occurs due to the aerobic mineralisation of the organic matter, while the remaining 1/3 is due to nitrification. In contrast, for the larger POC deposition fluxes, sedimentary $O_2$ consumption occurs primarily (50% or more) due to oxidation of the free sulphide produced by sulphate reduction, while the contributions of aerobic mineralisation and nitrification are lower.
toselect <- c("TotalAero", "TotalNit", "TotalSox", "O2.SWI.Flux", "O2.Deep.Flux") fac <- 365 # to molO2/m2/yr BUDGET <- std[toselect] BUDGET1 <- std1[toselect] BUDGET2 <- std2[toselect] O2budget <- rbind(std=unlist(BUDGET), std1=unlist(BUDGET1), std2=unlist(BUDGET2))*fac O2BUDG <- data.frame(O2FluxIn = O2budget[,"O2.SWI.Flux"], O2FluxOut = -O2budget[,"O2.Deep.Flux"], O2consByAeroMin = -O2budget[,"TotalAero"], O2consByNitri = -O2budget[,"TotalNit"]*2, # 2 moles per mole of NH3 nitrified O2consBySox = -O2budget[,"TotalSox"]*2) # 2 moles per mole of H2S oxidised knitr::kable(O2BUDG, digits = 2)
The output below demonstrates that the difference between the $O_2$ fluxes at the domain boundaries and the sum of all depth-integrated rates of $O_2$ consumption is equal to (nearly) zero, consistent with the mass conservation law.
rowSums(O2BUDG)
Sulphur budget
The following table ($S$ fluxes in $mol~S~m^{-2}~yr^{-1}$) shows that if sulphate reduction is significant, the total amount of $SO_4^{2-}$ consumed in the sediment by sulphate reduction is always greater than the amount of $SO_4^{2-}$ that enters the sediment from the overlying water. This is only possible if $SO_4^{2-}$ is additionally produced within the sediment, i.e., if the free sulphide produced by sulphate reduction is recycled back to sulphate. The process responsible for this recycling is the aerobic oxidation of sulphide, which, as noted in the previous section, is a major contributor to the sedimentary oxygen consumption.
toselect <- c("TotalSO4red", "TotalSox", "SO4.SWI.Flux", "SO4.Deep.Flux", "H2S.SWI.Flux", "H2S.Deep.Flux") fac <- 365 # to molS/m2/yr BUDGET <- std[toselect] BUDGET1 <- std1[toselect] BUDGET2 <- std2[toselect] SBUDG <- rbind(std=unlist(BUDGET), std1=unlist(BUDGET1), std2=unlist(BUDGET2))*fac SBUDG[,"TotalSO4red"] <- SBUDG[,"TotalSO4red"]*0.5 # 0.5 moles of SO4 per mole of C knitr::kable(SBUDG, digits = 4)
This conclusion is confirmed by the following output, which shows that the net flux of sulphate consumed in the sediment (which is equal to SO4.SWI.Flux - SO4.Deep.Flux) together with the gross flux of sulphide production (TotalSox) is equal to the gross (total) rate of sulphate reduction in the sediment (TotalSO4red).
Scheck <- cbind(SBUDG[,"SO4.SWI.Flux"]-SBUDG[,"SO4.Deep.Flux"] # net SO4 consumed + SBUDG[,"TotalSox"], # + H2S recycled SBUDG[,"TotalSO4red"]) # gross SO4 consumed Scheck
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.