Introduction

To begin, we will conduct our analysis using the SSP2-4.5 climate scenario.

First, load the packages you will need. These are a good starting point:

# Load libraries
library(matilda)
library(tidyverse)

Reminder that loading matilda will automatically load hector as long as it is installed. Loading tidyverse will load all packages in the Tidyverse family, including ggplot2 for plotting and dplyr for cleaning, organizing, and summarizing data.

TASK: Create a Hector core using the SSP2-4.5 initialization file. If you have forgotten how to code this, revisit previous tutorials for help. Name your Hector core object –> ssp245_core. Once initialized, the core object will appear in you global environment panel.

Producing CO2 Fertilization Parameters

Here we will create a perturbed parameter set of BETA values. BETA is the parameter in Hector that stand for the CO2 fertilization effect on the global terrestrial carbon system.

In an attempt to produce and identical result, set a seed number before producing you BETA parameter sample. This should make it possible for our random sample to match each others:

# set your seed number
set.seed(25)

# produce random sample of parameter values -- produce 50 draws
parameter_vals <- generate_params(ssp245_core, draws = 50)

This produces a set of 50 random parameter values for all parameters. Remember that we only want to vary the CO2 fertilization parameter (BETA).

Create a new dataframe (df) that includes only the BETA parameter:

# create BETA parameter df
beta_vals <- 
  parameter_vals %>% 
  select(BETA)

Function notes

%>% - is a piping function, it passes the object on the left of %>% to the function on the right of the %>%. The keyboard shortcut for this function is CMD + SHIFT + M on Mac and CNTRL + SHIFT + M on PC.*

select - Keeps (or drops) variables in a dataframe. In the code above, we are using select to create a new df object that only contains one column BETA. All other columns we produced in parameter_vals are dropped.

if you are not using dplyr you can perform the same task with base R code:

beta_vals <- 
  parameter_vals["BETA"]

Make a note of the mean BETA value form your random sample:

print(mean(beta_vals$BETA))
## [1] 0.6286004

This value should be pretty close to the default BETA value in Hector, and will converge on the default value as we increase the number of random samples we draw from the prior distribution:

fetchvars(ssp245_core, dates = NA, vars = "beta")
##   scenario year variable value      units
## 1 SSP2-4.5   NA     beta  0.65 (unitless)

Running Hector-Matilda with CO2 Fertilization Uncertainty

Let’s run the model with our BETA parameter set.

TASK: Using the set of random BETA parameters, run Hector with Matilda for the date range: 1850:2100 and extract the following output variables: VEG_C SOIL_C NPP RH NBP CONCENTRATIONS_CO2 and GMST.

You should expect the model run to take about 25-30 seconds to run with each of the 50 BETA values:

# variables to save
vars = c(VEG_C(), SOIL_C(), NBP(), NPP(), RH(), CONCENTRATIONS_CO2(), GMST())

# run the model with beta uncertainty
mod_result <- iterate_model(core = ssp245_core,
                            params = beta_vals,
                            save_years = 1850:2100, 
                            save_vars = vars)

Lets check to make sure our model run delivered all the variables we wanted:

unique(mod_result$variable)
## [1] "veg_c"             "soil_c"            "NBP"              
## [4] "NPP"               "RH"                "CO2_concentration"
## [7] "gmst"

Function Notes

unique - a useful function that prints all the unique values in a specified column. In order for this to work you should use $ to specify which column in the dataframe (df) you want the unique values for. Here, we are asking R to print all the unique variable names in the variable column of the mod_result df object.

Preparing Data for Plotting

Lets plot this early result to see how each of the output variables are affect by variations in the CO2 fertilization factor.

Because we will want to visualize the model runs across the spectrum of BETA values and we want each model run to be color mapped to match the BETA value used to run the model, we first need to create a df with the BETA values added to the model results. dplyr will help with this:

# add new df with run_number column
beta_vals_rn <- 
  beta_vals %>% 
  mutate(run_number = row_number()) # creates a new column from existing information in beta_vals df

# create an object that joins BETA to model results based on run_number
mod_beta_plot <- 
  beta_vals_rn %>% 
  left_join(mod_result, by = "run_number")

Function notes

left_join - Adds columns from one df to another df. Here we are adding the BETA column from beta_vals_rn to mod_result based on the run_number key. The result will be a new df called mod_beta_plot that has a BETA column that connects eachBETA value to the model result it produced.

Introductory visualization

TASK: Produce a multi-paneled plot that shows how varying BETA affects the output variables of interest. Color each model run by the BETA value used to run the model.

A couple of things to remember:

  1. Create a line plot grouping the data by run_number.
  2. Make sure you use facet_wrap to create the panels, using scales to allow the y-axis scale to vary between panels.
  3. Add labels to the X and Y axis and a title.

Consult the Basic Matilda Run Rpubs page to help with code composition if needed.

Your figure should look something like this:

Provide a basic analysis of this plot

Based on the definitions of these variables below, answer the following questions:

Variable definitions

  • VEG_C: vegetation carbon stock. How much carbon is stored in living plant material.
  • SOIL_C: soil carbon stock. How much carbon is stored in the soil.
  • NPP: Net primary production. The rate at which plants add carbon to plant biomass. Interpret as carbon entering the land carbon pool from the atmosphere from plant growth via photosynthesis.
  • RH: Heterotrophic respiration. Carbon released from soil organic matter as microbes break down organic matter. Compared to NPP, RH is carbon leaving the land carbon pool and returning to the atmosphere through decomposition.
  • NBP: Net biome production. Net annual carbon movement from the atmosphere to the land. In Hector, a positive NBP means land is taking up carbon (carbon sink), negative NBP means the land is giving off carbon (carbon source). If NBP is declining, but still positive, the land is still taking up carbon, just more slowly.
  • CO2 Concentrations: Atmospheric CO2 concentration. How much CO2 is in the atmosphere.
  • GMST: Global mean surface temperature anomaly. Global temperature change relative to some baseline time period. For default Hector, this baseline year is 1745.

Analysis Questions

  1. How does higher CO2 fertilization factor (BETA) change how plants respond to CO2 (use only NPP as evidence)?

  2. Based on your answer to question 1, how does NPP affect the amount of carbon stock in vegetation and soil over time?

  3. How does increasing the terrestrial carbon pool (stocks of veg and soil carbon) affect respiration (RH)? Why?

  4. If the CO2 fertilization factor is high (> 0.7), how does that affect the CO2 concentrations overtime?

  5. Relate the strength of the CO2 fertilization factor to the potential global temperatures at the end of the century (2100).

  6. In ~2075 (near the end of the century), the movement of carbon into the terrestrial carbon pool being to slow (NBP). What might explain this slowing flux? To help consider: NBP = NPP - RH - land-use change emissions

What’s Next

Next, we will discuss the importance of normalizing certain variables to a reference period, such as 1850–1900. This makes our projections comparable to a preindustrial baseline, allowing us to measure how much a variable has changed relative to that period. The 1850–1900 baseline is commonly used in major global climate analyses. Another important baseline is 1995–2014, often used as a “recent past” reference period, which helps us describe changes relative to more recent climate conditions.

We will also discuss how to calculate summary statistics from the Hector output. These statistics will help us describe how uncertainty in the BETA parameter affects different Hector projections.