Science questions

  1. How much do tree stem GHG fluxes vary over the course of the year?
  2. …

Read data

# Data file from the TEMPEST repository
message("Reading data...")
## Reading data...
dat <- read_parquet("tempest_tree_ghg_fluxes.parquet")

message(nrow(dat), " rows and ", ncol(dat), " columns")
## 3361 rows and 12 columns

Soil data

# Read soil temperature and moisture data
files <- list.files(path = "L2_data/", pattern = "parquet$", full.names = TRUE)
lapply(files, read_parquet) %>% 
  bind_rows() %>% 
  # compute a 'Date' column...
  mutate(Date = as.Date(TIMESTAMP)) %>% 
  # ...and average variables for each date
  group_by(Site, Plot, Date, research_name) %>% 
  summarise(Value = mean(Value), .groups = "drop") %>% 
  # reshape 
  pivot_wider(names_from = "research_name", values_from = "Value", names_repair = "universal") ->
  soildat
## New names:
## • `soil-temp-5cm` -> `soil.temp.5cm`
## • `soil-vwc-5cm` -> `soil.vwc.5cm`
# TODO: plot soil data

Quick summary and visualization of data

glimpse(dat)
## Rows: 3,361
## Columns: 12
## Groups: Year, Date, Plot, Timepoint, Species, ID [3,361]
## $ Year                  <dbl> 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, …
## $ Date                  <date> 2021-10-08, 2021-10-08, 2021-10-08, 2021-10-08,…
## $ Plot                  <chr> "Control", "Control", "Control", "Control", "Con…
## $ Timepoint             <chr> "(none)", "(none)", "(none)", "(none)", "(none)"…
## $ Species               <chr> "Tulip Poplar", "Red Maple", "Red Maple", "Beech…
## $ ID                    <chr> "C1", "C10", "C12", "C14", "C16", "C17", "C2", "…
## $ CO2_lin_flux.estimate <dbl> 0.42493599, 0.22701875, 0.08174378, 0.29282769, …
## $ CO2_lin_r.squared     <dbl> 0.9980241, 0.9926840, 0.9514107, 0.9952041, 0.84…
## $ CO2_rob_flux.estimate <dbl> 0.42512176, 0.22678729, 0.08119638, 0.29299455, …
## $ CH4_lin_flux.estimate <dbl> 0.1516792079, -0.0052212805, 0.1701407219, 0.000…
## $ CH4_lin_r.squared     <dbl> 0.999189922, 0.602961702, 0.999108777, 0.0170388…
## $ CH4_rob_flux.estimate <dbl> 0.1517037352, -0.0052248314, 0.1698198573, 0.000…
# 1. Print summary
summary(dat)
##       Year           Date                Plot            Timepoint        
##  Min.   :2021   Min.   :2021-10-08   Length:3361        Length:3361       
##  1st Qu.:2022   1st Qu.:2022-07-07   Class :character   Class :character  
##  Median :2023   Median :2023-06-06   Mode  :character   Mode  :character  
##  Mean   :2023   Mean   :2023-06-28                                        
##  3rd Qu.:2024   3rd Qu.:2024-06-12                                        
##  Max.   :2025   Max.   :2025-09-12                                        
##    Species               ID            CO2_lin_flux.estimate
##  Length:3361        Length:3361        Min.   :-0.08419     
##  Class :character   Class :character   1st Qu.: 0.11837     
##  Mode  :character   Mode  :character   Median : 0.30497     
##                                        Mean   : 0.38523     
##                                        3rd Qu.: 0.58891     
##                                        Max.   : 3.49284     
##  CO2_lin_r.squared   CO2_rob_flux.estimate CH4_lin_flux.estimate
##  Min.   :5.367e-05   Min.   :-0.09171      Min.   :-0.0970870   
##  1st Qu.:9.682e-01   1st Qu.: 0.11806      1st Qu.: 0.0007698   
##  Median :9.937e-01   Median : 0.30576      Median : 0.0049256   
##  Mean   :9.508e-01   Mean   : 0.38510      Mean   : 0.0259362   
##  3rd Qu.:9.985e-01   3rd Qu.: 0.58729      3rd Qu.: 0.0168636   
##  Max.   :1.000e+00   Max.   : 3.48976      Max.   : 0.9723704   
##  CH4_lin_r.squared   CH4_rob_flux.estimate
##  Min.   :1.300e-07   Min.   :-0.0892519   
##  1st Qu.:1.048e-01   1st Qu.: 0.0007614   
##  Median :5.729e-01   Median : 0.0049023   
##  Mean   :5.274e-01   Mean   : 0.0258854   
##  3rd Qu.:9.364e-01   3rd Qu.: 0.0167526   
##  Max.   :1.000e+00   Max.   : 0.9740040
# 2. ggplot(dat, ...)

Join with species data

# 3. Drop existing species column from dat
# 4. Read in TEMPEST_Trees_Instruments.csv

tree_inventory <- read.csv("TEMPEST_Trees_Instruments.csv") %>% 
  select(Sapflux_ID, spp) # don't need these

#
dat %>%
  left_join(tree_inventory,
            # in dat the column is "ID", while in the tree inventory it's
            # "Sapflux_ID", so we tell left_join that these two columns 
            # correspond to each other 
            by = c("ID" = "Sapflux_ID"), 
            # this is optional but good 'defensive programming': we specify
            # that we expect a many-to-one relationship (i.e., there shouldn't
            # be multiple matches for a given sapflux tree!)
            relationship = "many-to-one") ->
        tree_dat

# To help readability, generate the common name from the species code
tree_dat %>%
  mutate(Species = case_when(
    spp == "ACRU" ~ "Red maple",
    spp == "LITU" ~ "Tulip poplar",
    spp == "FAGR" ~ "American beech"
  )) -> 
  tree_dat

# Make a nice summary table to check things
tree_dat %>% 
  # for each plot and species...
  group_by(Plot, Species) %>% 
  # ...count the number of distinct IDs...
  summarise(value = n_distinct(ID), .groups = "drop") %>% 
  # ...and then reshape so that plot names run across the top
  pivot_wider(names_from = "Plot") %>% 
  knitr::kable()
Species Control Freshwater Seawater
American beech 6 6 6
Red maple 8 6 6
Tulip poplar 6 6 6

Plot annual cycle - raw data

# Filter out the TEMPEST flood measurements
tree_dat %>% 
  filter(Timepoint %in% c("(none)", "Pre-Treatment")) ->
  tree_dat_ac

message("Annual cycle dataset has ", nrow(tree_dat_ac), " rows and ",
        ncol(tree_dat_ac), " columns")
## Annual cycle dataset has 2184 rows and 13 columns
# EVAN:
# 1. What data frame do we want to use to plot?
# 2. This plots each tree as a separate color. How useful is this? What are alternatives?
ggplot(data = tree_dat, 
       mapping = aes(x = Date, y = CO2_lin_flux.estimate, color = Species))+
  geom_point() + # Add a layer of points
  labs(title = "CO2 flux by species by date")  # Add labels and a title

ggplot(data = tree_dat, 
       mapping = aes(x = Date, y = CH4_lin_flux.estimate, color = Species))+
  geom_point() + # Add a layer of points
  labs(title = "CH4 flux by species by date")  # Add labels and a title

# Summarize by …?

The trees are experimental units but what we’re actually interested in is…

# EVAN: how do we want to summarize and then plot the annual cycle data?

library(dplyr)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:arrow':
## 
##     duration
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
#calculate average CO2 flux per year per species

tree_dat %>%
  mutate(Year = year(Date)) %>%  # Extract the year from the Date column
  group_by(Year, Species) %>%  # Group by Year and Species
  summarise(avg_CO2_flux = mean(CO2_lin_flux.estimate, na.rm = TRUE), .groups = "drop") -> 
  avg_CO2flux_per_year_species

# View the resulting table
print(avg_CO2flux_per_year_species)
## # A tibble: 15 × 3
##     Year Species        avg_CO2_flux
##    <dbl> <chr>                 <dbl>
##  1  2021 American beech       0.0893
##  2  2021 Red maple            0.142 
##  3  2021 Tulip poplar         0.148 
##  4  2022 American beech       0.383 
##  5  2022 Red maple            0.310 
##  6  2022 Tulip poplar         0.545 
##  7  2023 American beech       0.417 
##  8  2023 Red maple            0.303 
##  9  2023 Tulip poplar         0.505 
## 10  2024 American beech       0.386 
## 11  2024 Red maple            0.269 
## 12  2024 Tulip poplar         0.524 
## 13  2025 American beech       0.281 
## 14  2025 Red maple            0.294 
## 15  2025 Tulip poplar         0.419
# Plotting the average CO2 flux by year and species

ggplot(avg_CO2flux_per_year_species, aes(x = Year, y = avg_CO2_flux, color = Species)) +
  geom_line(size = 1) +  # Add lines to connect points
  geom_point(size = 3) +  # Add points at each data point
  labs(
    title = "Average Annual CO2 Flux by Species", 
    x = "Year", 
    y = "Average CO2 Flux (g CO2 / m² / day)",
    color = "Species"
  ) +
  theme_minimal() +  # Clean minimal theme
  theme(
    legend.position = "top",  # Position the legend at the top
    axis.text.x = element_text(angle = 45, hjust = 1)  # Angle x-axis labels for better readability
  )
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#calculate average CH4 flux per year per species

tree_dat %>%
  mutate(Year = year(Date)) %>%  # Extract the year from the Date column
  group_by(Year, Species) %>%  # Group by Year and Species
  summarise(avg_CH4_flux = mean(CH4_lin_flux.estimate, na.rm = TRUE), .groups = "drop") -> 
  avg_CH4flux_per_year_species

# View the resulting table
print(avg_CH4flux_per_year_species)
## # A tibble: 15 × 3
##     Year Species        avg_CH4_flux
##    <dbl> <chr>                 <dbl>
##  1  2021 American beech      0.00259
##  2  2021 Red maple           0.0465 
##  3  2021 Tulip poplar        0.0119 
##  4  2022 American beech      0.0173 
##  5  2022 Red maple           0.0554 
##  6  2022 Tulip poplar        0.0141 
##  7  2023 American beech      0.0147 
##  8  2023 Red maple           0.0469 
##  9  2023 Tulip poplar        0.0165 
## 10  2024 American beech      0.0148 
## 11  2024 Red maple           0.0412 
## 12  2024 Tulip poplar        0.0133 
## 13  2025 American beech      0.0148 
## 14  2025 Red maple           0.0471 
## 15  2025 Tulip poplar        0.0130
# Plotting the average CH4 flux by year and species
ggplot(avg_CH4flux_per_year_species, aes(x = Year, y = avg_CH4_flux, color = Species)) +
  geom_line(size = 1) +  # Add lines to connect points
  geom_point(size = 3) +  # Add points at each data point
  labs(
    title = "Average Annual CH4 Flux by Species", 
    x = "Year", 
    y = "Average CO2 Flux (g CH4 / m² / day)",
    color = "Species"
  ) +
  theme_minimal() +  # Clean minimal theme
  theme(
    legend.position = "top",  # Position the legend at the top
    axis.text.x = element_text(angle = 45, hjust = 1)  # Angle x-axis labels for better readability
  )

# Change “Seawater” to “Saltwater”


tree_dat <- tree_dat %>%
  mutate(Plot = recode(Plot, "Seawater" = "Saltwater"))

Looking at average CO2 and CH4 by date instead of year

#calculate average CO2 flux per date per species

# Summarizing by specific Date and Species
avg_CO2flux_per_day_species <- tree_dat %>%
  group_by(Date, Species) %>%  # Group by full Date and Species
  summarise(avg_CO2_flux = mean(CO2_lin_flux.estimate, na.rm = TRUE), .groups = "drop")

# Ensure the Date column is in Date format (if needed)
avg_CO2flux_per_day_species$Date <- as.Date(avg_CO2flux_per_day_species$Date)


# Plotting the average flux by Date and Species
ggplot(avg_CO2flux_per_day_species, aes(x = Date, y = avg_CO2_flux, color = Species)) +
  geom_line(size = 1) +  # Add lines to connect points
  geom_point(size = 3) +  # Add points at each data point
  labs(
    title = "Average CO2 Flux by Date and Species", 
    x = "Date", 
    y = "Average CO2 Flux (g CO2 / m² / day)",
    color = "Species"
  ) +
  theme_minimal() +  # Clean minimal theme
  theme(
    legend.position = "top",  # Position the legend at the top
    axis.text.x = element_text(angle = 45, hjust = 1)  # Rotate x-axis labels for readability
  )


  #calculate average CH4 flux per date per species

# Summarizing by specific Date and Species
avg_CH4flux_per_day_species <- tree_dat %>%
  group_by(Date, Species) %>%  # Group by full Date and Species
  summarise(avg_CH4_flux = mean(CH4_lin_flux.estimate, na.rm = TRUE), .groups = "drop")

# Ensure the Date column is in Date format (if needed)
avg_CH4flux_per_day_species$Date <- as.Date(avg_CH4flux_per_day_species$Date)


# Plotting the average CH4 flux by Date and Species
ggplot(avg_CH4flux_per_day_species, aes(x = Date, y = avg_CH4_flux, color = Species)) +
  geom_line(size = 1) +  # Add lines to connect points
  geom_point(size = 3) +  # Add points at each data point
  labs(
    title = "Average CH4 Flux by Date and Species", 
    x = "Date", 
    y = "Average CH4 Flux (g CO2 / m² / day)",
    color = "Species"
  ) +
  theme_minimal() +  # Clean minimal theme
  theme(
    legend.position = "top",  # Position the legend at the top
    axis.text.x = element_text(angle = 45, hjust = 1)  # Rotate x-axis labels for readability
  )

Remove s7

# Specify the tree ID you want to remove (S7)
tree_id_to_remove <- "S7"

# Create a new dataset excluding the specified tree ID (S7)
tree_dat_no_S7 <- tree_dat %>%
  filter(ID != tree_id_to_remove)  # Exclude the tree with ID = "S7"

# Check column names to ensure 'ID' is correct
colnames(tree_dat)
##  [1] "Year"                  "Date"                  "Plot"                 
##  [4] "Timepoint"             "Species"               "ID"                   
##  [7] "CO2_lin_flux.estimate" "CO2_lin_r.squared"     "CO2_rob_flux.estimate"
## [10] "CH4_lin_flux.estimate" "CH4_lin_r.squared"     "CH4_rob_flux.estimate"
## [13] "spp"
# Specify the tree ID you want to remove (S7)
tree_id_to_remove <- "S7"

# Calculate the average CH4 flux per date and species excluding tree "S7"
avg_CH4flux_no_S7 <- tree_dat_no_S7 %>%
  group_by(Date, Species) %>%  # Group by Date and Species
  summarise(avg_CH4_flux = mean(CH4_lin_flux.estimate, na.rm = TRUE), .groups = "drop")

# Ensure the Date column is in Date format (if needed)
avg_CH4flux_no_S7$Date <- as.Date(avg_CH4flux_no_S7$Date)

# Plotting the average CH4 flux by Date and Species (excluding tree "S7")
ggplot(avg_CH4flux_no_S7, aes(x = Date, y = avg_CH4_flux, color = Species)) +
  geom_line(size = 1) +  # Add lines to connect points
  geom_point(size = 3) +  # Add points at each data point
  labs(
    title = "Average CH4 Flux by Year and Species (Excluding Tree S7)", 
    x = "Date", 
    y = "Average CH4 Flux (g CH4 / m² / day)",
    color = "Species"
  ) +
  theme_minimal() +  # Clean minimal theme
  theme(
    legend.position = "top",  # Position the legend at the top
    axis.text.x = element_text(angle = 45, hjust = 1)  # Rotate x-axis labels for readability
  )

Add error bars and remove 2021 for CH4 data

avg_CH4flux_stats_filtered <- tree_dat %>%
  filter(ID != "S7") %>%             # remove tree S7
  mutate(Year = year(Date)) %>%      # extract year
  filter(Year != 2021) %>%           # exclude 2021
  group_by(Year, Species) %>%        # group by Year for annual average
  summarise(
    avg_CH4_flux = mean(CH4_lin_flux.estimate, na.rm = TRUE),
    sd_CH4_flux  = sd(CH4_lin_flux.estimate, na.rm = TRUE),
    .groups = "drop"
  )

# Plot using the filtered dataset
ggplot(avg_CH4flux_stats_filtered, aes(x = Year, y = avg_CH4_flux, color = Species)) +
  geom_line(size = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(
    ymin = avg_CH4_flux - sd_CH4_flux,
    ymax = avg_CH4_flux + sd_CH4_flux
  ), width = 0.2) +
  labs(
    title = "Average Annual CH4 Flux by Species (Excluding Tree S7 and 2021)",
    x = "Year",
    y = "Average CH4 Flux (g CH4 / m² / day)",
    color = "Species"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

Incorporate Plots (SW, FW, Control) into figures

avg_CH4flux_stats_filtered <- tree_dat %>%
  filter(ID != "S7") %>%
  mutate(Year = lubridate::year(Date)) %>%
  filter(Year != 2021) %>%
  group_by(Year, Species, Plot) %>%   # <- include Plot here
  summarise(
    avg_CH4_flux = mean(CH4_lin_flux.estimate, na.rm = TRUE),
    sd_CH4_flux = sd(CH4_lin_flux.estimate, na.rm = TRUE),
    .groups = "drop"
  )
  
  ggplot(avg_CH4flux_stats_filtered,
       aes(x = Year, y = avg_CH4_flux, color = Species)) +
  geom_line() +
  geom_point() +
  geom_errorbar(aes(
    ymin = avg_CH4_flux - sd_CH4_flux,
    ymax = avg_CH4_flux + sd_CH4_flux
  ), width = 0.2) +
  facet_wrap(~ Plot) +   # <- split panels by plot
  theme_minimal()

Incorporate soildat into CH4 data

tree_soil_dat <- tree_dat %>%
  mutate(
    Date = as.Date(Date),
    Plot = recode(Plot, "C" = "Control", "F" = "Freshwater", "S" = "Seawater")
  ) %>%
  left_join(
    soildat %>%
      mutate(Date = as.Date(Date)) %>%
      group_by(Plot, Date) %>%
      summarise(soil.temp.5cm = mean(soil.temp.5cm, na.rm = TRUE), .groups = "drop"),
    by = c("Plot", "Date"),
    relationship = "many-to-one"
  )
  

  
  # Ensure Date class matches
tree_dat$Date <- as.Date(tree_dat$Date)
soildat$Date <- as.Date(soildat$Date)

# Optional: recode plots to full names
tree_dat <- tree_dat %>%
  mutate(Plot = recode(Plot, "C" = "Control", "F" = "Freshwater", "S" = "Seawater"))

soildat <- soildat %>%
  mutate(Plot = recode(Plot, "C" = "Control", "F" = "Freshwater", "S" = "Seawater"))

# Join soil temperature
tree_soil_dat <- tree_dat %>%
  left_join(
    soildat %>%
      group_by(Plot, Date) %>%
      summarise(soil.temp.5cm = mean(soil.temp.5cm, na.rm = TRUE), .groups = "drop"),
    by = c("Plot", "Date")
  )

# Check for unmatched rows
summary(tree_soil_dat$soil.temp.5cm)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.578  12.309  18.594  16.279  20.800  24.753     679
  # Prepare data for plotting
tree_soil_plot <- tree_soil_dat %>%
  filter(
    !is.na(soil.temp.5cm),
    !is.nan(soil.temp.5cm),
    !is.na(CH4_lin_flux.estimate),
    !is.nan(CH4_lin_flux.estimate)
    )
    
ggplot(tree_soil_plot,
       aes(x = soil.temp.5cm, y = CH4_lin_flux.estimate, color = Species)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  facet_wrap(~ Plot) +
  theme_minimal() +
  labs(
    title = "CH4 Flux vs Soil Temperature by Plot",
    x = "Soil Temperature (5 cm, °C)",
    y = "CH4 Flux (g CH4 / m² / day)",
    color = "Species"
  )
## `geom_smooth()` using formula = 'y ~ x'

Remove S7

# Prepare the dataset excluding tree S7
tree_soil_plot_no_S7 <- tree_dat %>%
  filter(ID != "S7") %>%                              # remove tree S7
  mutate(
    Date = as.Date(Date),
    Plot = recode(Plot, "C" = "Control", "F" = "Freshwater", "S" = "Seawater")
  ) %>%
  left_join(
    soildat %>%
      mutate(Date = as.Date(Date)) %>%
      group_by(Plot, Date) %>%
      summarise(soil.temp.5cm = mean(soil.temp.5cm, na.rm = TRUE), .groups = "drop"),
    by = c("Plot", "Date"),
    relationship = "many-to-one"
  ) %>%
  filter(
    !is.na(soil.temp.5cm),
    !is.nan(soil.temp.5cm),
    !is.na(CH4_lin_flux.estimate),
    !is.nan(CH4_lin_flux.estimate)
  )

# Plot using the new dataset

ggplot(tree_soil_plot_no_S7,
       aes(x = soil.temp.5cm, y = CH4_lin_flux.estimate, color = Species)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  facet_wrap(~ Plot) +
  theme_minimal() +
  labs(
    x = "Soil Temperature (5 cm, °C)",
    y = "CH4 Flux (g CH4 / m² / day)",
    title = "CH4 Flux vs Soil Temperature by Plot (Excluding Tree S7)"
  )
## `geom_smooth()` using formula = 'y ~ x'

Doing the same for CO2

# Prepare the dataset including S7 for CO2
tree_soil_plot_CO2 <- tree_dat %>%
  mutate(
    Date = as.Date(Date),
    Plot = recode(Plot, "C" = "Control", "F" = "Freshwater", "S" = "Seawater")
  ) %>%
  left_join(
    soildat %>%
      mutate(Date = as.Date(Date)) %>%
      group_by(Plot, Date) %>%
      summarise(soil.temp.5cm = mean(soil.temp.5cm, na.rm = TRUE), .groups = "drop"),
    by = c("Plot", "Date"),
    relationship = "many-to-one"
  ) %>%
  filter(
    !is.na(soil.temp.5cm),
    !is.nan(soil.temp.5cm),
    !is.na(CO2_lin_flux.estimate),
    !is.nan(CO2_lin_flux.estimate)
  )

# Plot using the new dataset
ggplot(tree_soil_plot_CO2,
       aes(x = soil.temp.5cm, y = CO2_lin_flux.estimate, color = Species)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  facet_wrap(~ Plot) +
  theme_minimal() +
  labs(
    x = "Soil Temperature (5 cm, °C)",
    y = "CO2 Flux (g CO2 / m² / day)",
    title = "CO2 Flux vs Soil Temperature by Plot (Including Tree S7)"
  )
## `geom_smooth()` using formula = 'y ~ x'