library(animation)
library(plotrix)
library(ggplot2)

1 Introduction

1.1 Overview

Forests play a crucial role in maintaining ecological balance, supporting biodiversity, and mitigating climate change by sequestering carbon dioxide. However, forest ecosystems are increasingly threatened by extreme weather events, including droughts, fires, and heavy rainfall caused by phenomena such as La Nina. This project aims to simulate the growth of forests in the Blue Mountains, New South Wales, Australia, over a period of 5 years, considering the impact of these extreme weather events.

##Motivation

The primary motivation behind this project is to understand how various climatic factors and events affect forest growth, mortality, and biomass. By simulating these scenarios, we can gain insights into how forests might respond to changing climatic conditions, which is critical for developing effective forest management and conservation strategies.

1.2 Parameters and Data Sources

The parameters used in this simulation are based on scientific literature and historical climate data for New South Wales, Australia.

1.3 General Parameters

Number of Trees (maxTrees): 100,000. This common upper limit is based on models like SORTIE and JABOWA (SORTIE, JABOWA models).

Forest Size (maxx, maxy): 200 x 200. This size creates a manageable yet realistic forest area and is used in models like JABOWA (SORTIE, JABOWA models).

Maximum Radius of a Tree (maxr): 10 meters. Models like SORTIE use a maximum tree radius to limit growth based on empirical observations (SORTIE, JABOWA models).

Growth Rate (default_gamma): 0.2. This rate accounts for the typical growth of trees (Botkin, Janak, & Wallis, 1972).

Mortality Probability (default_mu): 0.02. This probability is derived from long-term forest monitoring data (Shugart, 1984).

Reproduction Probability (lambda): 0.1. This probability is derived from seed dispersal studies used by SORTIE and JABOWA (Clarke et al., 1998).

1.4 Event Simulation Parameters

La Nina Probability (laNinaProb): 0.2. This is based on the frequency of La Nina events in Australia from the Bureau of Meteorology.

Drought Probability (droughtProb): 0.1. This probability is based on a 10% annual probability in Australia as stated by the Bureau of Meteorology.

Fire Probability (fireProb): 1/10. This represents the typical rate of fire in eucalyptus forests (Bradstock & Gill, 2001).

Fire Mortality Rate (fireMortalityRate): 0.8. This rate is based on studies on the impact of fire on tree mortality (Bradstock & Gill, 2001).

# All parameters bundled into a list
p <- list()
p$maxTrees <- 10^5     # Number of trees (alive or dead) that the simulation can handle
p$maxx <- 200          # x-size of forest
p$maxy <- 200          # y-size of forest
p$maxr <- 10           # Maximum radius of a tree
p$default_gamma <- 0.2 # Default growth per time step
p$default_mu <- 0.02   # Default mortality probability per time step
p$lambda <- 0.1        # Probability of producing surviving offspring after reaching maturity, per time step

# New parameters for extensions:
p$laNinaProb <- 0.2       # Probability of La Nina event in a given year (20% of years)
p$droughtProb <- 0.1      # Probability of drought in a given year (10% of years)
p$fireProb <- 1 / 10      # Probability of fire event in a given year (every 10 years)
p$fireMortalityRate <- 0.8 # Tree mortality rate during fire events (80%)

# Set default growth and mortality rates
p$gamma <- p$default_gamma
p$mu <- p$default_mu

tmax <- 60 # Number of time steps to simulate (5 years with 12 steps per year)
nini <- 10  # Initial number of trees

# Initial random forest
iniForest <- data.frame(ID = 1:nini,                             # ID for each tree 
                        alive = TRUE,                            # Living or dead?
                        x = runif(nini, min = 0, max = p$maxx),  # x position
                        y = runif(nini, min = 0, max = p$maxy),  # y position
                        r = runif(nini, min = 0, max = p$maxr),  # Radius
                        a = NA)                                  # Age
iniForest$a <- iniForest$r / p$gamma # Expected age based on their size
# Core simulation function with extended parameters:
simulateForest <- function(iniForest, tmax, p) {
  ftt <- as.vector(rep(NA, tmax + 1), "list") # "forest through time"
  ftt[[1]] <- iniForest
  nextID <- max(ftt[[1]]$ID, na.rm = TRUE) + 1
  
  for (i in 1:tmax) {
    # Reset growth and mortality rates to default
    p$gamma <- p$default_gamma
    p$mu <- p$default_mu
    
    # Apply climate effects (La Nina, drought)
    if (runif(1) < p$laNinaProb) {
      p$gamma <- p$gamma * 1.3  # Increase growth rate by 30% during La Nina
    }
    if (runif(1) < p$droughtProb) {
      p$gamma <- p$gamma * 0.5  # Decrease growth rate by 50% during drought
      p$mu <- p$mu * 1.5        # Increase mortality rate during drought
    }
    
    # Simulate fire events
    if (runif(1) < p$fireProb) {
      fireImpactRadius <- runif(1, min = 10, max = 50)  # Random fire impact radius
      fireCenterX <- runif(1, min = 0, max = p$maxx)
      fireCenterY <- runif(1, min = 0, max = p$maxy)
      affected <- sqrt((ftt[[i]]$x - fireCenterX)^2 + (ftt[[i]]$y - fireCenterY)^2) < fireImpactRadius
      ftt[[i]]$alive[affected] <- runif(sum(affected)) > p$fireMortalityRate
    }
    
    ftt[[i + 1]] <- data.frame(ID = rep(NA, p$maxTrees), alive = NA,                 
                               x = NA, y = NA, r = NA, a = NA)                  
    # Forest in new time step copied over from previous one:
    ftt[[i + 1]][1:nrow(ftt[[i]]),] <- ftt[[i]] 
    # Indices of living trees:
    alive <- which(ftt[[i + 1]]$alive)
    # Aging all trees:
    ftt[[i + 1]]$a[alive] <- ftt[[i + 1]]$a[alive] + 1
    # Growing trees:
    canGrow <- which(ftt[[i + 1]]$alive & ftt[[i + 1]]$r < p$maxr)
    ftt[[i + 1]]$r[canGrow] <- ftt[[i + 1]]$r[canGrow] + p$gamma
    # Make sure they are not growing beyond maxr:
    ftt[[i + 1]]$r[ftt[[i + 1]]$r > p$maxr] <- p$maxr
    # Death:
    ftt[[i + 1]]$alive[alive] <- (runif(length(alive)) > p$mu)
    # Reproduction:
    emptySlot <- which(is.na(ftt[[i + 1]]$ID))[1] # Next empty row in dataframe
    for (j in alive) {
      if (identical(ftt[[i + 1]]$r[j], p$maxr) && (runif(1) < p$lambda)) { # Reproduction happening?
        ftt[[i + 1]]$ID[emptySlot] <- nextID
        ftt[[i + 1]]$alive[emptySlot] <- TRUE
        ftt[[i + 1]]$x[emptySlot] <- runif(1, min = 0, max = p$maxx)
        ftt[[i + 1]]$y[emptySlot] <- runif(1, min = 0, max = p$maxy)
        ftt[[i + 1]]$r[emptySlot] <- 0
        ftt[[i + 1]]$a[emptySlot] <- 0
        nextID <- nextID + 1
        emptySlot <- emptySlot + 1
      }
    }
    # Get rid of empty rows:
    ftt[[i + 1]] <- ftt[[i + 1]][!is.na(ftt[[i + 1]]$ID),]
  }
  return(ftt)
}

# Plotting the forest at a given time:
plotForest <- function(forest, p) {
  if (!is.data.frame(forest)) {
    stop("Forest is not a data frame")
  }
  if (!"alive" %in% colnames(forest)) {
    stop("Column 'alive' not found in the forest data frame")
  }
  forest <- forest[order(forest$alive, decreasing = TRUE),]
  plot(c(0, p$maxx), c(0, p$maxy), type = 'n', axes = FALSE, ann = FALSE)
  cols <- ifelse(forest$alive, hsv(0.15 + (forest$ID %% 11) / 30, 
                                   0.8 + (forest$ID %% 5) / 40, 
                                   0.4 + (forest$ID %% 7) / 30), "grey")
  for (i in 1:nrow(forest)) {
    draw.circle(forest$x[i], forest$y[i], forest$r[i], col = cols[i], border = cols[i])
  }
}

# Plotting forest growth through time as a video:
forestMovie <- function(ftt, p, interval = 0.1, movie.name) {
  ani.options(interval = interval)
  saveGIF({
    for (i in 1:length(ftt)) {
      dev.hold()
      plotForest(ftt[[i]], p)
      text(x = 0, y = 1, adj = c(0, 1), label = paste0("time=", i), col = "red", cex = 1.5)
      ani.pause()
    }
  }, movie.name = movie.name)
  invisible(NULL)
}

# Function to plot results side by side
plotResults <- function(simResults, scenario, p, tmax) {
  par(mfrow=c(1, 3), mar=c(4, 4, 2, 1))
  for (i in 1:length(simResults)) {
    plotForest(simResults[[i]][[tmax + 1]], p)
    title(main = paste(scenario, "Replicate", i), cex.main = 1.5)
  }
  par(mfrow=c(1, 1))
}

# Function to plot time series data with custom x-axis labels
plotTimeSeries <- function(data, title, ylab) {
  par(mfrow=c(1, 3), mar=c(4, 4, 2, 1))
  for (i in 1:ncol(data)) {
    plot(data[,i], type='l', xlab='Time (years)', ylab=ylab, main=paste0(title, " Replicate ", i), xaxt='n')
    axis(1, at = seq(1, length(data[,i]), by = 12), labels = 2024:(2024 + length(seq(1, length(data[,i]), by = 12)) - 1), las = 2)
  }
  par(mfrow=c(1, 1))
}

# Function to run simulations and collect data for plotting
runSimulation <- function(simName, p, iniForest, tmax, output_dir, n_replicates = 3) {
  results <- list()
  treeDensity <- matrix(NA, nrow = tmax + 1, ncol = n_replicates)
  avgTreeRadius <- matrix(NA, nrow = tmax + 1, ncol = n_replicates)
  totalBiomass <- matrix(NA, nrow = tmax + 1, ncol = n_replicates)
  mortalityRate <- matrix(NA, nrow = tmax + 1, ncol = n_replicates)
  
  for (rep in 1:n_replicates) {
    simResult <- simulateForest(iniForest, tmax, p)
    results[[rep]] <- simResult
    
    treeDensity[,rep] <- sapply(simResult, function(forest) sum(forest$alive))
    avgTreeRadius[,rep] <- sapply(simResult, function(forest) mean(forest$r[forest$alive], na.rm = TRUE))
    totalBiomass[,rep] <- sapply(simResult, function(forest) sum(forest$r[forest$alive], na.rm = TRUE))
    mortalityRate[2:(tmax + 1),rep] <- sapply(2:length(simResult), function(i) sum(!simResult[[i]]$alive & simResult[[i-1]]$alive))
  }
  
  # Plot results
  plotResults(results, simName, p, tmax)
  
  # Create GIF for each replicate
  gif_paths <- c()
  for (rep in 1:n_replicates) {
    gif_path <- file.path(output_dir, paste0(simName, "_replicate", rep, "_forest_growth.gif"))
    forestMovie(results[[rep]], p, interval = 0.1, movie.name = gif_path)
    gif_paths <- c(gif_paths, gif_path)
  }
  
  # Plot time series data
  plotTimeSeries(treeDensity, paste0(simName, ": Tree Density"), "Tree Density")
  plotTimeSeries(avgTreeRadius, paste0(simName, ": Average Tree Radius"), "Average Tree Radius")
  plotTimeSeries(totalBiomass, paste0(simName, ": Total Forest Biomass"), "Total Forest Biomass")
  plotTimeSeries(mortalityRate, paste0(simName, ": Mortality Rate"), "Mortality Rate")
  
  return(gif_paths)
}

2 Baseline Simulation

2.1 Simulation and Results

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
# Set default growth and mortality rates
p$gamma <- p$default_gamma
p$mu <- p$default_mu

# Set the output directory
output_dir <- "/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024"

# Run Baseline Simulation
p$laNinaProb <- 0
p$droughtProb <- 0
p$fireProb <- 0
gif_paths <- runSimulation("baseline", p, iniForest, tmax, output_dir)
## Output at: /Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/baseline_replicate1_forest_growth.gif
## Output at: /Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/baseline_replicate2_forest_growth.gif
## Output at: /Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/baseline_replicate3_forest_growth.gif

for (gif_path in gif_paths) {
  cat(sprintf('<img src="%s" height="400px" width="400px" autoplay loop>', gif_path))
}
## <img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/baseline_replicate1_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/baseline_replicate2_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/baseline_replicate3_forest_growth.gif" height="400px" width="400px" autoplay loop>

2.2 Explanation and Inferences

In the baseline simulation, it is observed that tree density as well as tree biomass experience fluctuations but recover well over the 5 year period. On the contrary tree radii seem to have brief periods of increase, but always followed by an inevitable period of decrease in radii. Finally, mortality rates appear to be steady with some small fluctuations, but do not indicate of external stressors trending towards forest decline.

These results are informative because they present the current dynamics of this stochastic model of a “made-up” forest at baseline. The overarching statement regarding this simulation, is that the forest appears to be healthy. This is firstly indicated by the growth in tree biomass coupled with the proportional increase of tree density, implying that the two are positively correlated. It can be seen that increase in one corresponds to increase in the other. This could imply that the new forest is young and that there is a high reproduction rate, which contributes to new trees growing rather than older trees growing larger.

The above notion is further seconded by an additional observed relationship. Visually again, a relationship can be seen between tree density and the average tree radius, where there appears to be an inverse relationship, such that a peak in tree density corresponds to a trough in the average of tree radius. This informs us that as tree density increases, the average radius begins to decline, implying that tree density may correspond to an increase in competition for resources amongst proximal trees, which would slow down their growth.

3 Impact of La Nina

3.1 Simulation and Results

# Impact of La Nina
p$laNinaProb <- 0.2
p$droughtProb <- 0
p$fireProb <- 0

gif_paths <- runSimulation("laNina", p, iniForest, tmax, output_dir)

for (gif_path in gif_paths) {
  cat(sprintf('<img src="%s" height="400px" width="400px" autoplay loop>', gif_path))
}
## <img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/laNina_replicate1_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/laNina_replicate2_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/laNina_replicate3_forest_growth.gif" height="400px" width="400px" autoplay loop>

3.2 Explanation and Inferences

During La Nina events, the growth rate of trees increases by 30%, leading to a noticeable boost in tree density, average tree radius, and total forest biomass. The enhanced growth rate accelerates forest development, but the overall mortality rate remains similar to the baseline scenario, suggesting that La Nina has a predominantly positive impact on forest growth.

Despite the favorable growth conditions, a trend similar to the baseline simulation is seen in the relationship of tree density and average tree radius. Specifically, both are negatively correlated, such that an increase in tree density leads to a decrease in tree radius, indicating that very high density could still dampen the radius growth due to competition. It could be stipulated that tree growth occurs in a manner described in the following graph.

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
# Create a data frame for the plot
resources <- seq(0, 10, by = 0.1)
tree_growth <- 1 / (1 + exp(-resources + 5))

# Plot the graph
plot(resources, tree_growth, type = 'l', lwd = 2, col = 'blue', xlab = 'Resources', ylab = 'Tree Growth', main = 'Tree Growth vs Resources')

# Adding some aesthetic features similar to the provided image
lines(resources, tree_growth, col = "blue", lwd = 2)

The sigmoidal curve above, summarizes the potential dynamics of a forest, such that the presence of resources can rapidly accelerate forest growth in one or more metrics, e.g. tree radius or tree density. Usually, the types of resources present would favor one or the other, in the case of our model, that being the rapid increase of tree density. This hypothetical graph might also explain why although La Nina provides favorable growth conditions, we are not seeing an extremely higher increase in tree density or tree radius, which could be due to the pre-existing soil conditions, for example, providing enough resources to pre-saturate its growth. As such, even La Nina might not provide enough of an additive effect to see a marked difference.

4 Impact of Drought

4.1 Simulation and Results

p$laNinaProb <- 0
p$droughtProb <- 0.1
p$fireProb <- 0
gif_paths <- runSimulation("drought", p, iniForest, tmax, output_dir)

for (gif_path in gif_paths) {
  cat(sprintf('<img src="%s" height="400px" width="400px" autoplay loop>', gif_path))
}
## <img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/drought_replicate1_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/drought_replicate2_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/drought_replicate3_forest_growth.gif" height="400px" width="400px" autoplay loop>

4.2 Explanation and Inferences

Due to the adverse conditions, typically drought events reduce the growth rate of tree density and you see a gradual but slow increase in tree radius implying the limited resources reduce forest regeneration and existing trees grow larger and older in accordance with the resources present. From the graphs presented, although a clear trend cannot be seen, there is a higher degree of overall fluctuations, indicating that drought drastically changes the dynamics of the forest, switching from trend towards growth, to a trend towards survival.

5 Impact of Fire Events

5.1 Simulation and Results

p$laNinaProb <- 0
p$droughtProb <- 0
p$fireProb <- 1 / 10
gif_paths <- runSimulation("fire", p, iniForest, tmax, output_dir)

for (gif_path in gif_paths) {
  cat(sprintf('<img src="%s" height="400px" width="400px" autoplay loop>', gif_path))
}
## <img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/fire_replicate1_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/fire_replicate2_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/fire_replicate3_forest_growth.gif" height="400px" width="400px" autoplay loop>

5.2 Explanation and Inferences

Fire events typically cause a significant increase in tree mortality. There is a sharp decline in tree density and total biomass following each fire event. The average tree radius also decreases due to the loss of mature trees. The data illustrates the severe impact of fires on forest ecosystems. More specifically, the sharp reduction in tree density directly impacts the average tree radius, as fires predominantly kill larger, mature trees. Finally, a strong relationship is expected between density and total biomass, as with fire events would cause a decrease in both metrics.

6 Combined Impact of La Nina, Drought, and Fire

6.1 Simulation and Results

p$laNinaProb <- 0.2
p$droughtProb <- 0.1
p$fireProb <- 1 / 10
gif_paths <- runSimulation("combined", p, iniForest, tmax, output_dir)

for (gif_path in gif_paths) {
  cat(sprintf('<img src="%s" height="400px" width="400px" autoplay loop>', gif_path))
}
## <img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/combined_replicate1_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/combined_replicate2_forest_growth.gif" height="400px" width="400px" autoplay loop><img src="/Users/nicholaskarasavvidis/Desktop/UQ Assignment Project/UQ 23052024/combined_replicate3_forest_growth.gif" height="400px" width="400px" autoplay loop>

6.2 Explanation and Inferences

The combined simulation includes all weather effects of La Nina, drought, and fire events. The positive impact of La Nina is seen to be counterbalanced by the negative effects of drought and fire. While in isolation La Nina boosts growth, droughts and fires cause substantial tree mortality and therefore cause the reduction of biomass. This scenario provides a realistic representation of forest growth under fluctuating and dynamic climatic conditions. More specifically, the combined scenario shows a balance, with both tree density and average tree radius increasing during favorable conditions (La Nina) and decreasing during adverse events (drought and fire). Finally, the total biomass reflects the forest’s overall health, with fluctuations indicating the forest’s resilience and vulnerability to multiple climatic factors.

7 Conclusion

7.1 Summary

This study simulated forest growth in the Blue Mountains, New South Wales, Australia, considering various extreme weather events. The simulations revealed the following key points:

Baseline Scenario: Demonstrated natural forest growth and stable conditions.

La Nina: Increased growth rates and enhanced forest development.

Drought: Reduced growth rates and increased tree mortality.

Fire: Caused significant tree mortality and loss of biomass.

Combined Scenario: Showed a balance between growth and mortality, reflecting the complexity of real-world forest dynamics.

These findings underscore the importance of understanding and mitigating the impacts of extreme weather events on forest ecosystems to ensure their resilience and sustainability.

8 References

  1. Botkin, D. B., Janak, J. F., & Wallis, J. R. (1972). Some ecological consequences of a computer model of forest growth. Journal of Ecology, 60(3), 849-872. https://doi.org/10.2307/2258570

  2. Bradstock, R. A., & Gill, A. M. (2001). Living with fire and biodiversity. Australian Journal of Botany, 49(3), 279-289. https://doi.org/10.1071/BT01001

  3. Bureau of Meteorology, Australia. (n.d.). Retrieved from http://www.bom.gov.au

  4. Clark, J. S., Fastie, C., Hurtt, G., Jackson, S. T., Johnson, C., King, G. A., … & Wyckoff, P. (1998). Reid’s Paradox of Rapid Plant Migration: Dispersal Theory and Interpretation of Paleoecological Records. BioScience, 48(1), 13-24. https://doi.org/10.2307/1313224

  5. ForestGEO. (n.d.). Understanding forest diversity and dynamics through a global network. Retrieved from https://www.forestgeo.si.edu/

  6. MDPI. (2020). Forest Stand Dynamics and Its Applications. Retrieved from https://www.mdpi.com/

  7. Shugart, H. H. (1984). A Theory of Forest Dynamics: The Ecological Implications of Forest Succession Models. Springer-Verlag. https://doi.org/10.1007/978-1-4612-5190-6

  8. U.S. Geological Survey. (2017). Improving Understanding of Forest and Carbon Dynamics. Retrieved from https://www.usgs.gov/