20  Parcel Calculations

This section calculates the mean carbon stock in total harvested biomass and the mean carbon stock in extracted timber (merchantable timber that leaves the forest) for all parcels and years.

Due to a discrepency in the gazetted project area and the actual project area, it was necessary to adjust the strata areas based on this difference - therefore all strata areas were reduced by 1.67% to align with the legal definition of the site area.

Due to the use of spatially continuous timber volume maps, as opposed to an area averaged strata approach, the parcel level calculations have been simplified so that for each parcel the sum of the overlapping portion of timber volume maps considered.

Mean carbon stock of harvested biomass per unit area (\(C_{HB}\)) was calculated, in line with Eq. 3 of the VM0010 methods, as presented in Equation 20.1 \[ C_{HB} = V_{EX} \times BCEF_R \times CF \tag{20.1}\]

Where \(V_{EX}\) is the extracted timber volume, \(BCEF_R\) is a Biomass conversion and expansion factor applicable to wood removals in the project area (for which we use a value of 1.89), and \(C_F\) is carbon fraction for which use a value of 0.47 in line with Martin and Thomas (2011).

Then the extracted carbon stock \(C_{EX}\) was calculated, in line with Eq. 4 of the VM0010 methods, as presented in Equation 20.2 where \(W_D\) is the mean wood density of all trees measured in the forest plots - this was calculated as 0.507 \(g \space{cm^{-3}}\).

\[ C_{EX} = V_{EX} \times WD \times CF \tag{20.2}\]

The change in carbon stock of dead wood as logging slash resulting from timber harvest per unit area in stratum (\(\Delta C_{DW Slash}\)) was calculated, in line with Eq. 5 of the VM0010 methods, as presented in as in Equation 20.3:

\[ \Delta C_{DW Slash} = (C_{HB}- C_{EX}) + C_{RSD} \tag{20.3}\]

Mean carbon stock in timber from residual stand damage per unit area for species (\(C_{RSD}\)) was then calculated, in line with Eq. 6 of the VM0010 methods, as presented in Equation 20.4:

\[ C_{RSD} = C_{EX} \times F_{RSD} \tag{20.4}\]

where \(F_{RSD}\) is a dimensionless factor for residual stand damage which we define as 1.8. This value was selected to align with observations from Pinard and Putz (1996) for the Kuamut site. This is more conservative than the suggested value of 2.3 on Page 56 of the VM0010 methods which was considered too high based on the local forest structure.

Carbon stock of extracted timber that is assumed to be emitted immediately at the time of harvest (\(C_{WP0}\)) was calculated, in line with Eq. 9 of the VM0010 methods, as presented in Equation 20.5:

\[ C_{WP0} = \sum_k C_{EX} \times (WW_k \cdot P_k + SLF_k\cdot P_k) \tag{20.5}\]

Where \(k\) is the wood product (for which we include: Sawn Wood, Woodbase Panels, Other Industrial Roundwood, and Paper & Paperboard), \(WW_k\) is the fraction of biomass carbon from wood waste that is assumed to be emitted to the atmosphere immediately at the time of harvest for wood product \(k\), and \(SLF_k\) is the fraction of biomass carbon from the shortlived wood product pool that is assumed to be emitted to the atmosphere immediately at the time of harvest for wood product \(k\), and \(P_k\) is the proportion of the total extracted timber that is made up of wood product \(k\). See Table 20.1 for the values used for \(WW_k\), \(SLF_k\), and \(P_k\).

Table 20.1: Wood product parameters used in the calculation of \(C_{WP0}\), \(C_{WP}\), and \(C_{WP100}\)

The carbon stock of extracted timber from stratum i that is assumed to enter the wood products pool that is not immediately emitted at the time of harvest (\(C_{WP}\)) was calculated, in line with Eq. 10 of the VM0010 methods, as presented in Equation 20.6

\[ C_{WP} = \sum_k C_{EX} - C_{WP0} \tag{20.6}\]

The Cabron stored in wood products that are assumed to be retired between 3 - 100 years after harvest from strata (\(C_{WP100}\)) was calculated, in line with Eq. 11 of the VM0010 methods, as presented in Equation 20.7:

\[ C_{WP100} = C_{WP} \times (OF_k \cdot P_k) \tag{20.7}\]

Where \(OF_k\) is the fraction of biomass carbon for wood product type \(k\) that is assumed to be emitted to the atmosphere between 3 and 100 years following timber harvest.

Table 20.2: Natural Forest Management Timber Volume estimaes for the project, by year

20.1 Code

#' Plot level estimates from Tree-level estimates of biomass
#'
#' This file computes Plot level estimates from Tree-level estimates of biomass,
#' volume and carbon  (see Kuamut_TreeAllometrics.R). Combines the THPs into
#' combined parcels, reduce the total area to match Gazetted extent
#'
#' @param V_EX_31cmDBH_csv Path to the V_EX_31cmDBH.csv
#' (coupe timber volume > 31cm DBH)
#' @param V_EX_10cmDBH_csv Path to the V_EX_10cmDBH.csv
#' (coupe timber volume > 10cm DBH)
#' @param V_EX_ACD_10cm_csv Path to the V_EX_ACD_10cm.csv
#' (Aboveground carbon density for trees > 10cm DBH)
#' @param NFM_40cmDBH_csv Path to the NFM_40cmDBH.csv
#' (NFM timber volume > 40cm DBH)
#' @param NFM_10cmDBH_csv Path to the NFM_10cmDBH.csv
#' (NFM timber volume > 10cm DBH)
#' @param NFM_ACD_10cm Path to the NFM_ACD_10cm.csv (Aboveground carbon
#' density for trees > 10cm DBH)
#' @param NoGo_src Path to the NoGo spatial vector source
#' @param BCEF_R Biomass conversion and expansion factor applicable to wood
#' removals in the project area. Default is 1.89 (see IPCC 2006)
#' @param CF default is 0.47. See IPCC 2006 (also McGroddy et al., 2004)
#' @param WWk dafault is 0.24. Fraction of biomass carbon from wood waste that
#' is assumed to be emitted to the atmosphere immediately at the time of harvest
#' for wood product k, dimensionless 24% for developing countries.
#' Winjum et al 1998 (see also page 68 in VM0010v1.3)
#' @param WD The average wood density across all plots - use the result of
#' `tree_to_plot` in `KuamutTree2Plot.R`
#' @param F_RSD Default is 1.8. For Residual Damage, we aim to match Pinard and
#' Putz (1996) and therefore this value is more conservative than the
#' methodology's suggestion of 2.3, which was considered high for Kuamut.
#' @param parent_dir The parent directory to save the outputs.
#' @return A list containing the calculated values for each parcel.
#' @details
#' Runs After Timber Harvest Plan
#' Fits with Approved VCS Methodology VM0010
#' Version 1.3, 26 April 2016, Sectoral Scope 14
#' Methodology for Improved Forest Management: Conversion from Logged to
#' Protected Forest
#'
#' # 8.1.1 Calculation of carbon stocks in commercial timber volumes
#' This section calculates CHB,j,i|BSL, the mean carbon stock in total harvested
#' biomass in tC ha-1 and CEX,j,i|BSL, the mean carbon stock in extracted timber
#' (merchantable timber that leaves the forest) in tC<U+00B7>ha-1.
#' CHB =the mean carbon stock in total harvested biomass
#' CEX =the mean carbon stock in extracted timber (merchantable timber that
#' leaves the forest)
parcel_calculations <- function(
  V_EX_31cmDBH_csv,
  V_EX_10cmDBH_csv,
  V_EX_ACD_10cm_csv,
  NFM_40cmDBH_csv,
  NFM_10cmDBH_csv,
  NFM_ACD_10cm,
  NoGo_src,
  WD,
  BCEF_R = 1.89,
  CF = 0.47,
  WWk = 0.24,
  F_RSD = 1.8,
  parent_dir = set_path("Routputs"),
  parcels_out_parent = set_path("RoutputProducts")
) {
  # save option
  to_save <- getOption("kuamut.save.outs", default = TRUE)
  # set up table directories.
  tab_path <- file.path(parent_dir, "Tables")
  if (!dir.exists(parent_dir)) {
    dir.create(parent_dir)
  }
  if (!dir.exists(tab_path)) {
    dir.create(tab_path)
  }

  Mosaic_Parcels <- read.table(V_EX_31cmDBH_csv)

  Mosaic10cm <- read.table(V_EX_10cmDBH_csv)
  Mosaic_Parcels$V10cm <- Mosaic10cm[, 4]

  MosaicACD <- read.table(V_EX_ACD_10cm_csv)
  Mosaic_Parcels$ACD <- MosaicACD[, 4]

  ## NFM Parcels
  NFM_Parcels <- read.table(NFM_40cmDBH_csv)
  NFM_Parcels$V_EX <- NFM_Parcels$V_EX #

  NFM_10cm <- read.table(NFM_10cmDBH_csv)
  NFM_Parcels$V10cm <- NFM_10cm[, 4]

  NFM_ACD <- read.table(NFM_ACD_10cm)
  NFM_Parcels$ACD <- NFM_ACD[, 4]

  Parcels <- rbind(Mosaic_Parcels, NFM_Parcels)

  ### adjust area to match gazette
  StrataArea <- dplyr::filter(Parcels, Rotation == "First") |>
    group_by(Strata) |>
    dplyr::summarise(Area_ha = sum(Area_ha, na.rm = TRUE))

  # Add no-go area
  NoGo <- sf::read_sf(NoGo_src)

  NoGoArea <- sf::st_area(NoGo) |>
    units::set_units(ha) |>
    as.numeric()
  #  add to table
  StrataArea[4, 1] <- "No-Go"
  StrataArea[4, 2] <- NoGoArea
  # StrataArea

  # compare to total area
  Total <- sum(StrataArea$Area_ha)

  GazettedArea <- (69454 + 13927)
  AreaOver <- Total - GazettedArea

  ## percentage we will take to balance out
  TakePercent <- GazettedArea / Total

  # each strata proportion of the total area
  StrataArea$Area_ha_corrected <- StrataArea$Area_ha * TakePercent

  # each strata proportion of the total area
  StrataArea$ProportionOfArea <- StrataArea$Area_ha_corrected / GazettedArea

  # -------------------------------------------------------------------------- #
  ## Save StrataArea table for PD in word doc format

  # Check table here
  colnames(StrataArea) <- c(
    "Strata",
    "Area (ha)",
    "Corrected Area (ha)",
    "Proportion Of Project Area"
  )
  StrataArea.doc <- file.path(tab_path, "StrataArea.docx")
  if (to_save) {
    gt_tab(StrataArea, .filename = StrataArea.doc, .digits = 5)
  }
  # -------------------------------------------------------------------------- #

  ## Correct all parcels accordingly (reduce by 1.67%)
  Parcels$Area_ha <- Parcels$Area_ha * TakePercent

  # double check area matches
  StrataAreaCheck <- dplyr::filter(Parcels, Rotation == "First") |>
    group_by(Strata) |>
    dplyr::summarise(Area_ha = sum(Area_ha, na.rm = TRUE))
  sum(StrataAreaCheck$Area_ha) + (NoGoArea * TakePercent)

  # -------------------------------------------------------------------------- #
  ## Save V_EX table for PD in word doc format

  # order of ... is important - last var must be the data for pivot wider.
  df_to_word_table <- function(
    df,
    ...,
    .outpath,
    .round = 2,
    .daterange = 2016:2045,
    .tabnames = c(c(
      "Year of Logging",
      "Area (ha)",
      "Mosaic plantings",
      "Mosaic unplanted",
      "NFM"
    ))
  ) {
    vars <- enquos(..., .named = TRUE)
    df <- df |>
      tibble() |>
      dplyr::filter(LogYear %in% .daterange) |>
      select(!!!vars) |>
      tidyr::pivot_wider(
        names_from = Strata,
        values_from = tail(names(vars), n = 1)
      ) |>
      dplyr::mutate(
        across(where(is.numeric), ~ round(.x, .round)),
        across(
          everything(),
          ~ dplyr::case_when(
            is.na(.x) ~ " ",
            TRUE ~ as.character(.x)
          )
        )
      )

    colnames(df) <- .tabnames

    gt_tab(df, .filename = .outpath, .digits = 6)
  }

  tab_V_EX_2045.doc <- file.path(tab_path, "tab_V_EX_2045.docx")
  if (to_save) {
    df_to_word_table(
      Parcels,
      LogYear,
      Strata,
      Area_ha,
      V_EX,
      .outpath = tab_V_EX_2045.doc
    )
  }

  # -------------------------------------------------------------------------- #
  # Equation 3 # (Page 20) in VM0010v1.3

  Parcels$C_HB <- with(
    Parcels,
    V_EX * BCEF_R * CF
  )

  # -------------------------------------------------------------------------- #
  ## Save C_HB table for PD in word doc format

  tab_C_HB_2045.doc <- file.path(tab_path, "tab_C_HB_2045.docx")
  if (to_save) {
    df_to_word_table(
      Parcels,
      LogYear,
      Strata,
      Area_ha,
      C_HB,
      .outpath = tab_C_HB_2045.doc
    )
  }

  # -------------------------------------------------------------------------- #
  # Equation 4 # (Page 21) in VM0010v1.3
  # the mean carbon stock of extracted timber per unit area (so the carbon that
  # left the forest)
  # (note is already all species together so effectively Equation 4 and Equation
  # 8 (Page 24))

  Parcels$C_EX <- with(
    Parcels,
    V_EX * WD * CF
  )

  # -------------------------------------------------------------------------- #
  ## Save C_EX table for PD in word doc format

  tab_C_EX_2045.doc <- file.path(tab_path, "tab_C_EX_2045.docx")
  if (to_save) {
    df_to_word_table(
      Parcels,
      LogYear,
      Strata,
      Area_ha,
      C_EX,
      .outpath = tab_C_EX_2045.doc
    )
  }

  # -------------------------------------------------------------------------- #

  # Equation 5 # (Page 22) in VM0010v1.3
  ## IS below because we need to estimate parameters from Eq 6 and Eq7 first!

  # Equation 6 # (Page 23) in VM0010v1.3
  # 1/0.41 pretty much matches pinard abstract

  # F_RSD Methodology suggest 2.3 - but we estimate a lower, more conservative
  # value. Instead we use  F_RSD = 1.8 from Pinard and Putz (1996) we made this
  # more conservative as methodology suggestion of 2.3 was high in our situation
  Parcels$C_RSD <- with(
    Parcels,
    C_EX * F_RSD
  )

  # -------------------------------------------------------------------------- #
  ## Save C_RSD table for PD in word doc format

  tab_C_RSD_2045.doc <- file.path(tab_path, "tab_C_RSD_2045.docx")
  if (to_save) {
    df_to_word_table(
      Parcels,
      LogYear,
      Strata,
      Area_ha,
      C_RSD,
      .outpath = tab_C_RSD_2045.doc
    )
  }

  # -------------------------------------------------------------------------- #
  # Equation 7 # (Page 23) in VM0010v1.3
  # V_INF is not used, we therefor do not estimate C_INF

  # Equation 5 # (Page 22) in VM0010v1.3
  Parcels$dC_DW_Slash <- with(
    Parcels,
    (C_HB - C_EX) + C_RSD
  )

  # -------------------------------------------------------------------------- #
  ## Save C_RSD table for PD in word doc format

  tab_dC_DW_Slash_2045.doc <- file.path(tab_path, "tab_dC_DW_Slash_2045.docx")
  if (to_save) {
    df_to_word_table(
      Parcels,
      LogYear,
      Strata,
      Area_ha,
      dC_DW_Slash,
      .outpath = tab_dC_DW_Slash_2045.doc
    )
  }

  # -------------------------------------------------------------------------- #
  # Equation 8 # (Page 24) in VM0010v1.3
  # is already all species together
  Parcels$C_EX <- Parcels$C_EX

  # -------------------------------------------------------------------------- #
  # Parameters for Equation 9 # (Page 25) in VM0010v1.3 AND Equation 11, page 26

  # # fraction of biomass carbon from wood waste that is assumed to be emitted
  # to the atmosphere immediately at the time of harvest for wood product k,
  # dimensionless
  # 24% for developing countries. Winjum et al 1998
  # (see also page 68 in VM0010v1.3)

  # fraction of biomass carbon from the short lived wood product pool that is
  # assumed to be emitted to the atmosphere immediately at the time of harvest
  # for wood product k, dimensionless

  # Oxidized Fraction (OF)
  # Winjum et al 1998 (see also page 69 in VM0010v1.3)
  # Fraction of wood products that will be emitted to the atmosphere between 3
  # and 100 years after production;

  # We need to assume PROPORTION how much goes into each ("Sawnwood",
  # "WoodbasePanels", "OtherIndRoundwood", "PaperPaperboard")
  ## CURRENTLY ASSUMING 1/3rd OF EACH WOOD PRODUCT

  # DefineParamater
  WP <- data.frame(
    k = c("Sawnwood", "WoodbasePanels", "OtherIndRoundwood", "PaperPaperboard"),
    Proportion = c(0.19, 0.56, 0.25, 0),
    SLFk = c(0.12, 0.06, 0.18, 0.24),
    OF = c(0.86, 0.98, 0.99, 0.99),
    WWk = WWk
  )

  # -------------------------------------------------------------------------- #
  ## Save WP table for PD in word doc format
  # Names written out here for neatness in PD
  tab_WP <- data.frame(
    k = c(
      "Sawn Wood",
      "Woodbase Panels",
      "Other Industrial Roundwood",
      "Paper & Paperboard"
    ),
    round(WP[, c(2:5)], 2)
  )

  tab_WP.doc <- file.path(tab_path, "tab_WP.docx")
  if (to_save) {
    gt_tab(tab_WP, .filename = tab_WP.doc, .digits = 2)
  }
  # ------- equations performed for each 'k' -----------

  # Equation 9 # (Page 25) in VM0010 v1.3  (but they miss proportion out of the
  # formula... )
  c_wp0_eq <- function(.x) {
    tibble(
      !!(.x$k) := Parcels$C_EX *
        ((.x$Proportion * .x$WWk) + (.x$Proportion * .x$SLFk))
    )
  }

  # Equation 11 (page 29 in VM0010v1.3)
  c_wp100 <- function(.x) {
    tibble(!!(.x$k) := Parcels$C_WP * ((.x$Proportion * .x$OF)))
  }

  wood_prod_carbon <- function(.name, .fun) {
    WP |>
      group_split(k) |>
      purrr::map(~ .fun(.x)) |>
      bind_cols(
        LogYear = Parcels$LogYear,
        Strata = Parcels$Strata,
        Rotation = Parcels$Rotation
      ) |>
      relocate(
        LogYear,
        Strata,
        Sawnwood,
        WoodbasePanels,
        OtherIndRoundwood,
        PaperPaperboard,
        Rotation
      ) |>
      rowwise() |>
      mutate(!!.name := sum(c_across(Sawnwood:PaperPaperboard))) |>
      ungroup()
  }

  C_WP0 <- wood_prod_carbon("C_WP0", c_wp0_eq)

  # Set for parcels too.
  Parcels$C_WP0 <- rowSums(C_WP0[, c(3:6)])

  # -------------------------------------------------------------------------- #
  ## Save WP0 table for PD in word doc format
  # Names written out here for neatness in PD

  tab_C_WP0 <- dplyr::filter(C_WP0, LogYear %in% 2016:2045)
  tab_C_WP0 <- tab_C_WP0 |>
    magrittr::set_names(c(
      "Year of Logging",
      "Stratum",
      "Sawn Wood",
      "Woodbase Panels",
      "Other Industrial Roundwood",
      "Paper & Paperboard",
      "Rotation",
      "C_WP0"
    ))

  tab_WP0_2045.doc <- file.path(tab_path, "tab_WP0_2045.docx")
  if (to_save) {
    gt_tab(tab_C_WP0, .filename = tab_WP0_2045.doc, .digits = 2)
  }
  # -------------------------------------------------------------------------- #

  # Equation 10 (page 25, VM0010v1.3)
  # The amount of extracted carbon stock that is assumed to enter the wood
  # products pool that is not immediately emitted at harvest is calculated as
  # per equation 8 below:
  Parcels$C_WP <- with(
    Parcels,
    C_EX - C_WP0
  )

  # is calculated as Equation 11 (page 29 in VM0010v1.3):
  C_WP100 <- wood_prod_carbon("C_WP100", c_wp100)

  # Therefore, the carbon stock of wood products assumed to be retired
  # between 3-100 years following harvest

  # Add together each wood product
  ## (add up each of those rows)
  Parcels$C_WP100 <- rowSums(C_WP100[, 3:6])
  # Equation 11 Completed (page 26,VM0010v1.3)

  # -------------------------------------------------------------------------- #
  ## Save WP0 table for PD in word doc format
  tab_C_WP100 <- dplyr::filter(C_WP100, LogYear %in% 2016:2045)
  colnames(tab_C_WP100) <- c(
    "Year of Logging",
    "Stratum",
    "Sawn Wood",
    "Woodbase Panels",
    "Other Industrial Roundwood",
    "Paper & Paperboard",
    "Rotation",
    "C_WP100"
  )

  tab_C_WP100_2045.doc <- file.path(tab_path, "tab_C_WP100_2045.docx")
  if (to_save) {
    gt_tab(tab_C_WP100, .filename = tab_C_WP100_2045.doc, .digits = 2)
  }
  # -------------------------------------------------------------------------- #
  # Write completed Parcel calculations next step of annualised Calculations
  parcel_strata_calcs <- file.path(
    parcels_out_parent,
    "CalulatedParcelsStrata_VM10v1.3.csv"
  )
  write.csv(Parcels, parcel_strata_calcs, row.names = FALSE)

  return(list(
    ParcelStrata = parcel_strata_calcs,
    table_exp_list = list(
      tab_V_EX_2045.doc,
      tab_C_HB_2045.doc,
      tab_C_EX_2045.doc,
      tab_C_RSD_2045.doc,
      tab_dC_DW_Slash_2045.doc,
      tab_WP.doc,
      tab_WP0_2045.doc,
      tab_C_WP100_2045.doc,
      StrataArea.doc
    ),
    wood_product_params = WP
  ))
}

Martin, Adam R, and Sean C Thomas. 2011. “A Reassessment of Carbon Content in Tropical Trees.” PloS One 6 (8): e23533.

Pinard, Michelle A., and Francis E. Putz. 1996. “Retaining Forest Biomass by Reducing Logging Damage.” Biotropica 28 (3): 278–95. https://doi.org/10.2307/2389193.