Cutting out the middleman: calibrating and validating a dynamic vegetation model (ED2-PROSPECT5) using remotely sensed surface reflectance

Canopy radiative transfer is the primary mechanism by which models relate vegetation composition and state to the surface energy balance, which is important to light- and temperature-sensitive plant processes as well as understanding land–atmosphere feedbacks. In addition, certain parameters (e.g., specific leaf area, SLA) that have an outsized influence on vegetation model behavior can be constrained by observations of shortwave reflectance, thus reducing model predictive uncertainty. Importantly, calibrating against radiative transfer outputs allows models to directly use remote sensing reflectance products without relying on highly derived products (such as MODIS leaf area index) whose assumptions may be incompatible with the target vegetation model and whose uncertainties are usually not well quantified. Here, we created the EDR model by coupling the two-stream representation of canopy radiative transfer in the Ecosystem Demography model version 2 (ED2) with a leaf radiative transfer model (PROSPECT-5) and a simple soil reflectance model to predict full-range, high-spectral-resolution surface reflectance that is dependent on the underlying ED2 model state. We then calibrated this model against estimates of hemispherical reflectance (corrected for directional effects) from the NASA Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) and survey data from 54 temperate forest plots in the northeastern United States. The calibration significantly reduced uncertainty in model parameters related to leaf biochemistry and morphology and canopy structure for five plant functional types. Using a single common set of parameters across all sites, the calibrated model was able to accurately reproduce surface reflectance for sites with highly varied forest composition and structure. However, the calibrated model's predictions of leaf area index (LAI) were less robust, capturing only 46 % of the variability in the observations. Comparing the ED2 radiative transfer model with another two-stream soil–leaf–canopy radiative transfer model commonly used in remote sensing studies (PRO4SAIL) illustrated structural errors in the ED2 representation of direct radiation backscatter that resulted in systematic underestimation of reflectance. In addition, we also highlight that, to directly compare with a two-stream radiative transfer model like EDR, we had to perform an additional processing step to convert the directional reflectance estimates of AVIRIS to hemispherical reflectance (also known as “albedo”). In future work, we recommend that vegetation models add the capability to predict directional reflectance, to allow them to more directly assimilate a wide range of airborne and satellite reflectance products. We ultimately conclude that despite these challenges, using dynamic vegetation models to predict surface reflectance is a promising avenue for model calibration and validation using remote sensing data.

1 Introduction

Dynamic vegetation models play a vital role in modern terrestrial ecology and Earth science more generally. The terrestrial carbon cycle is a major biogeochemical feedback in the global climate system (Heinze et al., 2019), and accurate predictions of terrestrial carbon cycling rely on accurate representations of vegetation dynamics (Pacala and Deutschman, 1995). Vegetation also plays an important role in the water cycle and surface energy balance, with major climate implications (Bonan, 2008). In addition, the distribution of tree species, the structure of plant canopies, and many other variables simulated by dynamic vegetation models are also important predictors of biodiversity, making vegetation models an important tool for conservation management (McMahon et al., 2011). Robust calibration and validation of model projections is therefore of broad concern.

Past efforts to calibrate or constrain dynamic vegetation model parameters and states used a variety of data streams. Among these data streams, remote sensing is particularly promising due to its consistent measurement methodology and largely uninterrupted global coverage in recent decades. Data products derived from remote sensing observations have been used to inform, among others, vegetation phenology (Knorr et al., 2010; Viskari et al., 2015) and absorbed photosynthetically active radiation (Peylin et al., 2016; Schürmann et al., 2016; Zobitz et al., 2014). However, there are issues with using derived remote sensing products to calibrate vegetation models. The relationships between remotely sensed surface reflectance and vegetation structure and function are complex and multifaceted. Simple polynomial relationships between spectral indices (e.g., normalized difference vegetation index, NDVI; enhanced vegetation index, EVI) and vegetation properties (e.g., leaf area index, LAI) are often confounded by other ecosystem characteristics, including soil (Myneni and Williams, 1994) and snow (Zhang et al., 2020), or sensor configuration (Fensholt et al., 2004). More sophisticated approaches for estimating vegetation properties based on physically based radiative transfer models (RTMs) face issues of equifinality, whereby many different combinations of vegetation and soil properties can ultimately produce the same modeled surface reflectance (Combal et al., 2003; Lewis and Disney, 2007). Meanwhile, estimating quantities with more indirect relationships to surface reflectance, such as rates of primary productivity, requires a number of assumptions about resource use efficiency and other factors (Running et al., 2004) that can introduce considerable uncertainty and bias into the estimates. Collectively, these issues help explain the large differences in estimates of surface characteristics across different remote sensing instruments (Liu et al., 2018). Robust, pixel-level uncertainty estimates for remote sensing data products would help alleviate some of these concerns, but such estimates are not widely available for most data products. One way to overcome these limitations of derived remote sensing data products while still leveraging the capabilities of remote sensing is to work with lower-level surface reflectance products. Although generating these reflectance products still requires multiple processing steps, such as atmospheric correction, orthorectification, and correction for Sun-sensor geometry effects, all of these processing steps involve significantly fewer assumptions about the relationship between the remotely sensed signal and the surface property or phenomenon of interest than derived products. This can be accomplished by coupling dynamic vegetation models with leaf and canopy radiative transfer models that simulate surface reflectance as a function of known surface characteristics (Knorr and Lakshmi, 2001; Nouvellon et al., 2001; Quaife et al., 2008). Such an approach draws on decades of research on simulation of vegetation optical properties given their structural and biochemical characteristics (Dickinson, 1983; Sellers, 1985; Verhoef, 1984; Lewis and Disney, 2007; Jacquemoud et al., 2009; Pinty et al., 2004; Widlowski et al., 2007, 2015; Hogan et al., 2018) while avoiding the computational and conceptual challenges of inverse parameter estimation in radiative transfer modeling (Combal et al., 2003; Lewis and Disney, 2007).

Instead of coupling a dynamic vegetation model to an external canopy radiative transfer model, we propose a more direct approach of using a vegetation model's own internal representation of canopy radiative transfer. Radiative transfer models have long been an important component of land surface models (Dickinson, 1983; Sellers, 1985). Canopy radiative transfer is the primary mechanism by which models relate vegetation composition and state to the surface energy balance. This is important to both the plants themselves – as many plant processes (including evaporation and enzyme kinetics) are sensitive to temperature (Serbin et al., 2012) – and to the impact of plants on local, regional, and global climate (Bonan, 2008). Canopy radiative transfer also controls how much light is available to plants for photosynthesis and is therefore a first-order driver of plant function (Hikosaka and Terashima, 1995; Robakowski et al., 2004; Niinemets, 2016; Keenan and Niinemets, 2016). Canopy radiative transfer is particularly important to the current generation of demographically enabled dynamic vegetation models, where differences in canopy radiative transfer representations and parameters have major impacts on predicted community composition and biogeochemistry (Loew et al., 2014; Fisher et al., 2018; Viskari et al., 2019). Finally, parameters to which vegetation models are known to be highly sensitive – namely, those related to leaf biochemistry and canopy structure (Dietze et al., 2014; Raczka et al., 2018; Shiklomanov et al., 2020a) – play an important role in canopy radiative transfer. Therefore, calibration and validation against radiative transfer outputs can be an important source of constraint on a variety of model processes.

Our previous work demonstrated that predictions of carbon cycling and community composition by the Ecosystem Demography model version 2 (ED2; Medvigy et al., 2009) are highly sensitive to changes in parameters related to canopy structure and radiative transfer (Viskari et al., 2019). In this study, we build on this work by calibrating and validating the ED2 model using remotely sensed surface reflectance. First, we couple the internal ED2 canopy radiative transfer model to the PROSPECT-5 leaf radiative transfer model (Feret et al., 2008) and the Hapke soil reflectance model (Verhoef and Bach, 2007) to allow ED2 to predict surface reflectance spectra at 1 nm resolution across the complete visible–shortwave infrared (VSWIR) spectral region (400 to 2500 nm). Second, we jointly calibrate this model at 54 sites in the US Midwest and Northeast where coincident vegetation survey data and NASA Airborne Visible/Infrared Imaging Spectrometer – Classic (AVIRIS-Classic) surface reflectance observations are available. We hypothesize that, with known stand composition and informative priors on foliar biochemistry, calibration against airborne imaging spectroscopy will significantly constrain model parameters related to canopy structure. Although the scope of our study is limited to the ED2 model, both the underlying size-and-age structure approximation of ED2 as well as many aspects of its canopy radiative transfer (e.g., two-stream approximation, treatment of leaf angles) are common to other land surface models (e.g., FATES; Koven et al., 2020), meaning the insights from this work are more broadly applicable in model vegetation modeling.

2 Methods

2.1 ED2 model description

The ED2 model simulates plot-level vegetation dynamics and biogeochemistry (Moorcroft et al., 2001; Medvigy et al., 2009; Longo et al., 2019). ED2 has a fundamentally hierarchical structure: the fundamental unit of analysis is a plant cohort – a group of individual plants of similar size, age, and species composition (grouped into plant functional types, PFTs). A group of cohorts with a common disturbance history (i.e., time since last disturbance) constitutes a patch, and a group of co-located patches experiencing the same meteorological conditions constitutes a site. At the spatial scale of this work (60 × 60 m plots; see Sect. 2.3), we assume one patch per site.

Relevant to this work, ED2 includes a multi-layer canopy radiative transfer model that is a generalization of the two-layer, two-stream radiative transfer scheme in the Community Land Model (CLM) v4.5 (Oleson et al., 2013), which in turn is derived from Sellers (1985). At every time step, this model simulates the bi-hemispherical reflectance (BHR, also known as “intrinsic surface albedo” or “blue-sky albedo”; see Schaepman-Strub et al., 2006) as a function of that time step's vegetation composition and canopy structure. A complete description of the model derivation is provided in the supplementary information of Longo et al. (2019), but for completeness, we provide an abbreviated description below.

2.1.1 Radiative transfer parameters

Two-stream radiative transfer theory (Meador and Weaver, 1980) defines the change in radiative flux through a medium in terms of hemispherically integrated upward ($F_i^{\uparrow }$) and downward ($F_i^{\downarrow }$) radiative fluxes via the following system of differential equations (adapting the notation of Yuan et al., 2017):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{d}{F}_i^{\uparrow }}{\mathrm{d}x}}=-\underset{\text{Interception}}{\underbrace{(a_i+{\mathit{\gamma }}_i)F_i^{\uparrow }}}+\underset{\text{Diffuse scatter}}{\underbrace{{\mathit{\gamma }}_i{F}_i^{\downarrow }}}+\underset{\text{Direct backscatter}}{\underbrace{s_i{F}_i^{\odot }}}\text{(2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{F}_i^{\downarrow }}{\mathrm{d}x}}=-\underset{\text{Interception}}{\underbrace{(a_i+{\mathit{\gamma }}_i)F_i^{\downarrow }}}+\underset{\text{Diffuse scatter}}{\underbrace{{\mathit{\gamma }}_i{F}_i^{\uparrow }}}+\underset{\text{Direct scatter}}{\underbrace{s_i^{\prime }{F}_i^{\odot }}},\end{array}$$

where dx represents the vertical change in the total plant area index (combined area of leaves and woody elements), a describes absorption of diffuse radiation, γ describes scattering of diffuse radiation, s and s^′ describe the upward and downward scattering of direct (“beam”) radiation, and $F_i^{\odot }$ is the incident direct (or “beam”) radiative flux at canopy layer i.

Following Sellers (1985), the coefficients above are defined as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{a}_i+{\mathit{\gamma }}_i=\left[\mathrm{1}-\left(\mathrm{1}-{\mathit{\beta }}_i\right)\right]{\displaystyle \frac{\mathrm{1}}{{\overline{\mathit{\mu }}}_i}}\text{(4)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_i={\mathit{\beta }}_i{\mathit{\omega }}_i{\displaystyle \frac{\mathrm{1}}{{\overline{\mathit{\mu }}}_i}}\text{(5)}& {\displaystyle }& {\displaystyle }{s}_i={\displaystyle \frac{\mathrm{1}}{{\mathit{\mu }}_i^{\odot }}}{\mathit{\omega }}_i{\mathit{\beta }}_{\mathrm{0},i}\text{(6)}& {\displaystyle }& {\displaystyle }{s}_i^{\prime }={\displaystyle \frac{\mathrm{1}}{{\mathit{\mu }}_i^{\odot }}}{\mathit{\omega }}_i\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{0},i}\right),\end{array}$$

where ${\overline{\mathit{\mu }}}_i$ is the optical depth per unit plant area index for diffuse radiation, β_i is the backscattering coefficient for diffuse radiation, ω_i is the scattering coefficient for diffuse radiation, ${\mathit{\mu }}_i^{\odot }$ is the optical depth per unit plant area index for direct radiation, and β₀ is the backscattering coefficient for direct radiation.

For a given incident radiation (i.e., solar zenith) angle θ and leaf orientation angle φ, the optical depth is defined as

$$\begin{array}{}\text{(7)}& \mathit{\mu }\left(\mathit{\theta },\mathit{\phi }\right)={\displaystyle \frac{\cos \left(\mathit{\theta }\right)}{G(\mathit{\theta },\mathit{\phi })}},\end{array}$$

where G(θ,φ) is a function describing the projected leaf area. Following Goudriaan (1977), this function can be approximated as

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }G(\mathit{\theta },\mathit{\phi })\approx G^{*}(\mathit{\theta },{\mathit{\chi }}_i)={\mathit{\varphi }}_{\mathrm{1},i}+{\mathit{\varphi }}_{\mathrm{2},i}\cos \left(\mathit{\theta }\right)\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1},i}=\mathrm{0.5}-\mathrm{0.633}{\mathit{\chi }}_i-\mathrm{0.33}{\mathit{\chi }}_i^{\mathrm{2}}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{2},i}=\mathrm{0.877}-\left(\mathrm{1}-\mathrm{2}{\mathit{\varphi }}_{\mathrm{1},i}\right),\end{array}$$

where χ is the leaf orientation factor – a PFT-specific parameter whose values theoretically range from −1 (vertical leaves) to 0 (randomly distributed leaf angles) to 1 (horizontal leaves) but in practice are restricted to $-\mathrm{0.4}\le \mathit{\chi }\le \mathrm{0.6}$.

For a given solar zenith angle (θ_s), the optical depth for direct radiation, ${\mathit{\mu }}_i^{\odot }$, is defined as

$$\begin{array}{}\text{(11)}& {\mathit{\mu }}_i^{\odot }=\mathit{\mu }({\mathit{\theta }}_{\mathrm{s}},{\mathit{\chi }}_i)={\displaystyle \frac{\cos \left({\mathit{\theta }}_{\mathrm{s}}\right)}{G^{*}({\mathit{\theta }}_{\mathrm{s}},{\mathit{\chi }}_i)}}.\end{array}$$

The optical depth for diffuse radiation, ${\overline{\mathit{\mu }}}_i$, is defined as the integral of Eq. (7) over all zenith angles:

$$\begin{array}{}\text{(12)}& {\overline{\mathit{\mu }}}_i=\underset{\mathrm{0}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}{\displaystyle \frac{\cos \left(\mathit{\theta }\right)}{G^{*}(\mathit{\theta },{\mathit{\chi }}_i)}}\mathrm{d}\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{{\mathit{\varphi }}_{\mathrm{2},i}}}\left[\mathrm{1}+{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{1},i}}{{\mathit{\varphi }}_{\mathrm{2},i}}}\ln {\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{1},i}}{{\mathit{\varphi }}_{\mathrm{1},i}+{\mathit{\varphi }}_{\mathrm{2},i}}}\right].\end{array}$$

Following Sellers (1985), diffuse scattering (ω) and backscattering (β) coefficients of canopy elements (leaves or stems) are defined as a function of those elements' reflectance (R) and transmittance (T; wood transmittance is assumed to be zero). (We use index p to refer to PFT and p(i) to refer to the PFT of cohort i.)

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{\mathit{\omega }}_{i,\mathrm{leaf}}=R_{p\left(i\right),\mathrm{leaf}}+T_{p\left(i\right),\mathrm{leaf}}\text{(14)}& {\displaystyle }& {\displaystyle }{\mathit{\omega }}_{i,\mathrm{wood}}=R_{p\left(i\right),\mathrm{wood}}\text{(15)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{i,\mathrm{leaf}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\omega }}_i}}\left[R_{i,\mathrm{leaf}}+T_{i,\mathrm{leaf}}+\left(R_{i,\mathrm{leaf}}-T_{i,\mathrm{leaf}}\right)J\left({\mathit{\chi }}_i\right)\right]\text{(16)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{i,\mathrm{wood}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\omega }}_i}}\left[R_{i,\mathrm{wood}}+R_{i,\mathrm{wood}}J\left({\mathit{\chi }}_i\right)\right],\end{array}$$

where J(χ_i) captures the effect of leaf and branch inclination and is approximated as (similarly to Oleson et al., 2013)

$$\begin{array}{}\text{(17)}& J\left({\mathit{\chi }}_i\right)={\displaystyle \frac{\mathrm{1}+{\mathit{\chi }}_i}{\mathrm{2}}}.\end{array}$$

Both ω and β are calculated independently for leaves and wood and then averaged based on the relative effective area of leaves (L_i) and wood (W_i) within a canopy layer.

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }{\mathit{\omega }}_i={\mathit{\omega }}_{i,\mathrm{leaf}}{\displaystyle \frac{L_i}{L_i+W_i}}+{\mathit{\omega }}_{i,\mathrm{wood}}\left(\mathrm{1}-{\displaystyle \frac{L_i}{L_i+W_i}}\right)\text{(19)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_i={\mathit{\beta }}_{i,\mathrm{leaf}}{\displaystyle \frac{L_i}{L_i+W_i}}+{\mathit{\beta }}_{\mathrm{wood}}\left(\mathrm{1}-{\displaystyle \frac{L_i}{L_i+W_i}}\right)\end{array}$$

To account for non-uniform distribution of leaves within a canopy, ED2 has a PFT-specific clumping factor (q) parameter that serves as a scaling factor on leaf area index. Therefore, the effective leaf area index (L) is related to the true leaf area index (LAI) by

$$\begin{array}{}\text{(20)}& L_i={\mathrm{LAI}}_i\times q_{p\left(i\right)}.\end{array}$$

The leaf area of a cohort (LAI_i) is calculated as a function of leaf biomass (B_leaf,i, kgC plant⁻¹), specific leaf area (SLA_p, m² kgC⁻¹), and stem density (n_plant, plants m⁻²):

$$\begin{array}{}\text{(21)}& {\mathrm{LAI}}_i=n_{\mathrm{plant},i}{B}_{\mathrm{leaf},i}{\mathrm{SLA}}_{p\left(i\right)}.\end{array}$$

In turn, B_leaf,i is calculated from cohort diameter at breast height (DBH_i, cm) according to the following allometric equations:

$$\begin{array}{}\text{(22)}& B_{\mathrm{leaf},i}=\mathrm{b}\mathrm{1}{\mathrm{Bl}}_{p\left(i\right)}{\mathrm{DBH}}_i^{\mathrm{b}\mathrm{2}{\mathrm{Bl}}_{p\left(i\right)}},\end{array}$$

where b1Bl_p(i) and b2Bl_p(i) are PFT-specific parameters. The wood area of a cohort (WAI_i) is calculated directly from DBH according to a similar allometric equation:

$$\begin{array}{}\text{(23)}& {\mathrm{WAI}}_i=n_{\mathrm{plant},i}\mathrm{b}\mathrm{1}{\mathrm{Bw}}_{p\left(i\right)}{\mathrm{DBH}}_i^{\mathrm{b}\mathrm{2}{\mathrm{Bw}}_{p\left(i\right)}},\end{array}$$

where b1Bw_p(i) and b2Bw_p(i) are PFT-specific parameters.

Backscattering of direct radiation (${\mathit{\beta }}_i^{\odot }$) is defined as a function of single scattering albedo (α_s(θ_s)):

$$\begin{array}{}\text{(24)}& {\mathit{\beta }}_i^{\odot }={\displaystyle \frac{{\overline{\mathit{\mu }}}_i+{\mathit{\mu }}_i^{\odot }}{{\overline{\mathit{\mu }}}_i}}{\mathit{\alpha }}_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{s}}\right).\end{array}$$

Single scattering albedo (α_s(θ_s)) is in turn defined as an integral over all illumination angles (ϑ), following Sellers (1985):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{s}}\right)=\mathit{\omega }\underset{\mathrm{0}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\theta }}_{\mathrm{s}},\mathit{\vartheta })\cos \left(\mathit{\vartheta }\right)}{G({\mathit{\theta }}_{\mathrm{s}},\mathit{\phi })\cos \left({\mathit{\theta }}_{\mathrm{s}}\right)+G(\mathit{\vartheta },\mathit{\phi })\cos \left(\mathit{\vartheta }\right)}}\\ \text{(25)}& {\displaystyle }& {\displaystyle }\sin \left(\mathit{\vartheta }\right)\mathrm{d}\mathit{\vartheta }\text{(26)}& {\displaystyle }& {\displaystyle }\mathrm{\Gamma }(\mathit{\theta },\mathit{\vartheta })=G(\mathit{\theta },\mathit{\phi })G(\mathit{\vartheta },\mathit{\phi })P(\mathit{\theta },\mathit{\vartheta })\text{(27)}& {\displaystyle }& {\displaystyle }\underset{-\frac{\mathit{\pi }}{\mathrm{2}}}{\overset{\frac{\mathit{\pi }}{\mathrm{2}}}{\int }}P(\mathit{\theta },\mathit{\vartheta })G(\mathit{\vartheta },\mathit{\phi })\sin \left(\mathit{\vartheta }\right)\mathrm{d}\mathit{\vartheta }=\mathrm{1},\end{array}$$

where P(θ,ϑ) is the scattering phase function that defines the relative fraction of scattered flux in any direction relative to the projected leaf area in that direction (Dickinson, 1983). Substituting $G=G^{*}$ (Eq. 8), and assuming uniform scattering (i.e., $P(\mathit{\theta },\mathit{\vartheta })=\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}$, and therefore $\mathrm{\Gamma }(\mathit{\theta },\mathit{\vartheta })=\frac{G(\mathit{\theta },\mathit{\phi })}{\mathrm{2}}$; see a detailed discussion of these assumptions in Yuan et al., 2017) gives the following analytical solution to this integral:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}& {\mathit{\alpha }}_{\mathrm{s},i}\left(\mathit{\theta }\right)={\displaystyle \frac{{\mathit{\omega }}_i}{\mathrm{2}\left(\mathrm{1}+{\mathit{\varphi }}_{\mathrm{2},i}{\mathit{\mu }}_i^{\odot }\right)}}\\ & \left[\mathrm{1}-{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{1},i}{\mathit{\mu }}_i^{\odot }}{\mathrm{1}+{\mathit{\varphi }}_{\mathrm{2},i}{\mathit{\mu }}_i^{\odot }}}\ln \left({\displaystyle \frac{\mathrm{1}+({\mathit{\varphi }}_{\mathrm{1},i}+{\mathit{\varphi }}_{\mathrm{2},i}){\mathit{\mu }}_i^{\odot }}{{\mathit{\varphi }}_{\mathrm{1},i}{\mathit{\mu }}_i^{\odot }}}\right)\right].\end{array}\end{array}$$

2.1.2 Solution for the multi-layer canopy

In ED2, the direct radiation profile, $F_i^{\odot }$, is governed by exponential decay, following Beer's law:

$$\begin{array}{}\text{(29)}& {\displaystyle }& {\displaystyle }{F}_i^{\odot }=F_{i+\mathrm{1}}^{\odot }\exp \left(-{\displaystyle \frac{{\mathrm{TAI}}_i}{{\mathit{\mu }}_i^{\odot }}}\right)\text{(30)}& {\displaystyle }& {\displaystyle }{F}_{n+\mathrm{1}}^{\odot }=F_{\mathrm{sky}}^{\odot },\end{array}$$

where $F_{\mathrm{sky}}^{\odot }$ is the incident direct (“beam”) radiation from the atmosphere, a prescribed input; and TAI_i is the total plant area index, defined as the sum of effective leaf area index (L_i; Eq. 20) and wood area index (WAI).

For n cohorts, the full diffuse canopy radiation profile in ED2 is defined by a vector A of size 2n+2 that contains the upward ($F_i^{\uparrow }$) and downward ($F_i^{\downarrow }$) radiative fluxes for every “interface” immediately below cohort i; therefore, $F_{(n+\mathrm{1})}^{\uparrow }$ refers to the upward diffuse radiative flux from the top of the canopy towards the atmosphere (the quantity used to calculate the albedo), $F_{(n+\mathrm{1})}^{\downarrow }$ refers to the downward diffuse radiative flux from the atmosphere into the top of the canopy, $F_{\mathrm{1}}^{\uparrow }$ refers to the upward diffuse radiative flux from the ground into the canopy layer of the shortest cohort, and $F_{\mathrm{1}}^{\downarrow }$ refers to the downward diffuse radiative flux from the canopy layer of the shortest cohort towards the ground. To derive each of these F terms, ED2 uses the following analytical solution for Eqs. (2) and (1) (see Longo et al., 2019, Sect. S12, for a full derivation):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{F}_i^{\downarrow }=x_{(\mathrm{2}i-\mathrm{1})}{\mathit{\gamma }}_i^{+}\exp \left(-{\mathit{\lambda }}_i\mathrm{TAI}\right)+x_{\left(\mathrm{2}i\right)}{\mathit{\gamma }}_i^{-}\exp \left(+{\mathit{\lambda }}_i\mathrm{TAI}\right)\\ \text{(31)}& {\displaystyle }& {\displaystyle }+{\mathit{\delta }}^{+}\exp \left(-{\displaystyle \frac{{\mathrm{TAI}}_i}{{\mathit{\mu }}_i^{\odot }}}\right){\displaystyle }& {\displaystyle }{F}_i^{\uparrow }=x_{(\mathrm{2}i-\mathrm{1})}{\mathit{\gamma }}_i^{-}\exp \left(-{\mathit{\lambda }}_i\mathrm{TAI}\right)+x_{\left(\mathrm{2}i\right)}{\mathit{\gamma }}_i^{+}\exp \left(+{\mathit{\lambda }}_i\mathrm{TAI}\right)\\ \text{(32)}& {\displaystyle }& {\displaystyle }+{\mathit{\delta }}^{-}\exp \left(-{\displaystyle \frac{{\mathrm{TAI}}_i}{{\mathit{\mu }}_i^{\odot }}}\right),\end{array}$$

where x_i is a vector of cohort-specific unknowns, and the remaining terms are

$$\begin{array}{}\text{(33)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_i^{\pm }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\pm \sqrt{{\displaystyle \frac{\mathrm{1}-{\mathit{\omega }}_i}{\mathrm{1}-(\mathrm{1}-\mathrm{2}{\mathit{\beta }}_i){\mathit{\omega }}_i}}}\right)\text{(34)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}_i^{\pm }={\displaystyle \frac{\left({\mathit{\kappa }}^{+}\pm {\mathit{\kappa }}^{-}\right){\mathit{\mu }}_i^{\odot \mathrm{2}}}{\mathrm{2}\left(\mathrm{1}-{\mathit{\lambda }}_i^{\mathrm{2}}{\mathit{\mu }}_i^{\odot \mathrm{2}}\right)}}\text{(35)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_i^{\mathrm{2}}={\displaystyle \frac{\left[\mathrm{1}-(\mathrm{1}-\mathrm{2}{\mathit{\beta }}_i){\mathit{\omega }}_i\right]\left(\mathrm{1}-{\mathit{\omega }}_i\right)}{{\overline{\mathit{\mu }}}_i^{\mathrm{2}}}}\text{(36)}& {\displaystyle }& {\displaystyle }{\mathit{\kappa }}_i^{+}=-\left[{\displaystyle \frac{\mathrm{1}-\left(\mathrm{1}-\mathrm{2}{\mathit{\beta }}_i\right){\mathit{\omega }}_i}{{\overline{\mathit{\mu }}}_i}}+{\displaystyle \frac{\mathrm{1}-\mathrm{2}{\mathit{\beta }}_i}{{\mathit{\mu }}_i^{\odot }}}\right]{\displaystyle \frac{{\mathit{\omega }}_i{F}_{(i+\mathrm{1})}^{\odot }}{{\mathit{\mu }}_k^{\odot }}}\text{(37)}& {\displaystyle }& {\displaystyle }{\mathit{\kappa }}_i^{-}=-\left[{\displaystyle \frac{\left(\mathrm{1}-\mathrm{2}{\mathit{\beta }}_i\right)\left(\mathrm{1}-{\mathit{\omega }}_i\right)}{{\overline{\mathit{\mu }}}_i}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\mu }}_i^{\odot }}}\right]{\displaystyle \frac{{\mathit{\omega }}_i{F}_{(i+\mathrm{1})}^{\odot }}{{\mathit{\mu }}_k^{\odot }}}.\end{array}$$

The problem of solving for x_i in Eqs. (31) and (32) can be written as a matrix equation:

$$\begin{array}{}\text{(38)}& \mathbf{S}\mathit{x}=\mathit{Y},\end{array}$$

where $\mathit{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_{\mathrm{2}n+\mathrm{1}},x_{\mathrm{2}n+\mathrm{2}}\right)$, $\mathit{y}=\left(y_{\mathrm{1}},y_{\mathrm{2}},\mathrm{\dots },y_{(\mathrm{2}n+\mathrm{1})},y_{(\mathrm{2}n+\mathrm{2})}\right)$, and S is a $(\mathrm{2}n+\mathrm{2})\times (\mathrm{2}n+\mathrm{2})$ tridiagonal matrix. To solve this matrix equation, ED2 defines the following boundary conditions: at the top of the canopy ($i=n+\mathrm{1}$), $F_{(n+\mathrm{1})}^{\downarrow }\equiv F_{\mathrm{sky}}^{\downarrow }$, the incident diffuse flux from the atmosphere (a prescribed input); ${\mathrm{TAI}}_{(n+\mathrm{1})}=\mathrm{0}$; ${\overline{\mathit{\mu }}}_{(n+\mathrm{1})}=\mathrm{1}$; ${\mathit{\omega }}_{(n+\mathrm{1})}=\mathrm{1}$ (no absorption); and ${\mathit{\beta }}_{(n+\mathrm{1})}={\mathit{\beta }}_{(n+\mathrm{1})}^{\odot }=\mathrm{0}$ (no scattering; all radiance is transmitted). At the bottom of the canopy (i=1), we redefine the ground scattering, ω_g, based on a soil radiative transfer model (Sect. 2.2). With these boundary conditions, the elements of S are given by

$$\begin{array}{}\text{(39)}& \begin{array}{rl}& S_{\mathrm{1},\mathrm{1}}=\left({\mathit{\gamma }}_{\mathrm{1}}^{-}-{\mathit{\omega }}_g{\mathit{\gamma }}^{+}\right)\exp (-{\mathit{\lambda }}_{\mathrm{1}}{\mathrm{TAI}}_{\mathrm{1}})\\ & S_{\mathrm{1},\mathrm{2}}=\left({\mathit{\gamma }}_{\mathrm{1}}^{+}-{\mathit{\omega }}_g{\mathit{\gamma }}^{-}\right)\exp (+{\mathit{\lambda }}_{\mathrm{1}}{\mathrm{TAI}}_{\mathrm{1}})\\ & S_{(\mathrm{2}i,\mathrm{2}i-\mathrm{1})}={\mathit{\gamma }}_i^{+}\\ & S_{(\mathrm{2}i,\mathrm{2}i)}={\mathit{\gamma }}_i^{-}\\ & S_{(\mathrm{2}i,\mathrm{2}i+\mathrm{1})}=-{\mathit{\gamma }}_{(i+\mathrm{1})}^{+}\exp (-{\mathit{\lambda }}_{(i+\mathrm{1})}{\mathrm{TAI}}_{(i+\mathrm{1})})\\ & S_{(\mathrm{2}i,\mathrm{2}i+\mathrm{2})}=-{\mathit{\gamma }}_{(i+\mathrm{1})}^{-}\exp (+{\mathit{\lambda }}_{(i+\mathrm{1})}{\mathrm{TAI}}_{(i+\mathrm{1})})\\ & S_{(\mathrm{2}i+\mathrm{1},\mathrm{2}i-\mathrm{1})}={\mathit{\gamma }}_i^{-}\\ & S_{(\mathrm{2}i+\mathrm{1},\mathrm{2}i)}={\mathit{\gamma }}_i^{+}\\ & S_{(\mathrm{2}n,\mathrm{2}n+\mathrm{1})}=-{\mathit{\gamma }}_{(n+\mathrm{1})}^{+}\exp (-{\mathit{\lambda }}_{(n+\mathrm{1})}{\mathrm{TAI}}_{(n+\mathrm{1})})\\ & S_{(\mathrm{2}n,\mathrm{2}n+\mathrm{2})}=-{\mathit{\gamma }}_{(n+\mathrm{1})}^{-}\exp (+{\mathit{\lambda }}_{(n+\mathrm{1})}{\mathrm{TAI}}_{(n+\mathrm{1})}),\end{array}\end{array}$$

and the elements of y are given by

$$\begin{array}{}\text{(40)}& \begin{array}{rl}& y_{\mathrm{1}}={\mathit{\omega }}_{\mathrm{0}}{F}_{\mathrm{1}}^{\odot }-\left({\mathit{\delta }}_{\mathrm{1}}^{-}-{\mathit{\omega }}_g{\mathit{\delta }}_{\mathrm{1}}^{+}\right)\exp \left(-{\displaystyle \frac{\mathrm{TAI}}{{\mathit{\mu }}_{\mathrm{1}}^{\odot }}}\right)\\ & y_{\left(\mathrm{2}i\right)}={\mathit{\delta }}_{(i+\mathrm{1})}^{+}\exp \left(-{\displaystyle \frac{{\mathrm{TAI}}_{(i+\mathrm{1})}}{{\mathit{\mu }}_{(i+\mathrm{1})}^{\odot }}}\right)-{\mathit{\delta }}_i^{+}\\ & y_{(\mathrm{2}i+\mathrm{1})}={\mathit{\delta }}_{(i+\mathrm{1})}^{-}\exp \left(-{\displaystyle \frac{{\mathrm{TAI}}_{(i+\mathrm{1})}}{{\mathit{\mu }}_{(i+\mathrm{1})}^{\odot }}}\right)-{\mathit{\delta }}_i^{-}\\ & y_{(\mathrm{2}n+\mathrm{2})}=F_{\mathrm{sky}}^{\downarrow }-{\mathit{\delta }}_{(n+\mathrm{1})}^{+}.\end{array}\end{array}$$

Finally, the surface albedo (ρ) is defined as the fraction of the total radiative flux incident on the canopy ($F_{\mathrm{sky}}^{\odot }+F_{\mathrm{sky}}^{\downarrow }$) that is reflected:

$$\begin{array}{}\text{(41)}& \mathit{\rho }={\displaystyle \frac{F_{(n+\mathrm{1})}^{\uparrow }}{F_{\mathrm{sky}}^{\odot }+F_{\mathrm{sky}}^{\downarrow }}}.\end{array}$$

2.2 ED2-PROSPECT coupling

By default, ED2 performs the canopy shortwave radiative transfer calculations described in two broad spectral bands: visible (400–700 nm) and near infrared (700–2500 nm). For each of these regions, ED2 has user-defined prescribed, PFT-specific leaf and wood reflectance and transmittance values and calculates soil reflectance as the average of constant wet and dry soil reflectance values weighted by the relative soil moisture (0 indicates fully dry; 1 indicates fully wet). In this study, we modified ED2 to perform the same canopy radiative transfer calculations but in 1 nm increments across the range 400–2500 nm. We then simulated leaf reflectance and transmittance using the PROSPECT-5 leaf RTM, which has the following five parameters: effective number of leaf mesophyll layers (N, unitless, >=1), total chlorophyll content (Cab, µg cm⁻²), total carotenoid content (Car, µg cm⁻²), water content (Cw, g cm⁻²), and dry matter content (Cm, g cm⁻²) (Feret et al., 2008). For wood reflectance, we used a single representative spectrum – the mean of all wood spectra from Asner (1998), resampled to 1 nm resolution – for all PFTs. For soil scattering (ω_g), we used the simple Hapke soil submodel used in the soil–leaf–canopy RTM (Verhoef and Bach, 2007), whereby soil reflectance is the average of prescribed wet and dry soil reflectance spectra weighted by a relative soil moisture parameter (ϱ_soil, unitless, 0–1). The final coupled ED2-PROSPECT canopy radiative transfer model (hereafter known as “EDR”) has 12 parameters for each PFT – five parameters for PROSPECT, specific leaf area, two parameters each for the leaf and wood allometries, and clumping (q) and leaf orientation (χ) factors – and one site-specific parameter – the relative soil moisture (Table 1).

EDR shares many assumptions and internal coefficients with SAIL (Verhoef, 1984; Verhoef and Bach, 2007), a canopy radiative transfer model that is popular in the optical remote sensing community due to its ability to simulate both hemispherical and directional reflectance. Unlike EDR's vertically heterogeneous canopy, SAIL takes only a single homogenous canopy layer as an input, which precludes a valid comparison of the two models' simulations for real heterogeneous sites. Nevertheless, to help identify possible structural issues with EDR and to explore differences between hemispherical and directional reflectance streams, we compared the sensitivities of EDR and SAIL to LAI and solar zenith angles for a single-layer homogeneous canopy.

2.3 Site and data description

For model calibration, we selected 54 sites from the NASA forest functional type (FFT) field campaign that contained plot-level inventory data coincident with observations of NASA's AVIRIS-Classic. A full description of this dataset is provided in Singh et al. (2015). Briefly, each site consisted of a 60 × 60 m transect within which forest inventory data (stem density, species identity, and diameter at breast height, DBH) were collected. These sites are located in the United States upper Midwest, northern New York, and western Maryland (Fig. 1), and include stands dominated by either evergreen or deciduous trees and spanning a wide range of structures, from dense groups of saplings to sparse groups of large trees (Fig. 2).

Figure 1 Map of sites used in this analysis. Sites shown in Fig. 6 are labeled.

For this study, because our goal was only to calibrate the ED2 canopy radiative transfer parameters and not to evaluate ED2 predictions of vegetation dynamics, we prescribed the vegetation composition at each site based on the inventory data described above. We grouped the tree species in these sites into five different PFTs as defined by ED2: early successional hardwood, northern mid-successional hardwood, late successional hardwood, northern pine, and late successional conifer. The mappings of tree species onto these PFTs are provided as a CSV-formatted table in the file inst/pfts-species.csv in the source code repository for this project (see Code and Data Availability section).

Figure 2 Stand structure and composition characteristics of sites selected for analysis. The dashed line is a forest self-thinning curve (see Zeide, 2010) parameterized based on an analysis of US Forest Service Forest Inventory and Analysis (FIA) data (Travis Andrews, unpublished data). Colors indicate dominant PFT, calculated as the PFT contributing most to site total stem-density-weighted DBH. Sites shown in Fig. 6 are labeled.

AVIRIS-Classic measures the directional radiance of a surface from 365 to 2500 nm at approximately 10 nm increments. Atmospheric correction routines use this level-1 radiance product to estimate the surface reflectance (technically, hemispherical–directional reflectance factor, HDRF, sensu Schaepman-Strub et al., 2006) – a quantity that can be more directly related to intrinsic physical properties of the surface. For this study, in addition to the standard atmospheric correction and orthorectification conducted by NASA Jet Propulsion Laboratory (JPL), the AVIRIS data were also cross-track illumination corrected and bidirectional reflectance distribution function (BRDF) corrected, following the procedure of Lucht et al. (2000). Briefly, this BRDF correction estimates “intrinsic surface albedo” – the quantity that is simulated by EDR – from directional reflectance data by fitting a polynomial approximation to the Ross–Li semi-empirical BRDF model and then integrating this model over all angles. The full AVIRIS processing pipeline for the AVIRIS data (including the BRDF approximation) we used is described in Singh et al. (2015).

Because of unrealistic values in the shortwave infrared spectral region (>1300 nm) in the AVIRIS observations (likely caused by faulty atmospheric correction), we only used observations from 400 to 1300 nm for model calibration and validation. Following Shiklomanov et al. (2016), we used the relative spectral response functions of AVIRIS-Classic to relate the 1 nm EDR predictions to the 10 nm AVIRIS-Classic measurements.

For each AVIRIS-Classic observation, we retrieved the solar zenith angle directly from the flight-line metadata, where available, and calculated it based on the local time and position if not. In addition, we retrieved the relative fraction of diffuse vs. direct incident radiation from the hourly Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) meteorological reanalysis (Gelaro et al., 2017) for each observation's location and time (rounded to the nearest hour).

2.4 Model calibration

To estimate EDR parameters from AVIRIS observations, we used a Bayesian approach that builds on our previous work at the leaf scale (Shiklomanov et al., 2016). For a parameter vector Θ and matrix of observations X, the typical form of Bayes' rule is given by

$$\begin{array}{}\text{(42)}& \underset{\text{Posterior}}{\underbrace{P\left(\mathbf{\Theta }\right|\mathbf{X})}}\sim \underset{\text{Likelihood}}{\underbrace{P\left(\mathbf{X}\right|\mathbf{\Theta })}}\underset{\text{Prior}}{\underbrace{P\left(\mathbf{\Theta }\right)}}.\end{array}$$

Rather than performing a separate calibration at each site, we performed a single joint calibration across all sites. Therefore, our overall likelihood (P(X|Θ)) was the product of the likelihood at each site (P(X_s|Θ), for site s):

$$\begin{array}{}\text{(43)}& P\left(\mathbf{X}\right|\mathbf{\Theta })=\prod _{\mathrm{s}}P({\mathbf{X}}_{\mathrm{s}}\left|\mathbf{\Theta }\right).\end{array}$$

The likelihood at each site s is based on how well EDR-predicted albedo (R_pred,s) matches that site's observed AVIRIS albedo (X_s) given the known forest composition at that site (comp_s) and the current estimate of the overall parameter vector. Similar to Shiklomanov et al. (2016), we assumed the residual error between predicted and observed reflectance followed a multivariate normal distribution (MvNormal):

$$\begin{array}{}\text{(44)}& P\left({\mathbf{X}}_{\mathrm{s}}\right|\mathbf{\Theta })=\mathrm{MvNormal}({\mathbf{X}}_{\mathrm{s}}|R_{\mathrm{pred},\mathrm{s}},{\mathbf{\Sigma }}_{\mathrm{s}}),\end{array}$$

where Σ is the residual variance–covariance matrix. Shiklomanov et al. (2016) assumed Σ was a diagonal matrix with the same residual variance for all elements. For this study, we made two important changes to this methodology. First, to account for the large differences in the range of feasible reflectance values in different wavelength regions (for vegetation, reflectance in the 400–700 nm range is typically much lower than in the 700–1400 nm range), we used a heteroscedastic error model where the residual standard deviation (σ_s) was a linear function of the predicted reflectance R_(pred,s) with slope m and intercept b. Second, to account for autocorrelation in hyperspectral bands, we replaced the diagonal residual covariance matrix with an order-1 autoregressive (AR-1) covariance matrix. Collectively, these two changes produce the following calculation for Σ:

$$\begin{array}{}\text{(45)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{s}}=m{R}_{\mathrm{pred},\mathrm{s}}+b\text{(46)}& {\displaystyle }& {\displaystyle }{\mathbf{\Sigma }}_{\mathrm{s}}={\mathit{\sigma }}_{\mathrm{s}}{\mathit{\varrho }}^{\mathbf{H}}{\mathit{\sigma }}_{\mathrm{s}},\end{array}$$

where H is a matrix describing the distance between bands (0 on the diagonal, increasing regularly toward the corners) and ϱ is the AR-1 autocorrelation parameter. To simplify the inversion procedure, we first performed the inversion using a diagonal covariance matrix (i.e., ϱ=0), then calculated the mean ϱ from the residuals of this fit, and then used this average value (0.700) for our final inversion.

In addition, to mitigate sampling issues related to EDR's saturating response to increasing total LAI (Fig. A1), we added an additional term to our likelihood that assigns a uniform probability distribution over the range 0 to 10 to the EDR-predicted LAI for a given site (LAI_pred,s). In practice, this term causes any parameters resulting in total LAI greater than 10 to be immediately rejected but has no effect on parameters with LAI values less than 10. The maximum value of 10 was selected as a reasonable upper bound on temperate deciduous and evergreen forests in our study region. By comparison, the global maximum of MODIS LAI estimates is between 6 and 7, depending on collection (with most values less than 5; Fang et al., 2012; Yan et al., 2016). A global database of field LAI measurements (Iio et al., 2014) contains values as high as 23.5 for evergreen conifer trees and 12.1 for deciduous broadleaf trees, but these are extreme values, and our maximum of 10 is at least 3 standard deviations away from the mean value for evergreen conifer and deciduous broadleaf trees.

The final expression for the site-specific likelihood is therefore

$$\begin{array}{ll}\text{(47)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{pred},\mathrm{s}},{\mathrm{LAI}}_{\mathrm{pred},\mathrm{s}}=\mathrm{EDR}\left(\mathbf{\Theta }\right|{\mathrm{comp}}_{\mathrm{s}}){\displaystyle }& {\displaystyle }P\left({\mathbf{X}}_{\mathrm{s}}\right|\mathbf{\Theta })=\mathrm{MvNormal}({\mathbf{X}}_{\mathrm{s}}|R_{\mathrm{pred},\mathrm{s}},\mathbf{\Sigma })\\ \text{(48)}& {\displaystyle }& {\displaystyle }\mathrm{Uniform}\left({\mathrm{LAI}}_{\mathrm{pred},\mathrm{s}}\right|\mathrm{0},\mathrm{10}).\end{array}$$

Therefore, our parameter vector Θ consists of the following (summarized in Table 1): 10 EDR parameters per PFT – five parameters for the PROSPECT-5 model (N, Cab, Car, Cw, Cm) and five EDR parameters related to canopy structure (q, χ, SLA, b1Bl, b1Bw) – one parameter per site (relative soil moisture, ψ_s), and the residual slope (m) and intercept (b). With five PFTs and 54 sites, this means that Θ has length $(\mathrm{10}\times \mathrm{5})+\mathrm{54}+\mathrm{2}=\mathrm{106}$.

For priors on the PROSPECT-5 parameters and SLA, we performed a hierarchical multivariate analysis (Shiklomanov et al., 2020b) on PROSPECT-5 parameters and direct SLA measurements from (Shiklomanov, 2018, chap. 3). For priors on the leaf biomass allometry parameters, we fit a multivariate normal distribution to allometry coefficients from Jenkins et al. (2003, 2004) using the PEcAn.allometry package (https://github.com/pecanproject/pecan/tree/develop/modules/allometry, last access: 6 May 2021). For the clumping factor, we used a uniform prior across its full range (0 to 1), and for the leaf orientation factor, we used a weakly informative beta distribution rescaled to the range $(-\mathrm{1},\mathrm{1})$ and centered on 0.5 (Table 1).

To alleviate issues with strong collinearity between the allometry parameters and the specific leaf area, we fixed the allometry exponent parameters (b2Bl and b2Bw) to their prior means for each PFT. Doing so dramatically improved the stability of the inversion algorithm and the accuracy of the results.

We fit this model using the differential evolution with Snooker update (DEzs) Markov chain Monte Carlo (MCMC) sampling algorithm (ter Braak and Vrugt, 2008) as implemented in the R package BayesianTools (Hartig et al., 2019). We ran the algorithm using three independent chains for as many iterations as required to achieve convergence, assessed according to a Gelman–Rubin potential scale reduction factor (PSRF) diagnostic value of less than 1.1 for all parameters (Gelman and Rubin, 1992).

Table 1 EDR parameters and prior distributions.

¹ PFT-specific multivariate normal distribution fit to PROSPECT parameters and SLA from Shiklomanov (2018), chap. 3. ² PFT-specific results from Bayesian fits of allometric equations to allometry data from Jenkins et al. (2003, 2004) using the PEcAn.allometry package.

2.5 Analysis

To assess the extent to which AVIRIS-Classic observations were able to constrain parameter estimates, we compared the prior and posterior distributions for all parameters. To evaluate the performance of the calibrated model, we compared the posterior credible and predicted 95 % intervals of EDR-predicted spectra against the AVIRIS observations at each site. We examined the residuals between EDR-predicted and AVIRIS-observed reflectance across all sites pooled together and evaluated whether residuals varied systematically with site composition or structure by separating sites based on the dominant PFT (calculated as the PFT with the largest ${\sum }_i{\mathrm{DBH}}_i{n}_{\mathrm{plant},i}$ at each site), mean DBH, or mean stem density. We also compared the EDR-predicted LAI against field observations at each site, both across all sites together and within the above site groups based on composition and structure. To evaluate goodness of fit and additive and multiplicative biases, we used an ordinary least-squares regression of mean observed vs. posterior mean predicted LAI.

Figure 3  Marginal prior (pre-calibration; gray) and posterior (post-calibration; black) distributions of PFT-specific parameters related to leaf biochemistry and canopy structure. Distributions are shown as violin plots (rotated and mirrored kernel density plots). PFTs are abbreviated as follows: EH: early hardwood; MH: north mid-hardwood; LH: late hardwood; NP: northern pine; LC: late conifer. Leaf and wood biomass allometry panels are clipped at 0.2 to facilitate differentiation of posterior distributions.

3 Results

Model calibration improved the precision of most PFT-specific parameter estimates, including estimates of leaf parameters whose prior distributions were already independently constrained by an earlier analysis (Fig. 3). Across all PFT-specific parameters, the posterior 95 % credible interval (CI) was, on average, 10 % the size of the prior credible interval. The most constrained parameters on average were EDR canopy structure parameters – namely the wood biomass allometry (<1 % of prior CI), leaf biomass allometry (1 %), leaf orientation factor (8 %), and clumping factor (9 %) – while the least constrained parameters were those related to leaf morphology and biochemistry – namely, effective number of leaf layers (19 %), total chlorophyll content (16 %), total carotenoid content (15 %), specific leaf area (13 %), dry matter content (11 %), and leaf water content (11 %). For PFTs, the largest average relative constraint was for early hardwood (7 %) and the smallest relative constraint was for late hardwood (14 % of prior CI).

For leaf traits, PFT rankings of the posterior estimates largely followed the relative positions of the priors, though there were a few exceptions. In both the prior and posterior, the estimated effective number of leaf mesophyll layers (a.k.a., PROSPECT N parameter) was higher for needleleaf than broadleaf PFTs, with the highest value for northern pine and the lowest value for mid-hardwood. Similarly, specific leaf area (SLA) was lower in conifer than broadleaf PFTs, with the lowest value for late conifer, a higher value for northern pine (despite a similar prior), and higher values still in mid-hardwood and late hardwood and the highest value for early hardwood. Estimated total chlorophyll content (Cab) was similarly high for all hardwood PFTs in both the prior and posterior, but posterior estimates for late conifer and northern pine were lower. Posterior estimates of total carotenoid content (Car) were lower in early and mid-hardwood and northern pine and higher in mid-hardwood and late conifer. Posterior estimates of leaf water content (Cw) were low for early hardwood and northern pine and high for mid- and late hardwood and late conifer; these differences were despite strongly overlapping priors across all PFTs. Posterior estimates of leaf dry matter content (Cm) were lowest for mid- and late hardwood, higher for early hardwood and northern pine, and highest for late conifer, again despite a strongly overlapping prior across all PFTs.

Compared to leaf traits, canopy structural traits had less informative (and PFT-agnostic) priors, and the posterior distributions exhibited some differences across PFTs. Posterior leaf biomass allometry (b1Bl) estimates were lowest in early hardwood and northern pine, higher for late conifer, and highest for late and mid-hardwood. Posterior wood biomass allometry (b1Bw) estimates were lowest for early hardwood and late conifer, slightly higher in mid- and late hardwood, and highest for northern pine. Posterior canopy clumping factor (q) estimates were clustered at or near its upper limit of 1 (i.e., exhibited “edge-hitting behavior”) for early and mid-hardwood, were slightly lower in late hardwood and northern pine, and lowest in late conifer. Posterior leaf orientation factor (χ) estimates were lowest (near zero, indicating randomly distributed leaf angles) for northern pine, higher (more horizontal leaves) for early and mid-hardwoods and late conifer, and highest (at the upper limit of 0.6) for late hardwood. Finally, the calibration was able to constrain site-specific soil optical properties across all sites (Fig. A4).

Figure 4  Differences between AVIRIS-observed and EDR posterior predictive mean surface reflectance by site. Each thin gray line is a site-specific AVIRIS observation. Top left panel shows results across all sites, and remaining panels group sites according to dominant PFT. Within each panel, blue shading shows the 25 %–75 % quantile range and the thick black line is the median by wavelength for that specific site.

Figure 5  Differences between AVIRIS-observed and EDR posterior predictive mean surface reflectance by site, averaged across wavelength regions. Sites are grouped by dominant PFT, as in Fig. 4. Note the differences in the y-axis scale across panels.

Figure 6  Left: comparison between AVIRIS-observed (black) and EDR-predicted (mean prediction in green; 95 % posterior predictive interval in gray) surface reflectance for a geographically (Fig. 1), compositionally, and structurally (Fig. 2) representative sample of sites used in the calibration. Right: histogram of stem DBH by PFT at the corresponding site.

The accuracy and precision of EDR simulated spectra relative to AVIRIS observations varied across sites (Figs. 4, 5, 6, and A5). On average, EDR tended to accurately (within 0.01) reproduce reflectance in the 400–750 nm range, underpredict AVIRIS reflectance by ∼ 0.02 in the 750–1100 nm range, and overpredict AVIRIS reflectance by ∼ 0.03 in the 1100–1300 nm range. However, only the latter behavior was consistent across sites; below 1100 nm, sites dominated by any PFT could have low, accurate, or high estimates relative to the AVIRIS observations. The only consistent bias (expressed as interquartile range in bias not overlapping zero) we observed with respect to PFT was an underestimate of reflectance in the 750–1100 nm range for northern pine; otherwise, we did not observe any consistent patterns in mismatch between AVIRIS-observed and EDR-predicted reflectance with respect to tree size, stem density, or composition (Figs. 5 and A6–A15). The EDR posterior predictive interval overlapped AVIRIS observations in all but a few individual cases (Fig. A5), suggesting that our estimates of model uncertainty are realistic.

Figure 7 EDR predictions of site-specific true LAI compared to observed values. Horizontal error bars are posterior 95 % predictive intervals. Vertical error bars are mean ±1 standard deviation of the observed values. Dashed line shows the 1:1 relationship, and solid line is a least-squares predicted vs. observed regression with the equation marked in the upper left corner. Points are colored according to dominant PFT, calculated as in Fig. 4.

Figure 8 Bias in LAI predictions from calibrated EDR relative to observations, as a function of site mean DBH. Sites are grouped according to dominant PFT, same as in Fig. 4.

Figure 9 Same as above but as a function of mean site stem density.

Leaf area index predicted from calibrated EDR parameters captured 46 % of the variability in the observations (Fig. 7). The observed vs. predicted line had a slope of 0.35 and an intercept of 3.09, indicating that EDR calibration underpredicted true LAI on average but exaggerated LAI variability across sites. In general, EDR tended to underpredict LAI at high-density sites with low mean DBH and overpredict LAI at low-density sites with high mean DBH (Figs. 8 and 9). The trend with mean DBH was generally true across all PFTs but was most pronounced for early-hardwood- and late-conifer-dominated sites (Fig. 8), while the trend with stand density was most pronounced for late-conifer-dominated sites (Fig. 9).

Figure 10  Comparison of EDR and PRO4SAIL (labeled as “SAIL” for conciseness) predictions of reflectance for identical, single-cohort canopies as a function of solar zenith angle (θ_s). These simulations use identical PROSPECT and canopy structure parameters, a nadir-viewing sensor, and LAI of 3. For EDR, we use a single cohort with LAI of 3 prescribed directly. “SAIL: HR” is the “hemispherical reflectance” (“blue-sky albedo”), calculated (as in EDR) as the average of directional–hemispherical reflectance (DHR) and BHR streams weighted by the direct sky fraction (0.9; same value used for EDR). Similarly, “SAIL: DR” is the “directional reflectance”, calculated analogously from the bi-directional reflectance (BDR) and hemispherical–directional reflectance (HDR) streams. Individual SAIL fluxes are shown with dotted lines.

For identical canopies, EDR consistently predicted lower hemispherical reflectance than PRO4SAIL (Fig. 10). This difference was most pronounced when the Sun was directly overhead (${\mathit{\theta }}_{\mathrm{s}}={\mathrm{0}}^{\circ }$; cos (θ_s)=1) and declined with increasing solar zenith angle. For solar zenith angles typical of our study (${\mathit{\theta }}_{\mathrm{s}}\approx {\mathrm{30}}^{\circ }$; cos (θ_s)=0.85), EDR hemispherical reflectance predictions were very close to PRO4SAIL directional reflectance predictions over a wide range of LAI values (Fig. A16).

4 Discussion

Calibrating and validating vegetation models using optical remote sensing data have typically involved derived data products (e.g., MODIS gross primary production, GPP) that rely on their own models, in other words, “bringing the observations closer to the models”. In this study, we presented an alternative approach whereby we “bring the models closer to the observations” by training a vegetation model to simulate full-range hyperspectral surface reflectance that is closer to the measurements made by optical remote sensing instruments. We argue this is a more generalized approach, as many dynamic vegetation models already contain their own internal representations of canopy radiative transfer and thus can be modified to provide outputs that can mimic optical remote sensing (i.e., can be used as “satellite simulators” to connect underlying model state to emergent reflectance). We then demonstrated how this approach could be used to calibrate the model against airborne imaging spectroscopy data from AVIRIS-Classic. We found that such calibration reduced uncertainties in parameters related to leaf biochemistry and canopy structure, even for parameters with well-informed priors (Fig. 3). Moreover, we found that the calibrated model was able to reproduce observed surface albedo (Figs. 4, 5, 6, and A5) reasonably well across a large number of geographically (Fig. 1), structurally, and compositionally (Fig. 2) diverse sites. However, the calibrated model underpredicted LAI at sites with mostly small trees and overpredicted LAI at sites with mostly large trees (Figs. 7 and 8).

In this study, the vegetation composition at each site (including the PFT distribution and size–age structure) was prescribed in detail based on inventory data. This allowed us to focus the calibration on model-parameter-related canopy radiative transfer model parameters. However, ED2 is a dynamic vegetation model whose core purpose is to predict how vegetation composition and structure evolve through time. An important future direction of this work is to evaluate such dynamic ED2 simulations where vegetation composition and structure are predicted with some uncertainty. In ED2, shortwave canopy radiative transfer is already linked (through shared parameters and state variables) to other important model processes, including thermal radiative transfer, micrometeorology, photosynthesis, respiration, and competition (Longo et al., 2019), and therefore changes in canopy radiative transfer parameters have profound consequences for ED2 predictions of ecosystem fluxes and composition (Viskari et al., 2019). Future work could further strengthen this link by embedding the PROSPECT coupling demonstrated in this study into ED2 itself, replacing ED2's currently prescribed leaf optical properties with simulated optical properties that change with leaf morphology and biochemistry. For example, the PROSPECT leaf water content parameter (Cw) provides a physical link between leaf optical properties and hydraulics, so such a configuration could allow surface reflectance information to constrain ED2's recently developed dynamic hydraulics module (Xu et al., 2021).

With 54 sites in our calibration, any single site represents <2 % of the data, and for a joint calibration without site random effects, we have every reason to believe that the calibration is not overfitting to any individual site. Trying to fit any one site well would cause others to do worse (especially given the large observed variability in forest structure; Fig. 2) unless the EDR model structure was reasonable and the parameters chosen were genuinely good choices. We therefore did not perform an in-sample cross validation, as we believed the benefit of doing so would be low relative to the high computational cost. That said, we recognize that out-of-sample validation is a useful test of model performance, and we recommend performing out-of-sample validation in similar studies in the future.

The canopy radiative transfer model in ED2 is derived from the two-stream model of Sellers (1985) and adapted to a multi-level canopy. Similar versions of this two-stream formulation are present in other land surface models, including CLM (Oleson et al., 2013), SiB (Baker et al., 2008), Noah (Niu et al., 2011), tRIBS-VEGGIE (Ivanov et al., 2008), IMOGEN (Huntingford et al., 2008), and JULES (Best et al., 2011). Although the exact parameterization and implementation differs across these models, the similarity of the underlying conceptual framework and radiative transfer coefficients means that our approach should be directly transferable to all of these models.

Nevertheless, our analysis echoed some known challenges in canopy radiative transfer modeling. One challenge is equifinality in the contributions of leaf biochemistry, leaf morphology, and different aspects of canopy structure to canopy albedo, which means that multiple variable and parameter combinations can produce very similar canopy albedo responses (Lewis and Disney, 2007; Figs. A1–A3). We mitigated the equifinality between leaf traits and canopy structure by using informative priors on leaf traits from an independent data source (Shiklomanov, 2018). However, there is additional equifinality in the effects of the EDR canopy structure parameters. For example, because the effective LAI used in EDR's actual radiative transfer calculations is defined as the product of “true” LAI and clumping factor (Eq. 20), and because LAI is, in turn, derived from multiple parameters (leaf biomass allometry, specific leaf area; Eq. 21), these parameters collectively cannot be independently determined from reflectance data alone. At the same time, increasing the leaf orientation factor (more horizontal, or “planophile”, leaf orientation) has a similar (although not identical) effect to increasing LAI and clumping factor – namely, increasing canopy reflectance, especially in the near-infrared spectral region (Fig. A3). Collectively, these issues may help explain some of the edge-hitting behavior (parameter distributions clustered at the ends of the distribution) observed in our posterior estimates (Fig. 3) and some of the bias in our LAI estimates (Fig. 7). In future work, we suggest combining our approach with additional kinds of remote sensing measurements capable of directly constraining these structural parameters such as waveform lidar (which can provide a robust constraint on the canopy structural profile; Ferraz et al., 2020) to reduce equifinality.

That being said, one major advantage of the Bayesian calibration approach is that its output is a joint posterior distribution that includes not only fully quantified uncertainties for each parameter but also the variance–covariance matrix across the full set of parameters. Equifinality in parameters would manifest as strong pairwise correlation between parameters in the posterior distribution. Examining this correlation matrix shows that there are some parameter pairs with strong correlations, such as the hypothesized negative correlations between leaf allometry and clumping factor across some PFTs (Fig. A17). However, these correlations do not occur in all parameters that exhibited edge-hitting behavior. For instance, clumping factor exhibited edge-hitting behavior only for early- and mid-successional hardwood PFTs (Fig. 3), but the corresponding correlation coefficients were positive and near zero, respectively, while strong negative correlations for the other PFTs did not result in edge-hitting posteriors (Fig. A17). Similarly, the edge-hitting leaf orientation factor posterior for late hardwood (Fig. 3) had near-zero (or, in one case, positive) correlations with all other parameters (Fig. A17). Strong correlations also occurred among some of the PROSPECT parameters and between PROSPECT and structural parameters but contributed little to equifinality because the strong constraints on PROSPECT led to overall small covariance terms (results not shown). Finally, because our calibration captured all of these covariances, the presence of moderate equifinality did not preclude ecologically meaningful parameter constraints or accurate predictions because these covariances are being propagated into predictions. This is directly analogous to how a linear regression can have a tight confidence interval, despite high correlations between the slope and intercept, with that equifinality driving the characteristic hourglass shape of a regression confidence interval.

EDR tended to underpredict LAI at high-density sites with low mean DBH and overpredict LAI at low-density sites with high mean DBH (Figs. 7, 8, and 9). The relationship between DBH and LAI is controlled primarily by the leaf biomass allometry, which in EDR is fixed at the PFT level (Eq. 22). This fixed relationship neglects the known inter- and intra-specific variability in tree allocation strategies (Forrester et al., 2017; Dolezal et al., 2021). For example, Forrester et al. (2017) showed that the relationship between DBH and foliar biomass is modulated by tree age, stand density, and climate variables, none of which are accounted for in the ED2 allometry routines. This variability can be incorporated directly into ED2 by making allometry parameters dynamic functions of some of the aforementioned covariates or indirectly via hierarchical calibration whereby model parameters vary across sites (Dietze et al., 2008). Overall, this analysis reiterates the importance of evaluating models against multiple distinct variables – after all, none of these biases would have been apparent from looking at the reflectance simulations alone.

Our analysis also revealed some structural issues with EDR itself. EDR consistently predicted lower hemispherical reflectance than SAIL (Figs. 10 and A16). This difference can be attributed primarily to differences in how each model defines the direct radiation backscatter coefficient in the radiative transfer equation. A detailed description of the discrepancy is provided in Yuan et al. (2017). Briefly, EDR (and the Sellers, 1985 model from which EDR is derived) defines direct radiation backscatter as a function of the single-scattering albedo (Eq. 24), which in turn is a challenging integral involving the leaf scattering phase function and leaf-projected area function (Eq. 25). The analytical solution to this integral in EDR (Eq. 28) assumes a uniform scattering phase function, which is appropriate for point scatterers but less so for horizontal surfaces like leaves. The practical consequence of this assumption is a lower value of the direct radiation backscatter and therefore a lower albedo, which is consistent with the results of our sensitivity analysis. This underestimation of albedo described above may also help explain the edge-hitting behavior in our posterior distributions (Fig. 3) as well as the relatively low accuracy of our LAI estimates (Fig. 7). Specifically, our EDR calibration may be trying to compensate for underestimated albedos via a tendency to prefer higher effective LAI values (which results in higher values of the leaf biomass allometry and clumping factor for some PFTs; e.g., early and mid-hardwood and northern pine in Fig. 3) and more horizontal leaf distributions (i.e., higher leaf orientation factor; e.g., late hardwood in Fig. 3), both of which increase albedo (Figs. A1–A3).

Meanwhile, the SAIL definition of direct backscatter is a more simple function of leaf scattering, leaf angle distribution, and canopy optical depth that also produces a more accurate albedo estimate (Yuan et al., 2017). Given that underestimating albedo can have significant consequences for the biological, ecological, and physical predictions of ED2 (Viskari et al., 2019), incorporating this fix into the ED2 canopy RTM is an important future direction of our work. However, doing so is outside the scope of this work because it would require propagating the different coefficients through the complex, multiple-canopy-layer solution of EDR (Sect. 2.1.2) – a non-trivial task.

A significant caveat to the broader application of our approach is that there is a subtle but significant difference between the physical quantity EDR predicts and the quantities typically measured by optical remote sensing. Specifically, EDR predicts the hemispherical reflectance – the ratio of total radiation leaving the surface to the total radiation incident upon the surface, integrated over all viewing angles (also known as “blue-sky albedo”). On the other hand, optical remote sensing platforms typically measure the directional reflectance factor – the ratio of actual radiation reflected by a surface to the radiation reflected from an ideal diffuse scattering source (e.g., a Spectralon calibration panel) subject to the same illumination, in a specific viewing direction (Schaepman-Strub et al., 2006). These two quantities are numerically equivalent only for ideal Lambertian surfaces or, for non-Lambertian surfaces, under specific Sun-sensor geometries. However, vegetation canopies – the focus of this study – are known to exhibit reflectance with very strong angular dependence, so a comparison of canopy hemispherical reflectance with a remotely sensed directional reflectance factor is invalid.

Our specific analysis is valid because – as described in the Methods section (Sect. 2.3) – we used AVIRIS data that were also BRDF corrected, whereby the directional reflectance estimates from the atmospheric correction process were further converted to estimates of hemispherical reflectance via a polynomial approximation of the Ross–Li semi-empirical BRDF model (Lucht et al., 2000). Another dataset that would have been valid for our analysis (albeit, one with much lower spatial and spectral resolution) is the MODIS albedo product (MOD43), which takes advantage of the angular sampling of the MODIS instrument to quantify the surface BRDF characteristics and therefore more precisely estimate the albedo (Wang et al., 2004; Schaaf and Wang, 2015). However, our approach as described here would not be valid for the surface reflectance products produced by nadir-viewing instruments such as Landsat, Sentinel-2, or most airborne platforms, at least without additional processing steps on the data, or, preferably, modification of the underlying radiative transfer models to allow for the prediction of directional reflectance. Fortunately, the assumptions and parameters that comprise two-stream radiative transfer models like EDR and its parent model (Sellers, 1985) are readily adaptable to prediction of directional reflectance. For example, the SAIL model (Verhoef, 1984; Verhoef and Bach, 2007) – which predicts both hemispherical and directional reflectance, and which has a long history of successful application to remote sensing – makes the same general assumptions as EDR and even shares many of the underlying coefficients (Yuan et al., 2017). Alternatively, land surface models can take advantage of recent advances in radiative transfer theory to improve their accuracy without significant computational penalty (e.g., Hogan et al., 2018).

A related issue is the missing or simplistic treatment of two- and three-dimensional heterogeneity in canopy structure in EDR. For one, the treatment of leaves as infinitely small elements randomly distributed through the canopy space – a common feature of all two-stream approximations – neglects complex realities of the canopy light environment such as gaps and self-shading. In EDR, self-shading is handled via the clumping factor parameter, which functions as a scalar correction on the leaf area index (Eq. 20). A key feature of EDR design is its representation of multiple co-existing plant cohorts competing for light within a single patch; however, horizontal heterogeneity and competition between these cohorts is ignored. Improved representation of lateral energy transfer can improve the accuracy of simulations of the canopy light environment, and recent theoretical advances show that this can be accomplished without a significant loss in computational performance (Hogan et al., 2018). Treatment of horizontal competition also plays an important role in the outcomes of competition for light between different plants (Fisher et al., 2018). A useful avenue for development and parameterization of these models is comparison to more sophisticated and realistic three-dimensional representations of radiative transfer (e.g., Widlowski et al., 2007), which are themselves too computationally demanding to be coupled to ecosystem models, but from which empirical distributions and response functions could be derived and against which the behavior of simpler models could be evaluated.

5 Conclusions

Remote sensing observations are unrivaled in their spatial completeness and extent, notably extending to regions like the tropics and high latitudes that are relatively undersampled but have a disproportionate impact on the global climate system (Schimel et al., 2015) and/or global biodiversity (Jetz et al., 2016). At the same time, satellite time series provide multidecadal records with relatively high temporal frequency, which have tremendous utility for calibrating model projections of past ecological dynamics (Kennedy et al., 2014; Pasquarella et al., 2016). Used in combination with other emerging data sources, including global trait databases and eddy covariance measurements, remote sensing can be a transformative force in ecosystem ecology.

In this paper, we showed that using a vegetation model to directly simulate surface reflectance is a promising approach for calibrating and validating models against remotely sensed observations. To do this, we modified the ED2 dynamic vegetation model to predict full-range hyperspectral hemispherical surface reflectance and then calibrated this modified model against airborne, BRDF-corrected imaging spectroscopy data. The calibration successfully reduced uncertainties in model parameters related to canopy structure and leaf biogeochemistry for five plant functional types characteristic of temperate forests of the northeastern United States. The calibrated model was able to accurately reproduce observed surface reflectance for sites with highly varied forest composition and structure using a single common set of parameters (i.e., not site-specific parameters). However, the calibrated model predicted leaf area index values that did not agree well with observations and had parameter estimates that exhibited edge-hitting behavior, both of which suggest structural issues in the model. Comparison against a canopy radiative transfer model commonly used in the remote sensing community (PRO4SAIL; Verhoef and Bach, 2007) suggested that our model may be systematically underpredicting surface albedo. Given the direct role albedo plays in the canopy light and thermal environment in ED2, this bias could have significant downstream consequences for ED2 predictions of physiological and ecological processes. We therefore recommend structural changes to the ED2 canopy radiative transfer model to resolve this bias and recommend calibrating the updated model against remotely sensed surface reflectance, as we demonstrated here. We note that the basic structure and assumptions of the ED2 canopy radiative transfer scheme are shared by many other vegetation models, so we expect that both this issue and our recommendations for resolving it are highly transferable within the vegetation modeling community. More generally, we recommend the development of additional “observation operators” similar to ours for other classes of remote sensing data, such as thermal, microwave, and lidar, in ED2 and other dynamic vegetation models, to allow these models to take full advantage of remote sensing observations.

Appendix A: Supplementary figures

Figure A1  Sensitivity of EDR-predicted hemispherical reflectance to true LAI. These simulations assume a single-cohort canopy with effective number of mesophyll layers N=1.4, total chlorophyll content Cab =40, total carotenoid content Car =10, leaf water content Cw =0.01, leaf dry matter content Cm =0.01, clumping factor q=1, leaf orientation factor χ=0, cos (θ_s)=0.85, and soil moisture fraction ψ=0.5.

Figure A2  Same as above but with LAI fixed to 3 and varying clumping factor (q).

Figure A3  Same as above but instead varying leaf orientation factor (χ).

Figure A4  Site-specific relative soil moisture (0 indicates dry; 1 indicates wet) posterior estimates. Sites are sorted in order of increasing weighted evergreen fraction.

Figure A5  Comparison between AVIRIS-observed (black) and surface reflectance for each site used in the calibration. Sites are sorted in order of decreasing mean difference between observed and EDR-predicted reflectance (largest underestimates first; largest overestimates last).

Figure A6  Mean reflectance bias (EDR predicted − observed) for each by spectral region and dominant PFT as a function of site stem density. PFTs are abbreviated as follows: EH: early hardwood; MH: north mid-hardwood; LH: late hardwood; NP: northern pine; LC: late conifer.

Figure A7  EDR-predicted vs. observed spectra and species composition for the first quartile of sites by DBH.

Figure A8  As above but for the second quartile of sites by DBH.

Figure A9  As above but for the third quartile of sites by DBH.

Figure A10  As above but for the fourth quartile of sites by DBH.

Figure A11  As above but for sites where early hardwood trees had the largest mean DBH.

Figure A12  As above but for sites where mid-hardwood trees had the largest mean DBH.

Figure A13  As above but for sites where late hardwood trees had the largest mean DBH.

Figure A14  As above but for sites where pine trees had the largest mean DBH.

Figure A15  As above but for sites where late conifer trees had the largest mean DBH.

Figure A16  Same as Fig. 10 but varying LAI and fixing cos (θ_s)=0.85, a typical value for our study.

Figure A17  Posterior correlation matrix for PFT-specific parameters. Note that only correlations among parameters within the same PFT are shown – the full 106 × 106 dimensional correlation matrix is far too large to display in this format.