Restoration Treatments Enhance Tree Growth and Alter Climatic Constraints During Extreme Drought

Authors

Kyle C. Rodman, John B. Bradford, Alicia M. Formanack, Peter Z. Fulé, David W. Huffman, Thomas E. Kolb, Ana T. Miller-ter Kuile, Donald P. Normandin, Kiona Ogle, Rory J. Pedersen, Daniel R. Schlaepfer, Michael T. Stoddard, Amy E.M. Waltz

Overview

This archive folder includes tabular datasets, statistical models, and model outputs from a study of tree-level radial growth in five ponderosa pine-dominated experimental sites in Arizona, USA. At each site, restoration treatments (i.e., thinning and burning) were applied after an initial forest inventory. Subsequent inventories occurred at each site after treatment implementation, and increment cores were collected from a subset of trees in treated and untreated portions of each site in 2019 to supplement other available data. Descriptions of each dataset are provided below, but further information is available in Rodman et al. (2024). This folder was compressed in .zip format using Windows 11 v. 23H2x64

File Structure

This data package is split into four subfolders, as described below. To replicate our analyses, the R project file “Rodman_MultisiteCoreData.Rproj” should be opened in a new session of R Studio software. This allows the use of relative file paths in the “here” package to re-run scripts described below.

Code

This folder contains all R and Google Earth Engine (javascript) code used to process, analyze, and visualize data in this project. The three subfolders relate to each of these steps.

Code/Step1-DataProcessing

This folder includes the following scripts used to format data and create spatial datasets used in subsequent analyses. All files with the extension “.txt” were written in Java for use in cloud-based computing via Google Earth Engine. All files with the extension “.r” were written in R, for use in local processing.

GEE_GetClimate.txt

Obtains data on average annual climate at each field site from the Daymet dataset (Thornton et al. 2021), and exports them as a .csv file and raster format (geotiff; .tif), for use in making Fig. 1b of Rodman et al. (2024).

Step1a-CoreSubsetting.R

Reads in list of all trees on the sites in 2019, to develop a stratified random sample of increment cores to be collected and dated.

Step1b-FormattingForSOILWAT2.R

Formats field and spatial data as inputs for a daily water balance model (Schlaepfer et al. 2012, Bradford et al. 2014), used to estimate available water capacity - a predictor of tree growth.

Step1c-FormattingCoreData.R

Reads in annual tree ring width information, cleans it, converts it to annual basal area increment, and organizes it in an analysis-ready format. These data are used as a response variable in subsequent analyses.

Step1d-FormattingClimateData.R

Reads in daily weather data for vapor pressure and available water capacity and converts them to seasonal means.

Step1e-FormattingBAData.R

Takes tree lists, calculates live plot-level tree basal area in each field inventory year, and formats these data for use with tree ring data

Step1f-PrepDataForJAGS.R

Takes response data (annual basal area increment) and covariates (basal area, weather, tree size) and merges them into a named list object in .rds format for use as input data in JAGS software.

Step1g-SettingInitials.R

Fits a simplified version of our model, as a frequentist generalized linear mixed model, to create initial values for the Bayesian model in JAGS.

Code/Step2-Analysis

This folder includes the following files/scripts used to analyze data and evaluate fitted models.

ForMonsoon/Growth_v2_4000its/growthData_JAGS.rds

A named list, produced in the “Step1g-PrepDataForJAGS.R” script, and used as an input into the initial JAGS model (4,000 iterations in each of three chains).

ForMonsoon/Growth_v2_4000its/growthInitials.rds

A named list, produced in the “Step1h-SettingInitials.R” script, which gives starting parameter values in the initial JAGS model (4,000 iterations in each of three chains).

ForMonsoon/Growth_v2_4000its/main_job.sh

This is a batch script used to initiate a processing task in the Monsoon high performance computing cluster (with SLURM task manager) at Northern Arizona University to run the initial Bayesian model of tree growth (4,000 iterations across each of three chains).

ForMonsoon/Growth_v2_4000its/model.r

Used to define the structure of the initial Bayesian model of tree growth (4,000 iterations across each of three chains). It is called within the Monsoon cluster using the “main_job.sh” and “script.r” files

ForMonsoon/Growth_v2_4000its/script.r

Runs the initial Bayesian model of tree growth (4,000 iterations across each of three chains) using JAGS and jagsUI. It is called in the Monsoon cluster using the “main_job.sh” file

ForMonsoon/Growth_v3_Final/growthData_JAGS.rds

This object is a named list, produced in the “Step1f-PrepDataForJAGS.R” script, and used as an input into the final JAGS model (20,000 iterations in each of three chains).

ForMonsoon/Growth_v3_Final/growthInitials.rds

This object is a named list, produced in the “Step2b-GrowthModel_FinalInitsAndSettings.R” script, which gives starting parameter values in the final JAGS model (20,000 iterations in each of three chains).

ForMonsoon/Growth_v3_Final/main_job.sh

This is a batch script used to initiate a processing task in the Monsoon high performance computing cluster (with SLURM task manager) at Northern Arizona University, to run the final Bayesian model of tree growth (20,000 iterations across each of three chains).

ForMonsoon/Growth_v3_Final/model.r

Used to define the structure of the final Bayesian model of tree growth (20,000 iterations across each of three chains). It is called within the Monsoon cluster using the “main_job.sh” and “script.r” files

ForMonsoon/Growth_v3_Final/script.r

Runs the final Bayesian model of tree growth (20,000 iterations across each of three chains) using JAGS and jagsUI. It is called in the Monsoon cluster using the “main_job.sh” file

Step2b-GrowthModel_FinalInitsAndSettings.r

Reads in the initial Bayesian model of tree growth (4,000 iterations across each of three chains), and uses it to set initial values and determine the number of iterations needed for convergence in the final model.

Step2c-GrowthModel_EvaluateFit.r

Evaluates the final bayesian model of radial growth (“growthMod_full.rds”; full model object is not included in this archive because of large file sizes - 1.4gb) to assess convergence and summarize/interpret fitted parameters.

Step2d-GrowthModel_IntensityAndTimeSinceTrt.r

Compares time since thinning, time since burning, and treatment intensity with model residuals from hierarchical bayesian regression of tree growth, using only the data from trees that experienced treatment. This evaluates whether treatment effects might have waned throughout the coarse of the study, and creates Figs. S2 and S3 in Rodman et al.

Code/Step3-Visualization

This folder includes the following scripts used to visualize data and the results of the final Bayesian model of tree growth.

Step3a-SiteClimateSummary.r

Takes climate data outputs from “GEE_GetClimate.txt” and other locations in this data package, compares them to the overall distribution of Rocky Mountain ponderosa pine, and plots time series of precipitation and temperature at each of the five experimental sites. This script makes Fig. 1b and 1c in Rodman et al., as well as Fig. S4 (a summary of modeled soil moisture over time)

Step3b-GrowthDescriptive.r

Plots the mean +- 1.96*SE (95% interval) of annual tree basal area increment values at each of the five sites, as well as the timing of thinning and burning. This script makes Fig. 2 in Rodman et al.

Step3c-GrowthModelPlots.r

Extracts posterior summary statistics from the final Bayesian model of tree growth, and creates some visual summaries of the model results - covariate effects, marginal effect plots, and summaries of antecedent climate terms (i.e., lags and seasons that most influence growth). This script makes Figs. 3-6 in Rodman et al.

Data

This subfolder includes raw and summarized data used in Rodman et al., including tree ring data, forest inventories, water balance model outputs, and spatial data extractions (e.g., annual canopy cover, vapor pressure deficit). Individual files are split into additional subfolders, which are described below.

AnalysisReady

This folder includes files that have been cleaned and organized for use in the Bayesian hierarchical model of tree growth, and associated diagnostics.

growthData_merged.csv

This file is the final version of the data, used in statistical analyses for Rodman et al. (2024). Rows represent a single year of growth for an individual tree. Columns are as follows:

  • CoreID: A unique identifier for each tree, created by merging Site, Block, Treatment, Plot, and Tree Number variables from field inventories. Character. 670 unique values
  • Year: A numeric identifier for the year of tree growth in a given row (i.e., “1985”” includes all growth for a given tree in the 1985 calendar year). Integer. Ranges 1985-2018
  • bai_sq_mm: Annual basal area increment (in sq. mm/yr), calculated using the “bai.out” function in the dplr package in R. BAI represents the increase in cross-sectional area of a tree (i.e., CoreID column) in a given year (Year column) at 40 cm above ground level (i.e., coring height). Inputs to the bai.out function include measured ring widths and the bark-free diameter of a tree in 2019 (the year in which an increment core was collected). Numeric, floating point. Ranges 0-12248
  • DBH: Diameter at breast height (in cm; 1.37 m above ground level) for the corresponding tree and year. DBH values for each year were reconstructed from 2019 DBH measurements and measured ring widths. Numeric, floating point. Ranges 0.03-92.37
  • PlotID: A unique identifier for each plot (several trees may be present on a given plot), created by merging Site, Block, Treatment, and Plot variables. Character. 359 unique values
  • Site: Two-letter code identifying the site on which a tree was located. AS = Apache-Sitgreaves, CF = Centennial Forest, FV = Fort Valley, GV = Grandview, MT = Mt. Trumbull. Character. Five unique values.
  • Block: Numeric code identifying the replicated experimental block within a given site, on which a tree was located. Block codes are not unique across sites, and must be merged with the site variable to create a unique identifier. Integer. Unique values = 1, 2, 3, 4
  • Treatment: Numeric code identifying the treatment unit (within a given site and block) on which a tree was located. 1 = Control, 2 = Treated (i.e., thinned and burned), 10 = New Control (i.e., a control unit that was established after the project began). Treatment unit IDs are not unique across sites, and must be merged with the site and block variables to create a unique identifier. Integer. Unique values = 1, 2, 10
  • Plot: Numeric code identifying the plot (within a given site, block, and treatment unit) on which a tree was located. Plot #s are not unique across sites, and the “PlotID” column should be used instead if a unique identifier is needed. Integer. Ranges 1-25
  • lag_bai_sq_mm: Lagged annual basal area increment (in sq. mm/yr) from the previous year, calculated using the same methods as the “bai_sq_mm” column. Used in the modeling process to help account for temporal autocorrelation. Numeric, floating point. Ranges 0-12248
  • BA_ha: Annual live tree basal area estimates in sq. meters/hectare. Calculated using field surveys during inventory years. Values were extrapolated outside of field inventory years using constant values, and linearly interpolated between inventories. Numeric, floating point. Ranges 0.95-84.9
  • Thinned: Binary value of whether a treatment unit on which a tree was located had been thinned prior to the year of growth. If thinning occurred in the spring, then the same year was considered treated, but if it occurred in the summer or fall, then the following year was considered treated. Thus, treated units can be considered thinned or unthinned depending on the year, whereas untreated units are always considered untreated. 0 = Untreated, 1 = Thinned. Integer. Unique values = 0, 1
  • Burned: Binary value of whether a treatment unit on which a tree was located had been burned prior to the year of growth. If prescribed burning occurred in the spring, then the same year was considered treated, but if it occurred in the summer or fall, then the following year was considered treated. Thus, treated units can be considered burned or unburned depending on the year, whereas untreated units are always considered untreated. 0 = Untreated, 1 = Burned. Integer. Unique values = 0, 1
  • meanVPD_hPa_[season]: This description applies to five separate columns in the data, each representing values from a given [season] during the year of tree ring formation. Mean daily vapor pressure deficit, calculated using daily vapor pressure from Daymet (Thornton et al. 2021), and saturation vapor pressure for that day based on improved Magnus equations (Huang 2018). Original values were calculated in units of hPa, but were scaled and centered (i.e., converted to z-scores) within a given site and season. Thus, values represent the standardized spatial and temporal variation in meanVPD within a given site and season. Seasons were defined a priori based on the timing of photosynthesis and radial growth as follows: Winter (prior December to current February), Spring (March and April), Dry Summer (a.k.a. Early Summer; May and June), Monsoon Summer (a.k.a. Late summer; July to September), and Fall (October and November). Numeric, floating point. Ranges -2.81-3.64
  • SWA_mm_000to150_cm_[season]: This description applies to five separate columns in the data, each representing values from a given [season] during the year of tree ring formation. Mean daily available soil water in the top 150 cm of the soil profile (held at > - 3.0 mPa) in a given [season] during the year of ring formation. Values were calculated using SOILWAT2, a daily water balance model that incorporates daily weather, soils information, topography, forest cover, and understory plant cover (Schlaepfer et al. 2012, Bradford et al. 2014). Original values were calculated in units of mm, but were scaled and centered (i.e., converted to z-scores) within a given site and season. Thus, values represent the standardized spatial and temporal variation in available soil water within a given site and season. Seasons were defined a priori based on the timing of photosynthesis and radial growth as follows: Winter (prior December to current February), Spring (March and April), Dry Summer (a.k.a. Early Summer; May and June), Monsoon Summer (a.k.a. Late summer; July to September), and Fall (October and November). Numeric, floating point. Ranges -3.09-4.22
  • meanVPD_hPa_[season]_[lagYear]: This description applies to twenty separate columns in the data, each representing values from a given [season] during one of the four years prior to tree ring formation. [lagYear] identifes the number of years prior to ring formation, where lag1 is the year prior to ring formation and lag4 is four years prior. Units and methods of calculation are equivalent to “meanVPD_hPa_[season]” described above. Numeric, floating point. Ranges -2.81-3.64
  • SWA_mm_000to150_cm_[season]_[lagYear]: This description applies to twenty separate columns in the data, each representing values from a given [season] during one of the four years prior to tree ring formation. [lagYear] identifes the number of years prior to ring formation, where lag1 is the year prior to ring formation and lag4 is four years prior. Units and methods of calculation are equivalent to “SWA_mm_000to150_cm_[season]” described above. Numeric, floating point. Ranges -3.09-4.22
  • TreatmentID: Re-coded version of the “Thinned” column described above. Used for indexing within the Bayesian model of tree growth. 1 = untreated, 2 = treated. Integer. Unique values = 1, 2
  • PlotNum: Re-coded version of the “PlotID” and “Year” columns described above. Gives a number to each unique combination of plot and year, used for indexing within the Bayesian model of tree growth. Integer. Ranges 1-11790
  • BlockNum: Re-coded version of the “Site”, “Block”, and “Year” columns described above. Gives a number to each unique combination of block and year, used for indexing within the Bayesian model of tree growth. Integer. Ranges 1-462
  • BlockTreatNum: Re-coded version of the “Site”, “Block”, “Thinned”, and “Year” columns described above. Gives a number to each unique combination of block, treatment ID, and year, used for indexing within the Bayesian model of tree growth. Integer. Ranges 1-924
  • treeNum: Re-coded version of the “CoreID” column described above. Gives the unique identity of each tree, used for indexing within the Bayesian model of tree growth. Integer. Ranges 1-670
  • treeTrtNum: Re-coded version of the “CoreID” and “Thinned” columns described above. Gives a number to each unique combination of CoreID and treatment, used for indexing within the Bayesian model of tree growth. Integer. Ranges 1-918
growthData_timeSinceTrt.csv

This file is a subset of the original data, used to evaluate the effect of time since thinning and time since burning on residuals of the hierarchical Bayesian model of tree growth. This analysis, used to create Fig. S2 in Rodman et al. (2024), was used to determine if the effects of treatment declined throughout the study period. Descriptions of rows and columns and rows are identical to those in “growthData_merged.csv” described above, but with the exception that rows are restricted to only trees and years that were considered “treated”, so pre-treatment years and untreated units are removed. One newly added column is as follows:

  • resid: Posterior median residuals from the hierarchical Bayesian model of tree growth. Values give the difference between predicted and observed basal area increment (power-transformed; bai^0.2436712) for a given tree and year, with positive values indicating overpredictions, and negative values indicating underpredictions. Positive values in early years after treatment, and negative values at 10+ years after treatment would be indicative of treatment effects that wane over time. We observed no clear relationship between time since thinning/time since burning and model residuals, indicating a persistent treatment effect throughout the course of the study. Numeric, floating point. Ranges -6.62-3.24

Climate

This folder includes information about seasonal vapor pressure deficit or available soil water at the five experimental sites. Raw daily outputs are not included here due to large file sizes (> 10 gb for SOILWAT2 outputs and daily VPD data), but are available upon request. Please contact Kyle Rodman () with requests for these data

climateSummaries_seasonal.csv

This file includes summarized seasonal climate data, used as as covariates in the hierarchical Bayesian model of tree growth. Rows represent each combination of PlotID and year columns (identical to those described in “growthData_merged.csv” above). Columns are as follows:

  • PlotID: A unique identifier for each plot (several trees may be present on a given plot), created by merging Site, Block, Treatment, and Plot variables. Character. 359 unique values
  • Thinned: Binary value of whether a treatment unit on which a plot was located had been thinned prior to the year of a row. If thinning occurred in the spring, then the same year was considered treated, but if it occurred in the summer or fall, then the following year was considered treated. Thus, treated units can be considered thinned or unthinned depending on the year, whereas untreated units are always considered untreated. 0 = Untreated, 1 = Thinned. Integer. Unique values = 0, 1
  • Burned: Binary value of whether a treatment unit on which a plot was located had been burned prior to the year of a row. If prescribed burning occurred in the spring, then the same year was considered treated, but if it occurred in the summer or fall, then the following year was considered treated. Thus, treated units can be considered burned or unburned depending on the year, whereas untreated units are always considered untreated. 0 = Untreated, 1 = Burned. Integer. Unique values = 0, 1
  • Year: A numeric identifier for the year of climate data in a given row (i.e., “1985” includes weather for a given plot in the 1985 calendar year). Integer. Ranges 1980-2022
  • meanVPD_hPa_[season]: This description applies to five separate columns in the data, each representing values from a given [season] during the year identified in the “Year” column. Mean daily vapor pressure deficit, calculated using daily vapor pressure from Daymet (Thornton et al. 2021), and saturation vapor pressure for that day based on improved Magnus equations (Huang 2018). Values are given in hPa (hectopascals). Seasons were defined a priori based on the timing of photosynthesis and radial growth as follows: Winter (prior December to current February), Spring (March and April), Dry Summer (a.k.a. Early Summer; May and June), Monsoon Summer (a.k.a. Late summer; July to September), and Fall (October and November). Numeric, floating point. Ranges 2.1-22.7
  • SWA_mm_000to150_cm_[season]: This description applies to five separate columns in the data, each representing values from a given [season] during the year identified in the “Year” column. Mean daily available soil water in the top 150 cm of the soil profile (held at > - 3.0 mPa) in a given [season] during the year identified in the “Year” column. Values were calculated using SOILWAT2, a daily water balance model that incorporates daily weather, soils information, topography, forest cover, and understory plant cover (Schlaepfer et al. 2012, Bradford et al. 2014). Values are given in units of mm. Seasons were defined a priori based on the timing of photosynthesis and radial growth as follows: Winter (prior December to current February), Spring (March and April), Dry Summer (a.k.a. Early Summer; May and June), Monsoon Summer (a.k.a. Late summer; July to September), and Fall (October and November). Numeric, floating point. Ranges -3.09-4.22

Field

This subfolder includes raw and summarized information on tree growth, forest structure, and the timing of thinning and burning at each of five experimental sites in Arizona, USA, included in Rodman et al.

Cores/CoreSelection

This folder includes files used to create a stratified sample of increment cores for processing and dating in the laboratory. Files are created or used in “Step1a-CoreSubsetting.R” (see notes in this script for additional information) but are not described in detail here.

Cores/RawCoreData/[Site]

This description applies to five subfolders, one for each [Site]. Each subfolder has raw tree ring data, outputs of COFECHA software (Holmes 1983), and various notes from processing for each corresponding site. [Site] refers to AS = Apache-Sitgreaves, CF = Centennial Forest, FV = Fort Valley, GV = Grandview, MT = Mt. Trumbull. Files in each subfolder are as follows:

  • [Site]_COFECHA.txt: Outputs of COFECHA software, verifying the accuracy of dating.
  • [Site]_COFECHA_Scan and Notes: Scanned copy of COFECHA outputs, with notes about individual trees which were checked for accuracy based on COFECHA runs.
  • [Site]_Master.txt: Annual ring width measurements in Tucson format (i.e., fixed-width tabular data with ten ring width measurements per row). Series end dates are identified as “-9999” and values outside of the dated range are blank. Can be read into R use the “read.tucson” function in the dplR package. These files are processed and merged in the “Step1c-FormattingCoreData.r” script described above
Cores/RawCoreData/barkThicknessMod.rds

This object is a generalized linear mixed model, fit using the lme4 package in R, which can be read into R by calling “library(lme4)” and using the “readRDS” function. It was developed to predict bark thickness as a function of tree diameter and other covariates in Rodman et al. 2021. We used it to help develop predictions of basal area increment and reconstructed diameter at breast height. This model is used in the “Step1c-FormattingCoreData.r” script described above

Cores/TreeData/treeData_allSites.csv

This file includes forest inventory data from cored trees at each of five experimental sites in Arizona, USA. These data were used to help select a stratified random sample of cores to be dated and processed, and to link increment cores to tree-level data (e.g., diameter at breast height). Rows represent individual trees included in 2019 surveys that were cored (not a comprehensive list of trees on all sites). Columns are as follows:

  • Site: Two-letter code identifying the site on which a tree was located. AS = Apache-Sitgreaves, CF = Centennial Forest, FV = Fort Valley, GV = Grandview, MT = Mt. Trumbull. Character. Five unique values
  • Year: The year in which a tree was last surveyed. All trees in this dataset were last surveyed in 2019. Integer. One unique value = 2019
  • Block: Numeric code identifying the replicated experimental block within a given site, on which a tree was located. Block codes are not unique across sites, and must be merged with the site variable to create a unique identifier. Integer. Unique values = 1, 2, 3, 4
  • Trt: Numeric code identifying the treatment unit (within a given site and block) on which a tree was located. 1 = Control, 2 = Treated (i.e., thinned and burned), 10 = New Control (i.e., a control unit that was established after the project began). Treatment unit IDs are not unique across sites, and must be merged with the site and block variables to create a unique identifier. Integer. Unique values = 1, 2, 10
  • Plot: Numeric code identifying the plot (within a given site, block, and treatment unit) on which a tree was located. Plot #s are not unique across sites, and the “PlotID” column should be used instead if a unique identifier is needed. Integer. Ranges 1-25
  • Tree_Number: Numeric code identifying the number of a tree (within a given site, block, treatment unit, and plot). Tree Numbers are not unique across sites, and should be merged with Site, Block, Trt, and Plot, to join with tree ring data. Integer. Ranges 1-237
  • Sp_code: Four-letter code (first two letters of a genus and first two letters of a species) identifying the species of a given tree. All values in this dataset are PIPO (ponderosa pine), as this project was restricted to just this species for analysis. Character. One unique value (PIPO)
  • Cond: Condition class for a given tree. 1 = live and relatively healthy and 2 = declining/near death. Values greater than this (dead and in various stages of decomposition) are excluded from this dataset. Integer. Two unique values = 1, 2
  • DBH_x: Field-measured diameter at breast height (1.37m above ground level), in cm, for a given tree in the 2019 field inventory. Numeric, floating point. Range = 0-97
  • TotHt_m: Field-measured height (recorded using a sonar hypsometer), in m, for the top of a given tree in the 2019 field inventory. Numeric, floating point. Range = 1.37-42.7
  • CrownBaseHt_m: Field-measured height (recorded using a sonar hypsometer), in m, for the base of the live crown of given tree in the 2019 field inventory. Numeric, floating point. Range = 0-20.6
  • Inner_Ring: The year of the innermost ring that could be dated for a given increment core. Integer. Range 1466-1992
  • Estimate: An estimate of the number of rings between the pith and “Inner_Ring” of a core. Used to correct for offset from pith to obtain an age estimate for a tree. Values of -999 indicate that a core was not close to the pith and a center date could not be estimated. Integer. Range 0-45
  • Cdate: The estimated center date (i.e., year of pith at coring height; 40 cm above ground level) of a given tree. “Inner_Ring”” minus “Estimate”” columns. Integer. Range 1464-1985
Survival/treeList.csv

This file includes information about all trees taller than breast height (1.37 m) on each of five experimental sites in Arizona, USA. Unlike the “treeData_allSites.csv” file, this list is comprehensive, and includes information about all trees on a site, including uncored trees and those of species other than ponderosa pine. Each row represents a tree in a specific field inventory. Thus, most trees are represented multiple times, as they may have been surveyed in each field inventory. Dead trees were typically not surveyed again after the first inventory in which they were identified as “dead”. New ingrowth (i.e., trees that grew above 1.37 m during the course of the study) was included in only later surveys in which the individuals were taller than 1.37 m. Columns are as follows:

  • Site: Two-letter code identifying the site on which a tree was located. AS = Apache-Sitgreaves, CF = Centennial Forest, FV = Fort Valley, GV = Grandview, MT = Mt. Trumbull. Character. Five unique values
  • Year: The inventory year for the corresponding tree and row in the data. Integer. Range 1997-2019
  • Block: Numeric code identifying the replicated experimental block within a given site, on which a tree was located. Block codes are not unique across sites, and must be merged with the site variable to create a unique identifier. Integer. Unique values = 1, 2, 3, 4
  • Trt: Numeric code identifying the treatment unit (within a given site and block) on which a tree was located. 1 = Control, 2 = Treated (i.e., thinned and burned), 10 = New Control (i.e., a control unit that was established after the project began). Treatment unit IDs are not unique across sites, and must be merged with the site and block variables to create a unique identifier. Integer. Unique values = 1, 2, 10
  • Plot: Numeric code identifying the plot (within a given site, block, and treatment unit) on which a tree was located. Plot #s are not unique across sites, and the “PlotID” column should be used instead if a unique identifier is needed. Integer. Ranges 1-25
  • Tree_Number: Numeric code identifying the number of a tree (within a given site, block, treatment unit, and plot). Tree Numbers are not unique across sites, and should be merged with Site, Block, Trt, and Plot, to join with tree ring data. Integer. Ranges 1-999
  • Sp_code: Four-letter code (first two letters of a genus and first two letters of a species) identifying the species of a given tree. Character. Values are as follows:
    • JUDE: Juniperus deppeana
    • JUOS: Juniperus osteosperma
    • JUSC: Juniperus scopulorum
    • PIED: Pinus edulis
    • PIFL: Pinus flexilis or Pinus strobiformis
    • PIPO: Pinus ponderosa
    • QUGA: Quercus gambelii
    • RONE: Robinia neomexicana
  • Cond: Condition/decay class for a given tree in a given field inventory. Integer. Values are as follows:
    • 1 = live and relatively healthy
    • 2 = declining/near death
    • 3 = recently dead, with most fine branches and tight bark
    • 4 = recently dead, with few fine branches and loose bark
    • 5 = clean snag without many branches or bark
    • 6 = snag that is broken above 1.37 m
    • 7 = snag that is broken below 1.37 m and heavily decomposed
    • 8 = downed log
    • 9 = cut stump
    • 11 = stump hole (evidence that a tree was there at one time, with no woody debris)
  • DBH: Field-measured diameter at breast height (1.37m above ground level), in cm, for a given tree in the corresponding field inventory. Numeric, floating point. Range = 0-105.7
TreatmentDates/LEARN_treatmentYears.csv

This table gives the timing of thinning and burning at each site. Values give the first year in which the treated unit at a given site and block might be considered “Treated” by thinning and/or burning. As described in datasets above, this does not necessarily correspond to the calendar year of treatment - e.g., a treatment year of 2019 could indicate that a site was treated in the summer/fall of the previous year, or the spring of the current year. Numbers do not apply to untreated units at each site, where thinning and burning were never implemented. Rows represent each site and experimental block. Columns are as follows:

  • Site: Two-letter code identifying the site on which a tree was located. AS = Apache-Sitgreaves, CF = Centennial Forest, FV = Fort Valley, GV = Grandview, MT = Mt. Trumbull. Character. Five unique values
  • Block: Numeric code identifying the replicated experimental block within a given site, on which a tree was located. Block codes are not unique across sites, and must be merged with the site variable to create a unique identifier. Integer. Unique values = 1, 2, 3, 4
  • ThinningYear: Year in which a treated unit could first be considered “thinned”. Integer. Range 1999-2007
  • BurningYear1: Year in which a treated unit could first be considered “burned” at least once. NA values indicate that a block was never burned. Integer. Range 2000-2012
  • BurningYear2: Year in which a treated unit could first be considered “burned” at least twice. NA values indicate that a block was not burned twice. Integer. Range 2007-2008
  • BurningYear3: Year in which a treated unit could first be considered “burned” at least three times. NA values indicate that a block was not burned three times. Integer. Range 2015-2015
  • BurningYear4: Year in which a treated unit could first be considered “burned” at least four times. NA values indicate that a block was not burned four times. Integer. Range 2020-2020
bai_sq_mm.csv

This dataset has cleaned and formatted data on tree basal area increment (in sq. mm/year) and DBH (diameter in centimeters at 1.4 m above ground level) for all trees and years included in hierarchical Bayesian models of tree growth. Rows identify individual growth years for a given tree from 1950 to 2018. Data were restricted to 1985 to 2018 for subsequent analyses based on the availability of covariates. This file was produced using the “Step1c-FormattingCoreData.r” script described above. Columns are as follows:

  • Year: A numeric identifier for the year of tree growth in a given row (i.e., “1985”” includes all growth for a given tree in the 1985 calendar year). Integer. Ranges 1950-2018
  • CoreID: A unique identifier for each tree, created by merging Site, Block, Treatment, Plot, and Tree Number variables from field inventories. Character. 670 unique values
  • bai_sq_mm: Annual basal area increment (in sq. mm/yr), calculated using the “bai.out” function in the dplr package in R. BAI represents the increase in cross-sectional area of a tree (i.e., CoreID column) in a given year (Year column) at 40 cm above ground level (i.e., coring height). Inputs to the bai.out function include measured ring widths and the bark-free diameter of a tree in 2019 (the year in which an increment core was collected). Numeric, floating point. Ranges 0-12248
  • DBH: Diameter at breast height (in cm; 1.37 m above ground level) for the corresponding tree and year. DBH values for each year were reconstructed from 2019 DBH measurements and measured ring widths. Numeric, floating point. Ranges 0.03-92.37
ringWidths_mm.csv

This dataset has cleaned and formatted data on annual tree ring width for all trees and years included in hierarchical Bayesian models of tree growth. Rows identify individual growth years for a given tree from 1950 to 2018. Data were restricted to 1985 to 2018 for subsequent analyses based on the availability of covariates. This file was produced using the “Step1c-FormattingCoreData.r” script described above. Columns are as follows:

  • Year: A numeric identifier for the year of tree growth in a given row (i.e., “1985”” includes all growth for a given tree in the 1985 calendar year). Integer. Ranges 1950-2018
  • CoreID: A unique identifier for each tree, created by merging Site, Block, Treatment, Plot, and Tree Number variables from field inventories. Character. 670 unique values
  • ringWidth_mm: Annual ring width (in mm/yr), recorded in the laboratory using a measuring stage (Velmex, Bloomfield, NY) and Measure J2X software (VoorTech Consulting, Holderness, NH). Numeric, floating point. Ranges 0-13.77
  • DBH: Diameter at breast height (in cm; 1.37 m above ground level) for the corresponding tree and year. DBH values for each year were reconstructed from 2019 DBH measurements and measured ring widths. Numeric, floating point. Ranges 0.03-92.37

Spatial

This folder includes several files that provide the locations of each field plot/site, or extractions of spatial data (e.g., annual Landsat-derived canopy cover) at these locations, to be used as covariates in the hierarchical Bayesian model of tree growth. Individual files are further described below.

SiteLocations/LEARN_plotCoords.csv

This file provides the locations of all field plots at these five experimental sites. This may be useful for future users to extract spatial datasets at the locations of field inventories. Many plots that were not included in our analyses (because they contained no cored trees or were on blocks that experienced fire or another disturbance during the course of the study) are included here. Rows represent each field plot. Columns are as follows:

  • Long: Longitude of the center of a field plot. Units are in decimal degrees, and in the WGS84 coordinate system (EPSG:4326). Numeric, floating point. Range -113.21 to -109.62
  • Lat: Latitude of the center of a field plot. Units are in decimal degrees, and in the WGS84 coordinate system (EPSG:4326). Numeric, floating point. Range 34.15 to 36.38
  • Site: Two-letter code identifying the site on which a tree was located. AS = Apache-Sitgreaves, CF = Centennial Forest, FV = Fort Valley, GV = Grandview, MT = Mt. Trumbull. Character. Five unique values
  • Block: Numeric code identifying the replicated experimental block within a given site, on which a tree was located. Block codes are not unique across sites, and must be merged with the site variable to create a unique identifier. Integer. Unique values = 1, 2, 3, 4
  • Treatment: Numeric code identifying the treatment unit (within a given site and block) on which a tree was located. 1 = Control, 2 = Treated (i.e., thinned and burned), 10 = New Control (i.e., a control unit that was established after the project began). Treatment unit IDs are not unique across sites, and must be merged with the site and block variables to create a unique identifier. Integer. Unique values = 1, 2, 10
  • Plot: Numeric code identifying the plot (within a given site, block, and treatment unit) on which a tree was located. Plot #s are not unique across sites, and the “PlotID” column should be used instead if a unique identifier is needed. Integer. Ranges 1-25

Figures

This subfolder is intentionally empty, but included in this data package to permit reproducibility of R scripts in “Code/Step3-Visualization”.

ModelOutputs

This subfolder contains JAGS (Just Another Gibbs Sampler; v. 4.3.0) model outputs and associated model diagnostics created during this project. All .rds files can be opened in R using the “readRDS” function.

ModelOutputs/growthV2_bayesian4000its

Within the ModelOutputs folder, this subfolder contains a hierarchical Bayesian model used to predict tree growth as a function of treatment (i.e., thinned/burned or not treated), plot-level basal area, tree size, and interannual climate variability. The folder also includes associated diagnostics. This model and set of diagnostics are from the initial model (4,000 iterations in each of three chains), used to develop initial values and the number of iterations used in the final model (contained in “growthV3_bayesianFinal” described below). Code and input data used to develop this initial model are in “Code/Step2-Analysis/ForMonsoon/Growth_v2_4000its”

mcmcplots

This subfolder contains traceplots (created with the mcmcplots package in R) for tracked parameters in the the initial model, used to help assess convergence. Individual traceplots are in .png data formats. To view all of these plots in one place, open the “mcmcPlots_Growth4000its.html” file

growthMod_4000.rds

This model object, fit using JAGS v. 4.3.0 and the jagsUI package in R, was created on the Monsoon high performance computing cluster at Northern Arizona University. The model object can be read into R using the ‘readRDS’ function in R to summarize information about posterior distributions of each parameter.

Raftery_4000.rds

This object gives estimated time to convergence (following Raftery and Lewis 1992) for the initial model (4,000 iterations across each of three chains). The object can be read into R using the ‘readRDS’ function in R.

ModelOutputs/growthV3_bayesianFinal

Within the ModelOutputs folder, this subfolder contains a hierarchical Bayesian model used to predict tree growth as a function of treatment (i.e., thinned/burned or not treated), local canopy cover, tree size, and interannual climate. The folder also includes associated diagnostics. These models and diagnostics are from the final model (20,000 iterations in each of three chains), presented in Rodman et al. Code and input data used to develop this final model are in “Code/Step2-Analysis/ForMonsoon/Growth_v3_Final”

mcmcplots_full

This subfolder contains traceplots (created with the mcmcplots package in R) for tracked parameters in the the final model, used to help assess convergence. Individual traceplots are in .png data formats. To view all of these plots in one place, open the “mcmcPlots_GrowthAllits.html” file

Coefs_full.rds

This object gives summarized information on posterior distributions (median; 95% CI) of each model covariate in the final 20,000 iteration model of tree growth. These values are plotted in Fig. 3 of Rodman et al. (2024). The object can be read into R using the ‘readRDS’ function in R.

GOFSamples_full.rds

This object gives an additional 1,000 iterations in each of three chains (using final values from the final model with 20,000 iterations), used to calculate goodness of fit statistics. The object can be read into R using the ‘readRDS’ function in R.

PostSamples_full.rds

This object gives values of all 20,000 iterations in each of three chains for tracked parameters in the final Bayesian model of tree growth, used to create traceplots and evaluate model fit. The object can be read into R using the ‘readRDS’ function in R.

Rhat_growth_jagsui.rds

This object gives convergence statistics (r hat values following Gelman and Rubin 1992) for the final model (20,000 iterations across each of three chains). R hat values close to 1 indicate stability of a tracked parameter across the three chains. The object can be read into R using the ‘readRDS’ function in R.

References

  • Bradford, J. B., Schlaepfer, D. R., & Lauenroth, W. K. (2014). Ecohydrology of Adjacent Sagebrush and Lodgepole Pine Ecosystems: The Consequences of Climate Change and Disturbance. Ecosystems, 17(4), 590–605. https://doi.org/10.1007/s10021-013-9745-1
  • Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical science, 7(4), 457-472.
  • Holmes, R. L. (1983). Computer-assisted quality control in tree-ring dating and measurement. Tree-Ring Bulletin, 43, 69–78.
  • Huang, J. (2018). A Simple Accurate Formula for Calculating Saturation Vapor Pressure of Water and Ice. Journal of Applied Meteorology and Climatology, 57(6), 1265–1272. https://doi.org/10.1175/JAMC-D-17-0334.1
  • Raftery, A. E., & Lewis, S. M. (1992). One Long Run With Diagnostics: Implementation Strategies for Markov Chain Monte Carlo. Statistical Science, 7(4), 493–497.
  • Rodman, K. C., Bradford, J. B., Formanack, A. M., Fule, P. Z., Huffman, D. W., Kolb, T. E., Miller-ter Kuile, A. T., Normandin, D. P., Ogle, K., Pedersen, R. J., Schlaepfer, D. R., Stoddard, M. T., & Waltz, A. E. M. (2024). Restoration treatments enhance tree growth and alter climatic constraints during extreme drought. Ecological Applications.
  • Rodman, K. C., Veblen, T. T., Andrus, R. A., Enright, N. J., Fontaine, J. B., Gonzalez, A. D., Redmond, M. D., & Wion, A. P. (2021). A trait-based approach to assessing resistance and resilience to wildfire in two iconic North American conifers. Journal of Ecology, 109(1), 313–326. https://doi.org/10.1111/1365-2745.13480
  • Schlaepfer, D. R., Lauenroth, W. K., & Bradford, J. B. (2012). Ecohydrological niche of sagebrush ecosystems. Ecohydrology, 5, 453–466. https://doi.org/10.1002/eco.23
  • Thornton, P. E., Shrestha, R., Thornton, M., Kao, S. C., Wei, Y., & Wilson, B. E. (2021). Gridded daily weather data for North America with comprehensive uncertainty quantification. Scientific Data, 8(1), 1–17. https://doi.org/10.1038/s41597-021-00973-0