R/scenarioEmissions.R

In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis

#' Scenario biomass and emissions
#'
#' Function \code{scenarioLiveCarbonStock} estimates the amount of living biomass per time step.
#' Functions \code{scenarioDeadWood} and \code{scenarioProducts} estimate the dead and product biomass generated in each time step.
#' Finally, function \code{scenarioEmissions} estimates the carbon source and sink effects for each time step of the simulated scenario.
#'
#' @param scenario The result of function \code{\link{IFNscenario}}.
#' @param productDestination Data frame of product destination by species and diameter class
#' @param variable Either 'biomass' (default) or 'volume'
#' @param ... Additional parameters to  \code{\link{IFNproducts}}.
#'
#' @name scenarioEmissions
#'
#' @return Functions \code{scenarioLiveCarbonStock} and \code{scenarioDeadWood} return a list of two data frames with
#' biomass (dry weight or CO2, depending on \code{as.CO2}). The first data data frame contains biomass values in Mg for the target area, whereas the second contains
#' biomass values in Mg/ha.
#'
#' Function \code{scenarioProducts} can estimate both biomass or volumes (by calling function \code{\link{IFNproducts}}). Hence, and compared to the previous functions, the two data frames can correspond to volumes (m3 and m3/ha, respectively) instead of biomass values.
#'
#' Function \code{scenarioEmissions} returns a list with the following elements:
#' \itemize{
#'   \item{\code{initialDead}: Initial values of dead matter components (in Mg/ha).}
#'   \item{\code{initialProducts}: Initial values of products (in Mg/ha).}
#'   \item{\code{liveWoodChange}: Changes in live biomass components (in Mg/ha) per time step.}
#'   \item{\code{deadWoodGeneration}: Generation of dead biomass components (in Mg/ha) per time step.}
#'   \item{\code{productWoodGeneration}: Generation of products (in Mg/ha) per time step.}
#'   \item{\code{carbonBalance}: List of two frames with carbon balance details (including net emissions) in Mg/yr or Mg/ha/yr, respectively.}
#'   \item{\code{finalDead}: Final values of dead matter components (in Mg/ha).}
#'   \item{\code{finalProducts}: Final values of products (in Mg/ha).}
#' }
#' @export
#'
scenarioProducts<-function(scenario, productDestination, variable = "biomass",
                           ...) {

  variable = match.arg(variable, c("biomass", "volume"))

  products = names(productDestination)[-c(1:4)]
  if(variable=="biomass") destination = c("stumps","coarse_slash","medium_slash", "fine_slash", products)
  else destination = c("slash", products)

  treeDataCut = scenario$treeDataCut
  plotArea = scenario$plotDataSuppl$Area
  totalArea = sum(plotArea, na.rm=T)
  names(plotArea) = row.names(scenario$plotDataSuppl)
  y = IFNproducts(x= treeDataCut,
                  productDestination = productDestination,
                  variable = variable,
                  ...)
  y$Area = plotArea[y$ID]
  if(variable=="biomass") {
    for(var in destination) y[[var]] = (y[[var]]/1000)*(y$Area) # Transform to Mg
  } else {
    for(var in destination) y[[var]] = (y[[var]])*(y$Area) # Transform to m3
  }
  res = y %>% group_by(Step) %>%
    summarize_at(destination, sum, na.rm=T)
  res = as.data.frame(res)
  row.names(res) = res$Step
  res = res[,-1]
  res_ha = res
  res_ha[,destination] = res_ha[,destination]/totalArea
  if(variable=="biomass") return(list(Mg = res, Mg_ha = res_ha))
  return(list(m3 = res, m3_ha = res_ha))
}

#' @rdname scenarioEmissions
#' @param as.CO2 Flag to indicate output as kg of CO2 / ha instead of kg of dry weight / ha. Percentage of carbon per dry weight biomass by species are
#'              taken from Montero et al. (2005) (in turn, from Ibáñez et al. 2002).
#' @param area Either 'Atlantic' or 'Mediterranean' to specify allometric equations specific to the area (for Pinus pinaster)
scenarioDeadWood<-function(scenario, as.CO2 = FALSE, area = NA) {
  treeDataDead = scenario$treeDataDead
  plotArea = scenario$plotDataSuppl$Area
  names(plotArea) = row.names(scenario$plotDataSuppl)


  totalArea = sum(plotArea)

  destination = c("roots_dead", "coarse_dead", "medium_dead", "fine_dead")

  translateToDestiny<-function(Bs) {
    res = data.frame(matrix(0, nrow=nrow(Bs), ncol = length(destination)))
    names(res)<- destination
    row.names(res)<-1:nrow(res)
    #Set missing values to zero
    Bs$Roots[is.na(Bs$Roots)] = 0
    Bs$Needles[is.na(Bs$Needles)] = 0
    Bs$Leaves[is.na(Bs$Leaves)] = 0
    Bs$Bark[is.na(Bs$Bark)] = 0
    Bs$Branches[is.na(Bs$Branches)] = 0
    Bs$Stem[is.na(Bs$Stem)] = 0
    #assign
    res$roots_dead =Bs$Roots
    res$fine_dead = Bs$Leaves + Bs$Needles + Bs$Bark
    res$medium_dead = Bs$Branches
    res$coarse_dead = Bs$Stem
    return(res)
  }


  isSteps = "Step" %in% names(treeDataDead)
  if(isSteps) {
    steps = unique(treeDataDead$Step)
    numSteps = length(steps)
    res = data.frame(matrix(0, nrow = numSteps, ncol = length(destination)))
    names(res) =destination
    row.names(res) = steps
    for(i in 1:length(steps)) {
      Bs = IFNbiomass(treeDataDead[treeDataDead$Step==steps[i],],
                      as.CO2 = as.CO2, area = area)
      Bs$plotArea = plotArea[Bs$ID]
      Bs$Roots = Bs$Roots*Bs$plotArea/1000 # Translate from kg/ha to Mg
      Bs$Stem = Bs$Stem*Bs$plotArea/1000 # Translate from kg/ha to Mg
      Bs$Branches = Bs$Branches*Bs$plotArea/1000 # Translate from kg/ha to Mg
      Bs$Leaves = Bs$Leaves*Bs$plotArea/1000 # Translate from kg/ha to Mg
      Bs$Needles = Bs$Needles*Bs$plotArea/1000 # Translate from kg/ha to Mg
      Bs$Bark = Bs$Bark*Bs$plotArea/1000 # Translate from kg/ha to Mg
      prodi = translateToDestiny(Bs)
      res[i,destination] = colSums(prodi, na.rm=T)
    }
    return(list(Mg = res, Mg_ha = res/totalArea))
  } else {
    Bs = IFNbiomass(treeDataCut, DBHclasses = DBHclasslimits,
                    as.CO2 = as.CO2, area = area)
    Bs$plotArea = plotArea[Bs$ID]
    Bs$Roots = Bs$Roots*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Stem = Bs$Stem*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Branches = Bs$Branches*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Leaves = Bs$Leaves*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Needles = Bs$Needles*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Bark = Bs$Bark*Bs$plotArea/1000 # Translate from kg/ha to Mg
    prod = translateToDestiny(Bs)
    res = colSums(prod, na.rm=T)
    return(list(Mg = res, Mg_ha = res/totalArea))
  }
}

#' @rdname scenarioEmissions
scenarioLiveCarbonStock<-function(scenario, as.CO2 = FALSE, area = NA) {
  treeDataSequence = scenario$treeDataSequence
  plotArea = scenario$plotDataSuppl$Area
  names(plotArea) = row.names(scenario$plotDataSuppl)
  totalArea = sum(plotArea)

  destination = c("roots_live", "coarse_live", "medium_live", "fine_live")

  translateToDestiny<-function(Bs) {
    res = data.frame(matrix(0, nrow=nrow(Bs), ncol = length(destination)))
    names(res)<- destination
    row.names(res)<-1:nrow(res)
    #Set missing values to zero
    Bs$Roots[is.na(Bs$Roots)] = 0
    Bs$Needles[is.na(Bs$Needles)] = 0
    Bs$Leaves[is.na(Bs$Leaves)] = 0
    Bs$Bark[is.na(Bs$Bark)] = 0
    Bs$Branches[is.na(Bs$Branches)] = 0
    Bs$Stem[is.na(Bs$Stem)] = 0
    #assign
    res$roots_live =Bs$Roots
    res$fine_live = Bs$Leaves + Bs$Needles + Bs$Bark
    res$medium_live = Bs$Branches
    res$coarse_live = Bs$Stem
    return(res)
  }

  isSteps = "Step" %in% names(treeDataSequence)
  steps = unique(treeDataSequence$Step)
  numSteps = length(steps)
  res = data.frame(matrix(0, nrow = numSteps, ncol = length(destination)))
  names(res) =destination
  row.names(res) = steps
  for(i in 1:length(steps)) {
    Bs = IFNbiomass(treeDataSequence[treeDataSequence$Step==steps[i],],
                    as.CO2 = as.CO2, area = area)
    Bs$plotArea = plotArea[Bs$ID]
    Bs$Roots = Bs$Roots*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Stem = Bs$Stem*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Branches = Bs$Branches*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Leaves = Bs$Leaves*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Needles = Bs$Needles*Bs$plotArea/1000 # Translate from kg/ha to Mg
    Bs$Bark = Bs$Bark*Bs$plotArea/1000 # Translate from kg/ha to Mg
    prodi = translateToDestiny(Bs)
    res[i,destination] = colSums(prodi, na.rm=T)
  }
  return(list(Mg = res, Mg_ha = res/totalArea))
}

#' @rdname scenarioEmissions
#' @param initialDead Initial values of dead matter components (in Mg/ha)
#' @param initialProducts Initial values of product stocks (in Mg/ha)
#' @param deadHalfLives Half-live (in yrs) of dead matter components.
#' @param productHalfLives Half-live (in yrs) of products.
scenarioEmissions<-function(scenario, productDestination,
                            initialDead = c("roots_dead" = 5, "coarse_dead" = 3, "medium_dead" = 3, "fine_dead" = 1),
                            initialProducts = c("trituration" = 0.5, "poles" = 0.1, "boards" = 0.5, "construction" = 0.5, "firewood" = 0.1, "veneer" = 0.1),
                            deadHalfLives = c("roots_dead" = 20, "coarse_dead" = 10, "medium_dead" = 5, "fine_dead" = 2),
                            productHalfLives = c("trituration" = 2, "poles" = 25, "boards" = 15, "construction" = 35, "firewood" = 2, "veneer" = 25),
                            as.CO2 = FALSE, area = NA) {

  products = names(productDestination)[-c(1:4)]
  slash = c("stumps", "coarse_slash", "medium_slash", "fine_slash")
  dead =  c("roots_dead", "coarse_dead", "medium_dead", "fine_dead")
  totalArea = scenario$totalArea
  steps = scenario$steps

  cat(paste0("Calculating live carbon stocks...\n"))

  scenLiveStock = scenarioLiveCarbonStock(scenario = scenario,
                                          as.CO2 = as.CO2, area = area)

  scenStockChange = scenLiveStock$Mg_ha[2:(steps+1),] - scenLiveStock$Mg_ha[1:steps,]

  cat(paste0("Calculating dead wood...\n"))

  scenDeadWood = scenarioDeadWood(scenario = scenario,
                                  as.CO2 = as.CO2, area = area)

  cat(paste0("Calculating products...\n"))

  scenProd = scenarioProducts(scenario = scenario,
                              productDestination = productDestination,
                              as.CO2 = as.CO2, area = area)

  slashProd = scenProd$Mg_ha[,slash]
  deadWoodGeneration = scenDeadWood$Mg_ha[,dead]
  for(i in 1:nrow(slashProd)) {
    rn = row.names(slashProd)[i]
    if(rn %in% row.names(deadWoodGeneration)) {
      deadWoodGeneration[rn,] = deadWoodGeneration[rn,] + slashProd[rn,]
    }
  }
  productWoodGeneration = data.frame(matrix(0, nrow=nrow(deadWoodGeneration), ncol=length(products)))
  names(productWoodGeneration) = products
  for(i in 1:nrow(scenProd$Mg_ha)) {
    rn = row.names(scenProd$Mg_ha)[i]
    if(rn %in% row.names(productWoodGeneration)) {
      productWoodGeneration[rn,] = scenProd$Mg_ha[rn,products]
    }
  }

  deadComp = initialDead
  prodComp = initialProducts

  deadK = log(2)/deadHalfLives
  productK = log(2)/productHalfLives

  decompositionEmissions = rep(0, steps*10)
  productEndlifeEmissions = rep(0, steps*10)
  livestockVariation = rep(0, steps*10)
  deadstockGeneration = rep(0, steps*10)
  productstockGeneration = rep(0, steps*10)
  NPPSink = rep(0, steps*10)
  yearvec = 1:(steps*10)
  stepvec = rep(0, steps*10)
  for(step in 1:steps) {
    for(i in 1:10) {
      t = ((step-1)*10)+i
      stepvec[t] = step
      NPPSink[t] = sum(scenStockChange[step,]/10) + sum(as.numeric(deadWoodGeneration[step, ])/10) + sum(as.numeric(productWoodGeneration[step, ])/10)
      livestockVariation[t] = sum(scenStockChange[step,]/10)
      deadstockGeneration[t] = sum(as.numeric(deadWoodGeneration[step,])/10)
      productstockGeneration[t] = sum(as.numeric(productWoodGeneration[step,])/10)
      # if(smoothSink) {
      #   if((i<=5) && (step>1)) {
      #     p = (5+i)/10
      #     NPPSink[t] = NPPSink[t]*p + sum(scenStockChange[step-1,]/10)*(1-p)
      #   } else if((i>=6) && (step<steps)) {
      #     p = (i-6)/10
      #     NPPSink[t] = NPPSink[t]*(1-p) + sum(scenStockChange[step+1,]/10)*p
      #   }
      # }
      decompositionEmissions[t] = sum((1 - exp(-deadK))*deadComp)
      productEndlifeEmissions[t] = sum((1 - exp(-productK))*prodComp)
      deadComp = exp(-deadK)*deadComp + (as.numeric(deadWoodGeneration[step, ])/10)
      prodComp = exp(-productK)*prodComp + (as.numeric(productWoodGeneration[step,])/10)
    }
  }
  res_ha_yr = data.frame(Step = stepvec,
                   Years = yearvec,
                   LiveWoodChange = livestockVariation,
                   DeadWoodGeneration = deadstockGeneration,
                   ProductWoodGeneration = productstockGeneration,
                   NPPSink = NPPSink,
                   DecompositionEmissions = decompositionEmissions,
                   ProductEndlifeEmissions = productEndlifeEmissions,
                   EmissionNetEffect = decompositionEmissions + productEndlifeEmissions - NPPSink)
  res_yr = data.frame(Step = stepvec,
                     Years = yearvec,
                     LiveWoodChange = livestockVariation*totalArea,
                     DeadWoodGeneration = deadstockGeneration*totalArea,
                     ProductWoodGeneration = productstockGeneration*totalArea,
                     NPPSink = NPPSink*totalArea,
                     DecompositionEmissions = decompositionEmissions*totalArea,
                     ProductEndlifeEmissions = productEndlifeEmissions*totalArea,
                     EmissionNetEffect = (decompositionEmissions + productEndlifeEmissions - NPPSink)*totalArea)
  carbonBalance = list(Mg_yr = res_yr, Mg_ha_yr = res_ha_yr)

  finalDead = deadComp
  finalProducts = prodComp
  return(list(initialDead = initialDead, initialProducts = initialProducts,
              liveWoodChange = scenStockChange,
              deadWoodGeneration = deadWoodGeneration,
              productWoodGeneration = productWoodGeneration,
              carbonBalance = carbonBalance,
              finalDead = finalDead, finalProducts = finalProducts))
}