The study area comprises 31,209 ha inside the Willamette National Forest, OR, United States covered by different remote sensing datasets that include a recent lidar acquisition and a 30 m resolution map with above-ground biomass (AGB) predictions (Figure 1). Details on these datasets are provided in sections “Case 1: Fay–Herriot Models for New Inventories” and “Case 2: Fay–Herriot Models to Update Inventories.” Elevations range from 450 to 1700 m above sea level. Two forks of the Santiam river cross the study area from East to West and have numerous tributaries that form a complex drainage network where slopes do not have a dominant orientation. Conifers dominate vegetation with Douglas-fir, Pseudotsuga menziesii (Mirb) Franco, the most abundant species, and other conifers such as noble fir, Abies procera Rehder, silver fir, Abies amabilis Douglas ex J.Forbes, western hemlock, Tsuga heterophylla (Raf.) Sarg., and western red cedar, Thuja plicata Donn ex D.Don, as secondary species with a much lower abundance. Hardwood species have a minor presence with red alder, Alnus rubra Bong, and golden chinquapin, Chrysolepis chrysophylla (Douglas ex Hook.) Hjelmq., as the most important species in this group.

The study area contains 1616 stands with different management goals. Stand boundaries are the result of a continuous effort performed by forest managers in the study area and is based on the management history, structure, and composition of the forest. Stands are classified according to their management objectives into “Natural” and “Managed” (Figure 2). While the terms natural and managed can be the subject of lengthy discussions, we will keep this terminology as it is used in the FSVeg database. There are 696 natural stands and 920 managed stands. Natural stands occupy approximately two-thirds of the area. They are typically of larger size (i.e., average size = 31.04 ha, median size = 11.61 ha) than managed stands (i.e., average size = 9.90 ha, median size = 7.94 ha) (Figure 2). Managed stands have a past history of silviculture entry and often include artificial regeneration. In most cases, even-aged structures are subject to thinning and logging operations. Natural stands are subject to less intense management and tend to have larger dimensions and a larger internal variability in forest structure and ages (Figure 2).

For both Case 1 and Case 2, we considered estimating, for every stand in the study area, the total of AGB, and merchantable volume (V), for the year 2018, both expressed on a per unit area basis. Thus, for every stand, the unknown parameter of interest was

when considering AGB and

when considering V. In equations 1, 2 AGB_ti and V_ti are the AGB and V of the t-th tree in the i-th stand, and N_i and A_i are, respectively, number of trees and the area of the i-th stand. It is important to note that while μ_i, AGB_ti, V_ti, and V_ti were all unknown quantities, the stand area was known for every stand in the study area.

We used the Forest Vegetation Simulator (FVS), to compute an estimate, μ^g⁢i⁢j, of each parameter of interest for each VRP using the Horvitz–Thompson (HT), estimator

In equation 3, y_ijt represents either AGB_tij or V_tij for the t-th tree measured in j-th VRP in the i-th sampled stand and EF_ijt represents their respective expansion factors. The number of measured trees in the j-th VRP in the i-th stand is n_ij and the subindex g in μ^g⁢i⁢j indicates that it is a direct estimate based on the ground data.

For each sampled stand, VRP estimates μ^g⁢i⁢j from the n_i plots measured in the stand were averaged to produce a final direct ground estimate μ^g⁢i⁣* of AGB and V

A summary of the stand estimates based on the VRP data is presented in Table 1.

The HT estimator is unbiased and each VRP is assumed to provide an independent sample drawn under a sampling design that remains constant for all VRP in the stand. Thus, for a given stand, all μ^g⁢i⁢j were considered to be realizations of a random variable with mean μ_i and unknown variance σe⁢i⁢02. The final direct estimate for the stand, μ^g⁢i⁣*, is the average of n_i independent and identically distributed random variables. Therefore, μ^g⁢i⁣* is also a random variable with mean μ_i and its variance, σe⁢i2 equals σe⁢i⁢02ni. This allows establishing the following relation, equation 5, between the stand estimate μ^g⁢i⁣*, the unknown parameter of interest μ_i and the sampling error e_i

For any two stands, sampling errors are assumed to be independent of each other. Furthermore, due to the unbiasedness property of the HT estimator, errors are assumed to be distributed with zero mean and variance σe⁢i2. The variance σe⁢i⁢02 is unknown, but an unbiased estimator can be obtained pooling together the estimates of all VRP in a given stand as

Based on equation 6 the variance and the mean square error of μ^g⁢i⁣* is estimated using

Estimators in equations 4, 7 are typically used in stand-level inventories using only ground data when reporting estimates and measures of uncertainty for sampled stands.

Different patterns were observed regarding the selected models for Case 1 and Case 2 and for natural and managed stands. Regardless of the case and variable of interest, residuals for both cases tended to be centered around zero, and no significant departures with respect to the model assumptions were observed. The variability of these residuals was substantially larger for natural stands than for managed stands. Models for Case 2 tended to provide a shorter range of predicted values when compared with the models for Case 1 (Figure 4). In general, models for natural stands had a smaller number of predictors. This result was expected. The number of sampled stands is substantially smaller for natural stands than for managed stands, therefore estimated model coefficients for natural stands tend to have larger standard errors and less coefficients appeared as significant in the fitted models. The intercepts for models for natural stands for Case 1 were not significant and were removed. For both managed and natural stands, models for Case 1 had lower values of σ^v2. For managed stands, σ^v2 for AGB and V of the FH models for Case 2 were 26 and 42% larger than σ^v2 for Case 1, respectively. For natural stands, we observed a 4.14 and 3.27-fold increase in σ^v2 for AGB and V when comparing the values obtained for Case 2 with those obtained for Case 1 (Tables 2, 3). This indicates FH models based on a recent lidar acquisition explain more variance than the models for Case 2 (Tables 2, 3). When comparing models for managed and natural stands obtained under a given case scenario, we observed that for Case 1, σ^v2 for AGB and V in natural stands was 2.35 and 3.48 times larger than in managed stands, respectively. For Case 2, models for natural stands explained almost no variance, and for AGB and V, σ^v2 was 7.73 and 8.02 times larger than the one obtained in managed stands.

Selected models for Case 1 and Case 2 were used to obtain stand-level estimates and their associated mean squared errors for sampled and unsampled stands (Figure 5). For both cases and response variables, estimates based on FH models for managed stands tended to be smaller than estimates for natural stands and the same pattern was observed for the corresponding rmse (Figure 5). When considering rrmse, managed stands showed larger relative uncertainties. This is partly caused by the fact that managed stands stock substantially less AGB and V.

For both AGB and V, estimates and rrmse obtained for Case 1 tended to agree spatially with estimates and rrmse for Case 2 (Figure 5). For natural stands, Spearman rank correlation between estimates for Case 1 and Case 2 was 0.22 (p-value = 1.04 × 10^–8) for AGB and 0.12 (p-value = 1.39 × 10^–3) for V. The low agreement for natural stands seems to be caused by the low explanatory power of the models for Case 2 in natural stands. For managed stands, Spearman rank correlation between estimates for Case 1 and Case 2 were 0.69 (p-value < 10^–6) for AGB and 0.50 (p-value < 10^–6) for V, and when each map was grouped into 10 deciles, these categories tended to coincide. The same occurred with the estimated rmse maps (see Supplementary Figure 1).

For both types of stands, FH models for Case 1 and Case 2 provided improvements for direct ground estimates for AGB and V. Estimates from Case 1 were consistently more precise than those from Case 2 (Figure 6). This result was expected after observing the values obtained for σ^v2 for the different models. In general, improvements in efficiency and differences between cases were larger for natural stands (Figure 6). For natural stands, improvements in efficiency for Case 1 had an average △_effi of 37.36% for AGB and 33.10% for V (Table 4). For Case 2, the average of △_effi was 20.19% for AGB and 19.25% for V. For managed stands, the average of △_effi for Case 1 was 20.29% for AGB and 19.91% for V, and for Case 2, the average of △_effi was 17.55% for AGB and 16.05% for V (Table 4). The smaller values of △_effi in managed stands is explained by the larger homogeneity of this type of stands for which many of the direct ground estimates were already precise, leaving little room for improvements to the FH models. Differences in △_effi between cases for managed stands were smaller than the differences for natural stands. This seems to be the consequence of both the low explanatory power of the auxiliary variables for Case 2 in natural stands and the smaller room for improvements in managed stands.

When both scenarios were compared, estimates for Case 1 had, in general, lower errors than estimates from Case 2. The differences between cases were more important for natural stands than for managed stands. Multiple studies have shown that forest structural attributes correlate better with lidar auxiliary information than auxiliary variables from optical sensors. Auxiliary information for Case 1 proceeded from lidar. For Case 2, we used a previous remote sensing-based inventory the CMS1-AGB map, which heavily relies on metrics derived from 30-year climate normals, topographic, and Landsat variables, to which we added proxies for disturbances directly derived from Landsat images. This explains that Case 1 outperforms Case 2 for all response variables and stand types. Nevertheless, estimates from Case 2 are more efficient than direct ground estimates and regardless of the case, the rank correlations between estimates for Case 1 and Case 2 for managed stands indicated that both methods agree in the way the sort stands according to the predicted AGB or V. These results indicate that Case 2 is also useful for managed stands, and that certain management decisions, for example, concentrating harvest activities in the 10% of the managed stands with more volume, would tend to coincide regardless of which map (i.e., Case 1 or Case 2) is used to inform those decisions (see Supplementary Figure 1).

Two remarks should be made about the initial map for Case 2. On the one hand, the CMS1-AGB map has spatial and temporal coverage that cannot be matched by previous maps based only on lidar. Thus, this map can be used to develop similar stand-level inventories anywhere in the western United States. Furthermore, the multitemporal component of this map allow for possible applications of FH models to estimate changes and monitor vegetation dynamics that are not an option using single date lidar data. On the other hand, the CMS1-AGB map is expected to provide predictions with more noise than similar maps based on the ABA method and lidar data. This implies that results obtained in this study for Case 2 might improve substantially when the previous inventory is an ABA lidar-based inventory. Many countries have developed nationwide lidar acquisitions or are on the verge of completing such data collections, and national forest inventories can provide the necessary fixed area plots to use the ABA to develop maps based on lidar at national or regional scales. The effort required to develop these maps is large, and re-mapping is not expected to happen with a high frequency. This indicates that a potential niche of application of FH models and Case 2 is updating national or regional level ABA maps.

When comparing natural and managed stands, we observed that the former had larger estimated model variances, resulting in stand-level estimates with larger uncertainties in absolute terms. These differences are explained by the fact that natural stands are inherently more complex and variable than managed stands. Part of that complexity is not captured by predictors computed at the stand level. Relative uncertainties (i.e., rmse) were lower for natural stands than for managed stands. The higher stocking levels cause that in natural stands. Improvements in efficiency for natural stands were larger than those observed for managed stands, especially for Case 1. Finally, the differences in △_effi for Case 1 and Case 2 were relatively small for managed stands (i.e., about 3% difference between average values of △_effi) but large for natural stands (i.e., approximately 15% difference between average values of △_effi). The interaction of three different factors can explain this differentiated behavior. The first is that managed stands are relatively homogeneous units, and their direct ground estimates were more reliable than those obtained in natural stands. Thus, the potential for improved managed stands was more limited and made differences between cases smaller. Another factor is that the larger stocking levels are frequently associated with remote sensing predictions with larger uncertainties (Magnussen et al., 2014; Mauro et al., 2016; Breidenbach et al., 2018). Large uncertainties in the previous inventory must result in a poorer characterization of the initial state of the stands, which partially explains the performance drop for Case 2 in natural stands. The auxiliary information for Case 2 is primarily based on metrics derived from 30-year climate normals to capture steep AGB (and V) gradients, which largely compensated for signal saturation of Landsat variables (Hudak et al., 2020), albeit without sensitivity to local variation in stand structure. On the other hand, the lidar data used for Case 1 neither saturates in forested areas with closed canopies nor is insensitive to structure variation between or within stands, thus elevating the performance of the FH models for Case 1 compared to Case 2.

For all cases, response variables and stand types analyzed in this study, FH models allowed for gains in efficiency with respect to direct ground estimates. These results are concurrent with previous research using FH models in forest inventories (Goerndt et al., 2011; Magnussen et al., 2017; Mauro et al., 2017; Breidenbach et al., 2018; Ver Planck et al., 2018) and confirm that FH models: (1) allow using ground measurements that are easier to obtain than those used in ABA approaches and (2) results in efficiency improvements when compared to methods based only on ground data. Thus, while research efforts on using FH models are still necessary, there is substantial evidence that these models can play an essential role in operational forest inventory applications.

To the best of our knowledge, FH models have not been used in any operational forest inventory, and it somehow surprises how little attention FH models have received in the literature. While some research applications of FH models exist, the study of this type of model has been negligible compared to applications using conventional ABA approaches. The dominance of the traditional ABA approach can be explained by (1) its ability to produce high resolution maps with predictions of forest attributes and (2) its typically better predictive performance than FH models (Mauro et al., 2017; Breidenbach et al., 2018; Green et al., 2019). However, developing ABA models is not always a possibility. There are many scenarios where FH models can be a very appealing alternative; for example, only stand-level inventory data may be available. During the last decades, forest inventories have been consistently less constrained by the availability of useful auxiliary information, but the costs associated with ground data collection have increased, or at least not decreased at comparable rates. Simply put, ground data are too valuable to ignore, and FH models allow for an effective combination of those valuable datasets with different sources of auxiliary information.

A critical difference between traditional ABA models and FH models is the auxiliary information used by each technique. Workflows for preprocessing auxiliary information for traditional ABA models are well established and documented, with multiple tools available to implement these processing steps (i.e., Mc Gaughey, 2019; Roussel et al., 2020); this is not the case for FH models. In this study, we used as predictors stand summaries of: (1) gridded products (i.e., gridmetrics rasters) generated with FUSION (Mc Gaughey, 2019), (2) previously mapped estimates of forest attributes, and (3) changes in Landsat imagery from LCMS or computed using Google earth engine (Gorelick et al., 2017). Summarizing the entire point clouds within the stands under analysis is an alternative used in previous studies to compute lidar-based predictors (Ver Planck et al., 2018). Both options are valid from a methodological perspective as they provide standardized ways to compute auxiliary variables. Their effectiveness can differ if one preprocessing technique provided auxiliary variables that correlated better with the target responses than the other. However, as far as we know, no study to date has analyzed the differences in performance and tradeoffs of these two methods to generate stand-level predictors for FH models. Thus, this is an area where future research can help in establishing standardized processing workflows for lidar-assisted forest inventories using FH models.

This study presents two case scenarios in which basic FH models are used with VRP and demonstrates that FH models are a suitable alternative to use available auxiliary information to improve the efficiency of the estimation process. Our analysis presents a baseline for stand-level FH models and could be improved in different ways. One way is developing models that account for spatial correlations like those developed by Ver Planck et al. (2018). Another option is to use FH variants where the model variance is not constant (Breidenbach et al., 2018). A third option is to use multivariate FH models where correlations between different response variables can be considered to improve the results of univariate models (Benavent and Morales, 2016; Frank, 2020). In all cases, one factor that must be constantly considered is that estimates obtained from FH models are always based on a model (i.e., “model-based”). Thus, extrapolations entail high risks of producing biased estimates and model validation steps are critical to ensure that the fitted models correctly describe the populations under study.

Based on our findings and previous results (Goerndt et al., 2011; Breidenbach et al., 2018; Ver Planck et al., 2018), we envision that niche of application of stand-level FH models is not a replacement of traditional ABA methods but a complement for situations in which the time and resources available for ground data collection are limited or fixed radius plots with precise locations to develop ABA models are otherwise unavailable. This niche is larger than it might seem a priori for several reasons. One reason is that obtaining accurate coordinates for the ground measurements is not a constraint for FH models. For example, only the identifier of the stand where each ground observation was taken was necessary to develop this study using FSVeg data. Another reason, and probably the most compelling one, is that FH models can be applied with data from VRP or other sampling techniques such as sector plots (Iles and Smith, 2006) or transects (Warren and Olsen, 1964; Woodall and Monleon, 2008). This flexibility indicates numerous applications for fast inventories and monitoring problems in which FH models can be the preferred alternative. These applications include, but are not limited to, annual inventories for timber sales or fast updates of inventories after events like floods or wildland fires.

Improved AGB and V benefit both private companies and public land management agencies. Given the extremely high cost of establishing ground plots and the increasing demand for accurate biomass and carbon stock assessment, the inventory solution will require the innovative use of combined sources of remotely sensed and other auxiliary data. FH models based on VRP allow using remotely sensed information combined with an operative ground truth data collection and enable cost-effectively estimating forest attributes. In this study, the FH models have shown to be a viable and flexible option to estimate AGB or V and maximize the utility of both the ground inventory and environmental datasets. Moreover, different information for forest management planning is required at different levels or scales. For tactical planning, reasonably precise and unbiased estimates of forest variables for individual stands or polygons are already obtained using VRP because of its low cost and sampling efficiency at the stand level. Thus, FH models are an alternative for many established inventory programs to integrate their VRP data with lidar or other remote-sensing datasets to obtain more efficient and better information for sustainable forest management.