In dynamic-R/RTM: Reaction Transport Modelling: Course Tutorials and Exercises

Answers: The NPZD model

State variables in this model are expressed in $mol~N~m^{-3}$.

In the DIN uptake flux the worker is the PHYTOplankton, while both DIN and Light are rate limiting resources.
In the Grazing flux, the ZOOplankton does the work while the PHYToplankton is the limiting resource.
Faeces can only be produced if there is grazing; a fixed part of what is ingested (pFaeces) is expelled as faeces. That is, only the part $1-pFaeces$ leads to an increase in the ZOOplankton biomass.
$Grazing - FaecesProduction$ is used for the growth of the zooplankton.
For the mortality of the zooplankton, a second-order kinetics is used ($mortalityRate \times ZOO \times ZOO$). This is often done for small organisms that are predated upon by other small organisms. It is then implicitly assumed that the predator fluctuates together with the prey. Thus, if the prey concentration is low, there will be few predators and the predation pressure will also be low. If there are many prey, there will also be many predators, so the predation pressure will be very high.

Model implementation

require(deSolve)  # package with solution methods

# state variables, units = molN/m3
state.ini  <- c(DIN = 0.010, PHYTO = 0.0005, ZOO = 0.0003, DET = 0.005) 

# parameters
pars <- c(
  depth           = 10,        # [m] depth of the bay
  rUptake         = 1.0,       # [/day]
  ksPAR           = 140,       # [uEinst/m2/s]
  ksDIN           = 1.e-3,     # [molN/m3]
  rGrazing        = 1.0,       # [/day]
  ksGrazing       = 1.e-3,     # [molN/m3]
  pFaeces         = 0.3,       # [-]
  rExcretion      = 0.1,       # [/day]
  rMortality      = 400,       # [/(molN/m3)/day]
  rMineralisation = 0.05       # [/day]
)

#=============================================================================
# Model formulation
#=============================================================================

NPZD <- function(t, state, parms) {
 with(as.list(c(state, parms)),{

    # Forcing function = Light a sine function
    # light = (540+440*sin(2*pi*t/365-1.4)), 50% of light is PAR 
    # spring starts on day 81 (22 March)
    # We calculate the PAR in the middle of the water column; extinction coefficient = 0.05/m
    PAR <- 0.5*(540+440*sin(2*pi*(t-81)/365))*exp(-0.05*depth/2)

    # Rate expressions - all in units of [molN/m3/day]
    DINuptake      <- rUptake * PAR/(PAR+ksPAR) * DIN/(DIN+ksDIN)*PHYTO
    Grazing        <- rGrazing * PHYTO/(PHYTO+ksGrazing)*ZOO
    Faeces         <- pFaeces * Grazing
    ZooGrowth      <- (1-pFaeces) * Grazing
    Excretion      <- rExcretion * ZOO
    Mortality      <- rMortality * ZOO * ZOO
    Mineralisation <- rMineralisation * DET

    # Mass balances [molN/m3/day]
    dDIN      <- Mineralisation + Excretion - DINuptake  
    dPHYTO    <- DINuptake - Grazing
    dZOO      <- ZooGrowth - Excretion - Mortality
    dDET      <- Mortality - Mineralisation + Faeces
    TotalN    <- DIN+PHYTO+ZOO+DET              # [molN/m3]

    return (list(c(dDIN, dPHYTO, dZOO, dDET),   # the derivatives
                   TotalN = TotalN, PAR = PAR)  # ordinary output variables
           )
    })
}  # end of model equations

Model run

We run the model for 2 years.

outtimes <- seq(from = 0, to = 2*365, length.out = 100)  
out <- ode(y = state.ini, parms = pars, func = NPZD, times = outtimes)  # solution
plot(out, mfrow=c(2,3))