A leaf energy-balance model from LiDAR and flux tower data

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. There has been renewed interest in the temperature of forest canopies to quantify how environmental factors drive these processes. However, most studies only consider the canopy top. The vertical heterogeneity of forests drives vertical gradients in leaf temperature, however, these gradients are not well-quantified.

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.

Site selection

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

Data sources

We acquired aerial LiDAR point clouds and ancillary NEON products through the neonUtilities R package. For each site, we downloaded the most recent LiDAR acquisition (product ID DP1.30003.001) within a 50 m buffer centered on the flux tower. We then normalized each point cloud with k-nearest neighbor inverse distance weighting in the lidR R package . Points less than 1 m above the ground surface were discarded for compatibility with leaf area index calculation (see below).

We acquired flux tower micrometeorology data through the amerifluxR R package . This package provides an API to download the BASE-BADM product for all sites in the NEON network. Although these data are not processed through the ONEFlux pipeline, we found the BASE-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.

Model overview

Energy balance models estimate leaf temperature by equating incoming and outgoing sources of radiation. We adopted the leaf-level model in , parameterized for needle leaves. We acknowledge a scale mismatch with this choice: the model is developed for individual leaves, but we apply it to describe collections of leaves in a forest canopy. However, this model offers two key advantages that justify this choice. First, the model is computationally efficient in that it solves for canopy temperature directly without numerical optimization. Second, the model provides an additional output in the decoupling coefficient \(\Omega\), which gives further context on drivers of latent heat flux. For brevity, we do not repeat the entire model specification here. Given that needle leaves are often very well coupled with air temperature, it is instructive to focus on the leaf-air temperature difference (\(\Delta T\)):

\[\Delta T \propto R_n - \lambda E\]

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\).

Net radiation submodel

\(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 . In this approach, the LiDAR point cloud is partitioned into horizontal slices of constant thickness (in our case 1 m). For each slice, leaf area density (LAD) is calculated as

\[\mathrm{LAD}_z = \ln \frac{S_{z,in}}{S_{z,out}} \frac{1}{k \Delta z}\]

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.

Model calibration

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 . Precise fisheye lens parameters are not provided in NEON documentation, but based on mention of brand names we assumed that a Nikon Nikkor 10.5 mm fisheye lens was an appropriate assumption supported by 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.

Latent heat flux submodel

\(\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 for complete details on calculation of these variables. Generally, leaf temperature was well-coupled with air (\(\Omega \approx 0\)), so \(\lambda E \approx \lambda E_{imp}\). This endpoint is calculated as

\[\lambda E_{imp} = g_{tot} \frac{D}{P_a}\]

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 bigleaf to calculate canopy-scale \(g_s\) with the function surface.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 to estimate leaf-level \(g_s\) as a function of air temperature, downwelling shortwave radiation, and vapor pressure deficit. This approach enabled us to effectively downscale canopy-scale conductance values to vertical strata within the canopy.

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})\]

Future work

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 . This presents a new challenge in attributing canopy-scale photosynthesis to vertical strata.

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.