Exercise 5: Benthic Population Model of a River
ETH Zurich Course 701-0426-00L: Modelling Aquatic Ecosystems (Schuwirth)
April 30, 2025
Goals:
Understand the concepts and process formulations of the coupling of benthic population dynamics with nutrient and dissolved oxygen dynamics in the water column as described in section 11.6.
Understand the implementation of this model based on the R packages
stoichcalcandecosim.Understand the behaviour (time courses of the state variables and overall mass fluxes of this model).
Perform some simple sensitivity analyses to deepen the understanding of the behaviour of the model.
Task 1: Model formulation
Study and try to understand the model formulation given in section 11.6. Note that it may be useful to study first the ‘’intermediately complex’’ model described in section 11.5.
Here are the process table (Table 1) and process rates (Table 2):
Table 1: Process table of a model for the benthic population and oxygen and nutrient dynamics in a river.
\[ \begin{array}{|l|c|c|c|c|c|c|c|c|c|c|l|} \hline \textbf{Process} & \text{HPO}_4^{2-} & \text{NH}_4^+ & \text{NO}_3^- & \text{O}_2 & \text{ALG} & \text{ZOO} & \text{POMD} & \text{POMI} & \text{SPOMD} & \text{SPOMI} & \textbf{Rate} \\ & \text{(gP)} & \text{(gN)} & \text{(gN)} & \text{(gO)} & \text{(gDM)} & \text{(gDM)} & \text{(gDM)} & \text{(gDM)} & \text{(gDM)} & \text{(gDM)} & \\ \hline \text{Growth of algae NO}_3^- & - & & - & + & 1 & & & & & & \rho_{\text{gro,ALG,NO}_3^-} \\ \text{Growth of algae NH}_4^+ & - & - & & + & 1 & & & & & & \rho_{\text{gro,ALG,NH}_4^+} \\ \text{Respiration of algae} & + & + & & - & -1 & & & & & & \rho_{\text{resp,ALG}} \\ \text{Death of algae} & 0/+ & 0/+ & & 0/+ & -1 & & (1 - f_1)Y_{\text{ALG,death}} & f_1Y_{\text{ALG,death}} & & & \rho_{\text{death,ALG}} \\ \text{Growth of zooplankton} & + & + & & - & -\frac{1}{Y_{\text{ZOO}}} & 1 & \frac{(1 - f_1)f_e}{Y_{\text{ZOO}}} & \frac{f_1f_e}{Y_{\text{ZOO}}} & & & \rho_{\text{gro,ZOO}} \\ \text{Respiration of zoopl.} & & + & & - & & -1 & & & & & \rho_{\text{resp,ZOO}} \\ \text{Death of zooplankton} & 0/+ & 0/+ & & 0/+ & & -1 & (1 - f_1)Y_{\text{ZOO,death}} & f_1Y_{\text{ZOO,death}} & & & \rho_{\text{death,ZOO}} \\ \text{Nitrification} & & -1 & + & & & & & & & & \rho_{\text{nitri}} \\ \text{Oxic mineral. of org. part.} & + & + & & + & & & -1 & & & & \rho_{\text{miner,ox,POMD}} \\ \text{Ox. min. in sediment} & + & + & & & & & & & -1 & & \rho_{\text{miner,ox,SPOMD}} \\ \text{Anox. min. in sediment} & + & + & & - & & & & & & -1 & \rho_{\text{miner,anox,SPOMD}} \\ \text{Sed. of deg. org. part.} & & & & & & & -1 & & 1 & & \rho_{\text{sed,POMD}} \\ \text{Sed. of inert org. part.} & & & & & & & & -1 & & 1 & \rho_{\text{sed,POMI}} \\ \hline \end{array} \]
Table 2: Process rates of the model for benthic population and oxygen and nutrient dynamics in a river \[ \scriptsize \begin{array}{|l|l|} \hline \textbf{Rate} & \textbf{Rate expression} \\ \hline \rho_{\text{gro,SALG,NH}_4^+} & \begin{aligned} &k_{\text{gro,ALG},T_0} \cdot \exp\left(\beta_{\text{ALG}}(T - T_0)\right) \cdot \frac{I_0 \exp(-\lambda h)}{K_I + I_0 \exp(-\lambda h)} \cdot \\ &\min\left( \frac{C_{\text{HPO}_4^{2-}}}{K_{\text{HPO}_4^{2-},\text{ALG}} + C_{\text{HPO}_4^{2-}}}, \frac{C_{\text{NH}_4^+} + C_{\text{NO}_3^-}}{K_{N,\text{ALG}} + C_{\text{NH}_4^+} + C_{\text{NO}_3^-}} \right) \cdot \\ &\frac{p_{\text{NH}_4^+,\text{ALG}} C_{\text{NH}_4^+}}{p_{\text{NH}_4^+,\text{ALG}} C_{\text{NH}_4^+} + K_{\text{sha,SALG}}} \cdot D_{\text{SALG}} \end{aligned} \\ \hline \rho_{\text{gro,SALG,NO}_3^-} & \begin{aligned} &k_{\text{gro,ALG},T_0} \cdot \exp\left(\beta_{\text{ALG}}(T - T_0)\right) \cdot \frac{I_0 \exp(-\lambda h)}{K_I + I_0 \exp(-\lambda h)} \cdot \\ &\min\left( \frac{C_{\text{HPO}_4^{2-}}}{K_{\text{HPO}_4^{2-},\text{ALG}} + C_{\text{HPO}_4^{2-}}}, \frac{C_{\text{NH}_4^+} + C_{\text{NO}_3^-}}{K_{N,\text{ALG}} + C_{\text{NH}_4^+} + C_{\text{NO}_3^-}} \right) \cdot \\ &\frac{C_{\text{NO}_3^-}}{p_{\text{NH}_4^+,\text{ALG}} C_{\text{NH}_4^+} + C_{\text{NO}_3^-} + K_{\text{sha,SALG}}} \cdot D_{\text{SALG}} \end{aligned} \\ \hline \rho_{\text{resp,SALG}} & k_{\text{resp,ALG},T_0} \cdot \exp\left(\beta_{\text{ALG}}(T - T_0)\right) \cdot \frac{C_{\text{O}_2}}{K_{\text{O}_2,\text{ALG}} + C_{\text{O}_2}} \cdot D_{\text{SALG}} \\ \hline \rho_{\text{death,SALG}} & k_{\text{death,ALG}} \cdot D_{\text{SALG}} \\ \hline \rho_{\text{gro,SHET,NH}_4^+} & \begin{aligned} &k_{\text{gro,HET},T_0} \cdot \exp\left(\beta_{\text{HET}}(T - T_0)\right) \cdot \frac{p_{\text{NH}_4^+,\text{HET}} C_{\text{NH}_4^+}}{p_{\text{NH}_4^+,\text{HET}} C_{\text{NH}_4^+} + C_{\text{NO}_3^-}} \cdot \\ &\min\left( \frac{C_{\text{DOM}}}{K_{\text{DOM,HET}} + C_{\text{DOM}}}, \frac{C_{\text{O}_2}}{K_{\text{O}_2,\text{HET}} + C_{\text{O}_2}}, \frac{C_{\text{HPO}_4^{2-}}}{K_{\text{HPO}_4^{2-},\text{HET}} + C_{\text{HPO}_4^{2-}}}, \frac{C_{\text{NH}_4^+} + C_{\text{NO}_3^-}}{K_{N,\text{HET}} + C_{\text{NH}_4^+} + C_{\text{NO}_3^-}} \right) \cdot \\ &\frac{K_{\text{lim,SHET}}}{K_{\text{lim,SHET}} + D_{\text{SHET}}} \cdot D_{\text{SHET}} \end{aligned} \\ \hline \rho_{\text{gro,SHET,NO}_3^-} & \begin{aligned} &k_{\text{gro,HET},T_0} \cdot \exp\left(\beta_{\text{HET}}(T - T_0)\right) \cdot \frac{C_{\text{NO}_3^-}}{p_{\text{NH}_4^+,\text{HET}} C_{\text{NH}_4^+} + C_{\text{NO}_3^-}} \cdot \\ &\min\left( \frac{C_{\text{DOM}}}{K_{\text{DOM,HET}} + C_{\text{DOM}}}, \frac{C_{\text{O}_2}}{K_{\text{O}_2,\text{HET}} + C_{\text{O}_2}}, \frac{C_{\text{HPO}_4^{2-}}}{K_{\text{HPO}_4^{2-},\text{HET}} + C_{\text{HPO}_4^{2-}}}, \frac{C_{\text{NH}_4^+} + C_{\text{NO}_3^-}}{K_{N,\text{HET}} + C_{\text{NH}_4^+} + C_{\text{NO}_3^-}} \right) \cdot \\ &\frac{K_{\text{lim,SHET}}}{K_{\text{lim,SHET}} + D_{\text{SHET}}} \cdot D_{\text{SHET}} \end{aligned} \\ \hline \rho_{\text{resp,SHET}} & k_{\text{resp,HET},T_0} \cdot \exp\left(\beta_{\text{HET}}(T - T_0)\right) \cdot \frac{C_{\text{O}_2}}{K_{\text{O}_2,\text{HET}} + C_{\text{O}_2}} \cdot D_{\text{SHET}} \\ \hline \rho_{\text{death,SHET}} & k_{\text{death,HET}} \cdot D_{\text{SHET}} \\ \hline \rho_{\text{gro,SN1}} & \begin{aligned} &k_{\text{gro,N1},T_0} \cdot \exp\left(\beta_{N1}(T - T_0)\right) \cdot \\ &\min\left( \frac{C_{\text{NH}_4^+}}{K_{\text{NH}_4^+,\text{nitri}} + C_{\text{NH}_4^+}}, \frac{C_{\text{O}_2}}{K_{\text{O}_2,\text{nitri}} + C_{\text{O}_2}}, \frac{C_{\text{HPO}_4^{2-}}}{K_{\text{HPO}_4^{2-},\text{nitri}} + C_{\text{HPO}_4^{2-}}} \right) \cdot \\ &\frac{K_{\text{lim,SN}}}{K_{\text{lim,SN}} + D_{\text{SN1}}} \cdot D_{\text{SN1}} \end{aligned} \\ \hline \rho_{\text{resp,SN1}} & k_{\text{resp,N1},T_0} \cdot \exp\left(\beta_{N1}(T - T_0)\right) \cdot \frac{C_{\text{O}_2}}{K_{\text{O}_2,\text{nitri}} + C_{\text{O}_2}} \cdot D_{\text{SN1}} \\ \hline \rho_{\text{death,SN1}} & k_{\text{death,N1}} \cdot D_{\text{SN1}} \\ \hline \rho_{\text{gro,SN2}} & \begin{aligned} &k_{\text{gro,N2},T_0} \cdot \exp\left(\beta_{N2}(T - T_0)\right) \cdot \\ &\min\left( \frac{C_{\text{NO}_3^-}}{K_{\text{NO}_3^-,\text{nitri}} + C_{\text{NO}_3^-}}, \frac{C_{\text{O}_2}}{K_{\text{O}_2,\text{nitri}} + C_{\text{O}_2}}, \frac{C_{\text{HPO}_4^{2-}}}{K_{\text{HPO}_4^{2-},\text{nitri}} + C_{\text{HPO}_4^{2-}}} \right) \cdot \\ &\frac{K_{\text{lim,SN}}}{K_{\text{lim,SN}} + D_{\text{SN2}}} \cdot D_{\text{SN2}} \end{aligned} \\ \hline \rho_{\text{resp,SN2}} & k_{\text{resp,N2},T_0} \cdot \exp\left(\beta_{N2}(T - T_0)\right) \cdot \frac{C_{\text{O}_2}}{K_{\text{O}_2,\text{nitri}} + C_{\text{O}_2}} \cdot D_{\text{SN2}} \\ \hline \rho_{\text{death,SN2}} & k_{\text{death,N2}} \cdot D_{\text{SN2}} \\ \hline \rho_{\text{hyd}} & k_{\text{hyd,SPOM},T_0} \cdot \exp\left(\beta_{\text{hyd}}(T - T_0)\right) \cdot D_{\text{SPOM}} \\ \hline \end{array} \] ### Task 2: Model implementation
Study, complete and run the implementation of the model given below.
Load the package ecosim (and install if it is not yet
installed).
# load required packages:
if ( !require("ecosim") ) { install.packages("ecosim"); library(ecosim) }## Loading required package: ecosim## Loading required package: deSolve## Loading required package: stoichcalcNote that the packages deSolve, to numerially solve
systems of ordinary differential equations, and stoichcalc,
to calculate process stoichiometry were loaded automatically.
Define parameter values for stoichiometry, kinetics, river geometry input and environmental conditions:
# Definition of parameters:
param <- list(
# stoichiometric parameters
alpha.O.ALG = 0.50, # gO/gALG
alpha.H.ALG = 0.07, # gH/gALG
alpha.N.ALG = 0.06, # gN/gALG
alpha.P.ALG = 0.01, # gP/gALG
alpha.O.BAC = 0.30, # gO/gBAC
alpha.H.BAC = 0.08, # gH/gBAC
alpha.N.BAC = 0.10, # gN/gBAC
alpha.P.BAC = 0.02, # gP/gBAC
alpha.O.POM = 0.26, # gO/gPOM
alpha.H.POM = 0.07, # gH/gPOM
alpha.N.POM = 0.04, # gN/gPOM
alpha.P.POM = 0.01, # gP/gPOM
alpha.O.DOM = 0.26, # gO/gDOM
alpha.H.DOM = 0.07, # gH/gDOM
alpha.N.DOM = 0.04, # gN/gDOM
alpha.P.DOM = 0.01, # gP/gDOM
Y.N1 = 0.13, # gN1/gNH4-N
Y.N2 = 0.03, # gN2/gNO2-N
Y.HET = 0.6, # gHET/gDOM
Y.hyd = 1.0, # gDOM/gPOM
# kinetic parameters
k.gro.ALG = 1.5, # 1/d
k.gro.HET = 1.5, # 1/d
k.gro.N1 = 0.8, # 1/d
k.gro.N2 = 1.1, # 1/d
k.resp.ALG = 0.10, # 1/d
k.resp.HET = 0.20, # 1/d
k.resp.N1 = 0.10, # 1/d
k.resp.N2 = 0.10, # 1/d
k.death.ALG = 0.10, # 1/d
k.death.HET = 0.20, # 1/d
k.death.N1 = 0.10, # 1/d
k.death.N2 = 0.10, # 1/d
k.hyd.POM = 0.5, # 1/d
K.HPO4.ALG = 0.01, # gP/m3
K.HPO4.HET = 0.01, # gP/m3
K.HPO4.nitri = 0.02, # gP/m3
K.N.ALG = 0.04, # gN/m3
K.N.HET = 0.04, # gN/m3
p.NH4.ALG = 5, # -
p.NH4.HET = 5, # -
K.NH4.nitri = 0.5, # gN/m3
K.NO2.nitri = 0.5, # gN/m3
K.O2.nitri = 0.5, # gO/m3
K.O2.ALG = 0.5, # gO/m3
K.O2.HET = 0.5, # gO/m3
K.DOM.HET = 0.5, # gDOM/m3
beta.ALG = 0.046, # 1/degC
beta.HET = 0.046, # 1/degC
beta.N1 = 0.098, # 1/degC
beta.N2 = 0.069, # 1/degC
beta.hyd = 0.046, # 1/degC
K.I = 200, # W/m2
lambda = 0.1, # 1/m
K2.O2 = 10, # 1/d
K.shadow.ALG = 100, # gDM/m2
K.limit.HET = 10, # gDM/m2
K.limit.nitri = 5, # gDM/m2
# River geometry
L = 2000, # m
w = 20, # m
h = 0.5, # m
# Input and initial conditions
Q.in = 4, # m3/s
D.ALG.ini = 50, # gDM/m2
D.HET.ini = 20, # gDM/m2
D.N1.ini = 2, # gDM/m2
D.N2.ini = 1, # gDM/m2
D.POM.ini = 50, # gDM/m2
C.HPO4.ini = 0.4, # gP/m3
C.NH4.ini = 0.4, # gN/m3
C.NO3.ini = 4, # gN/m3
C.NO2.ini = 0, # gN/m3
C.O2.ini = 10, # gO/m3
C.DOM.ini = 3, # gDOM/m3
C.HPO4.in = 0.4, # gP/m3
C.NH4.in = 0.4, # gN/m3
C.NO3.in = 4, # gN/m3
C.O2.in = 10, # gO/m3
C.DOM.in = 3, # gDOM/m3
# environmental conditions
T0 = 20, # degC
T.min = 15, # degC
T.max = 25, # degC
I0.max = 600, # W/m2
t.max.I = 0.5, # d
t.max.T = 0.6, # d
p = 101325) # PaThis set of parameters contains the elemental mass fractions of nitrogen, phosphorus, oxygen and hydrogen. The largest mass fraction is then calculated as the difference to unity:
# choose carbon fractions to guarantee that the fractions sum to unity:
param$alpha.C.ALG <- 1 - (param$alpha.O.ALG + param$alpha.H.ALG +
param$alpha.N.ALG + param$alpha.P.ALG)
param$alpha.C.BAC <- 1 - (param$alpha.O.BAC + param$alpha.H.BAC +
param$alpha.N.BAC + param$alpha.P.BAC)
param$alpha.C.POM <- 1 - (param$alpha.O.POM + param$alpha.H.POM +
param$alpha.N.POM + param$alpha.P.POM)
param$alpha.C.DOM <- 1 - (param$alpha.O.DOM + param$alpha.H.DOM +
param$alpha.N.DOM + param$alpha.P.DOM)We define the ‘’yields’’ of algae and bacteria death processes as in exercise 4:
# choose yield of death to guarantee that no nutrients are required
# (oxygen content of POM was reduced to avoid need of oxygen):
param$Y.ALG.death <- min(1,
param$alpha.N.ALG/param$alpha.N.POM,
param$alpha.P.ALG/param$alpha.P.POM,
param$alpha.C.ALG/param$alpha.C.POM)
param$Y.BAC.death <- min(1,
param$alpha.N.BAC/param$alpha.N.POM,
param$alpha.P.BAC/param$alpha.P.POM,
param$alpha.C.BAC/param$alpha.C.POM)
print(param$Y.ALG.death)## [1] 0.5806452print(param$Y.BAC.death)## [1] 0.8064516Once all substances are defined, we can construct a composition
matrix by passing the list of substances to the function
calc.comp.matrix of the package
stoichcalc:
# Construction of composition matrix:
NH4 <- c(H = 4*1/14, # gH/gNH4-N
N = 1, # gN/gNH4-N
charge = 1/14) # chargeunits/gNH4-N
NO2 <- c(O = 2*16/14, # gO/gNO2-N
N = 1, # gN/gNO2-N
charge = -1/14) # chargeunits/gNO2-N
NO3 <- c(O = 3*16/14, # gO/gNO3-N
N = 1, # gN/gNO3-N
charge = -1/14) # chargeunits/gNO3-N
HPO4 <- c(O = 4*16/31, # gO/gHPO4-P
H = 1*1/31, # gH/gHPO4-P
P = 1, # gP/gHPO4-P
charge = -2/31) # chargeunits/gHPO4-P
HCO3 <- c(C = 1, # gC/gHCO3-C
O = 3*16/12, # gO/gHCO3-C
H = 1*1/12, # gH/gHCO3-C
charge = -1/12) # chargeunits/gHCO3-C
O2 <- c(O = 1) # gO/gO2-O
H <- c(H = 1, # gH/molH
charge = 1) # chargeunits/molH
H2O <- c(O = 1*16, # gO/molH2O
H = 2*1) # gH/molH2O
ALG <- c(C = param$alpha.C.ALG, # gC/gALG
O = param$alpha.O.ALG, # gO/gALG
H = param$alpha.H.ALG, # gH/gALG
N = param$alpha.N.ALG, # gN/gALG
P = param$alpha.P.ALG) # gP/gALG
HET <- c(C = param$alpha.C.BAC, # gC/gBAC
O = param$alpha.O.BAC, # gO/gBAC
H = param$alpha.H.BAC, # gH/gBAC
N = param$alpha.N.BAC, # gN/gBAC
P = param$alpha.P.BAC) # gP/gBAC
N1 <- c(C = param$alpha.C.BAC, # gC/gBAC
O = param$alpha.O.BAC, # gO/gBAC
H = param$alpha.H.BAC, # gH/gBAC
N = param$alpha.N.BAC, # gN/gBAC
P = param$alpha.P.BAC) # gP/gBAC
N2 <- c(C = param$alpha.C.BAC, # gC/gBAC
O = param$alpha.O.BAC, # gO/gBAC
H = param$alpha.H.BAC, # gH/gBAC
N = param$alpha.N.BAC, # gN/gBAC
P = param$alpha.P.BAC) # gP/gBAC
POM <- c(C = param$alpha.C.POM, # gC/gPOM
O = param$alpha.O.POM, # gO/gPOM
H = param$alpha.H.POM, # gH/gPOM
N = param$alpha.N.POM, # gN/gPOM
P = param$alpha.P.POM) # gP/gPOM
DOM <- c(C = param$alpha.C.DOM, # gC/gDOM
O = param$alpha.O.DOM, # gO/gDOM
H = param$alpha.H.DOM, # gH/gDOM
N = param$alpha.N.DOM, # gN/gDOM
P = param$alpha.P.DOM) # gP/gDOM
subst.comp <- list(C.NH4 = NH4,
C.NO2 = NO2,
C.NO3 = NO3,
C.HPO4 = HPO4,
C.HCO3 = HCO3,
C.O2 = O2,
C.DOM = DOM,
C.H = H,
C.H2O = H2O,
D.ALG = ALG,
D.HET = HET,
D.N1 = N1,
D.N2 = N2,
D.POM = POM)
alpha <- calc.comp.matrix(subst.comp)## [1] "Composition matrix (6x14) successfully constructed:"
## [1] "el. const.: H,N,charge,O,P,C"
## [1] "substances: C.NH4,C.NO2,C.NO3,C.HPO4,C.HCO3,C.O2,C.DOM,C.H,C.H2O,D.ALG,D.HET,D.N1,D.N2,D.POM"print(round(alpha,3))## C.NH4 C.NO2 C.NO3 C.HPO4 C.HCO3 C.O2 C.DOM C.H C.H2O D.ALG D.HET D.N1
## H 0.286 0.000 0.000 0.032 0.083 0 0.07 1 2 0.07 0.08 0.08
## N 1.000 1.000 1.000 0.000 0.000 0 0.04 0 0 0.06 0.10 0.10
## charge 0.071 -0.071 -0.071 -0.065 -0.083 0 0.00 1 0 0.00 0.00 0.00
## O 0.000 2.286 3.429 2.065 4.000 1 0.26 0 16 0.50 0.30 0.30
## P 0.000 0.000 0.000 1.000 0.000 0 0.01 0 0 0.01 0.02 0.02
## C 0.000 0.000 0.000 0.000 1.000 0 0.62 0 0 0.36 0.50 0.50
## D.N2 D.POM
## H 0.08 0.07
## N 0.10 0.04
## charge 0.00 0.00
## O 0.30 0.26
## P 0.02 0.01
## C 0.50 0.62Having defined the composition matrix and the stoichiometric parameters, we can define the stoichiometries of the processes to be considered by the model:
# Derivation of Process Stoichiometry:
# Growth of algae on ammonium:
nu.gro.ALG.NH4 <-
calc.stoich.coef(alpha = alpha,
name = "gro.ALG.NH4",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.ALG"),
subst.norm = "D.ALG",
nu.norm = 1)## [1] "Number of substances: 7"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 0"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Growth of algae on nitrate:
nu.gro.ALG.NO3 <-
calc.stoich.coef(alpha = alpha,
name = "gro.ALG.NO3",
subst = c("C.NO3","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.ALG"),
subst.norm = "D.ALG",
nu.norm = 1)## [1] "Number of substances: 7"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 0"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Respiration of algae:
nu.resp.ALG <-
calc.stoich.coef(alpha = alpha,
name = "resp.ALG",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.ALG"),
subst.norm = "D.ALG",
nu.norm = -1)## [1] "Number of substances: 7"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 0"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Death of algae:
nu.death.ALG <-
calc.stoich.coef(alpha = alpha,
name = "death.ALG",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.ALG","D.POM"),
subst.norm = "D.ALG",
nu.norm = -1,
constraints = list(c("D.ALG" = param$Y.ALG.death,
"D.POM" = 1)))## [1] "Number of substances: 8"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 1"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Growth of heterotrophic bacteria using ammonium:
nu.gro.HET.NH4 <-
calc.stoich.coef(alpha = alpha,
name = "gro.HET.NH4",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.DOM","C.H","C.H2O","D.HET"),
subst.norm = "D.HET",
nu.norm = 1,
constraints = list(c("C.DOM" = param$Y.HET,
"D.HET" = 1)))## [1] "Number of substances: 8"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 1"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Growth of heterotrophic bacteria using nitrate:
nu.gro.HET.NO3 <-
calc.stoich.coef(alpha = alpha,
name = "gro.HET.NO3",
subst = c("C.NO3","C.HPO4","C.HCO3","C.O2",
"C.DOM","C.H","C.H2O","D.HET"),
subst.norm = "D.HET",
nu.norm = 1,
constraints = list(c("C.DOM" = param$Y.HET,
"D.HET" = 1)))## [1] "Number of substances: 8"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 1"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Respiration of heterotrophic bacteria:
nu.resp.HET <-
calc.stoich.coef(alpha = alpha,
name = "resp.HET",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.HET"),
subst.norm = "D.HET",
nu.norm = -1)## [1] "Number of substances: 7"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 0"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Death of heterotrophic bacteria:
nu.death.HET <-
calc.stoich.coef(alpha = alpha,
name = "death.HET",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.HET","D.POM"),
subst.norm = "D.HET",
nu.norm = -1,
constraints = list(c("D.HET" = param$Y.BAC.death,
"D.POM" = 1)))## [1] "Number of substances: 8"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 1"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Growth of first stage nitrifiers:
nu.gro.N1 <-
calc.stoich.coef(alpha = alpha,
name = "gro.N1",
subst = c("C.NH4","C.NO2","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.N1"),
subst.norm = "D.N1",
nu.norm = 1,
constraints = list(c("C.NH4" = param$Y.N1,
"D.N1" = 1)))## [1] "Number of substances: 8"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 1"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Respiration of first stage nitrifiers:
nu.resp.N1 <-
calc.stoich.coef(alpha = alpha,
name = "resp.N1",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.N1"),
subst.norm = "D.N1",
nu.norm = -1)## [1] "Number of substances: 7"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 0"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Death of first stage nitrifiers:
nu.death.N1 <-
calc.stoich.coef(alpha = alpha,
name = "death.N1",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.N1","D.POM"),
subst.norm = "D.N1",
nu.norm = -1,
constraints = list(c("D.N1" = param$Y.BAC.death,
"D.POM" = 1)))## [1] "Number of substances: 8"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 1"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Growth of second stage nitrifiers:
nu.gro.N2 <-
calc.stoich.coef(alpha = alpha,
name = "gro.N2",
subst = c("C.NO2","C.NO3","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.N2"),
subst.norm = "D.N2",
nu.norm = 1,
constraints = list(c("C.NO2" = param$Y.N2,
"D.N2" = 1)))## [1] "Number of substances: 8"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 1"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Respiration of second stage nitrifiers:
nu.resp.N2 <-
calc.stoich.coef(alpha = alpha,
name = "resp.N2",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.N2"),
subst.norm = "D.N2",
nu.norm = -1)## [1] "Number of substances: 7"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 0"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Death of second stage nitrifiers:
nu.death.N2 <-
calc.stoich.coef(alpha = alpha,
name = "death.N2",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.H","C.H2O","D.N2","D.POM"),
subst.norm = "D.N2",
nu.norm = -1,
constraints = list(c("D.N2" = param$Y.BAC.death,
"D.POM" = 1)))## [1] "Number of substances: 8"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 1"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."# Hydrolysis of organic particles:
nu.hyd.POM <-
calc.stoich.coef(alpha = alpha,
name = "hyd.POM",
subst = c("C.NH4","C.HPO4","C.HCO3","C.O2",
"C.DOM","C.H","C.H2O","D.POM"),
subst.norm = "D.POM",
nu.norm = -1,
constraints = list(c("D.POM" = param$Y.hyd,
"C.DOM" = 1)))## [1] "Number of substances: 8"
## [1] "Number of elementary constituents: 6"
## [1] "Number of constraints: 1"
## [1] "Number of independent processes: 1"
## [1] "Stoichiometric coefficients successfully calculated."Finally, we bind the stoichiometries together for the full stoichiometric matrix:
# Combine process stoichiometries to stoichiometric matrix:
nu <- rbind(nu.gro.ALG.NH4,
nu.gro.ALG.NO3,
nu.resp.ALG,
nu.death.ALG,
nu.gro.HET.NH4,
nu.gro.HET.NO3,
nu.resp.HET,
nu.death.HET,
nu.gro.N1,
nu.resp.N1,
nu.death.N1,
nu.gro.N2,
nu.resp.N2,
nu.death.N2,
nu.hyd.POM)
print(round(nu,3))## C.NH4 C.NO2 C.NO3 C.HPO4 C.HCO3 C.O2 C.DOM C.H C.H2O
## gro.ALG.NH4 -0.060 0.000 0.000 -0.010 -0.360 0.930 0.000 -0.026 0.002
## gro.ALG.NO3 0.000 0.000 -0.060 -0.010 -0.360 1.204 0.000 -0.035 -0.002
## resp.ALG 0.060 0.000 0.000 0.010 0.360 -0.930 0.000 0.026 -0.002
## death.ALG 0.037 0.000 0.000 0.004 0.000 0.172 0.000 -0.002 0.011
## gro.HET.NH4 -0.033 0.000 0.000 -0.003 0.533 -1.635 -1.667 0.047 -0.022
## gro.HET.NO3 0.000 0.000 -0.033 -0.003 0.533 -1.483 -1.667 0.042 -0.025
## resp.HET 0.100 0.000 0.000 0.020 0.500 -1.528 0.000 0.036 -0.013
## death.HET 0.068 0.000 0.000 0.012 0.000 0.003 0.000 -0.004 0.004
## gro.N1 -7.692 7.592 0.000 -0.020 -0.500 -24.503 0.000 1.049 0.556
## resp.N1 0.100 0.000 0.000 0.020 0.500 -1.528 0.000 0.036 -0.013
## death.N1 0.068 0.000 0.000 0.012 0.000 0.003 0.000 -0.004 0.004
## gro.N2 0.000 -33.333 33.233 -0.020 -0.500 -36.110 0.000 -0.050 0.006
## resp.N2 0.100 0.000 0.000 0.020 0.500 -1.528 0.000 0.036 -0.013
## death.N2 0.068 0.000 0.000 0.012 0.000 0.003 0.000 -0.004 0.004
## hyd.POM 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000
## D.ALG D.HET D.N1 D.N2 D.POM
## gro.ALG.NH4 1 0 0 0 0.000
## gro.ALG.NO3 1 0 0 0 0.000
## resp.ALG -1 0 0 0 0.000
## death.ALG -1 0 0 0 0.581
## gro.HET.NH4 0 1 0 0 0.000
## gro.HET.NO3 0 1 0 0 0.000
## resp.HET 0 -1 0 0 0.000
## death.HET 0 -1 0 0 0.806
## gro.N1 0 0 1 0 0.000
## resp.N1 0 0 -1 0 0.000
## death.N1 0 0 -1 0 0.806
## gro.N2 0 0 0 1 0.000
## resp.N2 0 0 0 -1 0.000
## death.N2 0 0 0 -1 0.806
## hyd.POM 0 0 0 0 -1.000# write.table(nu,file="river2_nu.dat",sep="\t",col.names=NA)Now, we can proceed with the definition of transformation processes:
# Definition of transformation processes:
# Growth of algae on ammonium:
gro.ALG.NH4 <-
new(Class = "process",
name = "gro.ALG.NH4",
rate = expression(k.gro.ALG
*exp(beta.ALG*(T-T0))
*I0*exp(-lambda*h)/(K.I+I0*exp(-lambda*h))
*min(C.HPO4/(K.HPO4.ALG+C.HPO4),
(C.NH4+C.NO3)/(K.N.ALG+C.NH4+C.NO3))
*(p.NH4.ALG*C.NH4/(p.NH4.ALG*C.NH4+C.NO3))
*D.ALG*K.shadow.ALG/(K.shadow.ALG+D.ALG)),
stoich = as.list(nu["gro.ALG.NH4",]),
pervol = F)
# Growth of algae on nitrate:
gro.ALG.NO3 <-
new(Class = "process",
name = "gro.ALG.NO3",
rate = expression(k.gro.ALG
*exp(beta.ALG*(T-T0))
*I0*exp(-lambda*h)/(K.I+I0*exp(-lambda*h))
*min(C.HPO4/(K.HPO4.ALG+C.HPO4),
(C.NH4+C.NO3)/(K.N.ALG+C.NH4+C.NO3))
*(C.NO3/(p.NH4.ALG*C.NH4+C.NO3))
*D.ALG*K.shadow.ALG/(K.shadow.ALG+D.ALG)),
stoich = as.list(nu["gro.ALG.NO3",]),
pervol = F)
# Respiration of algae:
resp.ALG <-
new(Class = "process",
name = "resp.ALG",
rate = expression(k.resp.ALG*exp(beta.ALG*(T-T0))
*C.O2/(K.O2.ALG+C.O2)*D.ALG),
stoich = as.list(nu["resp.ALG",]),
pervol = F)
# Death of algae:
death.ALG <-
new(Class = "process",
name = "death.ALG",
rate = expression(k.death.ALG*D.ALG),
stoich = as.list(nu["death.ALG",]),
pervol = F)
# Growth of heterotrophic bacteria with ammonium:
gro.HET.NH4 <-
new(Class = "process",
name = "gro.HET.NH4",
rate = expression(k.gro.HET
*exp(beta.HET*(T-T0))
*min(C.DOM/(K.DOM.HET+C.DOM),
C.O2/(K.O2.HET+C.O2),
C.HPO4/(K.HPO4.HET+C.HPO4),
(C.NH4+C.NO3)/(K.N.HET+C.NH4+C.NO3))
*(p.NH4.HET*C.NH4/(p.NH4.HET*C.NH4+C.NO3))
*D.HET*K.limit.HET/(K.limit.HET+D.HET)),
stoich = as.list(nu["gro.HET.NH4",]),
pervol = F)
# Growth of heterotrophic bacteria with nitrate:
gro.HET.NO3 <-
new(Class = "process",
name = "gro.HET.NO3",
rate = expression(k.gro.HET
*exp(beta.HET*(T-T0))
*min(C.DOM/(K.DOM.HET+C.DOM),
C.O2/(K.O2.HET+C.O2),
C.HPO4/(K.HPO4.HET+C.HPO4),
(C.NH4+C.NO3)/(K.N.HET+C.NH4+C.NO3))
*(C.NO3/(p.NH4.HET*C.NH4+C.NO3))
*D.HET*K.limit.HET/(K.limit.HET+D.HET)),
stoich = as.list(nu["gro.HET.NO3",]),
pervol = F)
# Respiration of heterotrophic bacteria:
resp.HET <-
new(Class = "process",
name = "resp.HET",
rate = expression(k.resp.HET*exp(beta.HET*(T-T0))
*C.O2/(K.O2.HET+C.O2)*D.HET),
stoich = as.list(nu["resp.HET",]),
pervol = F)
# Death of heterotrophic bacteria:
death.HET <-
new(Class = "process",
name = "death.HET",
rate = expression(k.death.HET*D.HET),
stoich = as.list(nu["death.HET",]),
pervol = F)
# Growth of first stage nitrifiers:
gro.N1 <-
new(Class = "process",
name = "gro.N1",
rate = expression(k.gro.N1
*exp(beta.N1*(T-T0))
*min(C.HPO4/(K.HPO4.nitri+C.HPO4),
C.NH4/(K.NH4.nitri+C.NH4),
C.O2/(K.O2.nitri+C.O2))
*D.N1*K.limit.nitri/(K.limit.nitri+D.N1)),
stoich = as.list(nu["gro.N1",]),
pervol = F)
# Respiration of first stage nitrifiers:
resp.N1 <-
new(Class = "process",
name = "resp.N1",
rate = expression(k.resp.N1*exp(beta.N1*(T-T0))
*C.O2/(K.O2.nitri+C.O2)*D.N1),
stoich = as.list(nu["resp.N1",]),
pervol = F)
# Death of first stage nitrifiers:
death.N1 <-
new(Class = "process",
name = "death.N1",
rate = expression(k.death.N1*D.N1),
stoich = as.list(nu["death.N1",]),
pervol = F)
# Growth of second stage nitrifiers:
gro.N2 <-
new(Class = "process",
name = "gro.N2",
rate = expression(k.gro.N2
*exp(beta.N2*(T-T0))
*min(C.HPO4/(K.HPO4.nitri+C.HPO4),
C.NO2/(K.NO2.nitri+C.NO2),
C.O2/(K.O2.nitri+C.O2))
*D.N2*K.limit.nitri/(K.limit.nitri+D.N2)),
stoich = as.list(nu["gro.N2",]),
pervol = F)
# Respiration of second stage nitrifiers:
resp.N2 <-
new(Class = "process",
name = "resp.N2",
rate = expression(k.resp.N2*exp(beta.N2*(T-T0))
*C.O2/(K.O2.nitri+C.O2)*D.N2),
stoich = as.list(nu["resp.N2",]),
pervol = F)
# Death of second stage nitrifiers:
death.N2 <-
new(Class = "process",
name = "death.N2",
rate = expression(k.death.N2*D.N2),
stoich = as.list(nu["death.N2",]),
pervol = F)
# Hydrolysis of POM:
hyd.POM <-
new(Class = "process",
name = "nu.hyd.POM",
rate = expression(k.hyd.POM*exp(beta.hyd*(T-T0))
*D.POM),
stoich = as.list(nu["hyd.POM",]),
pervol = F)We now define the seasonally varying light intensity and temperature and, as a consequence of temperature variation, the seasonally varying saturation concentration of dissolved oxygen:
# Definition of environmental conditions:
cond <- list(I0 = expression(I0.max*0.5*(sign(cos(2*pi/1.0*(t-t.max.I)))+1)
*cos(2*pi/1.0*(t-t.max.I))^2),
T = expression( 0.5*(T.min+T.max)
+0.5*(T.max-T.min)
*cos(2*pi/1.0*(t-t.max.T))),
C.O2.sat = expression(exp(7.7117-1.31403*log(T+45.93))
*p/101325))We then define three reactors to approximate advective transport and roughly resolve longitudinal concentration gradients:
# Definition of reactor:
# River Sections:
reach1 <-
new(Class = "reactor",
name = "R1",
volume.ini = expression(L*w*h),
area = expression(L*w),
conc.pervol.ini = list(C.HPO4 = expression(C.HPO4.ini), # gP/m3
C.NH4 = expression(C.NH4.ini), # gN/m3
C.NO2 = expression(C.NO2.ini), # gN/m3
C.NO3 = expression(C.NO3.ini), # gN/m3
C.O2 = expression(C.O2.ini), # gN/m3
C.DOM = expression(C.DOM.ini)), # gDOM/m3
conc.perarea.ini = list(D.ALG = expression(D.ALG.ini), # gALG/m2
D.HET = expression(D.HET.ini), # gHET/m2
D.N1 = expression(D.N1.ini), # gN1/m2
D.N2 = expression(D.N2.ini), # gN2/m2
D.POM = expression(D.POM.ini)), # gPOM/m2
input = list(C.O2 = expression(K2.O2*L*w*h
*(C.O2.sat-C.O2))), # gas exchange
inflow = expression(Q.in*86400), # m3/d
inflow.conc = list(C.HPO4 = expression(C.HPO4.in),
C.NH4 = expression(C.NH4.in),
C.NO3 = expression(C.NO3.in),
C.O2 = expression(C.O2.sat),
C.DOM = expression(C.DOM.in)),
processes = list(gro.ALG.NH4,gro.ALG.NO3,resp.ALG,death.ALG,
gro.HET.NH4,gro.HET.NO3,resp.HET,death.HET,
gro.N1,resp.N1,death.N1,gro.N2,resp.N2,death.N2,
hyd.POM))
reach2 <-
new(Class = "reactor",
name = "R2",
volume.ini = expression(L*w*h),
area = expression(L*w),
conc.pervol.ini = list(C.HPO4 = expression(C.HPO4.ini), # gP/m3
C.NH4 = expression(C.NH4.ini), # gN/m3
C.NO2 = expression(C.NO2.ini), # gN/m3
C.NO3 = expression(C.NO3.ini), # gN/m3
C.O2 = expression(C.O2.ini), # gN/m3
C.DOM = expression(C.DOM.ini)), # gDOM/m3
conc.perarea.ini = list(D.ALG = expression(D.ALG.ini), # gALG/m2
D.HET = expression(D.HET.ini), # gHET/m2
D.N1 = expression(D.N1.ini), # gN1/m2
D.N2 = expression(D.N2.ini), # gN2/m2
D.POM = expression(D.POM.ini)), # gPOM/m2
input = list(C.O2 = expression(K2.O2*L*w*h
*(C.O2.sat-C.O2))), # gas exchange
processes = list(gro.ALG.NH4,gro.ALG.NO3,resp.ALG,death.ALG,
gro.HET.NH4,gro.HET.NO3,resp.HET,death.HET,
gro.N1,resp.N1,death.N1,gro.N2,resp.N2,death.N2,
hyd.POM))
reach3 <-
new(Class = "reactor",
name = "R3",
volume.ini = expression(L*w*h),
area = expression(L*w),
conc.pervol.ini = list(C.HPO4 = expression(C.HPO4.ini), # gP/m3
C.NH4 = expression(C.NH4.ini), # gN/m3
C.NO2 = expression(C.NO2.ini), # gN/m3
C.NO3 = expression(C.NO3.ini), # gN/m3
C.O2 = expression(C.O2.ini), # gN/m3
C.DOM = expression(C.DOM.ini)), # gDOM/m3
conc.perarea.ini = list(D.ALG = expression(D.ALG.ini), # gALG/m2
D.HET = expression(D.HET.ini), # gHET/m2
D.N1 = expression(D.N1.ini), # gN1/m2
D.N2 = expression(D.N2.ini), # gN2/m2
D.POM = expression(D.POM.ini)), # gPOM/m2
input = list(C.O2 = expression(K2.O2*L*w*h
*(C.O2.sat-C.O2))), # gas ex.
outflow = expression(Q.in*86400),
processes = list(gro.ALG.NH4,gro.ALG.NO3,resp.ALG,death.ALG,
gro.HET.NH4,gro.HET.NO3,resp.HET,death.HET,
gro.N1,resp.N1,death.N1,gro.N2,resp.N2,death.N2,
hyd.POM))We then define advective links to combine the three reactors:
# Definition of links:
reach1.reach2 <-
new(Class = "link",
name = "R1R2",
from = "R1",
to = "R2",
flow = expression(Q.in*86400))
reach2.reach3 <-
new(Class = "link",
name = "R2R3",
from = "R2",
to = "R3",
flow = expression(Q.in*86400))The definition of the model is concluded with the system definition that contains the reactors, links, environmental conditions, parameters and output time points for simulations:
# Definition of system:
# River system:
system.11.6 <- new(Class = "system",
name = "River",
reactors = list(reach1,reach2,reach3),
links = list(reach1.reach2,reach2.reach3),
cond = cond,
param = param,
t.out = seq(0,3,by=0.02)) #t.out=seq(0,10,by=0.02),Task 3: Model results
Peform a simulation of the model and try to understand the time courses of the state variables and the overall mass fluxes of phosphorus and nitrogen compounds.
The code below starts the simulations.
# Perform simulation:
res.11.6 <- calcres(system.11.6)We then plot the results to the screen and write them also with
different options to files. Please note that you have to adjust the size
of the pdf files with the width and height
arguments to have reasonable sizes of the legend and the descriptions of
the axes. You may also want to plot multiple lines into a single plot by
specifying the argument colnames of the function
plotres (you may type ?plotres to get a
description of the arguments of this function).
# Plot results:
#plotres(res.11.6) # default plots
# plot to screen:
plotres(res.11.6,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")))
# plot to file:
plotres(res.11.6,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")),
file = "exercise_5_results.pdf",
width = 12,
height = 9)## png
## 2
Task 4: Homework: Model sensitivity
What happens in the model if a herbicide enters the water and inhibits the growth of algae? What happens if a bactericide enters the water and inhibits the growth of a specific bacteria group (heterotrophic bacteria, step one nitrifiers, N1, or step two nitrifiers, N2)? Answer these questions by modifying the related kinetic parameters and interpret the changes in the simulation results.
Sensitivity analysis by manually setting parameters to zero:
# Sensitivity analysis
# manually setting parameters to zero:
system.11.6@param$k.gro.ALG <- 0
res.11.6_k.gro.ALG <- calcres(system.11.6)
system.11.6@param <- param
plotres(res.11.6_k.gro.ALG,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")))
plotres(res.11.6_k.gro.ALG,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")),
file = "exercise_5_results_kgroALG.pdf",
width = 12,
height = 9)## png
## 2system.11.6@param$k.gro.HET <- 0
res.11.6_k.gro.HET <- calcres(system.11.6)
system.11.6@param <- param
plotres(res.11.6_k.gro.HET,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")))
plotres(res.11.6_k.gro.HET,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")),
file = "exercise_5_results_kgroHET.pdf",
width = 12,
height = 9)## png
## 2system.11.6@param$k.gro.N1 <- 0
res.11.6_k.gro.N1 <- calcres(system.11.6)
system.11.6@param <- param
plotres(res.11.6_k.gro.N1,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")))
plotres(res.11.6_k.gro.N1,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")),
file = "exercise_5_results_kgroN1.pdf",
width = 12,
height = 9)## png
## 2system.11.6@param$k.gro.N2 <- 0
res.11.6_k.gro.N2 <- calcres(system.11.6)
system.11.6@param <- param
plotres(res.11.6_k.gro.N2,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")))
plotres(res.11.6_k.gro.N2,colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")),
file = "exercise_5_results_kgroN2.pdf",
width = 12,
height = 9)## png
## 2
Sensitivity analysis done by using the sensitivity analysis function
calcsens (note that the scaling factor is now set to 1 and
0 to get the base simulation and the simulations with the parameters set
to zero one after the other):
# use of sensitivity analysis function:
res.96.sens <- calcsens(system = system.11.6,
param.sens = c("k.gro.ALG",
"k.gro.HET",
"k.gro.N1",
"k.gro.N2"),
scaling.factors = c(1,0))
plotres(res = res.96.sens,
colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")))
plotres(res = res.96.sens,
colnames=list(c("C.NH4.R1", "C.NH4.R2", "C.NH4.R3"),
c("C.NO2.R1", "C.NO2.R2", "C.NO2.R3"),
c("C.NO3.R1", "C.NO3.R2", "C.NO3.R3"),
c("C.HPO4.R1","C.HPO4.R2","C.HPO4.R3"),
c("C.O2.R1", "C.O2.R2", "C.O2.R3"),
c("C.DOM.R1", "C.DOM.R2", "C.DOM.R3"),
c("D.ALG.R1", "D.ALG.R2", "D.ALG.R3"),
c("D.HET.R1", "D.HET.R2", "D.HET.R3"),
c("D.N1.R1", "D.N1.R2", "D.N1.R3"),
c("D.N2.R1", "D.N2.R2", "D.N2.R3"),
c("D.POM.R1", "D.POM.R2", "D.POM.R3")),
file = "exercise_5_results_sens.pdf",
width = 16,
height = 12)## png
## 2
Analysis of relevant time scales:
# time scale analysis:
print(paste("hydraulic retention time: ",round(param$L*param$w*param$h/param$Q.in/86400,3),"d"))
print(paste("inverse maximum specific growth rate ALG: ",round(1/param$k.gro.ALG,2),"d"))
print(paste("inverse maximum specific growth rate HET: ",round(1/param$k.gro.HET,2),"d"))
print(paste("inverse maximum specific growth rate N1: ",round(1/param$k.gro.N1,2),"d"))
print(paste("inverse maximum specific growth rate N2: ",round(1/param$k.gro.N2,2),"d"))
print(paste("inverse maximum specific respiration rate ALG:",round(1/param$k.resp.ALG,2),"d"))
print(paste("inverse maximum specific respiration rate HET:",round(1/param$k.resp.HET,2),"d"))
print(paste("inverse maximum specific respiration rate N1: ",round(1/param$k.resp.N1,2),"d"))
print(paste("inverse maximum specific respiration rate N2: ",round(1/param$k.resp.N2,2),"d"))
print(paste("inverse maximum specific death rate ALG: ",round(1/param$k.death.ALG,2),"d"))
print(paste("inverse maximum specific death rate HET: ",round(1/param$k.death.HET,2),"d"))
print(paste("inverse maximum specific death rate N1: ",round(1/param$k.death.N1,2),"d"))
print(paste("inverse maximum specific death rate N2: ",round(1/param$k.death.N2,2),"d"))
print(paste("inverse maximum specific hydrolysis rate: ",round(1/param$k.hyd.POM,2),"d"))## [1] "hydraulic retention time: 0.058 d"
## [1] "inverse maximum specific growth rate ALG: 0.67 d"
## [1] "inverse maximum specific growth rate HET: 0.67 d"
## [1] "inverse maximum specific growth rate N1: 1.25 d"
## [1] "inverse maximum specific growth rate N2: 0.91 d"
## [1] "inverse maximum specific respiration rate ALG: 10 d"
## [1] "inverse maximum specific respiration rate HET: 5 d"
## [1] "inverse maximum specific respiration rate N1: 10 d"
## [1] "inverse maximum specific respiration rate N2: 10 d"
## [1] "inverse maximum specific death rate ALG: 10 d"
## [1] "inverse maximum specific death rate HET: 5 d"
## [1] "inverse maximum specific death rate N1: 10 d"
## [1] "inverse maximum specific death rate N2: 10 d"
## [1] "inverse maximum specific hydrolysis rate: 2 d"It is interesting to note how the hydraulic retention time is one order of magnitude smaller than all the other rates lading to an advection dominated dynamics rather than a process dominated one.
Theory questions
What is the spatial arrangement of the reactors in the river model compared to the previous lake model? How is the advective flow implemented?
What is the difference between planktonic primary production (modelled in the previous lake models) and benthic primary production (modelled in this river model)?
If you have to develop a river water quality model, how do you decide if you model nitrification as a one step process or a two step process (as in section 8.6) or if you model the nitrifying bacteria explicitly (as in section 8.8.2)?
What do you have to consider when deciding about the number and length of the boxes to discretize a river in a model like this one? Hint: Have a look at transport processes described in section 6.1.2.4.
What is the relative importance of advective transport and transformation processes for concentrations in the reactors?