This workflow gives an overview of the temperature and moisture modifiers implemented in SoilR.

List of functions

Here is a list of functions that are commonly used in biogeochemical models to modify decomposition rates. This list was extracted from Table 2 of Sierra et al. (2015, JAMES 7:335-356).

\(f(X)\) Model or author name SoilR name
\(f(T) =\)
\(\left( \frac{T_{\textrm{max}}-T}{T_{\textrm{max}}-T_{\textrm{opt}}} \right)^{0.2} \exp \left( \frac{0.2}{2.63} \left( 1- \left( \frac{T_{\textrm{max}}-T}{T_{\textrm{max}}-T_{\textrm{opt}}} \right)^{2.63} \right) \right)\) Century fT.Century1
\(3.439 \exp \left( \frac{0.2}{2.63} \left( 1- \left( \frac{T_{\textrm{max}}-T}{T_{\textrm{max}}-T_{\textrm{opt}}} \right)^{2.63} \right) \left( \frac{T_{\textrm{max}}-T}{T_{\textrm{max}}-T_{\textrm{opt}}} \right)^{0.2} \right)\) Century fT.Century2
\(0.8 \exp (0.095 T_s)\) Daycent fT.Daycent1
\(0.56+(1.46 \arctan(\pi 0.0309 (T_s - 15.7)))/\pi\) Daycent fT.Daycent2
\(Q_{10}^{(T-10)/10}\) Q10 Q2, Q1.4
\(\exp \left(308.56 \left( \frac{1}{56.02}-\frac{1}{(T+273)-227.13} \right) \right)\) Lloyd and Taylor fT.LandT
\(\exp(-3.764+0.204 T (1-0.5 T/36.9))\) Kirschbaum fT.KB
\(\exp((\ln(Q_{10})/10)(T-20))\) Demeter fT.Demeter
\(\exp(-(T/(T_{\textrm{opt}}+T_{\textrm{lag}}))^{T_{\textrm{shape}}})Q_{10}^{(T-10)/10}\) Standcarb fT.Standcarb
————- ——-
\(f(W) =\)
\(\frac{1}{1+30 \exp(-8.5 W_i)}\) Century fW.Century
\(\left( \frac{W_i-b}{a-b} \right)^{d((b-a)/(a-c))} \left( \frac{W_i-c}{a-c} \right)^d\) Daycent fW.Daycent1
\(0.25+0.75(W_i)\) Demeter fW.Demeter
\((1-\exp(-(3/W_{\textrm{min}})(W_i+a)))^{b} \ \exp(-(W_i/(M_{\textrm{max}}+c))^d)\) Standcarb fW.Standacarb
\(4 W_i (1-W_i)\) if \(W_i \leq 0.5\); 1 if \(W_i >0.5\) Candy fW.Candy
\(\exp(-\exp(a-b \ W_i))\) Gompertz fW.Gompertz
\(b W_i +(1-b) W_i^2\) Myers NA
\(a W_i -b W_i^2\) Moyano fW.Moyano
\(\min[ \alpha W_i^f, \, \beta (1-W_i)^g]\) Skopp fW.Skopp

Generally, the functions are combined to produce a decomposition modifier \(\xi(t)\) in models of the form \[ \frac{d{\bf C}(t)}{dt} = {\bf I}(t) + \xi (t) \cdot {\bf A} \cdot {\bf C}(t), \] with \[ \xi(t) = \prod_i f(X_i(t)), \] where \(f(X_i)\) are functions that depend on environmental variables \(X_i\) such as temperature and moisture.

Functions in SoilR

Temperature functions

The list of temperature functions above are implemented in the SoilR package and can be used for comparing their behavior.

library(SoilR)
Temp=seq(-18,44,0.1) #Temperature range

fT1=fT.Century1(Temp)
fT2=fT.Century2(Temp)
fT3=fT.Daycent1(Temp)
fT4=fT.Daycent2(Temp)
fT5=fT.Q10(Temp)
fT6=fT.Q10(Temp,Q10=1.4)
fT7=fT.LandT(Temp)
fT8=fT.KB(Temp)
fT9=fT.Demeter(Temp)
fT10=fT.Standcarb(Temp)

To plot them

library(RColorBrewer) #Colors for plotting
pal=brewer.pal(5,"Set1") #Color palette
par(mar=c(4,4,1,0.5))
plot(Temp,fT1,type="l",xlab="Temperature (Celcius)", ylab="f(T)",ylim=c(0,10),col=pal[1])
lines(Temp,fT2,col=pal[1],lty=2)
lines(Temp,fT3,col=pal[2])
lines(Temp,fT4,col=pal[2],lty=2)
lines(Temp,fT5,col=pal[3])
lines(Temp,fT6,col=pal[3],lty=2)
lines(Temp,fT7,col=pal[4])
lines(Temp,fT8,col=pal[4],lty=2)
lines(Temp,fT9,col=pal[5])
lines(Temp,fT10,col=pal[5],lty=2)
abline(h=1,lty=2)
abline(h=0,lty=2)
legend("topleft",c("Century1","Century2","Daycent1","Daycent2","Q10 = 2.0","Q10 = 1.4","Loyd and Taylor","Kirschbaum","Demeter","Standcarb"),lty=rep(c(1,2),5),col=rep(pal,each=2),bty="n")

plot of chunk unnamed-chunk-2

Moisture functions

The moisture dependend functions are more tricky to compare because the original authors express moisture using different metrics. Some use soil water content, others the aridity index, and others soil water potential. To deal with this, we can take the range of these metrics and expressed in a range from 0 to 1. In the JAMES paper we called this metric the moisture index \(W_i\). Some moisture functions implemented in SoilR follow the original convention reported in the original papers, so we need to re-implement them here following the \(W_i\) convention:

# Moisture Functions
W1=seq(0.01,1,0.01)

#Standarize some functions to 'Moisture Index'.
fW.Century1=function (MoistIndex) 
{
  1/(1 + 30 * exp(-8.5 * (MoistIndex)))
}

#Need to modify the Daycent function so it doesn't calculate water filled pored space. It is replaced by the Moisture Index (MoistIndex \in 0,1)
fW.Daycent1.1=function (MoistIndex, a = 0.6, b = 1.27, c = 0.0012, d = 2.84) 
{
  wfps = MoistIndex
  wfunc = (((wfps - b)/(a - b))^(d * ((b - a)/(a - c))))*((wfps - c)/(a - c))^d
  return(wfunc)
}

fW.Standcarb1=function (Moist, MatricShape = 5, MatricLag = 0, MoistMin = 30, 
                        MoistMax = 350, DiffuseShape = 15, DiffuseLag = 4) 
{
  IncreaseRate = 3/MoistMin
  MatricLimit = (1 - exp(-IncreaseRate * (Moist + MatricLag)))^MatricShape
  DiffuseLimit = exp(-1*(Moist/(MoistMax + DiffuseLag))^DiffuseShape)
  MoistDecayIndex = MatricLimit * DiffuseLimit
  return(MoistDecayIndex)
}

fW.Candy=function(MoistIndex){
  Mi=MoistIndex # Moisture index. Original equation arrives here dividing VWC by pore volume
  fw=ifelse(Mi<=0.5, 4*Mi*(1-Mi),1)
  return(fw)
}

fW.Myers=function(MoistIndex,b=2){
  b*MoistIndex+((1-b)*MoistIndex^2)
}

#Run the functions
fW1=fW.Century1(W1)
fW2=fW.Daycent1.1(W1)
fW3=fW.Demeter(W1,Msat=1)
fW4=fW.Standcarb1(W1*100,MoistMin=15,MoistMax=90)
fW5=fW.Candy(W1)
fW6=fW.Gompertz(W1)
fW7=fW.Myers(W1)
fW8=fW.Moyano(W1)
fW9=fW.Skopp(W1)

And now we can plot them

plot(W1,fW1,type="l",xlab="Moisture index",ylab="f(W)",ylim=c(0,2),col=pal[1])
lines(W1,fW2,col=pal[2])
lines(W1,fW3,col=pal[3])
lines(W1,fW4,col=pal[4])
lines(W1,fW5,col=pal[5])
lines(W1,fW6,col=pal[1],lty=2)
lines(W1,fW7,col=pal[2],lty=2)
lines(W1,fW8,col=pal[3],lty=2)
lines(W1,fW9,col=pal[4],lty=2)
legend("topleft",c("Century","Daycent","Demeter","Standcarb","Candy"),lty=rep(1,5),col=pal[1:5],bty="n")
legend("topright",c("Gompertz","Myers","Moyano","Skopp"),lty=rep(2,4),col=pal[1:4],bty="n")

plot of chunk unnamed-chunk-4

Implement your own function

If you need a different modifier not included in SoilR, you can implement your own as a normal R function. For example, let's say you want to implement an exponential function as a temperature modifier with the form

\[ f(T) = A \cdot e^{k \cdot T} \]

Then, you implment it as

fT.Exponential = function(Temp, A, k){
    A*exp(k*Temp)
}