Based on the conceptual diagram, the mass balance equations are:
$$\frac{dPOC}{dt} = FluxPOC-HydrolysisPOC+BactMortality$$ $$\frac{dHMWC}{dt} = HydrolysisPOC - HydrolysisHMWC$$ $$\frac{dLMWC}{dt} = HydrolysisHMWC - BactUptake$$ $$\frac{dBACT}{dt} = BactUptake-GrowthRespiration-BasalRespiration-BactMortality$$
For the rate expressions, it is important to realise that bacteria are doing ALL the work to dissolve POC into smaller (HMWC) and smaller (LMWC) chunks, until they eventually take up the resource. Thus they are the worker for HydrolysisPOC, HydrolysisHMWC and BactUptake:
$$HydrolysisPOC = rPOC \times \frac{POC} {(POC + ksPOC)} \times BACT$$ $$HydrolysisHMWC = rHMWC \times \frac{HMWC} {(HMWC + ksHMWC)} \times BACT$$ $$BactUptake = rUptake \times \frac{LMWC} {(LMWC + ksLMWC)} \times BACT$$
Not all LMWC taken up by the bacteria yields new bacterial biomass, though. Part of this uptake is used to generate the energy required for biomass synthesis, which is done by respiration. This is called "activity respiration", and it is represented by the rate GrowthRespiration: $$GrowthRespiration = fLoss \times BactUptake$$ In addition, there is basal respiration (also called "maintenance respiration"), which provides energy for the basic body functions but does not lead to growth of the bacterial biomass: $$BasalRespiration = rBasal \times BACT$$
The parameter that defines the bacterial mortality has units of $(mol~C~m^{-3})^{-1}~d^{-1}$, while the bacterial mortality should have units of $mol~C~m^{-3}~d^{-1}$. This means that we need to multiply the mortality parameter with a squared concentration to obtain the correct units for the flux. Thus, the rate expression for the mortality is:\footnote{The description of the parameter as a "quadratic mortality constant" also gives this idea away.} $$BactMortality = rMort \times BACT^2$$
Modelers often use a quadratic mortality for small organisms that grow very fast and whose predator also grows very fast, so that the predator has a high concentration if the prey has a high concentration and vice versa.
Here is the R-code for the model:
require(deSolve) # package with solution methods # state variables, units = mmolC/m3 state.ini <- c(POC = 100, HMWC = 5, LMWC = 0.15, BACT = 5) # parameters pars <- c( rPOC = 0.75, # [d-1] rHMWC = 0.5, # [d-1] rUptake = 2, # [d-1] ksPOC = 100, # [mmolC m-3] ksHMWC = 5, # [mmolC m-3] ksLMWC = 0.5, # [mmolC m-3] rBasal = 0.1, # [d-1] fLoss = 0.5, # [-] rMort = 0.05, # [(mmolC m-3)-1 d-1] fluxPOC = 0.5 # [mmolC m-3 d-1] ) # Model function DetBact <- function(t, state, parms) { with (as.list(c(state, parms)), { # Rate expressions - all in units of [mmolC/m3/day] HydrolysisPOC <- rPOC * POC / (POC + ksPOC) * BACT HydrolysisHMWC <- rHMWC * HMWC / (HMWC + ksHMWC) * BACT UptakeLMWC <- rUptake * LMWC / (LMWC + ksLMWC) * BACT BasalResp <- rBasal * BACT ActivityResp <- UptakeLMWC * fLoss MortBact <- rMort * BACT^2 # Mass balances [molC/m3/day] dPOC.dt <- fluxPOC - HydrolysisPOC + MortBact dHMWC.dt <- HydrolysisPOC - HydrolysisHMWC dLMWC.dt <- HydrolysisHMWC - UptakeLMWC dBACT.dt <- UptakeLMWC - BasalResp - ActivityResp - MortBact return(list(c(dPOC.dt, dHMWC.dt, dLMWC.dt, dBACT.dt), totalC = POC + HMWC + LMWC + BACT)) }) }
And the model solution:
# output times outtimes <- seq(from = 0, to = 360, length.out = 100) # run the model for 360 days # solve this model, using the ode function from deSolve out <- ode(y = state.ini, parms = pars, func = DetBact, times = outtimes) # solution # visualise output plot(out, xlab="time (d)", ylab="mmolC/m3")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.