Next Article in Journal
Genomic Prediction of Tree Height, Wood Stiffness, and Male Flower Quantity Traits across Two Generations in Selected Individuals of Cryptomeria japonica D. Don (Japanese Cedar)
Next Article in Special Issue
Tree Biomass Modeling Based on the Exploration of Regression and Artificial Neural Networks Approaches
Previous Article in Journal
LIDAR-Based Forest Biomass Remote Sensing: A Review of Metrics, Methods, and Assessment Criteria for the Selection of Allometric Equations
Previous Article in Special Issue
Formulating Equations for Estimating Forest Stand Carbon Stock for Various Tree Species Groups in Northern Thailand
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Forests
Volume 14
Issue 10
10.3390/f14102096
forests-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Thinning Combined with Prescribed Burn Created Spatially Heterogeneous Overstory Structures in Contemporary Dry Forests: A Comparison Using LiDAR (2016) and Field Inventory (1934) Data

by
Sushil Nepal
1,*,
Bianca N. I. Eskelson
1,*,
Martin W. Ritchie
2 and
Sarah E. Gergel
3
1
Department of Forest Resources Management, The University of British Columbia, 2424 Main Mall, Vancouver, BC V6T 1Z4, Canada
2
USDA, Forest Service Pacific Southwest Research Station, 3644 Avtech Parkway, Redding, CA 96002, USA
3
Department of Forest and Conservation Sciences, The University of British Columbia, 2424 Main Mall, Vancouver, BC V6T 1Z4, Canada
*
Authors to whom correspondence should be addressed.
Forests 2023, 14(10), 2096; https://doi.org/10.3390/f14102096
Submission received: 30 August 2023 / Revised: 27 September 2023 / Accepted: 7 October 2023 / Published: 19 October 2023
(This article belongs to the Special Issue Modeling Aboveground Forest Biomass: New Developments)
Browse Figures
Review Reports Versions Notes

Abstract

:
Restoring current ponderosa pine (Pinus ponderosa Dougl. Ex P. and C. Laws)-dominated forests (also known as “dry forests”) to spatially resilient stand structures requires an adequate understanding of the overstory spatial variation of forests least impacted by Euro-American settlers (also known as “reference conditions”) and how much contemporary forests (2016) deviate from reference conditions. Because of increased tree density, dry forests are more spatially homogeneous in contemporary conditions compared to reference conditions, forests minimally impacted by Euro-American settlers. Little information is available that can be used by managers to accurately depict the spatial variation of reference conditions and the differences between reference and contemporary conditions. Especially, forest managers need this information as they are continuously designing management treatments to promote contemporary dry forest resiliency against fire, disease, and insects. To fill this knowledge gap, our study utilized field inventory data from reference conditions (1934) along with light detection and ranging and ground-truthing data from contemporary conditions (2016) associated with various research units of Blacks Mountain Experimental Forest, California, USA. Our results showed that in reference conditions, above-ground biomass—a component of overstory stand structure—was more spatially heterogeneous compared to contemporary forests. Based on semivariogram analyses, the 1934 conditions exhibited spatial variation at a spatial scale < 50 m and showed spatial autocorrelation at shorter ranges (150–200 m) compared to those observed in contemporary conditions (>250 m). In contemporary conditions, prescribed burn with high structural diversity treatment enhanced spatial heterogeneity as indicated by a greater number of peaks in the correlograms compared to the low structural diversity treatment. High structural diversity treatment units exhibited small patches of above-ground biomass at shorter ranges (~120 to 440 m) compared to the low structural diversity treatment units (~165 to 599 m). Understanding how spatial variation in contemporary conditions deviates from reference conditions and identifying specific management treatments that can be used to restore spatial variation observed in reference conditions will help managers to promote spatial variation in stand structure that has been resilient to wildfire, insects, and disease.

1. Introduction

Overstory stand structure of ponderosa pine-dominated forests—hereafter “dry forests” [1]—has undergone substantive changes since the Euro-American settlement of California, USA. Such changes have included shifts in species composition, an increase in the density of small trees [2], and a decline in the density of large trees [3,4,5]. An increase in densities of small trees and shifts in species composition are attributed to a wide range of factors, including fire suppression, logging, grazing, and climate change [2,3].
In California, tree-ring reconstructions indicate that spatial variation in the overstory of dry forests impacted by indigenous land use but prior to harvesting activities following Euro-American settlement (hereafter termed “reference conditions”) consisted of a mosaic of individual trees, tree groups, and gaps at small scales (<0.4 ha or 40–70 m). Such spatial variation in overstory stand structure resulted in open forests with sparse, large individuals of fire-resistant species [6,7]. Spatial variations in stand structure are the result of interactions among prior processes such as fire, regeneration, competition, and mortality [8,9]. However, spatial variations in stand structure such as tree density, basal-area, and above-ground biomass (AGB) have substantially changed following Euro-American settlement [10,11].
Within the dry forests of California, following Euro-American settlement, spatial variation in overstory stand structure has shifted from fire resilient to fire prone forests due to increased tree density of shade-tolerant species such as white fir (Abies concolor (Grod. & Glend)) and incense-cedar (Calocedrus decurrens (Torr.)) [10,11,12]. Frequent fires with return intervals of about 25 years are among the most important drivers of spatial variation in dry forests [13,14]. However, it is unclear how overstory stand structures may have been altered spatially due to the exceptionally long fire-free period after Euro-American settlement [3,15]. Therefore, quantitative descriptions of spatial changes in overstory stand structures are crucial for land managers as they can be used to assess the potential of wildfire due to increased fuel loads in the forest overstory [16,17].
The spatial variations in overstory stand structure over time have rarely been studied, which is potentially a result of both a general lack of spatially-explicit data collected prior to active management [3] and the fact that such data collection is extremely labor intensive and expensive. Hence, it is practically impossible to obtain census data over large areas [18]. Remote sensing technology, especially light detection and ranging (LiDAR), provides the capacity to obtain spatially explicit data over large areas in a timely and cost-effective manner [18]. Although LiDAR data are being extensively used in enhanced forest inventory, ground data are still required [19]. LiDAR-derived variables such as height metrics can be utilized with ground-based observations such as basal area (m2/ha), volume (m3/ha), and AGB (Mg/ha) to predict forest overstory stand structure at different times [20]. Thus, LiDAR data complement and can be used in conjunction with ground-based inventories to identify spatial changes in the forest overstory structure [21].
Furthermore, many studies examining the spatial variation of forest structure in California, USA, have been restricted to particular elevation ranges [2,22], management units [23,24], certain functional types, and small study areas [3,25,26]. Numerous studies have provided a general description of dry forests prior to active management following Euro-American settlement [27,28]. However, a description of the spatial variation of overstory stand structure over time is generally lacking. Many studies have investigated spatial variation in tree density and tree size (e.g., [9,29]), and at different scales based on averages (e.g., [30,31]). However, they have failed to capture the spatial variation of overstory over a larger landscape utilizing different metrics such as AGB as a measure of overstory stand structure. AGB is closely related to forest productivity and can be more accurately predicted and modeled from ground-measured biomass and LiDAR height metrics compared to trees per hectare or basal area per hectare [20]. Metrics such as trees per hectare and basal area per hectare are sensitive to the inability of LiDAR to capture small and understory trees very accurately [31].
This study was conducted in the dry forests of Blacks Mountain Experimental Forest (BMEF) in northern California, USA. BMEF remained untouched by selective logging operations from the late 1800s through 1940s [32], but may have been impacted by indigenous land use for hunting mule deer [33]. Within BMEF, our study focused on the section defined as the Blacks Mountain Ecological Research Project (BMERP), initiated in 1991. Within BMERP, two structural diversity thinning treatments with or without prescribed burn were implemented [32]. A low structural diversity treatment (LOD) was designed to produce a single canopy layer by removing dominant and large-sized trees [32]. A high structural diversity treatment (HID) was designed to leave all the dominant trees, abundant snags, and multiple canopy and forest openings [32]. To quantify the spatial variation of AGB, we utilized three types of data collected for the BMERP project area: ground inventory data from 1934 as reference data, as well as contemporary 2015 LiDAR data, and 2016 ground verification data. The specific objectives of our study were to:
  • compare the spatial variation of overstory AGB between contemporary forests and reference conditions using 2015 LiDAR data and 1934 census data;
  • compare the spatial variation of overstory AGB among various structural diversity treatments in contemporary conditions.
The variogram analyses implemented in this study, provide information on the inherent patch size and spatial heterogeneity for overstory AGB in both contemporary forests and reference conditions. Understanding the spatial differences in overstory AGB at two points in time across a landscape has various applications in forest management policies such as fuel classification, fire spread prediction, and post-disturbance vegetation changes [34,35,36]. Such information can improve our understanding of how trees and forests respond and will continue to respond towards changes in disturbance regimes [37] and help silviculturists design restoration treatments that move forests towards more resilient conditions similar to reference conditions.

2. Materials and Methods

2.1. Study Area

The BMEF, managed by the USDA Forest Service Pacific Southwest Research Station, is located in northeastern California (Figure 1) (40°40′ N, 10 121°10′ W), northeast of Lassen Volcanic National Park. The elevation of BMEF ranges between 1700 and 2100 m [32]. Slopes rarely exceed 30 percent [32]. Aspects are primarily west- and south-facing. At lower elevations, stands are dominated by ponderosa pine (Pinus ponderosa Dougl. Ex P. and C. Laws) with occasional occurrence of some Jeffrey pine (Pinus jeffreyi (Grev. And Balf.) [32]. At higher elevations, white fir (Abies concolor (Gord. And Glend.) Lindl.) and incense-cedar (Calocedrus decurrens (Torr.) Florin) dominate the stands. Classified as an interior ponderosa pine forest type [38], the 4358 ha forest has a wide range of stand conditions as a result of past research and management activities, as well as disturbance events [39].
For this study, we focused on ten of the twelve research units in the BMERP and four research natural areas (RNA, Figure 1C). BMERP was initiated as an interdisciplinary large-scale, long-term ecological research project at BMEF in 1991 [32]. The goals of BMERP were to: (a) understand the effects of forest structural complexity on the health and vigor of ponderosa pine ecosystems, (b) quantify the ecosystem’s resilience to natural and human-caused disturbances, and (c) determine how these ecosystems can be managed for sustained resource values [32]. The forests in the ten BMERP research units were subjected to two different types of treatments (Figure 1C, for details about treatments see Appendix A, Table A1). The first treatment consisted of three stand structures: low structural diversity (LOD), high structural diversity (HID), and research natural areas (RNAs) (Figure 1C; [32]). LOD and HID treatments had been randomly assigned to the ten BMERP research units ranging in size from 77 to 144 ha (Figure 1C). Each research unit was then split in half with one randomly assigned half receiving prescribed burn treatments (hereafter “burned”), whereas the other half did not receive the prescribed burn treatment (hereafter “unburned”) (Figure 1C). Collectively, the twelve units consisted of a total of twenty-four stands in our study: four RNAs, ten LOD stands, and ten HID stands (Figure 1C and Table A1). LOD stands were thinned to a uniformly spaced density of ~40 trees ha−1 retaining a single canopy layer of intermediate trees with height ranging from 12 to 30 m and a crown ratio > 50% [32]. In contrast, thinning in HID stands was conducted to retain all canopy layers that represented an overstory stand structure of a forest with multiple age classes and varying crown structures [32]. All large old trees were maintained with one smaller tree retained within the larger tree’s crown circumference [32]. As a part of the prescription, within the HID units, caches of high-density and small-diameter conifers were left [32,38].
The four RNAs (~40 ha each; Figure 1C) were set aside to serve as unmanaged, qualitative controls representative of the interior ponderosa pine type [32]. These RNAs have never received mechanical treatments, but fire exclusion has greatly increased their understory tree densities [32]. Two of the four RNA stands (RNA-B and RNA-C) received one application of prescribed burn in the late 1990s [32].

2.2. Data

2.2.1. 1934 Survey Data

During the fall of 1934, BMEF was divided into a total of 4074 rectangular plots of 1 ha (hereafter referred to as “1 ha plots”) in size to conduct a census survey ([40], Figure 1C). Within each 1 ha plot, all live overstory trees > 8.9 cm (3.5 inches) were tallied by species group and diameter class [40]. Diameter classes were 5.08 cm (2 inches) wide and labeled by even inches (i.e., 4, 6, 8, 10, etc.). Tree species were recorded as pine (Jeffrey or ponderosa pine), white fir, and incense-cedar [40].

2.2.2. 2015 LiDAR Data Acquisition and Processing

An airborne Light Detection and Ranging (LiDAR) dataset was acquired during the summer of 2015 using a Leica ALS50 PHASE II laser system (near-infrared) discrete return sensor mounted on a fixed-wing aircraft [41]. The aircraft was flown at an altitude of 900 m using an opposing-flight line with a side lap of at least 50% [41]. The scanning angle of the sensor was ±14° with an average return of 6.9 points per m2 and a standard deviation of 5.9 points per m2 (see [41]). LiDAR point cloud was processed and LiDAR metrics were calculated using the ‘lidR’ package in R ([42], see Appendix B for a full description of the process).

2.2.3. 2016 Field Data

The ten research units and four RNAs each have a permanent 100 m lattice grid (Figure 1D; [32]). The grid serves as the center points for most of the plot-level research conducted in BMEF [32]. In the summer of 2016, at every other grid point in all diagonal directions (282 m spacing, Figure 1D), 16 m radius plots (804 m2, hereafter referred to as “circular plots”) were measured within each treatment unit [43]. Standing live and dead trees ≥ 9 cm diameter at breast height (DBH) were stem mapped from the plot center and measured for total height and DBH [43]. These ground-truthing data were collected in a total of 154 circular plots: 65 LOD, 69 HID, and 20 RNA plots (Figure 1C; blue dots).

2.2.4. Overstory Above-Ground Biomass for 1935 and 2016

Using the 1934 survey data, I calculated biomass (Mg) for foliage, branch, and bole of individual trees > 9 cm DBH using species-dependent equations that were developed locally at BMEF [44]. Then, the above-ground biomass (AGB, Mg) for each tree was calculated by summing the individual biomass from foliage, bole, and branch. The height values used in the equations were estimated from local height–diameter equations developed by Dolph et al. [45]. For white fir and incense-cedar overstory AGB, we used equations and parameters suggested by Jenkins et al. [46]. The total plot overstory AGB was calculated by adding the overstory AGB of each species. We then converted plot AGB to AGB per hectare (Mg/ha) based on the given plot sizes.
For the 2016 field data, we calculated the AGB for individual trees using the same species equation parameters that were used for the 1934 biomass calculations. The total overstory AGB for each circular plot was calculated by taking the sum of the AGB calculated for individual trees of each species and converted into per ha values.

2.2.5. Biomass Model to Link LiDAR Metrics to 1 ha Plots for 2016

The distribution of the response variable—overstory above-ground biomass (AGB, Mg/ha)—was skewed to the right with values > 0 (Figure A1). Therefore, we fit generalized linear mixed effect models (glm) with a ‘gamma’ distribution and a ‘log link’ function to ensure that the expected value of AGB is always positive [47]. Modeling was performed in R using the ‘lme4’ package ([48]; see Appendix C for full details on the model building process and sumary results).
We divided the study area within BMERP into cells of size 32 × 32 m to match the resolution of the ground-truthing plots (Figure A4). LiDAR metrics were extracted at the cell level following the method described in Mauro et al. [43]. There was a total of 12,647 cells. The selected biomass model was then applied to the cells to predict overstory AGB. To calculate the overstory AGB for 1 ha plots in 2016, first, the predicted values of overstory AGB at cell level were created as a map (Figure A4). Then, the map was overlaid with the 1 ha plot map with the BMERP units. The intersect function in ArcMap 10.4.1 was used to extract the cells that fell within 1 ha plots [49]. If a cell straddled multiple 1 ha plots, overstory AGB was weighted based on the area of the cell, and area weighted values of cells that fell within the plots were summed to calculate the overstory AGB (Mg/ha) for 1 ha plots.

2.3. Data Analysis

We computed semivariograms at the scale of each BMERP unit and Moran’s I correlograms as a measure of spatial variation to evaluate our research objectives.

2.3.1. Comparing Semivariogram Components and Moran’s Correlogram between 1934 and 2016

Semivariograms—To meet our first research objective, we used semivariogram models following the methods of Rossi et al. [50] with the assumption of stationarity, i.e., that the variance in above-ground biomass (AGB) is the function of separation distance only [51]. The semivariance was calculated for a pair of observations of overstory AGB as a function of the separation distance (hereafter “lag distance”) between the sampled locations [51]. Using the ‘gstat’ spatial package in R version 4.1.1 [52], semivariogram models were built for all fourteen units using a 1 ha plot-level overstory AGB for 1934 and 2016. For the 2016 study units, we did not differentiate between the burned and unburned halves; hence, variograms represent the pooled variogram for each unit. Following the method of Fry and Stephens [53], we fit three semivariogram models—exponential, Gaussian, and spherical—for each unit in 1934 and 2016. We also checked the assumption of stationarity in our exploratory analysis using directional variograms and found that ranges were not different; hence, semivariance was only distance dependent and not direction dependent [53]. Furthermore, to allow for comparisons at a common scale, all variograms were standardized by dividing the semivariance by the overall sample variance [50]. For both 1934 and 2016, the best semivariogram models were selected based on the minimized root mean squared error (RMSE, [53]). Following the method suggested by Fry and Stephens [53], we used the selected semivariograms from both points in time and compared the values of range, nugget, and sill for all fourteen units using dot plots. The range describes the distance up to which overstory AGB values exhibited spatial autocorrelation and provides information on the inherent patch size and spatial heterogeneity for overstory AGB [54]. The sill values were compared to understand whether the amount of spatially dependent variance within a given range in overstory AGB differed between 1934 and 2016 [53,55]. Nugget values were compared to understand if there was spatial variation at a scale smaller than 50 m (hereafter referred as “fine-scale spatial variation”), which was the shortest distance between adjacent 1 ha plots in overstory AGB [55].
Moran’s I correlogram—For a more local measure of spatial autocorrelation, we followed the method of Jaquette et al. [56] and calculated Moran’s I over the range of 20 lags at an interval of 50 m between lags, which was the shortest possible distance between 1 ha plot centers. From our exploratory analysis, we found that Moran’s I could be calculated up to a distance of 1000 m for all units except RNA-D because the number of observations was <5% of total observations for the given unit and the spatial variation could not be interpreted easily with so few observations [57]. RNA-D was an exception and only allowed 800 m for the Moran’s I calculations. We also tested if Moran’s I at each lag was significantly different from 0 (alpha = 0.05) using Monte Carlo simulations with 1000 permutations [56]. Using the results from the 1000th permutation, all directional Moran’s I correlograms were constructed with lag distance (m) on the x-axis and Moran’s I value on the y-axis [57,58] for all units. We used the Moran’s I correlograms for 1934 and 2016 to compare spatial variation in overstory AGB in terms of the differences in:
(a)
fine-scale spatial variation—We evaluated whether the magnitude of Moran’s I at the first lag was significantly different from 0 (p < 0.05) for each unit in both years. A Moran’s I value at the first lag that is significantly different from 0 indicates a lack of fine-scale spatial variation within study units, i.e., patches of overstory AGB < 50 m do not exist [57]. If Moran’s I was not significantly different from 0, patches of AGB < 50 m exist in the study units.
(b)
periodicity in Moran’s I correlogram—Moran’s I values that are significantly different from 0 at different lag distances result in peaks (positive Moran’s I) and troughs (negative Moran’s I) at different lag intervals creating a periodic correlogram for overstory AGB [57,58]. A greater number of peaks and troughs in the correlogram indicates greater spatial variation, whereas fewer peaks in the correlogram suggest less spatial variation in AGB [57,58]. We visually compared if the number of peaks and troughs combined (collectively referred to as “peaks”) in the correlograms for the units were different between years as an indicator of greater or less spatial variation in AGB [57,58].

2.3.2. Effect of Management Treatments on Spatial Variation of Above-Ground Biomass in the Contemporary Forests

To compare the effect of management treatments on spatial variation of overstory AGB, we used semivariogram models built for each of the burned and unburned halves of the HID and LOD treatment units (Table A1). A total of 20 semivariograms—10 each for the burned and unburned halves of HID and LOD treatment units—were selected based on low RMSE values and used to compare the range, nugget, and sill. We also used Moran’s I correlogram constructed for both burned and unburned halves of HID and LOD treatment units to compare the fine-scale spatial variation and periodicity in spatial variation in overstory AGB as described in Section 2.3.1. From our exploratory analysis, we found that maximum distance for which we could construct Moran’s I correlograms was 1000 m for all the burned and unburned halves of the units.

3. Results

3.1. Spatial Variation between 1934 and 2016

3.1.1. Spatial Autocorrelation in Above-Ground Biomass Exhibited at Larger Ranges in 2016

Within research natural areas (RNAs), Gaussian and spherical models for various units exhibited low RMSE; thus, they were selected for comparison of spatial variation in overstory above-ground biomass (AGB) between 1934 and 2016 (Table 1, Figure A5). The models that did not converge for either year are not presented in the results. Within RNAs, the values of range did not differ substantially between RNA-A and RAN-B for both years (Figure 2 and Figure A5). However, within RNA-C and RNA-D, we observed larger ranges of ~400 m and 250 m, respectively, in 2016 compared to 1934 (~243 m for RNA-C and 167 m for RNA-D, Figure 2 and Figure A5). RNAs exhibited larger nugget values in 1934 compared to extremely small values in 2016 indicating the presence of fine-scale spatial variation in 1934 compared to 2016 (Figure 3 and Table 1). In 2016, the sill values were generally similar to those observed in 1934. We only observed higher sill values in RNA-D in 2016 compared to 1934, without much difference in other RNAs (Table 1).
Within units other than RNAs, spherical and Gaussian models exhibited low RMSE values and were selected (Table 1 and Figure A6). All units except 41 and 44 exhibited short ranges in 1934 compared to longer ranges in 2016 (Figure 2, Figure A6 and Figure A7). There was a very small difference in range values between both years in units 41, 43, 44, and 45 (Figure 2). Nugget values observed in 1934 were larger than those observed in 2016 within 14 units of the study, which indicated the presence of fine-scale spatial variation (Figure 3). The amount of spatial autocorrelation indicated by the sill values was greater in 1934 compared to 2016 in most of the units (Table 1 and Figure A6 and Figure A7).

3.1.2. Moran’s I Correlogram Exhibited Differences in Fine-Scale Spatial Variation and Periodicity between Two Years

Moran’s I correlogram suggested that fine-scale spatial variation was more pronounced in 1934 conditions compared to 2016, as indicated by the Moran’s I values—not significantly different from 0—at the first lag in most of the units (Table 2). For example, in 1934, nine out of fourteen units exhibited a Moran’s I value at the first lag that was not significantly different from 0, but in 2016, ten out of fourteen units had Moran’s I values at the first lag that were significantly different from 0 (Table 2). This was in agreement with all the higher nugget values found in 1934 as compared to low nugget values in 2016 (Figure 3). Thus, only four out of fourteen units exhibited fine-scale spatial variation in 2016 (Table 2).
In general, periodicity did not vary substantially between 2016 and 1934, as indicated by similar numbers of peaks between both years (Table 2, Figure A8). However, the presence of periodicity in the Moran’s I correlograms in both years indicated that overstory AGB occurred in heterogeneous patches (Figure A8). Exceptions were RNA-A and RNA-C, where the spatial variation in overstory AGB was more pronounced in 2016, as indicated by a larger number of peaks in 2016 compared to 1934 (Table 2 and Figure 4). The periodicity in the correlogram occurred at larger lag intervals (>250 m) for units 39, 43, and 48 in 2016, which is an indication of the presence of bigger patches of overstory AGB (Figure A8).

3.2. 2016 Spatial Variation in Treatment Units

3.2.1. HID Burned Halves Exhibited Spatial Autocorrelation at Short Ranges

Variogram models with the lowest RMSE for both burned and unburned halves of HID and LOD treatments were selected for comparing range, nugget, and sill (Table A5). Irrespective of the burned and unburned halves, all the HID treatment units exhibited spatial autocorrelation for overstory AGB at shorter ranges (~120 to 440 m) compared to the LOD treatment units (~165 to 599 m, Figure 5, Figure A9 and Figure A10). All the burned halves of HID showed larger nugget values compared to unburned halves and most of the LOD treatments, indicating the presence of fine-scale spatial variation in overstory AGB within burned halves (Figure 6).
Within LOD treatment units, most of the burned halves exhibited spatial autocorrelation of overstory AGB at a smaller range compared to unburned halves, but the difference was very small (Figure 5). Exceptions were units 43 and 44, where spatial autocorrelation in overstory AGB was exhibited at shorter ranges in unburned halves compared to burned halves (Figure 5). Within HID treatment units, all the unburned halves exhibited spatial autocorrelation in overstory AGB at larger ranges than the burned halves, except for units 38 and 48, where the differences in range between burned and unburned halves was very small (Figure 5 and Figure A10). The presence of fine-scale spatial variation in overstory AGB was indicated by larger nugget values in burned halves compared to unburned halves, except for units 41 and 42, which exhibited no nuggets (Figure 6).

3.2.2. HID Burned Halves Exhibited Greater Spatial Variation Indicated by Periodicity at Shorter Lag Distances

Irrespective of treatments, most of the burned and unburned halves exhibited Moran’s I values at the first lag that were significantly different from 0, which indicated the absence of spatial variation at distances below 50 m (Table 3, Figure A11). We found a greater number of peaks that were significantly different from 0 at shorter lag distances within burned halves of HID treatments than within unburned halves of HID and LOD treatments (Figure A11). For example, the burned half of unit 38, showed eight peaks at the interval of 100–150 m lag distances compared to only two peaks in the unburned half of unit 39, and three peaks in the burned half of unit 39, all of which occurred at larger lag distances between 300 and 400 m (Table 3, Figure A11A–F).
Correlograms for LOD treatments were more uniform with few significant peaks (p = 0.0009) of Moran’s I at longer lag distances compared to HID treatments (Figure A11, LOD vs. HID). However, between burned and unburned halves of LOD treatment units, we observed a greater number of significant (p = 0.0009) negative and positive Moran’s I values at short lag distances for the burned halves (Table 3, Figure A11, LOD). In addition, between unburned halves of HID and LOD, HID unburned halves for units 41, 42, 47, and 48 exhibited a greater number of peaks in Moran’s I values significantly different from 0 at shorter lag distances (Figure A11). Hence, both HID burned and unburned halves exhibited more spatial variation in AGB compared to LOD treatments (Figure A11, LOD vs. HID).

4. Discussion

Results showed that ranges, fine-scale variation, and periodicity in spatial variation of overstory above-ground biomass (AGB) were different among units in BMEF based on management and the point in time. Forests of 1934 were heterogeneous without marked human-induced management. In contrast, by 2016 these forests had been subjected to some form of thinning, burning, or a combination of both at various points in time [32,59]. Hence, spatial variation in overstory AGB for 2016 was likely influenced by these treatments. Our results showed that spatial variation in 2016 was more pronounced in prescribed burns.

4.1. The 1934 Forest Was Spatially Heterogeneous Particularly at Fine Scales

Most forest units had greater spatial variation in overstory AGB biomass in 1934 compared to 2016, as indicated by the short value of ranges (<200 m) from semivariogram, periodicity of Moran’s I correlogram and Moran’s I value at first lag not significantly different from 0 in RNAs [57,58]. The spatial variation in overstory AGB in BMEF is consistent with several other studies describing small and frequent patch size distributions (ranging from 60 to 150 m) of the overstory in old growth ponderosa pine forests in northern California [9,60] and Washington [1].
Further, our results identifying fine-scale variation in most of the units were consistent with other studies, which attributed the presence of fine-scale variation to frequent fire [61] and fine-scale variation in growing environmental conditions [62]. Exceptions occurred in units such as RNA-A and units 45 and 48, which lacked fine-scale (<50 m) spatial variation. The absence of fine-scale spatial variation is likely due to the absence of frequent fire after around 1880s at BMEF (unpublished data from Carl Skinner). In the absence of frequent fire, stand density possibly increased due to the infilling of gaps with shade-tolerant species such as incense-cedar and white fir [9] in units such as RNA-A and units 45 and 48, which occurred at higher elevations with abundant shade-tolerant species.

4.2. Spatial Autocorrelation at Long Ranges and Less Spatial Varaitons in 2016 except When Burned

The long ranges (>250 m) of the variograms suggested that overstory AGB in most units in 2016 consisted of small patch size and spatial variations, possibly due to increases in tree densities within 1934 forest gaps [3,63]. With the long-term (century-long) absence of frequent fire, heterogeneity at small scales (<50 m) could have decreased via a loss of fine-scale patches along with expansion of existing patches into gaps [3]. The loss of fine-scale patchiness in BMEF was evidenced by the decrease in nugget values and significant Moran’s I at first lag in 2016 compared to 1934, consistent with [61]. On the other hand, no difference in periodicity between 1934 and 2016 within prescribed burn RNA units indicates that at localized scales, prescribed burns likely emulated much of the 1934 spatial structure by killing trees in patches [64].

4.3. Prescribed Burns Enhanced Spatial Variation in Both LOD and HID Treatments

The high structural diversity (HID) treatment for BMEF was designed to approximately emulate a heterogeneous overstory stand structure, with numerous small openings and patches of large old trees [59]. Therefore, our results regarding spatial variation in overstory AGB—presence of fine-scale spatial variation, periodicity in patchiness at short lags, and short ranges—for HID treatments suggest that treatments were effective in creating greater spatial variation of overstory stand structure [58]. Furthermore, when prescribed burns were coupled with the HID treatment, spatial variation in overstory AGB was further enhanced as evidenced by greater periodicity in overstory AGB as compared to the unburned half of the HID units.
Prescribed burns can be patchy and localized due to a slow rate of spread and thus can kill small trees in groups and create openings [64,65]. In addition, second-order post-fire mortality of trees can also remove patches of trees [65,66]. For example, low-severity prescribed burns can leave fire-injured conifers that are receptive hosts for bark beetles [66]. The subsequent mortality of weakened trees—especially large old ponderosa pine—due to beetle attacks can produce openings and gaps within stands [66,67], therefore creating a heterogeneous overstory similar to our HID burned halves. Our results in HID burned units are consistent with Dodson and Peterson [68], who found that thinning increased spatial aggregation of residual trees at fine scales (50–150 m) and that prescribed fire of different burning intensity further promoted a mosaic of gaps with burned and unburned trees followed by beetle-related mortality.
The LOD treatment was designed to create open, even-aged stand conditions with a single-layer canopy consisting of evenly spaced trees [59]. The largest and smallest trees were removed from these stands for a unimodal diameter distribution only leaving intermediate-sized trees [32]. Especially within the unburned halves of the treatments, the long ranges from the semivariograms and the uniform correlograms with few peaks at greater lag distances suggested that LOD unburned halves consist of a more spatially homogeneous overstory AGB compared to HID and LOD burned units. Our results from unburned halves of the LOD units are consistent with [69], who found that after uniform removal of trees, mixed-conifer stands showed little spatial variation at broad scales. Similarly little spatial variation in spruce-fir forests was found by Kuehne et al. [70] when using thinning treatments involving the removal of small and large dominant trees. Irrespective of both LOD and HID treatments, burned halves consistently showed more periodicity in overstory AGB, indicating that fire is a key component in enhancing spatial heterogeneity in overstory AGB.

4.4. Management Implications

Our findings about enhanced spatial variation after post-treatment prescribed burns are similar to what other studies have reported for dry forests across California (e.g., [7,71]). Our comparison of spatial variation between reference and contemporary forests provides a multi-faceted, quantitative approach for evaluating forests over different time periods with a variety of potential implications. First, the spatial variation of overstory AGB that we detected in the 1934 reference forests helps portray the types of variation associated with very few forest management activities (i.e., in the absence of harvest or prescribed burning) and thus can serve as a useful point of reference for managers. Then, in utilizing the values of ranges from the semivariograms across time, managers can understand the distance or scale over which spatial variation in the overstory in contemporary forests have potentially departed from that of reference conditions. Secondly, it is important to remember that these reference conditions represent the accumulated overstory stand structures of dry forests that have been influenced by multiple past disturbances, rather than a condition at a single point in time. Nonetheless, when examined over time, changes in spatial heterogeneity can be used to understand and disentangle the fundamental environmental factors driving forest overstory variation. We were unable to quantify changes in species composition, because our data from 1934 were not spatially explicit at the individual tree level. However, in future work the evaluation of species composition changes over time may reveal different spatial variation than we have detected here in both reference and contemporary dry forests.
In addition, our comparison of reference and contemporary overstory spatial variation can help managers to explore the impact of fire exclusion on driving overstory spatial variation in dry forests [56]. Furthermore, our deeper comparison between reference conditions and contemporary forests that subsequently received different management treatments provides information on which treatments are useful in dry conifer forests such as BMEF for emulating overstory spatial variation in reference conditions. For example, in managing contemporary forests, treatments similar to HID with prescribed burns that preserve and/or create patchy spatial variation similar to reference conditions, may be warranted. Similarly, LOD with prescribed burn may be warranted whenever the goal of management is to reduce canopy fuel to promote resiliency against fire. In addition, comparing spatial variation among various management treatments in contemporary forests will allow managers to understand the effectiveness of restoration treatments such as HID with prescribed burns in creating desirable spatial variation for contemporary dry forests.

5. Conclusions

This study provides a unique opportunity to examine the spatial variation of overstory stand structure in forest reference conditions using field-collected data as opposed to using tree size and density variables derived from dendro-chronological reconstructions. Furthermore, we were able to fill a gap in the literature related to understanding the spatial variation of overstory above-ground biomass between contemporary and reference forests.
The spatial variation of above-ground biomass in reference conditions described in this study comes from forests that have minimum human impact following Euro-American settlement and hence represent a relatively intact dry forest ecosystem. Therefore, knowledge about spatial variation in such intact forests in reference conditions provides insights into stand development, tree interactions with each other, regeneration, and mortality. Such information can help guide sustainable forest management in the face of growing natural environmental disturbances. Our results further indicate that low-severity fire seems to be key for emulating the mosaic of alternating patches of biomass at regular intervals throughout these dry forests. Therefore, managers seeking to enhance ecological resilience are advised to use prescribed burning alone or in combination with treatments that include some degree of thinning.

Author Contributions

Conceptualization, S.N., B.N.I.E. and M.W.R.; Data curation, S.N. and M.W.R.; Formal analysis, S.N.; Investigation, S.N., B.N.I.E. and M.W.R.; Methodology, S.N.; Project administration, B.N.I.E.; Resources, B.N.I.E. and M.W.R.; Supervision, B.N.I.E.; Visualization, S.N.; Writing—original draft, S.N.; Writing—review & editing, B.N.I.E., S.E.G. and M.W.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the USDA Forest Service, Pacific Southwest Research Station (grant number: 15-IJ-11272139-016) and UBC four-year fellowship (4YF) program.

Data Availability Statement

The 1934 data have been submitted to U.S. Forest Service Research Data Archive and will be publicly available there.

Acknowledgments

We thank Lassen National Forest and Pacific Southwest Research Station for providing the data for this study. Special thanks go to Ethan Hammett and Brian Wing for providing the LiDAR and GIS data layers and support on any related questions and Carl Skinner for the unpublished fire data. We would also like to thank Paul Hacker and Francisco Mauro for valuable insights during LiDAR data processing.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Data Summary

Table A1. Summary of the applied combination of structural diversity treatments and prescribed burn in various units of BMERP.
Table A1. Summary of the applied combination of structural diversity treatments and prescribed burn in various units of BMERP.
UnitStructural Treatments *Treatment YearPrescribed BurnBurn Year# of 1 ha Plots
RNA-ANo treatment-Unburned-55
RNA-BNo treatment-Burned 199763
RNA-CNo treatment-Burned 199776
RNA-DNo treatment-Unburned-35
UNIT-39LOD1996Burned 199789
UNIT-39LOD1996Unburned-85
UNIT-40LOD1998Burned 200087
UNIT40LOD1998Unburned-79
UNIT-43LOD1996Burned 199768
UNIT-43LOD1996Unburned-90
UNIT-44LOD1997Burned 199956
UNIT-44LOD1997Unburned-49
UNIT-45LOD1997Burned 199957
UNIT-45LOD1997Unburned-75
UNIT-38HID1996Burned 199787
UNIT-38HID1996Unburned-93
UNIT-41HID1996Burned 199782
UNIT-41HID1996Unburned-79
UNIT-42HID1997Burned 199983
UNIT-42HID1997Unburned-101
UNIT-47HID1997Burned 199959
UNIT-47HID1997Unburned-48
UNIT-48HID1998Burned 200081
UNIT-48HID1998Unburned-76
Total 1720
* HID: high structural diversity treatment; LOD: low structural diversity treatment; RNA: research natural areas; # = number.
Forests 14 02096 g0a1
Figure A1. Histogram for above-ground biomass (Mg/ha) for 2016 ground data. The red and blue lines indicate the mean and median, respectively.
Figure A1. Histogram for above-ground biomass (Mg/ha) for 2016 ground data. The red and blue lines indicate the mean and median, respectively.
Forests 14 02096 g0a1

Appendix B. LiDAR Point Cloud Processing and Model Building Process

The LiDAR point cloud was processed using the ‘lidR’ package in R [42]. We filtered duplicates for the overlapping region of the point clouds and then classified the ground returns [42]. All the ground returns were used to obtain a digital terrain model of 1 × 1 m grid cell size [42]. The point cloud was normalized by subtracting the digital terrain model from all returns to remove the influence of terrain on above-ground returns [42]. From the normalized point cloud, we used the area-based approach [42] to calculate the LiDAR metrics (Table A2)—hereafter referred to as auxiliary variables—using the ‘cloudmetrics’ function related to the 154 circular ground plots. We also extracted average elevation (Table A2) for each ground plot in 2016 from a digital elevation model at 1 m resolution provided by the Pacific Southwest Research Station, Redding, USA, using the zonal statistics function in ArcMap 10.4.4. [72].
Table A2. LiDAR-derived and topographic auxiliary variables were used for biomass modeling.
Table A2. LiDAR-derived and topographic auxiliary variables were used for biomass modeling.
GroupsLiDAR Variables Description
Group A
zmeanMean elevation for all first returns above 2 m
zq5Elevation for first return in 5th percentile
zq10Elevation for first return in 10th percentile
zq15Elevation for first return in 15th percentile
zq20Elevation for first return in 20th percentile
zq25Elevation for first return in 25th percentile
Group B
zq30Elevation for first return in 30th percentile
zq35Elevation for first return in 35th percentile
zq40Elevation for first return in 40th percentile
zq45Elevation for first return in 45th percentile
zq50Elevation for first return in 50th percentile
Group C
zq55Elevation for first return in 55th percentile
zq60Elevation for first return in 60th percentile
zq65Elevation for first return in 65th percentile
zq70Elevation for first return in 70th percentile
zq75Elevation for first return in 75th percentile
Group D
zq80Elevation for first return in 80th percentile
zq85Elevation for first return in 85th percentile
zq90Elevation for first return in 90th percentile
zq95Elevation for first return in 95th percentile
pzmeanPercentage of first return above mean
Elev (m)Elevation from digital elevation model @1 m
Treatment × prescribed burn

Appendix C. Model Building, Selection Process, and Biomass Prediction Summary

To avoid over fitting of the model and to have about 10–15 observations per auxiliary variable, we divided the auxiliary variables into four groups (Table A2, [43,73]). Within each variable group, following the top-down model-building strategy for glm models described by Zuur et al. [47], the LiDAR-derived variables and the six treatments—three structural diversity treatments crossed with two prescribed burn treatments—were used as fixed effects. A random unit effect was also included in the model [43]. The variables that were not statistically significant (p-value > 0.05) and had a variance inflation factor (VIF) of 10 or greater were dropped sequentially [47,74]. The variables that were retained from each auxiliary variable group were combined into a single model and were retained if they were statistically significant (p < 0.05) and had VIF < 10 [47,74]. Variables that were not statistically significant during group-wise selection were brought back into the model to check if the performance of the combined model improved. The model goodness-of-fit at each model building step was assessed using a graph showing observed vs. predicted values (Figure A2) and graphs of the residual deviance [47]. We selected the biomass model with the lowest AIC and lowest residual deviance (Table A3) as our final model. Our final model consisted of the mean height for all first returns above 2 m (zmean, p = 1.4 × 10−7), the height for the first return in the 5th percentile (zq5, p = 0.05), crossed treatments as fixed effects (p = 0.05), and unit level random effects (Table A3).
The range of predicted values of overstory AGB for the cells was within the range of the observed values for the circular ground plots (Table A4, Figure A3). The mean and median predicted values were very close to the mean and median observed values (Table A4, Figure A3). Overstory AGB predicted within the low structural diversity (LOD) units was low and consistent with observed overstory AGB (Figure A4). Predicted overstory AGB in high structural diversity (HID) units and research natural areas (RNAs) was high compared to LOD but within the range of observed overstory AGB for HID and RNA units (Figure A4).
Table A3. Summary of the best-selected models from each auxiliary variable group. Highlighted in bold is the final model selected for above-ground biomass (Mg/ha) prediction.
Table A3. Summary of the best-selected models from each auxiliary variable group. Highlighted in bold is the final model selected for above-ground biomass (Mg/ha) prediction.
Variable GroupModelAICResidual Deviancep-Value
Group AFixed effects = zmean + zq5 + Treatment × prescribed burn
Random effect = units		1577.91557.9zmean
(p < 0.0001)
zq5
(p = 0.00087)
Group BFixed effects = zq50+ Treatment × prescribed burn
Random effect = units		1590.011572.0zq50
(p = 0.00034)
Group CFixed effects = zq75+ Treatment × prescribed burn
Random effect = units		1582.91564.1zq75
(p < 0.0001)
Group DFixed effects = zq95+ Treatment × prescribed burn
Random effect = units		1581.91561.9zq95
(p < 0.0001)
FinalFixed effects = zmean + zq5 + Elev Random effects = Treatments × prescribed fire15,916.61577.6Zmean (<0.001)
Zq5 (<0.001)
Elev (<0.05)
Forests 14 02096 g0a2
Figure A2. Observed vs. predicted above-ground biomass (Mg/ha). The models shown above represent the best model from each of the variable groups along with the final model. The red line is the best fitted line.
Figure A2. Observed vs. predicted above-ground biomass (Mg/ha). The models shown above represent the best model from each of the variable groups along with the final model. The red line is the best fitted line.
Forests 14 02096 g0a2
Table A4. Summary statistics of the observed and predicted overstory above-ground biomass (Mg/ha) values for circular plots (n = 153) and cells (n = 12,647).
Table A4. Summary statistics of the observed and predicted overstory above-ground biomass (Mg/ha) values for circular plots (n = 153) and cells (n = 12,647).
Above-Ground Biomass (Mg/ha)MinMaxMeanMedian
Observed (Circular plots)4.33327.53101.7183.48
Predicted
(32 × 32 m cells)		10.33381.65101.0289.03
Forests 14 02096 g0a3
Figure A3. Predicted above-ground biomass (AGB, Mg/ha). The predictions were made at the cell level for overstory AGB and compared with the overstory AGB from the circular ground-truthing plots. The blue and red lines represent the median and mean of overstory AGB, respectively.
Figure A3. Predicted above-ground biomass (AGB, Mg/ha). The predictions were made at the cell level for overstory AGB and compared with the overstory AGB from the circular ground-truthing plots. The blue and red lines represent the median and mean of overstory AGB, respectively.
Forests 14 02096 g0a3
Forests 14 02096 g0a4
Figure A4. Predicted above-ground biomass (Mg/ha). The predictions were made using the final model and LiDAR derived metrics in 2016 at 32 m × 32 m pixel level for high structural diversity treatments (HID), low structural diversity treatments (LOD), and research natural areas (RNA) within Blacks Mountain Experimental Forest (BMEF).
Figure A4. Predicted above-ground biomass (Mg/ha). The predictions were made using the final model and LiDAR derived metrics in 2016 at 32 m × 32 m pixel level for high structural diversity treatments (HID), low structural diversity treatments (LOD), and research natural areas (RNA) within Blacks Mountain Experimental Forest (BMEF).
Forests 14 02096 g0a4

Appendix D. Additional Table and Figures for the Results

Table A5. Summary of root mean squared error (RMSE) and values for the spherical, exponential, and Gaussian models selected for burned and unburned halves of HID and LOD units.
Table A5. Summary of root mean squared error (RMSE) and values for the spherical, exponential, and Gaussian models selected for burned and unburned halves of HID and LOD units.
UnitFireModelRMSESill
Low structural Diversity (LOD)
UNIT-39BurnedGaussian0.540.86
UNIT-39UnburnedGaussian0.711.04
UNIT-40BurnedGaussian0.720.93
UNIT-40UnburnedGaussian0.731.05
UNIT-43BurnedSpherical0.751.03
UNIT-43UnburnedSpherical0.741.00
UNIT-44BurnedGaussian0.641.15
UNIT-44UnburnedGaussian0.490.92
UNIT-45BurnedGaussian0.581.14
UNIT-45UnburnedGaussian0.651.07
High structural diversity (HID)
UNIT-38BurnedGaussian0.650.92
UNIT-38UnburnedGaussian0.590.87
UNIT-41BurnedSpherical0.740.84
UNIT-41UnburnedSpherical0.760.98
UNIT-42BurnedExponential0.741.13
UNIT-42UnburnedExponential0.731.32
UNIT-47BurnedGaussian0.630.90
UNIT-47UnburnedGaussian0.682.81
UNIT-48BurnedGaussian0.700.92
UNIT-48UnburnedGaussian0.740.88
Forests 14 02096 g0a5
Figure A5. Variogram models for research natural areas (RNA) in 1934 and 2016.
Figure A5. Variogram models for research natural areas (RNA) in 1934 and 2016.
Forests 14 02096 g0a5
Forests 14 02096 g0a6
Figure A6. Variogram models for low structural diversity units (LOD) in 1934 and 2016.
Figure A6. Variogram models for low structural diversity units (LOD) in 1934 and 2016.
Forests 14 02096 g0a6
Forests 14 02096 g0a7
Figure A7. Variogram models for high structural diversity units (HID) in 1934 and 2016.
Figure A7. Variogram models for high structural diversity units (HID) in 1934 and 2016.
Forests 14 02096 g0a7
Forests 14 02096 g0a8
Figure A8. Moran’s I correlogram for 1934 vs. 2016 within different units under study in Blacks Mountain Experimental Forest. Dotted lines represent the lower and upper limit of the non-significant spatial autocorrelation using 95% confidence envelopes for 1934 (red) and 2016 (blue) from Monte Carlo simulations [56]. Points connected with solid lines indicate Moran’s I at a given lag distance (m). Points above the dotted lines in the upper part of the envelope indicate a positive Moran’s I that is significantly different from 0, whereas points below the dotted lines in the lower part of the envelope show a negative Moran’s I that is significantly different from 0. Points within the dotted envelope show Moran’s I values that are not significantly different from 0.
Figure A8. Moran’s I correlogram for 1934 vs. 2016 within different units under study in Blacks Mountain Experimental Forest. Dotted lines represent the lower and upper limit of the non-significant spatial autocorrelation using 95% confidence envelopes for 1934 (red) and 2016 (blue) from Monte Carlo simulations [56]. Points connected with solid lines indicate Moran’s I at a given lag distance (m). Points above the dotted lines in the upper part of the envelope indicate a positive Moran’s I that is significantly different from 0, whereas points below the dotted lines in the lower part of the envelope show a negative Moran’s I that is significantly different from 0. Points within the dotted envelope show Moran’s I values that are not significantly different from 0.
Forests 14 02096 g0a8
Forests 14 02096 g0a9
Figure A9. Variogram models for burned and unburned halves of low structural diversity units (LOD) in 2016.
Figure A9. Variogram models for burned and unburned halves of low structural diversity units (LOD) in 2016.
Forests 14 02096 g0a9
Forests 14 02096 g0a10
Figure A10. Variogram models for burned and unburned halves of high structural diversity units (HID) in 2016.
Figure A10. Variogram models for burned and unburned halves of high structural diversity units (HID) in 2016.
Forests 14 02096 g0a10
Forests 14 02096 g0a11
Figure A11. Moran’s I correlograms for burned and unburned halves within different HID and LOD units in Blacks Mountain Experimental Forest. Dotted lines represent the lower and upper limit of the non-significant spatial autocorrelation using 95% confidence envelopes for Burned (red) and Unburned (blue) from Monte Carlo simulations [56]. Points connected with solid lines indicate Moran’s I at a given lag distance (m). Points above the dotted lines in the upper part of the envelope indicate a positive Moran’s I that is significantly different from 0, whereas points below the dotted lines in the lower part of the envelope show a negative Moran’s I that is significantly different from 0. Points within the dotted envelope show Moran’s I values that are not significantly different from 0.
Figure A11. Moran’s I correlograms for burned and unburned halves within different HID and LOD units in Blacks Mountain Experimental Forest. Dotted lines represent the lower and upper limit of the non-significant spatial autocorrelation using 95% confidence envelopes for Burned (red) and Unburned (blue) from Monte Carlo simulations [56]. Points connected with solid lines indicate Moran’s I at a given lag distance (m). Points above the dotted lines in the upper part of the envelope indicate a positive Moran’s I that is significantly different from 0, whereas points below the dotted lines in the lower part of the envelope show a negative Moran’s I that is significantly different from 0. Points within the dotted envelope show Moran’s I values that are not significantly different from 0.
Forests 14 02096 g0a11

References

  1. Harrod, R.J.; Gaines, W.L.; Hartl, W.E.; Camp, A. Estimating Historical Snag Density in Dry Forests East of the Cascade Range; U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station: Portland, OR, USA, 1998; p. 116. [Google Scholar]
  2. Dolanc, C.R.; Safford, H.D.; Thorne, J.H.; Dobrowski, S.Z. Changing Forest Structure Across the Landscape of the Sierra Nevada, CA, USA, since the 1930s. Ecosphere 2014, 5, 1–26. [Google Scholar] [CrossRef]
  3. Lydersen, J.M.; North, M.P.; Knapp, E.E.; Collins, B.M. Quantifying Spatial Patterns of Tree Groups and Gaps in Mixed-Conifer Forests: Reference Conditions and Long-Term Changes Following Fire Suppression and Logging. For. Ecol. Manag. 2013, 304, 370–382. [Google Scholar] [CrossRef]
  4. North, M.; Stine, P.; O’Hara, K.; Zielinski, W.; Stephens, S. An Ecosystem Management Strategy for Sierran Mixed-Conifer Forests; U.S. Department of Agriculture, Forest Service, Pacific Southwest Research Station: Albany, CA, USA, 2009. [Google Scholar]
  5. van Mantgem, P.J.; Stephenson, N.L. Apparent Climatically Induced Increase of Tree Mortality Rates in a Temperate Forest. Ecol. Lett. 2007, 10, 909–916. [Google Scholar] [CrossRef] [PubMed]
  6. Abella, S.R.; Denton, C.W. Spatial Variation in Reference Conditions: Historical Tree Density and Pattern on a Pinus Ponderosa Landscape. Can. J. For. Res. 2009, 39, 2391–2403. [Google Scholar] [CrossRef]
  7. Churchill, D.J.; Larson, A.J.; Dahlgreen, M.C.; Franklin, J.F.; Hessburg, P.F.; Lutz, J.A. Restoring Forest Resilience: From Reference Spatial Patterns to Silvicultural Prescriptions and Monitoring. For. Ecol. Manag. 2013, 291, 442–457. [Google Scholar] [CrossRef]
  8. Moore, D.A.; Carpenter, T.E. Spatial Analytical Methods and Geographic Information Systems: Use in Health Research and Epidemiology. Epidemiol. Rev. 1999, 21, 143–161. [Google Scholar] [CrossRef] [PubMed]
  9. Youngblood, A.; Max, T.; Coe, K. Stand Structure in Eastside Old-Growth Ponderosa Pine Forests of Oregon and Northern California. For. Ecol. Manag. 2004, 199, 191–217. [Google Scholar] [CrossRef]
  10. Clyatt, K.A.; Keyes, C.R.; Hood, S.M. Long-Term Effects of Fuel Treatments on Aboveground Biomass Accumulation in Ponderosa Pine Forests of the Northern Rocky Mountains. For. Ecol. Manag. 2017, 400, 587–599. [Google Scholar] [CrossRef]
  11. Collins, B.M.; Fry, D.L.; Lydersen, J.M.; Everett, R.; Stephens, S.L. Impacts of Different Land Management Histories on Forest Change. Ecol. Appl. 2017, 27, 2475–2486. [Google Scholar] [CrossRef]
  12. Oswalt, S.N.; Smith, W.B.; Miles, P.D.; Pugh, S.A. Forest Resources of the United States, 2012: A Technical Document Supporting the Forest Service Update of the 2010 RPA Assessment; U.S. Department of Agriculture, Forest Service, Washington Office: Washington, DC, USA, 2014; 218p. [Google Scholar]
  13. Kilgore, B.M. The Ecological Role of Fire in Sierran Conifer Forests. Its Application to National Park Management. Quat. Res. 1973, 3, 496–513. [Google Scholar] [CrossRef]
  14. Agee, J.K.; Skinner, C.N. Basic Principles of Forest Fuel Reduction Treatments. For. Ecol. Manag. 2005, 211, 83–96. [Google Scholar] [CrossRef]
  15. Taylor, A.H. Fire Disturbance and Forest Structure in an Old-Growth Pinus Ponderosa Forest, Southern Cascades, USA. J. Veg. Sci. 2010, 21, 561–572. [Google Scholar] [CrossRef]
  16. Chen, J.; Song, B.; Rudnicki, M.; Moeur, M.; Bible, K.; North, M.; Shaw, D.C.; Franklin, J.F.; Braun, D.M. Spatial Relationship of Biomass and Species Distribution in an Old-Growth. For. Sci. 2004, 50, 364–375. [Google Scholar]
  17. Kashian, D.M.; Romme, W.H.; Tinker, D.B.; Turner, M.G.; Ryan, M.G. Carbon Storage on Landscapes with Stand-Replacing Fires. Bioscience 2006, 56, 598–606. [Google Scholar] [CrossRef]
  18. Hall, S.A.; Burke, I.C.; Box, D.O.; Kaufmann, M.R.; Stoker, J.M. Estimating Stand Structure Using Discrete-Return Lidar: An Example from Low Density, Fire Prone Ponderosa Pine Forests. For. Ecol. Manag. 2005, 208, 189–209. [Google Scholar] [CrossRef]
  19. McRoberts, R.E.; Tomppo, E.O.; Næsset, E. Advances and Emerging Issues in National Forest Inventories. Scand. J. For. Res. 2010, 25, 368–381. [Google Scholar] [CrossRef]
  20. Lefsky, M.A.; Cohen, W.B.; Acker, S.A.; Parker, G.G.; Spies, T.A.; Harding, D. Lidar Remote Sensing of the Canopy Structure and Biophysical Properties of Douglas-Fir Western Hemlock Forests. Remote Sens. Environ. 1999, 70, 339–361. [Google Scholar] [CrossRef]
  21. Skowronski, N.S.; Lister, A.J. The Utility of Lidar for Large Area Forest Inventory Applications; U.S. Department of Agriculture, Forest Service, Northern Research Station: Newtown Square, PA, USA, 2012. [Google Scholar]
  22. Millar, C.I.; Stephenson, N.L.; Stephens, S.L. Climate Change and Forests of the Future: Managing in the Face of Uncertainity. Ecol. Appl. 2007, 17, 2145–2151. [Google Scholar] [CrossRef]
  23. Vankat, J.L.; Major, J. Vegetation Changes in Sequoia National Park, California. J. Biogeogr. 1978, 5, 377–402. [Google Scholar] [CrossRef]
  24. Beaty, R.M.; Taylor, A.H. Fire History and the Structure and Dynamics of a Mixed Conifer Forest Landscape in the Northern Sierra Nevada, Lake Tahoe Basin, California, USA. For. Ecol. Manag. 2008, 255, 707–719. [Google Scholar] [CrossRef]
  25. Ansley, J.S.; Battles, J.J. Forest Composition, Structure, and Change in an Old-Growth Mixed Confer Forest in the Northern Sierra Nevada. J. Torrey Bot. Soc. 1998, 125, 297–308. [Google Scholar] [CrossRef]
  26. Smith, T.F.; Rizzo, D.M.; North, M.P. Patterns of Mortality in an Old-Growth Mixed-Conifer Forest of the Southern Sierra Nevada, Californi. For. Sci. 2005, 51, 266–275. [Google Scholar]
  27. Knight, C.A.; Cogbill, C.V.; Potts, M.D.; Wanket, J.A.; Battles, J.J. Settlement-Era Forest Structure and Composition in the Klamath Mountains: Reconstructing a Historical Baseline. Ecosphere 2020, 11, e03250. [Google Scholar] [CrossRef]
  28. Fry, D.L.; Stephens, S.L.; Collins, B.M.; North, M.P.; Franco-Vizcaíno, E.; Gill, S.J. Contrasting Spatial Patterns in Active-Fire and Fire-Suppressed Mediterranean Climate Old-Growth Mixed Conifer Forests. PLoS ONE 2014, 9, e88985. [Google Scholar] [CrossRef] [PubMed]
  29. Sanchez Meador, A.J.; Parysow, P.F.; Moore, M.M. Historical Stem-Mapped Permanent Plots Increase Precision of Reconstructed Reference Data in Ponderosa Pine Forests of Northern Arizona. Restor. Ecol. 2010, 18, 224–234. [Google Scholar] [CrossRef]
  30. Sheridan, R.D.; Popescu, S.C.; Gatziolis, D.; Morgan, C.L.S.; Ku, N.W. Modeling Forest Aboveground Biomass and Volume Using Airborne LiDAR Metrics and Forest Inventory and Analysis Data in the Pacific Northwest. Remote Sens. 2015, 7, 229–255. [Google Scholar] [CrossRef]
  31. Richardson, J.J.; Moskal, L.M. Strengths and Limitations of Assessing Forest Density and Spatial Configuration with Aerial LiDAR. Remote Sens. Environ. 2011, 115, 2640–2651. [Google Scholar] [CrossRef]
  32. Oliver, W. Ecological Research at the Blacks Mountain Experimental Forest in Northeastern California; PSW-GTR-179; Pacific Southwest Research Station, Forest Service, U.S. Department of Agriculture: Albany, CA, USA, 2000; 66p. [Google Scholar]
  33. Gordon, D.T. History of the Blacks Mountain Experimental Forest, 1933 through 1981; U.S. Department of Agriculture, Forest Service, Pacific Southwest Forest and Range Experiment Station: Berkeley, CA, USA, 1981; pp. 1–209. [Google Scholar]
  34. Jia, G.J.; Burke, I.C.; Goetz, A.F.H.; Kaufmann, M.R.; Kindel, B.C. Assessing Spatial Patterns of Forest Fuel Using AVIRIS Data. Remote Sens. Environ. 2006, 102, 318–327. [Google Scholar] [CrossRef]
  35. Keane, R.E.; Nurgan, R.E.; Wagtendonk, J.W.V. Mapping Wildland Fuels for Fire Management across Multiple Scales: Integrating Remote Sensing, GIS, and Biophysical. Int. J. Wild. Fire 2001, 10, 301–319. [Google Scholar] [CrossRef]
  36. Rollins, M.G.; Keane, R.E.; Parson, R.A. Mapping Fuels and Fire Regimes Using Remote Sensing, Ecosystem Simulation, and Gradient Modeling. Ecol. Appl. 2004, 14, 75–95. [Google Scholar] [CrossRef]
  37. Swetnam, T.W. Fire History and Climate Chage in Giant Sequoia Groves. Science 1993, 262, 885–889. [Google Scholar] [CrossRef] [PubMed]
  38. Eyre, F.H. Forest Cover Types of the United States and Canada; U.S. Depratment of Agriculture, Forest Service, Washington Office: Washington, DC, USA, 1980. [Google Scholar]
  39. Ritchie, M.W. Multi-Scale Reference Conditions in an Interior Pine-Dominated Landscape in Northeastern California. For. Ecol. Manag. 2016, 378, 233–243. [Google Scholar] [CrossRef]
  40. Hasel, A.A. Instruction for Type Map and Inventory of Experimental Forest; Pacific Southwest Research Station: Berkeley, CA, USA, 1935; 9p. [Google Scholar]
  41. Wing, B.M.; Ritchie, M.W.; Boston, K.; Cohen, W.B.; Gitelman, A.; Olsen, M.J. Prediction of Understory Vegetation Cover with Airborne Lidar in an Interior Ponderosa Pine Forest. Remote Sens. Environ. 2012, 124, 730–741. [Google Scholar] [CrossRef]
  42. Roussel, J.R.; Auty, D.; Coops, N.C.; Tompalski, P.; Goodbody, T.R.H.; Meador, A.S.; Bourdon, J.F.; de Boissieu, F.; Achim, A. LidR: An R Package for Analysis of Airborne Laser Scanning (ALS) Data. Remote Sens. Environ. 2020, 251, 112061. [Google Scholar] [CrossRef]
  43. Mauro, F.; Ritchie, M.; Wing, B.; Frank, B.; Monleon, V.; Temesgen, H.; Hudak, A. Estimation of Changes of Forest Structural Attributes at Three Di Ff Erent Spatial Aggregation Levels in Northern California Using Multitemporal LiDAR. Remote Sens. 2019, 11, 923. [Google Scholar] [CrossRef]
  44. Ritchie, M.W.; Zhang, J.; Hamilton, T.A. Aboveground Tree Biomass for Pinus Ponderosa in Northeastern California. Forests 2013, 4, 179–196. [Google Scholar] [CrossRef]
  45. Leroy Dolph, K.; Mori, S.R.; Oliver, W.W. Height-Diameter Relationship for Conifer Species on the Blacks Mountain Experimental Forest; Research Paper; U.S. Department of Agriculture, Forest Service: Washington, DC, USA, 1995; p. 5. [Google Scholar]
  46. Jenkins, J.C.; Chojnacky, D.C.; Heath, L.S.; Birdsey, R.A. Comprehensive Database of Diameter-Based Biomass Resgressions for North American Tree Species; U.S. Department of Agriculture, Forest Service, Northeastern Research Station: Newtown Square, PA, USA, 2003. [Google Scholar]
  47. Zuur, A.F.; Ieno, E.N.; Walker, N.J.; Saveliev, A.A.; Smith, G.M. Mixed Effects Models and Extensions in Ecology with R; Springer: New York, NY, USA, 2009; p. 579. [Google Scholar]
  48. Bates, D.; Mächler, M.; Bolker, B.; Walker, S. Fitting Linear Mixed-Effects Models Using Lme4. J. Stat. Softw. 2015, 67, 1–48. [Google Scholar] [CrossRef]
  49. Reed, T.W.; Gulland, E.; West, G.; Mcmeekin, D.A.; Moncrieff, S. Geographic Metadata Searching with Semantic and Spatial Filtering Methods. In Proceedings of the GEOProcessing 2016: The Eighth International Conference on Advanced Geographic Information Systems, Applications, and Services, Venice, Italy, 24–28 April 2016. [Google Scholar]
  50. Rossi, R.E.; Mulla, D.J.; Journel, A.G.; Franz, E.H. Geostatistical Tools for Modeling and Interpreting Ecological Spatial Dependence. Ecol. Monogr. 1992, 62, 277–314. [Google Scholar] [CrossRef]
  51. Isaaks, E.; Srivastava, R.M. An Introduction to Applied Geostatistics, 1st ed.; Oxford University Press: New York, NY, USA, 1989; Volume 17. [Google Scholar]
  52. Pebesma, E.J. Multivariable Geostatistics in S: The Gstat Package. Comput. Geosci. 2004, 30, 683–691. [Google Scholar] [CrossRef]
  53. Fry, D.L.; Stephens, S.L. Stand-Level Spatial Dependence in an Old-Growth Jeffrey Pine-Mixed Conifer Forest, Sierra San Pedro Mártir, Mexico. Can. J. For. Res. 2010, 40, 1803–1814. [Google Scholar] [CrossRef]
  54. Keane, R.E.; Gray, K.; Bacciu, V.; Leirfallom, S. Spatial Scaling of Wildland Fuels for Six Forest and Rangeland Ecosystems of the Northern Rocky Mountains, USA. Landsc. Ecol. 2012, 27, 1213–1234. [Google Scholar] [CrossRef]
  55. Zawadzki, J.; Cieszewski, C.J.; Zasada, M.; Lowe, R.C. Applying Geostatistics for Investigations of Forest Ecosystems Using Remote Sensing Imagery. Silva Fenn. Monogr. 2005, 39, 599–617. [Google Scholar] [CrossRef]
  56. Jaquette, M.; Sánchez Meador, A.J.; Huffman, D.W.; Bowker, M.A. Mid-Scale Drivers of Variability in Dry Mixed-Conifer Forests of the Mogollon Rim, Arizona. Forests 2021, 12, 622. [Google Scholar] [CrossRef]
  57. Legendre, P.; Fortin, M.-J. Spatial Pattern and Ecological Analysis. Vegetatio 1989, 80, 107–138. [Google Scholar] [CrossRef]
  58. Radeloff, V.C.; Miller, T.F.; He, H.S.; Mladenoff, D.J. Periodicity in Spatial Data and Geostatistical Models: Autocorrelation between Patches. Ecography 2000, 23, 81–91. [Google Scholar] [CrossRef]
  59. Crotteau, J.S.; Ritchie, M.W. Long-Term Stand Growth of Interior Ponderosa Pine Stands in Response to Structural Modifications and Burning Treatments in Northeastern California. J. For. 2014, 112, 412–423. [Google Scholar] [CrossRef]
  60. Arthur, G.; Janet, F. Second-Order Neighborhood Analysis of Mapped Point Patterns. Ecology 1987, 68, 473–477. [Google Scholar] [CrossRef]
  61. Freeman, E.A.; Moisen, G.G. Evaluating Kriging as a Tool to Improve Moderate Resolution Maps of Forest Biomass. Environ. Monit. Assess. 2007, 128, 395–410. [Google Scholar] [CrossRef]
  62. Ziegler, J.P.; Hoffman, C.M.; Fornwalt, P.J.; Sieg, C.H.; Battaglia, M.A.; Chambers, M.E.; Iniguez, J.M. Tree Regeneration Spatial Patterns in Ponderosa Pine Forests Following Stand-Replacing Fire: Influence of Topography and Neighbors. Forests 2017, 8, 391. [Google Scholar] [CrossRef]
  63. North, M.P.; Tompkins, R.E.; Bernal, A.A.; Collins, B.M.; Stephens, S.L.; York, R.A. Forest Ecology and Management Operational Resilience in Western US Frequent-Fire Forests. For. Ecol. Manag. 2022, 507, 120004. [Google Scholar] [CrossRef]
  64. Kerns, B.K.; Thies, W.G.; Niwa, C.G. Season and Severity of Prescribed Burn in Ponderosa Pine Forests: Implications for Understory Native and Exotic Plants. Ecoscience 2006, 13, 44–55. [Google Scholar] [CrossRef]
  65. Sackett, S.S.; Haase, S.M.; Harrington, M.G. Prescribed Burning in Southwestern Ponderosa Pine. Eff. Fire Madrean Prov. Ecosyst. A Symp. Proc. 1996, 289, 178–186. [Google Scholar]
  66. Hood, S.M.; Smith, S.L.; Cluck, D.R. Predicting Mortality for Five California Conifers Following Wildfire. For. Ecol. Manag. 2010, 260, 750–762. [Google Scholar] [CrossRef]
  67. Davis, R.S.; Hood, S.; Bentz, B.J. Fire-Injured Ponderosa Pine Provide a Pulsed Resource for Bark Beetles. Can. J. For. Res. 2012, 42, 2022–2036. [Google Scholar] [CrossRef]
  68. Dodson, E.K.; Peterson, D.W. Dry Coniferous Forest Restoration and Understory Plant Diversity: The Importance of Community Heterogeneity and the Scale of Observation. For. Ecol. Manag. 2010, 260, 1702–1707. [Google Scholar] [CrossRef]
  69. Larson, A.J.; Stover, K.C.; Keyes, C.R. Effects of Restoration Thinning on Spatial Heterogeneity in Mixed-Conifer Forest. Can. J. For. Res. 2012, 42, 1505–1517. [Google Scholar] [CrossRef]
  70. Kuehne, C.; Weiskittel, A.; Pommerening, A.; Wagner, R.G. Evaluation of 10-Year Temporal and Spatial Variability in Structure and Growth across Contrasting Commercial Thinning Treatments in Spruce-Fir Forests of Northern Maine, USA. Ann. For. Sci. 2018, 75, 20. [Google Scholar] [CrossRef]
  71. Knapp, E.E.; Lydersen, J.M.; North, M.P.; Collins, B.M. Efficacy of Variable Density Thinning and Prescribed Fire for Restoring Forest Heterogeneity to Mixed-Conifer Forest in the Central Sierra Nevada, CA. For. Ecol. Manag. 2017, 406, 228–241. [Google Scholar] [CrossRef]
  72. Wilson, J.P. GIScience Research at the Thirty-Second Annual Esri International User Conference. Trans. GIS 2012, 16, 267–269. [Google Scholar] [CrossRef]
  73. Harrell, F.E.J.; Lee, K.L.; Califf, R.M.; Pryor, D.B.; Rosati, R.A. Regression Modelling Strategies for Improved Prognostic Prediction. Stat. Med. 1984, 3, 143–152. [Google Scholar] [CrossRef]
  74. Sheskin, D.J. Handbook of Parametric and Nonparametric Statistical Procedures, 5th ed.; Chapman & Hall: NewYork, NY, USA, 2000; p. 1928. [Google Scholar] [CrossRef]
Forests 14 02096 g001
Figure 1. Study area map. Panel (A) represents the United States of America with California being highlighted in grey. California with its county boundaries is shown in Panel (B). Blacks Mountain Experimental Forest is represented in Panel (C) with the 1 ha plots used for the 1934 census survey. Panel (C) also includes the polygons for the Blacks Mountain Ecological Research project (BMERP). Inside the BMERP polygons, the blue dots represent the permanent grid system within BMEF as shown in panel (D). The ground-truthing data were collected using a staggered pattern in 2016, as shown in Panel (D).
Figure 1. Study area map. Panel (A) represents the United States of America with California being highlighted in grey. California with its county boundaries is shown in Panel (B). Blacks Mountain Experimental Forest is represented in Panel (C) with the 1 ha plots used for the 1934 census survey. Panel (C) also includes the polygons for the Blacks Mountain Ecological Research project (BMERP). Inside the BMERP polygons, the blue dots represent the permanent grid system within BMEF as shown in panel (D). The ground-truthing data were collected using a staggered pattern in 2016, as shown in Panel (D).
Forests 14 02096 g001
Forests 14 02096 g002
Figure 2. Range values for above-ground biomass for the selected variogram models between 1934 vs. 2016. The top panel shows the range for research natural areas (RNA), the middle panel shows the range for the low structural diversity (LOD) units, and the bottom panel shows the range for the high structural diversity (HID) units in Blacks Mountain Ecological Research Project.
Figure 2. Range values for above-ground biomass for the selected variogram models between 1934 vs. 2016. The top panel shows the range for research natural areas (RNA), the middle panel shows the range for the low structural diversity (LOD) units, and the bottom panel shows the range for the high structural diversity (HID) units in Blacks Mountain Ecological Research Project.
Forests 14 02096 g002
Forests 14 02096 g003
Figure 3. Nugget values for above-ground biomass for the selected variogram models for 1934 and 2016. The top panel shows the nuggets for research natural areas (RNA), the middle panel shows the nuggets for the low structural diversity (LOD) units, and the bottom panel shows the nuggets for the high structural diversity (HID) units in Blacks Mountain Ecological Research Project.
Figure 3. Nugget values for above-ground biomass for the selected variogram models for 1934 and 2016. The top panel shows the nuggets for research natural areas (RNA), the middle panel shows the nuggets for the low structural diversity (LOD) units, and the bottom panel shows the nuggets for the high structural diversity (HID) units in Blacks Mountain Ecological Research Project.
Forests 14 02096 g003
Forests 14 02096 g004
Figure 4. Moran’s I correlogram for 1934 vs. 2016 within the four research natural areas (RNAs) in Blacks Mountain Ecological Research Project. RNA-B and RNA-C received prescribed burn treatments in 1997 and 1998, respectively. Dotted lines represent the lower and upper limit of the non-significant spatial autocorrelation using 95% confidence envelopes for 1934 (red) and 2016 (blue) from Monte Carlo simulations [56]. Points connected with solid lines indicate Moran’s I at a given lag distance (m). Points above the dotted lines in the upper part of the envelope indicate a positive Moran’s I that is significantly different from 0, whereas points below the dotted lines in the lower part of the envelope show a negative Moran’s I that is significantly different from 0. Points within the dotted envelope show Moran’s I values that are not significantly different from 0.
Figure 4. Moran’s I correlogram for 1934 vs. 2016 within the four research natural areas (RNAs) in Blacks Mountain Ecological Research Project. RNA-B and RNA-C received prescribed burn treatments in 1997 and 1998, respectively. Dotted lines represent the lower and upper limit of the non-significant spatial autocorrelation using 95% confidence envelopes for 1934 (red) and 2016 (blue) from Monte Carlo simulations [56]. Points connected with solid lines indicate Moran’s I at a given lag distance (m). Points above the dotted lines in the upper part of the envelope indicate a positive Moran’s I that is significantly different from 0, whereas points below the dotted lines in the lower part of the envelope show a negative Moran’s I that is significantly different from 0. Points within the dotted envelope show Moran’s I values that are not significantly different from 0.
Forests 14 02096 g004
Forests 14 02096 g005
Figure 5. Differences in range values for above-ground biomass for the selected variogram models in burned and unburned halves. The top panel shows the range values for low structural diversity (LOD) units, whereas the bottom panel shows the range values for high structural diversity (HID) units in the Blacks Mountain Ecological Research Project.
Figure 5. Differences in range values for above-ground biomass for the selected variogram models in burned and unburned halves. The top panel shows the range values for low structural diversity (LOD) units, whereas the bottom panel shows the range values for high structural diversity (HID) units in the Blacks Mountain Ecological Research Project.
Forests 14 02096 g005
Forests 14 02096 g006
Figure 6. Differences in sill values for above-ground biomass for the selected variogram models in burned and unburned halves. The top panel shows the sill values for low structural diversity (LOD) units, whereas the bottom panel shows the sill values for high structural diversity (HID) units in the Blacks Mountain Ecological Research Project.
Figure 6. Differences in sill values for above-ground biomass for the selected variogram models in burned and unburned halves. The top panel shows the sill values for low structural diversity (LOD) units, whereas the bottom panel shows the sill values for high structural diversity (HID) units in the Blacks Mountain Ecological Research Project.
Forests 14 02096 g006
Table 1. Summary of the root mean square error (RMSE), and nugget values for the models that converged. The selected variogram models for each unit are highlighted with bold font.
Table 1. Summary of the root mean square error (RMSE), and nugget values for the models that converged. The selected variogram models for each unit are highlighted with bold font.
UnitsModelsRMSESill
Research natural areas (RNA) 1934201619342016
RNA-ASpherical0.810.851.151.09
Exponential0.820.841.911.75
Gaussian0.740.801.010.99
RNA-BSpherical0.830.961.211.28
Exponential0.820.952.693.84
Gaussian0.790.941.090.99
RNA-CSpherical0.630.800.931.11
Exponential0.720.821.011.49
RNA-DGaussian0.500.651.011.56
Low structural diversity (LOD)
UNIT-39Spherical0.520.740.991.06
Exponential0.570.741.941.92
Gaussian0.490.700.950.98
UNIT-40Spherical0.760.731.041.11
Exponential0.750.751.921.94
UNIT-43Spherical0.650.701.030.91
Exponential0.690.711.161.06
Gaussian0.630.670.980.87
UNIT-44Spherical0.300.511.420.96
Exponential0.310.532.621.47
Gaussian0.280.480.800.76
UNIT-45Spherical0.690.802.501.12
Gaussian0.640.801.101.27
High structural diversity (HID)
UNIT-38Spherical0.480.670.810.95
Exponential0.520.690.851.07
Gaussian0.470.640.790.89
UNIT-41Spherical0.720.711.030.94
Exponential0.750.711.201.35
UNIT-42Spherical0.630.710.980.95
Exponential0.860.841.542.52
UNIT-47Spherical0.850.841.141.18
Exponential0.840.860.921.55
UNIT-48Spherical0.640.750.831.04
Exponential0.660.750.000.00
Table 2. Summary of the number of peaks and troughs of the Moran’s I correlograms and p-values for the Moran’s I at the first lag for research natural areas (RNAs) and other units in 1934 and 2016.
Table 2. Summary of the number of peaks and troughs of the Moran’s I correlograms and p-values for the Moran’s I at the first lag for research natural areas (RNAs) and other units in 1934 and 2016.
Unit First Lag
1934201619342016
# Peaks# Peaks
Research natural area (RNA)
RNA-A53p = 0.009 *p = 0.009 *
RNA-B33p = 0.23 nsp = 0.35 ns
RNA-C3 5p = 0.13 nsp = 0.06 *
RNA-D13p = 0.16 nsp = 0.47 ns
Low structural diversity (LOD)
UNIT-3911p = 0.009 *p = 0.009 *
UNIT-4022p = 0.06 nsp = 0.02 *
UNIT-4324p = 0.13 nsp = 0.009 *
UNIT-4402p = 0.16 nsp = 0.03 *
UNIT-4511p = 0.009 *p = 0.009 *
High Structural diversity (HID)
UNIT-3822p = 0.009 *p = 0.009 *
UNIT-4124p = 0.06 nsp = 0.17 ns
UNIT-4202p = 0.72 nsp = 0.93 ns
UNIT-4711p = 0.75 nsp = 0.05 *
UNIT-4822p = 0.009 *p = 0.009 *
* indicates significant at alpha = 0.05; ns = non-significant; # = number.
Table 3. Summary of the number of peaks and troughs of Moran’s I correlogram and p-values for Moran’s I at first lag for burned and unburned halves of low structural diversity treatment (LOD) and high structural diversity treatment (HID) units.
Table 3. Summary of the number of peaks and troughs of Moran’s I correlogram and p-values for Moran’s I at first lag for burned and unburned halves of low structural diversity treatment (LOD) and high structural diversity treatment (HID) units.
Unit First Lag
BurnedUnburned BurnedUnburned
# Peak# Peaks
Low structural diversity
UNIT-3932p = 0.009 *p = 0.009 *
UNIT-4054p = 0.11 nsp = 0.009 *
UNIT-4352p = 0.009 *p = 0.009 *
UNIT-4432p = 0.005 *p = 0.009 *
UNIT-4565p = 0.009 *p = 0.009 *
High structural diversity
UNIT-3883p = 0.009 *p = 0.002 *
UNIT-4143p = 0.007 *p = 0.27 ns
UNIT-4244p = 0.009 *p = 0.009 *
UNIT-4753p = 0.009 *p = 0.009 *
UNIT-4853p = 0.009 *p = 0.009 *
* indicates significant at alpha = 0.05; ns = non-significant; # = number.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Nepal, S.; Eskelson, B.N.I.; Ritchie, M.W.; Gergel, S.E. Thinning Combined with Prescribed Burn Created Spatially Heterogeneous Overstory Structures in Contemporary Dry Forests: A Comparison Using LiDAR (2016) and Field Inventory (1934) Data. Forests 2023, 14, 2096. https://0-doi-org.brum.beds.ac.uk/10.3390/f14102096

AMA Style

Nepal S, Eskelson BNI, Ritchie MW, Gergel SE. Thinning Combined with Prescribed Burn Created Spatially Heterogeneous Overstory Structures in Contemporary Dry Forests: A Comparison Using LiDAR (2016) and Field Inventory (1934) Data. Forests. 2023; 14(10):2096. https://0-doi-org.brum.beds.ac.uk/10.3390/f14102096

Chicago/Turabian Style

Nepal, Sushil, Bianca N. I. Eskelson, Martin W. Ritchie, and Sarah E. Gergel. 2023. "Thinning Combined with Prescribed Burn Created Spatially Heterogeneous Overstory Structures in Contemporary Dry Forests: A Comparison Using LiDAR (2016) and Field Inventory (1934) Data" Forests 14, no. 10: 2096. https://0-doi-org.brum.beds.ac.uk/10.3390/f14102096

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop