Carbon Sequestration

Carbon sequestration at the Laurieston

A simulation of carbon sequestration at the Laurieston site in Glasgow. The simulation is done in R over 300 Montecarlo simulations.

library(dplyr)
library(ggplot2)
library(tidyr)
library(raster)

CODE

RESULTS

tree_population_seeds <- seed_list %>% 
  dplyr::bind_rows()

  tree_population_seeds <- tree_population_seeds %>%
    dplyr::group_by(seed,ID) %>%
    dplyr::mutate(Sequestration = Carbon_storage - lag(Carbon_storage),
                  Sequestration = ifelse(Sequestration<0, 0, Sequestration),
                  'Carbon storage in standing trees' = ifelse(tree_condition>0, Carbon_storage, 0),
                  'Standing trees' = ifelse(tree_condition>0, 1, 0)) 
  
  tree_population_agr <- tree_population_seeds %>%
    dplyr::group_by(seed, Years) %>%
    summarise(across(everything(), sum)) %>%
    dplyr::mutate(Net_sequestration = Sequestration-Emission)
    
  data <- tree_population_agr %>% 
    as.data.frame() %>%
    dplyr::filter(seed==SEEDS[1]) %>%
    dplyr::select(Years, Carbon_storage, 'Carbon storage in standing trees') %>%
    dplyr::rename('Carbon storage'=Carbon_storage) %>%
    gather(key = "Indicators", value = "value", -Years) %>%
    drop_na()
  head(data, 3)

##   Years     Indicators    value
## 1     0 Carbon storage 3342.344
## 2     1 Carbon storage 3557.619
## 3     2 Carbon storage 3799.600

  ggplot(data, aes(x = Years, y = value)) + 
    geom_line(aes(color = Indicators), size = 1) +
    scale_color_manual('Indicators', values = c('red', 'blue')) +
    theme_minimal() +
    ggtitle("Carbon storage (one seed)") +
    ylab("(C)kg")

  data <- tree_population_agr %>% 
    as.data.frame() %>%
    dplyr::filter(seed==SEEDS[1]) %>%
    dplyr::select(Years, Emission, Sequestration) %>%
    gather(key = "Indicators", value = "value", -Years) %>%
    drop_na()
  head(data, 3)

##   Years Indicators    value
## 1     0   Emission 50.15033
## 2     1   Emission  3.93696
## 3     2   Emission 10.55044

  ggplot(data, aes(x = Years, y = value)) + 
    geom_line(aes(color = Indicators), size = 1) +
    scale_color_manual(values = c('red', 'blue')) +
    theme_minimal() +
    ggtitle("Sequestration vs emission (one seed)") +
    ylab("(C)kg")

  data <- tree_population_agr %>% 
    as.data.frame() %>%
    dplyr::filter(seed==SEEDS[1]) %>%
    dplyr::select(Years, 'Standing trees') %>%
    gather(key = "Indicators", value = "value", -Years) %>%
  rename(Populaiton=value)%>%
    drop_na()
  head(data, 3)

##   Years     Indicators Populaiton
## 1     0 Standing trees        101
## 2     1 Standing trees        100
## 3     2 Standing trees         96

  ggplot(data, aes(x = Years, y = Populaiton)) + 
    geom_line(aes(color = Indicators), size = 1) +
    scale_color_manual(values = c('red')) +
    theme_minimal() +
    ylab("Population (number of trees)") +
    ggtitle("Number of standing trees (one seed)")

  data <- tree_population_agr %>% 
    as.data.frame() %>%
    dplyr::filter(seed==SEEDS[1]) %>%
    dplyr::select(Years, Net_sequestration, Stored_decomposition) %>%
    gather(key = "Indicators", value = "value", -Years) %>%
    drop_na()
  head(data, 3)

##   Years        Indicators    value
## 1     1 Net_sequestration 215.2756
## 2     2 Net_sequestration 241.9807
## 3     3 Net_sequestration 185.8812

  ggplot(data, aes(x = Years, y = value)) + 
    geom_line(aes(color = Indicators), size = 1) +
    scale_color_manual(values = c('red','blue'), labels = c("Net sequestration", "Stored decomposition")) +
    theme_minimal() +
    ggtitle("Net sequestration and Stored decomposition (one seed)") +
    ylab("(C)kg")

  # Net sequestration over seeds
  data_stat <- tree_population_agr %>%
    dplyr::select(ID, Years, seed, Net_sequestration) %>%
    dplyr::group_by(Years) %>%
    dplyr::summarise(mean=mean(Net_sequestration), sd = sd(Net_sequestration), 
              cv = cv(Net_sequestration)) %>%
    drop_na()
  
  # Standard deviation band
  data <- data_stat %>%
    dplyr::mutate(Low = mean-sd,
                  High = mean+sd) %>%
    dplyr::select(Years, mean, Low, High)
    
  ggplot(data, aes(x = Years, y = mean, group = 1)) + 
    geom_line(col='red') + 
    geom_ribbon(aes(ymin = Low, ymax = High), alpha = 0.2) +
    ggtitle("Net sequestration - average over 300 seeds \n(standard deviation band)") +
    ylab("(C)kg")

  # Coefficient of variation
  
  data <- data_stat %>%
    dplyr::select(Years, cv) 
  
  ggplot(data, aes(x = Years, y = cv, group = 1)) + 
    geom_line(col='red') +
  ggtitle("Coeficient of variation of net sequestration over 300 seeds") +
    ylab("CoV(net sequestration)%")

  # Standard deviation
  data <- data_stat %>%
    dplyr::select(Years, sd)
  
  ggplot(data, aes(x = Years, y = sd, group = 1)) + 
    geom_line(col='red') +
    ggtitle("Standard deviation of net sequestration over 300 seeds") +
    ylab("SD(net sequestration)")

  # Deciles band
  mean_cl_quantile <- function(x, q = c(0.1, 0.9), na.rm = TRUE){
    dat <- data.frame(y = mean(x, na.rm = na.rm),
                      ymin = quantile(x, probs = q[1], na.rm = na.rm),
                      ymax = quantile(x, probs = q[2], na.rm = na.rm))
    return(dat)
  }
  
  data <- tree_population_agr %>%
    ungroup() %>%
    dplyr::select(Years, Net_sequestration) %>%
    drop_na() %>%
    gather(key = "Indicators", value = "value", -Years) 
  
  ggplot(data,aes(Years, value)) +
    stat_summary(geom = "line", fun = mean) +
    stat_summary(geom = "ribbon", fun.data = mean_cl_quantile, alpha = 0.2)  +
    ggtitle("Net sequestration - average over 300 seeds \n(quantile band)") +
    ylab("(C)kg")