Soil Carbon model (Allison et al. 2010)

Enzyme model of soil carbon

allison <- function(t,y,p) {
  with( as.list(c(y, p)), {
    
  death <- rdeath * M
  eprod <- renz.p * M
  eloss <- renz.l * E
  
  Vmaxu <- Vmaxu.0 * exp(- Eau / (gasconst * (temp+273) ) )
  kmu <- kmu.slope * temp + kmu.0
  assim <- Vmaxu * M * D / (kmu + D)
  
  cue <- cue.slope * temp + cue.0
  CO2 <- assim * (1-cue)
  
  Vmax <- Vmax.0 * exp( - Ea / (gasconst * (temp+273)) )
  km <- km.slope * temp + km.0
  decomp <- Vmax * E * S / (km + S)
  
  dM.dt <- assim * cue - death - eprod
  dE.dt <- eprod - eloss
  dS.dt <- inputS + death * MtoS - decomp
  dD.dt <- inputD + death * (1-MtoS) + decomp + eloss - assim
  
  return( list( c(dS.dt, dD.dt, dM.dt, dE.dt), CO2 = CO2, Total = M+E+S+D ) )
  } )
}
  
endtime <- 262800  # 30 y by hours
endtime / 24 / 365

## [1] 30

interval <- 8760 # 1 y
interval / 24

## [1] 365

t <- seq(0, endtime, by = interval)
y <- c(S = 111.876, D = 0.00144928, M = 2.19159, E = 0.0109579 )
p <- c(temp = 20, inputS = 0.0005, inputD = 0.0005, rdeath = 0.0002,
       renz.p = 5e-6, renz.l = 0.001, MtoS = 0.5, cue.0 = 0.63, cue.slope = -0.016,
       Vmax.0 = 1e8, Vmaxu.0 = 1e8, km.0 = 500, km.slope = 5,
       kmu.0 = 0.1, kmu.slope = 0.01, Ea = 47, Eau = 47, gasconst = 0.008314)
out <- ode(y, t, allison, parms=p) 

# plot in ggplot2
## restructure data
om <- melt(as.data.frame(out[,1:6]), id=1, measured = 2:5)
head(om)

##    time variable    value
## 1     0        S 111.8760
## 2  8760        S 111.8757
## 3 17520        S 111.8755
## 4 26280        S 111.8754
## 5 35040        S 111.8754
## 6 43800        S 111.8754

ggplot(om, aes(time, value) ) + geom_line() + 
  facet_wrap(~ variable, scales="free") + scale_y_continuous(expand = c(1,0) )