Study the implementation of the model described in section 11.2 for constant driving forces. Run the model implementation below and interpret the behavior of the solutions.
# Model with constant driving forces
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# load required packages:
if ( !require("deSolve") ) {install.packages("deSolve"); library("deSolve") }
if ( !require("ecosim") ) {install.packages("ecosim"); library("ecosim") }
# definition of model parameters:
param <- list(k.gro.ALG = 0.5, # 1/d
k.gro.ZOO = 0.4, # m3/gDM/d
k.death.ALG = 0.1, # 1/d
k.death.ZOO = 0.05, # 1/d
K.HPO4 = 0.002, # gP/m3
Y.ZOO = 0.2, # gDM/gDM
alpha.P.ALG = 0.003, # gP/gDM
A = 5e+006, # m2
h.epi = 5, # m
Q.in = 5, # m3/s
C.HPO4.in = 0.04, # gP/m3
C.HPO4.ini = 0.04, # gP/m3
C.ALG.ini = 0.1, # gDM/m3
C.ZOO.ini = 0.1) # gDM/m3# definition of transformation processes
# growth of algae:
gro.ALG <- new(Class = "process",
name = "Growth of algae",
rate = expression(k.gro.ALG
*C.HPO4/(K.HPO4+C.HPO4)
*C.ALG),
stoich = list(C.ALG = expression(1), # gDM/gDM
C.HPO4 = expression(-alpha.P.ALG))) # gP/gDM
# death of algae: TO BE COMPLETED
death.ALG <- new(Class = ...,
name = ...,
rate = ...,
stoich = ...) # gDM/gDM
# growth of zooplankton:
gro.ZOO <- new(Class = "process",
name = "Growth of zooplankton",
rate = expression(k.gro.ZOO
*C.ALG
*C.ZOO),
stoich = list(C.ZOO = expression(1), # gDM/gDM
C.ALG = expression(-1/Y.ZOO))) # gP/gDM
# death of zooplankton: TO BE COMPLETED
death.ZOO <- new(Class = ...,
name = ...,
rate = ...,
stoich = ...) # gDM/gDM# definition of reactor to describe the epilimnion of the lake: TO BE COMPLETED
epilimnion <-
new(Class = ...,
name = ...,
volume.ini = ...,
conc.pervol.ini = list(C.HPO4 = ..., # gP/m3
C.ALG = ..., # gDM/m3
C.ZOO = ...), # gDM/m3
inflow = ..., # m3/d
inflow.conc = list(C.HPO4 = expression(C.HPO4.in),
C.ALG = 0,
C.ZOO = 0),
outflow = ...,
processes = list(...))
# definition of the system consisting of a single reactor:
# observe the time frame of the simulation
system.11.2.a <- new(Class = "system",
name = "Lake",
reactors = list(epilimnion),
param = param,
t.out = seq(0,2*365,by=1)) # perform simulation:
res.11.2.a <- calcres(system.11.2.a)
# plot results with default options:
plotres(res.11.2.a)
# variables in a vector 'c()' are plotted in the same graph
plotres(res=res.11.2.a,colnames=c("C.ALG", "C.HPO4","C.ZOO"))
# variables in a list 'list()' are plotted in different graphs
plotres(res=res.11.2.a,colnames=list("C.ALG", "C.HPO4", "C.ZOO"))
# combination of the two
plotres(res=res.11.2.a,colnames=list("C.HPO4",c("C.ALG", "C.ZOO")))
# plot and save as pdf
plotres(res = res.11.2.a,
colnames = list("C.HPO4",c("C.ALG", "C.ZOO")),
file = "exercise_2_results_a.pdf",
width = 10,
height = 5)Study the implementation of the model extension to seasonally varying driving forces. First, run the model implementation below.
# Model with seasonally varying driving forces
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# copy the previous system definition:
system.11.2.b <- system.11.2.a
# extend model parameters:
param <- c(param,
list(beta.ALG = 0.046, # 1/degC
beta.ZOO = 0.08, # 1/degC
T0 = 20, # degC
K.I = 30, # W/m2
lambda.1 = 0.10, # 1/m
lambda.2 = 0.10, # m2/gDM
t.max = 230, # d
I0.min = 25, # W/m2
I0.max = 225, # W/m2
T.min = 5, # degC
T.max = 25)) # degC
# extend growth of algae and zooplankton by considering environmental factors:
gro.ALG.ext <-
new(Class = "process",
name = "Growth of algae extended",
rate = expression(k.gro.ALG
*exp(beta.ALG*(T-T0))
*C.HPO4/(K.HPO4+C.HPO4)
*log((K.I+I0)
/(K.I+I0*exp(-(lambda.1+lambda.2*C.ALG)*h.epi)))
/((lambda.1+lambda.2*C.ALG)*h.epi)
*C.ALG),
stoich = list(C.ALG = 1, # gDM/gDM
C.HPO4 = expression(-alpha.P.ALG))) # gP/gDM
gro.ZOO.ext <-
new(Class = "process",
name = "Growth of zooplankton",
rate = expression(k.gro.ZOO
*exp(beta.ZOO*(T-T0))
*C.ALG
*C.ZOO),
stoich = list(C.ZOO = expression(1), # gDM/gDM
C.ALG = expression(-1/Y.ZOO))) # gP/gDM
# re-define processes in the reactor "epilimnion":
epilimnion@processes <- list(gro.ALG.ext,death.ALG,gro.ZOO.ext,death.ZOO)
# make environmental conditions (light and temperature) time dependent:
epilimnion@cond <- list(I0 = expression(0.5*(I0.min+I0.max)+
0.5*(I0.max-I0.min)*
cos(2*pi/365.25*(t-t.max))), # W/m2
T = expression(0.5*(T.min+T.max)+
0.5*(T.max-T.min)*
cos(2*pi/365.25*(t-t.max)))) # degC
# plot the environmental conditions
t <- seq(1,2*365) # for two years
I0 <- numeric(0)
T <- numeric(0)
for(i in 1:length(t))
{
I0[i] <- eval(epilimnion@cond$I0, envir=c(param, t=t[i]))
T[i] <- eval(epilimnion@cond$T, envir=c(param, t=t[i]))
}
par(mfrow=c(1,2),xaxs="i",yaxs="i",mar=c(4.5,4.5,2,1.5)+0.1)
plot(t, I0, type="l")
plot(t, T, type="l")
# re-define the reactor "epilimnion" in the system definition:
system.11.2.b@reactors <- list(epilimnion)
# increase algal growth rate to compensate for new limitations:
param$k.gro.ALG <- 0.8
# replace parameters in the system definition:
system.11.2.b@param <- paramSecond, redo the simulations and plot and interpret the behavior of the solutions.
# redo simulations and plot results:
res.11.2.b <- calcres(system.11.2.b)
plotres(res=res.11.2.b, colnames=list("C.HPO4",c("C.ALG", "C.ZOO")))
plotres(res=res.11.2.b[1:365,], colnames=list("C.HPO4",c("C.ALG", "C.ZOO")))
plotres(res = res.11.2.b, # plot to pdf file
colnames = list("C.HPO4",c("C.ALG","C.ZOO")),
file = "exercise_2_results_b.pdf",
width = 10,
height = 5)We compare the two simulations.
# comparison of the two simulations:
plotres(res = list(const=res.11.2.a,dyn=res.11.2.b),
colnames = list("C.HPO4","C.ALG","C.ZOO"))
plotres(res = list(const=res.11.2.a,dyn=res.11.2.b),
colnames = list("C.HPO4","C.ALG","C.ZOO"),
file = "exercise_2_results_b.pdf",
width = 10,
height = 8)Fill in the missing terms and perform a sensitivity analysis of the model results to the parameters
\(C_{in, HPO^{2-}_{4}}\), \(Q_{in}\), \(k_{gro,ALG}\), \(k_{death,ALG}\), \(k_{gro,ZOO}\), \(k_{death,ZOO}\), and \(Y_{ZOO}\)
by modifying their values by factors of 2 and 1/2 for the model under constant driving forces using the function calcsens. Interpret the results. Hint: look at the solution of the Task 5 of Exercise 1.
# Sensitivity analysis with constant driving forces
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# calculate results of sensitivity analysis: TO BE COMPLETED
res.11.2.a.sens <- calcsens(...,
param.sens=...)
# plot results:
plotres(res = res.11.2.a.sens,
colnames = list("C.HPO4","C.ALG", "C.ZOO"))
plotres(res = res.11.2.a.sens,
colnames = list("C.HPO4","C.ALG", "C.ZOO"),
file = "exercise_2_results_a_sens.pdf",
width = 10,
height = 8)Do a sensitivity analysis for seasonally varying driving forces and discuss the differences to the case with constant driving forces.
# TO BE COMPLETEDAre the algae concentrations controlled bottom-up (by phosphate limitation) or top-down (by grazing of zooplankton)?
What is the reason for oscillating concentrations under constant driving forces? What happens when you introduce periodic driving forces?
What are the main deficits of the model compared to a real lake?
What is your expectation regarding the response of the model to the change in each parameter, does the result match your expectation and can you explain the observed changes?