Classes and methods for implementing aquatic ecosystem models,
for running these models, and for visualizing their results.

Models are built by constructing objects of the classes

`process-class`

,

`reactor-class`

,

`link-class`

,

`system-class`

.

A transformation processes (`process-class`

) is defined by a process rate (expression describing the dependence of the rate on substance or organism concentrations and external influence factors) and stoichiometric coefficients that describe how the rate affects different substances or organisms.
It is recommended to calculate the stoichiometric coefficients with the function `calc.stoich.coef`

of the package `stoichcalc`

from substance and organism compositions.
The output of this function can directly be used for the process definition.
A reactor (`reactor-class`

) describes a well-mixed compartment of the environment (or of a laboratory system).
For each reactor, inflow, outflow, substance and organism input and transformation processes can be defined.
A link (`link-class`

) describes advective and/or diffusive exchange of substances and/or organisms between well-mixed reactors.
Finally, a system (`system-class`

) consists of a single reactor or a set of isolated or linked reactors and can be used to describe a community or meta-community model of an ecosystem and the biogeochemical cycles.

Once a model is described by an object of the class system (`system-class`

), simulations can be performed using the member function

`calcres`

,

This function integrates the system of ordinary differential equations numerically using the function `ode`

of the package `deSolve`

and produces time series of the volumes and substance and organisms concentrations as a R matrix.
The results can be visualized with arbitrary R functions or a summary of all results can be produced with the function

`plotres`

.

To propagate stochasticity and uncertainty to the results, stochastic parameter time series can be generated with the function

`randou`

and parameter samples can be sampled with the function

`randnorm`

to get a sample from the predictive distribution by Monte Carlo simulation.

Package: | ecosim |

Type: | Package |

Version: | 1.2 |

Date: | 2014-02-11 |

License: | GPL (>= 2) |

Depends: | deSolve, stoichcalc |

The following demos are available:

lakemodel_simple

lakemodel_intermediate

lakemodel_complex

rivermodel_simple

rivermodel_complex

Peter Reichert

Maintainer: Peter Reichert <peter.reichert@eawag.ch>

Omlin, M., Reichert, P. and Forster, R., Biogeochemical model of lake Zurich: Model equations and results, Ecological Modelling 141(1-3), 77-103, 2001.

Reichert, P., Borchardt, D., Henze, M., Rauch, W., Shanahan, P., Somlyody, L. and Vanrolleghem, P., River Water Quality Model no. 1 (RWQM1): II. Biochemical process equations, Water Sci. Tech. 43(5), 11-30, 2001.

Reichert, P. and Schuwirth, N., A generic framework for deriving process stoichiometry in environmental models, Environmental Modelling & Software, 25, 1241-1251, 2010.

Soetaert, K., Petzoldt, T., and Woodrow Setzer, R. Solving differential equations in R: Package deSolve. Journal of Statistical Software, 33(9), 2010.

Soetaert, K., Cash, J., and Mazzia, F. Solving Differential Equations in R. Springer, Heidelberg, Germany. 2012.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 | ```
# Definition of parameters:
# =========================
param <- list(k.gro.ALG = 1, # 1/d
k.gro.ZOO = 0.8, # m3/gDM/d
k.death.ALG = 0.4, # 1/d
k.death.ZOO = 0.08, # 1/d
K.HPO4 = 0.002, # gP/m3
Y.ZOO = 0.2, # gDM/gDM
alpha.P.ALG = 0.002, # gP/gDM
A = 8.5e+006, # m2
h.epi = 4, # m
Q.in = 4, # m3/s
C.ALG.ini = 0.05, # gDM/m3
C.ZOO.ini = 0.1, # gDM/m3
C.HPO4.ini = 0.02, # gP/m3
C.HPO4.in = 0.04) # gP/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:
# ---------------
death.ALG <- new(Class = "process",
name = "Death of algae",
rate = expression(k.death.ALG*C.ALG),
stoich = list(C.ALG = expression(-1))) # 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:
# ---------------------
death.ZOO <- new(Class = "process",
name = "Death of zooplankton",
rate = expression(k.death.ZOO*C.ZOO),
stoich = list(C.ZOO = expression(-1))) # gDM/gDM
# Definition of reactor:
# ======================
# Epilimnion:
# -----------
epilimnion <-
new(Class = "reactor",
name = "Epilimnion",
volume.ini = expression(A*h.epi),
conc.pervol.ini = list(C.HPO4 = expression(C.HPO4.ini), # gP/m3
C.ALG = expression(C.ALG.ini), # gDM/m3
C.ZOO = expression(C.ZOO.ini)), # gDM/m3
inflow = expression(Q.in*86400), # m3/d
inflow.conc = list(C.HPO4 = expression(C.HPO4.in),
C.ALG = 0,
C.ZOO = 0),
outflow = expression(Q.in*86400),
processes = list(gro.ALG,death.ALG,gro.ZOO,death.ZOO))
# Definition of system:
# =====================
# Lake system:
# ------------
system <- new(Class = "system",
name = "Lake",
reactors = list(epilimnion),
param = param,
t.out = seq(0,365,by=1))
# Perform simulation:
# ===================
res <- calcres(system)
# Plot results:
# =============
plotres(res) # plot to screen
plotres(res,file="ecosim_example_plot1.pdf") # plot to pdf file
plotres(res, colnames=c("C.ALG", "C.ZOO")) # plot selected variables
plotres(res, colnames=list("C.HPO4",c("C.ALG", "C.ZOO")))
plotres(res[1:100,], colnames=list("C.HPO4",c("C.ALG", "C.ZOO"))) # plot selected time steps
plotres(res = res, # plot to pdf file
colnames = list("C.HPO4",c("C.ALG","C.ZOO")),
file = "ecosim_example_plot2.pdf",
width = 8,
height = 4)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.