Use of Multi-Temporal LiDAR to Quantify Fertilization Effects on Stand Volume and Biomass in Late-Rotation Coastal Douglas-Fir Forests

: Forest fertilization is common in coastal British Columbia as a means to increase wood production and potentially enhance carbon sequestration. Generally, the effects of fertilization are determined by measuring sample plots pre- and post-treatment, resulting in fertilization effects being determined for a limited portion of the treatment area. Applications of remote sensing-based enhanced forest inventories have allowed for estimations to expand to the wider forested area. However, these applications have not focused on monitoring the effects of silvicultural treatments. The objective of this research was to examine if a multi-temporal application of the LiDAR area-based method can be used to detect the fertilization effects on volume, biomass, and height in a second-growth Douglas-ﬁr ( Pseudotsuga menziesii ) stand. The study area on Vancouver Island was fertilized in January 2007, and sample plots were established in 2011. LiDAR acquisitions were made in 2004, prior to fertilization, and in 2008, 2011, and 2016, covering both treated and untreated areas. A total of 29 paired LiDAR blocks, comprised of four 20 m resolution raster cells, were selected on either side of the fertilization boundary for analysis of the effects across several different stand types differing in the percentage of Douglas-ﬁr, site index, and age. Random forest (RF) plot-level models were developed to estimate total stem volume and total stem biomass for each year of LiDAR acquisition using an area-based approach. Plot level results showed an increase in stem volume by 13% fertilized over control from 2005 to 2011, which was similar to a 14% increase in above-ground carbon stocks estimated using a tree-ring stand reconstruction approach. Plot-level RF models showed R 2 values of 0.86 (volume) and 0.92 (biomass) with relative cross-validated root mean square errors of 12.5% (volume) and 11.9% (biomass). For both the sample plots and LiDAR blocks, statistical results indicated no signiﬁcant differences in volume or biomass between treatments. However, signiﬁcant differences in height increments were detected between treatments in LiDAR blocks. The results from this research highlight the promising potential for the use of enhanced forest inventory methods to rapidly expand the assessment of treatment effects beyond sample plots to the stand, block, or landscape level.


Introduction
Forest fertilization as a silvicultural practice has increased in recent decades due to its observed positive impacts on tree growth [1][2][3][4]. Forest fertilization is the process of using sample plots, the accuracy of estimates is highly sensitive to conditions such as the local disturbance history [28] and can result in the overestimation of changes due to location bias of the sample plots [29]. Therefore, the level of accuracy for the wider area of analysis is directly linked to the number and location of sampling sites.
A method that can be applied for monitoring forests over an entire stand is tracking the carbon exchanges with Eddy covariance flux towers [30]. These towers are a common method to monitor the net exchange of carbon dioxide between forest canopies and the atmosphere and use cumulative high-frequency measurements of wind speed, direction, light, and CO 2 concentration from meteorological equipment to determine net ecosystem carbon production (NEP) [31]. While these towers can monitor NEP over time for larger areas than plots, the footprint of the contributing canopy area can be complex to determine and interpret [32]. Additionally, analysis conducted using flux towers lacks the ability to compare fertilization treatments against controls and thus can only compare the preand post-fertilization periods of the treated stand [30]. While both the flux tower and sample plot methods offer advantages for forest monitoring, in the context of carbon project verification, the ability to monitor effects over a larger operational area with a sample plot design's versatility would be ideal.
Advancements in remote sensing technologies can leverage observations made from sample plots to build models that are applied to the broader area, providing estimates of attributes for the entire area of data coverage [33]. One of the leading modelling methods, called the Area-Based (AB) approach, utilizes a suite of predictor metrics calculated from LiDAR data to develop a statistical model at the plot level for a forest attribute of interest, such as volume or biomass [34]. This statistical model can then be applied to similar stand types across the broader area for which there is observed data collected [35]. The advantage of this method is that it relies on collecting the same sample plot data that is usually collected in fertilization studies but allows for evaluating effects across a broader range of locations or even entire stands. Research using the AB approach has shown the ability to estimate variables at the stand level, such as volume [33] and biomass [36], with a high degree of accuracy while also showing the local within-stand variability of the estimates. In general applications of the AB approach, sample plots are established in the region of interest while sampling across the range of stand structures, primarily species mixtures and site qualities, to create a generalized model [37]. If applying an established AB model or sample plot network to detect fertilization impacts, if the treatment is not explicitly captured in the sampling design, the analysis may not capture the fertilization impacts. While the collection of LiDAR data and the use of AB-style forest inventories are becoming more common in forest management, to our knowledge, these area-based methods have not yet been used in a multi-temporal approach to evaluate the impacts of fertilization.
The objective of this research was to examine if a multi-temporal application of the LiDAR area-based method can be used to detect the fertilization effects on height, volume, and biomass, in a second-growth Douglas-fir (Pseudotsuga menziesii) stand. Area-based method results for one, four, and nine years post-fertilization were compared against observed sample plot results, as well as results from tree-ring stand reconstruction and carbon budget model predictions.

Study Area
The research study area was located in the Oyster River drainage northwest of Courtenay on eastern Vancouver Island, British Columbia, Canada ( Figure 1). The study area, located within the dry Coastal Western Hemlock biogeoclimatic subzone [38], was part of the Fluxnet Canada research network's coastal BC station [31]. The second-growth stands in the study area were mostly composed of Douglas-fir (Pseudotsuga menziesii) with components of western hemlock (Tsuga heterophylla), western redcedar (Thuja plicata), and red alder (Alnus rubra) [39]. The planting of Douglas-fir seedlings in 1949 established the main stand (DF49), with other species coming from natural regeneration after harvesting (1937, 1938, and 1943) and broadcast burning operations (1939 and 1943) [40]. The 2004 stand density index was 998 stems ha −1 and the stand had an estimated site index of 34 m [22]. In January 2007, a late rotation fertilization treatment was applied aerially to a central portion of the study area ( Figure 1) at a rate of 200 kg urea N ha −1 , a typical operational practice in this region [22]. located within the dry Coastal Western Hemlock biogeoclimatic subzone [38], was part of the Fluxnet Canada research network's coastal BC station [31]. The second-growth stands in the study area were mostly composed of Douglas-fir (Pseudotsuga menziesii) with components of western hemlock (Tsuga heterophylla), western redcedar (Thuja plicata), and red alder (Alnus rubra) [39]. The planting of Douglas-fir seedlings in 1949 established the main stand (DF49), with other species coming from natural regeneration after harvesting (1937, 1938, and 1943) and broadcast burning operations (1939 and 1943) [40]. The 2004 stand density index was 998 stems ha −1 and the stand had an estimated site index of 34 m [22]. In January 2007, a late rotation fertilization treatment was applied aerially to a central portion of the study area ( Figure 1) at a rate of 200 kg urea N ha −1 , a typical operational practice in this region [22].

Plot Establishment and Measurement
In 2011, three pairs of control and fertilized 0.04 ha circular plots (11.28 m radius) were established in a forest area (DF49a) on either side of the 2007 fertilization boundary ( Figure 1). Species and DBH were recorded for all trees above 9.0 cm DBH, heights measured for a subset of these trees, and a sample of 3 crop trees in each plot felled and stem disks cut to measure ring widths and wood properties [22]. These tree measurements were used to calculate stem volume using equations from Kozak [41] and biomass using equations from Canada's National Forest Inventory (https://nfi.nfis.org/en/biomass_calc, accessed on 15 August 2020 [42]). In 2017, all six plots were remeasured to record species, DBH, height, and other attributes following Canada's National Forest Inventory [43] sampling guidelines. Since several crop trees were felled in 2011, there were large decreases in plot-level volume and biomass between 2011 and 2017. To account for this anomaly through the whole time series, the values of these felled trees were removed from the plot-level volume and biomass totals in each year and were only included for a separate plot-level analysis of volume and biomass increments.

Tree-Ring Data Collection and Measurement
In 2013, all trees in the plots were increment cored, and tree rings analyzed to convert yearly ring width into DBH estimates for each tree's lifespan following procedures of Metsaranta et al. [27]. The estimated DBH values were then used with the Chapman-Richards (CR) equation to estimate tree heights. The CR equation is a common method for estimating tree heights from diameter measurements [44], and parameters were estimated from the pre-harvest plot data from the nearby DF49 stand [27]. For each tree's lifespan, merchantable volume was estimated from the DBH and height using regional taper equations from Kozak [41]. Tree volumes were then summed to the plot level and used to develop merchantable volume (m 3 ha −1 year −1 ) over age yield curves [27].

Biomass Increments Using Tree-Ring Volume Reconstructions and CBM-CFS3
The developed tree-ring yield curves were used as inputs to the CBM-CFS3 which used them to estimate above-and below-ground biomass dynamics and inputs to detrital and soil carbon pools [12]. The model was run for each of the plots using the ring curves from 1990 through 2011. The model was also run in a default mode using regional yield curves [10] for each plot's forest stand type and site index, without or with the model's fertilization disturbance effect [17] included. Annual above-ground biomass carbon increments were output from 2001 to 2011, which covered five years pre-and post-fertilization, similar to the periods examined by [22]. Comparisons of control and fertilized plot increments were made using two-way analysis of variance (ANOVA) for each of the pre-and post-fertilization periods.

LiDAR Data Aquisition
LiDAR acquisitions and remote sensing data for the area were obtained in 2004 [45], 2008 [46], 2011 [39], and 2016 (S. McLennan-Timberwest, personal communication, 8 December 2017) with various point densities and scan angle parameters. LiDAR points from each year were classified as ground or non-ground by the vendors using various versions of the Terrascan software (Terrasolid, Helsinki, Finland). Quality control to remove any duplicated points, air hits, or other errors was done on all point clouds through manual inspection and use of the lidR package in R [47].

Sample Plot Data for Area-Based Models
Matching the same years of LiDAR acquisition, the tree-ring-derived DBH and height estimates described above were subset for each tree into the 2004, 2008, and 2011 seasons. To represent the 2016 growing season, the measured attributes from the 2017 sampling were used. Individual tree volumes were then estimated using these data with the regional taper equations developed by Kozak [41], similar to Metsaranta et al. [27]. For each tree, the proportion that was considered as the stump, main stem, and treetop was determined by jurisdictional timber criteria. For these years, the British Columbia merchantability limits of a 30 cm stump height and a 10 cm top diameter were used [48]. Each tree's total stem volume consisted only of the tree's stem portions (i.e., those parts of the tree not including tops, stumps, and branches). Stem biomass was calculated using tree species, DBH, and height estimates with the online National Forest Inventory individual tree biomass calculator (https://nfi.nfis.org/en/biomass_calc, accessed on 15 August 2020 [42]) and included stem wood values only and not branches, bark, or foliage.
Individual tree volume and biomass within each plot was summed, and plot-level volume (m 3 ha −1 ) and biomass (Mg ha −1 ) were calculated for two groups. The first group (all trees) consisted of volume and biomass estimates from all trees within the plots in 2004, 2008, and 2011 and was only used in a plot-level analysis of treatment effects. The second group (subset trees) excluded estimates from the 2011 felled trees through all years and was used further to train the area-based models. Cumulative volume and biomass increments from 2004 (2008-2004; 2011-2004; 2016-2004) were calculated for each plot in addition to the periodic annual increments (2008-2004/4; 2011-2008/3; 2016-2011/5). A mixed-effects model with a fixed main effect of treatment and a random effect of plot [49] was used to determine if there were significant differences in volume or biomass cumulative or periodic increments in 2008, 2011, or 2016.

LiDAR Data Processing
As the LiDAR data was collected over several years by various vendors, there were structural differences between the acquisitions relating to point densities and scan angle ranges, with the 2016 data having almost ten times the average point density than that of the 2004 data. To address these differences in point cloud structure, all point clouds were thinned to a maximum point density of 3 points m −2 , and scan angles were limited to ±20 • . These limits were chosen based on the most constraining 2004 data. Thinning was accomplished using a moving window of 25 m 2 , in which points were randomly removed from the local window until the target density was reached, generating more consistent thinning results over areas of scan line overlap compared to using a global thinning target [47]. For the excluded trees removed from the field data, stem mapping information was used to geo-reference the stem locations from the plot center. Using an individual tree segmentation method [50], LiDAR points belonging to removed trees were identified in the 2004, 2008, and 2011 data and had their heights above ground set to 0, simulating their exclusion.
Multiple LiDAR-based predictor metrics were generated over the entire study area using the thinned point clouds for each year and a 20 m resolution raster with a comparable area as the sample plots (0.04 ha). Common predictor metrics used for estimating volume and biomass attributes include a combination of height, density, canopy cover, and statistical metrics [51]. In total, 47 point cloud metrics were generated using the lidR package in R [47] (Table S1 in Supplementary Materials). Additionally, point clouds were clipped to the sample plots' extent, and the same predictor metrics were generated at the plot level for model training [51].
Next, a series of paired "LiDAR blocks" was created to examine fertilization effects over the wider forested area around the fertilization boundary ( Figure 2). These blocks were created with a 20 m × 20 m fishnet grid aligned with the LiDAR metrics. This grid was intersected with forest inventory polygons and organized into analysis units (AUs) based on dominant species, stand age, and site index ( Table 1). Blocks of four cells were created on either side of the fertilization boundary with a 20 m buffer on the boundary's control side and a 10 m buffer on the fertilized side ( Figure 2). Pairs of blocks were selected if they shared a similar AU, and the set of blocks in the DF49a area were adjusted to ensure they contained the sample plots ( Figure 2). Ensuring the pairs were representative of the same area, the coefficient of variation (CoV) was calculated using the 95th percentile of LiDAR height, the mean LiDAR height, and the LiDAR point density for the four cells in a block as well as the eight cells that make up a pair. Groups of cells or block pairs with a CoV greater than or equal to 20% were removed. Additionally, pairs that belonged to AUs with less than three total pairs were also removed. In total, 29 paired blocks were selected across the four AUs (Table 1).
CoV greater than or equal to 20% were removed. Additionally, pairs that belonged to AUs with less than three total pairs were also removed. In total, 29 paired blocks were selected across the four AUs (Table 1).

Area-Based Modelling
Non-parametric random forest (RF) regression models were estimated using the plotlevel volume and biomass with the LiDAR predictor metrics. A Boruta variable selection method was employed to reduce the number of metrics used in model fitting. The Boruta method for feature selection is popularly used with machine learning applications in genetics, health sciences, and other disciplines where the number of predictor metrics is often much larger than the number of observations [52][53][54]. This method copies each metric and randomly permutes the data creating a set of noise metrics. The combined set of original and noise metrics are then used to fit many RF models where each metric's importance score is recorded. Metrics with an overall importance distribution significantly greater than the highest scoring noise variable make up the final selected metrics [55].

Area-Based Modelling
Non-parametric random forest (RF) regression models were estimated using the plotlevel volume and biomass with the LiDAR predictor metrics. A Boruta variable selection method was employed to reduce the number of metrics used in model fitting. The Boruta method for feature selection is popularly used with machine learning applications in genetics, health sciences, and other disciplines where the number of predictor metrics is often much larger than the number of observations [52][53][54]. This method copies each metric and randomly permutes the data creating a set of noise metrics. The combined set of original and noise metrics are then used to fit many RF models where each metric's importance score is recorded. Metrics with an overall importance distribution significantly greater than the highest scoring noise variable make up the final selected metrics [55]. Next, with the subset of metrics selected, RF models were trained to estimate total stem volume and total stem biomass using the field data in each year. The RF models' training was done using the caret and random forest packages in R [56,57]. A five-times repeated 10-fold cross-validation method was used while fitting 1000 regression trees [58]. This training method resulted in the removal of ten observations for each of the five repetitions resulting in 50 independent observations to assess model performance. Model performance was assessed using the independent observations and calculating the absolute and relative root mean square error (RMSE) and bias as described by White et al. [58]. The RF model for both the volume and biomass attributes was then applied across the study area using the suite of raster-based metrics in order to produce estimation surfaces for each year. Each year's estimated volume and biomass values were extracted from the surfaces of each cell of the LiDAR blocks, and block-level estimates for each year were calculated as the average value of the four cells in each treatment. As an additional check on model performance, the 95th percentile of LiDAR height in each cell was also extracted for increment analysis. Cumulative and periodic annual increments for volume, biomass, and height were calculated for the block-level estimates, similar to the sample plots. The block-level increments were used in mixed-effects models to test the significance of the fixed treatment effect with a LiDAR block's random effect. Tests were done each year for increments for AUs individually and combined across all AUs. Results from these tests indicated if the treatments in the LiDAR blocks produced a similar fertilization effect over controls as in the sample plots. The similarity of increment magnitudes between the sample plots and LiDAR blocks was assessed through visual comparison.

Sample Plots
In 2004, prior to the fertilization, average volumes in control plots (all trees) were slightly lower than in fertilized plots ( Table 2). Cumulative volume increments for both the all trees and subset trees groups were not significantly different between treatments in all years (p-values > 0.3, Table S2). Weighted by their 2004 volumes, control plots (subset trees) had an average cumulative increment from the 2004 volume ratio of 0.076, 0.175, and 0.320 by 2008, 2011, and 2016, respectively. In the fertilized plots, the weighted average ratio was 0.076, 0.181, and 0.303 by 2008, 2011, and 2016, respectively ( Table 2). By 2011 the un-weighted average cumulative volume increment in fertilized plots was 13% greater than in control plots. However, by 2016 the un-weighted average volume increment was only 3% greater in fertilized plots than in control plots ( Table 2). Periodic annual volume increments were not significantly different between treatments in any period (p-values > 0.3, Table S3).  (Table S4). In the 2012-2016 period, both treatments' overall rate decreased, with control plots having an average rate of 13.0 m 3 ha −1 yr −1 compared to a rate of 11.9 m 3 ha −1 yr −1 for the same period in the fertilized plots (Table S4) (Table 2). The biomass differences between treatment plots in 2004 were smaller than the volume differences. However, the control plots (all trees) had higher biomass values than the fertilized plots (Table 3). Similar to the plot-level volume increments, differences in cumulative biomass increments for both the all trees and subset trees groups between treatments were non-significant in all years (p-values > 0.3, Table S2). Fertilized plots continuously saw higher biomass gains with weighted ratios of 0.094, 0.229, and 0.372 in 2008, 2011, and 2016, respectively. In the same periods, control plots saw weighted ratios of 0.092, 0.216, and 0.293 (Table 3), respectively. By 2011 the un-weighted average cumulative biomass increment in fertilized plots was only 0.3% greater than in control plots. However, by 2016 the un-weighted average biomass increment was 20% greater in fertilized plots than in control plots (Table 3). For the period between 2005 and 2008, the periodic annual biomass increment in control plots was 3.0 Mg ha −1 yr −1 compared to an average rate of 2.8 Mg ha −1 yr −1 over the same period in fertilized plots (Table S5) (Table S5). However, similar to the volume periodic annual increments, there were no significant differences in rates between treatments over any period (p-values > 0.3, Table S3).

CBM-CFS3 Plot Biomass Increments Using Tree-Ring and Default Growth Curves
Annual biomass carbon increments (g C m −2 yr −1 ) using tree-ring curves varied more among control than treated plots across the 11-year period ( Figure 3). Prior to fertilization, annual biomass carbon increment in control and treated plots increased over time, and while the 5-year treated (242 SE 8) and control (269 SE 10) plot means were significantly different (p = 0.035), the 3-year means before fertilization were not (p = 0.122). Post-fertilization biomass increments in treated plots increased above control plots, with the 5-year means being significantly (p = 0.002) greater in treated (308 SE 9) than control plots (266 SE 10), a 16% increase. Cumulative biomass carbon increments (Mg C ha −1 ) from 2004 to 2008 and 2004 to 2011 were slightly greater in treated (11.8 SD 0.1, 21.8 SD 0.1) than control (11.0 SD 0.6, 19.2 SD 1.4) plots. Annual biomass carbon increments using default yield curves pre-fertilization ranged from 320 to 400 and were higher in plots to be fertilized (456 F) than control plots (123 C) (Figure 3). Post-fertilization, the model fertilizer effect predicted a 22% increase in 2007 and for the entire five-year period (Figure 3). means being significantly (p = 0.002) greater in treated (308 SE 9) than control plo SE 10), a 16% increase. Cumulative biomass carbon increments (Mg C ha −1 ) from 2008 and 2004 to 2011 were slightly greater in treated (11.8 SD 0.1, 21.8 SD 0.1) than (11.0 SD 0.6, 19.2 SD 1.4) plots. Annual biomass carbon increments using defau curves pre-fertilization ranged from 320 to 400 and were higher in plots to be fe (456 F) than control plots (123 C) (Figure 3). Post-fertilization, the model fertilize predicted a 22% increase in 2007 and for the entire five-year period (Figure 3).

Area-Based Models
The Boruta variable selection method results are shown in Table 4. The table p the metrics that had higher importance scores than the noise variables for each mod method took 3512 iterations for the volume model and resulted in the selection o DAR metrics, which consisted of both height and cover metrics ( Table 4). The method for the biomass model ran for the set maximum of 6000 iterations and s eleven total metrics. These metrics included similar height and cover metrics as t ume model and included several statistical metrics (Table 4).

Area-Based Models
The Boruta variable selection method results are shown in Table 4. The table provides the metrics that had higher importance scores than the noise variables for each model. The method took 3512 iterations for the volume model and resulted in the selection of six LiDAR metrics, which consisted of both height and cover metrics ( Table 4). The Boruta method for the biomass model ran for the set maximum of 6000 iterations and selected eleven total metrics. These metrics included similar height and cover metrics as the volume model and included several statistical metrics ( Table 4).
The fitted random forest (RF) volume model produced an R 2 value of 0.86, with a cross-validated RMSE of 67.33 m 3 ha −1 (12.55%) and a bias of 4.73 m 3 ha −1 (0.88%). Comparing the observed stem volumes to the modelled stem volumes showed a slight variable bias with no apparent patterns, except for the continual overestimating of volume in plot 8.3 ( Figure 4A). The fitted biomass RF model also reported a high R 2 value of 0.92 with a cross-validated RMSE of 17.15 Mg ha −1 (11.91%) and a bias of 0.42 Mg ha −1 (0.29%). Similar to the volume model, the biomass model showed some variable bias. However, the consistent over-predictions in plot 8.3 were not as large as in the volume model ( Figure 4B).  Table 4. LiDAR metrics selected and relative importance ranking for the volume (V, Vol) or biomass (B, Bio) models using Boruta variable selection method.   Figure 5A, and individual AU results are shown in Figure S1. Similar to the sample plots, across all AUs,  Figure 5A, and individual AU results are shown in Figure S1. Similar to the sample plots, across all AUs, there were no significant differences in volume increments in any year (p-values > 0.1, Table S6). Relative to their 2004 volumes, control blocks saw an average increase of 0.5, 0.9, and 6.  Table S7). Periodic annual increments averaged across all AUs were also not significantly different between treatment blocks during any period (p-values > 0.2, Table S8). However, for the AU that contained the sample plots (AU 146) for the period between 2009 and 2011, the fertilized blocks (2.6 m 3 ha −1 yr −1 ) increased at a significantly greater (p = 0.03) rate than in the control blocks (−0.7 m 3 ha 1 yr −1 ) (Table S9). time periods, fertilized blocks saw average increases of 2.0, 1.9, and 5.2 percent ov 2004 values ( Figure 5B). Similar to the volume results, there were no significant diff between treatments when increments were averaged across all AUs (Tables S6 a although the same individual AUs showed significant treatment differences in cum biomass increment (AU 445, Table S7) and periodic annual biomass increment (  Table S8).   Figure 5B). Similar to the volume results, there were no significant differences between treatments when increments were averaged across all AUs (Tables S6  and S8), although the same individual AUs showed significant treatment differences in cumulative biomass increment (AU 445, Table S7) and periodic annual biomass increment (AU146 , Table S8).

Predictor Class
In  Figure 6). Neither the average cumulative increment nor the average periodic annual increment across all AUs was significantly different (p = 0.08) between treatments by 2008 (Table S10). However, by 2011 the average cumulative height increment (2.61 m, SE 0.07) in fertilized blocks was significantly higher (p = 0.02) than in controls (2.39 m, SE 0.06), a 9% gain. Continuing this trend, by 2016, the average cumulative height increment (3.97 m SE 0.08) in fertilized blocks was significantly higher (p = 0.02) than in controls (3.56 m SE 0.06), for a 12% gain. Similar to the rates between 2005 and 2008, the average periodic annual increment between 2009 and 2011 was not significantly different (p = 0.05) between treatment blocks; however, for the period between 2012 and 2016, the rate in fertilized blocks (0.27 m yr −1 , SE 0.01) was significantly higher (p = 0.02) than in the control blocks (0.23 m yr −1 , SE 0.01). Height increment results for individual AUs can be seen in Figure S2 and Table S11. (0.27 m yr −1 , SE 0.01) was significantly higher (p = 0.02) than in the control blocks (0.23 m yr −1 , SE 0.01). Height increment results for individual AUs can be seen in Figure S2 and Table S11.

Discussion
This study investigated if area-based (AB) methods could be used with multi-temporal LiDAR data to detect fertilization effects across treated stands. Overall, the results show promise for LiDAR AB methods to be used to detect treatment impacts on stand volume, biomass, and height.

Method Comparisons
In the sample plots, the exclusion of felled trees appeared to have minimal impact on the analysis of treatment effects, as the all trees group showed from 2005 to 2011 that fertilized plots had an average un-weighted cumulative volume increment that was 13% greater than in control plots, which was similar to those seen in the subset trees group (Table 5). Weighted comparisons can be seen in Table S12. The differences between treatments had more variation in the biomass values whereby 2011, the subset trees group showed the fertilized cumulative biomass increments being less than the controls (Table  5). However, the differences between treatments were non-significant. Interestingly, CBM-CFS3 estimates using the tree-ring stand reconstruction approach showed that between 2005 and 2011, fertilized plots saw an increase in above-ground biomass (AGB) carbon stocks, compared to controls, that was more like the plot volume results than biomass values (Table 5). This difference in relative biomass increase suggests there may be differences between the biomass equations used in CBM-CFS3 and those used in the NFI online biomass calculator, or this could be due to the differences between AGB and stem wood biomass.
Plot-level periodic annual volume increments (subset trees) between 2005 and 2008 were 9% greater in fertilized over control plots, which increased to 16% greater in fertilized over control plots between 2009 and 2011; however, these differences were not significant (Tables S4 and 8). When examining the felled crop trees, Filipescu et al. [22] found that in fertilized plots, the volume 5-year average annual increment post-fertilization (2007-2011) was significantly greater than the average annual increment in control plots (Table 5). Average annual increments in AGB carbon from the stand reconstruction approach were similar, but slightly lower for the same post-treatment period as Filipescu et al. [22] (Table 5). For the plot-level biomass, periodic annual increments (subset trees) between 2005 and 2008 were −7% less in fertilized over control plots, which increased to 4% greater in fertilized over control plots between 2009 and 2011; however, these differences were not significant (Tables S5 and 8).

Discussion
This study investigated if area-based (AB) methods could be used with multi-temporal LiDAR data to detect fertilization effects across treated stands. Overall, the results show promise for LiDAR AB methods to be used to detect treatment impacts on stand volume, biomass, and height.

Method Comparisons
In the sample plots, the exclusion of felled trees appeared to have minimal impact on the analysis of treatment effects, as the all trees group showed from 2005 to 2011 that fertilized plots had an average un-weighted cumulative volume increment that was 13% greater than in control plots, which was similar to those seen in the subset trees group (Table 5). Weighted comparisons can be seen in Table S12. The differences between treatments had more variation in the biomass values whereby 2011, the subset trees group showed the fertilized cumulative biomass increments being less than the controls (Table 5). However, the differences between treatments were non-significant. Interestingly, CBM-CFS3 estimates using the tree-ring stand reconstruction approach showed that between 2005 and 2011, fertilized plots saw an increase in above-ground biomass (AGB) carbon stocks, compared to controls, that was more like the plot volume results than biomass values (Table 5). This difference in relative biomass increase suggests there may be differences between the biomass equations used in CBM-CFS3 and those used in the NFI online biomass calculator, or this could be due to the differences between AGB and stem wood biomass.
Plot-level periodic annual volume increments (subset trees) between 2005 and 2008 were 9% greater in fertilized over control plots, which increased to 16% greater in fertilized over control plots between 2009 and 2011; however, these differences were not significant (Tables S4 and S8). When examining the felled crop trees, Filipescu et al. [22] found that in fertilized plots, the volume 5-year average annual increment post-fertilization (2007)(2008)(2009)(2010)(2011) was significantly greater than the average annual increment in control plots (Table 5). Average annual increments in AGB carbon from the stand reconstruction approach were similar, but slightly lower for the same post-treatment period as Filipescu et al. [22] (Table 5). For the plot-level biomass, periodic annual increments (subset trees) between 2005 and 2008 were −7% less in fertilized over control plots, which increased to 4% greater in fertilized over control plots between 2009 and 2011; however, these differences were not significant (Tables S5 and S8). in pre-and post-treatment periods for stem volume (m 3 ha −1 ), stem biomass (Mg ha −1 ), and height (m) cumulative increment (gain) for plot subset trees and LiDAR blocks, and above-ground biomass carbon gain (Mg C ha −1 ) from a model reconstruction approach using plot tree-ring growth curves. Comparisons are also shown for periodic annual increment (PAI) of above-ground biomass carbon (Mg C ha −1 yr −1 ) from tree-ring curves for the same periods as volume growth (m 3 ha −1 yr −1 ) in [22]. Results from a flux tower in the main DF49 stand indicated that net ecosystem production (NEP) ranged from 268 to 410 (g C m −2 yr −1 ) from 1998 to 2006 due to interannual variability in weather [28]. For the 4-year post-fertilization period, the mean NEP from 2007 to 2010 (561 SD 38) increased by 64% compared to NEP for the 4-year pre-fertilization period, 2003-2006. However, Jassal et al. [30] found weather conditions differed between the periods and derived an empirical model for annual NEP and weather (annual air temperature, solar radiation, and 0-30 cm water content), which gave a lower mean estimate of NEP without fertilization from 2007 to 2010 (427 SD32), with fertilized being 31% greater than unfertilized. This was comparable to the fertilization effect noted by Filipescu et al. [22] (34%) but greater than that estimated using CBM-CFS3 with the tree-ring-based yield curves (14%) for the DF49a area. However, the effects noted by Filipescu et al. [22] were based on the extrapolation of crop tree increments to the rest of the plot, while the tree-ring yield curves were based on the individual tree increments.

Source
In the LiDAR blocks, cumulative volume (47.7 m 3 ha −1 F, 38.7 m 3 ha −1 C) and biomass Despite the differences in the magnitude of increments, the differences between treatments were largely non-significant in both the plot and LiDAR block-level results. However, the height increments in the LiDAR blocks showed significant differences between treatments in all years, other than 2008. Hilker et al. [46] found there were no noticeable differences in canopy height between treatments around the DF49 area in 2008. However, since the forests were fertilized in 2007, the height response might not be observed after one growing season. Littke et al. [59] found in paired single-tree fertilizations that fertilized tree height increments were greatest from 0 to 2 years and less from 2 to 4 years. They also reported that fertilized trees might show a basal area increment with little or no height increment. Hudak et al. [60] found using a bi-temporal LiDAR approach with random forest modelling that AGB increased at an average rate of 4.1 Mg ha −1 yr −1 over the six years between acquisitions in a complex conifer forest, which is much higher than increment rates seen at the LiDAR blocks in this study and more similar to the rates observed in the sample plots.

Sources of Uncertainty
While the results of this study clearly show the variation of treatment effects across methods and the treatment area, it is important to note that these results (as with all studies examining fertilization effects) are subject to the limitations of the applied methods. First, uncertainties exist related to the use of tree-ring analysis to produce plot or stand-level biomass estimates, including uncertainties introduced from the subsampling of trees within a plot and those from the use of allometric equations [61]. In this study, the influence of subsampling uncertainty was eliminated by utilizing the tree-ring analysis from all trees within the plots [27]. Uncertainty from allometric equations was also minimized in this study by using volume allometric equations that were species-and biogeoclimatic zonespecific [41] with parameters fit using plots from a harvested stand within the treatment area [27]. For the biomass calculations, regional equations were also used that specified the province (BC) and the ecozone (Pacific) of the data.
Another source of uncertainty comes from the removal of sampled trees from the LiDAR data. Stem mapping data often contain inaccuracies in absolute positioning due to GPS error [62], and similarly, LiDAR-identified treetops could have discrepancies in location between years due to differences in instrumentation, GPS error, and environmental conditions [51]. To limit this uncertainty, trees that were removed were manually matched with segmented tree crowns using both the location and height of field trees to identify matches from the segmented LiDAR point cloud. Finally, uncertainty was introduced into the analysis through the use of random forest (RF) models, which are limited by the inability to extrapolate beyond the set of training data [51]. The average 2004 values in LiDAR blocks were larger than the maximum of the 2004 values observed in sample plots, and suggest that the 2004 estimates approached the lower limit of the models. This would create lower yearly increments than reality if the actual 2004 values in these blocks were below the model's saturation point.
This study suggests that there is potential for the AB method to be used in a verification context. However, some challenges remain in optimizing the sampling design to describe changes over the entire treatment area. The sampling intensity of the sample plots used in this study appears to be appropriate given the low RMSE (12.55%/11.91%) and bias (0.88%/0.29%) values of the volume model, which showed a similar goodness of fit as other RF volume models applied in complex coniferous forests (RMSE of 26% [63]; RMSE of 33.24% [58]). However, in relation to the size of the plot level increments between years, the RMSE values were still high, and an increase in sample intensity might reduce the errors to an appropriate range for the increments. Additionally, the spatial distribution of the sample plots in future applications should be extended to include more of the variation within the treatment area. If available, plot locations could be determined by analyzing and stratifying the LiDAR point cloud structure utilizing a principal components analysis of LiDAR metrics [64].
As a method to augment sample plot estimates of treatment impacts, the AB approach offers the potential to increase the assessment's accuracy level and the spatial extent of the analysis. While more research is needed to determine the optimal sampling intensity and distribution for operational application, the fusion of the traditional sample design with a LiDAR area-based approach could offer the optimal design for a carbon verification context. As verification projects often balance the need for accuracy with project cost [21], the fusion of these two methods provides the flexibility to leverage plot-level observations to gain insights into the wider treatment area. This method should also become more accessible, particularly with the increased availability of multiple LiDAR acquisitions and the use of UAV-based LiDAR systems [65].

Conclusions
This study demonstrated that the impacts of fertilization observed in a ground plot sampling design do not necessarily translate across spatial scales. This was due to the differences between fertilization effects experienced in sample plots and those found in the broader treatment area. However, in terms of utilization for carbon project verification, the fusion of traditional sample plot methods and the LiDAR AB approach shows promise for describing fertilization changes across an entire operational treatment area. Future research on the fusion of these methods should focus on optimizing the sampling intensity and spatial distribution of the sampling design to balance a high level of accuracy with sampling costs.
Supplementary Materials: The following are available online: https://www.mdpi.com/article/10 .3390/f12050517/s1: Table S1. LiDAR predictor metrics used and selected by the Boruta variable selection method; Table S2. Mixed-effect model F and p-value results between treatments for volume and biomass cumulative increments with subset plot trees or all plot trees; Table S3. Mixed-effect model F and p-value results between treatments for volume and biomass periodic annual increments with all plot trees or subset plot trees; Table S4. Plot-level stem volumes (m 3 ha −1 ), periodic annual volume increments (m 3 ha −1 yr −1 ), and control (C) or fertilized (F) treatment averages for the subset trees in the plots, as well as all trees; Table S5. Plot-level stem biomass (Mg ha −1 ), periodic annual biomass increments (Mg ha −1 yr −1 ), and control (C) or fertilized (F) treatment averages for the subset trees in the plots, as well as all trees; Table S6. Mixed-effect model F and p-value results between treatments for volume and biomass increments in the LiDAR blocks across all AUs; Table S7. Mixed-effect model F and p-value results between treatments for volume and biomass increments in the LiDAR blocks for individual AUs; Table S8. Mixed-effect model F and p-value results between treatments for volume and biomass periodic annual increments in the LiDAR blocks across all AUs; Table S9. Mixed-effect model F and p-value results between treatments for volume and biomass periodic annual increments in the LiDAR blocks for individual AUs; Table S10. Mixed-effect model F and p-value results between treatments for 95th percentile of LiDAR height increments and periodic annual increments in the LiDAR blocks across all AUs; Table S11. Mixed-effect model F and p-value results between treatments for 95th percentile of LiDAR height increments and periodic annual increments in the LiDAR blocks for individual AUs; Figure S1. LiDAR block-level average (line) and standard error (shadow) of stem volume (left) and biomass (right) increments for control blocks