In Silviculturalist/forester: Forest Science and Forestry

knitr::opts_chunk$set(
  dev='pdf',
  collapse = TRUE,
  comment = "#>"
)

An example of a simple growth model in R.

#Carbonnier production table Oak Sweden example.
library(forester)
library(tidyverse)


#Inputs
total_age <- 30
increment <- 5
H100_Oak <- 28
stems <- 1300
silt_fraction <- 50

#empty results table
output1 <- matrix(nrow = 0, ncol = 14)
colnames(output1) <- c(
  "age",
  "dominant_height_m",
  "diameter_BA",
  "mean_height_L",
  "stems",
  "BA",
  "Volume",
  "thinned_diameter_BA",
  "thinned_stems",
  "thinned_BA",
  "thinned_volume",
  "total_thinned",
  "MAI",
  "BA_per_annum_rent"
)
output1 <- data.frame(output1)

#thinning programmes
# A thins every 5 years.
#thinning amounts subjective, copied from table XI.1

thin_BA <-
  c(
    2.2,
    4.5,
    3.9,
    3.7,
    3.2,
    2.7,
    2.3,
    2.1,
    2,
    1.9,
    1.9,
    1.9,
    1.8,
    1.8,
    1.7,
    1.7,
    1.6,
    1.5,
    1.5,
    1.4,
    1.3,
    1.3,
    0
  )

period <- 1

while (total_age < 140) {

  if (period == 1) {
    ## Initial variables.
    ## Inherited for later periods.

    #calculate height
    dominant_height_m <-
      Carbonnier_1975_Sweden_height_trajectories_Oak(
        total_age = 100,
        top_height_m = H100_Oak,
        age_2 = total_age,
        output = "Height"
      )

    #calculate mean diameter
    mean_diameter <-
      Carbonnier_1975_initial_stand_mean_diameter(stems = stems, dominant_height_m = dominant_height_m)

    #calculate Basal area from stems and mean diameter.
    BA <- stems * ((pi * ((
      mean_diameter / 100
    ) / 2) ^ 2))

    #calculate mean height
    mean_height <-
      Carbonnier_1975_mean_height_before_thinning_Oak(dominant_height_m = dominant_height_m, stems = stems)

    #calculate form factor
    form_factor <-
      Carbonnier_1975_form_factor_Oak(Lorey_mean_height = mean_height)

    #volume m3.
    volume <- mean_height * BA * form_factor



    percent_growth_BA <- 0
    total_thinned <- 0
  }

  if (period > 1) {
    if (period < 3) {
      last_thin <- thin_BA[period]
      second_last_thin <- thin_BA[period - 1]
      third_last_thin <- 0
    } else {
      last_thin <- thin_BA[period]
      second_last_thin <- thin_BA[period - 1]
      third_last_thin <- thin_BA[period - 2]
    }

    percent_growth_BA <-
      Carbonnier_1975_compound_basal_area_growth_Oak(
        SIH100 = H100_Oak,
        silt_fraction = silt_fraction,
        total_age = total_age,
        BA_before_thinning = BA,
        BA_before_thinning_understory = 0,
        mean_basal_area_diameter_after_thinning = mean_diameter2,
        mean_height_Lorey_after_thinning = mean_height2,
        stems_after_thinning = stems2,
        removal_last_thinning_BA = last_thin,
        removal_last_thinning_BA_understory = 0,
        removal_second_last_thinning_BA = second_last_thin,
        removal_second_last_thinning_BA_understory = 0,
        removal_third_last_thinning_BA = third_last_thin,
        removal_third_last_thinning_BA_understory = 0
      )


    #updating key variables if not period 1.
    stems <- stems2
    BA <- BA2 * (1 + (percent_growth_BA / 100)) ^ increment
    volume <- volume2
    mean_diameter <- mean_diameter2
    total_age <- total_age + increment

    #new dominant height
    dominant_height_m <-
      Carbonnier_1975_Sweden_height_trajectories_Oak(
        total_age = 100,
        top_height_m = H100_Oak,
        age_2 = total_age,
        output = "Height"
      )
  }

  #mean height before thinning
  mean_height <-
    Carbonnier_1975_mean_height_before_thinning_Oak(dominant_height_m = dominant_height_m,
                                                    stems = stems)

  #form factor before thinning
  form_factor <- Carbonnier_1975_form_factor_Oak(mean_height)

  #volume before thinning
  volume <- mean_height * BA * form_factor

  #mean diameter before thinning
  mean_diameter <- 2 * sqrt((BA / stems) / pi) * 100


  #thinning
  BA2 <- BA - thin_BA[period]

  BA_thinned <- thin_BA[period]

  #quotient 0.92 recommended.
  thinned_mean_diameter <- ifelse(thin_BA[period] > 0,
                                  0.92 * mean_diameter,
                                  0)

  #calculate thinned stems
  thinned_stems <- ifelse(thin_BA[period]>0,
    (BA_thinned / ((pi * ((
      thinned_mean_diameter / 100
    ) / 2) ^ 2))),
    0)

  #calculate remaining stems
  stems2 <- stems - thinned_stems

  #calculate diameter corresponding to mean basal area
  mean_diameter2 <- 2 * sqrt((BA2 / stems2) / pi) * 100

  #mean height after thinning
  mean_height2 <-
    Carbonnier_1975_mean_height_after_thinning_Oak(
      Loreys_mean_height_before_thinning = mean_height,
      basal_area_weighted_mean_diameter_cm_before_thinning = mean_diameter,
      basal_area_weighted_mean_diameter_cm_after_thinning = mean_diameter2
    )

  #form factor after thinning
  form_factor2 <- Carbonnier_1975_form_factor_Oak(mean_height2)

  #volume after thinning
  volume2 <- mean_height2 * BA2 * form_factor2

  #removed volume
  thinned_volume <- volume - volume2

  #total thinned
  total_thinned <- total_thinned+thinned_volume

  #MAI total volume (since start).
  MAI <- (volume+total_thinned) / total_age




  output2 <- c(
    "age" = total_age,
    "dominant_height_m" = dominant_height_m,
    "diameter_BA" = mean_diameter2,
    "mean_height_L" = mean_height2,
    "stems" = stems2,
    "BA" = BA2,
    "Volume" = volume2,
    "thinned_diameter_BA" = thinned_mean_diameter,
    "thinned_stems" = thinned_stems,
    "thinned_BA" = BA_thinned,
    "thinned_volume" = thinned_volume,
    "total_thinned"=total_thinned,
    "MAI" = MAI,
    "BA_percent_growth" = percent_growth_BA
  )

  output1 <- if (period > 1) {
    dplyr::bind_rows(output1, output2)
  } else {
    output2
  }

  #update period
  period <- period + 1
}

table1 <- round(output1,1)

table1

Plotting

Basal Area Annual Increment

output1 %>% ggplot(aes(x=BA,y = BA_percent_growth))+geom_point()

Basal Area

## Timeline simulator.
before <- output1$BA+output1$thinned_BA
age <- output1$age
after <- output1$BA
before<- cbind("BA"=before,"age"=age)

after <- cbind("BA"=after,"age"=(age+0.001))

BA <- data.frame(rbind(before,after))
colnames(BA) <- c("BA","age")

#timeline table XI.1
age <- c(seq(30,140,5))
age2 <- age+0.001
BA_tableXI1 <- c(12.8,13.2,13.5,13.8,14.2,14.8,15.5,16.2,16.8,17.4,17.8,18.2,18.4,18.6,18.8,18.9,19,19,19,19,19,19,20.2)
volume_before <- BA_tableXI1+thin_BA

BA_after <- data.frame(age=age2,BA=BA_tableXI1+thin_BA)
BA_before <- data.frame(age=age,BA=BA_tableXI1)

BA2 <- bind_rows(BA_before,BA_after)

#plot
BA %>% ggplot(aes(x=age,y=BA,linetype="Simulator"),color="grey30",alpha=0.9)+geom_line()+geom_line(data=BA2,aes(x=age,y=BA,linetype="Table"))

Volume

## Timeline simulator.
before <- output1$Volume+output1$thinned_volume
age <- output1$age
after <- output1$Volume
before<- cbind("volume"=before,"age"=age)

after <- cbind("volume"=after,"age"=(age+0.001))

volume <- data.frame(rbind(before,after))
colnames(volume) <- c("volume","age")

volume %>% ggplot(aes(x=age,y=volume))+
  geom_line()+
  geom_line(data=output1,aes(x=age,y=diameter_BA*3),linetype=2)+
  theme_classic()+
  labs(
  title="Oak 28 Stand Development",
  subtitle="Southern Sweden\nSoil particle content<0.06mm 50%\nNegligible understory < 6-9 m2/ha.\nThinning repeated every 5 years",
  x="Total Stand Age",
  y="Standing Volume, cu.m.",
  caption="After growth functions by Carbonnier (1975)"
  )+
  scale_y_continuous(sec.axis=sec_axis(~./3,name="Diameter corresponding to mean Basal Area, cm"))