Forecasting Forest Inventory Using Imputed Tree Lists for LiDAR Grid Cells and a Tree-List Growth Model

A method to forecast forest inventory variables derived from light detection and ranging (LiDAR) would increase the usefulness of such data in future forest management. We evaluated the accuracy of forecasted inventory from imputed tree lists for LiDAR grid cells (20 × 20 m) in spruce (Picea sp.) plantations and tree growth predicted using a locally calibrated tree-list growth model. Tree lists were imputed by matching measurements from a library of sample plots with grid cells based on planted species and the smallest sum of squared difference between six inventory variables. Total and merchantable basal area, total and merchantable volume, Lorey’s height, and quadratic mean diameter increments predicted using imputed tree lists were highly correlated (0.75–0.86) with those from measured tree lists in 98 validation plots. Percent root mean squared error ranged from 12.8–49.0% but was much lower (4.9–13.5%) for plots with ≤10% LiDAR-derived error for all plot-matched variables. When compared with volumes from 15 blocks harvested 3–5 years after LiDAR acquisition, average forecasted volume differed by only 1.5%. To demonstrate the novel application of this method for operational management decisions, annual commercial thinning was planned at grid-cell resolution from 2018–2020 using forecasted inventory variables and commercial thinning eligibility rules.


Introduction
The use of light detection and ranging (LiDAR) to produce forest inventories has significantly changed forest management [1,2].Around the time of LiDAR acquisition, tree-level data (e.g., species, diameter at breast height (DBH), total height) are collected from calibration ground plots and used to calculate plot-level variables of interest (e.g., average tree height, total volume).These plots are selected to represent the range of structural variability within the area of interest [3].LiDAR point cloud statistics (e.g., mean return height, standard deviation of return heights, height percentiles) are calculated for the corresponding plot location and area.Using parametric or non-parametric methods, the statistical relationship between the variables of interest and point-cloud statistics is determined and used to predict variable values in non-sampled areas at grid-cell resolutions matching the size of the calibration plots (e.g., 20 × 20 m) [3,4].Such model-based estimates derived from LiDAR have been shown to provide accurate forest-and stand-level estimates of average tree DBH [5,6], average tree height [7,8], crown base height [9,10], basal area [11,12], and volume [13,14], with good spatial characterization of within-stand variability.Traditional forest inventory methods (i.e., sampling stands from broad species composition groups at various stages of development and management types to construct average age-based yield curves that are then applied to photo-interpreted forest inventory polygons), while accurate at the forest-level [15,16], do not match the stand-level accuracies of LiDAR estimates without expensive, supplemental field sampling and do not offer any of the within-stand spatial variability [1,17] that is needed for precision forestry (i.e., the prescription of silvicultural and harvest treatments at high spatial resolution).
Commercial thinning is one silvicultural treatment that requires spatially accurate estimates of forest inventory.Commercial thinning is a partial harvest treatment, typically prescribed in immature managed stands, with the goal of removing merchantable stems to provide economic return and improve growing conditions for residual stems [18,19].While thinning does not typically increase annual stand volume accrual [20,21], it does increase light, soil moisture, and nutrient availability for the residual crop trees [22,23], increasing their size and value [21,24].Accurate timing of commercial thinning entries, with respect to stand development, is vital to ensure effective response from the treatment and to realize maximum economic benefits.Optimally-timed commercial thinning that is conducted once density-dependent mortality is imminent will provide the greatest number of harvested stems with the highest average tree volume [25].Commercial thinning too early will result in lower average tree volume and reduced economic return [18], while thinning too late will result in lower volume due to tree mortality and lower live crown ratio [21].A live crown ratio of 35-40% has been cited for a number of conifers as the minimum level before tree vigor and response to thinning is reduced [18,25].LiDAR can provide the inventory data required to delineate and prioritize thinning treatments [26], or to identify LiDAR grid cells (e.g., 20 × 20 m areas) that are eligible for thinning [27].
While LiDAR-derived forest inventory can provide inventory data to make sound commercial thinning decisions, any forest inventory becomes less accurate and useful as the years since data acquisition increase.To increase the longevity of LiDAR-derived inventory and reduce the costs associated with LiDAR acquisition, processing, and inventory predictions, simplified process-based and size-class growth models have been used to forecast stand-level variables at grid-cell resolution [28][29][30].In addition, tree-level variables have been forecasted using tree lists imputed for LiDAR grid cells and a tree-list growth model [31].Tree-list growth models (e.g., [32,33]) predict growth and yield of individual trees on an annual or periodic basis using species-specific diameter, height, mortality, and live crown ratio increment equations to derive forecasted inventory [34,35].Providing forest managers with forecasted tree-level inventory for each LiDAR grid cell would allow users to track changes in species-specific products and identify grid cells that meet harvest requirements.
A major constraint of using tree-list growth models to forecast LiDAR inventory is the tree-level data (i.e., species, DBH, and height) required as model inputs.Predicting tree-level inventory for LiDAR grid cells requires a substantial number of calibration plots covering the full range of inventory conditions [3,36].In studies by Peuhkurinen et al. [37] and Packalén and Maltamo [38], close to 500 calibration plots were needed to produce species-specific diameter distribution estimates for their respective 1200-and 2000-ha study areas.In similar stand types, species-independent diameter distributions can be produced with as few as 50-100 plots on much larger land bases [39,40].If tree-level inventory is produced for a land base, incorporating it into a tree-list growth model introduces additional challenges.As each LiDAR grid cell would inevitably contain its own unique combination of estimated tree-level inventory, each cell would need to be forecasted individually.This becomes computationally challenging for large land bases.
A common method used for imputing forest inventory data is k-nearest neighbor imputation.k-nearest neighbor is a non-parametric method for selecting the most similar k neighbors for target data from a library of reference data, typically ignoring any spatial distance between the target and reference data.Neighbors are determined based on absolute differences between indicator variables that are shared between both the target and reference data.This method was originally used with forest inventory to impute complex inventory data, such as basal area distribution, for forest stands [41,42] and more recently has proved to be a powerful tool to impute LiDAR-derived stand-level variables [4,43].When the k-nearest neighbor method is used with LiDAR data to impute inventory variable values, LiDAR point-cloud statistics are used as the indicator variables to determine the nearest neighbors.This method requires that all reference data have available LiDAR point cloud statistics, ultimately limiting the amount of reference data available for most studies.Recently, Lamb et al. [44] used non-parametric k-nearest neighbor methods to impute tree-level inventory for LiDAR cells in spruce (Picea sp.) plantations in New Brunswick, Canada.Six LiDAR-derived variables (i.e., total and merchantable basal area, total and merchantable volume, top height, and merchantable quadratic mean diameter) and known planted tree species were used to match stand structural variables estimated from LiDAR to those in a library of sample plots.Since matches were based on shared inventory variable values and not LiDAR point-cloud statistics, sample plots in the reference library did not need to be located within the LiDAR acquisition area.As a result, over 5500 sample plot measurements were used to impute tree lists for 20 × 20 m LiDAR grid cells across 83,000 ha of spruce plantations.Matches were determined based on planted species and minimum sum of squared differences between six inventory variables.With this method, only 100 LiDAR calibration plots were required to predict the LiDAR-derived inventory variables for the 83,000 ha of spruce plantations.LiDAR grid cells were matched with only one plantation plot measurement (k = 1), which limited the number of unique tree lists to the 2870 spruce plantation plot measurements that were matched.This provides a framework for forecasting LiDAR grid cells at the tree-level without the need for a large number of calibration plots, while minimizing computational challenges associated with forecasting tree lists for cells.
The objective of this study was to assess the accuracy of forecasted inventory using tree-level inventory imputed for LIDAR grid cells in New Brunswick spruce plantations based on the methods of Lamb et al. [44].To do so, imputed tree lists were forecasted using a regional tree-list growth model [33] calibrated using growth and mortality data from permanent sample plots in spruce plantations from the study area.We hypothesized that inventory increments predicted using observed tree lists would be strongly correlated to increments predicted from imputed tree lists and the locally calibrated growth model.This would provide forest managers with a method to extend the longevity of LiDAR-derived inventory and allow use of the forecasted data to make operational, tactical, and strategic management decisions.

Study Area
The Black Brook District is a 210,000 ha forest in northwestern New Brunswick, Canada (47 • 9'51" N to 67 • 55'27" W) (Figure 1a), owned by J.D. Irving, Limited.The District is within Ecoregion 3 [45] and contains some of the most intensively managed lands in Canada [46].The forest is comprised primarily of uneven-aged, shade-tolerant hardwoods, managed using selection and group shelterwood harvests (65,000 ha) and spruce plantations, commercially thinned one to three times before final harvest (83,000 ha).Spruce plantations include 37,000 ha of black spruce (Picea mariana (Mill.)B.S.P.), 33,500 ha of white spruce (Picea glauca (Moench) Voss), 11,000 ha of Norway spruce (Picea abies (L.) Karst.), and 1500 ha of red spruce (Picea rubens Sarg.).Close to 400 permanent sample plots, each measured every 3-5 years since 1980, have been installed in these spruce plantations to monitor growth and mortality.

LiDAR Data and Calibration/Validation Plot Data
LiDAR data and calibration/validation plot data as described in Lamb et al. [44] were used in this study (Table 1).A private LiDAR contractor acquired high density LiDAR (6 pulses•m −2 , with up to 8 returns per pulse) in September 2011 for the northern portion (north of Highway 17, Figure 1b) of the District and in July 2013 for the southern portion.Scanning parameters described in Lamb et al. [44] were identical in both acquisitions.Ninety-five circular 400-m 2 plots were established in the summer of 2012 in the area flown in 2011 for use in calibrating LiDAR predictions (Figure 1b).An additional five plots were established in 2014 in the area flown in 2013 (Figure 1b).The plots were established in 41 black spruce, 45 white spruce, and 14 Norway spruce plantations to cover the range of inventory conditions observed across spruce plantations in the District [47].Plot centers were geo-referenced using an SX Blue II GPS system to obtain sub-meter accuracy.Species, status (live or dead), and DBH were measured for all trees with DBH ≥3 cm.Three randomly selected trees within each species and 5 cm DBH class were measured for height.For each plot, total and merchantable basal area (BA t , BA m ), gross total and merchantable volume (VOL t , VOL m ), top height (TopHt; average height of the 100 stems per hectare (4 stems per 400-m 2 plot) with the largest DBH), and merchantable quadratic mean diameter (QMD m ) were calculated using the formulae described in Lamb et al. [44], with 9.1 cm DBH as the lower merchantability limit.Volumes were calculated using taper equations from Li et al. [48] for commercial species and from Honer et al. [49] for non-commercial species.The calibration plots were used for validation as described in a subsequent section, however, two of the 100 calibration plots (one white spruce and one black spruce) partially fell outside of J.D. Irving, Limited forest resource inventory plantation polygon boundaries.Although these two calibration plots were included for calibration of LiDAR-derived forest inventory, they were not used for validation (n = 98).

LiDAR-Derived Forest Inventory Predictions
Following best practices for the use of LiDAR in area-based methods presented by Woods et al. [2] and White et al. [3], the LiDAR contractor clipped the point cloud data to correspond to the location and area of each circular 400-m 2 calibration ground plot.The LiDAR was flown to achieve a first return density of 6•m −2 with 8 returns per pulse, resulting in a nominal first return point density of 2400 and the potential for 19,200 returns per plot.Point cloud statistics were calculated for each plot from the clipped data using all LiDAR returns (i.e., no height threshold was used) with FUSION LiDAR Software [50].Statistics that were calculated included, but were not limited to, mean height of returns, standard deviation of return heights, return height percentiles, etc. [2,4].The LiDAR point cloud overlaying the District's spruce plantations was divided into 20 × 20 m grid cells and the same suite of point cloud statistics were calculated for each grid cell.Using R [51], the LiDAR contractor ran Random Forests TM [52] in regression mode to build ensemble decision trees to determine the statistical relationships between LiDAR point cloud statistics from the calibration plots and the corresponding forest inventory variables (BA t , BA m , VOL t , VOL m , TopHt, and QMD m ).Default parameters of 500 trees, one-third of the predictor variables used at each split, and a random two-thirds of the observations used to build any one tree were used for each Random Forests ensemble decision tree.The resulting decision trees were used to predict estimates of forest inventory at 20 × 20 m grid-cell resolution for all Black Brook District spruce plantations.These LiDAR-based forest inventory methods are used operationally in many jurisdictions in Canada and are based on more than 25 years of scientific research on the application of airborne laser scanning data to forest inventory [2][3][4].Upon completion of all inventory predictions by the contractor, the finalized LiDAR-derived forest inventory was delivered to J.D. Irving, Limited, and was then provided for our use in this study.Forest inventory predictions derived from point cloud statistics will herein be referred to as "LiDAR-derived".

Ancillary Forest Inventory Plot Data
The library of plots used in the plot-matching approach was described in detail by Lamb et al. [44] and only a summary will be given here.This included: (1) permanent fixed-area plots installed in spruce plantations by the New Brunswick Department of Energy and Resource Development (NB ERD) [53] and J.D. Irving, Limited (including MacLean et al. [46] and Omari and MacLean [54]), with all trees measured for DBH and height in each measurement; and (2) variable radius plantation temporary sample plots from the NB ERD Forest Development Survey database [55], which used a 3 m 2 •ha −1 basal area factor prism to sample trees for species and DBH and a 10 m 2 •ha −1 basal area factor prism to sample trees for height.Overall, this included over 65,000 variable-radius and permanent fixed-area plots installed in spruce plantations throughout New Brunswick.As with the LiDAR calibration/validation plots, stand-level inventory variables were calculated using the equations presented in Lamb et al. [44].
Since detailed species information was not available from the LiDAR-derived inventory and the only information available for each spruce plantation polygon was planted species, plots within the plot-matching library required species compositions representative of those observed throughout the District's spruce plantations.This was necessary, as the company's spruce plantations are intensively managed to retain high planted species percentages, compared to plantations on public lands, which can contain higher proportions of non-planted species (e.g., intolerant hardwoods and balsam fir (Abies balsamea (L.) Mill.)).Almost 400 Black Brook spruce plantation permanent sample plots with close to 1400 total measurements, as well as the LiDAR calibration plots, were used calculate the 95th percentile of species compositions by basal area.Only ancillary plots that fell within the range of observed species composition (i.e., ≥60% spruce, <40% balsam fir, <17% hardwood, and <5% non-commercial species) were included in the plot-matching library.This included 2983 black spruce, 1184 white spruce, 63 Norway spruce, and 216 red spruce dominated plantation plot measurements.Since white spruce and Norway spruce shared similar ranges of forest inventory and are both planted on similar site conditions, the 1184 white spruce plots were combined with the Norway spruce plots to increase the sample size of Norway spruce plots.Inventory variables were recalculated for these plots using Norway spruce height-diameter relationships and taper equations.This was deemed necessary, as over 270,000 LiDAR grid cells were located within Norway spruce plantations.In total, 5630 plot measurements were available for plot matching.The LiDAR calibration plots were not included in the plot-matching library, as these plots were to be used for validation, as described in a subsequent section.

Plot-Matching Approach
A key prerequisite for the successful application of the plot matching method was that plot measurements and LiDAR grid cells share the same planted/dominant spruce species.Planted species was identified for each LiDAR grid cell by intersecting the cell centroids with the J.D. Irving, Limited forest resource inventory polygons using ArcGIS (ESRI, Redlands, CA, USA), which contained planted species and year of establishment for each plantation polygon.This resulted in almost 2.1 million LiDAR grid cells in spruce plantations: 927,360 black spruce, 843,360 white spruce, 272,130 Norway spruce, and 41,810 red spruce.
A nearest-neighbor matching algorithm built within the FORUS program [33] was used to match LiDAR grid cells to one (k = 1) ancillary plot based on the smallest sum of squared difference between each inventory variable used, referred to as 'distance' (Equation ( 1)).The influence of different variable ranges (e.g., BA t vs. VOL t ) was removed from matching by normalizing each variable by dividing values by the maximum of observed plot values for each variable.Match distance was calculated [44] as: where w is variable weight, c is LiDAR cell value, p is ancillary plot value (i . . .n), for each variable v: planted species, BA t , BA m , VOL t , VOL m , TopHt, or QMD m .
Planted species was converted to a numeric value (e.g., black spruce = 1, white spruce = 2, etc.).To ensure that LiDAR grid cells and ancillary plot matches shared the same planted/dominant spruce, an arbitrarily high variable weight of 1000 was used for planted/dominant species, whereas all other variables (BA t , BA m , VOL t , VOL m , TopHt, and QMD m ) were given an equal weight of 1.This method of imputing tree-level inventory for LiDAR grid cells has been shown to result in highly correlated variable values (r = 0.91-0.99)and statistically equivalent basal area distributions (86% of the time (α = 0.05)) when compared with those measured in 98 validation plots [44].Variables extracted from matched plots are referred to herein as "plot-matched variables".

Forest Growth Model and Model Calibration
The growth model used to forecast forest inventory for LiDAR grid cells was Open Stand Model (OSM) [56].OSM is a tree-list growth simulation model calibrated for the Acadian Forest Region that predicts forest vegetation changes in response to stand development for a range of management types and site conditions.For each annual iteration of growth, forest inventory attributes for individual trees (DBH and height) are altered using species-specific DBH and height increment and tree survival equations.Species DBH and height measurements are used to calculate volume using regional taper equations [48,49].Growth rates are influenced by the Acadian biomass growth index (BGI) model, which takes into account climate, lithology, soils, and topography [57].The BGI is available as a 20 × 20 m spatial raster (http://www.forusresearch.com/bgi.php)and was used to determine BGI for each plot and grid cell location.
While OSM is calibrated for intensively managed plantations in the Acadian Forest, the Black Brook District is known to have above average productivity and plantation management practices compared to most other areas in this region, so it was necessary to ensure that growth derived from OSM would be representative of the District's plantations.To do this, the 399 permanent sample plots found in Black Brook spruce plantations were used to calibrate OSM.Plot measurements were used only if no thinning treatments occurred between two measurements to remove the influence of treatments on observed growth and yield.This resulted in 173 plots in plantations ranging from 13-45 years old, with 534 total measurements: 175 black spruce, 295 white spruce, 1 Norway spruce, and 63 red spruce.
Height increment calibration was done within the OSM software.Each of the 534 plot measurements was forecasted in OSM using its respective BGI to produce predicted height increments, and these were compared with observed increments.OSM contains methods for (1) validation of height and DBH increment models against local observations and (2) auto-adjustment of these models to correct bias through simple linear regression [33].Variables considered in model building included total plot basal area, DBH, and basal area of larger trees.Height increment models were built and used to calibrate growth for black spruce, white spruce, red spruce, and balsam fir.No height increment model was built for Norway spruce due to low sample size.DBH increment calibration models were not used, as they did not result in improved predictions of DBH increment.
In addition to height increment calibration, the need for mortality calibration was also assessed.The FORUS platform allows for the modification of background mortality rates by user-defined diameter classes, as well as the maximum relative density limit where self-thinning is to be applied.Plot-level observed and predicted annual mortality rates (stems•year −1 ) were calculated for five diameter classes (3-10 cm, 10.1-20 cm, 20.1-30 cm, 30.1-40 cm, and >40 cm) to compare mortality.Based on the predicted and observed mortality, mortality rates were reduced by 50% and 80% in the 10.1-20 cm and 20.1-30 cm DBH classes, respectively.Relative density where self-thinning was to be applied was increased from 85% to 100%.The adjustments to the mortality rates and self-thinning lines reduced the overprediction of mortality in the 10.1-20 cm and 20.1-30 cm DBH classes from 47% and 80% to only 11% and 12%, respectively.Mortality rates in all other DBH classes differed by no more than 6%.

Forecasting Forest Inventory for LiDAR Grid Cells
To obtain forecasted inventory for LiDAR grid cells, each plot used in matching was forecasted using OSM and a BGI value representative of the LiDAR grid cell to which it was matched.The BGI values for Black Brook spruce plantations ranged from 2251 kg•ha −1 •year −1 to 4894 kg•ha −1 •year −1 .To reduce the combinations of plot and BGI values, and thus the need for individual calculations, three BGI values were used for forecasts: 3000, 3750, and 4500 kg•ha −1 •year −1 .The 5630 inventory plots used in plot matching were forecasted for 10 years in 1-year time steps using each of these three BGI values.The BGI of each LiDAR grid cell was determined using ArcGIS (ESRI, Redlands, CA, USA) and the regional BGI raster (http://www.forusresearch.com)and rounded to the nearest of the three BGI values used for forecasting.Each LiDAR grid cell was then linked with its forecasted matched plot, creating 10-year forecasted BA t , BA m , VOL t , VOL m , basal area weighted average height (Lorey's height; LHt), total quadratic mean diameter (QMD t ), and live crown ratio for each grid cell.Mean annual height increment was calculated for each grid cell by dividing its LHt by age of plantation at year of forecast.

Plot-Level Comparison of Predicted Inventory Variable Increments
The 100 calibration/validation plots were used to test the accuracy of inventory variable increments predicted using tree lists imputed from plot matching.
Tree lists from the calibration/validation plots were input into OSM using the BGI for each plot location to obtain annual increments predicted from measured tree lists for six inventory variables (BA t , BA m , VOL t , VOL m , LHt, QMD t ).Since these plots were 400-m 2 circular areas, they covered more than one single square 400-m 2 LiDAR grid cell.Plot match-derived inventory variable increments were calculated for each plot by summing the product of the variable increment multiplied by the corresponding proportions of the grid cell intersecting a plot area.The annual increments of the six variables predicted from tree lists imputed by plot matching were plotted against annual increments predicted from actual, measured plot tree lists.Pearson correlation coefficient (r), mean error, and root-mean-square error (RMSE) were calculated for all six inventory variables.To compare across variables, mean error and RMSE were both calculated as a percentage by dividing by the mean of each respective variable and multiplying by 100.This comparison was done using 98 of the 100 calibration/validation plots, as two plots partially fell outside of spruce plantation polygon boundaries and were excluded from validation.

Plantation-Level Comparison of Harvest Volume and Forecasted Volume
Harvested volume data for 15 spruce plantation harvest blocks in the Black Brook District were compared with corresponding volumes forecast from tree lists imputed by plot matching.These blocks were located in six black spruce and nine white spruce plantations with areas ranging from 4 to 45 ha.All blocks were clearcut in 2016 with the same product specifications as used to calculate VOL m .Four blocks were located in the 2011 LiDAR acquisition area, while the remaining 11 were located in the 2013 LiDAR area, equating to 5-and 3-year plot match-derived forecasted VOL m , respectively.The 2016 plot match-derived VOL m was intersected with the block polygons using ArcGIS (ESRI, Redlands, CA, USA) to calculate forecasted VOL m and this was plotted against corresponding measured harvest volumes.In addition, VOL m derived from the original unforecasted LiDAR predictions were plotted against measured harvest volumes to assess improvements when using plot match-derived forecasted VOL m .Pearson correlation coefficients (r) and mean bias, both weighted by harvest block area, were calculated for each case.A paired t-test, weighted by block area, was used to determine if plot match-derived forecasted VOL m and LiDAR-derived VOL m were significantly different from measured harvest volumes.

Mapping Commercial Thinning Eligibility for Operational Harvest Planning
J.D. Irving, Limited commercial thinning eligibility rules were used to rank and map LiDAR grid cells that would qualify for commercial thinning in 2018 using plot match-derived forecasted inventory variables.These rules use "points" assigned to LiDAR-or plot match-derived inventory variables (i.e., VOL m , average merchantable tree volume (VOL m •stems per hectare −1 ), live crown ratio, and mean annual height increment (m•year −1 )) as well as minimum/maximum age and planted species derived from forest inventory polygons.The sum of these points was used to rank cells for commercial thinning eligibility.Using plot match-derived forecasted variable values for 2018 and forest resource inventory polygons, commercial thinning eligibility points were calculated for each grid cell.Cells were then categorized based on their commercial thinning eligibility as high ranking (points ≥ 9.5), medium ranking (7.5-9.4), and low ranking (0.1-7.4).High ranking cells are expected to have the greatest response to release and highest financial return, in comparison to medium and low ranking cells.
To demonstrate the process for prioritizing annual commercial thinning treatments using plot match-derived forecasted inventory, the commercial thinning ranking system and a measure of block-level competition were used to spatially schedule commercial thinning for 2018, 2019, and 2020.The objective was to target high ranking blocks that would provide the greatest financial return and prioritize blocks to capture density-dependent mortality.Using the dissolve tool in ArcGIS (ESRI, Redlands, CA, USA), all high-ranking cells sharing a border were dissolved together to create new contiguous block polygons.Block polygons less than 4 ha were removed from commercial thinning consideration, as this was the smallest harvest block size observed from the 2016 harvest data.To quantify block-level competition, basal area was used, as it is commonly used as a measure of competition [58,59].Mean BA t was calculated for each block polygon using the 2018 plot match-derived forecasted BA t .Block polygons were selected in order of decreasing mean BA t until the sum of total area selected was ≤3000 ha (average total annual area commercially thinned in the District since 2014 [60]).These polygons represented the proposed 2018 commercial thinning blocks.To select commercial thinning blocks for 2019 and account for proposed treatments in 2018, all 2019 plot match-derived grid cells that were coincident with the proposed 2018 commercial thinning blocks were removed using the erase tool in ArcGIS (ESRI, Redlands, CA, USA).Following the same process used for selecting 2018 commercial thinning blocks, ≤3000 ha of high ranking, high BA t areas were selected from the 2019 plot match-derived forecasted inventory.This same process was repeated for the 2020 plot match-derived forecasted inventory.This resulted in proposed 2018, 2019, and 2020 commercial thinning blocks that were mapped to show the distribution across the District.

Comparison of Inventory Variable Increments Predicted from Plot-Matched versus Field-Measured Tree Lists
Inventory variable increments predicted with OSM using tree lists imputed by plot matching were strongly correlated with those using measured tree lists from 98 calibration/validation plots.For the six inventory variables compared, r ranged from 0.75-0.86(Table 2 and Figure 2).Percentage mean error and percentage RMSE ranged from −3.6-0.0%and 12.8-36.6%for BA t , VOL t , VOL m , LHt, and QMD t , but were 11.3% and 49.0% for BA m , respectively.Table 2. Comparison between the annual increments of six inventory variables forecast for 5 years using Open Stand Model using tree lists imputed by plot matching versus forecasts using tree lists measured in 98 calibration/validation plots.Variables included basal area, merchantable basal area, volume, merchantable volume, Lorey's height, and quadratic mean diameter; comparison measures include Pearson correlation coefficient (r), mean error, and root-mean-square error (RMSE).

Inventory Variable r Mean Error Mean Error (%) RMSE RMSE (%)
Total basal area (m

Comparison of Harvested Volumes with Predicted Volumes Derived from Forecasted Plot Matches and LiDAR
Comparison of harvested volumes from 15 spruce plantations harvested in 2016 revealed that plot match-derived forecasted VOL m were highly accurate (Table 3 and Figure 3).Differences between plot match-derived forecasted VOL m and harvested volumes ranged from −12.0-14.4% (Table 3) with an overall mean difference of only 1.5% when weighted by block area.When year of LiDAR acquisition was considered, mean differences were −4.1% and 3.0% for harvest blocks in the 2011 and 2013 LiDAR acquisition areas, respectively.The greatest errors in this comparison resulted from the inaccuracy of the LiDAR-derived VOL m .For harvest blocks 11 and 12, with errors of 14.4% and 14.3%, respectively, LiDAR-derived VOL m was 11.8 and 7.0 m 3 •ha −1 greater, respectively, than volumes harvested 3 years after LiDAR-derived predictions.There was a stronger correlation and reduced mean bias between harvested volumes and plot match-derived forecasted VOL m (0.92 and 2.5 m 3 •ha −1 ) when compared to LiDAR-derived VOL m (0.83 and 11.0 m 3 •ha −1 ) (Figure 3).In addition, the results from the paired t-test indicated that plot match-derived forecasted VOL m was not significantly different from harvest volumes (p = 0.54), but that the unforecasted LiDAR-derived VOL m was significantly different (p = 0.032).

Mapping Commercial Thinning Ranking and Scheduling Annual Commercial Thinning Treatments for Spruce Plantations in the Black Brook District
Using the 2018 plot match-derived forecasted inventory, commercial thinning rankings were applied to individual grid cells using J.D. Irving, Limited commercial thinning eligibility rules.When grid cells were removed if they did not fall within the minimum and maximum age limits, 32,370 ha (809,250 cells) were "eligible" for commercial thinning (Figure 4).Ranking individual cells resulted in 9870 ha, 12,780 ha, and 9720 ha of high, medium, and low ranking area (Figure 4a).When bordering cells sharing the same commercial thinning rank were dissolved into larger polygons, the result was 49,760 individual polygons.This number was greatly reduced when minimum polygon size was considered, as there were only 2840 and 1470 polygons ≥1 ha and ≥4 ha, representing 27,920 ha and 25,160 ha, respectively.When annual commercial thinning was planned for high ranking blocks ≥4 ha using plot match-derived forecasted BA t as a competition index, the result was 2970 ha, 2990 ha, and 2970 ha for years 2018, 2019, and 2020, respectively (Figure 5).Minimum block size was 4 ha for all years, however, maximum size was 290 ha, 200 ha, and 100 ha, and mean size was 27 ha, 22 ha, and 19 ha for years 2018, 2019, and 2020, respectively.The number of blocks scheduled for years 2018, 2019, and 2020 was 112, 134, and 157, respectively.

Discussion
LiDAR-derived forest inventory introduced a major shift in how forest management is practiced in Canada [2,3] and abroad [61,62].With LiDAR, the planning of harvest and silviculture treatments no longer needs to be determined by photo-interpreted stand polygons and strata-based yield curves.Instead, within-and between-stand variability can be determined using this high-resolution inventory, allowing for prescription development prior to site visits.However, without a means of forecasting the inventory, forest managers would eventually be forced to return to stand-based yield curves for operational-, tactical-, and strategic-level planning between inventory refresh cycles.Using tree-list imputation methods [44], we have presented a framework for forecasting grid-cell inventory variables using imputed tree lists and a locally calibrated tree-list growth model.This extends the life of LiDAR-derived inventory and allows forest managers to continue high-resolution forest planning into the future.
Tree lists imputed for spruce plantation LiDAR grid cells from plot matches allowed for accurate inventory variable increment predictions when input into a locally calibrated growth model, as validated by comparing variable increments predicted using measured tree lists for 98 calibration/validation plots.Plot matches resulted in strong correlation, low mean bias, and low RMSE for BA t , BA m , VOL t , VOL m , LHt, and QMD t increment predictions, compared with increments predicted from measured trees.Much of the error associated with inventory variable increment predictions resulted from calibration/validation plots located in young plantations with merchantable values of 0. For example, of the 20 plots with plot match-derived BA m increments that differed from measured tree-list increments by more than ±0.5 m 2 •ha −1 •year −1 , 11 had LiDAR-derived merchantable values of 0. Without merchantable attributes, plot matches were based largely on BA t , VOL t , and TopHt.This introduced the possibility that matched plots may not be representative of actual conditions, as the combination of BA t , VOL t , and TopHt was not sufficient to accurately describe density or average diameter.The addition of QMD t as a plot-matched variable would ensure that plot matches are more representative of density and average diameters, however, this LiDAR-derived variable was not predicted by the LiDAR contractor so was not available for matching.Matching to grid cells in unmerchantable stands has additional challenges in that few plots were in young stands.Of the 5630 plot measurements used in matching, only 9% had no merchantable tree measurements, compared with 35% of all LiDAR grid cells with merchantable value predictions of 0 and 19% of LiDAR calibration plots.The lack of plot measurements in the plot-matching library with no merchantable tree measurements was likely the result of the forest inventory assessment design, which favors sampling merchantable stand types to assess current operable growing stock rather than unmerchantable stands.The relatively high percentage of LiDAR grid cells with merchantable value predictions of 0 is largely based on the management and age distribution of these plantations, which are intended for clearcut before the age of 40.Despite increased variance in the growth increment predictions of plots in young stands, predictions remained unbiased, suggesting that errors will average out when grid cells are aggregated.
When plantation harvest blocks were examined, plot match-derived forecasted VOL m showed strong agreement with harvested volumes, suggesting that forecasted plot matches may provide estimates of stand-level volume that are acceptable for operational planning.Although plot match-derived forecasted VOL m accuracy varied between blocks, overall, plot match-derived forecasted VOL m was only 1.5% greater than harvested volume and was found to not be significantly different.The slightly higher average VOL m may be the result of defects or cull not accounted for in the calculation of VOL m , however, this was beyond the scope of our study and was not assessed.The benefits of plot matching and the resulting forecasted VOL m were evident by the increased correlation and decreased bias when compared with 3-5-year-old LiDAR-derived VOL m .For fast growing stands, such as the spruce plantations in the Black Brook District, the addition of plot match-derived forecasted inventory with LiDAR-derived estimates presents a significant advantage for operational forest management.
The use of LiDAR-or plot match-derived inventory variables to spatially map areas of interest may be the most beneficial advantage of this high-resolution inventory.Quantifying variability in inventory variables is of greater importance than knowing stand-level averages, as variability directly impacts operational costs, response to treatment, wildlife habitat, etc. [2].This is particularly true for silvicultural treatments such as commercial thinning, which have very exacting stand-condition requirements for positive cost-benefit outcomes.Watt et al. [26] concluded that LiDAR-derived estimates of average tree height and base to live crown were sufficiently precise to delineate and schedule commercial thinning in Douglas-fir (Pseudotsuga menziesii (Mirb.)Franco) stands.However, as time since acquisition of LiDAR-derived inventory increases, it becomes less representative of actual conditions.We have demonstrated how plot match-derived inventory variables, forecast for 3-5 years with a tree-list growth model, and commercial thinning ranking rules can be used to delineate high-ranking grid cells beyond the year of LiDAR acquisition.In addition, we presented a novel method for planning annual commercial thinning treatments using the commercial thinning ranking and a measure of competition.Although BA t was used as the measure of competition, other competition indexes (e.g., stand density index or relative density [35]) can be calculated from imputed tree lists and used to prioritize thinning treatments.
Forecasting forest inventory is an important issue for all inventory methods.Although methods to obtain forecasted inventory at grid-cell resolution have been proposed [28][29][30], the variables that can be obtained from the growth models used in these studies are limited, as they did not use tree-level inventory.The advantage of the plot matching method over other LiDAR forecasting methods is the ability to generate tree lists and to calculate additional inventory variables that were not estimated from LiDAR, and to combine them to generate forecasts.Variables such as species-specific product volume [37], biomass and carbon [63], and wildlife habitat [64] have all been estimated using LiDAR, but each variable requires a unique set of calibration plots to cover the variability observed on the land base [3].Using the methods presented in this study, future research should focus on assessing the ability to impute such complex variables.
Despite the benefits that forecasted plot matches present, forecast inventory estimates are limited by several potential issues.The accuracies of plot match-derived inventory variables and subsequent forecasts are limited by the accuracy of the LiDAR-derived inventory [44].Even with accurate LiDAR-derived inventory, plot match-derived forecasted inventory may still be inaccurate for stand conditions not represented in the plot measurement database, so traditional age-based yield curves may be more accurate for rare forest cover types.Finally, the accuracy of the tree-list growth model used to forecast matches will influence results.Since our study area had a large number of permanent sample plot measurements available to calibrate height and mortality increments for OSM, this likely contributed to the accuracy of estimates of forecasted variable values.The fact that our study was conducted on even-aged spruce plantations also likely contributed to the high accuracy of estimates of forecasted variable values and should be taken into consideration when attempting to apply these methods to more complex stand types (e.g., uneven-aged mixedwood stands).

Conclusions
We imputed tree-level inventory data for 2.1 million LiDAR grid cells in spruce plantations using a nearest-neighbor matching algorithm and plot data from throughout the province of New Brunswick [44].This method uses the library of existing plot measurements and a small number of LiDAR-derived inventory variables to impute complex tree-level inventory data for each grid cell.Using a locally calibrated tree-list growth model and BGI value for the cell location, imputed tree lists were forecast 10 years in 1-year time steps to obtain forecasted inventory variables (i.e., BA t , VOL t , LHt, etc.) at 20 × 20 m grid-cell resolution.Increments of variables predicted using tree lists imputed from plot matching were strongly correlated with those predicted from measured tree lists, suggesting that imputed tree lists were sufficiently accurate for forecasting inventory variables.This was supported by the comparison of plot match-derived forecasted VOL m and harvested volume, which showed stronger correlation and lower bias when compared with 3-5-year-old LiDAR-derived VOL m .Using plot match-derived forecasted inventory variables (e.g., VOL m , mean annual height increment, live crown ratio, etc.) and commercial thinning eligibility rules, we demonstrated the ability to spatially plan annual silviculture treatments that were otherwise unobtainable with the older LiDAR-derived inventory.
The benefit of this approach is the use of a tree-list growth model to obtain high resolution forecasted inventory data.Typically, forecasting forest inventory for a land base at grid cell resolution with a tree-list growth model would be computationally infeasible.However, we have presented a novel and practical solution that allows forest managers to make detailed operational-, tactical-, and strategic-management plans, not otherwise available using traditional inventory methods.

Figure 1 .
Figure 1.Location of (a) Black Brook District in northwestern New Brunswick, and (b) light detection and ranging (LiDAR) calibration/validation plots in spruce plantations.LiDAR was acquired in 2011 and 2013 for areas north and south of Highway 17, respectively.

Figure 2 .
Figure 2. Comparison of the annual increments over 5 years predicted using Open Stand Model for six inventory variables (a-f): basal area, merchantable basal area, volume, merchantable volume, Lorey's height, and quadratic mean diameter for tree lists imputed by plot matching and measured in 98 calibration/validation plots.

Figure 3 .
Figure 3.Comparison of merchantable volume derived from tree lists imputed by plot matching and forecasted using Open Stand Model (green circles) and original unforecasted LiDAR-derived predictions (red squares) with actual harvest volumes from 15 harvest blocks harvested in 2016.Time elapsed between LiDAR acquisition and harvest ranged from 3-5 years with an average of 3.4 years.Volume difference between forecasted plot matches and original LiDAR ranged from −2.0-40.8m 3 •ha −1 with an average of 13.5 m 3 •ha −1 .Point sizes represent relative harvest block size (3.9-45.0ha).Bias is presented in variable units (m 3 •ha −1 ).

Figure 4 .
Figure 4. Commercial thinning (CT) ranking for LiDAR grid cells in (a) Black Brook District and (b) one sample stand using inventory variables derived from tree lists imputed by plot matching and forecasted to 2018 using Open Stand Model.Cell ranking was determined using J.D. Irving, Limited commercial thinning eligibility rules for first commercial thinning, which consider species, age, gross merchantable volume, mean tree gross merchantable volume, mean annual height increment, and live crown ratio.

Figure 5 .
Figure 5. Optimal year of first commercial thinning (CT) for 2018, 2019, and 2020 using inventory variables derived from tree lists imputed by plot matching and forecasted to 2018 using Open Stand Model.Proposed blocks include high ranking grid cells based on J.D. Irving, Limited commercial thinning eligibility rules as well as a block-level competition index (basal area).J.D. Irving, Limited commercial thinning eligibility rules consider species, age, gross merchantable volume, mean tree gross merchantable volume, mean annual height increment, and live crown ratio.Priority was placed on blocks with greatest total basal area to capture density dependent mortality.

Table 1 .
Summary of the inventory variable values from the 400-m 2 LiDAR calibration plots (n = 100). 2

Table 3 .
Actual harvested versus forecasted merchantable volume derived from tree lists imputed by plot matching using Open Stand Model for 15 blocks harvested in 2016 in Black Brook District.Harvest blocks located in the 2011 and 2013 LiDAR acquisition areas represent 5-and 3-year plot match-derived forecasted merchantable volume, respectively.Averages were weighted by harvest block area.All blocks were clearcut with the same product specifications as was used to calculate plot match-derived merchantable volume.