Detailed methods for my iPoster at AGU 2023.
This post describes the methods used to generate data presented at the 2023 AGU Fall Meeting. However, this project is a work in progress. Assumptions, simplifications, and approaches described here are not necessarily representative of later iterations of this work.
Temperature is the ultimate biophysical parameter. The thermal environment of a plant drives respiration, photosynthesis, and transpiration rates
In this project, we sought to answer two questions:
To answer these questions, we identified nine evergreen forests in the United States and used energy balance modeling to estimate canopy temperature from colocated aerial LiDAR and eddy covariance flux tower micrometeorology. We then validated our model against radiometers mounted at multiple vertical strata for each site. Our data indicate that across a range of climates, growth forms, and canopy structures, leaf temperature is warmer than air in the upper canopy and cooler than air in the lower canopy. Variation in canopy temperature is driven almost completely by air temperature, but variation in the leaf-air temperature difference is driven by stomatal conductance.
We selected study locations from the National Ecological Observatory Network (NEON) which met the following conditions:
This initial search yielded 11 sites. Two (YELL/US-xYE and SOAP/US-xSP) were excluded after analysis of radiometer position and view angle indicated that the radiometers at each site were viewing grasses, not tree canopies. The nine remaining sites are listed in the table below. The canopy heights below differ slightly from NEON metadata because we calculated canopy height from LiDAR returns in the immediate vicinity of each flux tower (see below).
| Site name | Ameriflux ID | NEON ID | Latitude | Longitude | Tower height (m) | Canopy height (m) |
|---|---|---|---|---|---|---|
| Abby Road | US-xAB | ABBY | 45.7624 | -122.3303 | 20 | 13.5 |
| Caribou Creek | US-xBN | BONA | 65.1540 | -147.5026 | 19.5 | 6.5 |
| Delta Junction | US-xDJ | DEJU | 63.8811 | -145.7514 | 22.5 | 10 |
| Jones Ecological Research Center | US-xJE | JERC | 31.1948 | -84.4686 | 47 | 42.5 |
| Rocky Mountain National Park | US-xRM | RMNP | 40.2759 | -105.5459 | 25.5 | 17 |
| Ordway-Swisher Biological Station | US-xSB | OSBS | 29.6893 | -81.9934 | 35.5 | 21 |
| Talladega National Forest | US-xTA | TALL | 32.9505 | -87.3933 | 35.5 | 23.5 |
| Lower Teakettle | US-xTE | TEAK | 37.0058 | -119.0060 | 59 | 26 |
| Wind River Experimental Forest | US-xWR | WREF | 45.8205 | -121.9519 | 74 | 58 |
We acquired aerial LiDAR point clouds and ancillary NEON products through the neonUtilities R packagelidR R package
We acquired flux tower micrometeorology data through the amerifluxR R package BASE-BADM product for all sites in the NEON network. Although these data are not processed through the ONEFlux pipelineBASE-BADM product a more convenient way to access all instrumentation on the tower to reconstruct vertical profiles of meteorological variables. Data were filtered using the default parameters in the amf_filter_base function. For each variable, we determined its measurement height using the table returned by amf_var_info as of 14 November 2023.
Our micrometeorology dataset served two purposes. First, observations of air temperature, wind speed, etc. were key driving variables in the leaf energy balance model. Second, observations of canopy temperature served as a validation dataset for our model. We note that tower observations of canopy temperature are measured by wide-angle radiometers that can be several meters away from foliage. Thus, non-leaf elements (branches, soil, grasses, metal infrastructure, etc.) can contaminate the radiometer signal. We mitigated this issue by discarding observations with a measurement height less than 1 m.
Energy balance models estimate leaf temperature by equating incoming and outgoing sources of radiation. We adopted the leaf-level model in
where \(R_n\) is net isothermal radiation and \(\lambda E\) is latent heat flux. Due to high boundary layer conductance, \(\lambda E\) is driven by leaf stomatal conductance (\(g_s\)). Overall energy balance is therefore driven by two submodels of \(R_n\) and \(g_s\).
\(R_n\) was calculated as
\[R_n = \alpha(\mathrm{SW \downarrow} - \mathrm{SW \uparrow}) + \epsilon(\mathrm{LW \downarrow} - \mathrm{LW \uparrow}) - \epsilon \sigma T_a^4\]where \(\mathrm{SW}\) and \(\mathrm{LW}\) are fluxes of shortwave and longwave radiation, respectively, arrows indicate direction of flux (downwelling vs. upwelling), \(\alpha\) is shortwave reflectance of the leaf surface, \(\epsilon\) is leaf emissivity, \(\sigma\) is the Stefan-Boltzmann constant, and \(T_a\) is air temperature. Our micrometeorology dataset provided top-of-canopy measurements for each radiation flux, a vertical profile of air temperature, and a vertical profile of photosynthetically-active radiation (PAR). We made the following assumptions:
We acknowledge the limitations of these assumptions and intend to revise these in future iterations of the model. Under these assumptions, we are tasked with calculating \(\mathrm{SW\downarrow}\) on a vertical profile from PAR sensors. We calculated the proportion of top-of-canopy shortwave radiation at vertical strata according to
where \(S_{in}\) and \(S_{out}\) are the number of laser pulses entering and exiting the slice, respectively, \(\Delta z\) is slice thickness, and \(k\) is a calibration constant such that the sum of LAD is equal to leaf area index. The proportion of top-of-canopy radiation at height \(z\) is then calculated as
\[\frac{I_z}{I_o} = \exp(-k_I \sum_{i=0}^z \mathrm{LAD}_i \Delta z)\]where \(k_I\) is an additional calibration constant and \(i=0\) corresponds to the canopy top. Note that \(\Delta z\) cancels out when the two equations above are combined. Therefore, once values for \(k\) and \(k_I\) are known, \(\frac{I_z}{I_o}\) can be calculated for any vertical position in the canopy.
Values of whole-canopy leaf area index are required to calibrate \(k\), while vertical profiles of light availability are required to calibrate \(k_I\). To determine leaf area index, (LAI) we utilized digital hemispherical photographs (DHPs) available through the NEON data portal (product ID DP1.10017.001). DHPs are collected in the footprint of the flux tower at roughly shoulder height. Thus, the LAI values we calculate are a site average instead and ignore understory foliage. For each site, we randomly selected 10 DHPs per year and estimated LAI with the hemispheR R package hemispheR. No trends in LAI were observed over observation years, so we calculated site LAI as the median of all images processed. With site LAI known, we then calibrated \(k\) such that the sum of LAD along a vertical profile was equal to site LAI.
We used photosynthetically-active radiation (PAR) sensors on each flux tower to determine an appropriate value for \(k_I\). Sensors are placed at multiple heights, allowing us to determine \(k_I\) as the slope of \(\log(I_z / I_o)\) vs. cumulative LAI. We calculated \(\log(I_z / I_o)\) as the geometric mean of all observations at each measurement height.
\(\lambda E\) was calculated as
\[\lambda E = \Omega \lambda E_{eq} + (1 - \Omega) \lambda E_{imp}.\]This parameterization models \(\lambda E\) as the combination of radiation-limited (\(\lambda E_{eq}\)) and stomatally-imposed (\(\lambda E_{imp}\)) endpoints, related by the decoupling coefficient \(\Omega\). We refer the reader to
where \(g_{tot}\) is the reciprocal sum of leaf boundary layer conductance (\(g_{bH}\)) and stomatal conductance (\(g_s\)), \(D\) is vapor pressure deficit, and \(P_a\) is air pressure. Leaf boundary layer conductance was typically much higher than stomatal conductance, so \(g_{tot} \approx g_s\). Air pressure and vapor pressure deficit could be calculated directly from tower instrumentation. Therefore, this submodel is primarily concerned with estimating \(g_s\).
We estimated \(g_s\) on a vertical profile by downscaling canopy-scale \(g_s\) with our micrometeorological data via a linear mixed effects model. To do so, we first developed a dataset of canopy-scale \(g_s\). Canopy-scale \(g_s\) is estimated from total latent heat flux, but this method is meaningful only when total evapotranspiration is driven by plant transpiration (i.e. \(T/ET \approx 1\)). We filtered latent heat flux observations to those without rainfall in the prior 72 hours to ensure the latent heat flux was dominated by transpiration. We then used the R package bigleafsurface.conductance. Then, we divided canopy-scale \(g_s\) by average site LAI to get an estimate of average leaf \(g_s\) throughout the canopy. We discarded the top 5% of \(g_s\) values at each site since these were physically unreasonable. We then fit a linear mixed effects model with the R package lme4
Selecting the appropriate fixed and random effects for a linear mixed effects model is complicated by large sample sizes. Typically, metrics such as the Akaike or Bayesian information criterion are used to rank different model forms. These criteria select for models that fit data well and against models with many parameters. In our situation, a large sample size (\(n = 25,579\) across all sites) meant that very small improvements in model fit outpaced the penalty associated with model complexity. Thus, information criteria were not useful in selecting a model form. Instead, we fit fixed slopes across all sites and random slopes and intercepts for VPD at each site. In R formula format, the final model took the following form.
\[g_s \sim T_a + \mathrm{SW\downarrow} + \mathrm{VPD} + (\mathrm{VPD} | \mathrm{site})\]Current data indicate that our approach to downscaling canopy-scale \(g_s\) does not accurately represent variation in latent heat flux throughout the canopy. With our current methods, cooling from latent heat is relatively constant on a vertical profile. One avenue to improve this limitation is to couple photosynthesis and stomatal conductance as prior authors have
Our approach to modeling longwave radiation is also quite limited. Longwave radiation data are only available at the top of the canopy, so we cannot derive extinction coefficients as we do for shortwave radiation. If an analytical or empirical relationship with LAI is not known, we may simply retrieve reasonable estimates from the literature.