Sensitivity Analysis of Forest Carbon Accumulation Using a Threshold-Based Growth Model

Evaluating Parameter Influence on Maximum Forest Biomass Under Climate-Driven Uncertainty

I. Background

This document explores a variance-based sensitivity analysis using the Sobol indices of a piecewise ordinary differential equation of forest growth.

a) Sobol Indices

Sobol indices are calculated to determine how sensitive the model output is to each of the parameter values. We can generate three indices (total-effect, first-order, and second-order) to help us understand how each parameter influences the model dynamics. The total-effect index (\(T_i\)) is a metric that accounts for interactions and tells us how much of the total variation in the model output can be attributed to each parameter. The first-order index (\(S_i\)) indicates how much each parameter alone impacts the output variation. The three scenarios which can arise when comparing these two indices are outlined below:

$$ \[\begin{align} \text{High}\ T_i \ \text{and Low}\ S_i &= \text{parameter interaction drives dynamics} \\ T_i \approx\ S_i &= \text{parameter acts independetly to drive dynamics} \\ \text{Low}\ T_i &= \text{parameter does not influence dynamics} \end{align}\] $$

Finally, the second-order index accounts for parameter-pair interactions and is most useful when an interaction is determined to be driving the model dynamics.

b) Forest Growth Model

In this document, we will use a piecewise forest model to simulate the maximum carbon stock over 300 years of growth.

$$

=

\[\begin{cases} r \cdot C, & \text{if } C < \text{threshold} \\ g \cdot \left(1 - \frac{C}{K} \right), & \text{if } C \geq \text{threshold} \end{cases}\]

$$

$$ \[\begin{align} \text{where:} \\ C &= \text{forest carbon stock/size (kg C)} \\ r &= \text{exponential growth rate before canopy closure (1/year)} \\ g &= \text{linear growth rate after canopy closure (kg C/year)} \\ K &= \text{carrying capacity of the forest (kg C)} \\ \text{threshold} &= \text{canopy closure level at which growth dynamics shift (kg C)} \\ t &= \text{time (years)} \end{align}\] $$

We will examine two simulations throughout this workflow. In the first simulation, we will evaluate our forest_growth function by using the average parameter values and simulating across 300 years of growth. Then, in simulation two, we will examine how uncertainty in four parameters (\(r\), \(g\), \(K\), \(threshold\)) impacts the variation in the maximum simulated forest size over time.

II. Set Up

We will begin this process by loading in the necessary libraries.

a) Load Libraries

#this clears out the environment
rm(list = ls())

#load the necessary libraries
library(calecopal)
library(knitr)
library(sensitivity)
library(tidyverse)
library(here)
library(kableExtra)
library(htmltools)
library(janitor)
library(deSolve)

b) Color Palette

Additionally, we can establish the color palette to be used throughout the document.

#create a color palette to use throughout the project
kelp1 <- cal_palette(name = "kelp1", n = 50, type = "continuous")

#create a named list of colors from the kelp1 palette
colors <- c(r = "#C88C2A",
            g = "#3D744F",
            K = "#D22600",
            threshold = "#038C9F")

III. Simulation One: No Parameter Uncertainty

Now we can begin with the first simulation. In this simulation, we will not include parameter uncertainty in the model. Instead we will use a predetermined mean value for each of the parameters.

$$ \[\begin{align} r &= \text{0.01}\ \text{year}^{-1} \\ g &= \text{2 kgC/year} \\ K &= \text{250 kgC} \\ threshold &= \text{50 kgC} \\ \end{align}\] $$

a) Forest Growth Function

First, we will load in the forest_growth function to examine the structure. This function is our ordinary differential equation (ODE) that we will use to calculate the rate of change of the forest over time.

#load the forest_growth R script
source(here("scripts", "forest_growth.R"))

#view the function
forest_growth

function (t, initial_state, params) 
{
    C <- initial_state["C"]
    r <- as.numeric(params["r"])
    g <- as.numeric(params["g"])
    K <- as.numeric(params["K"])
    threshold <- as.numeric(params["threshold"])
    if (C < threshold) {
        dCdt <- r * C
    }
    else {
        dCdt <- g * (1 - C/K)
    }
    return(list(c(dCdt)))
}

b) Initial Conditions

We can start by defining the initial conditions for the simulations. For both scenarios, we will start with an initial forest size of \(10\ kg C\) and simulate the forest growth over 300 years. Additionally, in this first scenario we will simulate the growth over 300 years and 1000 years.

#define the initial carbon size
C0 <- c(C = 10) #start with a forest of size of 10kgC

#create a times vector 
timesteps1 <- seq(from = 0, to = 300, by = 1) #run the simulation for 300 years

#create a times vector 
timesteps2 <- seq(from = 0, to = 1000, by = 1) #run the simulation for 1000 years

#define the parameter values
param_means <- c(r = 0.01, #exponential growth rate before canopy closure
                 g = 2.0, #linear growth rate after canopy closure 
                 K = 250, #carrying capacity
                 threshold = 50) #canopy closure limit 

c) Run Simulation

Then we can run the simulation and store the results in a dataframe.

#use the ode solver to run the simulation 
results1 <- ode(y = C0, #use the inital carbon stock for the initial state
               times = timesteps1,#use the timesteps for the time sequence 
               func = forest_growth, #function defined t, y, parms
               parms = param_means) #parameters to be passed to the function defined above

#store the results as a dataframe
results_df1 <- as.data.frame(results1)

#use the ode solver to run the simulation 
results2 <- ode(y = C0, #use the inital carbon stock for the initial state
               times = timesteps2,#use the timesteps for the time sequence 
               func = forest_growth, #function defined t, y, parms
               parms = param_means) #parameters to be passed to the function defined above

#store the results as a dataframe
results_df2 <- as.data.frame(results2)

d) Visualize Results

Now we can compare how the simulations look when we allow the forest to grow for either 300 years or 1000 years.

i) 300 Years

#plot the output of the forest carbon growth 
ggplot(results1, aes(x = time, y = C)) +
  geom_line(color = "forestgreen", size = 1.2) +
  geom_hline(yintercept = param_means["threshold"], linetype = "dashed", color = colors[["threshold"]], size = 0.8) +
  geom_hline(yintercept = param_means["K"], linetype = "dashed", color = colors[["K"]], size = 0.8) +
  coord_cartesian(ylim = c(0, 300)) +
  labs(title = "Forest Carbon Growth Trajectory",
       subtitle = "Simulation One: No Parameter Uncertainty", 
       x = "Time (years)",
       y = "Forest Carbon (kg C)") +
  annotate("text", x = 250, y = param_means["threshold"] + 10, 
           label = "Threshold", color = "black", size = 4, family = "Times New Roman") +
  annotate("text", x = 250, y = param_means["K"] + 10, 
           label = "Carrying Capacity (K)", color = "black", size = 4, family = "Times New Roman") +
  theme_classic() +
  theme(text = element_text(family = "Times New Roman"),
        plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
        plot.subtitle = element_text(size = 14, hjust = 0.5),
        axis.title = element_text(face = "bold", size = 14),
        axis.text = element_text(size = 10))

From this graph, we can determine that it takes over 150 years for the forest to reach the threshold of \(50\ kgC\) and switch from exponential growth to linear growth. We can also see that the forest only reaches a maximum size of less than \(200\ kgC\) by the end of 300 years. Therefore, it appears that the forest will likely not near the carrying capacity of \(250\ kgC\) in any of the runs in simulation two.

ii) 1000 Years

#plot the output of the forest carbon growth 
ggplot(results2, aes(x = time, y = C)) +
  geom_line(color = "forestgreen", size = 1.2) +
  geom_hline(yintercept = param_means["threshold"], linetype = "dashed", color = colors[["threshold"]], size = 0.8) +
  geom_hline(yintercept = param_means["K"], linetype = "dashed", color = colors[["K"]], size = 0.8) +
  coord_cartesian(ylim = c(0, 300)) +
  labs(title = "Forest Carbon Growth Trajectory",
       subtitle = "Simulation One: No Parameter Uncertainty", 
       x = "Time (years)",
       y = "Forest Carbon (kg C)") +
  annotate("text", x = 850, y = param_means["threshold"] + 10, 
           label = "Threshold", color = "black", size = 4, family = "Times New Roman") +
  annotate("text", x = 850, y = param_means["K"] + 10, 
           label = "Carrying Capacity (K)", color = "black", size = 4, family = "Times New Roman") +
  theme_classic() +
  theme(text = element_text(family = "Times New Roman"),
        plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
        plot.subtitle = element_text(size = 14, hjust = 0.5),
        axis.title = element_text(face = "bold", size = 14),
        axis.text = element_text(size = 10))

We can see that it takes the forest approximately 700 years to reach the carrying capacity of \(250\ kgC\). Once the forest reaches the canopy threshold (\(50\ kgC\)), there is an increase in the immediate growth but this slows down as the forest approaches the carrying capacity of \(250\ kgC\).

IV. Simulation Two: Parameter Uncertainty

Moving on to simulation two, we can now include parameter uncertainty in our analysis. This process involves generating parameter sample sets using the Sobol method. We will then use these Sobol sample sets to calculate the maximum size of the forest during each loop of the function.

$$ \[\begin{align} r &= \text{normally distributed with a mean of 0.01}\ \text{year}^{-1}\ \text{and a standard deviation of 0.0001}\ \text{year}^{-1} \\ g &= \text{normally distributed with a mean of 2 kgC/year and a standard deviation of 0.02 kgC/year} \\ K &= \text{normally distributed with a mean of 250 kgC and a standard deviation of 2.5 kgC} \\ threshold &= \text{normally distributed with a mean of 50 kgC and a standard deviation of 0.5 kgC} \\ \end{align}\] $$

a) Load Functions

First, we can load the three additional functions needed for simulation two.

i) Sobol Samples Function

This first function is used to transform the Sobol object into the true distributions.

#load in the sobol samples function 
source(here("scripts", "sobol_samples.R"))

#view the sobol_samples function 
sobol_samples

function (X_unit, parm_means, parm_sds) 
{
    parm_names <- names(parm_means)
    for (i in seq_along(parm_means)) {
        X_unit[, i] <- qnorm(pmax(X_unit[, i], 0), mean = parm_means[i], 
            sd = parm_sds[i])
    }
    colnames(X_unit) <- parm_names
    as.data.frame(X_unit)
    return((X_unit))
}

ii) Maximum Carbon Function

This second function is used to calculate the maximum size of the forest returned from each sample set simulation.

#load in the maximum carbon function 
source(here("scripts", "max_carbon.R"))

#view the max_carbon function 
max_carbon

function (result, threshold) 
{
    max_C <- max(result$C, na.rm = TRUE)
    idx <- which(result$C > threshold)[1]
    idx <- ifelse(is.na(idx), length(result$C), idx)
    return(list(maxC = max_C))
}

iii) Forest Wrapper Function

Thie final funciton is a wrapper function that is used to

#load in the forest wrapper function 
source(here("scripts", "forest_wrapper.R"))

#view the forest_wrapper function 
forest_wrapper 

function (r, g, K, threshold, initial_state, sim_times, ode_func, 
    max_func) 
{
    params <- list(r = r, g = g, K = K, threshold = threshold)
    state <- initial_state
    results <- ode(y = state, times = sim_times, func = ode_func, 
        parms = params)
    colnames(results) <- c("time", "C")
    results_df <- as.data.frame(results)
    metrics <- max_func(result = results_df, threshold = threshold)
    return(metrics)
}

b) Sobol Indices

Now we can use those function to calculate the Sobol indices for this forest_growth model.

i) Sampling Conditions

We will start by defining the sampling conditions for this simulation. Here we will use the same parameter distributions defined above which uses the value from simulation one as the mean and generates a normal distribution with a standard deviation of 10% of the mean.

#define the seed for reproducibility
set.seed(42)

#determine the number of samples
np <- 1000  

#determine the number of parameters (r, g, K, threshold)
k <- 4     

#calculate the parameter standard deviations (10% of mean)
param_sds <- param_means * 0.10

ii) Generate Sample Matrices

After establishing the sampling conditions, we can generate the unit-uniform sample matrix using the sobolSalt method.

#generate a sobol object using the sobolSalt method 
sobol_object <- sobolSalt(model = NULL, 
                          X1 =  matrix(runif(np * k), nrow = np, ncol = k),
                          X2 = matrix(runif(np * k), nrow = np, ncol = k),
                          scheme = "B",
                          nboot = 100)

#extract the sample matrix from the sobol object 
X_raw <- sobol_object$X

iii) Transform Samples to True Distributions

Now we will use the sobol_samples function to transform the sample sets to their true distributions.

#create a dataframe with sets of Sobol parameter values 
param_df <- sobol_samples(X_unit = X_raw, #input the raw X object 
                          parm_means = param_means, #use the param_means defined above 
                          parm_sds = param_sds) %>%  #use the param_sds defined above 
  as.data.frame(.)

iv) Maximum Forest Size

Additionally, we will now use the forest_wrapper, forest_growth, and max_carbon functions to calculate the maximum size of the forest for each sample set of parameters simulated over 300 years.

#create a list with the maximum size of the forest for each sample set simulated over 300 years
all_results <- param_df %>% #input the param_df 
  pmap(forest_wrapper, #for each row in the param_df, pass it through the forest_wrapper function 
       initial_state = C0, #define the initial state to be used
       sim_times = timesteps1, #input the time vector to be used
       ode_func = forest_growth, #call the ode function - the forest_growth function 
       max_func = max_carbon) #call the summary function used to calculate the max size of the forest for each simulation 

#create a summary dataframe from the results with the maximum carbon from each year
summary_df <- tibble(max_C = map_dbl(all_results, 1))

c) Calculate Sobol Indices

Now that we have finalized the dataframe containing the maximum size of the forest from each sample set simulation, we can calculate the sobol indices.

#update the matrix of the Sobol object to contain a full matrix of parameter sets
sobol_maxC <- tell(sobol_object, #method to estimate the sensitivity measures of the incomplete Sobol object
                   summary_df$max_C) #estimated using the forest size values output from the model 

#update the rownames for the total-effect and first-order indices
rownames(sobol_maxC$T) <- rownames(sobol_maxC$S) <- colnames(param_df)[1:4]

i) Total-Effect Index

The total-effect indices tell us the sum of the variance associated with a particular parameter and the interaction with other parameters.

#create a dataframe with the total-effect indices
sobol_T <- sobol_maxC$T %>%
  as.data.frame() %>%
  rownames_to_column("parameter") %>%
  arrange(desc(original)) %>%  #original holds total-effect index
  clean_names()  %>% 
  mutate(influential = case_when(min_c_i <= 0 & max_c_i >= 0 ~ "no", #CI includes zero so the parameter is NOT influential
                                 TRUE ~ "yes")) %>% #CI excludes zero so the parameter is influential
  select(-"std_error", -"bias")

#create a plot with the total-effect indices
ggplot(sobol_T, aes(x = parameter, y = original, fill = parameter)) +
  geom_col(alpha = 0.7) +
  geom_text(aes(label = round(original, 2)), 
            vjust = -0.5, 
            size = 4,
            family = "Times New Roman") +
  labs(title = "Total-Effect Sobol Indices (Max Carbon)",
       subtitle = "Simulation Two: Parameter Uncertainty", 
       y = "T Index",
       x = "Parameter") +
  scale_fill_manual(values = colors) +
  theme_classic() +
  theme(text = element_text(family = "Times New Roman"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        plot.subtitle = element_text(size = 14, hjust = 0.5), 
        axis.title = element_text(size = 14, face = "bold"),
        axis.text = element_text(size = 12),
        legend.position = "none")

From this analysis it appears that the linear growth rate (\(r\)) and the carrying capacity (\(K\)) are the dominant drivers of model dynamics. However, we can further examine how many of the parameters are influential as well:

#create a kable
t_kable <- kable(sobol_T,
                 col.names = c("Parameter", "Total-Effect Index", "Min C.I.", "Max C.I", "Influential ?"),
                 caption = htmltools::tags$div(style = "text-align: center; font-size: 1rem;",
                                               list(htmltools::tags$strong("Total-Effect Indices"),
                                                    htmltools::tags$br(),
                                                    htmltools::tags$span("Maximum Forest Size"))),
                  escape = FALSE) %>%
  kable_styling(full_width = FALSE, font_size = "0.8rem", position = "center") %>%
  row_spec(row = 0, bold = TRUE, background = "rgba(174, 209, 200, 0.4)") %>%
  row_spec(row = 0:4, align = "center") %>% 
  kable_classic(html_font = "Times New Roman") 

#view the kable
t_kable

Total-Effect Indices
Maximum Forest Size

Parameter	Total-Effect Index	Min C.I.	Max C.I	Influential ?
r	0.3723040	0.3212767	0.4059387	yes
K	0.3712474	0.3190687	0.4069773	yes
g	0.2285285	0.2038744	0.2498172	yes
threshold	0.0574407	0.0505793	0.0624075	yes

From this we can conclude that all of the parameters are influencing the model dynamics. To investigate if these parameters function independently or interact, we can examine the first-order indices as well.

ii) First-Order Index

The first-order indices indicate the variance associated with the single parameter alone, not including interactions. When a parameter has a high total-effect index and a low first-order index, this indicates that the parameter is mostly important through interactions. When the total-effect and first-order index are roughly equivalent, this means that the parameter acts independently more than interactive. And finally, if the total-effect index is low, this means the parameter is not very influential.

#create a dataframe with the first-order indices
sobol_S <- sobol_maxC$S %>%
  as.data.frame() %>%
  rownames_to_column("parameter") %>%
  arrange(desc(original)) %>%  #original holds total-effect index
  clean_names()  %>% 
  mutate(influential = case_when(min_c_i <= 0 & max_c_i >= 0 ~ "no", #CI includes zero so the parameter is NOT influential
                                 TRUE ~ "yes")) %>% #CI excludes zero so the parameter is influential
  select(-"std_error", -"bias")

#create a plot with the total-effect indices
ggplot(sobol_S, aes(x = parameter, y = original, fill = parameter)) +
  geom_col(alpha = 0.7) +
  geom_text(aes(label = round(original, 2)), 
            vjust = -0.5, 
            size = 4,
            family = "Times New Roman") +
  labs(title = "First-Order Sobol Indices (Max Carbon)",
       subtitle = "Simulation Two: Parameter Uncertainty",
       y = "S Index", 
       x = "Parameter") +
  scale_fill_manual(values = colors) +
  theme_classic() +
  theme(text = element_text(family = "Times New Roman"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        plot.subtitle = element_text(size = 14, hjust = 0.5), 
        axis.title = element_text(size = 14, face = "bold"),
        axis.text = element_text(size = 12),
        legend.position = "none")

We can once again confirm which parameters are influential independently by examining the maximum and minimum range of the confidence intervals.

#create a kable
s_kable <- kable(sobol_S,
                 col.names = c("Parameter", "First-Order Index", "Min C.I.", "Max C.I", "Influential ?"),
                 caption = htmltools::tags$div(style = "text-align: center; font-size: 1rem;",
                                               list(htmltools::tags$strong("First-Order Indices"),
                                                    htmltools::tags$br(),
                                                    htmltools::tags$span("Maximum Forest Size"))),
                  escape = FALSE) %>%
  kable_styling(full_width = FALSE, font_size = "0.8rem", position = "center") %>%
  row_spec(row = 0, bold = TRUE, background = "rgba(174, 209, 200, 0.4)") %>%
  row_spec(row = 0:4, align = "center") %>% 
  kable_classic(html_font = "Times New Roman") 

#view the kable
s_kable

First-Order Indices
Maximum Forest Size

Parameter	First-Order Index	Min C.I.	Max C.I	Influential ?
r	0.3536260	0.3099274	0.4147958	yes
K	0.3447861	0.3121313	0.3922842	yes
g	0.2114491	0.1543123	0.2977824	yes
threshold	0.0471592	-0.0111194	0.1326670	no

d) Visualize Forest Size

Now we can create a boxplot of the maximum forest sizes produced from each run of sample parameters.

#plot the maximum forest size 
ggplot(summary_df, aes(x = "", y = max_C)) +
  geom_boxplot(fill = "forestgreen", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Maximum Forest Carbon",
       y = "Maximum Carbon (kg C)", x = "") +
  theme_minimal() +
  theme(text = element_text(family = "Times New Roman"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        axis.title = element_text(size = 14, face = "bold"),
        axis.text = element_text(size = 12))

V. Discussion

The Sobol sensitivity analysis shows that parameters \(r\) (early exponential growth rate) and \(K\) (carrying capacity) have the largest influence on the maximum forest size, both with total-effect indices of 0.37, indicating they drive much of the output variability. Since the values of the first-order and total-effect indices are very similar, this suggests that interactions between parameters are minimal, and most of the influence comes from their independent effects. In contrast, \(threshold\) has little influence on model outcomes under the conditions simulated, implying uncertainty in this parameter contributes little to the final forest carbon stock prediction. Climate change could impact the early exponential growth rate (\(r\)) through a variety of stressors that may hinder sapling survival such an extreme changes in temperature, soil moisture, or nutrient availability. If climate change reduces the early exponential growth rate even further than it would likely take longer for the forest to reach threshold.