Allometric Scaling and Resource Limitations Model of Tree Heights: Part 1. Model Optimization and Testing over Continental USA

Abstract

A methodology to generate spatially continuous fields of tree heights with an optimized Allometric Scaling and Resource Limitations (ASRL) model is reported in this first of a multi-part series of articles. Model optimization is performed with the Geoscience Laser Altimeter System (GLAS) waveform data. This methodology is demonstrated by mapping tree heights over forested lands in the continental USA (CONUS) at 1 km spatial resolution. The study area is divided into 841 eco-climatic zones based on three forest types, annual total precipitation classes (30 mm intervals) and annual average temperature classes (2 °C intervals). Three model parameters (area of single leaf, α, exponent for canopy radius, η, and root absorption efficiency, γ) were selected for optimization, that is, to minimize the difference between actual and potential tree heights in each of the eco-climatic zones over the CONUS. Tree heights predicted by the optimized model were evaluated against GLAS heights using a two-fold cross validation approach (R² = 0.59; RMSE = 3.31 m). Comparison at the pixel level between GLAS heights (mean = 30.6 m; standard deviation = 10.7) and model predictions (mean = 30.8 m; std. = 8.4) were also performed. Further, the model predictions were compared to existing satellite-based forest height maps. The optimized ASRL model satisfactorily reproduced the pattern of tree heights over the CONUS. Subsequent articles in this series will document further improvements with the ultimate goal of mapping tree heights and forest biomass globally.

1. Introduction

Several recent articles have reported generating spatially continuous maps of forest canopy heights and/or biomass using a combination of remote sensing data, in-situ measurements and non-physical/non-physiological or statistical scaling approaches (e.g., [1–8]). Tree height estimation, and potentially biomass, is now possible with altimeter data from terrestrial, airborne, and satellite lidar (e.g., [1–4,9–13]). Lidar waveform data from the Geoscience Laser Altimeter System (GLAS) instrument onboard the Ice, Cloud and land Elevation Satellite (ICESat) have been used to map global and regional forest heights [1,2] and live aboveground biomass [3,4,13]. However, discrete distributions of tree heights retrieved from GLAS data should be extrapolated to generate continuous maps of forest heights or biomass [1,2,4]. This “black-box” type of extrapolation has the obvious limitation that it is often done using non-physical/non-physiological procedures in conjunction with spatially continuous remote sensing and climate data.

Physical/physiological models for mapping tree heights or biomass rely on mechanisms governing plant growth. The Allometric Scaling and Resource Limitations (ASRL) model [14] is one such physical/physiological model. This predicts local maximum tree heights. The ASRL model integrates allometric scaling laws of trees and energy budgets limited by local resources such as water, air temperature, sunlight, and wind [14]. Some researchers (e.g., [15–18]) however doubt the relevance of plant allometric scaling laws given the high variability observed in actual forests [19,20]. Other studies have demonstrated the applicability of scaling laws for quantifying forest structure and dynamics (e.g., [21,22]) and estimating live biomass in forest stands (e.g., [23,24]). The ASRL model implements the steady-state allometric approach based on the assumption that physiological traits of trees generally follow allometric scaling rules [14]. Nevertheless, the allometric coefficients and scaling exponents of the ASRL model are assumed constant across different eco-climatic zones and forest types of varying age classes. This often results in disparities between measurements and model predictions. Here, a significant progress in mapping tree heights and biomass is possible if the power of allometric scaling laws, local energy budgets and resource limitations can be incorporated with the advancements of remote sensing altimetry (i.e., GLAS data) for scaling purposes.

Generating continuous fields of tree heights and biomass is the larger objective of this multi-part series of articles. In this first article, we focus on how the ASRL model can be used with GLAS data to map actual tree heights over the continental USA (CONUS) at 1 km spatial resolution. The ASRL model is briefly explained in Section 2 together with key equations and parameters. Section 3 includes descriptions of input data for ASRL model and GLAS data preprocessing. Information of the model optimization and evaluation is provided in Section 4 followed by results and discussion (Section 5) and concluding remarks (Section 6). The second paper of this series [25] examines in detail how the same procedures work at a local scale, specifically at several FLUXNET sites. Future articles in the series will consolidate these results and extend them to biomass estimation.

2. The ASRL Model

The ASRL model [14] predicts potential maximum tree heights using a combination of allometric scaling laws and energy budgets constrained by local resource limitations, such as water, air temperature, solar radiation, and wind. The model incorporates estimates of parameters related to tree geometry (e.g., canopy radius and leaf area), light (e.g., soil reflectance, leaf absorptivity and deep canopy reflection coefficient), and water flow (e.g., root absorption coefficient, depth of a stomata, and an exponent for metabolism).

In the fundamental premises of the ASRL model, a tree obtains sufficient resources (water and nutrients) to meet its needs for the growth and the availability of local resources limits the maximum potential growth. This is expressed by an inequality equation of basal metabolic rates (Q_p ≥ Q_e ≥ Q₀), where Q_p is the available flow rate, Q_e refers to the potential evaporative flow rate, and Q₀ corresponds to the required flow rate of resources in a tree [14]. Q₀ is solely determined by allometric scaling rules, while Q_p and Q_e are additionally associated with local resources, such as water, air temperature, solar radiation, and wind. The maximum tree growth can be calculated given Q_p and Q_e. Key equations of the ASRL model are given in

Table 1. Key equations of the Allometric Scaling and Resource Limitations (ASRL) model [14]. Underlined variables (α, η, and γ) are selected parameters in the ASRL model optimization, further explained in Section 4.3.

Categories	Variables	Symbols	Key Equations	Sub-Variables
Equations of Basal Metabolic Rates	Available Flow Rate	Qp	Qp = γ π rroot2Pinc	γ = Root Absorption Efficiency; rroot = Radial Extent of Root System; Pinc = Incoming Precipitation Rate
	Evaporative Flow Rate	Qe	Qe = af Ecan μw ρw−1	af = Effective Area over the Latent Heat Flux Loss; Ecan = Evaporative Flux of Canopy; μw = Molar Mass of Water (= 1.80 × 10−2 kg·mol−1); ρw = Density of Water (= 1.0 × 103 kg·m−3)
	Required Flow Rate	Q0	Q0 = β2hη2	β2 = Proportionality Constant for Metabolism (≈9.2 × 10−7 L·day−1 cm·−η2); h = tree height; η2 = Exponent for Metabolism (≈2.7)
Sub-equations of Evaporative Flow Rate	Effective Area over the Latent Heat Flux Loss	af	af = 2 aL δs as	aL = Total One-sided Area of All Leaves on a Tree; δs = Density of Stomata on a Leaf (=220 stomata mm−2); as = Area of a Single Stomata (=235.1 μm2)
	Total One-sided Area of All Leaves on a Tree	aL	aL = α nN	α = Area of Single Leaf; n = Branching Parameter (=2); N = Number of Branching Generations
	Number of Branching Generations	N	N = 2 ln (r0/rN)/ln n	r0 = Maximum Stem Radius; rN = Radius of Terminal Branch (=0.4 mm)
	Evaporative Flux of Canopy	Ecan	Refers to [14]	Equation uses Rate of Absorbed Solar Radiation (Rabs) along with Canopy Radius (rcan) and Area of Single Leaf (α)
Sub-allometric Scaling Equations	Radial Extent of Root System	rroot	rroot = β31/4h	β3 = Root to Stem Mass Proportionality (≈0.423)
	Maximum Stem Radius	r0	r0 = 0.5 (β2/β1)1/η1hη2/η1	β1 = Proportionality Constant for Metabolism (=0.257 L·day−1 cm−η1); η1 = Exponent for Metabolism (=1.8)
	Canopy Radius	rcan	rcan = β5hη	β5 = Proportionality Constant for Canopy Radius (=35.24 cm·m−η); η = Exponent for Canopy Radius

[14].

The maximum tree growth varies depending on many factors (e.g., climatic and soil condition, forest types and stand ages), but the ASRL model implements consistent allometric scaling parameters and exponents across different eco-climatic regimes and forest types of varying age classes [14]. In this study, we test where the ASRL model prediction successes and fails. Our optimization process adjusts several allometric parameters to minimize the difference between actual observations and the model predictions. More details are explained in Section 4.3.

3. Data

3.1. Input Data for the ASRL Model

The key input climatic variables include annual total precipitation, annual average temperature, annual incoming solar radiation, annual average wind speed and annual average relative humidity. Additionally, two categories of ancillary input data are needed: (a) Digital Elevation (DEM) and Leaf Area Index (LAI) for initializing the ASRL model and (b) land cover and tree cover for delineating forested lands.

Table 2. Climatic and other ancillary variables used for ASRL model simulations.

Types	Required Input Variables	Units	Temporal Range	Spatial Resolution	Used Data Sets
Climatic Variables	Annual Total Precipitation	mm	1980–1997	1 km	DAYMET model [28] Annual Average Relative Humidity (%) was computed by the formula provided by World Meteorological Organization (WMO) [29]
	Annual Average Temperature	°C	1980–1997	1 km
	Annual Incoming Solar Radiation	W/m2	1980–1997	1 km
	Annual Average Vapor Pressure	hPa	1980–1997	1 km
	Annual Average Wind Speed	m/s	2000–2008	32 km	North American Regional Reanalysis (NARR) data [30]
Ancillary Variables I	Digital Elevation (DEM)	m	2009	30 m	National Elevation Dataset (NED) [31]
	Growing Season Average Leaf Area Index (LAI)	N/A	2003–2006 Jun–Sep	1 km	Post-processed Moderate Resolution Imaging Spectroradiometer (MODIS) LAI products [32]
Ancillary Variables II	Land cover	N/A	2006	30 m	National Land Cover Database (NLCD) [33]
	Percentage of Tree Cover	%	2005	250 m	MODIS Vegetation Continuous Fields (VCF) Collection 5 [34]

lists input data (climate data and ancillary inputs) required for the ASRL model. Finer gridded data (e.g., 30 m or 250 m) were resampled to 1 km resolution in this study using the majority principle for categorical values and cubic convolution for numerical values [26]. Wind speed data at 32 km resolution were spatially interpolated to 1 km resolution using an Inversed Distance Weighting (IDW) method [27].

3.1.1. Climate Data

The ASRL model was mostly driven by climatic variables derived from the DAYMET model [28]. The DAYMET model uses daily weather observations (1980–1997) to produce wall-to-wall climate grids of annual total precipitation, annual average temperature, annual incoming solar radiation and annual average vapor pressure over CONUS. Annual relative humidity, an input parameter of the ASRL model, was derived from annual average temperature and annual vapor pressure using functional relationships provided by the World Meteorological Organization (WMO) [29] (Section S1.1 of Supplementary Information). Wind speed was derived from the North American Regional Reanalysis (NARR) data [30]. The NARR provides monthly mean values from 1979 till present at a spatial resolution of 0.3-degree (∼32 km). Monthly mean values from years 2000 to 2008 were averaged to obtain annual average values and used as input to the ASRL model.

3.1.2. Ancillary Data

The first set of ancillary data (DEM and LAI) is used for the initial ASRL predictions of potential tree heights. The growing season (June to September) average LAI data were calculated from a refined version of the standard Moderate Resolution Imaging Spectroradiometer (MODIS) LAI products (1 km grids) for the time period from 2003 to 2006 [32] (Section S1.2).

The second set of ancillary data is required during the ASRL model prediction and the parametric optimization to identify forested lands. The model simulations were conducted on spatial regions categorized into three forest classes—deciduous, evergreen, and mixed forests—each with percent tree cover ≥50 percent based on the MODIS Vegetation Continuous Field (VCF) product (Figure 1).

3.2. GLAS Tree Heights

The GLAS laser altimetry data provide information related to land elevation and vegetation height at a spatial resolution of ∼70 m (ellipsoidal footprints) and at ∼170 m spaced intervals [35,36]. The latest release (Release-33) of the standard GLAS product corresponding to the GLAS Level-2 Land Surface Altimetry (GLA14; L2 Land Surface Altimetry) for the period 2003 to 2006 was obtained from the National Snow and Ice Data Center (NSIDC) for this study. The GLA14 product was used to estimate forest canopy heights within each footprint (e.g., [1,3]) using geolocation information and waveform parameters, such as signal beginning and echo energy peaks [37]. There are notable heterogeneities in the dimension and shape of the individual GLAS footprints. The GLAS instrument was designed to have a fixed footprint size, but the dimensions of footprints are significantly changed depending on the laser periods (e.g., orbits and spans of campaigns) [38]. To simplify, we assumed that all GLAS footprints have a circular diameter of 70 m [39]. Data from May to October of each year were used, as these come from the growing season and correspond to the MODIS LAI product.

GLAS waveform data are affected by three degrading factors: (a) atmospheric forward scattering and signal saturation, (b) background noise (low cloud) and (c) slope gradient effects. Additionally, GLAS footprints over non-forest and/or bare ground must be filtered from analysis. Five screening steps were applied to remove invalid GLAS waveform data prior to retrieval of tree heights (

Table 3. Five screening steps to remove invalid Geoscience Laser Altimeter System (GLAS) footprints over the CONUS. Final valid GLAS footprints are 126693 in this study.

Screening Steps	Description	Number of Valid GLAS Footprints	References
1. Atmospheric Forward Scattering and Signal Saturation Filter	- Cloud-free and saturation-free GLAS waveform data; - Internal flag of GLAS data “FRir_qaFlag = 15” and “satNdx = 0”	1,822,739	[3]
2. NLCD and VCF Filters	- GLAS footprints over forested lands; - Geolocation of NLCD and VCF pixels (pixels nearest to the center of a GLAS footprint); - Deciduous, evergreen, and mixed forests with greater than 50% of the tree cover	1,659,061	-
3. Background Noise Level (Low Cloud) Correction Filter	- No background noise level in GLAS waveform data; - Absolute difference (≤ 50 m) between the NED DEM and the internal elevation (“i_elev”) of GLAS waveform data	161,533	[36,40]
4. Slope Gradient Correction Filter	- GLAS footprint over non- high topographic condition; - Slope value < 20 ° of the nearest pixel from GLAS data; - Additionally correction of the potential bias (= footprint size × tan (slope))	129,705	[25,40,41]
5. Removal of Remaining Outliers	- Using two standard deviations from the mean of GLAS tree heights (5 m < HGLAS ≤ 100 m)	126,693 (Final)	-

). Note that we removed any remaining outliers using two standard deviations from the mean of GLAS tree heights (5 m < H_GLAS ≤ 100 m) in this analysis. Final valid GLAS footprints were intersected with the pixels (=1 km) over forested lands. We averaged tree heights derived from the GLAS footprints falling in a pixel. This generates a raster distribution of GLAS heights (Figure S1).

There are two approaches of retrieving tree heights from GLAS waveform data [40]: (a) the “statistical analysis for examining full GLAS waveform extents” [41–44] and (b) the “decomposition of GLAS waveforms into multiple Gaussian distribution curves” [1,3,6,40,45–47]. We considered only the second approach in this study.

We used two standard altimetry variables of GLA14 product (signal begin range increment, SigEndOff; centroid range increment for the last Gaussian Peak, gpCntRngOff 1). Theoretically, gpCntRngOff 1 is assumed to represent the ground level elevation within a GLAS field-of-view, while SigBegOff refers to the highest point of a surface. Amongst five possible GLAS height metrics representative of tree heights [25], we used the best metric that closely resembled field-measured and Laser Vegetation Imaging Sensor tree heights (R² = 0.70; RMSE = 4.42 m; [25]): that is, “SigBegOff – gpCntRngOff 1” with correction of the potential bias [40] (Equation (1)).

$$H_{\mathit{GLAS}}=(\mathit{SigBegOf}f-\mathit{gpCntRngOff})-\frac{d\times \text{tan}\theta }{2}$$

(1)

Here, SigBegOff is the signal begin range increment and gpCntRngOff 1 refers to the last peak of Gaussian of GLAS waveform, d is the footprint diameter of GLAS data (∼70 m) and θ is the slope value nearest to the geolocation of GLAS footprint center.

4. Methods

4.1. Defining Climatic Zones

The forested area in the CONUS was categorized into 841 climatic zones based on three forest types, annual total precipitation (30 mm intervals) and annual average temperature (2 °C intervals). An empirical orthogonal panel [48,49] was used to identify the pattern of these two climatic variables (horizontal axis—annual total precipitation and vertical axis—annual average temperature), and to associate forested grids to climatic zones (Figure S2;

Table 4. Definition of climatic zones for grouping pixels within a forested area. Three forest types with fixed ranges of annual total precipitation (30 mm intervals) and annual average temperature (2 °C intervals) yield a total of 5805 segments (3 × 129 × 15). Of these, only 841 climatic zones have at least one GLAS observation.

Forest Types (Deciduous, Evergreen, and Mixed Forests)	Annual Total Precipitation (mm)			Annual Average Temperature (°C)			Climatic Zones
	Lower Limits	Upper Limits	Intervals	Lower Limits	Upper Limits	Intervals	Effective
3	300	4,170	30	−5	25	2	841

).

The reasons for defining climatic zones are twofold: First, a direct comparison of GLAS heights, which represent actual tree heights, with potential tree heights predicted by the ASRL model is not valid. Second, optimization of the ASRL model for every forested pixel (over 1.3 million pixels) is not computationally practical. Thus, the optimization was performed at the climatic zone level for each of the three forest types.

4.2. Initial ASRL Model Prediction of Potential Tree Heights

ASRL model simulations were performed over forested areas (Section 3.1.2). The ASRL model predicts potential tree heights at 1 km spatial resolution using input climatic and ancillary variables (Section 3.1). There are notable disparities between model predicted tree heights and actual observations—the reason being that the ASRL model includes constant scaling exponents and parameters across different climatic regimes and forest types.

4.3. Optimization of the ASRL Model

Optimization of the ASRL model was designed to simultaneously adjust multiple scaling parameters. This optimization was aimed to minimize the difference between actual tree heights derived from GLAS data and tree heights predicted by the ASRL model (Figure 2). The underlying theoretical framework is based on Powell’s optimization methodology [50] that results in finding the minima of a multidimensional function. This algorithm is efficient for generating the convergence of a function due to (a) its bi-directional search algorithm over the vector of multi-variables and (b) nonessential calculation of derivatives for each variable [50]. A function involving three variables (merit function as in Equation (2)) was formulated and implemented based on Press et al.[51] and Kuusk and Nilson [52].

Three parameters of the ASRL model—area of single leaf (α), exponent for canopy radius (η) and root absorption efficiency (γ)—were selected for optimization. Initial values of these three parameters (α, η, and γ) were set to 13 cm², 1.14, and 0.33, respectively. These values are comparable to the representative values (averages) from the TRY database [53] and also based on Kempes et al.[14]. Although there are other physiological traits available from the TRY database [54], this study was limited to optimization of these three parameters.

The collection of solar radiation for plant growth is associated with the coefficient for canopy transmissions. Here, α produces the total leaf area based on the branching generation theory [55]. In the ASRL model, the canopy-level budget is collected from the energy budget in a single leaf based on the allometric geometry of canopy [14]. The value η controls the scaling of canopy radius with tree height, which is related to the rate of absorbed solar radiation. Lastly, γ determines the available flow rate given the incoming rate of precipitation within the root capture area. The tallest tree takes γ (=1/3) on average and local γ varies across different soil type and hydrology [14].

Kempes et al.[14] have performed the sensitivity tests of several allometric scaling exponents, and η showed the second least sensitivity. The value γ was tested in the optimization procedures of Kempes et al.[14], generating clear improvement of the model predictions. In this study, α was additionally selected for optimization due to (a) its strong relationship to net absorbed radiation, sensible and latent heat fluxes (e.g., [56]) and (b) considerable variability of α across different climatic zones and forest types [57,58].

The optimization process stops the iterative adjustment of the three parameters when it finds the maximum likelihood estimates of each parameter that result in minimizing the merit function. To reach an optimal solution, we implemented variable ranges (lower and upper boundaries as in the TRY database) for each of the input parameters such that 1 cm² ≤ α < 100 cm², 0.8 ≤ η < 1.5, and 0.1 ≤ γ < 0.8. Equation (2) shows the merit function,

$$J(p_1,p_2,p_3)={\sum }_{k=1}^3\left[{\left(p_k-p_{\mathit{kb}}\right)}^4{w}_k^2+\frac{{\left(p_k-p_{k0}\right)}^2}{{\mathrm{\Delta }}_{p_k}^2}\right]+{\sum }_{i=1}^N\frac{{\left(H_{\mathit{ASRL},i}-H_{\mathit{GLAS},i}\right)}^2}{{\mathrm{\Delta }}_H^2}$$

(2)

Here p_k is the selected ASRL model parameter (k; 1 = area of single leaf, 2 = exponent for canopy radius, and 3 = root absorption efficiency), p_k0 refers to the initial values of each parameter in the original ASRL model, p_kb refers to the boundary limits ((lower limit + upper limit)/2) for each parameter, Δ_pk is the standard deviation associated with each parameter with respect to initial values, w_k is a scalar weight (w_k = 0 when p_k ∈ [p_{k_lower-limit}, p_{k_upper-limit}], otherwise w_k = 10), N is the total number of comparison sets (i) for GLAS heights (H_GLAS) and ASRL model predictions (H_ASRL) with given parameter values for each climatic zone and Δ_H is the standard deviation associated with H_GLAS and H_ASRL. The iterative adjustment process continues until it finds the minimum of the function J (p₁, p₂, p₃) for each of the climatic zones.

A noteworthy limitation of this optimization exercise is that forest stand ages are not directly involved in the optimization process. Tree heights and growth rates vary depending on forest types and sites due to different growing conditions. Those are clearly related to forest stand ages [59–61]. When a tree ages, its height increases along with decline in the rate of its vertical growth over time—young forest stands (∼10 years) grow in the southeastern region, while old forest stands (∼900 years) inhabit the western coasts in CONUS [62]. However, it does not necessarily mean that our methodology neglects forest stand ages in tree height estimations. GLAS waveform data indirectly brings age information of forests into the ASRL model to find appropriate scaling parameters, as actual heights are associated with forest stand ages.

Performance of the ASRL model was tested by comparing GLAS tree heights and model predicted heights (with and without optimization). The goal was to show the efficacy of the optimization process. We calculated R² and root-mean-square-error (RMSE) from relationships between GLAS tree heights and model predicted heights in each climatic zone (Equation (3)).

$${\mathit{RMSE}}_{\mathit{ASRL}}=\sqrt{\frac{{\sum }_{i=1}^n{\left({\overline{H}}_{\mathit{ASRL},i}-{\overline{H}}_{\mathit{GLAS},i}\right)}^2}{n}}$$

(3)

Here H̄_ASRL is the mean of tree heights from ASRL model predictions (with and without optimization) in each climatic zone, H̄_GLAS is the mean of tree heights derived from GLAS waveform data in each climatic zone, and i refers to the number of climatic zones (n = 841).

The training datasets of GLAS used in model optimization are identical to the test datasets of GLAS used for evaluation [63]—this was first done to assess whether the optimization scheme was correctly implemented or not. It was not meant to establish validity of the optimized ASRL model. The actual evaluation of the optimized ASRL model was performed as detailed below.

4.4. Evaluation of the Optimized ASRL Model Results

The prediction of the optimized ASRL model was evaluated in two parts: (a) two-fold cross validation approach and (b) two inter-comparisons of optimized ASRL model prediction (H_{opt ASRL}) with forest canopy heights produced by Simard et al.[1] (H_Simard) and Lefsky [2] (H_Lefsky). Each evaluation performs inter-comparisons at the climatic zone level and at the pixel level.

4.4.1. Two-Fold Cross Validation

A two-fold cross validation approach was performed: That is, we randomly divided the original sample input data into two sets of training and test data. The first half of the GLAS tree heights was used as a training data to optimize the ASRL model in each climatic zone. The test data was generated by averaging the remaining half of the GLAS tree heights in each climatic zone and used for model evaluation purposes (Equation (4)). In addition, pixel level comparisons were performed to evaluate model prediction errors (H_{opt ASRLtraining} − H_GLAStest). We selected spatially corresponding tree height values (the nearest pixels) in pixel level comparisons.

$${\mathit{RMSE}}_{\mathit{ASRL}\hspace{0.17em}\text{two-foldcross}\hspace{0.17em}\text{valid}.}=\sqrt{\frac{{\sum }_{i=1}^n{\left({\overline{H}}_{\mathit{opt}\hspace{0.17em}\mathit{ASRL}\hspace{0.17em}\text{training}i}-{\overline{H}}_{\mathit{GLAS}\hspace{0.17em}\text{test},i}\right)}^2}{n}}$$

(4)

Here H̄_{opt ASRLtraining} is the predicted height by the optimized ASRL model using the first set of GLAS training data for each climatic zone, H̄_GLAStest is the mean of tree heights computed from the second set of GLAS test data in each climatic zone, and i refers to the number of climatic zones (n = 245). In these climatic zones, the number of pixels with GLAS tree height data was more than 20.

4.4.2. Inter-comparison with Other Forest Height Maps

The optimized model evaluations were additionally performed by comparing model predicted heights with H_Simard and H_Lefsky. Linear regression analysis between model predicted tree heights and the two maps was performed for each of the climatic zones. Pixel level evaluations used differences in histograms that were differentiated by forest types: deciduous, evergreen, and mixed forests. Some caveats are in order regarding these inter-comparisons: (a) the metric of forest height map in Lefsky [2] is Lorey’s height—basal area weighted mean height, while Simard et al.[1] and our research used maximum canopy height, (b) the forest height map of Simard et al.[1] does not allow tree height values >40 m, (c) both Simard et al. and Lefsky differ in their definition of forested lands, and (d) final products of Simard et al.[1] and Lefsky [2] are at different spatial resolution (1 km and 500 m, respectively) and different map projection.

To facilitate inter-comparison with H_Simard, we resampled and reprojected the forest height map of Simard et al.[1] to match our map of model predicted tree heights. The comparison was then performed over pixels that spatially corresponded to our definition of forested lands.

On the other hand, a direct comparison between Lefsky [2]’s forest heights and our results was not feasible for the reason of different measures (Lorey’s height versus maximum tree heights). Therefore, we used Lefsky [2]’s input GLAS heights data, rather than the final product of Lefsky [2]. For the inter-comparison, we averaged Lefsky’s GLAS heights falling in a pixel (=1 km) of the forest lands.

There are certain limitations of our analysis in the ASRL model predictions and evaluations: (a) up- and down-scaling approaches of resampling may cause potential errors due to the aggregation of heterogeneity in finer grids and the neglect of discontinuity in coarser datasets [64,65] and (b) the reprojection possibly results in certain modification of true pixel values [66].

5. Results and Discussion

5.1. Initial ASRL Model Predictions of Potential Tree Heights

A continuous map of potential tree heights (H_{potential ASRL}) was generated with the unoptimized ASRL model at 1 km resolution (Figure 3(a)). Maximum potential tree heights were greater than 50 m in both the Northeastern Appalachian and Pacific Northwestern forest corridors. The model predicted lower values of potential tree heights (≤35 m) in the Southeast.

We noted discrepancies between model predictions and GLAS tree heights (H_GLAS; actual tree height; Figure S1). A low correlation was observed (Figure 3(b)) in each climatic zone (R² = 0.06; RMSE = 22.8 m). In addition, there was significant skewness in the histograms of actual (mean = 31.3 m; standard deviation = 11.5) and potential (mean = 45.5 m; std. = 23.6) tree heights at the pixel level (Figure 3(c)). Tree heights were overestimated especially in the northeastern forests as compared to H_GLAS (Figure S3(a)). A plausible reason could be that the ASRL model does not accurately reflect the spatial/temporal dynamics in the estimation of internal flow balances (metabolic flow requirement, available flow, and evaporative flow) across different eco-climatic regimes and forest types [14].

5.2. Optimized ASRL Model Predictions

The optimized model was then used to generate a spatially continuous map of tree heights (H_{opt ASRL}; Figure 4(a)). We noted a significant improvement in predictions of tree heights both at the climatic zone level (Figure 4(b)) and individual pixel level (Figure 4(c); Figure S3(b)): (a) the RMSE decreased from 22.8 m (without optimization) to 3.1 m (after optimization) with an increase in R² from 0.06 to 0.8 (P < 0.01); (b) the histograms show a better agreement between distributions of GLAS tree heights (mean = 31.3 m; std. = 11.5) and the optimized model predictions (mean = 30.4 m; std. = 8.5); and (c) relatively smaller model prediction errors over the Northeastern Appalachian and Pacific Northwestern forest corridors as compared to the unoptimized ASRL model predictions.

Figure S4 shows that the ASRL model prediction errors at the individual pixel level decreased from 15.10 m (H_GLAS − H_{potential ASRL}) to −0.80 m (H_GLAS − H_{opt ASRL}). However, the optimized ASRL model poorly predicted tree heights over complex terrains (e.g., ∼20 m underestimation for the redwood stands in the Pacific Northwestern mountains of California and Oregon; Figure S3b). Other GLAS-based models also reported relatively large prediction errors in the estimation of tree heights [1] and biomass [3] in these forests. Interpolation of annual precipitation (e.g., [67]) and temperature (e.g., [68]) may have produced large uncertainties in climatic variables that are sensitive to topographic features. Note that these are critical inputs to the ASRL model. Other plausible reasons for this discrepancy may be: (a) GLAS undersampling for some of the climatic zones (Figure S5) that resulted in fewer comparison sets in the merit function (Equation (2)) and (b) topographic influence on GLAS waveform data which could not perfectly be rectified by our slope gradient filter (Table 3).

The optimized parameters are shown in Figure 5. There are notable changes in the optimal area of single leaf (initial value α: 13.0 cm²) that ranged from 1.5 cm² to 90.0 cm². The root absorption efficiency (initial value γ: 0.33) converged to a relatively narrower range of values (from 0.05 to 0.65), while ∼80% of the optimized exponent for canopy radius (η) fell within the range of ±10% of its initial value (1.14). Kempes et al.[14] have also reported a stable median relative error against the percent change of a single scaling parameter (i.e., η).

The area of single leaf of deciduous forests (mean α = 19.3 cm²) was higher than that of evergreen forests (mean α = 9.1 cm²). The original ASRL model precludes inclusion of forest types. The optimization process allows combining allometric scaling laws with features that are representative of specific forest types. Optimized α values are well correlated with the variability in forest types, annual total precipitation and annual average temperature in each climatic zone. Warm (annual average temperature = ∼15 °C) and wet (annual total precipitation ≥ ∼1,500 mm) regions displayed a larger value of α for both deciduous and evergreen forests. In cold regions (annual average temperature = ∼5 °C), the optimized value of α for evergreen forests increased with annual total precipitation. These results are supported by other studies that examined relationships between leaf traits and environmental conditions [69–71].

Similar trends in the optimized γ values were observed in warm and wet regions. However, evergreen forests generally showed higher optimized γ values compared to deciduous forests in relatively dry regions. Water availability is spatially heterogeneous for an individual species within a location [72]. For example, evergreen and deciduous plants in dry regions have different root systems and water use efficiencies (evergreen > deciduous as in [73,74]). Kempes et al.[14] have demonstrated an improvement of the ASRL model based on optimization of γ that generated a lower variance in the model error.

5.3. Evaluation of the Optimized ASRL Model

5.3.1. Two-Fold Cross Validation Approach

Figure 6(a) shows the two-fold cross validation comparison (R² = 0.59; RMSE = 3.31 m; P < 0.01). Histograms comparing the test GLAS heights (mean = 30.8 m; std. = 10.7) and tree heights predicted by the optimized model (mean = 30.6 m; std. = 8.4) show considerable similarity (Figure 6(b)), even though it gives relatively less correlations than using all of valid GLAS tree heights. The satisfactorily low prediction errors (mean = −0.61 m; std. = 12.91) are shown in Figure 6(c). We achieved the stability of the optimized model predictions from the two-fold cross validation.

5.3.2. Inter-comparison with Other Forest Height Maps

Forest height maps from Simard et al.[1] and our study portray similar patterns of tree heights over the CONUS: (a) taller trees (> 40 m) in the Pacific Northwestern forests of California and Oregon, (b) relatively medium-to-tall trees (30 to 40 m) in the northeastern forested regions, and (c) smaller trees (∼20 to 30 m) along the Great Lake and the Mississippi River basin. It should be noted that the regression tree procedure described in Simard et al.[1] is based on the GLAS altimetry variables of the Gaussian decomposition (signal beginning, signal end, and last Gaussian peak). Simard et al.[1] also included a similar set of environmental layers related to elevation, temperature, and precipitation.

Figure 7(a) depicts a scenario where H_{opt ASRL} is relatively higher in the northwestern and northeastern forested regions as compared to H_Simard. At the scale of climatic zones (Figure 7(b)), the optimized ASRL model predictions are moderately correlated to height values derived by Simard et al.[1] (R² = 0.45; RMSE = 8.01 m; P < 0.01). Average values of H_{opt ASRL} for each of the climatic zones were usually higher. Figure 7c shows that the differences between these two maps are nearly independent of forest type. The differences in height values can likely be attributed to differences in definitions of forests—the map from Simard et al.[1] used forested areas corresponding to classes such as mosaic crops, open forest, and saline flooded forests. An added caveat, as noted in Simard et al.[1], was the inability of their regression tree model to simulate forest heights >40 m.

We also compared our forest height map with Lefsky’s [2] original GLAS-based tree heights. Figure 8(a) shows a one-to-one comparison between the average height values obtained from H_{opt ASRL} and H_Lefsky for each of the climate zones. Overall, there is a moderate correlation (R² = 0.41; RMSE = 6.72 m; P < 0.01) between H_{opt ASRL} and H_Lefsky. Mean values of H_{opt ASRL} are an underestimate. Figure 8b shows the pixel level difference (H_{opt ASRL} – H_Lefsky) histograms for three forest types—deciduous forests show higher differences (mean = −5.7 m; std. = 10.3) followed by the evergreen forests (mean = −4.2 m; std. = 12.6) and mixed forest types (mean = −4.5 m; std. = 9.5). A plausible reason could be that Lefsky [2] applied a different height retrieval procedure (statistical) based on full GLAS waveform extents, while our study used the standard Gaussian decomposition approach.

6. Concluding Remarks

An optimization of the Allometric Scaling and Resource Limitations (ASRL) model with Geoscience Laser Altimeter System (GLAS) waveform data was performed to generate a spatially continuous map of tree heights over the continental USA (CONUS) at 1 km resolution. The optimization is designed to minimize differences between actual heights (based on GLAS waveforms) and potential tree heights predicted by the ASRL model. This study covered all forested lands with over 50% tree cover. These were categorized into 841 climatic zones based on forest types (deciduous, evergreen, and mixed forests), fixed intervals of annual total precipitation (30 mm) and annual average temperature (2 °C). The optimization procedure simultaneously adjusted three model parameters (area of single leaf, α; exponent for canopy radius, η; and root absorption efficiency, γ) in each of the climatic zones.

After testing for correctly implementing the optimization technique, tree heights predicted by the optimized model were first evaluated using a two-fold cross validation approach. Regression analysis was used to assess the correlation between predictions of tree heights by the optimized model (H_{opt ASRL} training) and test GLAS tree heights (H_GLAStest) in all climatic zones. Mean values of H_{opt ASRL} explained 59% of the variability in H_GLAStest mean estimates in each of the climatic zones and, on average, showed an estimation error of 3.31 units of height. A similar evaluation of the optimized ASRL model was performed at FLUXNET sites—this is detailed in the second of this multi-article series [25]. A comparison at the pixel level to quantify the skewness between H_{opt ASRL} training (mean = 30.8 m; standard deviation = 10.7) and H_GLAS (mean = 30.6; std. = 8.4) was performed. Predicted tree heights by the optimized model agreed better with GLAS tree heights (mean = −0.6 m; std. = 12.9) in comparison to the estimates from the unoptimized ASRL model. However, the optimized ASRL model still poorly predicted tree heights over the Pacific Northwestern Mountains of California and Oregon.

Second, tree height predictions by the optimized ASRL model were compared with available forest height products derived independently but from the GLAS data—Simard et al.[1] (H_Simard) and Lefsky [2] (H_Lefsky). The results indicate moderate correlation between optimized ASRL model predicted heights and forest heights from Simard et al.[1] and Lefsky [2] for all climatic zones (R² = 0.45 and RMSE = 8.01 m for H_Simard; R² = 0.41 and RMSE = 6.72 m for H_Lefsky). H_{opt ASRL} was an overestimate compared to H_Simard and an underestimate compared to H_Lefsky and with significant skewness at the individual pixel level—these discrepancies can be attributed due to different definitions of heights and forested lands between these studies and certain inherent limitations of the various approaches.

Predictions of tree heights by the ASRL model were clearly improved by the optimization technique reported in this article. The optimization successfully compensated for certain limitations of the original ASRL model, which did not account for effects related to spatio-temporal variability in climatic-regimes and forest types. The results demonstrate the potential for a more generic applicability of the ASRL model for estimation of tree heights. Nevertheless, the optimized ASRL model still yields ambiguous results over complex terrains, possibly due to uncertainties in input climatic data and topographic effects in the GLAS waveform data. The optimization methodology reported in this article has certain limitations: e.g., (a) a limited number of scaling parameters (α, η, and γ) were explored in the model optimization, (b) stand age was not directly considered in the optimization, (c) soil conditions were neglected in the optimization and (d) we assumed that allometric scaling laws at individual tree level were applicable at larger scales. Also, our analysis could not take into account the uncertainties derived from resampling and reprojection of maps and data at different scales and projections. Alleviation of these limitations should be addressed in future articles in this series.

Supplementary Information

remotesensing-05-00284-s001.docx