---
title: 'Exercise 4: Biogeochemical-Ecological Lake Model'
author: "<Your Name>"
date: "April 17, 2024"
output:
  pdf_document: default
  word_document: default
  html_document: default
subtitle: 'ETH Zurich Course 701-0426-00L: Modelling Aquatic Ecosystems (Schuwirth)'
---

```{r setup, include=TRUE}
# to display chunk code along with its results
knitr::opts_chunk$set(echo = TRUE)
```


### Goals:

* Understand the concepts and process formulations of the coupled biogeochemical-ecological lake model described in section 11.4.

* Understand the implementation of this model based on the R packages `stoichcalc` and `ecosim`.

* 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.4. 
Note that it may be useful to study first the ''intermediately complex'' model described in section 11.3.

### 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).

```{r}

# load required packages:

if ( !require("ecosim") ) {install.packages("ecosim"); library(ecosim) }
```

Note that the packages `deSolve`, to numerially solve systems of ordinary differential equations, 
and `stoichcalc`, to calculate process stoichiometry were loaded automatically.

#### 2.1 Define parameters and calculate stoichiometric coefficients

Define parameter values for stoichiometry, kinetics, reactor configuration and input:

```{r}

# Define parameter values:

param <- list(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.005,        # gP/gALG
              alpha.O.ZOO          = 0.50,         # gO/gZOO
              alpha.H.ZOO          = 0.07,         # gH/gZOO
              alpha.N.ZOO          = 0.06,         # gN/gZOO
              alpha.P.ZOO          = 0.01,         # gP/gZOO
              alpha.O.POM          = 0.39,         # gO/gPOM
              alpha.H.POM          = 0.07,         # gH/gPOM
              alpha.N.POM          = 0.06,         # gN/gPOM
              alpha.P.POM          = 0.007,        # gP/gPOM
              Y.ZOO                = 0.2,          # gZOO/gALG
              f.e                  = 0.2,          # gPOM/gALG
              f.I                  = 0.2,          # gPOMI/(gPOMD+gPOMI)
              k.gro.ALG            = 0.8,          # 1/d
              k.gro.ZOO            = 0.4,          # m3/gDM/d
              k.resp.ALG           = 0.10,         # 1/d
              k.resp.ZOO           = 0.10,         # 1/d
              k.death.ALG          = 0.10,         # 1/d
              k.death.ZOO          = 0.05,         # 1/d
              k.nitri              = 0.1,          # gN/m3/d
              k.miner.ox.POM       = 0.02,         # 1/d
              k.miner.ox.POM.sed   = 5.0,          # gDM/m2/d
              k.miner.anox.POM.sed = 5.0,          # gDM/m2/d
              K.POM.miner.sed      = 10,           # gDM/m2
              K.HPO4               = 0.002,        # gP/m3
              K.N                  = 0.04,         # gN/m3
              p.NH4                = 5,            # - 
              K.O2.ZOO             = 0.2,          # gO/m3
              K.O2.resp            = 0.5,          # gO/m3
              K.O2.nitri           = 0.4,          # gO/m3
              K.O2.miner           = 0.5,          # gO/m3
              K.NO3.miner          = 0.1,          # gN/m3
              K.NH4.nitri          = 0.5,          # gN/m3
              A                    = 5e+006,       # m2
              h.epi                = 5,            # m
              h.hypo               = 10,           # m
              Q.in                 = 5,            # m3/s
              C.HPO4.in            = 0.04,         # gP/m3
              C.NO3.in             = 0.5,          # gN/m3
              C.O2.in              = 10,           # gO/m3
              C.ALG.ini            = 0.1,          # gDM/m3
              C.ZOO.ini            = 0.1,          # gDM/m3
              C.POMD.ini           = 0.0,          # gDM/m3
              C.POMI.ini           = 0.0,          # gDM/m3
              C.HPO4.ini           = 0.04,         # gP/m3
              C.NH4.ini            = 0.1,          # gN/m3
              C.NO3.ini            = 0.5,          # gN/m3
              C.O2.ini             = 10,           # gO/m3
              D.POMD.ini           = 0,            # gDM/m2
              D.POMI.ini           = 0,            # gDM/m2
              beta.ALG             = 0.046,        # 1/degC
              beta.ZOO             = 0.08,         # 1/degC
              beta.BAC             = 0.046,        # 1/degC
              T0                   = 20,           # degC
              K.I                  = 30,           # W/m2
              lambda.1             = 0.10,         # 1/m
              lambda.2             = 0.10,         # m2/gDM  
              v.ex.O2              = 1,            # m/d
              v.sed.POM            = 1,            # m/d
              Kz.summer            = 0.02,         # m2/d
              Kz.winter            = 20,           # m2/d
              h.meta               = 5,            # m
              t.max                = 230,          # d
              I0.min               = 25,           # W/m2
              I0.max               = 225,          # W/m2
              T.min                = 5,            # degC
              T.max                = 25,           # degC
              p                    = 101325)       # Pa

```

This set of parameters contains the elemental mass fractions of nitrogen, phosphorus, oxygen and hydrogen.
The largest mass fraction, of carbon, is then calculated as the difference to unity:

```{r}
# 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.ZOO <- 1 - (param$alpha.O.ZOO + param$alpha.H.ZOO + 
                        param$alpha.N.ZOO + param$alpha.P.ZOO)
param$alpha.C.POM <- 1 - (param$alpha.O.POM + param$alpha.H.POM + 
                        param$alpha.N.POM + param$alpha.P.POM)

```

Next we define the composition of all substances to be considered for the calculation of the stoichiometries 
of all processes considered in the model.
Note that we have to be careful with the units.
The natural unit to measure the amount of inorganic chemical compounds is moles.
This is not the case for organic compounds as they consist of a large number of different organic molecules 
and we just consider their mean composition. 
For this reason, we measure organic compounds as dry mass.
As we have to consider conservation of charge as well as conservation of the elements, we need to add charge to the composition vectors:

```{r}
# Construct composition matrix:

NH4    <- c(H      = 4*1/14,                # gH/gNH4-N
            N      = 1,                     # gN/gNH4-N
            charge = 1/14)                  # chargeunits/gNH4-N
NO3    <- c(O      = 3*16/14,               # gO/gNO3-N
            N      = 1,                     # gN/gNO3-N
            charge = -1/14)                 # chargeunits/gNO3-N
N2     <- c(N      = 1)                     # gN/gN2-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
ZOO    <- c(C      = param$alpha.C.ZOO,     # gC/gZOO
            O      = param$alpha.O.ZOO,     # gO/gZOO
            H      = param$alpha.H.ZOO,     # gH/gZOO
            N      = param$alpha.N.ZOO,     # gN/gZOO
            P      = param$alpha.P.ZOO)     # gP/gZOO
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
```

Once 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`:

```{r}
subst.comp <- list(C.NH4    = NH4,
                   C.NO3    = NO3,
                   C.N2     = N2,
                   C.HPO4   = HPO4,
                   C.HCO3   = HCO3,
                   C.O2     = O2,
                   C.H      = H,
                   C.H2O    = H2O,
                   C.ALG    = ALG,
                   C.ZOO    = ZOO,
                   C.POMD   = POM,
                   D.POMD   = POM,
                   C.POMI   = POM,
                   D.POMI   = POM)

alpha <- calc.comp.matrix(subst.comp)
print(round(alpha,3))
```

Death converts living algae or zooplankton into dead particulate organic matter.
The typical stoichiometry of this process would be -1 for algae or zooplankton and +1 for particulate organic matter.
However, we run into troubles if the composition of some of these organic components are different.
This could be resolved by introducing dead organic particles with all different compositions of living organisms.
To avoid this increase in the number of state variables (in particular if there would be more organisms to distinguish) 
and as most of the dead organic particles will be mineralized later anyway, we introduced a ``yield'' for death so that 
part of the organism becomes dead particulate organic matter and the remaining part is mineralized.
The idea is to keep the fraction that becomes particulate organic matter as large as possible and 
not to have negative stoichiometric coefficients of nutrients or oxygen.
The criteria we apply here are that we do not need an uptake of nitrogen or phosphorus for dying.
Oxygen will have to be checked in the resulting stoichiometry.
This leads to the following expressions for the yield:

```{r}
# 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.ZOO.death <-  min(1,
                        param$alpha.N.ZOO/param$alpha.N.POM,
                        param$alpha.P.ZOO/param$alpha.P.POM,
                        param$alpha.C.ZOO/param$alpha.C.POM)

print(param$Y.ALG.death)
print(param$Y.ZOO.death)
```

Having defined the composition matrix and the stoichiometric parameters, we can define the stoichiometries of the processes to be considered by the model.
Each process requires the composition matrix, it has a name, a list of substances to be considered, an indication of which coefficient should be chosen to make the stoichiometry unique, the value to which this coefficient is set (typically +1 or -1), and, optionally, the coefficients of additional constraints.
These coefficients have to be provided in the form $\sum_j \nu_j \gamma_j = 0$.
This means that e.g. the condition $\nu_\mathrm{ZOO}=- Y \nu_{ALG}$ indicating that one mass unit of algae is only converted into $Y$ (= yield) mass units of zooplankton has to be transformed into $\nu_\mathrm{ZOO} + Y \nu_{ALG} = 0$ to read the coefficients (see the first constraint below for the growth of zooplankton).

```{r}

# Calculate stoichiometric coefficients of selected processes:

# 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","C.ALG"),
                    subst.norm  = "C.ALG",
                    nu.norm     = 1)

# 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","C.ALG"),
                    subst.norm  = "C.ALG",
                    nu.norm     = 1)

# 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","C.ALG"),
                    subst.norm  = "C.ALG",
                    nu.norm     = -1)

# 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","C.ALG","C.POMD","C.POMI"),
                    subst.norm  = "C.ALG",
                    nu.norm     = -1,
                    constraints = list(c("C.ALG"  = param$Y.ALG.death,
                                         "C.POMD" = 1,
                                         "C.POMI" = 1),
                                       c("C.POMD" = -param$f.I,
                                         "C.POMI" = 1-param$f.I)))

# Growth of zooplankton:

nu.gro.ZOO <- 
   calc.stoich.coef(alpha       = alpha,
                    name        = "gro.ZOO",
                    subst       = c("C.NH4","C.HPO4","C.HCO3","C.O2","C.H",
                                    "C.H2O","C.ALG","C.ZOO","C.POMD","C.POMI"),
                    subst.norm  = "C.ZOO",
                    nu.norm     = 1,
                    constraints = list(c("C.ZOO"  = 1,
                                         "C.ALG"  = param$Y.ZOO),
                                       c("C.POMD" = 1,
                                         "C.POMI" = 1,
                                         "C.ALG"  = param$f.e),
                                       c("C.POMD" = -param$f.I,
                                         "C.POMI" = 1-param$f.I)))

# Respiration of zooplankton:

nu.resp.ZOO <- 
   calc.stoich.coef(alpha       = alpha,
                    name        = "resp.ZOO",
                    subst       = c("C.NH4","C.HPO4","C.HCO3","C.O2",
                                    "C.H","C.H2O","C.ZOO"),
                    subst.norm  = "C.ZOO",
                    nu.norm     = -1)

# Death of zooplankton:

nu.death.ZOO <- 
   calc.stoich.coef(alpha       = alpha,
                    name        = "death.ZOO",
                    subst       = c("C.NH4","C.HPO4","C.HCO3","C.O2",
                                    "C.H","C.H2O","C.ZOO","C.POMD","C.POMI"),
                    subst.norm  = "C.ZOO",
                    nu.norm     = -1,
                    constraints = list(c("C.ZOO"  = param$Y.ZOO.death,
                                         "C.POMD" = 1,
                                         "C.POMI" = 1),
                                       c("C.POMD" = -param$f.I,
                                         "C.POMI" = 1-param$f.I)))

# Nitrification:

nu.nitri <- 
   calc.stoich.coef(alpha       = alpha,
                    name        = "nitri",
                    subst       = c("C.NH4","C.NO3","C.O2","C.H","C.H2O"),
                    subst.norm  = "C.NH4",
                    nu.norm     = -1)

# Oxic mineralization of suspended organic particles:

nu.miner.ox.POM <- 
   calc.stoich.coef(alpha       = alpha,
                    name        = "miner.ox.POM",
                    subst       = c("C.NH4","C.HPO4","C.HCO3","C.O2",
                                    "C.H","C.H2O","C.POMD"),
                    subst.norm  = "C.POMD",
                    nu.norm     = -1)

# Oxic mineralization of sedimented organic particles:

nu.miner.ox.POM.sed <- 
   calc.stoich.coef(alpha       = alpha,
                    name        = "miner.ox.POM.sed",
                    subst       = c("C.NH4","C.HPO4","C.HCO3","C.O2",
                                    "C.H","C.H2O","D.POMD"),
                    subst.norm  = "D.POMD",
                    nu.norm     = -1)

# Anoxic mineralization of sedimented organic particles:

nu.miner.anox.POM.sed <- 
   calc.stoich.coef(alpha       = alpha,
                    name        = "miner.anox.POM.sed",
                    subst       = c("C.NH4","C.NO3","C.HPO4","C.HCO3","C.N2",
                                    "C.H","C.H2O","D.POMD"),
                    subst.norm  = "D.POMD",
                    nu.norm     = -1,
                    constraints = list(c("C.NO3" = 1,
                                         "C.N2"  = 1)))

# Sedimentation of degradable organic particles:

nu.sed.POMD <- 
   calc.stoich.coef(alpha       = alpha,
                    name        = "sed.POMD",
                    subst       = c("C.POMD","D.POMD"),
                    subst.norm  = "C.POMD",
                    nu.norm     = -1)

# Sedimentation of inert organic particles:

nu.sed.POMI <- 
   calc.stoich.coef(alpha       = alpha,
                    name        = "sed.POMI",
                    subst       = c("C.POMI","D.POMI"),
                    subst.norm  = "C.POMI",
                    nu.norm     = -1)
```

Finally, we bind the stoichiometries together for the full stoichiometric matrix:

```{r}

# Construct stoichiometric matrix:

nu <- rbind(nu.gro.ALG.NH4,
            nu.gro.ALG.NO3,
            nu.resp.ALG,
            nu.death.ALG,
            nu.gro.ZOO,
            nu.resp.ZOO,
            nu.death.ZOO,
            nu.nitri,
            nu.miner.ox.POM,
            nu.miner.ox.POM.sed,
            nu.miner.anox.POM.sed,
            nu.sed.POMD,
            nu.sed.POMI)
print(round(nu,3))
```

This stoichiometric matrix should be checked in particular for the signs of the stoichiometric coefficients.
As we chose the yields for death accordingly above, death of neither algae or zooplankton consumes phosphate or ammonium.
However, we also have to check that it does not consume oxygen, which is not so easy to avoid.
If it would consume oxygen, we would have to switch to the approach of using multiple fractions of dead organic matter (we do not have this problem if the compostion of dead organic matter is the same as that of the dying organism).

#### 2.2 Define the transformation processes and the final system

We can now proceed with the definition of transformation processes.
This is analogous as we did it in the exercises 1 and 2 with the exception that we can use the stoichiometries as calculated above with `stoichcalc`.
As we pass the numerical values of the stoichiometric coefficients that were evaluated above using the model parameters, we will not be able to do sensitivity analyses of the model with regard to parameters that affect the stoichiometry by just redoing simulations with modified parameter values (we can still do so for the parameters that characterize kinetics, input, etc.).
If we are interested in sensitivity regarding stoichiometric parameters, we have to re-evaluate stoichiometry with modified parameter values, pass the modified stoichiometry to the processes in the model and then redo the simulations.

```{r}

# 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))
                            *log((K.I+I0)/
                                (K.I+I0*exp(-(lambda.1+lambda.2*C.ALG)*h.epi)))
                                /((lambda.1+lambda.2*C.ALG)*h.epi)
                            *min(C.HPO4/(K.HPO4+C.HPO4),
                                (C.NH4+C.NO3)/(K.N+C.NH4+C.NO3))
                            *(p.NH4*C.NH4/(p.NH4*C.NH4+C.NO3))
                            *C.ALG),
       stoich = as.list(nu["gro.ALG.NH4",]))

# 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))
                           *log((K.I+I0)/
                                (K.I+I0*exp(-(lambda.1+lambda.2*C.ALG)*h.epi)))
                                /((lambda.1+lambda.2*C.ALG)*h.epi)
                           *min(C.HPO4/(K.HPO4+C.HPO4),
                                (C.NH4+C.NO3)/(K.N+C.NH4+C.NO3))
                           *(C.NO3/(p.NH4*C.NH4+C.NO3))
                           *C.ALG),
       stoich = as.list(nu["gro.ALG.NO3",]))
       
# 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.resp+C.O2))*C.ALG),
       stoich = as.list(nu["resp.ALG",]))

# Death of algae:

death.ALG <- 
   new(Class = "process",
       name   = "death.ALG",
       rate   = expression(k.death.ALG*C.ALG),
       stoich = as.list(nu["death.ALG",]))

# Growth of zooplankton:

gro.ZOO <- 
   new(Class  = "process",
       name   = "gro.ZOO",
       rate   = expression(k.gro.ZOO*exp(beta.ZOO*(T-T0))
                           *(C.O2/(K.O2.ZOO+C.O2))
                           *C.ALG
                           *C.ZOO),
       stoich = as.list(nu["gro.ZOO",]))

# Respiration of zooplankton:

resp.ZOO <- 
   new(Class  = "process",
       name   = "resp.ZOO",
       rate   = expression(k.resp.ZOO*exp(beta.ZOO*(T-T0))*
                             (C.O2/(K.O2.resp+C.O2))*C.ZOO),
       stoich = as.list(nu["resp.ZOO",]))

# Death of zooplankton:

death.ZOO <- 
   new(Class = "process",
       name   = "death.ZOO",
       rate   = expression(k.death.ZOO*C.ZOO),
       stoich = as.list(nu["death.ZOO",]))
```

In a second step, we define additional processes: 

* nitrification, 

* oxic mineralization of `POMD` and `SPOMD`, 

* anoxic mineralization of `SPOMD`. 

The missing equations for the expression of the rates are completed based on Table 11.10.

```{r}

# Nitrification:

nitri <- 
   new(Class = "process",
       name   = "nitri",
       rate   = expression(k.nitri
                           *exp(beta.BAC*(T-T0))
                           *min(C.NH4/(K.NH4.nitri+C.NH4),C.O2/
                                  (K.O2.nitri+C.O2))),
       stoich = as.list(nu["nitri",]))

# Oxic mineralization of suspended organic particles:

miner.ox.POM <- 
   new(Class  = "process",
       name   = "miner.ox.POM",
       rate   = expression(k.miner.ox.POM
                           *exp(beta.BAC*(T-T0))
                           *C.O2/(K.O2.miner+C.O2)
                           *C.POMD),
       stoich = as.list(nu["miner.ox.POM",]))

# Oxic mineralization of sedimented organic particles:

miner.ox.POM.sed <- 
   new(Class  = "process",
       name   = "miner.ox.POM.sed",
       rate   = expression(k.miner.ox.POM.sed
                           *exp(beta.BAC*(T-T0))
                           *C.O2/(K.O2.miner+C.O2)
                           *D.POMD/(K.POM.miner.sed+D.POMD)),
       stoich = as.list(nu["miner.ox.POM.sed",]),
       pervol = F)

# Anoxic mineralization of sedimented organic particles:

miner.anox.POM.sed <- 
   new(Class  = "process",
       name   = "miner.anox.POM.sed",
       rate   = expression(k.miner.anox.POM.sed
                           *exp(beta.BAC*(T-T0))
                           *C.NO3/(K.NO3.miner+C.NO3)
                           *(D.POMD/(K.POM.miner.sed+D.POMD))^2),
       stoich = as.list(nu["miner.anox.POM.sed",]),
       pervol = F)
```

The sediment modeled here is simplified. To model a real sediment, we would need different sediment layers, in-situ concentrations varying with the depth of the sediment, and diffusion. We could then model the oxygen diffusion limitation in the deeper sediment and implement the same process rates as in chapter 8.5 that would depend on the in-situ concentrations of oxygen, nitrate etc. in the sediment. With increasing sediment thickness, the mineralization rates would become constant. Here we do not have a layer for the sediment, so we model it with a state variable for the density of sedimented particulate organic matter, $D_{SPOMD}$, and the behavior of oxic and anoxic mineralization are differentiated with a Monod limitation term and a quadratic Monod limitation term, respectively (see Figure 4.2). This ensures that at small sediment thickness (i.e low density of `SPOM`), oxic mineralization is dominating, but after a certain density of `SPOM`, exactly $K_{SPOM,miner,sed}$, anoxic mineralization becomes more important than oxic mineralization.

As a last step, we define the processes of sedimentation of `POMD` and `POMI`.

```{r}
# Sedimentation of degradable organic particles:

sed.POMD <- 
   new(Class  = "process",
       name   = "sed.POMD",
       rate   = expression(v.sed.POM/h.hypo*C.POMD),
       stoich = as.list(nu["sed.POMD",]))

# Sedimentation of inert organic particles:

sed.POMI <- 
   new(Class  = "process",
       name   = "sed.POMI",
       rate   = expression(v.sed.POM/h.hypo*C.POMI),
       stoich = as.list(nu["sed.POMI",]))

```

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.
In addition, we define the turbulent diffusion coefficient, $K_z$, that quantifies mixing between the hypolimnion (deep water) and the epiliminion (upper part of the water column).
To describe, mixing during the winter and stratification in the summer, we use a large value in winter and a small value in summer.

```{r}
# Definition of environmental conditions:

cond.epi  <- list(I0       = expression( 0.5*(I0.min+I0.max)
                                        +0.5*(I0.max-I0.min)
                                            *cos(2*pi/365.25*(t-t.max))),
                  T        = expression( 0.5*(T.min+T.max)
                                        +0.5*(T.max-T.min)
                                            *cos(2*pi/365.25*(t-t.max))),
                  C.O2.sat = expression(exp(7.7117-1.31403*log(T+45.93))
                                        *p/101325))
                                          
cond.hypo <- list(I0 = 0,
                  T  = 5)

cond.gen  <- list(Kz=expression( 0.5*(Kz.summer+Kz.winter)
                                -0.5*(Kz.winter-Kz.summer)
                                    *sign(cos(2*pi/365.25*(t-t.max))+0.4)))
                                    
# Plot the environmental conditions

t <-  1:365 
  
Kz <- numeric(0)
I0 <- numeric(0)
T <- numeric(0)  
C.O2.sat <- numeric(0)      
                       
for(i in 1:length(t))
{
   Kz[i]       <- eval(cond.gen$Kz, envir=c(param, t=t[i]))
   I0[i]       <- eval(cond.epi$I0, envir=c(param, t=t[i]))
   T[i]        <- eval(cond.epi$T, envir=c(param, t=t[i]))
   C.O2.sat[i] <- eval(cond.epi$C.O2.sat, envir=c(param, T=T[i]))
}

#Plots
par.def <- par(no.readonly=TRUE)  # store current plotting parameter
par(mfrow=c(2,2),
    xaxs="i",
    yaxs="i",mar=c(4.5,4.5,2,1.5)+0.1) # modify plot parameters                            
plot(t, Kz , ylim=c(0,25))
plot(t, I0, type="l")
plot(t, T, type="l")
plot(t, C.O2.sat, type="l")
par(par.def)  # restore previous plotting parameters 
```

We now define the two boxes that allow us to describe mixed and stratified conditions of the lake.
Inflow and outflow are into and out of the upper box, called epilimnion.
In the lower box, called hypolimnion, there is no primary production due to the absence of light.
Also the sediment is assumed to be only relevant in the hypoliminion, which is typical for deep lakes with steep banks.

```{r}

# Definition of reactors:

# Epilimnion:

epilimnion <- 
   new(Class            = "reactor",
       name             = "Epi",
       volume.ini       = expression(A*h.epi),
       conc.pervol.ini  = list(C.HPO4 = expression(C.HPO4.ini),  # gP/m3
                               C.NH4  = expression(C.NH4.ini),   # gN/m3
                               C.NO3  = expression(C.NO3.ini),   # gN/m3
                               C.O2   = expression(C.O2.ini),    # gO/m3
                               C.ALG  = expression(C.ALG.ini),   # gDM/m3
                               C.ZOO  = expression(C.ZOO.ini),   # gDM/m3
                               C.POMD = expression(C.POMD.ini),  # gDM/m3
                               C.POMI = expression(C.POMI.ini)), # gDM/m3
       input            = list(C.O2   = expression(v.ex.O2*A
                                                   *(C.O2.sat-C.O2))), # gas ex.
       inflow           = expression(Q.in*86400),                # m3/d
       inflow.conc      = list(C.HPO4 = expression(C.HPO4.in),
                               C.NH4  = 0,
                               C.NO3  = expression(C.NO3.in),
                               C.O2   = expression(C.O2.in),
                               C.ALG  = 0,
                               C.ZOO  = 0,
                               C.POMD = 0,
                               C.POMI = 0),
       outflow          = expression(Q.in*86400),
       cond             = cond.epi,
       processes        = list(gro.ALG.NH4,gro.ALG.NO3,resp.ALG,death.ALG,
                               gro.ZOO,resp.ZOO,death.ZOO,
                               nitri,miner.ox.POM))

# Hypolimnion:

hypolimnion <- 
   new(Class            = "reactor",
       name             = "Hypo",
       volume.ini       = expression(A*h.hypo),
       area             = expression(A),
       conc.pervol.ini  = list(C.HPO4 = expression(C.HPO4.ini),  # gP/m3
                               C.NH4  = expression(C.NH4.ini),   # gN/m3
                               C.NO3  = expression(C.NO3.ini),   # gN/m3
                               C.O2   = expression(C.O2.ini),    # gO/m3
                               C.ALG  = expression(C.ALG.ini),   # gDM/m3
                               C.ZOO  = expression(C.ZOO.ini),   # gDM/m3
                               C.POMD = expression(C.POMD.ini),  # gDM/m3
                               C.POMI = expression(C.POMI.ini)), # gDM/m3
       conc.perarea.ini = list(D.POMD = expression(D.POMD.ini),  # gDM/m2
                               D.POMI = expression(D.POMI.ini)), # gDM/m2
       cond             = cond.hypo,
       processes        = list(resp.ALG,death.ALG,
                               gro.ZOO,resp.ZOO,death.ZOO,
                               nitri,
                               miner.ox.POM,miner.ox.POM.sed,miner.anox.POM.sed,
                               sed.POMD,sed.POMI))

```

The final element needed is the metalimnion that separates the epilimnion from the hypolimnion.
We neglect its volume (this means we expand the epilimnion and the hypolimnion) and specify directed flow of dead organic particles through sedimentation from the epilimnion to the hypolimnion and diffusive exchange of all substances with the turbulent diffusion coefficient $K_z$ that was defined above to be large in winter and small in summer.

```{r}

# Definition of links:

# Exchange between epilimnion and hypolimnion:

metalimnion <-
   new(Class      = "link",
       name       = "Metalimnion",
       from       = "Epi",
       to         = "Hypo",
       qadv.spec  = list(C.POMD = expression(v.sed.POM*A),
                         C.POMI = expression(v.sed.POM*A)),
       qdiff.gen  = expression(A/h.meta*Kz))
```

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:

```{r}

# Definition of the lake system:

system.11.4 <- new(Class    = "system",
                   name     = "Lake",
                   reactors = list(epilimnion,hypolimnion),
                   links    = list(metalimnion),
                   cond     = cond.gen,
                   param    = param,
                   t.out    = seq(0,730,by=1)) #number of days 365, 730, 1095 ,1460

```

#### 2.3 Perform simulations and plot results

Now, we start simulations.
The default option is to use an implicit integration scheme that is stable also for large time steps; as an alternative, we could apply an explicit Runge-Kutta scheme that is faster for each time step but needs more of them:

```{r}

# Perform simulation:

ptm <- proc.time()
res.11.4 <- calcres(system.11.4)
print(proc.time()-ptm)

ptm <- proc.time()
res.11.4.rk4 <- calcres(system.11.4,method="rk4")
print(proc.time()-ptm)
```

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).

```{r fig.height=10,fig.width=10}

# Plot results:

plotres(res      = res.11.4,
        colnames =list(c("C.NH4.Epi","C.NH4.Hypo"),
                       c("C.NO3.Epi","C.NO3.Hypo"),
                       c("C.HPO4.Epi","C.HPO4.Hypo"),
                       c("C.O2.Epi","C.O2.Hypo"),
                       c("C.ALG.Epi","C.ALG.Hypo"),
                       c("C.ZOO.Epi","C.ZOO.Hypo"),
                       c("C.POMD.Epi","C.POMD.Hypo"),
                       c("C.POMI.Epi","C.POMI.Hypo"),
                       c("D.POMD.Hypo","D.POMI.Hypo")))
                  
plotres(res      = res.11.4,
        colnames = list(c("C.NH4.Epi","C.NH4.Hypo"),
                        c("C.NO3.Epi","C.NO3.Hypo"),
                        c("C.HPO4.Epi","C.HPO4.Hypo"),
                        c("C.O2.Epi","C.O2.Hypo"),
                        c("C.ALG.Epi","C.ALG.Hypo"),
                        c("C.ZOO.Epi","C.ZOO.Hypo"),
                        c("C.POMD.Epi","C.POMD.Hypo"),
                        c("C.POMI.Epi","C.POMI.Hypo"),
                        c("D.POMD.Hypo","D.POMI.Hypo")),
        file     = "exercise_4_results.pdf",
        width    = 8,
        height   = 8)

# compare numerical solutions:
plotres(res      = list(lsoda=res.11.4,rk4=res.11.4.rk4),
        colnames = list(c("C.NH4.Epi","C.NH4.Hypo"),
                        c("C.NO3.Epi","C.NO3.Hypo"),
                        c("C.HPO4.Epi","C.HPO4.Hypo"),
                        c("C.O2.Epi","C.O2.Hypo"),
                        c("C.ALG.Epi","C.ALG.Hypo"),
                        c("C.ZOO.Epi","C.ZOO.Hypo"),
                        c("C.POMD.Epi","C.POMD.Hypo"),
                        c("C.POMI.Epi","C.POMI.Hypo"),
                        c("D.POMD.Hypo","D.POMI.Hypo")),
        file     = "exercise_4_results_rk4.pdf",
        width    = 8,
        height   = 8)

```

Overall mass balance of phosphorus over the whole simulation period:

```{r results="hold"}

# Phosphorus mass balance:

nr.days  <- (as.numeric(rownames(res.11.4)[nrow(res.11.4)])-
               as.numeric(rownames(res.11.4)[1]))
nr.steps <- (nrow(res.11.4)-1)

F.in.P  <- c(HPO4    = param$Q.in*param$C.HPO4.in*nr.days*86400/1e6)
F.out.P <- c(HPO4    = sum(param$Q.in*res.11.4[,"C.HPO4.Epi"])*
               nr.days/nr.steps*86400/1e6,
             ALG     = sum(param$Q.in*res.11.4[,"C.ALG.Epi"]*param$alpha.P.ALG)*
                         nr.days/nr.steps*86400/1e6,
             ZOO     = sum(param$Q.in*res.11.4[,"C.ZOO.Epi"]*param$alpha.P.ZOO)*
                         nr.days/nr.steps*86400/1e6,
             POMD    = sum(param$Q.in*res.11.4[,"C.POMD.Epi"]*param$alpha.P.POM)*
                         nr.days/nr.steps*86400/1e6,
             POMI    = sum(param$Q.in*res.11.4[,"C.POMI.Epi"]*param$alpha.P.POM)*
                         nr.days/nr.steps*86400/1e6)
Acc.P   <- c(HPO4    = param$A/1e6*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.HPO4.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.HPO4.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.HPO4.Epi"]+
                            param$h.hypo*res.11.4[1,"C.HPO4.Hypo"])),
             ALG     = param$A/1e6*param$alpha.P.ALG*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.ALG.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.ALG.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.ALG.Epi"]+
                            param$h.hypo*res.11.4[1,"C.ALG.Hypo"])),
             ZOO     = param$A/1e6*param$alpha.P.ZOO*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.ZOO.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.ZOO.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.ZOO.Epi"]+
                            param$h.hypo*res.11.4[1,"C.ZOO.Hypo"])),
             POMD    = param$A/1e6*param$alpha.P.POM*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.POMD.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.POMD.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.POMD.Epi"]+
                            param$h.hypo*res.11.4[1,"C.POMD.Hypo"])),
             POMI    = param$A/1e6*param$alpha.P.POM*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.POMI.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.POMI.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.POMI.Epi"]+
                            param$h.hypo*res.11.4[1,"C.POMI.Hypo"])),
             POMDsed = param$A/1e6*param$alpha.P.POM*
                         (res.11.4[nrow(res.11.4),"D.POMD.Hypo"]-
                            res.11.4[1,"D.POMD.Hypo"]),
             POMIsed = param$A/1e6*param$alpha.P.POM*
                         (res.11.4[nrow(res.11.4),"D.POMI.Hypo"]-
                            res.11.4[1,"D.POMI.Hypo"]))

print("Input (tP)")
print(signif(F.in.P,3))
print("Output (tP)")
print(signif(F.out.P,3))
print(signif(sum(F.out.P),3))
print("Accumulation (tP)")
print(signif(Acc.P,3))
print(signif(sum(Acc.P),3))

print("Input - Output - Accumulation:")
print(sum(F.in.P)-sum(F.out.P)-sum(Acc.P))

```

Overall mass balance of nitrogen over the whole simulation period:

```{r results="hold"}

# Nitrogen mass balance:

nr.days <- (as.numeric(rownames(res.11.4)[nrow(res.11.4)])-
              as.numeric(rownames(res.11.4)[1]))
nr.steps <- (nrow(res.11.4)-1)

F.in.N  <- c(NO3     = param$Q.in*param$C.NO3.in*nr.days*86400/1e6)
F.out.N <- c(NO3     = sum(param$Q.in*res.11.4[,"C.NO3.Epi"])*
               nr.days/nr.steps*86400/1e6,
             NH4     = sum(param$Q.in*res.11.4[,"C.NH4.Epi"])*
               nr.days/nr.steps*86400/1e6,
             ALG     = sum(param$Q.in*res.11.4[,"C.ALG.Epi"]*param$alpha.N.ALG)*
                         nr.days/nr.steps*86400/1e6,
             ZOO     = sum(param$Q.in*res.11.4[,"C.ZOO.Epi"]*param$alpha.N.ZOO)*
                         nr.days/nr.steps*86400/1e6,
             POMD    = sum(param$Q.in*res.11.4[,"C.POMD.Epi"]*param$alpha.N.POM)*
                         nr.days/nr.steps*86400/1e6,
             POMI    = sum(param$Q.in*res.11.4[,"C.POMI.Epi"]*param$alpha.N.POM)*
                         nr.days/nr.steps*86400/1e6)
Acc.N   <- c(NO3     = param$A/1e6*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.NO3.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.NO3.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.NO3.Epi"]+
                            param$h.hypo*res.11.4[1,"C.NO3.Hypo"])),
             NH4     = param$A/1e6*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.NH4.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.NH4.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.NH4.Epi"]+
                            param$h.hypo*res.11.4[1,"C.NH4.Hypo"])),
             ALG     = param$A/1e6*param$alpha.N.ALG*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.ALG.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.ALG.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.ALG.Epi"]+
                            param$h.hypo*res.11.4[1,"C.ALG.Hypo"])),
             ZOO     = param$A/1e6*param$alpha.N.ZOO*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.ZOO.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.ZOO.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.ZOO.Epi"]+
                            param$h.hypo*res.11.4[1,"C.ZOO.Hypo"])),
             POMD    = param$A/1e6*param$alpha.N.POM*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.POMD.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.POMD.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.POMD.Epi"]+
                            param$h.hypo*res.11.4[1,"C.POMD.Hypo"])),
             POMI    = param$A/1e6*param$alpha.N.POM*
                       ((param$h.epi*res.11.4[nrow(res.11.4),"C.POMI.Epi"]+
                            param$h.hypo*res.11.4[nrow(res.11.4),"C.POMI.Hypo"])-
                        (param$h.epi*res.11.4[1,"C.POMI.Epi"]+
                           param$h.hypo*res.11.4[1,"C.POMI.Hypo"])),
             POMDsed = param$A/1e6*param$alpha.N.POM*
                         (res.11.4[nrow(res.11.4),"D.POMD.Hypo"]-
                            res.11.4[1,"D.POMD.Hypo"]),
             POMIsed = param$A/1e6*param$alpha.N.POM*
                         (res.11.4[nrow(res.11.4),"D.POMI.Hypo"]-
                            res.11.4[1,"D.POMI.Hypo"]))

print("Input (tN)")
print(signif(F.in.N,3))
print("Output (tN)")
print(signif(F.out.N,3))
print(signif(sum(F.out.N),3))
print("Accumulation (tN)")
print(signif(Acc.N,3))
print(signif(sum(Acc.N),3))

print("Input - Output - Accumulation:")
print(sum(F.in.N)-sum(F.out.N)-sum(Acc.N))
```

### Task 3: Model results 

Study the time courses of the state variables and the overall fluxes of phosphorus and nitrogen resulting from running the code below and interpret the results.

### Task 4: Homework: Model sensitivity

Perform a sensitivity analysis with respect to the most important kinetic parameters and interpret the results.

```{r fig.height=10,fig.width=10}

# system.11.4@t.out = seq(0,2*365.25,by=1)   # analyze results over two years

system.11.4.sens <- calcsens(system.11.4,
                             param.sens=c("C.HPO4.in",
                                          "k.gro.ZOO",
                                          "k.death.ZOO",
                                          "k.gro.ALG",
                                          "k.death.ALG",
                                          "K.HPO4",
                                          "k.nitri",
                                          "Q.in"))

plotres(res      = system.11.4.sens,
        colnames = list("C.NH4.Epi",#"C.NH4.Hypo",
                        "C.NO3.Epi",#"C.NO3.Hypo",
                        "C.HPO4.Epi",#"C.HPO4.Hypo",
                        "C.O2.Epi","C.O2.Hypo",
                        "C.ALG.Epi",#"C.ALG.Hypo",
                        "C.ZOO.Epi",#"C.ZOO.Hypo",
                        #c("C.POMD.Epi","C.POMD.Hypo"),
                        #c("C.POMI.Epi","C.POMI.Hypo"),
                        "D.POMD.Hypo","D.POMI.Hypo"
                        ),
        col = c("black","red","blue"),
        ncol= 3)

plotres(res      = system.11.4.sens,
        colnames = list("C.NH4.Epi",#"C.NH4.Hypo",
                        "C.NO3.Epi",#"C.NO3.Hypo",
                        "C.HPO4.Epi",#"C.HPO4.Hypo",
                        "C.O2.Epi","C.O2.Hypo",
                        "C.ALG.Epi",#"C.ALG.Hypo",
                        "C.ZOO.Epi",#"C.ZOO.Hypo",
                        #c("C.POMD.Epi","C.POMD.Hypo"),
                        #c("C.POMI.Epi","C.POMI.Hypo"),
                        "D.POMD.Hypo","D.POMI.Hypo"
                        ),
        file     = "exercise_4_results_sens_a.pdf",
        width    = 12,
        height   = 12,
        col = c("black","red","blue"),
        ncol= 3)

plotres(res      = system.11.4.sens,
        colnames = list(c("C.NH4.Epi","C.NH4.Hypo"),
                        #c("C.NO3.Epi","C.NO3.Hypo"),
                        c("C.HPO4.Epi","C.HPO4.Hypo"),
                        c("C.O2.Epi","C.O2.Hypo"),
                        c("C.ALG.Epi","C.ALG.Hypo"),
                        c("C.ZOO.Epi","C.ZOO.Hypo"),
                        #c("C.POMD.Epi","C.POMD.Hypo"),
                        #c("C.POMI.Epi","C.POMI.Hypo"),
                        c("D.POMD.Hypo","D.POMI.Hypo")),
        file     = "exercise_4_results_sens_b.pdf",
        width    = 10,
        height   = 8,
        col = rainbow(6,end=0.8))

```

### Theory questions

1. Why is it important that some stoichiometric coefficients are indicated as 0/+?

2. How is the metalimnion represented in the model?

3. Look at the process rates of mineralization in the sediment. Why are they different from those in chapter 8.5?

4. Look at the mass balance for P and N. 
If there is a difference between input and output + accumulation, where does it come from? 
Hint: Have a look at the stoichiometric coefficients for anoxic mineralization.
