Estimating timber volume loss due to storm damage in Carinthia, Austria, using ALS/TLS and spatial regression models

A spatial regression model framework is presented to predict growing stock volume loss due to storm Adrian which caused heavy forest damage in the upper Gail valley in Carinthia, Austria, in October 2018. Model parameters were estimated using growing stock volume measured with a terrestrial laser scanner on 62 sample plots distributed across five sub-regions. Predictor variables were derived from high resolution vegetation height measurements collected during an airborne laser scanning campaign. Non-spatial and spatial candidate models were proposed and assessed based on fit to observed data and out-of-sample prediction. Spatial Gaussian processes associated model intercepts and regression coefficients were used to capture spatial dependence. Results show a spatially-varying coefficient model, which allowed the intercept and regression coefficients to vary spatially, yielded the best fit and prediction. Two approaches were considered for prediction over blowdown areas: 1) an areal approach that viewed each blowdown as a single prediction unit indexed by its centroid; and 2) a block approach where each blowdown was partitioned into smaller prediction units to better align with sample plots' spatial support. Joint prediction was used to acknowledge spatial dependence among block units. Results demonstrated the block approach is preferable as it mitigated change-of-support issues encountered in the areal approach. Despite the small sample size, predictions for 55% of the total 564 blowdown areas, accounting for 93% of the total loss, had a coefficient of variation less than 25%. Key advantages of the proposed regression framework are the ability to quantify uncertainty in spatial covariance parameters, propagate parameter uncertainty through to prediction, and provide statistically valid prediction point and interval estimates for individual blowdowns and collections of blowdowns.


Introduction
Analyses by Schelhaas et al. (2003) showed that storms were responsible for more than half the total damages in European forests for the period 1950-2000. The majority of this forest damage occurred in the Alpine zone of mountainous regions. In Austria, for the period 2002-2010, storms damaged 3.1 million m 3 annually (Thom et al., 2013), representing 0.26 % of the total growing stock as well as 12 % of the total annual fellings.
In Switzerland, the storm damage was 17 and 22 times higher in the period  than in the two preceding 50 years periods (Usbeck et al., 2010).
According to Usbeck et al. (2010), the possible explanations for such an increasing trend were manifold and include increased growing stocks, enlarged forested area, milder winters and hence tendency for wet and unfrozen soils, and higher recorded maximum gust wind speeds. Since 1990, 85 % of all primary damage in the Western, Central and Northern European regions was caused by catastrophic storms with maximum gust wind speed between 50-60 ms -1 (Gregow et al., 2017). The relevance of storm damage is likely to increase, as simulations with climate scenario data and forest stand projections suggested a higher probability of exceeding the critical wind speed and, hence, wind throw events (Olofsson and Blennow, 2008;Blennow et al., 2010). Such a trend could negatively impact the carbon balance, especially in Western and Central European regions (Lindroth et al., 2009).
Wind throw events have direct and indirect impact on long-term sustained yield planning. A common indirect cost stems from the fact that storm-felled trees provided a surplus of breeding material for bark beetles and promoted their rapid population increase causing extra timber losses in the subsequent years (Marini et al., 2016). Hence, it is recommended that fallen trees be removed within two years of the disturbance (Groot et al., 2018). Following a wind throw event and prior to harvesting the storm-felled trees, strategic salvage harvest planning is needed (Jirikowski and Pröll, 2003). Implementing such plans requires unplanned/unbudgeted logging road and site development efforts as well as interaction with external harvesting, transportation, and processing industry.
To inform strategic salvage harvest planning, a rapid and accurate estimate of the spatial extent and local severity of damage is needed. The task is to provide such estimates for individual and collections of discrete blowdowns. In most situations, a paucity of existing field inventory data within or adjacent to blowdowns precludes designbased estimation. Rather, field data from an existing set of sparsely sampled inventory plots or a small purposively selected set of sample plots, can be coupled with auxiliary data using a model to yield viable estimates. In such settings, the forest response variable measured on sample plots is regressed against meaningful auxiliary data. Commonly, such auxiliary data come as remotely sensed variables collected via satellite, aircraft, or unmanned aerial vehicle (UAV) based sensors. Remotely sensed data are routinely used to support forest inventories and have been shown to increase accuracy and precision of the estimates, and reduce field data collection efforts; see Köhl et al. (2006) and other references herein.
In practice, different techniques are used to predict forest variables using field measurements and remote sensing data. Nonparametric imputation methods, using some flavor of nearest neighbor (NN) interpolation (Moeur and Stage, 1995), achieved robust prediction for possibly multivariate response variable vectors (Tomppo and Halme, 2004;LeMay and Temesgen, 2005;Maltamo et al., 2006;Packalén andMaltamo, 2007). Hudak et al. (2008) demonstrated the NN prediction accuracy can be enhanced through resampling and classifications with "random forests" (Breiman, 2001). A key shortcoming of NN imputation methods, however, is the lack of statistically robust variance estimates, with the exception of some approximations presented and assessed in McRoberts et al. (2007McRoberts et al. ( , 2011 and Magnussen (2013).
Similar to NN imputation methods, geostatistical methods yield accurate and precise spatial prediction of forest variables. Additionally, geostatistical methods provide a solid theoretical foundation for probability-based uncertainty quantification, see, e.g., Ver Hoef and Temesgen (2013). In such settings, point-referenced observations are in-dexed by spatial locations, e.g., latitude and longitude, and predictive models build on classical kriging methods (Cressie, 1993;Schabenberger and Gotway, 2004;Chilès and Delfiner, 2012). Additional flexibility for specifying and estimating regression mean and covariance components comes from recasting these predictive models within a Bayesian inferential framework (Banerjee et al., 2014). Similar extensions and benefits have been demonstrated for data with observations indexed by areal units, particularly in smallarea estimation settings, many of which build on the Fay-Harriot model (Fay andHerriot, 1979), see, e.g., Ver Planck et al. (2017);Ver Planck et al. (2018) and references therein.
Following the study region description and data in Section 2.1, a general model that subsumes various sub-models is developed in Section 2.2 along with implementation details for parameter estimation, prediction, and model selection in Section 2.3. Here too, we offer two different approaches for prediction over the areal blowdown units.
Results presented in Section 3 focus first on model selection then comparison among blowdown predictions, which is followed by discussion in Section 4. Some final remarks and future direction are provided in Section 5.

Study region and model data
The study region was located in the southern region of the Austrian federal state of Carinthia, within the upper Gail valley and near the Dellach forest research and training center, which is jointly operated by the Institute of Forest Growth and the Institute of Forest Engineering of the University of Natural Resources and Life Sciences Vienna ( Figure 1(a)). On October 28, 2018, the storm Adrian formed over the western Mediterranean Sea and achieved wind gust speeds of 130 km/h throughout Carinthia. Within 72 hours, the storm was producing 627 l/m 2 of precipitation at the Plöckenpass meteorological station, located near the study region, and inflicting heavy damage on the region's forest (Zimmermann, 2018). ministrative district (Bezirk), which covers the study region. Using high-resolution aerial images provided by the Carinthian Forest Service, blowdown occurrences in the district were identified and delineated with high accuracy. As illustrated in Figure 1(b) the blowdowns were concentrated in five distinct sub-regions labeled Frohn, Laas, Liesing, Mauthen, and Ploecken. A total of 564 blowdowns were delineated, totalling 212.3 ha of affected forest. Table 1 provides the number of blowdown occurrences across the sub-regions and affected area characteristics. Forest inventory data were not available for the study region; therefore, a field measurement campaign was initiated in May 2020 (post-Adrian) to collect growing stock timber volume measurements suitable for estimating volume loss due to blowdown. A total of n=62 sample plots were installed in unaffected forest adjacent to blowdowns ( Figure 1(b)). Plot locations were chosen to characterize forest similar in structure and composition to the pre-blowdown forest. No plots were located within blowdowns. Plot measurements were conducted using the terrestrial laser scanning (TLS) GeoSLAM ZEB HORIZON system (GeoSLAM Ltd., Nottingham, UK). Position and diameter at breast height (DBH) for the approximately 5586 measurement trees were derived from 3D point clouds on a 20 m radius plot using fully automated routines demonstrated in Gollob et al. (2019Gollob et al. ( , 2020 and in Ritter et al. (2017Ritter et al. ( , 2020. Tree height was estimated using a Näslund function formulated as a mixed-effects model with plot-level random effects. Stem volume was calculated using a traditional stem-form function (Pollanschütz, 1965). For each plot, growing stock timber volume was expressed as m 3 /ha (i.e., computed as the sum of tree volume scaled by the 7.958 fixed-area plot tree expansion factor). Table 1 provides a summary of growing stock timber volume by sub-region, the values of which serve as the response variable observations in subsequent regression models. Due to the relatively small number of sample plots within any one of the five spatially disjoint sub-regions ( Figure 1), we aimed to pool the 62 sample plot measurements. An ideal pooled model would allow for intra-and inter-location specific relationships between the response and predictor variables. Such spatially varying relationships are particularly attractive because they accommodate potential impact of unobserved (and for the most part unobservable) spatially-explicit mediating factors such as disturbance history, genetics, and local growth environments. When cast within a regression framework, a spatially varying coefficients (SVC) model uses a smoothly-varying spatial process to pool information and avoid overfitting (Gelfand et al., 2003;Finley, 2011).
Given prediction is our primary goal, a preferred model would also estimate residual spatial correlation (i.e., spatial dependence among measurements not explained by the predictor variables) and use it to improve prediction performance. We anticipated the residual spatial correlation would decrease as distance between measurements increased.
Again, following Gelfand et al. (2003) and broader geostatistical literature cited in Section 1, residual spatial dependence is effectively estimated using a spatial process.
To accommodate the anticipated data features, assess varying levels of model complexity, and deliver statistically valid probabilistic uncertainty quantification we model response y(s) at generic spatial location s as where x j (s), for each j = 1, . . . , p, is the known value of a predictor variable at location s, β j is the regression coefficient corresponding to x j (s), β 0 is an intercept, and (s) follows a normal distribution with mean zero and variance τ 2 . Here, τ 2 is viewed as the measurement error variance. The quantities w 0 (s) and w j (s) are spatial random effects corresponding to the intercept and predictor variables, respectively, thereby yielding a spatially varying regression model. To allow for varying levels of model complexity, the δs in Equation (1) are binary indicator variables used to turn on and off the spatial random effects (i.e., a value of 1 means the given spatial random effect is included in the model and 0 otherwise). When δ = 1 the associated space-varying regression coefficient can be seen as having a global effect, i.e., β 0 or β j s, with local adjustments, i.e., w 0 (s) or w j (s)s.
Over n locations, a given spatial random effect w = (w(s 1 ), w(s 2 ), w(s 3 ), . . . , w(s n )) follows a multivariate normal distribution with a length n zero mean vector and n × n covariance matrix Σ with (i, j) th element given by C(s i , s j ; θ). Clearly, for any two generic locations s and s * locations within the study region the function used for C(s, s * ; θ) must result in a symmetric and positive definite matrix Σ. Such functions are known as positive definite functions, details of which can be found in Cressie (1993), Chilès and Delfiner (2012), and Banerjee et al. (2014), among others. Here we specify is a positive support correlation function with φ comprising one or more parameters that control the rate of correlation decay and smoothness of the process. The spatial process variance is given by σ 2 , i.e., Var(w(s)) = σ 2 . This covariance function yields a stationary and isotropic process, i.e., a process with a constant variance and a correlation depending only on the Euclidean distance separating locations. The Mátern correlation function is a flexible class of correlation functions with desirable theoretical properties (Stein, 1999) and is given by where ||s−s * || is the Euclidean distance between s and s * , φ = {φ, ν} with φ controlling the rate of correlation decay and ν controlling the process smoothness, Γ is the Gamma function, and K ν is a modified Bessel function of the third kind with order ν. While it is theoretically ideal to estimate both φ and ν, it is often useful from a computational standpoint to fix ν and estimate only φ. For our current analysis, such a concession is reasonable given there is likely little information gain in estimating both parameters.
In the subsequent analysis, we consider the following candidate models defined using the general model (1): 1. Non-spatial-all δ are set to zero. This is simply a multiple regression model.
2. Space-varying intercept (SVI)-δ 0 equals 1 and all other δs are set to zero. This model estimates a spatial process and associated parameters θ 0 for the intercept, but the impact of the covariates is assumed to be the same across the study region.
3. Space-varying coefficients (SVC)-all δs are set to 1. This model allows all regression coefficients to vary spatially over the study region. Each spatial process has its own parameters θ 0 and θ j for j = (1, 2, . . . , p).

Implementation and analysis
To facilitate uncertainty quantification for model parameters and subsequent prediction, the Non-spatial, SVI, and SVC candidate models defined in Section 2.2 were fit within a Bayesian inferential framework, see, e.g., Gelman et al. (2013)  Data and annotated R code needed to fully reproduce the analysis and results will be provided as Supplementary Material upon publication or if requested by reviewers.

Prediction
Our interest is in predicting the 2020 growing stock timber volume response variablẽ , and τ 2(l) are the l-th joint samples from the parameters' posterior distribution, X is the n 0 × p matrix of predictors at the n 0 prediction locations, and I is the n 0 × n 0 identity matrix. The multivariate Normal PPDs for the SVI and SVC models are given in Finley and Banerjee (2020b,a). Importantly, the SVI and SVC candidate models use joint composition sampling to acknowledge the spatial correlation among prediction locations. A given candidate model's PPD is evaluated using each of its parameter's M posterior samples, i.e., l = 1, 2, . . . , M , to generate M PPD samples for each of the n 0 prediction locations.
These PPD samples are summarized analogously to the parameters' posterior samples.
Prediction point estimates could include the PPD mean or median, and interval esti-mates can be built off the PPD standard deviation or directly expressed using a set of lower and upper percentiles in the form of a credible interval. For example, a point and dispersion estimate for the growing stock timber volume at a prediction location could be the mean and standard deviation over that location's M PPD samples.

Model selection
Our analysis had two separate model selection steps. First, in an effort to build parsimonious models and reduce possible issues arising from collinearity among the often highly correlated ALS predictor variables, we identified a common set of predictors to use in the three candidate models. Given our focus on prediction, predictor variable selection aimed to minimize prediction error using leave-one-out (LOO) cross-validation among the 62 sample plot measurements. Second, given the common set of predictors, model fit and prediction performance were used to select the "best" candidate model for subsequent blowdown area prediction. Here, again, candidate model prediction performance was assessed using LOO cross-validation among the 62 sample plot measurements.
For the first model selection step, the Non-spatial model was used to select the set of ALS predictor variables that minimized LOO mean squared prediction error (MSPE).
As described in Section 2.1, candidate variables were summaries of the ALS point cloud height distribution and included its mean, median, minimum, maximum, and standard deviation. A backward variable selection algorithm, implemented using the leaps (Lumley, 2020) and caret (Kuhn, 2020) R packages, identified the predictor set that minimized LOO MSPE.
The second model selection step used model fit and LOO prediction measures to identify the "best" model for subsequent blowdown area prediction. The deviance information criterion (DIC) (Spiegelhalter et al., 2002) and widely applicable information The three candidate models were also assessed based on LOO cross-validation prediction MSPE and continuous rank probability score (CRPS) (Gneiting and Raftery, 2007

Blowdown prediction
Given the sample plot dataset and posited model identified using methods in Section 2.3.2, two approaches were used to predict growing stock timber volume for blowdown areas. We refer to the approaches as areal and block prediction.
The areal approach views each blowdown as a single prediction unit indexed by its polygon centroid and predictor variables computed as an average of the 1 m × 1 m resolution ALS values over its extent. Figure 2(  stands within the population held at least two plots and the study goal was to improve stand-level point and interval estimates through a smoothing conditional autoregressive model (CAR) model. In our current setting, the sample data comprise response variable measurements collected on fixed-area plots with clearly defined spatial support over which the ALS predictor variables were also computed-the response and predictor variable measurements are spatially aligned at the plot-level. Also, no sample plots fall within the blowdown prediction units; hence, our inferential goal is squarely on out-of-sample prediction. Despite these differences we pursue the areal prediction and clearly acknowledge the change-of-support issue between the discrete fixed-area plot data used to fit the model and different spatial support of the prediction units (see, e.g., Schabenberger and Gotway (2004) for a thorough description of spatial change-ofsupport problems). While not appropriate from a statistical standpoint, we do see this approach used in practice and therefore include it here for comparison with the block approach that mitigates change-of-support issues by better aligning the spatial support of the measurement and prediction units.
For the areal approach, the l-th PPD sample of total growing stock volume (m 3 ) for a given blowdown isỹ (l) A = Aỹ(s) (l) , where A is the blowdown's area (ha) andỹ(s) (l) is the growing stock volume (m 3 /ha) predicted at the blowdown's centroid s.
The block approach partitions each blowdown into grid cells with the same area as the fixed-area sample plots (i.e., 0.126 ha). Each cell is indexed using its centroid, and predictor variables are computed as an average of the 1 m × 1 m resolution ALS values over its extent; see illustration in Figure 2(b). Akin to block kriging (Wackernagel, 2003), the response variable prediction for a given blowdown is an aggregate of multiple point predictions within the blowdown extent. More specifically, given a blowdown divided into n 0 cells and corresponding vector holding the l-th joint PPD sampleỹ (l) (m 3 /ha), the PPD sample of total growing stock volume (m 3 ) for the blowdown isỹ , where a(s i ) is the area of the i-th cell that falls within the blowdown (a i is at most 0.126 ha when the cell's extent is completely within the blowdown).
Composition sampling was again used to generate M samples fromỹ A 's andỹ B 's PPD for each blowdown. These PPD samples were summarized analogously to the parameters' posterior samples to yield prediction point, dispersion, and interval estimates for the 564 blowdowns. In addition to blowdown specific PPD summaries, extra composition sampling was used to estimate the total volume PPD by sub-region and region.

Candidate models
Following Section 2.3.2, the model that yielded minimum LOO cross-validation MSPE included only the ALS point cloud height distribution mean predictor variable. We refer to this predictor variable as mean canopy height and set it as x 1 in Equation (1) for all candidate models.
Parameter estimates for candidate models are given in Table 2. As expected, the β 1 estimates indicate a strong positive relationship between 2012 mean canopy height and 2020 growing stock volume. The addition of spatial random effects decreased the non-spatial residual variance τ 2 . The reapportionment of τ 2 to σ 2 0 and non-negligible σ 2 1 suggests a substantial portion of variance-not explained by the ALS predictor-had a spatial structure and the ALS predictor had space-varying impact.
The SVI model spatial decay parameter estimates suggest a fairly localized spatial structure. Given the exponential spatial correlation function and km map projection units, the distance d 0 at which the spatial correlation drops to 0.05 (an arbitrary, but commonly used, cutoff) is estimated by solving 0.05 = exp(−φd 0 ) for d 0 providing d 0 = − log(0.05)/φ. The distance d 0 is commonly referred to as the effective spatial range. Using the SVI model φ 0 estimates, the corresponding spatial range is 0.26 (0.10,

1.58) km.
Compared with the SVI model, the SVC model further reduced the non-spatial residual variance by taking into account the space-varying impact of x 1 . Relative to the effective spatial range of the intercept process, the process on x 1 had a long spatial range, i.e., 29.96 (0.88, 99.86) km. That is, the spatially varying slope coefficient for x 1 represented sub-region differences in the relationship between growing stock volume and mean canopy height that were probably caused by unmeasured species, genetic, or environmental factors.

Blowdown prediction
The SVC model provided the "best" fit and LOO predictive performance (Table 3) and therefore served as the prediction model for the blowdowns. Following methods in Section 2.3.3, areal and block growing stock volume PPD mean and standard deviation were computed for each blowdown, the results of which are plotted in Figure 3. Figure 3(a) shows negligible difference between areal and block PPD means. However, as shown in Figure 3(b), compared with the block approach, the areal prediction resulted in a consistently larger PPD standard deviation. Additionally, supplemental Figure S1 shows the areal PPD coefficient of variation (CV) is generally larger than the block PPD CV.  Table 4 provides sub-region growing stock volume loss totals and corresponding 95 % confidence intervals. Although total blowdown area in Frohn was larger than other subregions, its per unit area growing stock loss 510.68 m 3 /ha (32,086.1 m 3 /62.83 ha) was less than that in Laas (637.26 m 3 /ha), Liesing (818.21 m 3 /ha), Mauthen (784.76 m 3 /ha), and Ploecken (638.50 m 3 /ha). This disparity between Frohn and the other sub-regions is because the blowdowns in the Frohn sub-region were concentrated in relatively unproductive forests close to the alpine treeline zone. Maps of PPD summaries per blowdown were prepared for each sub-region. For example, Figure 4 shows blowdowns in the Frohn sub-region colored by growing stock volume PPD estimate mean, standard deviation, and coefficient of variation. Similar maps for the other sub-regions are provided in the Supplemental Material. Figure 5 shows the distribution of PPD CV by blowdown area.  Given the mean canopy height predictor, the SVC model showed consistently better fit and predictive performance compared to the Non-spatial and SVI models (Table 3), and was therefore selected as the model to generate areal and block prediction for the blowdowns.
In regression, and most other modeling contexts, it is assumed the predictor variables used to estimate model parameters are the same as those used in subsequent prediction. This assumption does not hold for the areal prediction approach described in Section 2.3.3, because the ALS variables computed over the 62 0.126 ha sample plots could have a different distribution than the ALS variables computed over the variable area blowdowns. Although it is reasonable to expect the mean of a given ALS variable to be similar between the sample plot and blowdown distributions, the dispersion of the two distributions could be quite different. That is, a change-of-support was immanent between the sample plots (model data) and the prediction units' areal extent. This fact alone should precluded areal prediction application; however, our analysis results further underscore the approach's shortcomings.
In the areal approach each blowdown was considered as a single prediction point, and its predicted volume per unit area was scaled by the blowdown's area to arrive at the blowdown's total volume. Such an approach might seem intuitive, but is problematic in the current setting. The issue is akin to how variance scales with forest inventory plot size. It is often the case that a forest variable's variance is inversely related to the area over which the variable is measured (Freese, 1962). Consequently, a timber volume variable shows more variability when measured on smaller plots than on larger plots, especially when applied to structurally complex forests. This is because the larger plots average over local scale structural variability. When using a single prediction at the blowdown's centroid, the areal prediction variance does not scale with blowdown area. In contrast, the block approach more accurately scales prediction variance with blowdown area as reflected in Figure 3(b). This is because its PPD is the result of an average over possibly multiple spatially correlated prediction units (grid cells) within the blowdown's boundary-the number of prediction units contributing to this average scales with blowdown area.
Finally, the PPDs rely on the spatial correlation between observed and prediction locations, e.g., using Equation 2, to appropriately weight the contribution of observed data for prediction. For the areal approach, this correlation is computed between observed locations and a blowdown's polygon centroid, which is the average location relative to the polygon's vertices. If the polygon's shape is highly irregular, with a large boundary length to area ratio, then its centroid might not represent well the distance between observed locations and the polygon's extent. As an extreme example, the centroid of a polygon in the shape of a letter "C" or "L" will fall outside the polygon boundary.
Again, such issues are circumvented using the block prediction approach.
A particularly useful quality of the Bayesian inferential paradigm is PPDs can be generated for any arbitrary set of prediction units from a single grid cell to sets of blowdowns that comprise a sub-region or entire study region. Inference proceeds from desired joint PPD samples. These PPD samples can be summarized using measures of central tendency, dispersion, intervals, and transformations, with results presented as maps or tables, e.g., Figure 4 and Table 4.

Conclusions
The study goal to estimate growing stock timber volume loss due to storm Adrian in the Austrian upper Gail valley in Carinthia, was met using ALS and TLS measurements coupled through a flexible spatial regression model cast within a Bayesian inferential framework. Limited data availability and its configuration in space and time presented several inferential constraints on how statistically robust estimates could be pursued.
Our proposed regression model was designed to leverage information from a small set of plot samples and aerial ALS as well as spatial autocorrelation among forest measurements and nonstationary relationships between response and predictor variables.
Three candidate models of varying complexity were assessed using model fit and outof-sample prediction. Performance metrics supported the SVC model which was then used to make areal and block predictions over the blowdowns. Our results showed that in contrast to the areal approach, the block approach mitigates issues with change-ofsupport by matching the prediction unit to the sample plot extent. Using this finer spatial scale prediction unit and a joint prediction algorithm that acknowledges spatial correlation among prediction units, the total growing stock volume PPD captures the correlation among prediction units within a given blowdown and scale appropriately with blowdown area. The block prediction approach facilitated statistically sound inference at various spatial scales, i.e., blowdown, sub-region, and region levels.
The proposed methodology and annotated code that yields fully reproducible results can be adapted to deliver damage assessment for other forest disturbance events in future periods and different geographic regions. While additional sample plot data would improve estimates in our current study, we were able to demonstrate a fairly high level of accuracy and precision is achievable using a limited sample size. This small sample, i.e., n=62 plots, was collected using a TLS which, compared with traditional individual tree measurements, allowed for time and effort efficient data collection. Based on our experience from this and other efforts, a field crew can collect ∼20-25 plots per day, even in difficult alpine terrain.
Our study illustrated a methodology to efficiently deliver information required for strategic salvage harvesting following storm and other disturbances. Future work focuses on augmenting the SVC model to incorporate large spatial datasets using Nearest Neighbor Gaussian Processes , automate remote sensing predictor SVC block prediction total growing stock coefficient of variation (%) SVC areal prediction total growing stock coefficient of variation (%) Figure S1: Summaries of each blowdown's total growing stock volume (m 3 ) posterior predictive distribution computed using areal and block prediction approach. Points represent blowdowns, colored by area, and broken down by sub-region with a one-to-one line.