High-resolution time-series of VOC fluxes from leaf litter

In an incubation experiment using tropical leaf litter at the Biosphere II, four chambers of litter and one blank were measured continuously for VOC fluxes, where the background air was from the Biosphere II. Air temperature and relative humidity were measured in each chamber with iButtons. 304 VOCs were identified by their mass-to-charge ratio (mz), with recognized compounds identified from the GLOVOCS library and Giovanni.

In this repository, we process the raw data through a pipeline that involves calibration -> diurnal detrending -> model fitting -> parameter extraction. The VOC responses were categorized into 9 classes based on their response to initial and secondary experimental wetting. The original and detrended timeseries, fitted models, and parameter comparison can for all 304 VOCs can be interactively explored in a Shiny app.

Data notes

Initial chamber flux measurements occurred at 2020-02-07 18:00:00. Some instrument error occurred on day 11, so data was trimmed to end on 2020-02-17.

15 mL of double distilled water added to four chambers and blank at the following times:

• chamber 1: 2020-02-14 15:11:00
• chamber 2: 2020-02-14 15:26:00
• chamber 3: 2020-02-14 15:41:00
• chamber 4: 2020-02-14 14:56:00
• blank: 2020-02-14 14:51:00

Fluxes were calculated from the difference in relative concentrations between blanks and chambers. These differences were multiplied by the molar flow (from ideal gas constant, pressure, temperature, and flow rate divided by chamber surface area) to obtain fluxes in units of micro-moles per meter squared per second ($\mu mol m^{-2} s^{-1}$). Positive fluxes indicate compound emission and negative fluxes indicate compound uptake.

Methods

Upon data exploration, some VOC timeseries exhibited a diurnal fluctuation seemingly unrelated to the initial incubation or the litter rewetting event. Furthermore, no consistent differences among the four replicate chambers were found. Thus, each VOC timeseries were summarized across the four chambers, divided into two temporal phases using the mean pulse addition time of 2020-02-14 15:00:00, then decomposed with STL (seasonal and trend decomposition using loess) using a 24 hour window. All subsequent analysis was performed withe decomposed mean fluxes (trend and remainder), after the diurnal pattern was removed (see 'Raw timeseries' tab of the Shiny app).

Visualization of the decomposed mean fluxes revealed three major categories of flux responses to either initial incubation or litter rewetting:

1. no response

2. exponential decline

3. delayed pulse and decline

Thus, we created an algorithm to fit exponential and Ricker models to both phases and developed criteria to select the model of best fit. Based on the results of the algorithm, each VOC was assigned to one of nine categories. Finally, common parameters were calculated for exponential and Ricker fits to allow for comparison across all 304 VOCs and 2 phases.

Exponential model

For each observation $i$, the decomposed flux means were given a normal likelihood:

$flux_i \sim Normal(\mu_i, \sigma_{p}^2)$

where $\mu$ is the expected value of each flux mean, $\sigma^2$ is the measurement error variance, and $p$ represents the phase associated with each observation $i$. The exponential model was defined as:

$\mu_i = c_p + a_p \cdot e^{-b_p \cdot t_i}$

where the parameters $a$, $b$, and $c$ are independently estimated for each phase and $t$ is the time in hours since the beginning of the phase associated with each observation $i$.

Parameters $a$ and $c$ were given relatively non-informative normal priors, while $b$ as assigned $Gamma(0.1, 0.1)$ prior. Measurement error precisions ($1/\sigma_p^2$) were assigned conjugate, relatively non-informative $Gamma(0.1, 0.1)$ priors. Posterior predictive loss was computed as the sum of squared differences between the observed and replicated data.

Ricker model

For each observation $i$, the decomposed flux means were given a normal likelihood:

$flux_i \sim Normal(\mu_i, \sigma_{p}^2)$

where $\mu$ is the expected value of each flux mean, $\sigma^2$ is the measurement error variance, and $p$ represents the phase associated with each observation $i$. The Ricker model was defined as:

$\mu_i = c_p + a_p \cdot t_i \cdot e^{-b_p \cdot t_i}$

where the parameters $a$, $b$, and $c$ are independently estimated for each phase and $t$ is the time in hours since the beginning of the phase associated with each observation $i$.

Parameter $c$ was given a relatively non-informative normal prior, while $a$ and $b$ were assigned $Gamma(0.1, 0.1)$ priors. Parameters $a$ and $c$ were restricted to the same direction (positive or negative). Measurement error precisions ($1/\sigma_p^2$) were assigned conjugate, relatively non-informative $Gamma(0.1, 0.1)$ priors. Posterior predictive loss was computed as the sum of squared differences between the observed and replicated data.

Algorithm implementation

The above models were implemented in JAGS 4.3.0 via the R package 'rjags'and R 4.1.1. Custom functions were written to run each model in three parallel Markov chain Monte Carlo with randomly generated initial values. After a burn-in period of 10000 iterations, chains were sampled for 15000 iterations thinned by 5, yielding a total of 3000 relatively independent samples. If the chains had not converged ($\hat{R} > 1.3$), the model was restarted with the final iteration of the previous model run until convergence was achieved. Monitored parameters and replicated data were saved for further evaluation.

To assess best fit for each VOC and phase, model statistics were calculated for both exponential and Ricker model outputs. Statistics included $\hat{R}$, posterior predictive loss, and $r^2$ between observed and replicated. The best fit was determined in two steps:

1. For each VOC (n = 304) and phase (n = 2), the model (exponential or Ricker) with the lowest posterior predictive loss was selected if $\hat{R} < 1.35$ and $r^2 > 0.4$. Otherwise, no model was assigned ("none").

2. For each VOC (n = 304), the best fits for phase I and phase II were combined into one of nine categories (e.g., "exponential, Ricker", "none, exponential", "none, none")

In each category, VOCs can be selected to visualize the observed and modeled timeseries and to compare the raw parameters between the phases (see 'Classification' tab of the Shiny app).

Comparing derived parameters

Additional parameters were derived to be comparable between exponential and Ricker model fits.

1. Maximum or minimum flux rate

2. Time (hours) until 20% of maximum or minimum was attained

3. Slope of exponential or Ricker decline

One or more categories can be selected to compare the derived parameters for each VOC and phase (see 'Parameters' tab of the Shiny app).

File structure

• app/ contains .Rdata files and Shiny script

  • app.R creates app with 3 tabs to show raw and decomposed data, VOCs classified into 9 categories with model fit, and comparison of extracted parameters of time of peak flux, initial flux slope, and maximum flux

• data_clean/ contains cleaned and/or summarized csvs produced by code in the scripts/ folder

  • env_chamber_hourly.csv temperature, RH, and VPD for each chamber rounded to hourly
  • env_hourly.csv mean temperature, RH, and VPD across chambers rounded to hourly
  • env_long.csv mean temperature, RH, and VPD across chambers rounded to nearest 5 minutes
  • flux_chamber_hourly.csv flux of all VOCs and environmental variables by chamber, rounded to hourly
  • flux_hourly.csv mean flux of all VOCs and environmental variables across chambers, rounded to hourly

• data_raw/ consists of original files shared by Lia Crocker

  • iButton_continuous.xlsx contains air temperature and relative humidity for each chamber as separate tabs

  • VOC_masses_ID.xlsx annotates the molecular weights by compount names and groups where available

  • chambers/

    • Chamber_3_Blank_VOC_Concentration.csv raw VOC concentrations, reported as averages and standard deviation for each compound
    • Chamber_X_Litter_VOC_Flux.csv VOC time series of average flux with datetime, date, and days since start of experiment.

  • Litter_VOC_daily_avg_flux.csv is average across days and chamber

• flux/ contains the raw data and calibrations to produce fluxes. See /flux/README.md for further documentation

• models/ contains trials of different modelling approaches

  • 01-model-exploration/ sandbox for testing a mixture model between diurnal variation and pulse response
    • control.R
  • 02-model-VPD/ attempt to model fluxes either a linear function of VPD or a pulse response

    • 01-run-model.R

    • 02-test-function.R

    • 03-categorize.R

    • fit_models.R

    • VPD_pulse_response.jags

    • VPD_response.jags

  • 03-model-Ricker/ attempt to model fluxes as either a linear function of VPD, a Ricker response, or a combination of both

    • 01-run-model.R

    • 02-test-function.R

    • 03-run-function.R

    • 04-classify.R

    • fit_models.R

    • Ricker_only.jags

    • VPD_only.jags

    • VPD_Ricker.jags

  • 04-model-decomposed/ model fluxes by pulse as either exponential decay or Ricker response

    • 01-test-function.R sandbox to test understanding of the Ricker parameters

    • 02-process-output.R combines and applies categorization criteria to output of parallel/singleIteration.R

    • 03-checking-misclassification.R sand box to verify why particular VOCs might be misclassified

    • 04-plotting-params.R gather model parameters to summarize and plot

    • 05-data-for-app.R prepare model data needed for app/

    • 06-create-functions.R sandbox to test how to find % of max in both exponential and Ricker models

    • exponential.jags

    • Ricker_only.jags

    • fit_models_exp.R

    • fit_models_rick.R

    • parallel/

      • singleIteration.R script to run for a single VOC, can be locally parallelized if run from terminal
  • 05-model-decomposed-chamber/ currently empty except for data inputs. Not considered a high priority because running hierarchical Ricker and exponential models will likely involve numerous convergence issues

• plots/ contains overall plot and four sub-folders containing collections of plots

  • fig_env.png mean and chamber-specific temp, RH, and VPD during experiment

  • heatmaps/ contains 5 subsets of VOC plotted as the temporal signature of fluxes

  • STL_decomp

  • STL_deomp_chamber

  • VOC contains 304 plots of flux and VPD over time\

• scripts/ consists of R scripts for cleaning, preparing, and plotting data

  • 01-organize-env.R reads in the iButton data, calculates VPD, and outputs with rounded 5-minute time stamps across chambers, rounded hourly time stamps by chamber and across chambers
  • 02-organize-flux.R reads in chamber specific raw fluxes, joins with hourly environmental data, and outputs chamber specific and averaged across-chamber fluxes
  • 02b-check-tz.R checks the timing of the moisture additions to verify an error in the original time zones
  • 03-plot-heatmap.R tests different clustering approaches to visualize the temporal fingerprint of fluxes across all VOCs, output in plots/heatmaps
  • 04-plot-env-ts.R plots timeseries and bivariate plots of mean and chamber-specific and VPD
  • 05-test-ts-models.R follows the fpp3 package and accompanying textbook for different timeseries approaches
  • 06-decompose-ts.R uses STL decomposition to remove the diurnal trend of average chamber and chamber-specific fluxes for each VOC, outputs are saved scripts/04-model-decomposed/ and app/
  • 07-data-for-app.R curate average chamber and chamber-specific fluxes for app/
  • pheatmap.R was sourced from Lia Crocker to use the test the pheatmap package

• source/ functions used elsewhere

  • round.POSIXct.R custom rounding function for POSIX objects

renv

This project uses renv to capture R package dependencies for the Shiny app in app.R. This helps to ensure that the deployed Shiny app is able to run and will use the exact versions of all packages recorded in renv.lock by renv::snapshot(). When you open this project you will be prompted to run renv::restore() to install necessary packages. You may need to "downgrade" your R version to get package installation to work. For managing multiple installations of R, check out rig. The .renvignore file is there to only make renv track dependencies in app/app.R, so if you'd like to run code in other files, you may need to install additional packages not installed by renv::restore().

---

Developed in collaboration with the University of Arizona CCT Data Science team