Parameters

The function returning a named vector with parameter values:

ParmsCO2fix <- function(
      ### External conditions
      IE = 200.0,  # IE external light intensity. If we consider the value of 200 and its
      # units of J PAR / (m2 s), this means for a day-length of 12 hours
      # 8.64 MJ / (m2 d). This is a reasonable value as an average over a
      # season, but the highest values during days of the season may reach
      # 30 MJ / (m2 d) or even more.
      CDGS = 200,  # CDGS days in a growing season, d/season, which is effectively d/y.
      CDLH = 12,   # CDLH sun hours in a day, hour/day = h/d.
  
      ### Tree specific parameters
      K = 0.7,     # K Light extinction coefficient (ranges from 0.5 - 0.9) m2 soil / (m2 leaf)
      SLA = 15.0,  # Specific Leaf Area (ranges from 10 - 20) m2 leaf / (kg leaf)
      FCL = 0.5,   # Fraction of C per mass of leaf: kgC / (kg of Leaf)
      LUE = 1.E-8, # Light use efficiency (ranges from 1.(08) - 4.3(2) gDW/MJ or
      # 0.001- 0.0043 kgDW/MJ), Sinclair and Muchow, 1999).
      # However, since the unit is here given in kg CO2 / (J PAR), the
      # range is between 1.8E-9 - 7.9E-9. The LUE is taken as 1.E-8,
      # which appears rather high.
      # Photosynthesis (PG), etc. are in terms of C(arbon), and the LUE should
      # also be given in gC/MJ. Assuming 50% of C per gram of dry matter
      # (DW) this would mean a value between 0.5 * (0.001- 0.0043 kgDW/MJ)
      # = 0.0005 - 0.0022 kgC/MJ.
      CF = 0.9,    # CF dimensionless curvature factor of the NRH or photosynthesis curve.
      PC = 1.E-6,  # PC photosynthesis capacity in kg CO2/(m2 s).
  
      # Allocation fractions: (kg assimilates to organ)/(kg assimilates produced by photosynthesis)
      NF = 0.2,      # constant allocation fraction of assimilates to foliage
      NS = 0.6,      # constant allocation fraction of assimilates to stem
      NR = 1 - (NF + NS), # The fraction to roots: 1 - fraction to (foliage+stem)
  
      # RF, RR and RS are (carbon) mass-based relative respiration rates for
      # F (foliage), R (roots) and S (shoot). Unit 1/year.
      RF = 0.1, 
      RR = 0.1,
      RS = 0.1,
      # relative death rates or turnover rates, for
      # F (foliage), R (Roots), S (stem), and transition from sapwood to heartwood Unit 1/year.
      LRF = 1.0, 
      LRR = 0.2,
      LRS = 0.05, # equal for sapwood and heartwood
      TSH = 0.1,
  
      ### Management related parameters
      NH = 0.1, # NH relative harvesting rate, 1/year, NP and NT: constant allocation
      # Unit NP       ; unit NT
      # (kgC of paper) (kgC of timber wood)
      # ------------- --------------------
      # (kgC of wood) (kgC of wood)
      NP = 0.2,   # fractions of harvested product to paper, NP, and timber, NT.
      NT = 1 - NP,# fractions of harvested product to timber
      LRP = 0.5,  # relative decay rates or turnover rates, for
      LRT = 0.02, # P (Paper) and T (Timber). Unit 1/year.
      LRW = 0.1   # loss rate waste
      # The inverse of these numbers gives the average residence time of the
      # product. If RGT=0.02, this means an average residence time of 50 years.
  ){
  
  # Constants, such as molecular weights. Constants must not be changed.
  MWC <- 12 # molar weight of C, g/mol
  MWO <- 16 # molar weight of O, g/mol
  # number of seconds per hour
  CSPH <- 3600
  
  # Gather parameters in a named vector
  y <- c(IE = IE, CDGS = CDGS, CDLH = CDLH, K = K, SLA = SLA,
         FCL = FCL, LUE = LUE, CF = CF, PC = PC, 
         NF = NF, NS = NS, NR = NR,
         RF = RF, RR = RR, RS = RS, LRF = LRF, LRR = LRR, LRS = LRS, TSH = TSH,
         NH = NH, NP = NP, NT = NT, LRP = LRP, LRT = LRT, LRW = LRW,
         MWC = MWC, MWO = MWO, CSPH = CSPH)
  
  # Return
  return(y)
}

Initial conditions

The function returning a named vector with initial state values (all in units kgC m^-2), here initialization according to the fractions that are allocated to the organs during the simulation:

InitCO2fix <- function(WF = 0.02, # WF, kgC in Foliage/m2
                       WR = 0.02, # WR, kgC in Roots/m2 
                       WS = 0.06, # WS, kgC in Stems/m2
                       WH = 0.0,  # WH, kgC in Stems/m2
                       WP = 0.0,  # WP  kgC in paper/ m2 
                       WT = 0.0,  # WT, kgC in timber/m2
                       WW = 0.0   # WW, kgC in waste/m2
                       ){
  # Gather initial conditions in a named vector; given names are names for the state variables in the model
  y <- c(WF = WF,
         WR = WR,
         WS = WS,
         WH = WH,
         WP = WP,
         WT = WT,
         WW = WW)
  
  # Return
  return(y)
}

Differential equations

The function computing the rates of change of the state variables with respect to time:

RatesCO2fix <- function(t, y, parms) {
  # use the with() function to be able to access the parameters and states easily
  # remember that the with() function is with(data, expr), where
  # - 'data' is ideally a list with named elements
  #   (thus combine 'y' and 'parms' first, using the function c(), then convert to list using as.list())
  # - 'expr' is an expression (i.e. code) that is evaluated
  #   this can span multiple lines of code when embraced with curly brackets {}
  with(
    as.list(c(y, parms)),{
      
      ### Optional: forcing functions: get value of the driving variables at time t (if any)
      
      ### Optional: auxiliary equations
      LAI <- WF * SLA / FCL
      IA <- IE *(1 - exp(-K * LAI ))
      
      PGS <- (1 / (2 * CF)) * (LUE * IA + PC -
                                 sqrt((LUE * IA + PC)^2 -4 * CF * LUE * IA* PC) )
      PG <- PGS * CSPH * CDLH * CDGS * (MWC / (MWC + 2 * MWO))
      PN <- PG - (RF * WF + RR * WR + RS * WS)
      
      ### Rate equations
      DWFDT <- PN * NF - LRF * WF - NH * WF
      DWRDT <- PN * NR - LRR * WR - NH * WR
      DWSDT <- PN * NS - LRS * WS - NH * WS - TSH*WS
      DWHDT <- TSH* WS - LRS * WH - NH * WH
      
      HT <- NH * (WF + WR + WS + WH) # total harvest
      HU <- NH * (WS + WH)           # useful harvest
      
      # The useful harvest, HU, is only the stem, so: (NH*DWS).
      # This is distributed over the two states paper and timber. These are
      # fractions of the useful harvest: paper NP; timber NT=(1-NP).
      DWPDT <- HU * NP - LRP * WP # paper
      DWTDT <- HU * NT - LRT * WT # Timber
      DWWDT <- HT - HU + LRF * WF + LRR * WR +LRS *(WS+WH) - LRW * WW # Waste
      
      ### Gather all rates of change in a vector
      # - the rates should be in the same order as the states (as specified in 'y')
      # - it can be a named vector, but does not need to be
      RATES <- c(DWFDT = DWFDT,
                 DWRDT = DWRDT,
                 DWSDT = DWSDT,
                 DWHDT = DWHDT,
                 DWPDT = DWPDT,
                 DWTDT = DWTDT,
                 DWWDT = DWWDT)
      
      ### Optional: get in/out flow used to compute mass balances (or set MB <- NULL)
      # not included here (thus here use MB <- NULL), see template 3
      MB <- NULL
      
      ### Optional: gather auxiliary variables that should be returned (or set AUX <- NULL)
      # - this should be a named vector or list!
      AUX <- NULL
      
      # Return result as a list
      # - the first element is a vector with the rates of change (in the same order as 'y')
      # - all other elements are (optional) extra output, which should be named
      outList <- list(c(RATES, # the rates of change of the state variables (same order as 'y'!)
                        MB),   # the rates of change of the mass balance terms (or NULL)
                      AUX) # optional additional output per time step
      return(outList)
    })
}

With these 3 functions now defined, the CO2Fix mini-model has been implemented in R code, and the model can now numerically be solved using the ode function from the deSolve package.

Here, let’s solve the model for 3 different values of parameter CDGS: first for its default value (200; solved model stored in states), then when increasing it by 10% (solved model stored in statesHI), and then when decreasing it with 10% (solved model stored in statesLO).

states <- ode(y = InitCO2fix(),
              times = seq(from = 0, to = 50, by = 1),
              func = RatesCO2fix,
              parms = ParmsCO2fix(CDGS = 200),
              method = "ode45")
statesHI <- ode(y = InitCO2fix(),
                times = seq(from = 0, to = 50, by = 1),
                func = RatesCO2fix,
                parms = ParmsCO2fix(CDGS = 1.1*200),
                method = "ode45")
statesLO <- ode(y = InitCO2fix(),
                times = seq(from = 0, to = 50, by = 1),
                func = RatesCO2fix,
                parms = ParmsCO2fix(CDGS = 0.9*200),
                method = "ode45")

We can show these different scenarios:

plot(states, statesHI, statesLO)

where the black line is the default value for CDGS, the red line is the high value (10% increase), and the green line shows the simulation for the value reduced by 10%.

We can now calculate the sensitivity index (SI) and elasticity index (EI) of the parameter CDGS on state total carbon (Wf + Ws + Wr + Wp + Wt + Ww) at the end of the simulation:

carbonColumns <- c("WF", "WS", "WR", "WP", "WT", "WW")
deltaState <- sum(statesHI[nrow(statesHI),carbonColumns]) - sum(statesLO[nrow(statesLO),carbonColumns]) 
deltaParm <- 0.1*200
SI <- deltaState / (2*deltaParm)
SI

## [1] 0.1532905

EI <- (deltaState/sum(states[nrow(states),carbonColumns])) / ((2*deltaParm) / 200)
EI

## [1] 1.00145

This approach can be repeated for the other parameters and both state variables. Repeating the same procedure for all parameters and all state variables yields the following elasticity indices:

# Get initial conditions, parameters and the simulation output times
inits <- InitCO2fix()
parms <- ParmsCO2fix()
times <- seq(from = 0, to = 50, by = 1)

# Create empty vectors to hold the differences in the state variable (stateDiff) and parameters (parmDiff)
stateDiff <- parms # this makes of copy of the named parameter vector
parmDiff  <- parms
stateDiff[] <- NA # This sets all elements within the vector to value NA
parmDiff[] <- NA
# the result is 2 empty (values: NA) named numeric vectors (with the names of the parameters)

# Set the fraction by which each parameter is decreased/increased during sensitivity analysis
changeFraction <- 0.1 # 10%

# Set the time at which we want to evaluate the values of the state variables (a value that occurs in 'times'!)
evalTime <- max(times) # we use the last time in the output, but it can be other times as well

# In a 'for'-loop with iterator "i":
# - create 2 copies of 'parms': called 'parmsLo' and 'parmsHi'
# - reduce the value of the i-th parameter in parmsLo with 'changeFraction'
# - increase the value of the i-th parameter in parmsHi with 'changeFraction'
# - solve the ODE model for both parameter sets (they are identical except for the i-th parameter)
# - retrieve the values of the state 'stateName' after 'evalTime' time units and store the difference in 'stateDiff'
# - store the difference between the elevated and reduced parameter value in 'parmDiff'
for(i in seq_along(parms)){
  # create copies of parms
  parmsLo <- parms
  parmsHi <- parms
  
  # Reduce/increase the value of the i-th parameter
  parmsLo[[i]] = (1 - changeFraction) * parms[[i]]
  parmsHi[[i]] = (1 + changeFraction) * parms[[i]]
  
  # Solve the ODE model for both sets of parameters
  states_lo <- ode(y = inits,
                   times = times,
                   parms = parmsLo,
                   func = RatesCO2fix,
                   method = "ode45")
  states_hi <- ode(y = inits,
                   times = times,
                   parms = parmsHi,
                   func = RatesCO2fix,
                   method = "ode45")
  
  # Retrieve the values of the state variable after evalTime time units and store as the i-th element in stateDiff
  states_hi <- subset(states_hi, time==evalTime)
  states_lo <- subset(states_lo, time==evalTime)
  stateDiff[[i]] <- sum(states_hi[carbonColumns]) - sum(states_lo[carbonColumns])
  
  # Store the difference between the elevated and reduced parameter value in the i-th element of parmDiff
  parmDiff[[i]] <- parmsHi[[i]] - parmsLo[[i]]
}

We can now compute the sensitivity index:

SI <- stateDiff / parmDiff
round(SI, 5)

##            IE          CDGS          CDLH             K           SLA           FCL           LUE            CF            PC            NF            NS 
##  2.003000e-02  1.532900e-01  2.554840e+00  8.990000e-02  4.200000e-03 -1.202800e-01  4.005913e+08  2.195191e+01  2.670798e+07  1.594245e+01  2.883252e+01 
##            NR            RF            RR            RS           LRF           LRR           LRS           TSH            NH            NP            NT 
##  1.413560e+01 -4.209980e+00 -1.513376e+01 -5.417006e+01  6.896000e-02  1.307830e+00 -4.405809e+01  5.137620e+00  4.870975e+01  1.309810e+00  1.913173e+01 
##           LRP           LRT           LRW           MWC           MWO          CSPH 
## -5.290700e-01 -2.922817e+02 -9.355903e+01  1.859430e+00 -1.400880e+00  8.520000e-03

We can also calculate the elasticity indices, here rounded to 5 digits and sorted on decreasing absolute value:

states <- ode(y = inits,
                  times = times,
                  parms = parms,
                  func = RatesCO2fix,
                  method = "ode45")
Es <- parms / sum(states[nrow(states),carbonColumns]) * SI
round(Es[order(abs(Es),decreasing=TRUE)],5)

##     CSPH     CDLH     CDGS       PC      MWO      MWC       CF       NS       NT      LRW      LRT       RS       NH      LUE       IE       NF       NR 
##  1.00145  1.00145  1.00145  0.87242 -0.73216  0.72886  0.64536  0.56509  0.49995 -0.30561 -0.19095 -0.17695  0.15911  0.13085  0.13085  0.10415  0.09235 
##      LRS       RR      TSH       RF      LRP       NP      LRR      LRF        K      SLA      FCL 
## -0.07196 -0.04943  0.01678 -0.01375 -0.00864  0.00856  0.00854  0.00225  0.00206  0.00206 -0.00196

We could use a barplot to show the elasticities:

barplot(Es, col = 1:8, beside = TRUE)
legend(x = "topright", horiz = FALSE, bty = "n", legend = names(Es), pch = 22, cex = 0.7, col = 1:8)