chlorophyllCycling

Code for the modeling and analyses behind:

A trigger may not be necessary to cause senescence in deciduous broadleaf forests Kathryn I. Wheeler and Michael C. Dietze

Preprint: biorxiv.org/content/10.1101/2023.06.07.544057v1

Model overview

Daily canopy greenness is modeled as a Bayesian state-space process: greenness at day t depends on greenness at day t-1 plus temperature- and photoperiod-driven chlorophyll synthesis, minus a constant proportional breakdown rate. The chlorophyll cycling (CC) model is fit using JAGS (via rjags/runjags) and compared against a historical-average ("climatology") null model.

Note: in the code, the fitted parameters b0, b3, and b4 correspond to b0, b1, and -1*b2 (i.e., b4 = -b2) in the manuscript's Equation 1.

Manuscript Fig. 1. Parts of the yearly greenness curve. Gray shading indicates data excluded as spring (1 January–30 June). Full green-down (GD) encompasses the summer gradual GD, senescence, and abscission; this study focuses on predicting greenness and the amount of inflection around the SOS.

Manuscript Fig. 2. Example time series of modeled 95% credible intervals (blue shading) versus calibration (black circles) and validation (red triangle) data at (a) Howland, ME and (b) Alligator River, NC, calibrated to the full autumn; (c,d) Lac Clair, QC, calibrated to 43 vs. 65 days/year; and (e) Willow Creek, WI, calibrated to 76 days/year. The model reproduces withheld greenness, including senescence, without being trained on any senescence-period data.

Data

The scripts expect the following input files (paths are set in generalVariables.R, see below):

• allPhenocamDBsitesComplete.csv — site metadata, with columns siteName (PhenoCam site name), Lat, Long, URL (PhenoCam data archive URL for the site/ROI), PFT (plant functional type, i.e. DB), TZ (offset from GMT, e.g. -5 for ET), startDate (first date of a full year of data, e.g. "2013-01-01"), URL2/URL3 (additional PhenoCam URLs if the camera/ROI changed), endDate
• allPhenocamDBsitesSummary.csv — summary site table used to generate Table S1
• phenocamSitesBadYears.csv — site-years visually assessed as poor quality due to large gaps of missing data

PhenoCam greenness data and ERA5 reanalysis meteorology are openly available at doi.org/10.1038/s41597-019-0229-9 and doi.org/10.1002/qj.3803, respectively.

Setup

1. Install the required R packages: rjags, runjags, doParallel, ecoforecastR, RColorBrewer, scoringRules, scales, phenopix, zoo, suncalc, ncdf4, reticulate, randomForest, mgcv, tidyverse, caTools. (ecoforecastR and phenopix are installed from GitHub rather than CRAN.)
2. Update the file paths in generalVariables.R (e.g. headFilePath, ERA5dataFolder) to match your machine. This file also defines shared variables (site lists, output folders, MCMC settings) and is sourced by nearly every other script.
3. downloadERA5Calibration.R calls the Copernicus Climate Data Store API via reticulate; you need CDS API credentials configured to run it.

Scripts

Prepping data

• downloadERA5Calibration.R — downloads ERA5 reanalysis data for the selected sites (requires CDS API access, see Setup)
• downloadPhenocam.R — downloads PhenoCam GCC data based on the URLs in the site data file
• load_ERA5.R — functions to load the downloaded ERA5 met data from the saved NetCDF files; sourced by createDataObjects.R
• createElmoreFitsForRescaling.R — fits the phenopix ElmoreFit curve to each site and saves the rescaling info as paste0(dataDirectory,siteName,"_phenopixOutputs.RData")
• createDataObjects.R — combines PhenoCam and ERA5 data (including photoperiod, calculated with suncalc) into a comprehensive data object per site, saved as paste0(dataDirectory,siteName,"_dataFinal.RData")

Calculating PhenoCam transition dates

• estimatePhenoCamTransitions_changePoint.R — fits a changepoint Bayesian model to PhenoCam data to estimate start-of-senescence (SOS) transition dates; saves fits to paste0(transitionEstimateOutputsFolder,siteName,"_",yrName,"_PhenoCam_changePointCurve_varBurn.RData")
• calculatePhenoCamTransitionsFromMean.R — combines the estimated PhenoCam SOS transition dates into one file, "phenocamTransitions_fromMeanFiltered.csv"

Fitting models

• runModelIterations.R — helper functions to iteratively run the JAGS models below until convergence and a large enough effective sample size are reached
• createModelCalibration_climatology.R — fits the historical-average null model, saved as paste0(climatologyModelOutputsFolder,siteName,"_climatology_forecast_calibration_varBurn.RData")
• createModelCalibrations_CCmodel.R — calibrates the chlorophyll cycling model for each site and a specified number of included autumn days, saved as paste0(CCmodelOutputsFolder,siteName,"_",n,"_ccModel_forecast_calibration_varBurn.RData")
• createModelCalibrations_triggerModel.R — calibrates an alternate cooling-degree-day (CDD) trigger model for comparison, saved to triggerModelOutputsFolder

Validation site predictions

• uncertaintyAnalysisHindcasts_allSites.R — generates predictions at validation sites from each calibration site's fitted model

Fit plots

• plotClimatologyFits.R — plots the historical-average (climatology) model fits
• plotHindcast_forecasts.R — plots the model fits for various amounts of included autumn data at specified sites
• plotPhenocamTransitions.R — plots the estimated PhenoCam transitions

Analyses

• investigateInflectionPoints.R — determines whether an inflection (SOS) was predicted for each calibration site/year/included-data combination, saving each as paste0(siteName,"_inflectionPointData_15.RData")
• investigateInflectionPoints_OOSsites.R — same, for out-of-sample (validation) site predictions, saving each as paste0(calSite,"_inflectionPointData_OOSsites_mean15_",n,".RData")
• createCRPSpercentageMatrix.R — builds a matrix of the percentage of site-years where the CC model outperforms climatology across differing amounts of included data (Fig. 3c), saving "crpsMat_includedVsDay_Complete.RData", "reorganizedDat_includedVsDay.RData", and "daysOffset_includedVsDay.RData"
• testingImportanceOfTiming.R — fits a random forest model to test whether environmental conditions after a transition are better predictors of that transition than conditions before it (H3)

Manuscript Fig. 5. Importance ranking of predictors in the fitted random forest model (testingImportanceOfTiming.R) used to predict whether the start-of-senescence inflection occurred. The mean product of temperature (T) and photoperiod (D) the week after a date was the most important predictor — evidence for conditions driving senescence rather than a pre-set trigger. Lat and MAT refer to latitude and mean annual temperature, respectively.

Creating figures

• createManuscriptFiguresFINAL.R — creates the manuscript figures (further edited outside of R)
• createCalibrationParameterTable.R — creates the supplementary table of fitted parameter values for each calibration site

Shared helper functions

• ciEnvelope.R — plots a credible-interval polygon on a figure
• calculateStart.R — determines the start of senescence from the second difference (inflection point) of a greenness time series
• generalVariables.R — shared file paths, site lists, and general settings; sourced by nearly every script above (see Setup)

Output file naming conventions

• ERA5 data: paste0(ERA5dataFolder,siteName,"_",start_date,"_",end_date,"_era5TemperatureMembers.nc")
• Rescaled parameters: paste0(dataDirectory,siteName,"_phenopixOutputs.RData")
• Data objects: paste0(dataDirectory,siteName,"_dataFinal.RData")
• Changepoint model outputs (SOS estimate): paste0(transitionEstimateOutputsFolder,siteName,"_",yrName,"_PhenoCam_changePointCurve_varBurn.RData")
• Climatology model outputs: paste0(climatologyModelOutputsFolder,siteName,"_climatology_forecast_calibration_varBurn.RData")
• CC model outputs: paste0(CCmodelOutputsFolder,siteName,"_",n,"_ccModel_forecast_calibration_varBurn.RData")
• CRPS values of out-of-sample predictions: paste0("outOfSampleSites_crps_",calSite,"_183.csv")
• Inflection data, calibration sites: paste0(siteName,"_inflectionPointData_15.RData")
• Inflection data, out-of-sample sites: paste0(calSite,"_inflectionPointData_OOSsites_mean15_",n,".RData")