Forestry, Faculty of
DSpace
UBCV
Nepal, Sushil
Doctor of Philosophy - PhD
University of British Columbia
Spatial patterns of forest overstory prior to timber harvest activities associated with Euro-American settlement can serve as helpful forest reference conditions. Yet, such historical reference patterns can be challenging to quantify due to a lack of appropriate spatially explicit field data and the limitations of dendrochronological techniques. As a result, many studies have investigated contemporary spatial patterns of tree density and tree size, but few have successfully described changes in spatial patterns. In this dissertation, I first used overstory census data collected in 1934, over an area of ~4,000 ha, to characterize forest overstory reference conditions for Blacks Mountain Experimental Forest (BMEF), California, USA. Specifically, I quantified the 1934 spatial patterns of stand density index by seven vigor classes for Ponderosa pine/Jeffrey pine and white fir and incense-cedar. The forest reference conditions for BMEF consisted of clusters of open forests dominated by high vigor ponderosa pine trees associated with soil moisture and elevation gradients. Second, I used both 1934 and 2018 regeneration data from BMEF to compare the stocking between reference and contemporary conditions. The stocking of the reference conditions was dominated by ponderosa pine seedlings and was associated with available water capacity at 150 cm and overstory basal area (m2 ha-1). Third, I utilized 2016 Light Detection and Ranging data and overstory ground truth data to compare among stands that received thinning treatments, prescribed burn treatments, and a combination of thinning with prescribed burn treatments. I then compared above-ground biomass (Mg ha-1) between contemporary (2016) and reference (1934 census data) conditions. Spatial patterns in above-ground biomass were more heterogeneous in reference conditions than in contemporary conditions. However, in contemporary conditions, prescribed burns combined with thinning were better able to emulate the spatial patterns of reference conditions than thinning treatments alone. The findings of this dissertation can be utilized by managers: a) as a reference point against which they can evaluate the success of restoration treatments, b) to understand the relationships between regeneration and environmental factors, which shape the overall stand structure, and c) to identify contemporary stands at-risk of insect, disease, and fire that may benefit from restoration treatments.
SPATIAL PATTERNS OF STAND STRUCTURE IN PINUS PONDEROSA DOMINATED FORESTS WITHIN THE BLACKS MOUNTAIN EXPERIMENTAL FOREST, CALIFORNIA, USA, IN 1930s AND CONTEMPORARY CONDITIONS by Sushil Nepal B.Sc., Institute of Agriculture and Animal Science, 2011 M.Sc., Institute of Agriculture and Animal Sciences, 2014 M.Sc., Forestry, Northern Arizona University, 2017 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (FORESTRY) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) June 2023 © Sushil Nepal, 2023 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled: Spatial patterns of stand structure in Pinus ponderosa dominated forests within the Blacks Mountain Experimental Forest, California, USA, in 1930s and contemporary conditions submitted by Sushil Nepal in partial fulfilment of the requirements for the degree of Doctor of Philosophy in Forestry Examining Committee: Bianca N.I. Eskelson, Associate Professor, Forest Resources Management, UBC Supervisor Martin W. Ritchie, Research Forester, United States Department of Agriculture, Forest Service, Pacific Southwest Research Station Supervisory Committee Member Nicholas C. Coops, Professor, Forest Resources Management, UBC Supervisory Committee Member Sarah E. Gergel, Professor, Landscape Ecology and Conservation, UBC Supervisory Committee Member Sean M. Smukler, Associate Professor, Faculty of Land and Food Systems, UBC University Examiner Jeanine M. Rhemtulla, Associate Professor, Forest and Conservation Sciences, UBC University Examiner iii Abstract Spatial patterns of forest overstory prior to timber harvest activities associated with Euro-American settlement can serve as helpful forest reference conditions. Yet, such historical reference patterns can be challenging to quantify due to a lack of appropriate spatially explicit field data and the limitations of dendrochronological techniques. As a result, many studies have investigated contemporary spatial patterns of tree density and tree size, but few have successfully described changes in spatial patterns. In this dissertation, I first used overstory census data collected in 1934, over an area of ~4,000 ha, to characterize forest overstory reference conditions for Blacks Mountain Experimental Forest (BMEF), California, USA. Specifically, I quantified the 1934 spatial patterns of stand density index by seven vigor classes for Ponderosa pine/Jeffrey pine and white fir and incense-cedar. The forest reference conditions for BMEF consisted of clusters of open forests dominated by high vigor ponderosa pine trees associated with soil moisture and elevation gradients. Second, I used both 1934 and 2018 regeneration data from BMEF to compare the stocking between reference and contemporary conditions. The stocking of the reference conditions was dominated by ponderosa pine seedlings and was associated with available water capacity at 150 cm and overstory basal area (m2 ha-1). Third, I utilized 2016 Light Detection and Ranging data and overstory ground truth data to compare among stands that received thinning treatments, prescribed burn treatments, and a combination of thinning with prescribed burn treatments. I then compared above-ground biomass (Mg ha-1) between contemporary (2016) and reference (1934 census data) conditions. Spatial patterns in above-ground biomass were more heterogeneous in reference conditions than in iv contemporary conditions. However, in contemporary conditions, prescribed burns combined with thinning were better able to emulate the spatial patterns of reference conditions than thinning treatments alone. The findings of this dissertation can be utilized by managers: a) as a reference point against which they can evaluate the success of restoration treatments, b) to understand the relationships between regeneration and environmental factors, which shape the overall stand structure, and c) to identify contemporary stands at-risk of insect, disease, and fire that may benefit from restoration treatments. v Lay Summary Forest managers aiming to restore resilience of contemporary forests to better reflect what existed prior to settlement timber harvest—aka reference conditions—want to utilize information on the spatial patterns of forests. In this dissertation, I quantified the spatial patterns of reference conditions using tree vigor and density utilizing 1934 data from Blacks Mountain Experimental Forest (BMEF), USA, California. I also examined the differences in regeneration between reference (1934) and contemporary forests (2018) utilizing field data collected at BMEF. Finally, I quantified and compared the differences in spatial patterns of above-ground biomass between reference and contemporary BMEF utilizing 1934 data and 2016 light detection and ranging and inventory data. This dissertation provides unique information on overstory spatial patterns of tree vigor and regeneration in reference conditions, and highlights differences in spatial patterns of overstory conditions between reference and contemporary forests. vi Preface The dissertation's research questions and objectives were originally conceived from discussions between myself and my supervisory committee: Drs. Bianca Eskelson, Martin Ritchie, Nicholas Coops, and Sarah Gergel. All the data used for this dissertation were provided by the United States Department of Agriculture, Forest Service, Pacific Southwest Research Station under the supervision of Dr. Martin Ritchie. The research chapters of this dissertation (Chapters 3-5) were written as a series of manuscripts for publishing in peer-reviewed journals. For each research chapter, the objectives, data organization, statistical analysis and preparation of draft manuscripts, including tables and figures, were performed by me with contributions and editorial advice from Dr. Bianca Eskelson and Dr. Martin Ritchie. Versions of chapters 3 and 4 have been published as: Nepal, S.N., Eskelson, B.N.I, Ritchie, M.W. Gergel, S.E., 2023. Spatial patterns of vigor by stand density across species groups and its drivers in a pre-harvest ponderosa pine-dominated landscape in northern California. Forest Ecology and Management, 534 (101972). https://doi.org/10.1016/j.foreco.2023.120867 • Sushil Nepal identified the research problem, designed research, performed the statistical analysis and wrote the manuscript for publication. • Bianca Eskelson and Martin Ritchie conceived the idea, advised on the research approach, and provided data and editorial support. • Sarah Gergel provided editorial support. vii Nepal, S.N., Eskelson, B.N.I, Ritchie, M.W. 2020. Difference in regeneration conditions in Pinus ponderosa dominated forest in northern California, USA, over an 83 years period. Forests. 11 (581). https://doi.org/10.3390/f11050581 • Sushil Nepal identified the research problem, performed the statistical analysis and wrote the manuscript for publication. • Martin Ritchie designed the sampling strategy for the 2018 data collection. • Bianca Eskelson and Martin Ritchie conceived the idea, advised on the research approach, and provided data and editorial support. A version of Chapter 5 is being prepared for submission as peer-reviewed publication. viii Table of Contents Abstract ......................................................................................................................................... iii Lay Summary .................................................................................................................................v Preface ........................................................................................................................................... vi Table of Contents ....................................................................................................................... viii List of Tables ................................................................................................................................xv List of Figures ............................................................................................................................ xvii List of Abbreviations ............................................................................................................... xxiii Acknowledgements .................................................................................................................. xxvi Dedication ................................................................................................................................ xxvii Chapter 1: Introduction ................................................................................................................1 1.1 Background ..................................................................................................................... 1 1.1.1 Spatial patterns in contemporary dry forests: a risk to resilience? ............................. 1 1.1.2 Limitations in approaches for understanding spatial patterns of overstory between reference and contemporary conditions .................................................................................. 4 1.2 Dissertation overview and research goals ....................................................................... 5 Chapter 2: Study site and data ...................................................................................................10 2.1 Blacks Mountain Experimental Forest.......................................................................... 10 2.1.1 Land use and fire regime of Blacks Mountain Experimental Forest ........................ 10 2.1.2 Research natural areas (RNAs) ................................................................................. 11 2.1.3 Blacks Mountain Ecological Research Project ......................................................... 12 2.2 Data Representing Reference Conditions ..................................................................... 15 ix 2.2.1 Forest inventory methods commonly used in the 1930s in the western United States ………………………………………………………………………………………15 2.2.2 1934 census data ....................................................................................................... 16 2.2.3 1934 regeneration data .............................................................................................. 18 2.3 Contemporary LiDAR data ........................................................................................... 18 2.4 Auxiliary data................................................................................................................ 19 Chapter 3: Spatial patterns of vigor by stand density across species groups and its drivers in a ponderosa pine-dominated landscape of the 1930s in northern California ....................21 3.1 Introduction ................................................................................................................... 21 3.2 Data ............................................................................................................................... 24 3.2.1 Study area.................................................................................................................. 24 3.2.2 Overstory data from 1934 census survey .................................................................. 26 3.2.3 Auxiliary soil and topography variables ................................................................... 26 3.2.3.1 Soil variables ..................................................................................................... 27 3.2.3.2 Topography variables........................................................................................ 27 3.3 Analysis......................................................................................................................... 28 3.3.1 Principle component analysis to reduce data dimension .......................................... 28 3.3.2 Cluster analysis ......................................................................................................... 29 3.3.2.1 Interpretation of cluster categories based on SDI values .................................. 29 3.3.3 Association between clusters and auxiliary variables ............................................... 30 3.3.3.1 Classification and regression trees .................................................................... 30 3.3.3.2 Multinomial logistic regression ........................................................................ 30 3.4 Results ........................................................................................................................... 31 x 3.4.1 Principal components 1 and 2 explained the variance in the data ............................ 31 3.4.2 Cluster analysis ......................................................................................................... 32 3.4.2.1 Cluster size k = 4 was internally consistent, spatially structured, and lacking hierarchical arrangement ................................................................................................... 32 3.4.2.2 Spatially aggregated clusters of different stand density index and tree vigor .. 34 3.4.3 Elevation, surface-soil depth, and sub-soil depth were associated with clusters ...... 38 3.5 Discussion ..................................................................................................................... 41 3.5.1 Low stand density was associated with tree vigor .................................................... 41 3.5.2 Spatial patterns of vigor in BMEF indicate a resilient forest ecosystem .................. 43 3.5.3 Elevation and soil characteristics influenced tree vigor ........................................... 43 3.5.4 Management implications ......................................................................................... 46 3.6 Conclusions ................................................................................................................... 47 Chapter 4: Difference in Regeneration Conditions in Pinus ponderosa Dominated Forests in Northern California, USA, over an 83 Year Period .................................................................49 4.1 Introduction ................................................................................................................... 49 4.2 Materials and Methods .................................................................................................. 53 4.2.1 Study area.................................................................................................................. 53 4.2.2 Management history.................................................................................................. 53 4.2.3 Regeneration data...................................................................................................... 57 4.2.3.1 1934 regeneration data ...................................................................................... 57 4.2.3.2 Contemporary regeneration data ....................................................................... 57 4.2.4 Data compilation ....................................................................................................... 57 4.2.4.1 Stocking ............................................................................................................ 57 xi 4.2.4.2 Explanatory variables for 1934 conditions ....................................................... 58 4.3 Statistical analysis ......................................................................................................... 61 4.3.1 Stocking in reference conditions ............................................................................... 61 4.3.2 Effect of overstory, soil, and topography on stocking .............................................. 61 4.3.3 Stocking in contemporary forests under different management treatments ............. 63 4.3.4 Stocking in reference vs. contemporary forests ........................................................ 63 4.4 Results ........................................................................................................................... 63 4.4.1 High stocking for ponderosa pine species in reference conditions ........................... 63 4.4.2 Overstory and soil variables were associated with stocking ..................................... 68 4.4.3 Differences in stocking among management treatments in 2018 ............................. 69 4.4.4 Stocking in reference conditions higher than contemporary stocking ...................... 69 4.5 Discussion ..................................................................................................................... 72 4.5.1 Stocking in reference conditions was dominated by ponderosa pine at low elevations ………………………………………………………………………………………72 4.5.2 Stocking in reference conditions was dominated by height classes 2 and 3 ............. 73 4.5.3 Available water capacity at 150 cm and overstory basal area were drivers of stocking in reference conditions ........................................................................................... 75 4.5.4 Contemporary stocking was low in thinned stands and differed across treatments . 77 4.5.5 Higher stocking in reference conditions than in managed, contemporary stands ..... 79 4.6 Conclusions ................................................................................................................... 79 Chapter 5: Quantifying changes in spatial patterns of Pinus ponderosa dominated overstory above-ground biomass across 80 years using LiDAR and inventory data ............81 5.1 Introduction ................................................................................................................... 81 xii 5.2 Methods......................................................................................................................... 84 5.2.1 Study area.................................................................................................................. 84 5.2.2 Data ........................................................................................................................... 86 5.2.2.1 1934 census data ............................................................................................... 86 5.2.2.2 2015 LiDAR data acquisition and processing .................................................. 86 5.2.2.3 2016 field data................................................................................................... 86 5.2.2.4 Overstory above-ground biomass for 1934 and 2016 ....................................... 87 5.2.2.5 Biomass model to link LiDAR metrics to 1 ha plots for 2016 ......................... 87 5.3 Data analysis ................................................................................................................. 88 5.3.1 Comparing semivariogram components and Moran’s correlogram between 1934 and 2016…………………………………………………………………………………………88 5.3.2 Effect of management treatments on spatial patterns of above-ground biomass in the contemporary forests ............................................................................................................. 91 5.4 Results ........................................................................................................................... 91 5.4.1 Differences in spatial patterns between 1934 and 2016 ............................................ 91 5.4.1.1 Spatial autocorrelation in above-ground biomass exhibited at larger range in 2016………………………………………………………………………………………91 5.4.1.2 Moran’s I correlogram exhibited differences in small scale spatial patterns and periodicity between two years .......................................................................................... 92 5.4.2 2016 spatial patterns in treatment units .................................................................... 98 5.4.2.1 HID Burned halves exhibited spatial autocorrelation at short ranges ............... 98 5.4.2.2 HID Burned halves exhibited heterogeneous spatial patterns indicated by periodicity at shorter lag-distances ................................................................................... 99 xiii 5.5 Discussion ................................................................................................................... 103 5.5.1 The 1934 forest was spatially heterogeneous particularly at small-scales ............. 103 5.5.2 Larger and more homogenous spatial patterns in 2016 except when burned ......... 104 5.5.3 Prescribed burn enhanced spatial heterogeneity in both LOD and HID treatments 105 5.5.4 Management implications ....................................................................................... 106 5.6 Conclusion .................................................................................................................. 108 Chapter 6: Conclusion ...............................................................................................................109 6.1 Overall summary and contributions ............................................................................ 109 6.2 Limitations .................................................................................................................. 111 6.2.1 Reference conditions spatial pattern chapter limitations ........................................ 111 6.2.2 Regeneration chapter limitations ............................................................................ 112 6.2.3 Reference vs. contemporary spatial patterns chapter limitation ............................. 113 6.3 Future directions and recommendations ..................................................................... 114 6.3.1 Utilize fire data to improve the understanding of spatial patterns in reference conditions ............................................................................................................................ 114 6.3.2 Use of spatially explicit regeneration data and climate variables ........................... 114 6.3.3 Fine scale spatial patterns and use of environmental covariates ............................. 115 6.4 Resiliency comes from heterogeneity ......................................................................... 116 Bibliography ...............................................................................................................................117 Appendices ..................................................................................................................................150 Appendix A Additional tables for Chapter 2 and 3 ................................................................ 150 Appendix B LiDAR point cloud processing and model building process .............................. 157 Appendix C Additional tables for Chapter 5 .......................................................................... 159 xiv Appendix D Additional figures for Chapter 5 ........................................................................ 163 xv List of Tables Table 3.1 Overstory and auxiliary variables. .................................................................................. 28 Table 3.2 Summary statistics for each cluster category given as mean SDI and RD with standard deviation in parentheses by species and vigor group for cluster categories in cluster ................. 37 Table 3.3 Summary of Akaike Information criteria (AIC) value, kappa coefficient, percentage of correctly classified plots (% correct) for multivariable models across the auxiliary variables group. ............................................................................................................................................ 40 Table 4.1 Treatment, compartment name, management description, compartment area, and number of 0.0016 ha quadrats collected in 1934 and 2018 for each of the 12 compartments for which the non-parametric ANOVA and t-tests were conducted. .................................................................. 56 Table 4.2 Compartment-level response and explanatory variables. ................................................ 60 Table 4.3 Summary statistics of stocking of all species, seedling density of all species, available water capacity at 150 cm (AWC150), overstory basal area per hectare for all species, and overstory density for all species for the 98 compartments in 1934 that were used in the model development. The overstory basal area was calculated for all trees >8.9 cm in diameter at breast height............................................................................................................................................. 61 Table 4.4 Pseudo-R2, Akaike Information Criterion (AIC), and precision parameter for selected beta regression models. Estimated coefficients with p-value in parentheses are given for each explanatory variable. The final beta regression model is highlighted in black font. .................... 68 Table 5.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. .......... 94 xvi Table 5.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. ....................................................................................................................................... 98 Table 5.3 Summary of the number of peaks and troughs of Moran’s I correlogram and p-values for the Moran’ I at first lag Burned and Unburned halves of low structural diversity treatment (LOD) and high structural diversity treatment (HID). ............................................................................ 103 Table A.1 Summary of the applied combination of structural diversity treatments and prescribed burn in various units of BMERP. ................................................................................................ 150 Table A.2 Brief description of each Dunning’s vigor classes (Dunning 1928). ........................... 151 Table A.3 Summary of observed and broken-stick model eigenvalues along with variance explained by each of the principal components. ......................................................................... 152 Table A.4 Summary of gamma statistics (γ) and Moran’s I for cluster sizes k = 2 to 10 ............. 153 Table A.5 Stand density index (SDI) by seven vigor classes for two species group: ponderosa pine (PP) and cedar-fir (CF) with standard deviation in parentheses. ................................................ 154 Table C.1 LiDAR-derived and topography auxiliary variables were used for biomass modeling………………………………………………………………………………………..159 Table C.2 Summary of the best-selected models from each auxiliary variable group. Highlighted in bold is the final model selected for above-ground biomass (Mgha-1) prediction…………...160 Table C.3 Summary statistics of the observed and predicted above-ground biomass (Mgha-1) values for circular plots (n=153) and cells (n=12,647)…………………………………………161 Table C.4 Summary of RMSE and values for the spherical, exponential, and Gaussian models selected for Burned and Unburned halves of HID and LOD units……………………………………………………………..…………………………........162 xvii List of Figures Figure 1.1 Schematic diagram showing links between different chapters, analysis, and data used. LOD: low structural diversity, HID: high structural diversity, RNA: research natural area, LiDAR: light detection and ranging, SDI: stand density index, PCA: principal component analysis. C: control, B: burn, TC-B: combination thinning with prescribe burn, TC: combination thinning, TB-B: thin from below with prescribe burn, TB: thin from below. ................................ 9 Figure 2.1 Study area map: A) United States of America, B) California, C) Blacks Mountain Experimental Forest showing sections, compartments, and 1 ha plots. The research natural areas (RNA), and the Burned and Unburned halves of the low structural diversity (LOD) and high structural diversity (HID) treatments in the figure are part of Blacks Mountain Ecological Research Project (BMERP) within the study area. ....................................................................... 14 Figure 3.1 The figure shows A) the sections and 1 ha plots of Blacks Mountain Experimental Forest in B) northern California in C) the United States of America. .......................................... 25 Figure 3.2 Results from PCA biplot and cluster analysis for cluster sizes A) k = 2, B) k = 3, C) k = 4, and D) k = 5. SDI: stand density index, PP = ponderosa pine, and Jeffrey pine group, CF = incense-cedar and white fir group, number in the subscript to the seven Dunning vigor classes, where 1, 2, 3, and 4 are high and 4, 5, and 7 are low vigor classes. ............................................. 33 Figure 3.3 Cluster analysis results for cluster sizes A) k = 2, B) k = 3, C) k = 4, and D) k = 5 with E) elevation (m) gradient and F) sub-soil depth (cm) gradient for BMEF. Cluster size k= 4 was chosen as it was internally consistent, spatially structured and lack hierarchical arrangement. .. 36 Figure 3.4 Decision tree to find the auxiliary variables predicting each of the four clusters in cluster size k = 4. The root nodes represent the value of variables above or below which the split xviii occurs. The number in parentheses represents the number of plots used in the classification for a given leaf node. ............................................................................................................................. 39 Figure 4.1 Location of Blacks Mountain Experimental Forest (BMEF) in northern California with 100 compartments and six management treatments. Treatments were applied at the compartment level. C: control, B: prescribed fire (1999/1997, 2011), TC: combination thinning (1996), TC-B: combination thinning (1996) with prescribed fire (1997), TB: thinning from below (2012), TB-B: thinning from below (2012) with prescribed fire (2017). ............................................................. 55 Figure 4.2 Observed stocking proportion of natural regeneration in 2018 (black points) with model fit of 0.0016ha stocking (dashed line) and bootstrapped 95% confidence interval (grey) from Ritchie (2020). .............................................................................................................................. 65 Figure 4.3 Histograms of 1934 stocking proportion calculated at 0.0016 ha for 100 compartments for a) all species, b) Ponderosa pine, c) white fir, and d) incense-cedar, with mean (solid black line) and median (dashed line). The scale on y-axis are different for both panels of figure. ....... 66 Figure 4.4 Histograms of 1934 stocking proportion at 0.0016 ha for 100 compartments by four height classes and three species:a) ponderosa pine, height class 1, b) incense-cedar, height class 1, c) white fir, height class 1, d) ponderosa pine, height class 2, e) incense-cedar, height class 2, f) white fir, height class 2, g) ponderosa pine, height class 3, h) incense-cedar, height class 3, i) white fir, height class 3, j) ponderosa pine, height class 4, k) incense-cedar, height class 4, l) white fir, height class 4. Height Class 1: 0–0.31 m, Class 2: 0.31–0.91 m, Class 3: 0.91–1.83 m and Class 4: >1.83 m for seedlings less than 8.89 cm diameter at breast height (DBH), with mean (solid black line) and median (dashed line). The data used for ponderosa pine includes all 100 compartments while the data for incense-cedar and white fir were from the 62 compartments that contained regeneration of both species or either of them. ............................................................ 67 xix Figure 4.5 Relationship between overstory basal area per hectare and stocking proportion calculated at 0.0016 ha for 98 compartments in 1934 (a). Relationship between stocking proportion by species and corresponding overstory basal area of that species for (b) ponderosa pine, (c) incense-cedar, and (d) white fir. The highlighted grey area represents the range of incense-cedar and white fir overstory basal area compared to ponderosa pine. ........................... 70 Figure 4.6 Stocking proportion calculated at 0.0016 ha in 1934 (red) and 2018 (blue) for 12 compartments by (a) years, (b) species, and (c) height class for various management treatments. Height Class 1: 0–0.31 m, Class 2: 0.31–0.91 m, Class 3: 0.91–1.83 m and Class 4: >1.83 m for seedlings less than 8.89 cm diameter at breast height (dbh). IC: incense-cedar, PP: ponderosa pine, WF: white fir. C: control, B: prescribed fire, TC: combination thinning, TC-B: combination thinning with prescribed fire, TB: thinning from below, TB-B: thinning from below with prescribed fire. .............................................................................................................................. 71 Figure 5.1 Study area map: A) United States of America, B) California, C) Blacks Mountain Experimental Forest, and D) permanent grid system within Blacks Mountain Ecological Research project (BMERP)........................................................................................................... 85 Figure 5.2 Range values for above-ground biomass as indicated by selected variogram models between 1934 vs. 2016 within research natural area (RNA) units (top panel), low structural diversity (LOD) units (middle panel), and high structural diversity (HID) units (bottom panel) in Blacks Mountain Experimental Forest.......................................................................................... 95 Figure 5.3 Nugget values for above-ground biomass as indicated by selected variogram models between 1934 vs. 2016 within research natural area (RNA) units (top panel), low structural diversity (LOD) units (middle panel), and high structural diversity (HID) units (bottom panel) in Blacks Mountain Experimental Forest.......................................................................................... 96 xx Figure 5.4 Moran’s I correlogram for 1934 vs. 2016 within the four research natural areas (RNAs) in Blacks Mountain Experimental Forest. 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 (Jaquette et al. 2021). 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, while 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. ................................................................................................. 97 Figure 5.5 Differences in range for above-ground biomass as indicated by selected variogram models in between Burned and Unburned halves of low structural diversity (LOD) units (upper panel) and high structural diversity (HID) units (lower panel) in Blacks Mountain Experimental Forest........................................................................................................................................... 101 Figure 5.6 Differences in sill for above-ground biomass as indicated by selected variogram models between Burned and Unburned halves within high structural diversity (HID) units (upper panel) and low structural diversity (LOD) units (lower panel) in Blacks Mountain Experimental Forest...................................................................................................................................................... 102 Figure D.1 Histogram for above-ground biomass (Mgha-1) for 2016 ground data. The red line indicates mean, and blue indicates median. The red line indicates mean while blue line median...................................................................................................................................................... 163 Figure D.2 Predicted above-ground biomass (Mgha-1) using the final model and LiDAR derived metrics in 2016 at 32 m X 32 m pixel level. HID: high structural diversity treatment; LOD: low xxi structural diversity treatment; RNA: research natural area; BMEF: Blacks Mountain Experimental Forest. ................................................................................................................... 164 Figure D.3 Observed vs. predicted above-ground biomass (Mgha-1) values for the best model from each of the variable groups along with the final model. ............................................................. 165 Figure D.4 Predicted overstory above-ground biomass (AGB, Mgha-1) for 1 A) predicted overstory AGB at cell level, and B) observed overstory AGB for circular ground truthing plots. The blue line represents the median and red line represents the mean for AGB. ....................... 166 Figure D.5 Observed overstory above-ground biomass (AGB, Mgha-1) for circular ground truthing plots in 2016. ............................................................................................................................... 167 Figure D.6 Variogram models for research natural areas (RNA) in 1934 and 2016. ................... 168 Figure D.7 Variogram models for low structural diversity units (LOD) in 1934 and 2016. ........ 169 Figure D.8 Variogram models for high structural diversity units (HID) in 1934 and 2016. ........ 170 Figure D.9 Differences in Moran’s I correlogram between 1934 vs. 2016 within different units under study in Blacks Mountain Experimental Forest. Dotted lines represent the lower and upper limit of the no significant spatial autocorrelation using 95% confidence envelopes for 1934 (red) and 2016 (blue) from Monte Carlo simulations. Points connected with solid lines indicate Moran’s I at a given lag-distance. Points above the dotted lines in upper part of envelope indicate a positive Moran’s I that is significantly different from 0, while points below the dotted lines in lower part of the envelope show negative Moran’s I significantly different from 0. Points within the dotted envelope shows no significant Moran’s I. ................................................................. 171 Figure D.10 Variogram models for Burned and Unburned halves of low structural diversity units (LOD) in 2016............................................................................................................................. 172 xxii Figure D.11 Variogram models for Burned and Unburned halves of high structural diversity units (HID) in 2016. ............................................................................................................................. 173 Figure D.12 Differences in Moran’s I correlogram between 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 no significant spatial autocorrelation using 95% confidence envelopes for Burned (red) and Unburned (blue) from Monte Carlo simulations. Points connected with solid lines indicate Moran’s I at a given lag-distance (m). Points above the dotted lines in upper part of envelope indicate a positive Moran’s I that is significantly different from 0, while points below the dotted lines in lower part of the envelope show negative Moran’s I significantly different from 0. Points within the dotted envelope shows no significant Moran’s I...................................................................................................................................................... 174 xxiii List of Abbreviations γ—Gamma statistic ABG—Above-ground biomass (Mg ha-1) AIC—Akaike information criterion Aka—Also known as ANOVA—Analysis of variance AWC150—Available water capacity at 150 cm soil depth AWC50— Available water capacity at 50 cm soil depth B—Burning BA—Basal area per hectare BMEF—Blacks Mountain Experimental Forest BMERP—Blacks Mountain Ecological Research Project C—Control CART—Classification and regression tree CF—Incense-cedar and white fir DBH—Diameter at breast height DEM—Digital elevation model ELV—Elevation EUI—Ecological unit inventory GIS—Geographic information system HID—High structural diversity k—Number of clusters LiDAR—Light detection and ranging xxiv LOD—Low structural diversity Mgha-1—Mega grams per hectare n—Sample size p—p-value PC—Principal component PCA—Principal component analysis phi—Precision parameter PP— Ponderosa pine/Jeffrey pine RD—Relative density RMSE—Root mean squared error RNA—Research natural areas SAGA—System of automated geoscientific analysis SBD—Sub-soil depth SDI—Stand density index SFD—Surface-soil depth SL—Slope ST—Surface stoniness TB—Thinning from below TB-B—Thinning from below with prescribed fire TC—Combination thinning TC-B—Combination thinning with prescribed fire TPH—Trees per hectare TWI—Topographic wetness index xxv USA—United States of America USDA—Unites States Department of Agriculture xxvi Acknowledgements First and foremost, I would like to express my sincere gratitude to my research supervisor, Associate Prof. Bianca Eskelson, for her encouragement, guidance, support, and patience with me that made all this work possible. Many thanks to my thesis committee, Dr. Martin Ritchie, Prof. Dr. Nicholas Coops, and Prof. Dr. Sarah Gergel. Thank you all for the assistance in analysis, interpretation of results, and for the help along the way. The four of you motivated me at every step on this journey and this project would simply not have been possible without the assistance of each of you. I also have tremendous gratitude towards Prof. Dr. Alan H. Taylor for acting as the external examiner for this dissertation. There are tons of people who helped me finish my studies. I would like to thank USDA PSW-Redding’s field crew for collecting regeneration data, and Dr. Earl Alexander for providing valuable insights on the Ecological Unit Inventory map. Special thanks go to Ethan Hammett and Brian Wing for providing the GIS data layers and support on any related questions. I would also like to thank Dr. Carl Skinner for sharing his fire scar data for BMEF. The Four-Year Fellowship provided by UBC, and funding provided by the USDA Forest Service, Pacific Southwest Research Station (grant number: 15-IJ-11272139-016), were essential to support the projects. Thanks to Dr. Adam Polinko and Dr. Andrew Sanchez Meador for helping me with the R code at times. I would also like to extend my gratitude to Dr. Paul Hacker and Dr. Francisco Mauro for providing valuable insights during LiDAR data processing. Gratitude to my beloved and dearest wife Dorothy Poudel for all those inspiring thoughts during my low times, you will be in my heart forever. And my dearest brother Bibek Adhikari for being there with me through my thick and thin. Lastly, thanks to my family and friends in Canada and Nepal for their love and support. xxvii Dedication This dissertation is dedicated to my late uncle Sushant Bagale, my late grandfather Kashinath Nepal, grandfather Dilli Ram Bagale, my father Shiva Prasad Nepal, my mom Manju Nepal, father-in-law Bhanu Poudyal, mother-in-law Manju Poudel, and the people who have supported me throughout my Ph.D. Thank you all for your support and encouragement. 1 Chapter 1: Introduction 1.1 Background 1.1.1 Spatial patterns in contemporary dry forests: a risk to resilience? Globally, dry forests have undergone major changes in stand structure, which are impacting the resilience—capacity to recover from disturbances such as wildfire, insects and disease—of forests to disturbance (Seidl et al., 2016). More specifically, overstory spatial patterns—mosaic of individual trees, tree groups, and gaps at a small-scale—can involve complex changes in the horizontal arrangement of species, tree density as well as tree vigor (Franklin & Pelt, 2004). The changes in overstory spatial patterns have arguably diminished the resiliency of contemporary dry forests against insects, disease, and wildfire spread. Overstory spatial patterns found in forests least impacted by timber harvesting following Euro-American settlement (hereafter referred as “reference conditions”) have been identified as supporting greater resilience to biotic and abiotic stresses, as compared to contemporary forests (Hessburg et al., 2019). For example, a forest with low tree density may consist of trees with high vigor that are resistant to insect, disease, fire, and drought (Merschel et al., 2021). Young trees may have high vigor (Dunning, 1928; Keen, 1936) as they secret more resins than old trees, resulting in stronger defense mechanism against insect attack (Perrakis et al., 2011). Young trees also have potential of faster recovery after insect or disease attack (Dunning, 1928; Kolb et al., 2007). For example, in reference conditions low tree density indicated minimum competition for resources, which resulted in clusters of high vigor and medium-sized trees (< 27.9–58.4 cm of diameter at breast height, DBH; Dunning, 1928; Keen, 1936) with few large (i.e., >54 cm DBH) fire-resistant, old trees (Churchill et al., 2018; Stephens et al., 2015). Collectively, such clusters 2 form heterogenous overstory stand structures, which are an indicator of forest function and as such, tracking changes in such indicators over time is imperative to research and management. Context of dry forests of California In western North America, many dry forests have undergone major changes in stand structure following the Euro-American settlement, which has created major differences in overstory spatial patterns between reference and contemporary forests (Boyden et al., 2005). Especially, ponderosa pine (Pinus ponderosa Laws.) dominated forests (also referred to as “dry forests” in this study) in the western United States of America (USA) have been particularly impacted. One important component of spatial patterns that has changed in contemporary dry forests is above-ground biomass (AGB), especially within California (Chen et al., 2004). AGB has substantially increased due to an increase in tree density of shade-tolerant species such as white fir (Abies concolor (Grod. And Glend)) and incense-cedar (Calocedrus decurrens (Torr.)) as a result of fire exclusion since the early 1900s (Collins et al., 2017). The increase in AGB has altered the fuel dynamics of frequent-fire forests, which has resulted in forest overstory stand structures and spatial patterns that are more prone to stand-replacing wildfire (Hurteau et al., 2016). Quantifying the differences in spatial patterns of AGB in dry forests is important for fuel classification, fire spread prediction, and post-disturbance vegetation changes (Jia et al., 2006; Keane, 2001; Rollins et al., 2004). Overstory stand structure in dry forests of the western United States of America (USA) are the product of complex interactions between forest growth, disturbance regimes (Boyden et al., 2005), as well as underlying terrain. Within dry forests of California, overstory spatial patterns reflect the availability of growing resources, which influence tree vigor (Patton et al., 2021). In dry forests in reference conditions, groups of medium-diameter trees, gaps, and 3 occasional large diameter trees of fire-resistant species formed heterogeneous overstory stand structures (Churchill et al., 2013; Larson & Churchill, 2012). These structural characteristics and spatial patterns of dry forests in reference conditions were the results of stands being substantially altered due to fire exclusion and livestock grazing around the early 1900s (Cooper, 1960; Stoddard et al., 2006). Historically, prior to the arrival of Euro-settlers, who changed forest harvest and fire, overstory spatial patterns in dry forests were associated with frequent fires at a mean return interval of 12-14 years (Cooper, 1960; Kilgore, 1973). While natural fires influenced tree density, basal area, and understory vegetation, topography also influenced forest species composition in reference conditions (Show & Kotok, 1924). Soil characteristics such as texture and moisture varied based on elevation, slope, and aspect (Stephenson, 1998). Because of such drivers, pine species were more dominant in the lower montane zone on xeric sites, while fir species dominated the mid to upper montane zones and mesic sites within northern California (Show & Kotok, 1924). Overstory spatial patterns also provide insights about how environmental processes in reference conditions, such as regeneration, mortality, and competition have shaped the current stand structure (Youngblood et al., 2004, Fry & Stephens, 2010; Taylor, 2010). For example, stand density can either promote or impede tree regeneration by manipulating microsite characteristics (Bonnet et al., 2005; Ouzts et al., 2015). Stand composition affects the regeneration rates of different co-occurring species by affecting the availability of trees to produce and distribute viable seeds under post-disturbance environments (Bailey & Covington, 2002). 4 Analyses of differences in spatial patterns between reference conditions and contemporary forest conditions (e.g., of the last 30 years) can increase our understanding of stand structure development pathways (Fortin & Dale, 2005; Levin, 1992). Analyzing the differences in spatial patterns will also allow forest managers to understand the effectiveness and potential of current management interventions such as prescribed fire and thinning that could restore the spatial patterns of reference conditions (Foster et al., 2003; Larson & Churchill, 2012; Reynolds et al., 2013). However, few studies describe the spatial perspective of differences in forest overstory stand structure between contemporary and reference conditions (Lydersen et al., 2013). 1.1.2 Limitations in approaches for understanding spatial patterns of overstory between reference and contemporary conditions Several studies have documented considerable differences in stand structure and composition over the past century for the frequent-fire dry forests of California (Scholl & Taylor, 2010; Collins et al., 2017). However, these studies based their findings on reference stand structures and their spatial patterns in dry forests using dendro-chronological reconstruction data (Bonnicksen & Stone, 1982; Youngblood et al., 2004). Reconstruction data used to quantify spatial patterns in California are available at small scales of tree- and stand-level (e.g., Lydersen et al., 2013; Pawlikowski et al., 2019; Youngblood et al., 2004), and therefore fail to address the relationship of spatial patterns to ranges of the environmental factors at landscape-level scales (Churchill et al., 2013). The small scale at which spatial patterns and their relationship to environmental processes are studied indicated a substantial lack of landscape-level scale, spatially-explicit data for both contemporary and reference conditions (Lydersen et al., 2013). 5 Often, wherever data are available for reference conditions, the sampling method is unknown, which hinders appropriate selection of statistical analysis methods (Stephens, 2000). It is also unclear from reconstruction studies how overstory spatial patterns may have changed over time as dry forests have experienced a long fire free interval post Euro-American settlement (Lydersen et al., 2013; Taylor, 2010). Most information that is available for differences in overstory stand structure between reference and contemporary conditions is quantified using stand-level attributes such as average tree size, tree density, basal area, and species composition (Boyden et al., 2005; Stephens, 2000). There is a lack of spatially explicit information on overstory stand structure for reference conditions. Such lack of information hinders the development of spatial models that can serve as guidelines for restoration efforts to help structure future forest developmental pathways (Shan et al., 2009; Youngblood et al., 2004). 1.2 Dissertation overview and research goals This dissertation addressed three research objectives around reference and contemporary conditions of dry forests in northern California, USA (Figure 1.1): 1) describe spatial patterns of overstory structures in dry forests in reference conditions of the 1930s, 2) compare stocking in contemporary and reference conditions, and 3) compare above-ground biomass in contemporary and reference conditions. I first introduce the study area of Blacks Mountain Experimental Forest (BMEF) in northern California and provide an overview of the data that were used for the reference and contemporary conditions analyzed (Chapter 2). BMEF makes a unique study site, because 1934 field survey on overstory and regeneration are available. The forest conditions in 1934 were minimally impacted by timber harvest, and I refer to them as “reference conditions”. Extensive field survey data for reference conditions as they are available for BMEF are very rarely available for the dry forests found across the USA. In addition, contemporary regeneration 6 and overstory data as well as light detection and ranging (LiDAR) data are available for BMEF, which allowed for the unique comparison of the differences in overstory and regeneration data in time. The details on data available for BMEF and their description are presented in Chapter 2. In Chapter 3, I used overstory data collected in 1934 that consisted of a complete census of live trees within 3,830 plots of ~1 ha size to identify the spatial patterns in overstory stand structure at a landscape-scale. I described the reference conditions of overstory spatial patterns in terms of stand density and vigor. Then, I linked the identified spatial patterns to soil and topography variables. The specific research objectives in this chapter were to: 1. Describe the spatial patterns of reference conditions of dry forests based on stand density index by seven vigor classes and two conifer species groups, and 2. Identify the soil and topography variables associated with the spatial patterns of forests in reference conditions of BMEF. In Chapter 4, I used a sample of regeneration data from 1934, as well as regeneration surveys from 2018 that allowed for comparisons of regeneration indicated by stocking. Then, I linked the 1934 stocking to soil, topography, and fire history variables along with overstory metrics such as basal area per hectare using beta regression models. In this study, I answered the following four research questions: 1. What was the stocking of ponderosa pine, incense-cedar, and white fir for four different height classes in 1934 prior to active forest management across BMEF? 2. Were reference conditions stocking in BMEF associated with overstory stand structure, soil, and topography? 3. Are there differences in stocking among contemporary (2018) forests that received several management treatments: control, prescribed fire, combination thinning, 7 combination thinning with prescribed fire, thinning from below, and thinning from below with prescribed fire? If so, which treatments differ? 4. How does reference condition stocking (1934) for the observed tree species differ from contemporary stocking (2018) for the 12 compartments that received management treatments prior to 2018? Chapters 3 and 4 complement each other by describing overstory and understory vegetation structures in references conditions and linking them to soil and topography information (Figure 1.1). Chapter 5 focused on comparing differences in spatial patterns in above-ground biomass (Mgha-1) of forests in reference (1934 census data) and contemporary (2016 field and 2015 LiDAR data) conditions. The specific research objectives in this chapter were to: 1. Use census data from 1934 and 2015 LiDAR data to compare the differences in spatial patterns of overstory above-ground biomass (Mgha-1) between reference and contemporary conditions, and 2. Compare the differences in spatial patterns of overstory AGB among various management treatments in contemporary conditions using LiDAR data. Chapters 4 and 5 complemented each other in terms of comparing reference and contemporary forests where the contemporary stands have received some of the same management treatments (Figure 1.1). Lastly, Chapter 6 provides a summary of the results with an overall conclusion. In this chapter, I highlight the contributions of this dissertation to forest research that investigates reference conditions. Chapter 6 also highlights how to utilize the information on overstory spatial patterns and the differences in overstory spatial patterns between reference and contemporary 8 conditions to inform management decision for contemporary forests. Further, Chapter 6 describes the limitations of this dissertation and future research directions. 9 Figure 1.1 Schematic diagram showing links between different chapters, analysis, and data used. LOD: low structural diversity, HID: high structural diversity, RNA: research natural area, LiDAR: light detection and ranging, SDI: stand density index, PCA: principal component analysis. C: control, B: burn, TC-B: combination thinning with prescribe burn, TC: combination thinning, TB-B: thin from below with prescribe burn, TB: thin from below. 10 Chapter 2: Study site and data 2.1 Blacks Mountain Experimental Forest The studies presented were conducted in Blacks Mountain Experimental Forest (BMEF, 40°40’ N, 121°10’ W), located in northern California, USA (Figure 2.1). BMEF was established in 1934 and is managed by the United States Department of Agriculture Forest Service, Pacific Southwest Research Station. BMEF covers an area of 4,000 ha within an elevation range between 1700 and 2100 m (Ritchie, 2016). BMEF is dominated by ponderosa pine in lower elevations. Higher elevations include white fir (Abies concolor (Gord. and Glend.) Lindl.) and incense-cedar (Calocedrus decurrens (Torr.) Florin) along with ponderosa pine (Wing et al., 2012). In lower elevations, Jeffrey pine (Pinus jeffreyi (Grev. and Balf)) is found within some stands (Oliver, 2000). Western juniper (Juniperus occidentalis Hook.) is sparsely distributed as individuals or in small groves in lower elevations of the forest (Oliver, 2000). The lower elevations are dominated by soil of subgroup Typic Argixerolls with 1–3 m deep basaltic bedrock, and mesic soil temperature regimes, while higher elevations are dominated by Andic Argixerolls with frigid soil temperature regimes (Zhang et al., 2008). Within BMEF, summers are warm and dry, while winters are cold and wet with average daily temperature ranging from -9 to 29.8 ˚C (Zhang et al., 2008). Annual precipitation is approximately 460 mm and falls primarily as snow from October through May (Zhang et al., 2008). 2.1.1 Land use and fire regime of Blacks Mountain Experimental Forest Before Euro-American settlers arrived (prior to the 1840s), the land of BMEF was used by the Native American Atsugewi tribe for mule deer hunting (Gordon, 1981). Later, between 1840-1850, the Lassen Trail that passes right through the middle of BMEF was mostly used by Euro-11 American settlers as a shortcut from the Oregon Trail to Sacramento Valley (Gordon, 1981). BMEF also experienced extensive summer cattle grazing that started in the 1850s. By the 1890s, due to reduced carrying capacity of ranges for cattle, livestock switched from cattle to sheep, which continued up until 1918, which was the last year that a grazing permit was issued (Gordon, 1981). The heavy sheep grazing destroyed the low vegetation, brush, and ground-hugging broadleaved herbs within BMEF and its surroundings, which reduced surface fuels for low-severity fire within BMEF (Gordon, 1981). Frequent low-severity fire, with a mean fire return interval between 5 to 17 years, once occurred in BMEF and has been effectively excluded from the experimental forest since the late 1800s (unpublished data from Carl Skinner). Many of the fires in BMEF are likely to have ignited due to lightning rather than human ignition (Gordon, 1981). The disruption of the fire regime due to widespread sheep grazing was later followed by aggressive wildfire suppression by the Forest Service in the early 1900s (Gordon, 1981; Ritchie, 2016). Until the Cone Fire of 2002, the last fire in BMEF found in the fire record was in 1899 (unpublished data from Carl Skinner). In the early 1900s, BMEF had not yet been impacted by timber harvesting due to lack of access. Harvesting and other active management activities in BMEF started in the late 1930s (Ritchie, 2016) after a transportation system had been established in about 1934. 2.1.2 Research natural areas (RNAs) In 1934, BMEF was divided into 100 compartments (Figure 2.1C). Of the 100 compartments, four were set aside as research natural areas (RNAs) during the establishment of the experimental forest in 1934 (Figure 2.1C). RNAs represent the remnants of the earlier structure of undisturbed and late successional forests, which until 1997 experienced fire exclusion (Camann et al., 2008; Oliver, 2000). However, after 1997 prescribed burns occurred as part of 12 management activities in two of the four RNAs, creating the following opportunities for comparison: • RNA A and D: No treatment, no management (Camann et al., 2008; Oliver, 2000); • RNA B: Prescribed fire burned with high secondary mortality 1999 (Camann et al., 2008); • RNA C: Two prescribed fires (1997, 2011) burned with low severity throughout the area (Camann et al., 2008). 2.1.3 Blacks Mountain Ecological Research Project The Blacks Mountain Ecological Research Project (BMERP) was initiated in 1991 as an interdisciplinary long-term ecological research project at BMEF (Oliver, 2000) with the goals 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 (Oliver, 2000). Ten research units of the BMERP were subjected to two different types of treatments (Figure 2.1C). The first treatment consisted of two stand structures: low structural diversity (LOD), and high structural diversity (HID), while four research natural areas (RNAs) units served as control (Figure 2.1C, Oliver, 2000). Five research units were each randomly assigned to both the LOD and HID treatments ranging in size from 77 to 144 ha (Figure 2.1C). Each research unit was then split in half with one randomly assigned half receiving prescribed burn treatments (hereafter “Burned”), while the other half did not receive the prescribed burn treatment (hereafter “Unburned”) (Figure 2.1C). Hereafter, I will refer to each Burned and Unburned half of the LOD, HID, and RNA units as stands (Appendix B: Table A.1). 13 Collectively, a total of 24 stands were used in this research: four RNAs (as described above), ten LOD stands, and ten HID stands (Figure 2.1C). LOD stands were thinned to a uniformly spaced density of ~40 trees ha−1 retaining a single canopy layer of intermediate trees—height ranging from 12 to 30 m and crown ratio > 50% (Oliver, 2000). In contrast, thinning in HID stands was conducted to retain all canopy layers that represented a feature of forest with multiple age classes and varying crown structures (Oliver, 2000). All large, old trees were maintained with one smaller tree retained within the larger tree's crown circumference (Oliver, 2000). As a part of the prescription, within the HID stands, small areas were left unharvested (Crotteau & Ritchie, 2014; Oliver, 2000). 14 Figure 2.1 Study area map: A) United States of America, B) California, C) Blacks Mountain Experimental Forest showing sections, compartments, and 1 ha plots. The research natural areas (RNA), and the Burned and Unburned halves of the low structural diversity (LOD) and high structural diversity (HID) treatments in the figure are part of Blacks Mountain Ecological Research Project (BMERP) within the study area. 15 2.2 Data Representing Reference Conditions I used data from different areas within BMEF for the data Chapters 3, 4, and 5. This section provides a general overview of inventory methods commonly used in the 1930s and a description of the methods used to collect the 1934 overstory and the 1934 regeneration data analyzed in this dissertation. 2.2.1 Forest inventory methods commonly used in the 1930s in the western United States Timber depletion following Euro-American settlement—especially addressed in the 1920 “Carper Report”—was a major concern during the early 1900s to the late 1930s. Hence, most of the forest inventory primarily focused on assessing timber volume (Frayer & Furnival, 1998; Lentz, 1931). Statisticians and researchers around the 1930s used systematic line-plot surveys to estimate timber volumes and realized that problems arose in estimating sampling errors, because of their non-random and systematic nature of sample plot, plot size, and sampling intensity (see e.g., Hasel, 1938; Hasel, 1942). In addition, mean timber volume of a certain area was estimated from a small percentage of timberland taken as sample that inevitably underestimated the true mean volume (Hasel, 1938). Hasel showed his interest in understanding the practical effects of random vs. systematic sampling, plot size, plot configuration, and sampling intensity on the estimation of sampling error (Hasesl, 1938; Hasel, 1942). With this objective in mind, Hasel conducted a complete census of the overstory trees > 8.9 cm diameter at breast height (DBH) in a virgin timber forest dominated by ponderosa and Jeffrey pines using contiguous 1 ha strips (2.5 by 10 chain or 50 x 200 m, hereafter referred to as “plots”) between the summer of 1933 and fall of 1934 (Hasel, 1934). This census (“referred to as 1934 census data”) conducted at BMEF facilitated the grouping of data into plots of various sizes at a variety of sampling intensities (Hasel, 1934). 16 Utilizing his complete census from BMEF, Hasel (1942) came up with a few recommendations such as: a) smaller sampling units had smaller sampling errors compared to larger ones, b) long narrow plots were more efficient than square plots of the same area, c) line-plots with two or more equally spaced plots of fixed length were as efficient as single-plot sampling units or contiguous strips, and d) dividing the area into blocks of equal size and sampling each block with the same intensity were efficient in accurately estimated sampling error. 2.2.2 1934 census data The 4,000 ha of BMEF were divided into sections (Figure 2.1C) to conduct a census survey between the summer of 1933 and fall of 1934 (Hasel, 1934). Within sections, a strip sampling method was employed within 1 ha plots (Hasel, 1934, Figure 3.1C). Within each plot, all live trees >8.9 cm (3.5 inch) were tallied by 5.08 cm (2.0 inch) wide diameter classes so each diameter class was labelled by even inches (4, 6, 8, 10…). The class width of the first diameter class is slightly narrower than the others with only 3.81 cm (2.0 inch), because the lower-class boundary was 8.9 cm (3.5 inch) rather than 7.63 cm (3.0 inch). Tree species were recorded as yellow pine (Jeffrey or ponderosa pine), white fir, and incense-cedar (Hasel, 1934). Yellow pine represented a grouping of Jeffrey pine (Pinus jeffreyi Grev. & Balf) and ponderosa pine, two closely related species. To simplify, I refer to this grouping as ponderosa pine throughout the dissertation. In addition, one of the seven Dunning’s vigor classes—represented by numbers 1 through 7—was assigned to the individual trees (Dunning, 1928; see Appendix A; Table A.2 for details). Vigor classes 1, 2, and 6 represent young trees (50-150 years in age), resistant to insects and diseases within dominant and co-dominant crown positions (Dunning, 1928). Vigor class 3 represents mature trees (150-300 years), resistant to insects and disease within a co-dominant crown position (Dunning, 1928). 17 Vigor classes 4, 5, and 7 represent mature trees (150-300 years), with greater susceptibility to insects and disease within the intermediate and the suppressed crown positions (Dunning, 1928). The 1934 census data were some of the first quantitative records available for the ponderosa pine/mixed-conifer forests of northern California. The 1934 census data also possibly recorded tree vigor information for the first time and tree vigor data opened up other areas of forestry research such as forest health and forest entomology beside just timber volume. In addition, within the areas where data prior to Euro-American settlement were absent, census data from timber cruising are some of the oldest data that have been used to understand the forest overstory stand structure of reference conditions in various regions of the USA (e.g., Hagmann et al., 2013; Scholl & Taylor, 2010). While using 1934 census data, I identified that some of the plots within BMEF fell along mountain ridges or bordered meadows, which resulted in the plot area being less than 1 ha in size. Any calculations that were used for these plots in Chapters 3, 4, and 5 were proportionately increased or reduced to the 1 ha basis. The plot locations within BMEF were defined by the coordinates of the center of the plot and were obtained from a plot map provided by the Pacific Southwest Research Station, Redding, California, USA. 18 2.2.3 1934 regeneration data I used 1934 regeneration data collected in 100 compartments (Figure 2.1C) within BMEF for Chapter 4. Regeneration data were systematically sampled across the entire BMEF in 1934. Along the length of each 1 ha plot, four 1 milacre quadrats were installed at every 50 m (164 feet) interval (PSW Research Station, 1934). Within each quadrat, seedling data were tallied by three species—ponderosa pine/Jeffrey pine, white fir, and incense-cedar. Data were also tallied by height classes of 0-0.3, 0.3-0.9, 0.9-1.82 and > 1.82 m for seedlings less than 8.89 cm (3.5 inches) diameter at breast height (DBH) (PSW Research Station, 1934). The data were aggregated to compartment-level (Figure 2.1), resulting in seedling density by species and height class for 100 compartments (PSW Research Station, 1934). Hasel’s 1934 regeneration data were one of the first detailed regeneration records that existed in the California region, which were collected using a well-designed sampling method. 2.3 Contemporary LiDAR data 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 (Wing et al., 2012). The aircraft was flown at an altitude of 900 m using an opposing-flight line with a side lap of at least 50% and the scanning angle of ±14˚ (Wing et al., 2012). The narrow beam setting of beam (~25 cm diameter on ground) resulted in a very low percentage of multiple returns and a very high percentage of single returns (Wing et al., 2012). The average return was 6.9 points per m2 with a standard deviation of 5.9 points per m2 (Wing et al., 2012). LiDAR data derived for ten of the twelve research units within the Blacks Mountain Ecological Research Project (BMERP, Figure 2.1C) within BMEF were used in Chapter 5. 19 2.4 Auxiliary data The soil and vegetation information for BMEF was available from an ecological unit inventory (EUI) map (Alexander, 1994). The EUI map was provided by the Pacific Southwest Research Station, Redding, California. The EUI map, developed for the study area in 1994, characterized the whole BMEF and identified different geological landforms, soils, and potential vegetation types (Alexander, 1994). Soil samples were taken within hand-dug “model pits” at least 6 feet in depth along a single transect running diagonally through BMEF from the southwestern section to the northeastern section of BMEF (Alexander, 1994). As the landforms and vegetation type that are described in the EUI map changed along the transect, new soil samples were taken by digging additional model pits (Alexander, 1994). In addition to model pits, soil samples were also analyzed from the road-cuts and borrow pits along the transect. Soil types were then described by USDA nomenclature (Soil Survey Staff, 1951). Using the method suggested by Alexander (1988), the available water capacity was estimated between 50 and 150 cm based on soil texture and water retention differences between 0.1 and 15 bars suction. For other areas of BMEF outside the transect, soil samples were taken randomly along each landform by classical auger method (D. Young, personal communication, June 21, 2023). The soil samples were then linked to one of the soil types identified within the model pits (D. Young, personal communication, June 21, 2023). The EUI map was used to derive the soil variables that were used in Chapters 4 and 5. The details on methods used to derive soil variables are discussed in the respective chapters. A digital elevation model (DEM, Wing et al., 2019) at 1 m resolution was also provided by the Pacific Southwest Research Station, Redding, USA. The DEM was used to derive various 20 topography variables used in Chapters 3, 4, and 5, and the methods used to derive these topography variables are discussed in the respective chapters. 21 Chapter 3: Spatial patterns of vigor by stand density across species groups and its drivers in a ponderosa pine-dominated landscape of the 1930s in northern California 3.1 Introduction Prior to Euro-American settlement, ponderosa pine (Pinus ponderosa Laws.) dominated forests (hereafter referred to as “dry forest”) had overstory stand structures, which facilitated resilience to fire, insect, and disease (Stephens et al., 2008). Selective timber harvesting in the early 1900s has since altered these original (hereafter referred to as “reference conditions”) stand conditions. Conventional descriptions of overstory stand structure in reference conditions are based on stand-level attributes such as average stand density (trees ha-1) and basal area (m2 ha-1) (Boyden et al., 2005). Largely missing from the study of reference condition forests is an understanding of tree vigor, which reflects the ability of trees to ward off insect and/or disease attacks (Dunning, 1928; Keen, 1936). Understanding reference condition forests, in terms of tree vigor and the spatial patterns of trees, where spatial patterns are defined as clusters of stands at various stand densities (Getzin et al., 2006), can be used by foresters to identify and prioritize contemporary stands in need of restoration (Patton et al., 2021). However, spatial patterns of trees and tree vigor in reference conditions are difficult to ascertain due to a paucity of spatially explicit information for pre-harvest forests across large landscapes (Lydersen et al., 2013). The reference condition ponderosa pine-dominated forests in northern California, United States of America (USA), consisted of open, park-like overstory stand structures with low stand density (50-130 trees ha-1) (Cooper, 1961; Norman & Taylor, 2003). In these forests, elevation and moisture gradients in conjunction with fire and fire frequency influenced the spatial patterns 22 of species and stand densities (Youngblood et al., 2004). Specifically, low to mid-elevation (<1700 m) and arid sites consisted of forests with low tree density (~36-88 trees ha-1) dominated by ponderosa pine and Jeffrey pine (Pinus jeffreyi Grev. & Balf.) (Collins et al., 2021; Hagmann et al., 2013). Low stand densities at lower elevation were maintained by high frequency fires with return intervals of 12-15 years (Norman & Taylor, 2005; Taylor, 2000). In contrast, white fir (Abies concolor (Grod. and Glend.) and incense-cedar (Calocedrus decurrens (Torr.) Florin) increase at higher elevation (1900-2200 m), mesic sites that experienced more infrequent fires at mean fire return interval > 20 years (Show & Kotok, 1924, Norman & Taylor, 2003), achieving abundance and tree densities >200 trees ha-1. In addition to elevation and moisture, the spatial patterns of stand density is also influenced by environmental processes such as regeneration, mortality, and competition (Stephens, 2000; Youngblood et al., 2004). Therefore, the spatial patterns of density provide important insights on patterns of tree esstablishment (Boyden et al., 2005), resource utilization, and growth at a given location (Olsen et al., 1995). Stand densities in ponderosa pine-dominated forests are often related to the availability of growing resources, which influence tree vigor (Patton et al., 2021). An open forest with low tree density has trees with high vigor, increasing its resistance to insect, disease, fire, and drought (Merschel et al., 2021). For example, in reference condition ponderosa pine-dominated forests, low tree density resulted in low competition for resources, which resulted in clusters of high vigor (Dunning, 1928; Keen, 1936) and medium-size trees (> 27.9–58.4 cm of diameter at breast height, DBH) with few large (i.e., >58.4 cm DBH) fire-resistant, old trees (Churchill et al., 2018; Stephens et al., 2015). Young trees with high vigor (Dunning, 1928; Keen, 1936) can produce more resins than old trees, providing a defense mechanism against insect attacks (Perrakis et al., 2011). Young trees also potentially recover faster after insect or disease attack (Dunning, 1928; 23 Kolb et al., 2007). In contrast, in ponderosa pine forests with high tree density, competition for resources can produce stands with low vigor, consisting of old and mature slowly-growing trees that are weak competitors for resources and easily killed by insects, disease, and drought (Dunning, 1928; Patton et al., 2021). However, older trees are more fire-resistant than younger trees (Dunning, 1928), provide wildlife habitat (Reynolds et al., 2013), and are living records of past natural disturbances such as fire and drought (Bailey & Idle, 2001). The information that currently exists about overstory stand density in reference conditions forests of California often comes from dendrochronological data (Everett, 2008; Taylor & Skinner, 2003) from only a few sites and trees (Stephens et al., 2003). Dendrochronological data often fail to capture the spatial patterns at landscape scales (Finley & Zhang, 2019; Stephens et al., 2003). Several studies capture landscape scale spatial patterns (e.g., Scholls & Taylor, 2010; Beaty & Taylor, 2008) using dendrochronological reconstructions. As reconstructions are derived from contemporary trees, dendrochronological approaches may not capture information from trees that existed during long ago pre-harvest periods (Zhang et al., 2019). Hence, a complete picture of the reference conditions stand density, species composition, and vigor cannot be depicted using dendrochronological data (Trotsiuk et al., 2018). A reliance on dendrochronology for studying forests of northern California is largely the result of a scarcity of reference conditions information regarding overstory stand structures at the landscape-scale (Stephens, 2000; Stephens et al., 2015). When such data exist, the exact location of study sites, poorly explained inventory methods, and limited spatial extent are major impediments to depicting landscape-scale spatial patterns of overstory stand structure (Stephens et al., 2015). Thus, there is a knowledge gap that requires field-based measurements of species-level tree vigor 24 information of pre-harvest forest condition to help managers understand how spatial patterns of overstory stand structures maintain ecosystem processes and function (Churchill et al., 2018). To fill this knowledge gap, I analyzed a landscape scale (~4,000 ha) census of live trees for the Blacks Mountain Experimental Forest (BMEF) in northern California, USA, containing information about tree vigor for four species: ponderosa and Jeffrey pine, incense-cedar, and white fir (Hasel, 1934). My goal was to characterize spatial patterns of tree vigor in the reference condition ponderosa pine-dominated forests and understand the drivers of such patterns. More specifically, we aimed to: 1) describe the spatial patterns of stand densities by seven vigor classes and two conifer species groups—ponderosa pine with Jeffrey pine and incense-cedar with white fir; and 2) identify the soil and topography variables associated with these reference conditions forest patterns. I expected to see heterogeneous spatial patterns in reference conditions that were associated with environmental soil and topography variables. Managers can use my results as baseline guide to identify, compare, and contrast stand densities in stands of different vigor classes and identify contemporary stands at-risk of insect, disease, and fire that may benefit from restoration treatments. Such baseline information will help managers design management prescriptions that achieve contemporary forest stand structures that more closely resemble the more open, reference condition forests that were more resilient to fire, insect, and disease. 3.2 Data 3.2.1 Study area I conducted this study in BMEF located in northeastern California, USA (Figure 3.1). A detailed description of the BMEF was given in Chapter 2 (see section 2.1). 25 Figure 3.1 The figure shows A) the sections and 1 ha plots of Blacks Mountain Experimental Forest in B) northern California in C) the United States of America. 26 3.2.2 Overstory data from 1934 census survey Using the 1934 overstory data (see section 2.2.2), I calculated the stand density index (SDI) by seven vigor classes for the two species groups, which resulted in 14 overstory variables for each 1 ha plot (Table 3.1). For the SDI calculations, I used the summation method proposed by Long & Daniel (1990) for uneven-aged and mixed-species forests. 𝑆𝐷𝐼𝑖𝑗 = ∑ 𝑇𝑃𝐻𝑖𝑗 (𝐷𝐵𝐻𝑖𝑗25)1.6 where 𝑆𝐷𝐼𝑖𝑗 represents the stand density index of trees in the ith vigor class and jth species group. 𝑇𝑃𝐻𝑖𝑗 and 𝐷𝐵𝐻𝑖𝑗 represent, respectively, the trees per hectare and the midpoint of the diameter class for the trees in ith vigor class and jth species group. The exponent of 1.6 was used in this study, because it is a commonly used exponent for ponderosa pine-dominated mixed-species forests in northern California, USA (Youngblood et al., 2004). Due to the diameter classes, SDI is bound to 9.1 cm (3.6 inches) at the lower end (Curtis, 2010). Some of the plots within BMEF fell along mountain ridges or bordered meadows, which resulted in the plot area being < 1 ha in size. In my exploratory analysis, I found that plots < 0.3 ha in size often had inflated tree density and basal area per hectare values due to the small plot size. Therefore, 244 plots < 0.3 ha in size were dropped, which resulted in a new total of 3830 plots used in my analysis. The plot locations within BMEF were defined by the coordinates of the center of the plot and were obtained from a plot map provided by the Pacific Southwest Research Station, Redding, California, USA. 3.2.3 Auxiliary soil and topography variables Auxiliary variables consisted of two groups—soil and topography (Table 3.1). 27 3.2.3.1 Soil variables A shape-file of an ecological unit inventory (EUI) map for BMEF (Alexander, 1994) had soil information stored as categorical variables (Table 3.1) for each EUI polygon. The plot map was overlaid with the EUI map and the intersect function in ArcMap 10.4.1 was used to extract the EUI polygons related to the plots (Reed et al., 2016). If a plot straddled multiple polygons from the EUI map, the midpoint of the soil variables from the EUI polygon that covered the maximum area within the plot was assigned to the plot. I did not use a weighted average because not every plot straddled multiple EUI polygons. Therefore, I did not have multiple soil variable values to work with, and hence the maximum area approach was the best choice. 3.2.3.2 Topography variables I extracted average elevation, slope, and aspect (Table 3.1) for each plot from a digital elevation model (DEM, Wing et al., 2019) at 1 m resolution provided by the Pacific Southwest Research Station, Redding, USA, using the zonal statistics function in ArcMap 10.4.4. (Wilson, 2012). Slope was converted to slope proportion. The topographic wetness index (TWI) was calculated based on slope, using the equations provided by Beven & Kirkby (1979). TWI indicates the effect of local topography on runoff flow direction and accumulation, which gives the moisture index of a given area (Beven & Kirkby, 1979). TWI is calculated by evaluating the flow direction, flow accumulation, and slope of an area using a DEM as follows: 𝐼 = 𝑙𝑛𝛼𝑡𝑎𝑛𝛽 where 𝐼 is the index value, 𝑙𝑛 is the natural logarithm, 𝛼 is the local upslope contributing area, and tan𝛽 is the local slope calculated using a 1 m DEM. 28 Table 3.1 Overstory and auxiliary variables. Variables Source Overstory variables SDIPP#: Stand density index (SDI) by seven Dunning vigor classes (# = 1, …, 7) for ponderosa and Jeffery pine (PP) 1934 Census Survey SDICF#: SDI by seven Dunning vigor classes (# = 1, …, 7) for incense-cedar and white fir (CF) 1934 Census Survey Auxillary variables Soil variables Available water capacity at 50 cm depth (AWC50, cm) Ecological Unit Inventory (EUI) Map (1991) Available water capacity at 150 cm depth (AWC150, cm) EUI Map (1991) Sub-soil depth (SBD, cm) EUI Map (1991) Surface-soil depth (SFD, cm) EUI Map (1991) Surface stoniness (ST, %) EUI Map (1991) Topography variables Aspect (degree) Digital Elevation Model (DEM at 1 m resolution) Slope proportion (SL) DEM at 1 m resolution Elevation (ELV, m) DEM at 1m resolution Topographic wetness index (TWI) DEM at 1 m resolution # = 1 to 7, corresponding to Dunning’s vigor classes (Dunning, 1928) 3.3 Analysis 3.3.1 Principle component analysis to reduce data dimension Principal component analysis (PCA; Hotelling, 1933) was used to reduce the dimension of our 14 overstory variables—SDI by two species groups—that were calculated for each of the 3830 29 plots. PCA reduced the 14 overstory variables into non-correlated principal components (PCs) using the “vegan” package (Oksanen et al., 2022) in R version 3.2.6 (R Core Team, 2020). PCs were retained if the observed eigenvalues (Kaiser, 1960) from the PCA exceeded the eigenvalues generated by the broken-stick distribution (Jackson, 1993; MacArthur, 1957). The broken-stick distribution assumes that the total variance in data—sum of eigen values—can be broken into pieces, where the length of each piece represents the variance explained by each principal component (Jackson, 1993). Therefore, if an observed eigenvalue from a PC exceeds the eigenvalue from the broken-stick distribution then all the PCs smaller than the observed eigenvalues are retained (Jackson, 1993). 3.3.2 Cluster analysis The retained PCs for the 3830 plots were used to perform k-means clustering (Jain & Dubes, 1988) using the function “k-means” in R version 3.2.6 (R Core Team, 2020). I set the cluster size from k = 2 to 10. I followed methods described in Plant (2012) and Roel & Plant (2004) to select a cluster size k, which I determined to be k=4 (see details in Appendix A.6). Based on cluster size k=4, I interpreted clusters based on their SDI values as described below: 3.3.2.1 Interpretation of cluster categories based on SDI values In this study, for the selected cluster size k, SDI values of clusters were compared and interpreted in relation to the maximum SDI (hereafter referred to as “SDImax”) of 1000 trees ha-1 for the Californian dry forests (Oliver, 1995). My study area also consisted of white fir and incense-cedar as minor species components, however, BMEF was ponderosa pine-dominated, hence I used 1000 as the max SDI. I converted the SDI of each cluster into relative density (RD) for easier interpretation. RD is calculated by dividing the given SDI by SDImax (Woodall et al., 2005). A cluster with a mean RD ≤ 0.25 was interpreted as open forest, and a cluster with a 30 mean RD between > 0.25 and < 0.35 was interpreted as a dense forest which is at the stage of tree competition initiation (Long, 1985; Long & Shaw, 2005). A cluster with RD ≥ 0.35 was an indication of a lower limit of full site occupancy and a dense forest with trees competing for resources (Long, 1985). 3.3.3 Association between clusters and auxiliary variables I used classification and regression tree analysis (CART, Breiman et al., 1984) and multinomial logistic regression (Hosmer & Lemeshow, 1989) to associate clusters with auxiliary variables using the “rpart” and “nnet” packages, respectively, in R version 3.2.6. (R Core Team, 2020). 3.3.3.1 Classification and regression trees CART analysis results in a decision tree with root nodes and leaf nodes that respectively represent the input auxiliary variables (Table 3.1) and cluster (Breiman et al., 1984; Roel & Plant, 2004). I used a complexity parameter of 0.01 to prune the tree and find the optimal tree size (Breiman et al., 1984) using the inbuilt ‘cp’ function in the ‘rpart’ package of R version 3.2.6. (R Core Team, 2020). 3.3.3.2 Multinomial logistic regression Prior to the model-building exercise, aspect and slope were incorporated into the model using the transformations presented by Stage (1976). The transformations allowed for the cosine shift on the optimum aspect and amplitude of the slope using the coefficients from the linear model that would allow for the interaction effect of aspect and slope (Stage, 1976). The model building process used in this study held ~30 % of the data for validation and ~70 % for training (Shoot et al., 2021). In my first step, using the training data, I fit multivariable multinomial logistic regression models separately for the soil and topography variable groups using a forward step-wise selection with cluster as the response (Hosmer & Lemeshow, 1989). From the two auxiliary 31 variable groups, variables that gave evidence of importance (p-value < 0.05), gave the lowest Akaike Information Criterion (AIC, Akaike, 1973), and a high percentage of correctly classified plots were retained (Hosmer & Lemeshow, 2000). I also calculated the kappa coefficient (Cohen, 1960) using the validation data to assess the degree of agreement between the predicted and the actual cluster for each model (kappa: -1 to 1; 1 perfect agreement, < 0 indicates no better agreement than expected by chance). Kappa was calculated using the following formula provided by (Kuhn, 2008): 𝑘𝑎𝑝𝑝𝑎 = 𝑂𝑏𝑠𝑒𝑟𝑣𝑒𝑑 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑐𝑎𝑡𝑒𝑔𝑜𝑟𝑦 − 𝐸𝑥𝑝𝑒𝑐𝑡𝑒𝑑 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑐𝑎𝑡𝑒𝑔𝑜𝑟𝑦1 − 𝐸𝑥𝑝𝑒𝑐𝑡𝑒𝑑 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑐𝑎𝑡𝑒𝑔𝑜𝑟𝑦 In my second step, I combined the variables that were retained from the soil and topography groups in a single multivariable model (Hosmer & Lemeshow, 2000). Variables that were not statistically significant previously were considered for entry back into the combined multivariable models to check for further improvement of the models (Hosmer & Lemeshow, 2000; Ferster et al., 2016 ). Interactions between auxiliary variables were evaluated for the combined multivariable model (Hosmer & Lemeshow, 2000). The model with the lowest AIC value (Akaike, 1973), a high percentage of correctly classified plots, and a high value for the kappa coefficient was selected as my final model. 3.4 Results 3.4.1 Principal components 1 and 2 explained the variance in the data PC1 and PC2, which were retained, explained 41 % of the variance in our data (Fig. 2) and had higher eigenvalues compared to eigenvalues from the broken-stick distribution (Appendix A: Table A.3). PC1 captured most information on the SDI and vigor for the PP group. PC1 was strongly correlated with the SDI values of the PP group in high vigor classes 3 and 6 along with 32 low vigor classes 4, 5, and 7 (Figure 3.2). PC2 captured most of the information on the SDI and vigor for the CF group (Figure 3.2). PC2 was strongly correlated with the SDI values of the CF group in high vigor classes 1 and 3 along with low vigor classes 4, 5, and 7 (Figure 3.2). 3.4.2 Cluster analysis 3.4.2.1 Cluster size k = 4 was internally consistent, spatially structured, and lacking hierarchical arrangement The gamma statistic (γ) was highest for cluster size k = 5 followed by cluster sizes k = 2 and k = 4 which indicated that clusters partitioned in each of the 10 replications for any given cluster size k were identical and internally consistent (Appendix A: Table A.4). Cluster size k = 2 was most spatially aggregated with plots of similar SDI, vigor, and species group next to each other followed by cluster sizes k = 3 and k = 4 as indicated by Moran’s I values (Appendix A: Table A.4). Cluster size k = 2 had higher γ and Moran’s I values than cluster size k = 4, however, two clusters could only capture the information on SDI by the two species groups (Figure 3.2A). In contrast, four clusters in cluster size k = 4 captured the information on SDI by vigor across the two species group (Figure 3.2 A vs. C). As the cluster size increased to k ≥ 5, both γ and Moran’s I decreased (Appendix A: Table A.4). In addition, the new cluster in cluster size k ≥ 5 was hierarchically arranged, being a subset of cluster 3 from cluster size k = 4. Therefore, I selected cluster size k = 4 for further CART and multinomial logistic regression analyses. 33 Figure 3.2 Results from PCA biplot and cluster analysis for cluster sizes A) k = 2, B) k = 3, C) k = 4, and D) k = 5. SDI: stand density index, PP = ponderosa pine, and Jeffrey pine group, CF = incense-cedar and white fir group, number in the subscript to the seven Dunning vigor classes, where 1, 2, 3, and 4 are high and 4, 5, and 7 are low vigor classes. 34 3.4.2.2 Spatially aggregated clusters of different stand density index and tree vigor Within cluster size k = 4, there were 4 distinctive clusters: ponderosa pine-high vigor (PP-High), ponderosa pine-low vigor (PP-Low), mixed species-high vigor (Mixed-High), and mixed species-low vigor (Mixed-Low) (Figure 3.2C and 3.3C). The PP-High cluster contained most of the plots (47.5%), followed by Mixed-Low (25.6%) (Figure 3.3C, Table 4). PP-High—The PP-High cluster had low RD (≤0.25) and mostly consisted of the PP species group in the high vigor classes 1, 2 and 6 with a substantial amount of vigor class 3 trees (Table 3.2, Appendix A: Table A.5). Within the PP-High cluster, 61% of plots had RD ≤ 0.25, while 6% of plots had RD ≥ 0.35 (Table 3.2). The PP-High cluster dominated the central and southwestern sections of BMEF that were in low elevation with greater sub-soil depth (Fig. 3C, E, and F). The presence of incense-cedar and white fir was low (Table 3.2, Appendix A: Table A.5). This cluster represented an open, high vigor forest at low elevation (Table 3.2). PP-Low—Within the PP-High cluster, mostly in the southern-central section of BMEF, were small, spatially aggregated patches of the PP-Low cluster (Figure 3.3 C, E, and F). The average RD of the PP-Low cluster was high compared to the PP-High cluster and mostly consisted of the 3.2PP species group in the low vigor classes 4, 5, and 7 (Table 3.2, Appendix A: Table A.5). Within, the PP-Low cluster, 49.7% of plots had RD ≥ 0.35, while only 4.2% of plots had RD ≤ 0.25 (Table 3.2). The presence of incense-cedar and white fir was low (Table 4, Appendix A: Table A.5). Since the majority of forest was > 0.25 RD, this cluster represented a dense, low vigor forest at low elevation compared to the PP-High cluster (Table 3.2). Mixed-High—The average RD of the Mixed-High cluster was comparable with the PP-Low (Table 3.2). However, the Mixed-High cluster consisted of trees from the two species groups while the PP-Low had mostly ponderosa pine species group (Table 3.2). Within the 35 Mixed-High cluster ~ 6% of plots had RD ≤ 0.25, while ~61% of plots had RD ≥ 0.35. The Mixed-High cluster dominated the northeastern, northwestern, and southeastern sections of BMEF that were high elevation areas with shallow sub-soil depth (Figure 3.3C, E, and F). High vigor trees in Dunning classes 1, 2, 3, and 6 in the CF group were mostly present in the Mixed-High cluster (Table 3.2, Appendix A: Table A.5). Mixed-Low—Within the Mixed-High cluster, throughout the northeastern, northwestern, and southeastern sections of BMEF, were small, spatially aggregated patches of the Mixed-Low cluster (Figure 3.3C, E, and F), which mostly included trees in the CF group, that were low vigor trees in Dunning classes 4, 5, and 7 (Table 3.2). Within the Mixed-Low clusters, ~17% of plots had RD ≤ 0.25, while ~32% of plots had RD ≥ 0.35 (Table 3.2). The trees in the PP group in both Mixed-High and Mixed-Low clusters were mostly low vigor trees in Dunning class 5 (Appendix A: Table A.5). 36 Figure 3.3 Cluster analysis results for cluster sizes A) k = 2, B) k = 3, C) k = 4, and D) k = 5 with E) elevation (m) gradient and F) sub-soil depth (cm) gradient for BMEF. Cluster size k= 4 was chosen as it was internally consistent, spatially structured and lack hierarchical arrangement.37 Table 3.2 Summary statistics for each cluster category given as mean SDI and RD with standard deviation in parentheses by species and vigor group for cluster categories in cluster. Clusters (n) Ponderosa pine and Jeffrey pine (PP) Incense-cedar and white fir (CF) Cluster’s mean SDI min SDI max SDI RD RD < 0.25 (% of plots) RD > 0.35 (% of plots) High vigor Low vigor High vigor Low vigor Mixed-high (304) 93.9 (62.7) 112.8 (46.2) 102.2 (35.4) 45.6 (35.3) 354.5 (74.8) 172.4 634.5 0.36 6.2 51.0 PP-low (723) 119.6 (50.4) 191.5 (64.3) 16.3 (15.8) 28.3 (26.1) 356.0 (71.0) 100.8 667.5 0.36 4.2 49.7 PP-high (1821) 119.2 (57.2) 99.3 (55.4) 2.7 (7.7) 2.3 (7.6) 223.6 (84.5) 164.7 723.4 0.22 61.0 6.1 Mixed-low (982) 83.9 (41.7) 225.9 (53.8) 43.8 (19.8) 47.0 (33.6) 316.8 (72.2) 61.8 643.6 0.32 16.8 32.1 Total: 3830 SDI=stand density index; PP = ponderosa pine group; CF = incense-cedar and white fir group; RD: Relative density of a given cluster category, min: minimum SDI value with in the given cluster category, max: maximum SDI value within the given cluster category; n: number of plots38 3.4.3 Elevation, surface-soil depth, and sub-soil depth were associated with clusters Based on both CART analysis and multinomial logistic regression, I found that elevation, surface-soil depth, and sub-soil depth were the important variables that were related to clusters in cluster size k = 4 (Figure 3.4, Table 3.4). The CART analysis also suggested that elevation, sub-soil depth, and surface-soil depth were related to clusters (Figure 3.4). For example, the PP-High cluster aggregated at lower elevation (<1806 m) at high surface-soil depth (Figure 3.4, node 1), while the Mixed-Low cluster aggregated at high elevation (>1806 m) at low sub-soil depth (Figure 3.4, node 6). In addition to the CART results, my final multinomial model suggested that surface stoniness, available water capacity at 150 cm depth, and its interaction with sub-soil depth were also related to clusters (Table 3.4, Model 15). My final model resulted in the highest classification accuracy, the highest kappa coefficient, and the lowest AIC value (Table 3.4, Model 15). My final model showed an interaction of sub-soil depth and available water capacity (AWC) at 150 cm, which was consistent with results from the CART analysis. Both CART results and interaction term in my final model suggested that at high elevations, areas where the sub-soil depth was high, the AWC at 150 m decreased resulting in low tree vigor clusters (Figure 3.4, Table 3.4). Including the slope-aspect transformation and the topographic wetness index in my models did not increase model performance substantially (Table 3.4, Model 14 vs. Models 17-18). 39 Figure 3.4 Decision tree to find the auxiliary variables predicting each of the four clusters in cluster size k = 4. The root nodes represent the value of variables above or below which the split occurs. The number in parentheses represents the number of plots used in the classification for a given leaf node. 40 Table 3.3 Summary of Akaike Information criteria (AIC) value, kappa coefficient, percentage of correctly classified plots (% correct) for multivariable models across the auxiliary variables group. Highlighted bold is the multinomial logistic regression model selected in each model building step, ELV: elevation (m), ST = surface stoniness, SBD = sub-soil depth (cm), SFD = surface-soil depth (cm), AWC150 = available water capacity at 150 cm soil depth (cm), TWI = topographic wetness index, SL = proportion of slope, cos(aspect) = cosine of aspect, sin(aspect) = sine of aspect, {SL+cos(aspect)+sin(aspect)} = Stage (1976) transformation. Auxiliary variables Model AIC % correct Kappa Multivariable models for soil ST+SBD 1 7962 57.3 0.31 AWC150 +SBD 2 7986 57.6 0.32 ST+SFD+SBD 3 7900 56.1 0.29 ST+SBD++SFD +AWC150 4 7850 62.3 0.41 ST+SBD+AWC150 5 7909 62.3 0.41 Multivariable models for topography ELV+TWI 6 7517 59.3 0.32 ELV+{SL+cos(aspect)+sin(aspect)} 7 7503 59.4 0.33 ELV+ TWI+{SL+cos(aspect)+sin(aspect)} 8 7497 59.3 0.32 Combined models for soil and topography ELV+{SL+cos(aspect)+sin(aspect)}+ST+SBD+TWI+AWC150 9 7148 61.9 0.38 ELV+{SL+cos(aspect)+sin(aspect)}+ST+SBD+TWI 10 7188 60.4 0.36 ELV+{SL+cos(aspect)+sin(aspect)}+ST+SBD+AWC150 11 7150 62.1 0.39 ELV+ST+SBD+TWI+AWC150 12 7152 61.6 0.38 ELV+ST+SBD+AWC150 13 7176 62.0 0.38 Combined models with interaction ELV+ ST+TWI+SBD+AWC150 +(SBD*AWC150) 14 6921 65.5 0.44 ELV+ST+SBD+AWC150 +(SBD*AWC150) 15 6752 66.3 0.46 ELV+ST+SBD+AWC150 +(SBD*AWC150)+(ELV*SBD) 16 6933 64.5 0.43 ELV+{SL+cos(aspect)+sin(aspect)}+ST+SBD+AWC150 +(SBD*AWC150) 17 6922 64.7 0.43 ELV+{SL+cos(aspect)+sin(aspect)}+ST+SBD+AWC150 +(SBD*AWC150)+(ELV*SBD) 18 6903 64.5 0.43 41 3.5 Discussion Before timber harvest activity in 1934, within the Blacks Mountain Experimental Forest (BMEF) nearly half of the plots consisted of low density stands, with distinct stand density patterns by tree vigor and species groups. The high vigor ponderosa pine stands (PP-High) were aggregated at lower elevations of the central and southwestern sections of BMEF and occurred in lower densities. In contrast, low vigor PP stands, and incense-cedar and white fir groups (CF-Low) were aggregated at high elevations in the southeastern, northeastern, and northwestern sections of BMEF, and occurred at higher densities. The spatial patterns of stand densities based on tree vigor and species groupings showed strong associations with elevation, soil moisture content, and sub-soil depth throughout the BMEF. 3.5.1 Low stand density was associated with tree vigor The low stand density (RD ≤ 0.25) for the PP-High cluster, found in the central and southwestern section of BMEF, supports previous assertions that pre-harvest ponderosa pine forests of California were open, park-like structures (Collins et al., 2021; North et al., 2022). Low-density stands are maintained by the repeated occurrence of low intensity and high frequency fire (10-12 years fire return interval; Norman & Taylor, 2003). Similar to Stephens et al. (2015), low density PP-High clusters within the southwestern section of BMEF could be the result of a frequent fire regime that existed before active fire suppression, which began in the early 1990s (Norman & Taylor, 2003; Norman & Taylor, 2005). The presence of trees in high vigor classes (classes 1, 2, 3, and 6) within the PP-High cluster in the southwestern section of BMEF supports the notion that at low stand density (RD ≤ 0.25), there is little tree competition within the forest (Oliver, 1995; Long, 1985), which likely enhances tree vigor (Cochran et al., 1994). The presence of high vigor trees in the PP-High 42 cluster in the southwestern and central section of BMEF is an indication of potentially resilient overstory stand structures, as a lack of suitable tree hosts for pine beetles to spread may support resilience (Oester et al., 2018; Olsen et al., 1995). Under bark beetle attack, the open overstory stand structure (RD ≤ 0.25) likely limits dispersal of insects among trees making a forest resilient against insects (Churchill et al., 2018). In addition, the southwestern section also contains mature ponderosa pine trees (vigor class 3) with thick bark (Dunning, 1928; Keen, 1936). Thick bark, which peels off upon burning, helps prevent the spread of fire into the canopy, creating a fire-resilient overstory (Fitzgerald, 2005; Norman & Taylor, 2003). Thick bark also acts as an insulator, preventing heat transfer into the cambium and thus minimizing scorching injuries (Fitzgerald, 2005; Norman & Taylor, 2003). High stand densities (RD ≥ 0.35) were found for the PP-Low, Mixed-High, and Mixed-Low clusters. Even though the mean relative density on these clusters seems similar, there were differences in percentage of the plots that were at RD ≤ 25 and the plots that were at RD ≥ 0.35. In addition, there were also differences in the species that dominate each cluster. For example, there was a substantial presence of low vigor trees (classes 4, 5, and 7) in the CF groups in the Mixed-High cluster compared to PP-Low cluster, which included low vigor PP groups. Low vigor trees in the CF group in Mixed-High and Mixed-Low clusters indicate competition-induced stress associated with high stand density (Larsson et al., 1983; Olsen et al., 1995). On the other hand, the density of the ponderosa pine species in vigor class 5 was relatively high in all clusters except PP-High, which indicated the existence of old and mature ponderosa pine trees within 1934 BMEF (Dunning, 1928). Their presence possibly was an indication of the fire-resistant nature of old and mature ponderosa pine trees (Fitzgerald, 2005). However, the low vigor in those old and mature ponderosa pine trees indicated that they are weak competitors of resources 43 and hence experience more competition stress compared to incense-cedar and white fir trees (Graham & Jain, 2005). 3.5.2 Spatial patterns of vigor in BMEF indicate a resilient forest ecosystem Unlike other studies (e.g., Churchill et al., 2017; Hessburg et al., 1994), the overall spatial pattern of the four types of forested clusters identified in 1934 BMEF represents a landscape with unique and heterogenous spatial structures in terms of stand density, tree vigor, and species composition as suggested by Merschel et al. (2021). In addition, where other studies (e.g., North et al., 2022) relied on SDI alone to make inferences about forest resiliency, one of the most unique attributes of my data set is the tree vigor information assessed during the 1934 inventory. The vigor information collected during the inventory allowed me to link the relative SDI with tree vigor to capture the spatial patterns of heterogeneous and resilient landscapes in accurate detail. This heterogeneous spatial pattern is indicative of forest stands that are resilient to insect and disease (e.g., Merschel et al., 2021). In general, heterogeneous spatial structures can reduce the spread of disease and insects when suitable tree hosts are only in patches dispersed throughout the landscape (Merschel et al., 2021). A patchy occurrence of suitable tree hosts lacks connectivity in both the roots and crowns, thus limiting the infestations to small localized tree groups or stands rather than the entire forest, hence, making the landscape resilient to insect and disease (Merschel et al., 2021). 3.5.3 Elevation and soil characteristics influenced tree vigor The influence of elevation on vegetation types within ponderosa pine-dominated landscapes in northern California has long been recognized (Show & Kotok, 1924), and the association of greater pine dominance at low elevations has also been observed (Fites-Kaufman et al., 2007, Lydersen & North, 2012, Skinner & Taylor, 2018). Within BMEF, the elevation gradient 44 influenced the species composition and stand density for the given spatial patterns of clusters. Similar to my results, within southern Sierra Nevada ponderosa pine forests (circa 1911), clumps of ponderosa pine forest occurred at elevations ≤1800 m, while mixed-species forest (of incense-cedar and white fir) were aggregated at elevations around ≥1850 m (Stephens et al., 2015). What is new however, is our understanding of tree vigor patterns, summarized next. Ponderosa pine dominance at lower elevations in pre-harvest forests of California has been observed elsewhere (Show & Kotok, 1924; Stephens et al., 2015). However, I further discerned relationships between soil characteristics and spatial patterns of vigor in pre-harvest forests. The amount of water at depth of 150 cm, sub-soil depth, and surface stoniness can vary with soil types (Hatfield et al., 2017) and influence the spatial patterns of tree species and vigor (McCullough & Wagner, 1987). For example, mollisols have high available water capacity at 150 cm, greater sub-soil depth, and less surface stoniness (Hatfield et al., 2017). Therefore, the southwestern section of BMEF is made up of deep mollisols or andisols (Alexander, 1994), which has a high moisture content, that likely favored the growth of the PP group observed in the PP-High cluster. In addition, ponderosa pine under deep soil conditions also develops a well distributed tap root system (Zouhar, 2001) that is likely related to their drought resistance and make them less susceptibel to insect and disease attack, hence high vigor (Merschel et al., 2021). On the contrary, species such as white fir and incense-cedar, which were intermixed within the Mixed-High and Mixed-Low clusters, are likely of low vigor compared to ponderosa pine, because of their shallow rooting habitat, which makes them susceptible to root diseases that additionally cause moisture stress and increase insect susceptibility (Zou et al., 2008; Zouhar, 2001). However, shallow rooting habitat, which is needed for plants to establish under frigid and shallow mollisols conditions of higher elevation within the northeastern, and southeastern 45 sections of BMEF (Alexander, 1994), could be the reason for the dominance of the CF group in these areas in my study (Zouhar, 2001). My findings are consistent with (North et al., 2002), who detected a cluster of incense-cedar and white fir within shallow outcrop areas in the southern Sierra Nevada, USA. The difference in soil moisture indicated by available water capacity at 150 cm pertaining to various soil characteristics and competition for resources such as light likely influenced the tree vigor in various clusters within BMEF. For example, low tree vigor in the PP-Low, Mixed-High, and Mixed-Low clusters compared to the PP-High cluster could likely be the result of soil moisture stress (McCullough & Wagner, 1987). In low soil moisture conditions, more carbohydrates are allocated for the repair of roots than towards the crown development, which substantially depletes the carbohydrate reserve and reduces tree vigor (Grulke et al., 2020; McCullough & Wagner, 1987). In addition, under the high stand density (Mixed-High, Mixed-Low, and PP-Low) at high elevation, competition for resources such as light likely resulted in a decrease in chlorophyll content and other photosynthetic pigments, eventually reducing photosynthesis and tree vigor (Kocher & Harris, 2007; Schubert, 1974). Also, at high elevation, where the sub-soil depth increases, moisture becomes the limiting factor as moisture exists in deeper soil depth (Zouhar, 2001). As moisture becomes available at greater depth, the shallow rooting habit of the incense-cedar and white fir species prevents them from drawing moisture from greater sub-soil depth causing moisture stress and improper growth and development (Zou et al., 2008; Zouhar, 2001). Therefore, my results from the CART analysis and the interaction term in the multinomial logistic regression model suggested that the Mixed-Low cluster was apparent in high elevation areas with higher sub-soil depth and low available water capacity at 150 cm depth. 46 3.5.4 Management implications My analysis of the pre-harvest forest inventory is unique in terms of the details captured regarding overstory stand structures at a landscape scale and provides one of the first notable records of pre-harvest conditions of ponderosa pine-dominated forests in northern California (Oliver, 2000). At the time of original data collection in 1934, BMEF had not yet been harvested (Oliver, 2000) and fire suppression had been implemented only since the late 1800s (Unpublished data from Carl Skinner). Hence, the spatial patterns I identified for 1934, represent a relatively intact forest, undisturbed by Euro-American settlers. Such forest still exhibited the overstory stand structures, species composition, and patterns of vigor structure expected in an ecosystem subjected to frequent fire. Thus, my results provide a reasonable baseline and context to guided monitoring of forest change within BMEF and similar ponderosa pine forests in northern California. The variation in the spatial patterns of overstory stand structures in my study conveys an important message to managers on landscape-scale heterogeneity that existed prior to any active human management. The entire BMEF landscape was heterogenous and non-uniform in terms of the spatial patterns of forest cluster types. Thus, my results on the pre-harvest spatial patterns of vigor by SDI can be used to improve management guidelines to develop spatially heterogeneous, resilient forests in the face of environmental change (Stephens et al., 2015). For example, within the ponderosa pine dominated forests, managers can aim at reducing the RD below 0.25 by retaining high vigor trees in clumps, like PP-High clusters in the southwestern section of my study area which is more resilient compared to the rest of the BMEF. In addition, quantifying the spatial patterns of tree vigor by SDI, prior to any management treatments, could help identify 47 dense low vigor stands at risk of insect, disease, and fire where management treatments can be prioritized. My study comes with several limitations. First, we cannot evaluate spatial patterns at a scale finer than 1 ha due to the design of the 1934 census. This could be important from a wildlife perspective (Larson & Churchill, 2012). Second, excluding trees < 9.1 cm (3.6 in) DBH leaves out the regeneration that was established around the early 1900s, which may have influenced tree vigor due to understory and overstory resources competition (Riegel et al., 1992). My study does not incorporate any data on disturbances such as fire and climate, which could have affected the spatial patterns in 1934 BMEF. Hence, my results provide an incomplete picture of the pre-harvest spatial patterns of the overstory stand structures of frequent fire ecosystems (Lydersen et al., 2013). Even in the areas where the fire regimes were still intact, due to the patchy nature of fire (Show & Kotok, 1924), small trees might still be present within BMEF, which might substantially alter the spatial patterns of the overstory. 3.6 Conclusions The management of current dry forest ecosystems depends on the comprehensive understanding of pre-harvest spatial patterns based on the three important components of overstory, which are stand density, tree species, and vigor. My study quantified the spatial patterns of stand density index (SDI) across vigor classes in two conifer species groups using landscape scale (~4000 ha) 1934 field data from BMEF. I found that high vigor forests with ponderosa pine as dominant species were spatially aggregated at low SDI, while low vigor forests were spatially aggregated with a mixture of incense-cedar, white fir, and ponderosa pine trees at high SDI. I also found that the spatial patterns of stand 48 density index by vigor and species groups were associated with elevation, soil moisture, and soil depth. Forest managers in California use pre-harvest forest conditions, which were resilient against insects, fire, and disease, as a guide to inform their management decisions. My results will help forest managers to understand that resiliency against insects, disease, and fire in the pre-harvest mixed-species dry forest comes from the spatial patterns of species vigor and was influenced by environmental factors. Forest managers can use my results as a baseline to make an informed decision regarding management prescriptions that will allow them to retain a tree species at various SDI and vigor combinations within different environmental gradients in a landscape. The results from my study are unique as they include the information on tree vigor, which is an underrepresented component of forest stand structure in the current body of literature. In addition, studies on pre-harvest spatial patterns often come from reconstruction and dendrochronological studies, and the current body of literature lacks research on spatial patterns based on field data. Hence, this research pioneers the description of pre-harvest spatial patterns using actual field data in a large landscape of 4,000 ha. 49 Chapter 4: Difference in Regeneration Conditions in Pinus ponderosa Dominated Forests in Northern California, USA, over an 83 Year Period 4.1 Introduction The resiliency of dry forests can be characterized by desired post-disturbance structure and composition, which are resistant to disturbances such as fire and insect attack. This requires the implementation of appropriate methods to support and enhance regeneration (Long et al., 2018). In dry forests of California, public agencies focus forest management efforts on reducing forest fuels by using thinning and prescribed fire treatments (McKelvey et al., 1996). However, a common secondary objective focuses on promoting more resilient future forests ideally dominated by fire-resistant and shade-intolerant ponderosa pine (Pinus ponderosa Laws.) (Fitzgerald, 2005). Meeting these two objectives requires an understanding of regeneration dynamics, which is a key component of silvicultural decision-making processes in sustainable forest management (Stein, 1992). Many studies on reference conditions have shown that conifer seed production, germination, subsequent establishment, and survival are associated with both biotic and abiotic factors (Puhlick et al., 2012). Abiotic factors such as fire influence the establishment and survival of seedlings (Allen et al., 2002; Brown & Sieg, 1996). By exposing mineral soil and preparing fine seed beds (Weaver, 1943), low- to moderate-severity fires prior to the early 1900s (Show & Kotok, 1924; Taylor, 2010) affected regeneration of ponderosa pine forests of the southern Cascades. Soil moisture (Schubert, 1974), and topography (Dodson & Root, 2013) are considered important limiting abiotic factors for tree regeneration and growth in northern California and the Sierra Nevada (Van Pelt & Franklin, 1999). Differences in regeneration 50 densities across elevation and moisture gradients have been reported (Dodson & Root, 2013). At higher elevation, regeneration densities are usually higher due to improved water availability (Dodson & Root, 2013). Successful regeneration also depends on stand structure (i.e., species composition, vertical and horizontal distribution of trees) (Angell et al., 2014) by modifying availability of water and nutrients (Oliver & Larson, 1996). Biotic factors such as disease or insect outbreaks, good seed crop (Cooper, 1960), and anthropogenic disturbances including prior silviculture treatments are capable of altering the trajectory of stand dynamics by creating establishment sites and resources for regeneration (Zald et al., 2008). Disentangling the relative importance of this suite of interacting biotic and abiotic forces with enough specificity to guide management prescriptions remains challenging, however. In ponderosa pine forests of western United States, prior to the late 1800s, regeneration densities were often low enough to maintain an open, park-like forest structure (Bailey & Covington, 2002). However, after the late 1800s, an absence of fire increased seedling densities producing a pulse of dense regeneration (McKelvey & Johnston, 1992). Subsequent increased tree densities resulting from fire exclusion provided a conducive environment for establishment of more shade-tolerant species such as white fir (Abies concolor) and incense-cedar (Calocedrus decurrens (Torr.) Florin) under the dense shade of large trees (McKelvey et al., 1996; Van Wagtendonk et al., 2006; Zald et al., 2008). As they began to exhibit more homogeneous structure and composition (Kilgore, 1973; Pearson, 1950), many ponderosa pine forests were at greater risk of high-severity, stand-replacing fires (Stephens et al., 2015). Thus, the future trajectory of ponderosa pine forests depends on regeneration responses to management as well as environmental conditions (Puhlick et al., 2012). High levels of stand density tend to suppress individual tree growth and may lead to stand stagnation (Zald et al., 2008). Besides fire, 51 topography plays an important role in conifer regeneration with aspect and slope interacting to influence solar radiation, temperature, and soil moisture availability for germinating seedlings (Dodson & Root, 2013; Farnden, 1994). Yet, the direct effects of topographic and moisture gradients are not well understood. Much of our knowledge about regeneration under reference conditions is based on reconstructions using tree-ring data which may fail to distinguish among the effects of fire, topography, and soil characteristics in regeneration (League & Veblen, 2006). Both the fragmentary nature of tree-ring evidence (subject to decay or destruction by subsequent fires) as well as the difficulty of obtaining an annual resolution of germination dates from samples of large trees (Savage et al., 1996) present limitations. Tree-ring reconstructions also fail to document the survival rate of the seedlings established post-fire (Harley et al., 2018). Furthermore, many recent regeneration studies aimed at understanding the mechanisms of post-fire regeneration in ponderosa pine forests only help reveal contemporary conditions (Fulé et al., 2012; Keyser et al., 2008). Thus, in general, studies of post-fire regeneration in appropriate reference conditions for ponderosa pine forests are generally lacking. Forest inventories based on field surveys can provide quantitative measures of regeneration such as seedling density (Stein, 1992) but can suffer from some limitations. First, natural tree regeneration is known to be highly spatially heterogeneous and thus direct measures of regeneration density can be misleading. Second, seedlings often occur in clumps and field survey data tend to be highly aggregated with high variability (Lowdermilk, 1927). In order to address issues of skewness and variability in the collected regeneration data, this study used a variation of the stocked-quadrat method to estimate stocking proportion, which was originally introduced by Lynch & Schumacher (1941). Stocking in this study was defined as percent of 52 0.0016 ha plots with at least one tree seedling or sapling, and was modelled from available seedling density data using the equation from Ritchie (2020). I conducted this study in Blacks Mountain Experimental Forest (BMEF), because highly detailed regeneration data are available for yellow pine, incense-cedar, and white fir by four height classes, in a unique dataset collected in 1934 and 2018. The dataset from 1934 reflects conditions prior to active forest management and thus serves as a valuable form of reference conditions for comparison to contemporary forest conditions observed in 2018. Forest stand conditions in 1934 had likely already been changed due to sheep grazing, which may have altered fire frequency (Skinner & Taylor, 2006) relative to forest conditions prior to European settlement. Nevertheless, I will be using ‘reference conditions’ for the purposes of this work to answer the following four questions: 1. What was the stocking of ponderosa pine, incense-cedar, and white fir for four different height classes in 1934 prior to active forest management across BMEF? 2. Was stocking in reference conditions in BMEF associated with overstory stand structure, soil, and topography? 3. Are there differences in stocking among contemporary (2018) forests that received several management treatments (thinning, burning and combination thinning)? If so, which treatments differ? 4. How does stocking in reference conditions (1934) for the observed tree species differ from contemporary stocking (2018) for 12 compartments that received management treatments prior to 2018? My study allowed detailed quantitative descriptions of reference (1934) and contemporary (2018) stocking of different species under various stand conditions that were dominated by 53 ponderosa pine. I was also able to quantify differences in contemporary stocking—a given stand density based on a management objective or a desirable density—under several management scenarios. Information on regeneration in reference conditions by species, height classes, and environmental factors associated with regeneration is critically important to managers for prescribing forest management treatments that will determine sustainability of forest values in the future. 4.2 Materials and Methods 4.2.1 Study area This study was conducted at the Blacks Mountain Experimental Forest (Figure 4.1). A detailed description of the BMEF was given in Chapter 2 (see section 2.1). 4.2.2 Management history Due to observed high rates of insect-caused mortality (Eaton, 1959; Hallin, 1959), between 1937 and 1960, most compartments of the forest were subject to sanitation/salvage cutting that removed many of the decadent and high risk trees in the forest, those of risk class 3 and 4 (Salman & Bongberg, 1942). These harvests were usually light, removing less than 20% of the standing volume. Remnants of the earlier structure remained in four research natural areas (RNAs; Table 4.1, Figure 4.1). The RNAs represent undisturbed and late successional forests (Camann et al., 2008; Oliver, 2000). During the 20th century, there was no fire at BMEF since 1899 (unpublished fire scar data from Carl Skinner, United States Forest Services, retired). However, after 1997 prescribed fire was applied in two of the four RNAs. The Bull thinning project (hereafter, thinning from below) was started in 2011 and concluded in 2012. Prescribed burns within this project were performed in fall of 2017. Thinning in a Low Structural Diversity treatment, removing all but the mid-canopy trees (hereafter, 54 combination thinning), was conducted in 1996 and prescribed burns were performed in 1997 (Crotteau & Ritchie, 2014). All observations on current conditions were conducted in 2018. Treatments are summarized as follows (see details in Table 4.1): • Control (C) • Prescribed fire (B, burn only)—prescribed fire in 1999 / 1997 and 2011 • Combination thinning (TC)—thinning in 1996 • Combination thinning with prescribed fire (TC-B)—thinning in 1996 and prescribed fire in 1997 • Thinning from below only (TB)—thinning in 2012 • Thinning from below with prescribed fire (TB-B)—thinning in 2012 and prescribed fire in 2017 55 Figure 4.1 Location of Blacks Mountain Experimental Forest (BMEF) in northern California with 100 compartments and six management treatments. Treatments were applied at the compartment level. C: control, B: prescribed fire (1999/1997, 2011), TC: combination thinning (1996), TC-B: combination thinning (1996) with prescribed fire (1997), TB: thinning from below (2012), TB-B: thinning from below (2012) with prescribed fire (2017). 56 Table 4.1 Treatment, compartment name, management description, compartment area, and number of 0.0016 ha quadrats collected in 1934 and 2018 for each of the 12 compartments for which the non-parametric ANOVA and t-tests were conducted. Treatment Compartment Management description Area (ha) # of 0.016 ha quadrats 1934 2018 RNA C B24-4 (RNA A) No treatment 35 296 140 C C4-21 (RNA D) No treatment 38 210 108 B P21-28 (RNA B) Prescribed fire in 1999 39 288 157 B G2-12 (RNA C) Prescribed fire in 1997 and 2011 57 444 225 LOD TC G14-3 Overstory above 50.8 cm removed, and remaining stands thinned from below to a residual density of 250 tree per hectare in 1996 45 350 175 TC B24-7 Overstory above 50.8 cm removed, and remaining stands thinned from below to a residual density of 250 tree per hectare in 1996 41 322 158 TC-B G17-6 Same as G14-3 followed by prescribed fire in 1997 46 367 193 TC-B B27-11 Same as B24-7 followed by prescribed fire in 1997 77 642 210 Bull thinning TB G17-1 Thinning from below up to 40.38 cm diameter at breast height (dbh) in 2012. For clumps of trees >40.38 cm dbh, half removed while largest and healthiest trees retained. Trees between 7.62–40.38 cm DBH and not in aggregates of large trees, a residual density stocking of 13.77 to 18.37 m2 ha−1 of basal area retained. Retain all snags >30.48 cm DBH except where they pose a hazard adjacent to roads or other improvements. 32 230 131 TB G10-1 Similar to G17-1 in 2012. But trees between 7.62–40.38 cm DBH and not in aggregates of large trees, a residual density stocking of 9.18 m2 ha−1 of basal area was retained. 52 341 213 TB-B G7-6 Thinning from below to a target of 18–23 m2 ha−1. Prescribed fire in fall of 2017. 72 469 291 TB-B G2-8 Thinning from below to a target of 18–23 m2 ha−1. Prescribed fire in fall of 2017. 38 343 149 LOD: Low structural diversity treatment; RNA: research natural area. 57 4.2.3 Regeneration data Data for the analysis came from two sources: (i) 1934 reference conditions regeneration data from 100 compartments (Figure 4.1) and (ii) 2018 regeneration data—hereafter referred to as contemporary regeneration data—for 12 compartments spanning a range of current conditions at BMEF. 4.2.3.1 1934 regeneration data For the detailed description of 1934 regeneration data collected in 100 compartments (Figure 4.1) please refer to section 2.2.3. 4.2.3.2 Contemporary regeneration data In 2018, the 1934 sampling design was imitated within 12 compartments, which spanned across six recent treatments with two replications each across the entire BMEF (Table 4.1, Figure 4.1). The 2018 field sampling was specifically designed and conducted in order to compare the contemporary data to the available 1934 data, because little quality data on regeneration were collected since 1934. While in 1934 square 0.0016 ha (4-milacre) quadrats were used, I used circular plots of radius 2.27 m (7.45 feet) split into four 0.000404 ha wedges for ease of implementation in the field. Plots were established at an interval of 50 m spacing among permanently marked grid points at 100 m square spacing. The 2018 seedling data, which were collected within the milacre quadrats, were then aggregated to the compartment level as was done with the 1934 data. 4.2.4 Data compilation 4.2.4.1 Stocking Stocking, as defined by Ginrich (1967), is a relative term that describes a given stand density based on a management objective or a desirable density. For 1934, only the seedling density per 58 hectare summaries calculated from the 0.0016 ha (4-milacre) plots were available rather than the raw inventory data. Therefore, I was unable to observe stocking, defined as percent of plots with at least one tree seedling or sapling per 0.0016 ha, directly for the 1934 data. Density measurements alone can be misleading because of the known spatial heterogeneity of natural tree regeneration (Gill, 1950; Lowdermilk, 1927). There have been several attempts at relating understory tree density and stocking (Ghent, 1969; Gill, 1950; Lynch & Schumacher, 1941), because stocking is a more reliable measure in assessing management objectives or desirable tree occupancy (Ginrich, 1967) and accounts for the spatial heterogeneity of natural regeneration. More recently Ritchie (2020), using data from the 2000 Storrie fire on the Lassen National Forest, conducted a validation of the Lynch & Schumacher (1941) probit model that predicts stocking as a function of understory tree density, and found that the Storrie fire data matched well with the Lynch & Schumacher (1941) model fit. The results from Ritchie (2020) were applied here to estimate stocking from the compiled, available seedling density data: Φ−1(𝑝) = −0.21307 + 1.22876 𝑙𝑜𝑔10(𝑥) (1) where Φ-1(p) is the inverse of the standard normal (probit) of proportion (p) of stocked plots and x is mean number of trees per plot. To evaluate the fit of the new equation I plotted my observed stocking in 2018 at 0.0016 ha plots over the predicted stocking 95% confidence interval from Ritchie (2020). 4.2.4.2 Explanatory variables for 1934 conditions Three groups of explanatory variables were used: overstory, soil, and topography. I compiled overstory variables—trees per hectare and overstory basal area per hectare (Table 4.2)—for each of 100 compartments using available census data from 1934. In 1934, all live trees >8.9 cm DBH were tallied by 5.08 cm wide diameter classes within 4055 plots (~1 ha in size) throughout 59 BMEF. If a 1 ha plot straddled two compartments, the overstory variables were weighted by the proportion of the plot that fell into a particular compartment. Two compartments that were primarily composed of meadows were dropped, and I used the remaining 98 compartments for any analyses that included the 1934 overstory data (Table 4.3). Topographic variables—aspect, slope, and elevation—were extracted from a digital elevation model (DEM) at a 1-m resolution (Table 4.2). The zonal statistic function from ArcMap 10.4.1 summarized the values of the DEM in terms of average elevation, slope, and aspect within the extent specified by the shape file for the 100 compartments (Wilson, 2012). Aspect and slope were incorporated into the model using the transformations suggested by Stage (1976). Vegetation information for BMEF was available from an ecological unit inventory (EUI) shapefile. This EUI map, developed for BMEF in 1994, identified different geological landforms, soils, and potential vegetation types (Alexander, 1994). The compartment map was overlaid with the EUI map and the intersect function in ArcMap 10.4.1 was used to extract the soil variables (Table 4.2) from the EUI map at compartment level (Reed et al., 2016). Multiple categories of variables could fall in a single compartment. If a compartment straddled multiple soil variable categories, an area-weighted average of the midpoints was assigned to the compartment. Hydrological processes and the control of topography over such processes can be quantified by the topographic wetness index (Beven & Kirkby, 1979). For the 100 compartments, topographic wetness index was calculated based on the slope from a 1 m resolution DEM, using the equations provided by Beven & Kirkby (1979). Similarly, solar insolation, which takes into account topographic characteristics such as slope, aspect, and elevation (Dubayah & Rich, 1995), was calculated for each of the 100 compartments using the DEM and Spatial Analyst tool from ArcMap 10.4.1. 60 Table 4.2 Compartment-level response and explanatory variables. Variables group Source Response variable Stocking Field data (1934) for 100 compartments Field data (2018) for 12 compartments Explanatory variables Overstory Trees per hectare (TPH) Field data (1934) Basal area per hectare (BA, m2 ha-1) Field data (1934) Topographic variables Aspect (degrees) DEM at 1 m resolution Slope (%) DEM at 1 m resolution Elevation (m) DEM at 1m resolution Soil variables AWC150 (cm) EUI Map (Alexander, 1994) AWC50 (cm) EUI Map (Alexander, 1994) Sub-soil depth (cm) EUI Map (Alexander, 1994) Surface soil depth (cm) EUI Map (Alexander, 1994) Surface-stoniness (%) EUI Map (Alexander, 1994) Topographic wetness index DEM at 1 m resolution Solar insolation DEM at 1 m resolution AWC50: available water capacity at 50 cm depth, AWC150: available water capacity at 150cm depth, DEM: digital elevation model. EUI = ecological unit inventory. 61 Table 4.3 Summary statistics of stocking of all species, seedling density of all species, available water capacity at 150 cm (AWC150), overstory basal area per hectare for all species, and overstory density for all species for the 98 compartments in 1934 that were used in the model development. The overstory basal area was calculated for all trees >8.9 cm in diameter at breast height. Variables Min First Quartile Median Mean Third Quartile Max Stocking 0.21 0.52 0.61 0.60 0.70 0.79 Seedlings per hectare 827 4097 6302 7070 10,051 16,772 AWC150 (cm) 4.5 12.5 15.0 14.7 17.7 21.6 Overstory basal area (m2 ha-1) 5.31 12.86 14.74 14.77 16.92 25.47 Overstory trees per hectare (TPH) 61 100 116 124 140 229 Min: minimum value, Max: maximum value. 4.3 Statistical analysis 4.3.1 Stocking in reference conditions I conducted exploratory data analyses for the 1934 regeneration data using histograms, with stocking in each compartment and frequency on the x- and y-axis, respectively. Similarly, I constructed histograms for stocking by height class and species for the reference conditions data from 1934. 4.3.2 Effect of overstory, soil, and topography on stocking Meadows were part of the reference conditions forest landscape. For two out of the 100 compartments, the area of meadows exceeded 70%, which made these compartments different 62 from the other compartments, for which the area of meadows was below 40%, with 77 compartments not having any meadows. Therefore, I dropped the two compartments that were dominated by meadows from my analysis and used the remaining 98 compartments in the model development. Stocking is continuous and restricted to the standard unit interval (0, 1). Therefore, I fit beta regression models that take the distribution of stocking into account (Ferrari & Cribari-Neto, 2004) using the ‘betareg’ package in R version 3.6.2 (core Team, 2018) with a logit link function (Ferrari & Cribari-Neto, 2004). I developed the models as follows: My first step was to identify the variables that were significantly associated with stocking. I fit univariable beta regression models for all explanatory variables (Table 4.2). Then, I selected variables with low p-value (p < 0.05) to be used in multivariable beta regression models across the three groups of explanatory variables (Table 4.2), keeping variables significant at α-level of 0.05 and dropping variables that had variance inflation factor greater than three (Zuur et al., 2009). This resulted in a model that contained available water capacity at 150 cm (AWC150) and overstory basal area per hectare as explanatory variables. However, in my exploratory data analysis, the overstory basal area per hectare showed a curvilinear relationship with stocking. In sorder to account for the trend observed in the data and to mimic this trend, I decided to use a range of exponents between 0.10 to 2 at the interval of 0.10 for the overstory basal area per hectare to identify the value of exponent that resulted in the best model. I selected the model that resulted in high pseudo-R2 and high precision parameter (phi) values (Ferrari & Cribari-Neto, 2004) as well as a low Akaike Information Criterion (AIC) value (Akaike, 1973). I also checked the assumption of constant variance and goodness-of-fit using residual vs. linear predictor and half normal plots (Moral et al., 2017), respectively, for all beta regression models (Beven & Kirkby, 1979; Dubayah & Rich, 1995). 63 4.3.3 Stocking in contemporary forests under different management treatments The effect of six different management treatments on contemporary stocking was investigated using a Kruskal–Wallis non-parametric one-way ANOVA (Kruskal & Wallis, 1952). The treatment factor had six levels (Table 4.1). When differences among treatment levels were detected, a pairwise Wilcoxon rank sum test (Wilcoxon, 1945) was used to determine which treatment levels had different stocking. 4.3.4 Stocking in reference vs. contemporary forests Stocking by species and height class was available for 12 compartments in both measurement years. Dot plots were constructed for both measurement years for overall stocking, stocking by three species, and stocking by four height classes to investigate if contemporary stocking was higher than stocking in reference conditions. 4.4 Results 4.4.1 High stocking for ponderosa pine species in reference conditions My results indicated that the stocking observed in 2018 for 0.0016 ha plots corresponded well, for observations below about 1000 seedlings per hectare, with the 95% confidence limits of predicted 0.0016 ha stocking from Ritchie's (2020) equation (Figure 4.2). Although there is evidence of some lack of fit at the highest levels of density where model predictions of stocking may be low by around 10%. However, since for very high levels of density the model approaches one, errors in estimation will actually start to decrease for the extreme density levels observed in 1934. In 1934, the majority of compartments exhibited higher seedling densities than observed in the 12 compartments for 2018 (Table 4.3). Therefore, it is possible that the stocking predicted from Ritchie’s (2020) equation may be a bit lower than actual stocking. However, I do not know for sure how the equation performs for high seedling densities over 1500 trees ha−1, because I did 64 not have any 2018 data in that range for validation (Figure 4.2). The estimated stocking in 1934 was dominated by ponderosa pine (Figure 4.3b). My data indicated that there were no white fir and incense-cedar seedlings in 38 compartments in 1934. Most of the incense-cedar and white fir regeneration was present in the compartments at higher elevations, located in the northeast portion of BMEF. Ponderosa pine had the highest stocking (mean proportion = 0.53 and median proportion = 0.54), followed by white fir (mean proportion = 0.10 and median proportion = 0.21) and incense-cedar (mean proportion = 0.10 and median proportion = 0.02) (Figure 4.3). The majority of the compartments were dominated by seedlings in height classes 2 and 3 irrespective of species (Figure 4.4d–i). However, among the three species ponderosa pine had the highest stocking in height classes 2 (mean = 0.33 and median = 33) and 3 (mean = 0.33 and median = 0.34) (Figure 4.4d, g). Stocking in height classes 1 and 4, were low compared to other height classes for all three species (Figure 4.4). 65 Figure 4.2 Observed stocking proportion of natural regeneration in 2018 (black points) with model fit of 0.0016ha stocking (dashed line) and bootstrapped 95% confidence interval (grey) from Ritchie (2020). 66 Figure 4.3 Histograms of 1934 stocking proportion calculated at 0.0016 ha for 100 compartments for a) all species, b) Ponderosa pine, c) white fir, and d) incense-cedar, with mean (solid black line) and median (dashed line). The scale on y-axis are different for both panels of figure. 67 Figure 4.4 Histograms of 1934 stocking proportion at 0.0016 ha for 100 compartments by four height classes and three species:a) ponderosa pine, height class 1, b) incense-cedar, height class 1, c) white fir, height class 1, d) ponderosa pine, height class 2, e) incense-cedar, height class 2, f) white fir, height class 2, g) ponderosa pine, height class 3, h) incense-cedar, height class 3, i) white fir, height class 3, j) ponderosa pine, height class 4, k) incense-cedar, height class 4, l) white fir, height class 4. Height Class 1: 0–0.31 m, Class 2: 0.31–0.91 m, Class 3: 0.91–1.83 m and Class 4: >1.83 m for seedlings less than 8.89 cm diameter at breast height (DBH), with mean (solid black line) and median (dashed line). The data used for ponderosa pine includes all 100 compartments while the data for incense-cedar and white fir were from the 62 compartments that contained regeneration of both species or either of them. 68 4.4.2 Overstory and soil variables were associated with stocking I chose the model with the lowest AIC value (−191.86) and highest pseudo-R2 (0.49) and precision parameter (29.27) compared to the other tested models (Table 4.4). However, given the small differences in AIC values among the models (Table 4.4), this model is not necessarily superior to the other models in terms of model fit. This final beta regression model showed that basal area per hectare raised to an exponent of 0.10 (p-value < 0.001) and available water capacity at 150 cm (p-value < 0.001) were significantly associated with reference conditions stocking (Table 4.4). Table 4.4 Pseudo-R2, Akaike Information Criterion (AIC), and precision parameter for selected beta regression models. Estimated coefficients with p-value in parentheses are given for each explanatory variable. The final beta regression model is highlighted in black font. Coefficients (p-value) Explanatory variables Pseudo R2 AIC Precision parameter (phi) AWC150 (cm) Overstory basal area (m2 ha-1) AWC150 + (Basal area)0.10 0.49 −191.86 29.27 0.059 (<0.0001) 11.05 (<0.0001) AWC150 + (Basal area)0.20 0.483 −190.46 28.85 0.059 (<0.0001) 4.251 (<0.0001) AWC150 + (Basal area)0.30 0.475 −189.00 28.43 0.059 (<0.0001) 2.178 (<0.0001) AWC150 + (Basal area)0.40 0.467 −187.51 28.00 0.059 (<0.0001) 1.253 (<0.0001) AWC150 + (Basal area)0.50 0.458 −185.98 27.56 0.058 (<0.0001) 0.768 (<0.0001) AWC150+ (Basal area) 0.413 −178.13 25.45 0.057 (<0.0001) 0.099 (<0.0001) AWC150+ (Basal area)1.50 0.365 −170.47 23.54 0.055 (<0.0001) 0.016 (<0.0001) AWC150+ (Basal area)2 0.319 −163.51 21.92 0.053 (<0.0001) 0.003 (<0.0001) AWC150: available water capacity at 150 cm depth, Basal area: overstory basal area in m2 ha−1. 69 Stocking appeared to increase with increasing overstory basal area per hectare up to about 12 m2 ha−1; above this level the relationship hits a plateau (Figure 4.5a). Available water capacity at 150 cm had a positive relationship with stocking. Stocking increased as more water was available to plants at depth of 150 cm (Table 4.4). The slope-aspect transformations (Stage, 1976) had p-values between 0.04 and 0.07. Because the inclusion of these model terms did not improve the model substantially but increased model complexity, I decided not to include them. None of the remaining topography variables—elevation, solar insolation, and topographic wetness index—were significant in the model. 4.4.3 Differences in stocking among management treatments in 2018 Stocking in 2018 varied across treatments (Kruskal–Wallis test p-value < 0.0001) (Figure 4.6a). Control (C) compartments had the highest stocking followed by burned (B) compartments (p-value < 0.0001; Figure 4.6a). The Wilcoxon rank sum test revealed that stocking in TC-B and TB were not significantly different from each other (p-value = 0.0816). The difference in stocking between TC-B and TB-B (p-value = 0.0028) and TC and TB-B (p-value = 0.0014) was significant. Significant differences in stocking (p < 0.0001) were found between all other treatment levels. 4.4.4 Stocking in reference conditions higher than contemporary stocking Estimated stocking in 1934 was greater than observed stocking in 2018 across all the compartments (Figure 4.6a). Control treatments had the highest stocking in 2018 compared to the other treatments. These two compartments also had the highest stocking in 1934 (Figure 4.6a). The majority of compartments were dominated by stocking of ponderosa pine followed by incense-cedar in both reference 1934 and contemporary 2018 conditions (Figure 4.6b). Stocking in 1934 was dominated by height classes 2 and 3 (Figure 4.6c). In 2018 stocking was dominated 70 by height class 4 in control compartments (Figure 4.6c), while compartments that received the burn treatment were mostly stocked with seedlings from height classes 1 and 4 (Figure 4.6c). The compartments that received thinning treatments had very low stocking, with height classes 1 and 2 being a bit more abundant than height classes 3 and 4 (see TC-B, TC and TB in Figure 4.6c). Figure 4.5 Relationship between overstory basal area per hectare and stocking proportion calculated at 0.0016 ha for 98 compartments in 1934 (a). Relationship between stocking proportion by species and corresponding overstory basal area of that species for (b) ponderosa pine, (c) incense-cedar, and (d) white fir. The highlighted grey area represents the range of incense-cedar and white fir overstory basal area compared to ponderosa pine. 71 Figure 4.6 Stocking proportion calculated at 0.0016 ha in 1934 (red) and 2018 (blue) for 12 compartments by (a) years, (b) species, and (c) height class for various management treatments. Height Class 1: 0–0.31 m, Class 2: 0.31–0.91 m, Class 3: 0.91–1.83 m and Class 4: >1.83 m for seedlings less than 8.89 cm diameter at breast height (dbh). IC: incense-cedar, PP: ponderosa pine, WF: white fir. C: control, B: prescribed fire, TC: combination thinning, TC-B: combination thinning with prescribed fire, TB: thinning from below, TB-B: thinning from below with prescribed fire. 72 4.5 Discussion My research was able to meet the following four objectives, as follows. First, I wanted to understand stocking densities for different species and height classes within reference forests. I determined that prior to the advent of forest management, reference stocking densities were dominated by ponderosa pine in height classes 2 and 3 throughout BMEF. Second, I aimed to determine the factors associated with stocking densities in reference forests and I determined the major drivers were overstory basal area and soil moisture. Third, I wished to understand whether contemporary stocking densities varied across various management treatments and I determined There were differences in stocking among management treatments, where combination thinning with prescribed burn reduced the stoking in contemporary forests. Lastly, I compared contemporary stocking densities with those of reference conditions and determined stocking was higher in reference forests than in contemporary forests. 4.5.1 Stocking in reference conditions was dominated by ponderosa pine at low elevations In the early 1930s, after sheep grazing (Skinner & Taylor, 2006) but prior to the onset of contemporary timber harvesting activities, low elevation sites of BMEF were dominated by ponderosa pine stands (Ritchie, 2016). The landscape in BMEF has gentle slopes with the majority of the forest in south or southwestern aspect and with an average slope of less than 15% (Ritchie, 2016). Such south-facing gentle slopes have been reported to favor seedling germination of ponderosa pine in northern California (Fowells & Stark, 1965). The observed dominance of ponderosa pine stocking in these low elevation sites of BMEF can further be explained by overstory basal area per hectare being a major driver of stocking in the reference conditions BMEF. Ponderosa pine dominated the overstory, and its presence not only facilitates 73 the establishment of regeneration by supplying ample seeds (Shive et al., 2018), but also ameliorates the harsher environmental conditions caused by high temperatures, solar radiation, and wind (Stein & Kimberling, 2003). In the reference conditions, 96% of BMEF fell within the moderate range for basal area (7 to 23 m2 ha−1), a range reported to be conducive for establishment of ponderosa pine seedlings (Flathers et al., 2016; Puhlick et al., 2012). More specifically, in 1934, BMEF overstory exhibited low to moderate basal area for ponderosa pine (0 to 25 m2 ha−1) and low basal area (0 to 7 m2 ha−1) for incense-cedar and white fir. The observed low to moderate basal area of the overstory maintained an open park-like structure in the majority of BMEF, which favored germination and growth of shade-intolerant ponderosa pine seedlings. The observed high stocking densities for incense-cedar and white fir in some compartments in the northeast portion of BMEF may be due to being situated in more mesic and higher elevation sites. Mesic and higher elevation sites have been reported to be more conducive to these species (Fowells & Stark, 1965). A mixture of all three species in the overstory may have created a dense overstory, which provided enough shade for the growth of shade tolerant incense-cedar and white fir (Oliver & Dolph, 1992). At the same time, the presence of mature incense-cedar and white fir trees provided the seeds necessary for germination (Oliver & Dolph, 1992). 4.5.2 Stocking in reference conditions was dominated by height classes 2 and 3 Land survey records across dry forests in Oregon, California, and parts of northern Arizona in the late 1800s showed that regeneration—seedlings and/or saplings—were observed over 35%–57% of the dry-forests (Baker & Williams, 2015). The absence of surface fires after Euro-American settlement failed to control their growth and establishment in these areas (Baker & 74 Williams, 2015). At the time of the 1934 regeneration data collection, BMEF had experienced over 35 fire-free years (unpublished fire scar data from Carl Skinner, United States Forest Services, retired). This fire free period may be attributed to heavy grazing by sheep in the area (Skinner & Taylor, 2006). Hence, seedlings that were established around the early 1900s due to the absence of fire could have stagnated in height growth throughout BMEF. Stand stagnation could explain the observed higher stocking in height classes 2 (0.31–0.91 m) and 3 (0.91–1.83 m) compared to height classes 1 and 4 throughout BMEF in 1934. Lacking age data for the 1934 forest conditions, I was unable to make any specific comments on height and age dynamics of seedlings within the ponderosa pine forest in the BMEF reference conditions. Oliver (1967) reported dense thickets of saplings that had an average height of 2.7 m at the age of 70 within BMEF as a result of stand stagnation. Similarly, Weaver (1943) reported stand stagnation after the establishment of regeneration in ponderosa pine forests in Washington and Oregon, USA, because of fire exclusion. Therefore, the seedlings/saplings that were observed in smaller height classes 1 (0–0.31 m) and 2 could be as old as seedlings/saplings in height classes 3 and 4 (>1.3 m) due to stagnation. Differences in stocking for various height classes could also be attributed to two causes. First, a lag in regeneration—a period of no regeneration—between two pulses of regeneration could have occurred (Cooper, 1960). Regeneration lags can occur due to the periodicity of cone production, whereby heavy cone production years are followed by years of few or no cone production (Shepperd et al., 2006). Ponderosa pine can experience regeneration lags from 15 to 70 years between establishment of seedlings (Lauvaux et al., 2016; Kolb et al., 2007) which can be influenced by unusual climatic episodes. Fluctuations between wetter fire-free and drier fire-prone periods can influence regeneration (Keyes & Manso González, 2015) as long as drier and 75 fire-prone periods will result in a longer gap in regeneration before another good seed year (Savage & Mast, 2005). Such gaps in regeneration might result in the differentiation of seedlings into various height classes. The differentiation of stocking in various height classes could also be attributed to variability in growth due to site effects (Oliver & Power, 1971) such as localized competition for nutrients and soil moisture (Schubert, 1974). For example, seedlings germinated at the same time in favorable sites might gain greater height than seedlings in poorer sites. In the western Sierra Nevada, southern Cascade Range, and Warner Mountains of northeastern California, Oliver & Power (1971) found that seedlings that germinated at the same time reached diameter at breast height at 13.5 years of age on poor-sites as compared to 4.4 years on better sites. 4.5.3 Available water capacity at 150 cm and overstory basal area were drivers of stocking in reference conditions Interactions of topographic variables—elevation, slope, and aspect—result in heterogeneous forest structure in the dry forests across the Sierra Nevada and Southern Cascades, usually through changing the patterns and amount of regeneration along topographic gradients (Lydersen & North, 2012). However, in my study, topographic variables such as elevation, slope, and aspect were not found to be major drivers of stocking. This is possibly due to the limited range in elevation—1700 to 2100 m—and a gentle slope of less than 15% on average (Skinner & Taylor, 2006). A good seed crop followed by a favorable moisture regime and receptive seed bed is needed for the establishment of ponderosa pine seedlings (Cooper, 1960). I observed that within BMEF available water capacity at the depth of 150 cm was one driver of stocking in 1934. A large number of seedlings possibly germinated due to ample moisture availability within their 76 rooting zone (Fowells & Stark, 1965; Haase, 1986; Harr, 1977). Ponderosa pine seedlings have been shown to survive in areas with low precipitation by maintaining high water use efficiencies and root systems—12 to 50 cm deep—that utilize moisture available in their rooting zone (Daubenmire, 1968; Haase, 1986). Stocking showed a curvilinear relationship with increasing basal area in the overstory (Figure 4.5a). Overstory basal area was dominated by ponderosa pine. However, the highest observed stocking and overstory densities in higher-elevation compartments in the northeast portion of BMEF could be due to the presence of all three species (Clark et al., 1999; Oliver & Dolph, 1992; Zhang et al., 2013). While the observed relationship between overstory basal area per hectare and stocking was stronger than the relationship between available water capacity at the depth of 150 cm, the overall model fit was not very good (low pseudo-R2 of 0.49). Stocking tended to be slightly overpredicted for lower ranges of overstory basal area, and variability in the data was quite high. The fitted beta regression model provides general insights about the variables that explain reference conditions stocking in the BMEF environment. My model, however, comes with the limitations that the soil variables used in the model were derived from an Ecological Unit Inventory (EUI) map produced in 1994 (Alexander, 1994) decades after 1934, when the reference conditions regeneration data were collected. In the absence of field-collected soil data, I have to make an assumption that the information on soil variables derived from the EUI map represents the soil conditions that existed in BMEF around 1934 and has not changed. Additional variability in stocking could possibly be explained by incorporating fire history variables such as time since fire, fire frequency, and fire return interval into the model (Taylor, 2010). However, in this study, availability of fire history variables was limited to 28 of the 100 compartments. The 77 fire history variables were statistically significant when I added them to the model using the 28 compartments. Considering the small sample size of 28, the model was overfit with these additional variables and I did not think it was meaningful to report the results in detail. Yet, this exploratory analysis confirms that fire history variables could explain some of the observed variability in stocking. 4.5.4 Contemporary stocking was low in thinned stands and differed across treatments In 2018, compartments that received thinning treatments—thinning from below and combination thinning—with or without prescribed fire had lower stocking than the compartments that received the burn only treatment. The mechanical thinning treatments mostly removed the trees up to a diameter of 40.38 cm including the seedlings and saplings up to a diameter of ~7 cm, which can explain the low observed stocking in these compartments (Hallin, 1959; Ritchie, 2016). In addition, within the compartments that received prescribed fires after the thinning treatments, fire might have killed additional seedlings resulting in lower stocking compared to treatments that received thinning only. Another possible reason for low stocking within the combination thinning and thinning from below could be the lack of seed sources near thinned sites (Clark et al., 1999; Zhang et al., 2013). Within the thinned sites, due to reductions in canopy cover, seedlings are likely to get exposed to high levels of direct solar radiation (Keeton & Franklin, 2005; McDonald, 1983). This may lead to seedling desiccation and death within thinned compartments (Fowells & Stark, 1965), where regeneration was already removed during the thinning treatments. In contrast, higher stocking in the burn only treatment compared to the treatments that received thinning can possibly be attributed to the fire having reduced heavy accumulation of litter and duff, which allowed seeds to reach mineral soil and extract adequate soil moisture for 78 germination (Haase, 1986). At the same time, the overstory that survived the prescribed fire, provided favorable environmental conditions at microsites for germination and growth of seedlings (Haase, 1986). The presented results for the burn only treatment should be taken with caution, as time since management might have played an important role in determining the amount of stocking across various management scenarios at BMEF. Burn only treatments were implemented in 1997 while the compartments that received thinning from below in 2012 were burned in 2017 (Crotteau & Ritchie, 2014). I collected the regeneration data in 2018. There was a gap of 20 years between application of treatments and data collection in the burn only treatment. During this 20-year window, seedlings might have established in the burn only treatments and resulted in the observed higher stocking. In the compartments that received thinning from below and combination thinning with or without fire, there was a gap of only one year between application of fire and data collection. This may explain the observed lower stocking in these treatments. Because time since treatment was different for the different treatments, the observed treatment effect may be confounded with time. The observed lower stocking in burned compartments compared to the control compartments may be explained by the mortality of the seedlings over time due to crown scorch and cambium wounding (Bailey & Covington, 2002). In the absence of any management activities and fire for many decades in the control compartments (Ritchie, 2016), the seedlings could have stagnated within the stands during various episodes of recruitment, which resulted in high stocking (Oliver, 1967). In addition, the control compartments were located at somewhat higher elevations and within mesic sites than the burn only treatments. These conditions are conducive in terms of available microsites for shade-tolerant white fir and incense-cedar (Zald et al., 2008), which was confirmed by the higher stocking observed for incense-cedar and white-fir 79 in the control compartments. Therefore, within these compartments higher stocking may also be explained by the higher abundance of seedlings of the shade-tolerant species and a more mesic environment, rather than by the fact that they did not receive any management. 4.5.5 Higher stocking in reference conditions than in managed, contemporary stands The observed lower stocking in 2018 compared to 1934 may indicate that the contemporary forest structure is not as conducive to regeneration as it used to be prior to 1934. Due to infilling of trees, there may be lack of growing space in contemporary forests (Knapp et al., 2013). In 1934, however, the forest structure was more open, maintained by frequent surface fires, which reduced shading and resource competition and prepared fine mineral seedbeds for seedling germination (Bailey & Covington, 2002). However, the differences in stocking between reference and contemporary forests could also simply be due to the management treatments that were applied in the contemporary forests, where the mechanical thinning treatments removed most of the regeneration and prescribed fires may have killed additional regeneration (Hallin, 1959; Ritchie, 2016). 4.6 Conclusions The trajectories of future forest stand structures in ponderosa pine forests of northern California depend on the regeneration response to management and environmental conditions. Regeneration studies in forests representing reference conditions allows managers to understand the relationships between regeneration and environmental factors, which in turn shape overall stand structure. My field inventory data-driven study showed dense regeneration in reference conditions, dominated by shade-intolerant ponderosa pine seedlings (in height classes 2 and 3) at the lower elevations of the forest. More mesic sites at higher elevations exhibited high stocking of shade-tolerant incense-cedar and white fir. Seedlings that were established around the early 80 1900s found available growing space in open park-like conditions favoring the germination and growth of ponderosa pine seedlings. Additionally, in absence of fire to check the growth, seedlings could have stagnated in height growth throughout BMEF, which resulted in high stocking in height classes 2 and 3. Overstory basal area per hectare and available water capacity at 150 cm depth were related to stocking in reference conditions. Stocking increased with both increasing overstory basal area as well as increasing available water capacity. Topography variables were not significant in the final beta regression model. The amount of regeneration was different across management treatments in contemporary conditions. The time since management played a key role in determining the differences in stocking within contemporary BMEF. Compartments that were recently thinned and burned had lower stocking compared to compartments that were burned and thinned 20 years prior. The stocking in the contemporary BMEF was lower than the stocking in the reference conditions forests. The primary reason for lower stocking in contemporary conditions was because most of the small trees, seedlings, and saplings were either removed by thinning or killed due to fire by recent management activities. The results of this study are unique because much of our current knowledge on regeneration in reference conditions has been based on dendro-chronological reconstructions, which cannot provide the kind of detailed information about regeneration in reference conditions that I have presented in this study. 81 Chapter 5: Quantifying changes in spatial patterns of Pinus ponderosa dominated overstory above-ground biomass across 80 years using LiDAR and inventory data 5.1 Introduction Overstory stand structure of ponderosa pine-dominated forests—hereafter “dry forests” (Harrod et al., 1998)—has undergone substantive changes since Euro-American settlement of California, USA. Such changes have included shifts in species composition, an increase in the density of small trees (Dolanc et al., 2014), and a decline in the density of large trees (Lydersen et al., 2013; North et al., 2009; van Mantgem & Stephenson, 2007). 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 (Dolanc et al., 2014; Lydersen et al., 2013). In California, tree-ring reconstructions indicate that heterogeneous spatial patterns in the overstory of dry forests prior to 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), resulting in open forests with sparse large individuals of fire-resistant species (Abella & Denton, 2009; Churchill et al., 2013). Such spatial patterns are the result of interactions among prior processes such as fire, regeneration, competition, and mortality (Moore & Carpenter, 1999; Youngblood et al., 2004). However, overstory spatial patterns, especially of above-ground biomass (AGB) have substantially changed following Euro-American settlement (Clyatt et al., 2017; Collins et al., 2017). Within dry forests of California, overstory AGB (Mg ha-1) has substantially increased following Euro-American settlement and fire exclusion with increases in tree density of shade-82 tolerant species such as white fir (Abies concolor (Grod. & Glend)) and incense-cedar (Calocedrus decurrens (Torr.)) (Clyatt et al., 2017; Collins et al., 2017; Oswalt et al., 2014). This increase in overstory AGB has altered the fuel dynamics of frequent-fire forests, in turn creating forest overstory stand structures and spatial patterns more prone to stand-replacing wildfires (Hurteau et al., 2016). Frequent fires with return intervals of less than 25 years are among the most important drivers of spatial patterns in dry forests (Agee & Skinner, 2005; Kilgore, 1973). However, it is unclear how overstory spatial patterns may have been altered due to the exceptionally long fire-free period after Euro-American settlement (Lydersen et al., 2013; Taylor, 2010). Therefore, quantitative descriptions of changes in spatial patterns of overstory AGB are crucial for land managers as they can be used to assess the potential of wildfire due to increased fuel load in overstory (Chen et al., 2004; Kashian et al., 2006). Shifts in spatial patterns over time have rarely been addressed in the literature, which is potentially a result of both a general lack of spatially-explicit data collected prior to active management (Lydersen et al., 2013) 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 (Hall et al., 2005). 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 (Hall et al., 2005). Although LiDAR data are being extensively used in enhanced forest inventory, ground data are still required (McRoberts et al., 2010). LiDAR derived variables such as height metrics can be utilized with ground-based observations such as basal area (m2 ha-1), volume (m3 ha-1), and AGB (Mg ha-1) to predict forest overstory stand structure at different times (Lefsky et al., 1999; Means et al., 2000). Thus, LiDAR data 83 complements and can be used in conjunction with ground-based inventories to identify changes in the spatial patterns of forest overstory structure (Skowronski & Lister, 2012). Furthermore, many studies examining shifts in spatial patterns of forest structure in California, USA, have been restricted to particular elevation ranges (Dolanc et al., 2014; Millar et al., 2007), management units (Beaty & Taylor, 2008; Vankat & Major, 1978), certain functional types (Bolsinger, 1988), and small study areas (Ansley & Battles, 1998; Smith et al., 2005). Numerous studies have provided a general description of dry forests prior to active management following Euro-American settlement, however, a description of the spatial patterns of overstory stand structure and their changes over time is lacking (e.g., Knight et al., 2020, Fry et al., 2014). Many studies have investigated spatial patterns of tree density and tree size (e.g., Sanchez Meador et al., 2010; Youngblood et al., 2004), but generally ignored the overstory AGB component. Some studies have addressed the AGB component of overstory stand structure in terms of AGB averages and distribution, however, these studies have not addressed the spatial components and differences in spatial patterns of overstory AGB at different points in time (e.g., Ritchie, 2016; Sheridan et al., 2015). My study was conducted in the dry forests of Blacks Mountain Experimental Forest (BMEF) in northern California, USA. BMEF remained untouched by selective logging operations of the late 1800s until the 1940s (Oliver, 2000). Within the BMEF, I focused on the section defined as the Blacks Mountain Ecological Research Project (BMERP), initiated in 1991. Within BMERP various thinning and prescribed burn treatments were implemented to create different structural diversity in forest stand structure (Oliver, 2000). To quantitatively describe differences in spatial patterns of AGB, I utilized three types of data collected for the BMERP project area: ground inventory data from 1934 (as reference data), as well as contemporary data 84 from LiDAR (2015) and ground verification data (2016). The specific objectives of our study were to: a) Compare spatial patterns of overstory AGB between contemporary forests and reference conditions using 2015 LiDAR data and 1934 census data, and b) Compare the differences in spatial patterns of overstory AGB among various management treatments in contemporary conditions. Understanding the differences in spatial patterns of 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 (Jia et al., 2006; Keane et al., 2001; Rollins et al., 2004). Such information can improve our understanding of how trees and forests respond and will continue to respond towards changes in disturbance regimes (Swetnam, 1993) and help silviculturists design restoration treatments that move forests towards more resilient reference conditions. 5.2 Methods 5.2.1 Study area For this study, I focused on ten of the twelve research units in the Blacks Mountain Ecological Research Project (BMERP) and four research natural areas (RNA, Figure 5.1C). For the detailed description of RNA and BMERP refer to Chapter 2 (Section 2.1.1 and section 2.1.2). 85 Figure 5.1 Study area map: A) United States of America, B) California, C) Blacks Mountain Experimental Forest, and D) permanent grid system within Blacks Mountain Ecological Research project (BMERP). 86 5.2.2 Data 5.2.2.1 1934 census data Please refer to Chapter 2 (section 2.2.2) for the details about 1934 census data. 5.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 (Wing et al., 2012). The aircraft was flown at an altitude of 900 m using an opposing-flight line with a side lap of at least 50% (Wing et al., 2012). 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, Wing et al., 2012). The LiDAR point cloud was processed and LiDAR metrics were calculated using the ‘lidR’ package in R (Roussel et al., 2020, see Appendix B1 for a full description of the process). 5.2.2.3 2016 field data The ten research units and four RNAs each have a permanent 100 m lattice grid (Figure 5.1D; Oliver, 2000). The grid serves as the center points for most of the plot level research conducted in BMEF (Oliver, 2000). In the summer of 2016, at every other grid point in all diagonal directions (282 m spacing, Figure 5.1D), 16 m radius plots (804 m2, hereafter referred to as “circular plots”) were measured within each treatment unit (Mauro et al., 2019). Standing live and dead trees ≥ 9 cm diameter at breast height (DBH) were stem mapped from the plot center and measured for total height, DBH, crown width, and height to live and dead crown (Mauro et al., 2019). These ground-truthing data were collected in a total of 154 circular plots: 65 LOD, 69 HOD, and 20 RNA plots (Figure 5.1C; blue dots). 87 5.2.2.4 Overstory above-ground biomass for 1934 and 2016 Using the 1934 survey data, I calculated biomass (Mg) for foliage, branch, and bole of individual trees > 9cm DBH using species dependent equations that were developed locally at BMEF (Ritchie et al., 2013). 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. (1995). For white-fir and incense-cedar overstory AGB, I used equations and parameters suggested by Jenkins et al. (2003). The total plot overstory AGB was calculated by adding the overstory AGB of each species. I then converted plot AGB to AGB per hectare (Mg ha-1) based on the given plot sizes. For the 2016 field data, I 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. 5.2.2.5 Biomass model to link LiDAR metrics to 1 ha plots for 2016 I fit generalized linear mixed effect models (glm) with a‘gamma’ distribution to ensure that the expected value of AGB is always positive (Zuur et al., 2009). A ‘log link’ function was used to account for the to the distribution of the response variable—overstory above-ground biomass (AGB, Mg ha-1)—which was skewed to the right with values > 0 (Zuur et al., 2009, Appendix D: Figure D.1). Modelling was done in R using the ‘lme4’ package (Bates et al., 2015, see Appendix B.2 for the full details on model building process and results). In order to match the scale of ground-truthing plots, I divided the study area within BMERP into cells of size 32 x 32 m (Appendix D: Figure D.2). LiDAR metrics were extracted at 88 the cell level following the method described in Mauro et al. (2019). 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 (Appendix D: Figure D.2). 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 (Reed et al., 2016). 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-1) for 1 ha plots (see Appendix C.3. for details for the summary of prediction results). 5.3 Data analysis I computed semivariograms at the scale of the study units and Moran’s I correlogram as a measure of spatial patterns to evaluate my research objectives. 5.3.1 Comparing semivariogram components and Moran’s correlogram between 1934 and 2016 Semivariograms—To compare the spatial patterns of overstory AGB between reference and contemporary conditions, I used semivariogram models following the methods of Rossi et al. (1992) with the assumption of stationarity that the variance in above-ground biomass (AGB) is the function of separation distance only (Isaaks & Srivastava, 1989). 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 (Isaaks & Srivastava, 1989). Using the ‘gstat’ spatial package in R version 4.1.1 (Pebesma, 2004), semivariogram models were built for all fourteen units using 1 ha plot-level overstory AGB for 1934 and 2016. For the 2016 study units, I did not differentiate between the Burned and Unburned halves, hence variograms 89 represent the pooled variogram for each unit. Following the method of Fry & Stephens (2010), I fit three semivariogram models—exponential, Gaussian, and spherical—for each unit in 1934 and 2016. I also checked the assumption of stationarity in my exploratory analysis using directional variograms and found that ranges were not different, hence semivariance was only distance dependent and not direction dependent (Fry & Stephens, 2010). Furthermore, to allow for comparisons at a common scale, all variograms were standardized by dividing the semivariance by the overall sample variance (Rossi et al., 1992). For both 1934 and 2016, the best semivariogram models were selected based on the minimized root mean squared error (RMSE, Fry & Stephens, 2010). Following the method suggested by Fry & Stephens (2010), I 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 (Keane et al., 2012). The sill was compared to understand whether the amount of spatially dependent variance in overstory AGB differed between 1934 and 2016 (Fry & Stephens, 2010; Zawadzki et al., 2005). Nugget values were compared to understand if there was fine-scale spatial-variation (< 1 ha) in overstory AGB (Zawadzki et al., 2005). Moran’s I correlogram—For a more local measure of spatial autocorrelation, I followed the method of Jaquette et al. (2021) 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 my exploratory analysis, I found that Moran’s I could be calculated up to a distance of 1000 m for all units, because the number of observations was < 5% of total observations for the given unit and the spatial patterns could not be interpreted easily with so few 90 observations (Legendre & Fortin, 1989). RNA-D was an exception and only allowed 800m for the Moran’s I calculations. I also tested if Moran’s I at each lag was significantly different from 0 (alpha = 0.05) using Monte Carlo simulations with 1000 permutations (Jaquette et al., 2021). 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 (Legendre & Fortin, 1989; Radeloff et al., 2000) for all units. I used the Moran’s I correlograms for 1934 and 2016 to compare spatial patterns in overstory AGB in terms of the differences in: a. fine-scale spatial patterns—I 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 patterns within study units, i.e., patches of overstory AGB < 50 m do not exist (Legendre & Fortin, 1989). 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 (Legendre & Fortin, 1989; Radeloff et al., 2000). A greater number of peaks and troughs in the correlogram indicates heterogeneous spatial patterns, while fewer peaks in the correlogram suggest homogeneous spatial patterns (Legendre & Fortin, 1989; Radeloff et al., 2000). I visually compared if the number of peaks and troughs in the correlograms for the units were different between years as an indicator of heterogeneous or homogeneous spatial patterns (Legendre & Fortin, 1989; Radeloff et al., 2000). 91 5.3.2 Effect of management treatments on spatial patterns of above-ground biomass in the contemporary forests To compare the effect of management treatments on spatial patterns of overstory AGB, I used semivariogram models built for each Burned and Unburned halves of the HID and LOD treatment units (Appendix C: Table C.6). 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 visually compare and interpret the range, nugget, and sill. I 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 patterns and periodicity in spatial patterns in overstory AGB as described in section 2.3.1. From my exploratory analysis I found that maximum distance for which I could construct Moran’s I correlograms was 1000 m for all the Burned and Unburned halves of the units. 5.4 Results 5.4.1 Differences in spatial patterns between 1934 and 2016 5.4.1.1 Spatial autocorrelation in above-ground biomass exhibited at larger range in 2016 Within research natural areas (RNAs), Gaussian and spherical models for various units exhibited low RMSE, thus, were selected for comparison of spatial patterns in overstory above-ground biomass between 1934 and 2016 (Table 5.1). The models that did not converge for either year are not presented in the results. The values of range from both years indicated that spatial autocorrelation in overstory AGB was exhibited at shorter range in 1934 compared to 2016, without big differences in RNA-A and RNA-B (Figure 5.2 & Appendix D: Figure D.6). Exceptions were RNA-C and RNA-D, with larger range of ~400 m and 250 m, respectively, in 92 2016 compared to 1934 (~243 m for RNA-C and 167 m for RNA-D, Figure 5.2 & Appendix D: Figure D.6). RNAs exhibited larger nugget values in 1934 compared to extremely small values in 2016 indicating the presence of fine-scale spatial patterns in 1934 compared to 2016 (Figure 5.3 & Table 5.1). I only observed higher sill values in RNA-D in 2016 compared to 1934, without much difference in other RNAs (Table 5.1). Within units other than RNAs, spherical and Gaussian models exhibited low RMSE values and were selected (Table 5.1). All units except 41 and 44 exhibited short ranges in 1934 compared to longer ranges in 2016 (Figure 5.2, Appendix D: Figure D.7-8). There was a very small difference in range between both years in units 41, 43, 44, and 45 (Figure 5.2). The larger nugget values for all units in 1934 compared to 2016 indicated the presence of fine-scale spatial patterns (Figure 5.3). However, exceptions were units 40, 41, 42, 47 and 48, where the nugget values were extremely small in 1934 and no nugget was present in 2016 (Figure 5.3 & Table 5.1). The amount of spatial autocorrelation indicated by the sill values was greater in 1934 compared to 2016 (Table 5.1 & Appendix D: Figure D.6-8), exceptions were units 45, 38, and 41, where the sill was higher in 2016 compared to 1934 (Appendix D: Figure D.7-8). 5.4.1.2 Moran’s I correlogram exhibited differences in small scale spatial patterns and periodicity between two years Moran’s I correlogram suggested that fine-scale spatial patterns were 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 5.2). For example, in 1934, nine out of 14 units exhibited a Moran’s I value at the first lag that was not significantly different from 0, while in 2016 ten out of 14 units had Moran’s I values at the first lag that were significantly different from 0 (Table 5.2). This was in agreement with all the higher nugget values found in 1934 as 93 compared to low nugget values in 2016 (Figure 5.3). Thus, only four out of 14 units exhibited fine-scale spatial patterns in 2016 (Table 5.2). In general, periodicity did not vary substantially between 2016 and 1934 as indicated by similar numbers of peaks and troughs between both years (Table 5.2 & Appendix D: Figure D.9). However, the presence of periodicity in the Moran’s I correlograms in both years indicated that overstory AGB occurred in heterogeneous patches (Appendix D: Figure D.9). Exceptions were RNA-A and RNA-C, where the heterogenous spatial pattern in overstory AGB was more pronounced in 2016 as indicated by a larger number of peaks and troughs in 2016 compared to 1934 (Figure 5.4 & Table 5.2). 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 presence of bigger patches of overstory AGB (Appendix D: Figure D.9). 94 Table 5.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. Units Models RMSE Sill Research natural areas (RNA) 1934 2016 1934 2016 RNA-A Spherical 0.81 0.85 1.15 1.09 Exponential 0.82 0.84 1.91 1.75 Gaussian 0.74 0.80 1.01 0.99 RNA-B Spherical 0.83 0.96 1.21 1.28 Exponential 0.82 0.95 2.69 3.84 Gaussian 0.79 0.94 1.09 0.99 RNA-C Spherical 0.63 0.80 0.93 1.11 Exponential 0.72 0.82 1.01 1.49 RNA-D Gaussian 0.50 0.65 1.01 1.56 Low structural diversity (LOD) UNIT-39 Spherical 0.52 0.74 0.99 1.06 Exponential 0.57 0.74 1.94 1.92 Gaussian 0.49 0.70 0.95 0.98 UNIT-40 Spherical 0.76 0.73 1.04 1.11 Exponential 0.75 0.75 1.92 1.94 UNIT-43 Spherical 0.65 0.70 1.03 0.91 Exponential 0.69 0.71 1.16 1.06 Gaussian 0.63 0.67 0.98 0.87 UNIT-44 Spherical 0.30 0.51 1.42 0.96 Exponential 0.31 0.53 2.62 1.47 Gaussian 0.28 0.48 0.80 0.76 UNIT-45 Spherical 0.69 0.80 2.50 1.12 Gaussian 0.64 0.80 1.10 1.27 High structural diversity (HID) UNIT-38 Spherical 0.48 0.67 0.81 0.95 Exponential 0.52 0.69 0.85 1.07 Gaussian 0.47 0.64 0.79 0.89 UNIT-41 Spherical 0.72 0.71 1.03 0.94 Exponential 0.75 0.71 1.20 1.35 UNIT-42 Spherical 0.63 0.71 0.98 0.95 Exponential 0.86 0.84 1.54 2.52 UNIT-47 Spherical 0.85 0.84 1.14 1.18 Exponential 0.84 0.86 0.92 1.55 UNIT-48 Spherical 0.64 0.75 0.83 1.04 Exponential 0.66 0.75 0.00 0.00 95 Figure 5.2 Range values for above-ground biomass as indicated by selected variogram models between 1934 vs. 2016 within research natural area (RNA) units (top panel), low structural diversity (LOD) units (middle panel), and high structural diversity (HID) units (bottom panel) in Blacks Mountain Experimental Forest. 96 Figure 5.3 Nugget values for above-ground biomass as indicated by selected variogram models between 1934 vs. 2016 within research natural area (RNA) units (top panel), low structural diversity (LOD) units (middle panel), and high structural diversity (HID) units (bottom panel) in Blacks Mountain Experimental Forest. 97 Figure 5.4 Moran’s I correlogram for 1934 vs. 2016 within the four research natural areas (RNAs) in Blacks Mountain Experimental Forest. 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 (Jaquette et al. 2021). 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, while 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. 98 Table 5.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 1934 1934 2016 2016 1934 2016 #Peaks #Troughs #Peaks #Troughs RNA-A 1 4 1 2 p=0.009* p=0.009* RNA-B 2 1 2 1 p=0.23ns p=0.35ns RNA-C 2 3 4 2 p=0.13ns p=0.06* RNA-D 0 1 2 1 p=0.16ns p=0.47ns UNIT-39 1 4 1 2 p=0.009* p=0.009* UNIT-40 2 1 2 1 p=0.06ns p=0.02* UNIT-43 2 3 4 2 p=0.13ns p=0.009* UNIT-44 0 1 2 1 p=0.16ns p=0.03* UNIT-45 1 4 1 2 p=0.009* p=0.009* UNIT-38 2 1 2 1 p=0.009* p=0.009* UNIT-41 2 3 4 2 p=0.06ns p=0.17ns UNIT-42 0 1 2 1 p=0.72ns p=0.93ns UNIT-47 1 4 1 2 p=0.75ns p=0.05* UNIT-48 2 1 2 1 p=0.009* p=0.009* *indicates significant at alpha=0.05; ns=non-significant 5.4.2 2016 spatial patterns in treatment units 5.4.2.1 HID Burned halves exhibited spatial autocorrelation at short ranges Variograms models with lowest RMSE for both Burned and Unburned halves of HID and LOD treatments were selected for comparing range, nugget, and sill (Appendix C: Table C.4). 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.5, Appendix D: Figure D.10 & 11). All the Burned halves of HID showed larger nuggets compared to Unburned halves and most of the LOD treatments, indicating presence of fine-scale spatial patterns in overstory AGB within Burned halves (Figure 5.6). 99 Within LOD treatment units, most of the Burned halves exhibited spatial autocorrelation of overstory AGB at smaller range compared to Unburned halves, but the difference was very small (Figure 5.5). Exceptions were units 43 and 44, where spatial autocorrelation in overstory AGB was exhibited at shorter range in Unburned halves compared to Burned halves (Figure 5.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.5 & Appendix D: Figure D.10). The presence of fine-scale spatial patterns in overstory AGB were indicated by larger nugget values in Burned halves compared to Unburned halves, except for units 41 and 42, which exhibited no nuggets (Figure 5.6). 5.4.2.2 HID Burned halves exhibited heterogeneous spatial patterns 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 patterns at distances below 50 m (Table 5.3 & Appendix D: Figure D.12). I found a greater number of peaks and troughs of positive and negative Moran’s I values 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 (Appendix D: Figure D.12). For example, the Burned half of unit 38, showed four peaks and five troughs at the interval of 100-150 m lag-distances compared to only one peak in the Unburned half of unit 39, and two peaks in the Burned half of unit 39, all of which occurred at larger lag distances between 300 and 400 m (Table 5.3 & Appendix D: Figure D.12 A & F). 100 Correlograms for LOD treatments were more uniform with few statistically significant peaks (p = 0.0009) of Moran’s I at longer lag-distances compared to HID treatments (Appendix D: Figure D.12, LOD vs. HID). However, between Burned and Unburned halves of LOD treatment units, I observed a greater number of statistically significant (p = 0.0009) negative and positive Moran’s I values at short lag-distances for the Burned halves (Table 5.3 & Appendix D: Figure D.12, 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 and troughs of positive and negative Moran’s I values that were significantly different from 0 at shorter lag-distances (Appendix D: Figure D.12). Hence, both HID Burned and Unburned halves exhibited more spatially heterogeneous spatial patterns for overstory AGB compared to LOD treatments (Appendix D: Figure D.12, LOD vs. HID). 101 Figure 5.5 Differences in range for above-ground biomass as indicated by selected variogram models in between Burned and Unburned halves of low structural diversity (LOD) units (upper panel) and high structural diversity (HID) units (lower panel) in Blacks Mountain Experimental Forest. 102 Figure 5.6 Differences in sill for above-ground biomass as indicated by selected variogram models between Burned and Unburned halves within high structural diversity (HID) units (upper panel) and low structural diversity (LOD) units (lower panel) in Blacks Mountain Experimental Forest. 103 Table 5.3 Summary of the number of peaks and troughs of Moran’s I correlogram and p-values for the Moran’ I at first lag Burned and Unburned halves of low structural diversity treatment (LOD) and high structural diversity treatment (HID). Unit First lag Burned Burned Unburned Unburned Burned Unburned #Peak #Troughs #Peaks #Troughs LOD UNIT-39 2 1 1 1 p=0.009* p=0.009* UNIT-40 3 2 2 2 p=0.11ns p=0.009* UNIT-43 2 3 1 1 p=0.009* p=0.009* UNIT-44 1 3 1 1 p=0.005* p=0.009* UNIT-45 2 4 3 2 p=0.009* p=0.009* HID UNIT-38 4 4 2 1 p=0.009* p=0.002* UNIT-41 2 2 2 1 p=0.007* p=0.27ns UNIT-42 2 2 2 2 p=0.009* p=0.009* UNIT-47 3 2 2 1 p=0.009* p=0.009* UNIT-48 3 2 2 1 p=0.009* p=0.009* *indicates significant at alpha=0.05; ns=non-significant 5.5 Discussion Results showed that ranges, fine-scale variation, and periodicity in spatial patterns of overstory above-ground biomass (AGB) varied among units in BMEF. Forests of 1934 were heterogeneous even without marked human induced management. In contrast, by 2016 these forests had been subjected to some form of thinning, burning, or combination of both at various points in time (Crotteau & Ritchie, 2014; Oliver, 2000). Hence, spatial patterns in overstory AGB for 2016 are likely influenced by these treatments. My results showed that spatial patterns in 2016 varied among treatments and were more pronounced in prescribed burns. 5.5.1 The 1934 forest was spatially heterogeneous particularly at small-scales Most forest units had heterogenous spatial patterns in overstory AGB biomass in 1934 as indicated by the short ranges (<200 m) of the semivariograms, periodicity of Moran’s I 104 correlograms, and Moran’s I value at first lag not significantly different from 0 in RNAs (Legendre & Fortin, 1989; Radeloff et al., 2000). The heterogeneous spatial patterns in overstory AGB in BMEF are consistent with several other studies describing small patch size distributions (ranging from 60-150 m) of the overstory in old growth ponderosa pine forests in northern California (Arthur & Janet, 1987; Youngblood et al., 2004), and Washington (Harrod et al., 1998).s Further, my results identifying fine-scale variation in most of the units were consistent with the study from Freeman & Moisen (2007), who attributed the presence of fine-scale variation in biomass to larger nugget values. Exceptions occurred in units such as RNA-A, 45, and 48, which lacked fine-scale (<50 m) spatial patterns. The absence of fine-scale spatial patterns is likely due to a lack of frequent fire that ended around the late 1880s at BMEF (unpublished data from Carl Skinner). In the absence of frequent fire, stand density possibly increased with the infilling of gaps with shade-tolerant species such as incense-cedar and white fir (Youngblood et al., 2004) in units such as RNA-A, 45, and 48, which occurred at higher elevations with abundant shade-tolerant species. 5.5.2 Larger and more homogenous spatial patterns in 2016 except when burned The long ranges (> 250 m) of variograms suggested that overstory AGB in most units in 2016 consisted of larger more homogeneous spatial patterns, possibly due to increases in tree densities within 1934 forest gaps (Lydersen et al., 2013; North et al., 2004). With the 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 (Lydersen et al., 2013). 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, which was consistent with Freeman & Mois 105 (2007). On the other hand, no difference in periodicity between 1934 and 2016 within prescribed burn RNA units indicates that at localized scales, prescribed fire likely emulated much of the 1934 spatial structure by killing trees in patches (Kerns et al., 2006). 5.5.3 Prescribed burn enhanced spatial heterogeneity in both LOD and HID treatments The high structural diversity (HOD) treatment for BMEF was designed to approximately emulate a heterogeneous overstory stand structure, with numerous small openings and patches of large old trees (Crotteau & Ritchie, 2014). Therefore, my results regarding overstory spatial patterns in AGB—presence of fine-scale patterns, periodicity in patchiness at short lags, and short ranges—for HID treatments suggest that treatments were effective in creating the heterogeneous spatial patterns of overstory stand structure (Radeloff et al., 2000). Furthermore, when prescribed burns were coupled with the HID treatment, spatial heterogeneity in overstory AGB biomass 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 kill small trees in groups and create openings (Kerns et al., 2006; Sackett et al., 1996). In addition, second-order post-fire mortality of trees can also remove patches of trees (Sackett et al., 1996, Hood et al., 2010). For example, low severity prescribed burns can leave fire-injured conifers that are receptive hosts for bark beetles (Hood et al., 2010). Subsequent mortality of weakened trees, especially large old ponderosa pine - due to beetle attack can produce openings and gaps within stands (Hood et al., 2010, Davis et al., 2012); therefore, creating a heterogeneous overstory similar to our HID Burned halves. My results in HID Burned units are consistent with Dodson & Peterson (2010), who found that thinning increased spatial aggregation of residual 106 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 few large gaps (Crotteau & Ritchie, 2014). The largest and smallest trees were removed from these stands for a unimodal diameter distribution only leaving intermediate-sized trees (Oliver, 2000). Especially within the Unburned halves of the treatments, the long ranges from the semivariograms, and 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. My results from Unburned halves of the LOD units are consistent with Larson et al. (2012), who found that after uniform removal of trees, mixed-conifer stands showed homogeneous spatial patterns at broad scales. Similar homogeneous overstory spatial patterns in spruce-fir forests were found by Kuehne et al. (2018), 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. 5.5.4 Management implications My findings about enhanced spatial heterogeneity after post-treatment prescribed burns are similar to what other studies have reported for dry forests across California (e.g., Churchill et al., 2013; Knapp et al., 2017). My comparison of spatial patterns 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 patterns of overstory AGB that I detected in the 1934 reference forests helps portray the types of patterns 107 associated with very little 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 semivariogram ranges across time, managers can understand the distance (or scale) over which spatial patterns 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 patterns. In addition, comparison of reference and contemporary overstory spatial patterns can help managers to explore the impact of fire exclusion on driving overstory spatial patterns in dry forests (Jaquette et al., 2021). Furthermore, my 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 patterns of reference conditions. For example, in managing contemporary forests, treatments similar to HID with prescribed burns that preserve and/or create patchy and heterogeneous spatial patterns similar to reference conditions, may be warranted. In addition, comparing changes in the spatial patterns of the contemporary forests will allow managers to understand the effectiveness of the restoration treatments—such as HID with prescribed burns—in creating desirable spatial patterns for contemporary dry forests. 108 5.6 Conclusion This study provides a unique opportunity to examine the spatial patterns 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, I was able to fill a gap in the literature related to understanding the spatial pattern of overstory above-ground biomass between contemporary and reference forests. The spatial patterns of above-ground biomass in reference conditions described in this study come from forests that have minimum human impact following Euro-American settlement, and hence represent a relatively intact dry forest ecosystem. Therefore, knowledge about spatial patterns of such intact forests in reference conditions provides insights on 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. My results further indicate that low severity fire seems to be key important 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. 109 Chapter 6: Conclusion 6.1 Overall summary and contributions Quantification of reference conditions can help enable identification of contemporary forests potentially at-risk of insects, disease, and fire and then allow managers to design restoration treatments that better emulate the processes that help create forest stand structures resilient to fire, insects, and disease. However, a near universal lack of spatially explicit stand structure data representing overstory reference conditions makes it challenging to understand long-term changes in spatial patterns in forests using contemporary spatial datasets—such as LiDAR. This lack of reference data is particularly problematic for managers of ponderosa pine dominated forests (or “dry forests”) who wish to understand forest reference conditions and use such baseline conditions to guide management. In this dissertation, I utilized distinctive census data from reference conditions (1934) collected in Blacks Mountain Experimental Forest (BMEF), in California, USA, over an area of ~ 4,000 ha to quantify the spatial patterns of overstory stand structure in historical reference conditions. More specifically, I first used stand density index (SDI) for tree species across seven vigor classes to identify overstory spatial patterns in forest reference conditions using 1934 census data (Chapter 3). Spatial clusters consisting of low tree density and high tree vigor within in 1934 reference forests indicated that forest conditions consisted of open stands with healthy (high vigor) trees. More specifically, my chapter 3 provided information on two components of forest resiliency in reference conditions: a) resiliency against insects, and b) resiliency against fire. For example, in 1934, BMEF consisted of more than 50% ponderosa pine high vigor forest: young, dominant and/or co-dominant and well-developed trees in vigor classes 1, 2, and 6 trees. Based on Dunning’s (1928) observations these trees were mostly insect and disease resistant and 110 high vigor, which made BMEF a resilient forest against both insects and disease. Furthermore, each cluster I identified had an ample amount of vigor class 5 trees, which were over mature (> 300 years old), thick barked, fire resistant, and dominant trees (Dunning, 1928). The presence of vigor class 5 trees throughout BMEF indicated that the forest structure in BMEF consisted substantially of fire-resistant old growth trees similar to the forests of the pre-settlement era. Because dendrochronological information cannot be used to describe a detailed spatial patterns of tree vigor in historical reference conditions, my results contribute to filling gaps in our knowledge regarding an empirical baseline for forest managers to identify forests that are at-risk of insects, disease, and fire. In Chapter 4, I utilized reference and contemporary regeneration data based on three species and four height classes from BMEF to quantify the differences in regeneration between reference and contemporary forest conditions. In addition, I also linked regeneration in reference forests with environmental and soil variables. Because natural tree regeneration is known to be highly spatially heterogeneous, direct measures of regeneration density can be misleading. Thus, I used stocking—a more expressive measure of assessing tree occupancy—which accounts for the spatial heterogeneity of natural regeneration. Not only did I find that stocking was dominated by ponderosa pine seedlings in reference conditions, it was also associated with available water capacity at 150 cm soil depth and overstory basal area. Results from Chapter 4 are critically important to help managers understand how regeneration patterns can be considered in forest management treatments. My research utilizes field collected regeneration data from 1934, which allowed me to depict an accurate and detailed picture of regeneration. In contrast, existing reconstruction studies fail to capture information such as seedlings established post-fire. 111 In Chapter 5, I used contemporary (2016) remotely sensed LiDAR data in conjunction with the 1934 census data specific to the Blacks Mountain Ecological Research Project (BMERP) to quantify the differences in spatial patterns in above-ground biomass (AGB) between reference and contemporary forest conditions. Differences in overstory spatial patterns are among the less well studied aspects of dry forests. Studies investigating spatial patterns of tree density and size (e.g., Sanchez Meador et al., 2010; Youngblood et al., 2004), have generally ignored overstory components such as above-ground biomass, which can provide important information on forest fuels. Thus, my study is one of the first to report spatial patterns of reference conditions using AGB (rather than tree density and size). My study of reference conditions pertains to a forest that was least impacted by human management. Thus, my comparison of spatial patterns between reference and contemporary conditions contributes methods suitable for identifying reference points against which departures from spatial patterns can be assessed. In addition, because many management units I examined received various forms of thinning and prescribed burn treatments, my results can serve as a guide for measuring the success of restoration treatments aimed at emulating the spatial patterns of reference conditions. 6.2 Limitations 6.2.1 Reference conditions spatial pattern chapter limitations The study in Chapter 3 comes with a few inherent limitations. First, the study does not provide insights into fine-scale spatial patterns of less than 1 ha due to the design of the 1934 census data collection. Spatial patterns at scale smaller than 1 ha can be important from a wildlife perspective (Kolb et al., 2007). Secondly, I excluded trees < 8.9 cm DBH, which leaves out the recent regeneration that was established around the early 1900s. Regeneration established around 1900s likely influenced tree vigor due to competition of understory and overstory for resources (Riegel 112 et al., 1992). Even within areas of BMEF where fire regimes were still intact in 1934, due to the patchy nature of fire, small trees (excluded from my study) must have escaped fire and their establishment could have influenced the spatial patterns of reference conditions (Taylor, 2010). In addition, standing dead trees possibly killed by fire were also not included in the 1934 data. For 1934 data, there was also uncertainty associated with the number of trees that were established after the start of fire suppression in 1905, which makes it difficult to associate the effect of fire suppression on spatial patterns (Scholl & Taylor, 2010; Taylor, 2010). Therefore, my results from Chapter 3 potentially provide an incomplete picture of spatial patterns of reference conditions in dry forests with frequent fires. There were some uncertainties associated with soil variables. For example, BMEF is a highly rocky area, which likely resulted in soil samples being restricted to shallower soil depth and a scenario where soil had to be resampled at a new location, which possibly could provide inaccurate classification of the soil types, texture, and available water capacity (D. Young, personal communication, June 21, 2023). 6.2.2 Regeneration chapter limitations Most of the limitations for the regeneration study presented in Chapter 4 came from the variables that were used to fit the beta regression models. Soil variables from the available EUI map were produced in 1994, decades after 1934. Therefore, my analysis assumes that the 1994 soil variables represent the soil conditions that existed in BMEF around 1934 and that there was no change in soil conditions between 1934 and 1994. Secondly, the fitted models were unable to explain all existing variability in stocking. The unexplained variability in stocking could be explained by incorporating fire history variables such as time since fire, fire frequency, and fire return interval into the model. However, availability of fire history variables was limited to 28 of the 100 compartments. The fire history variables were statistically significant when I added them 113 to the model using only 28 compartments. However, due to the small sample size of 28, the model was overfit, and I did not think it was meaningful to use them as auxiliary variables in my model. Yet, this exploratory analysis confirms that fire history variables could explain some of the observed variability in stocking. 6.2.3 Reference vs. contemporary spatial patterns chapter limitation In Chapter 5, limitations result from the scale at which I studied the spatial patterns. Because 1934 overstory data were collected at the scale of 1 ha, I could not utilize the full potential of LiDAR data in 2016 to derive the spatial patterns at the individual tree level. It is important to understand the tree level spatial patterns such as microsite conditions influence differences in competitive status of the trees (Britwum et al., 2023), which will highly alter the tree spatial patterns at different scales (Perry et al., 2009). Another limitation includes the limited availability of ground-truthing data in 2016 outside the treatment BMERP units, rendering it difficult to understand the spatial patterns of the overstory stand structure throughout Blacks Mountain Experimental Forest. Chapters 3 and 4 suggested that stand structure and spatial patterns varied based on elevation and soil gradient. However, I excluded large areas of BMEF from my study. The excluded areas can have variation in elevation and soil, which could influence the spatial patterns, therefore including those areas might provide a different spatial pattens than what I have depicted in my study. Hence, the spatial patterns in overstory above-ground biomass that I have discerned may be quite different than what currently exists if the entire scale of BMEF is taken into consideration. Therefore, the results from my study should be taken with caution and can only be extrapolated to ponderosa pine dominated ecosystems that received some form of combination treatments—thinning, prescribed burn, or a combination of both. 114 6.3 Future directions and recommendations 6.3.1 Utilize fire data to improve the understanding of spatial patterns in reference conditions In dry forests of western North America, studies of spatial patterns have led to an improved understanding of stand dynamics. Studies have revealed that relatively frequent, low- to mixed-severity fires in reference conditions shaped stand structure by creating and maintaining generally open, uneven-aged stands consisting of a mosaic of individual trees, tree groups, and openings (Larson & Churchill, 2012). The forest mosaic in dry forests have been characterized as aggregated at smaller scales (0.004 ha) with heterogeneous within group spatial patterns (Larson et al., 2012; Pearson, 1950). One of the important drivers of these mosaic aggregates is fire as fires reinforced these mosaic spatial patterns (Larson & Churchill, 2012). Elements of forest dynamics such as demography (Boyden et al., 2005), mortality (Boyden et al., 2005), tree growth, and regeneration (Fajardo et al., 2007) are regulated by fire. Therefore, one of the important drivers of spatial patterns to be included in future studies can be fire, as it effectively creates a heterogeneous mosaic of forest stands at different sizes and scales (Larson & Churchill, 2012); such heterogeneity can favor the resiliency of forest stands. 6.3.2 Use of spatially explicit regeneration data and climate variables My study on regeneration (Chapter 4) has provided some insights on the potential roles of topography and tree overstory in influencing regeneration. However, several research studies have suggested that the rate and pattern of tree regeneration in ponderosa pine forests is influenced by a number of complex interacting spatial phenomena such as the spatial variability in solar radiation, temperature, precipitation, and the behavior of wildlife (Bonnet et al., 2005; Chambers et al., 2016). My study did account for the environmental covariates, but lack of direct 115 measure of climate variables and spatial patterns of regeneration limits my ability to discern the exact mechanisms controlling regeneration patterns. In addition, ponderosa pine recruitment and establishment is highly dependent on inter-annual climate (Shepperd et al., 2006). Specifically, changes in moisture regimes following the fire event could lead to changes between facilitation and competition among species (Brown et al., 2015; Pielou, 1962). Therefore, one of the avenues to explore in future studies would be utilizing climate variables and spatially explicit regeneration data. Including climate variables and spatial information would result in increased understanding of regeneration patterns at smaller scales than what I have observed at compartment level in my study. 6.3.3 Fine scale spatial patterns and use of environmental covariates Spatial patterns of forest structure and their development are imprinted by ecological processes such as disturbances and management (Ziegler et al., 2017). Ecological processes also influence the spatial patterns of forest stand structure at multiple scales—more finer scales than coarse—and time (Hessburg et al., 1999). Thus, multiscale assessment of spatial pattern change is essential to understand change in forest ecosystem dynamics and function (Hessburg et al., 1999). As growth and structure of trees vary depending upon the microsite conditions, utilizing data from variables such as topographic wetness index, slope, aspect, and soil taken at scale < 1 ha can be included in discerning the spatial patterns in forest stand structure (Dungey et al., 2018). At places where data are absent on direct measures of environmental factors, an alternative approach would be to include variables such as presence of stumps and coarse woody material (Chambers et al., 2016). Coarse woody material can be used as proxy that could influence patterns of moisture and soil and hence fine scale spatial patterns. At a broader scale, one of the promising avenues to explore is the use of overstory variables and associated 116 environmental variables at finer scale, which will allow mangers to understand how small scale (< 1 ha) environmental conditions affect spatial patterns of overstory stand structure. Future research that aims to understand the influence of environmental conditions on spatial patterns should also look across scales to improve the understanding on how stand development occurs at different scales (Ziegler et al., 2017). 6.4 Resiliency comes from heterogeneity Disturbance regimes for dry forests in Western North America are continuously changing, especially exacerbated by interaction of land management and climate change (Fettig et al., 2013). An emerging goal for ecologists and forest managers is to build resilient dry forests for current and future climates in the event of climate change and catastrophic disturbances such as wildfire and landscape-scale insect and disease outbreaks (Miller et al., 2007). While thinning and burning can be common prescriptions, goals to promote resiliency may, however, require an in depth understanding of interaction of overstory components such as tree density, tree vigor, and regeneration in space and time. Such understanding can refine current restoration treatments, which eventually will result in heterogeneous overstory spatial patterns that can withstand and recover after any disturbance event. 117 Bibliography Abella, S. R., & Denton, C. W. (2009). Spatial variation in reference conditions: Historical tree density and pattern on a pinus ponderosa landscape. Canadian Journal of Forest Research, 39(12), 2391–2403. https://doi.org/10.1139/X09-146 Agee, J. K., & Skinner, C. N. (2005). Basic principles of forest fuel reduction treatments. Forest Ecology and Management, 211(1–2), 83–96. https://doi.org/10.1016/j.foreco.2005.01.034 Alexander, E. B. (1994). Experimental forests and ranges of the Pacific Southwest Research Station. Unpublished Report on File at the Pacific Southwest Research Station, USDA Forest Service, Redding, California. http://www.cabdirect.org/abstracts/19930670177.html Alexander, E. B. (1988). Morphology, fertility and classification of productive soils on serpentinized peridotite in California (U.S.A.). Geoderma, 41(3-4), 337-351.https://doi.org/10.1016/0016-7061(88)90069-9. Akaike, H. (1973). Information theory and the maximum likelihood principle. In: Petrov, B.N., Csaki, F. (Eds.), Proceedings of the 2nd International Symposium on Information Theory September 2-8, 1971; Budapest, Hungary, pp. 267-281.https://doi.org/10.1007/978-1-4612-1694-0_15 Allen, C. D., Savage, M., Falk, D. A., Suckling, K. F., Swetnam, T. W., Schulke, T., Stacey, P. B., Morgan, P., Hoffman, M., & Klingel, J. T. (2002). Ecological restoration of southwestern ponderosa pine ecosystems: A broad perspective. Ecological Applications, 12(5), 1418–1433. https://doi.org/10.1890/1051-0761(2002)012[1418:EROSPP]2.0.CO;2 Angell, N., Waring, K. M., & Graves, T. A. (2014). Predicting height growth of sugar pine regeneration using stand and individual tree characteristics. Forestry, 87(1), 85–98. https://doi.org/10.1093/forestry/cpt028 118 Ansley, J. S., & Battles, J. J. (1998). Forest composition, structure, and change in an old-growth mixed confer forest in the northern Sierra Nevada. Journal of the Torrey Botanical Society, 125(4), 297–308. https://doi.org/10.2307/2997243 Bailey, J. D., & Covington, W. W. (2002). Evaluating ponderosa pine regeneration rates following ecological restoration treatments in northern Arizona, USA. Forest Ecology and Management, 155(1–3), 271–278. https://doi.org/10.1016/S0378-1127(01)00564-3 Bailey, J., & Idle, D. (2001). Four corners regional forest assessment. School of forestry report. Northern Arizona University, Flagstaff, AZ. Baker, W. L., & Williams, M. A. (2015). Bet-hedging dry-forest resilience to climate-change threats in the western USA based on historical forest structure. Frontiers in Ecology and Evolution, 2, 1–7. https://doi.org/10.3389/fevo.2014.00088 Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1 SE-Articles), 1–48. https://doi.org/10.18637/jss.v067.i01 Beaty, R. M., & Taylor, A. H. (2008). Fire history and the structure and dynamics of a mixed conifer forest landscape in the northern Sierra Nevada, Lake Tahoe Basin, California, USA. Forest Ecology and Management, 255(3–4), 707–719. https://doi.org/10.1016/j.foreco.2007.09.044 Behrens, T., Viscarra Rossel, R. A., Kerry, R., MacMillan, R., Schmidt, K., Lee, J., Scholten, T., & Zhu, A. X. (2019). The relevant range of scales for multi-scale contextual spatial modelling. Scientific Reports, 9(1), 1–9. https://doi.org/10.1038/s41598-019-51395-3 Beven, K. J., & Kirkby, M. J. (1979). A physically based, variable contributing area model of basin hydrology. Hydrological Sciences Bulletin, 24(1), 43–69. 119 https://doi.org/10.1080/02626667909491834 Bonnet, V. H., Schoettle, A. W., & Shepperd, W. D. (2005). Postfire environmental conditions influence the spatial pattern of regeneration for Pinus ponderosa. Canadian Journal of Forest Research, 35(1), 37–47. https://doi.org/10.1139/x04-157 Bonnicksen, T. M. ., & Stone, E. C. (1982). Reconstrution of a presettlement giant sequoia-mixed conifer forest community using the aggregation approach. Ecology, 63(4), 1134–1148. https://doi.org/10.2307/1937250 Bolsinger, C. L. (1988). The hardwoods of California''s timberlands, woodlands, and savannas. Res. Bull. PNW-RB-148. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station. 148 p. https://doi.org/10.2737/PNW-RB-148 Boyden, S., Binkley, D., & Shepperd, W. (2005). Spatial and temporal patterns in structure, regeneration, and mortality of an old-growth ponderosa pine forest in the Colorado Front Range. Forest Ecology and Management, 219(1), 43–55. https://doi.org/10.1016/j.foreco.2005.08.041 Breiman, L., Friedman, H. J., Olshen, A. R., & Stone, J. C. (1984). Classification and regression trees. CHAPMAN & HALL/CRC. Boca Raton, New York. Britwum, S., Marshall, P. L., Eskelson, B. N. I., & Barbeito, I. (2023). Temporal changes in tree spatial patterns in uneven-aged interior Douglas-fir dominated stands managed under different thinning treatments. Forest Ecology and Management, 528, 120640. https://doi.org/10.1016/j.foreco.2022.120640 Brown, P. M., Battaglia, M. A., Fornwalt, P. J., Gannon, B., Huckaby, L. S., Julian, C., & Cheng, A. S. (2015). Historical (1860) forest structure in ponderosa pine forests of the northern front range, Colorado. Canadian Journal of Forest Research, 45(11), 1462–1473. 120 https://doi.org/10.1139/cjfr-2014-0387 Brown, P. M., & Sieg, C. H. (1996). Fire history in interior ponderosa pine communities of the Black Hills, South Dakota, USA. International Journal of Wildland Fire, 6(3), 97–105. https://doi.org/10.1071/WF9960097 Camann, M. A., Gillette, N. E., Lamoncha, K. L., & Mori, S. R. (2008). Response of forest soil Acari to prescribed fire following stand structure manipulation in the southern Cascade Range. Canadian Journal of Forest Research, 38(5), 956–968. https://doi.org/10.1139/X07-241 Chambers, M. E., Fornwalt, P. J., Malone, S. L., & Battaglia, M. A. (2016). Patterns of conifer regeneration following high severity wildfire in ponderosa pine – dominated forests of the Colorado Front Range q. Forest Ecology and Management, 378, 57–67. https://doi.org/10.1016/j.foreco.2016.07.001 Chen, J., Song, B., Rudnicki, M., Moeur, M., Bible, K., North, M., Shaw, D. C., Franklin, J. F., & Braun, D. M. (2004). Spatial relationship of biomass and species distribution in an old-growth. Forest Science, 50(3), 364–375. 10.1371/journal.pone.0186742 Churchill, D. J., Carnwath, G. C., Larson, A. J., & Jeronimo, S. A. (2017). Historical forest structure, composition, and spatial pattern in dry conifer forests of the western Blue Mountains, Oregon. Technical Report, U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station, Portland, OR: PNW-GTR-95(November), 1–93. Churchill, D. J., Larson, A. J., Dahlgreen, M. C., Franklin, J. F., Hessburg, P. F., & Lutz, J. A. (2013). Restoring forest resilience: From reference spatial patterns to silvicultural prescriptions and monitoring. Forest Ecology and Management, 291, 442–457. https://doi.org/10.1016/j.foreco.2012.11.007 121 Churchill, D. J., Seager, S. T., Larson, A. J., Schneider, E. E., Kemp, K. B., & Bienz, C. (2018). Ecological functions of spatial pattern in dry forests: implications for forest restoration. The Nature Conservancy Newsleter, January. Clark, J. S., Beckage, B., Camill, P., Cleveland, B., HilleRisLambers, J., Lichter, J., McLachlan, J., Mohan, J., & Wyckoff, P. (1999). Interpreting recruitment limitation in forests. American Journal of Botany, 86(1), 1–16. https://doi.org/10.2307/2656950 Clyatt, K. A., Keyes, C. R., & Hood, S. M. (2017). Long-term effects of fuel treatments on aboveground biomass accumulation in ponderosa pine forests of the northern Rocky Mountains. Forest Ecology and Management, 400, 587–599. https://doi.org/10.1016/j.foreco.2017.06.021 Cochran, P. H., Geist, J. M., Clemens, D. L., Clausnitzer, R. R., & Powell, D. C. (1994). Suggesting stocking levels for forest stands in northeastern Oregon and southeastern Washington. Research Paper - US Department of Agriculture, Forest Service, PNW-RN-513. Collins, B. M., Fry, D. L., Lydersen, J. M., Everett, R., & Stephens, S. L. (2017). Impacts of different land management histories on forest change. Ecological Applications, 27(8), 2475–2486. https://doi.org/10.1002/eap.1622 Collins, B. M., Bernal, A., York, R. A., Stevens, J. T., Juska, A., & Stephens, S. L. (2021). Mixed-conifer forest reference conditions for privately owned timberland in the southern Cascade Range. Ecological Applications, 31(7), 1–14. https://doi.org/10.1002/eap.2400 Cooper, C. F. (1960). Changes in vegetation, structure, and growth of southwestern pine forests since white settlement. Ecological Monographs, 30(2), 129–164. https://doi.org/10.2307/1948549 122 Cooper, C. F. (1961). Pattern in Ponderosa Pine Forests. Ecology, 42(3), 493–499. https://doi.org/10.2307/1932235 core Team, R. (2018). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Crotteau, J. S., & Ritchie, M. W. (2014). Long-term stand growth of interior ponderosa pine stands in response to structural modifications and burning treatments in northeastern California. Journal of Forestry, 112(5), 412–423. https://doi.org/10.5849/jof.13-090 Curtis, R.O. (2010). Effect of diameter limits and stand structure on relative density indices: A case study. Western Journal of Applied Forestry, 25(4):169-175. https://doi.org/10.1093/wjaf/25.4.169. Daubenmire, R. (1968). Soil moisture in relation to vegetation distribution in the mountains of northern Idaho. Ecology, 49(3), 431–438. Davis, R. S., Hood, S. M., & Bentz, B. J. (2012). Fire-injured ponderosa pine provide a pulsed resource for bark beetles. Canadian Journal of Forest Research, 42(12), 2022-2036. Dodson, Erich K., & Peterson, D. W. (2010). Dry coniferous forest restoration and understory plant diversity: The importance of community heterogeneity and the scale of observation. Forest Ecology and Management, 260(10), 1702–1707. https://doi.org/10.1016/j.foreco.2010.08.012 Dodson, E. K., & Root, H. T. (2013). Conifer regeneration following stand-replacing wildfire varies along an elevation gradient in a ponderosa pine forest, Oregon, USA. Forest Ecology and Management, 302, 163–170. https://doi.org/10.1016/j.foreco.2013.03.050 Dolanc, C. R., Safford, H. D., Thorne, J. H., & Dobrowski, S. Z. (2014). Changing forest structure across the landscape of the Sierra Nevada, CA, USA, since the 1930s. Ecosphere, 123 5(8), 1–26. https://doi.org/10.1890/ES14-00103.1 Dubayah, R., & Rich, P. M. (1995). Topographic solar radiation models for GIS. International Journal of Geographical Information Systems, 9(4), 405–419. https://doi.org/10.1080/02693799508902046 Dunning, D. (1928). A tree classification for the selection forests of the Sierra Nevada. Journal of Agricultural Research, 36(9), 755–771. Eaton, C. B. (1959). Insect-caused maortality in relation to methods of cutting in ponderosa pine on Blacks Mountain Experimental Forest. US Department of Agriculture, Forest Service, PSW-GTR-43, 36p. Everett, R. G. (2008). Dendrochronology-based fire history of mixed-conifer forests in the San Jacinto Mountains, California. Forest Ecology and Management, 256(11), 1805–1814. https://doi.org/10.1016/j.foreco.2008.04.036 Fajardo, A., Graham, J. M., Goodburn, J. M., & Fiedler, C. E. (2007). Ten-year responses of ponderosa pine growth, vigor, and recruitment to restoration treatments in the Bitterroot Mountains, Montana, USA. Forest Ecology and Management, 243(1), 50–60. https://doi.org/10.1016/j.foreco.2007.02.006 Farnden, C. (1994). Forest regeneration in the ESSF zone of north-central British Columbia. Canadian Forest Service, Pac. For. Res. Centre, Information Report BC-X-351, 31p. Ferrari, S. L. P., & Cribari-Neto, F. (2004). Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799–815. https://doi.org/10.1080/0266476042000214501 Ferster, C. J., Eskelson, B. N. I., Andison, D. W., & LeMay, V. M. (2016). Vegetation mortality within natural wildfire events in the Western Canadian boreal forest: What burns and why? 124 Forests, 7(9), 1–20. https://doi.org/10.3390/f7090187 Fettig, C. J., Reid, M. L., Bentz, B. J., Sevanto, S., Spittlehouse, D. L., & Wang, T. (2013). Changing climates, changing forests: A Western North American perspective. Journal of Forestry, 111, 214–228. https://doi.org/10.5849/jof.12-085 Finley, K., & Zhang, J. (2019). Climate effect on ponderosa pine radial growth varies with tree density and shrub removal. Forests, 10(6). https://doi.org/10.3390/f10060477 Fites-Kaufman, J., Rundel, P. W., Stephenson, N. L., & Weixelman, D. A. (2007). Montane and subalpine vegetation of the Sierra Nevada and Cascade Ranges. In: Barbour, M.G., Major, J. (Eds.), Terrestrial Vegetation of California. University of California Press, Berkeley, CA, USA, pp. 456–501. Fitzgerald, S. (2005). Fire ecology of ponderosa pine and the rebuilding of fire-resilient ponderosa pine. USDA Forest Service General Technical Report, PSW-GTR-198, 28p. Flathers, K. N., Kolb, T. E., Bradford, J. B., Waring, K. M., & Moser, W. K. (2016). Long-term thinning alters ponderosa pine reproduction in northern Arizona. Forest Ecology and Management, 374, 154–165. https://doi.org/10.1016/j.foreco.2016.04.053 Fortin, M. J., & Dale, M. R. T. (2005). Spatial analysis: A guide for ecologists. Cambridge University Press. https://doi.org/DOI: 10.1017/CBO9780511542039 Foster, D., Swanson, F., Aber, J., Burke, I., Brokaw, N., Tilman, D., & Knapp, A. (2003). The importance of land-use legacies to ecology and conservation. Foster_Bioscience_2003. 53(1).77-88. https://doi.org/10.1641/0006-3568(2003)053[0077:TIOLUL]2.0.CO;2 Fowells, H. A., & Stark, N. B. (1965). Natural regeneration in relation to environment in the mixed conifer forest type of California. U.S. Forest Service Research Paper PSW - 24, 14p. Franklin, J. F., & Pelt, R. Van. (2004). Spatial aspects of structural complexity in old-growth 125 forests. Journal of Forestry, May, 22–28. https://www.jstor.org/stable/40925512 Frayer, W. E., & Furnival, G. M. (1998). History of Forest Surveys Sampling Designs in the United States. Gen. Tech. Rep. NC-212. St. Paul, MN: U.S. Dept. of Agriculture, Forest Service, North Central Forest Experiment Station. 42-49. Freeman, E. A., & Moisen, G. G. (2007). Evaluating kriging as a tool to improve moderate resolution maps of forest biomass. Environmental Monitoring and Assessment, 128(1–3), 395–410. https://doi.org/10.1007/s10661-006-9322-6 Fry, D. L., & Stephens, S. L. (2010). Stand-level spatial dependence in an old-growth Jeffrey pine -mixed conifer forest, Sierra San Pedro Mártir, Mexico. Canadian Journal of Forest Research, 40(9), 1803–1814. https://doi.org/10.1139/X10-122 Fry, D. L., Stephens, S. L., Collins, B. M., North, M. P., Franco-Vizcaíno, E., & Gill, S. J. (2014). Contrasting spatial patterns in active-fire and fire-suppressed mediterranean climate old-growth mixed conifer forests. PLoS ONE, 9(2). https://doi.org/10.1371/journal.pone.0088985 Fulé, P. Z., Crouse, J. E., Roccaforte, J. P., & Kalies, E. L. (2012). Do thinning and/or burning treatments in western USA ponderosa or Jeffrey pine-dominated forests help restore natural fire behavior? Forest Ecology and Management, 269, 68–81. https://doi.org/10.1016/j.foreco.2011.12.025 Arthur, G., & Janet, F. (1987). Second-order neighborhood analysis of mapped point patterns Ecology, 68(3), 473–477. http://www.jstor.org/stable/1938452 Getzin, S., Dean, C., He, F., Trofymow, J. A., Wiegand, Ke., & Wiegand, T. (2006). Spatial patterns and competition of tree species in a Douglas‐fir chronosequence on Vancouver. Ecograph, 29, 671–682. https://doi.org/10.1111/j.2006.0906-7590.04675.x 126 Ghent, A. W. (1969). Studies of regeneration in forest stands devastated by the spruce budworm. Forest Science, 15(4), 417–429. https://doi.org/10.1093/forestscience/15.4.417 Gill, C. B. (1950). A study of the relationship between the number of trees per acre and dispersion. The Forestry Chronicle, 26(3), 186–196. https://doi.org/10.5558/tfc26186-3 Ginrich, S. (1967). Measuring and evaluating stocking and stand density in upland hardwood forests in the central states. Forest Science, 13(1), 38–53. https://doi.org/10.1093/forestscience/13.1.38 Gordon, T. D. (1981). History of the Blacks Mountain Experimental Forest, 1933 through 1981. Pacific Southwest Research Station, USDA Forest Service, PSW-GTR, 1–209. Graham, R. T., & Jain, T. B. (2005). Ponderosa pine ecosystems. General Technical Report - Pacific Southwest Research Station, USDA Forest Service, PSW-GTR-198, 1–32. Grulke, N., Maxfield, J., Riggan, P., & Schrader-Patton, C. (2020). Pre-emptive detection of mature pine drought stress using multispectral aerial imagery. Remote Sensing, 12(14), 1–23. https://doi.org/10.3390/rs12142338 Haase, S. M. (1986). Effect of prescribed burning on soil moisture and germination of southwestern ponderosa pine seed on basaltic soils. US Department of Agriculture, Forest Service, Research Note -462, 1–6. Hagmann, R. K., Franklin, J. F., & Johnson, K. N. (2013). Historical structure and composition of ponderosa pine and mixed-conifer forests in south-central oregon. Forest Ecology and Management, 304, 492–504. https://doi.org/10.1016/j.foreco.2013.04.005 Hall, S. A., Burke, I. C., Box, D. O., Kaufmann, M. R., & Stoker, J. M. (2005). Estimating stand structure using discrete-return lidar: An example from low density, fire prone ponderosa pine forests. Forest Ecology and Management, 208(1–3), 189–209. 127 https://doi.org/10.1016/j.foreco.2004.12.001 Hallin, W. E. (1959). The application of unit area control in the management of ponderosa-Jeffrey pine at Black Mountains Experimental Forest. US Department of Agriculture, Forest Service, Techinical Bulletin No. 1191, 96p. https://doi.org/10.1017/CBO9781107415324.004 Hancock, J.M., & Zvelebil, M.J. (2004). Principal components analysis (PCA), in: Dictionary of bioinformatics and computational biology. John Wiely & Sons Inc, New York. https://doi.org/10.1002/9780471650126.dob0562.pub2 Harley, G., Baisan, C., Brown, P., Falk, D., Flatley, W., Grissino-Mayer, H., Hessl, A., Heyerdahl, E., Kaye, M., Lafon, C., Margolis, E., Maxwell, R., Naito, A., Platt, W., Rother, M., Saladyga, T., Sherriff, R., Stachowiak, L., Stambaugh, M., & Taylor, A. H. (2018). Advancing dendrochronological studies of fire in the United States. Fire, 1(11), 1–6. https://doi.org/10.3390/fire1010011 Harr, R. D. (1977). Water flux in soil and subsoil on a steep forested slope. Journal of Hydrology, 33(1–2), 37–58. https://doi.org/10.1016/0022-1694(77)90097-X Harrell, F. E. J., Lee, K. L., Califf, R. M., Pryor, D. B., & Rosati, R. A. (1984). Regression modelling strategies for improved prognostic prediction. Statistics in Medicine, 3(2), 143–152. https://doi.org/10.1002/sim.4780030207 Harrod, R. J., Gaines, W. L., Hartl, W. E., & Camp, A. (1998). Estimating historical snag density in dry forests east of the Cascade Range. USDA Forest Service - General Technical Report PNW, PNW-GTR-428, 116. Hasel, A. A. (1934). Instruction for type map and inventory of experimental forest. PSW Reserach Station, 9p. 128 Hasel, A. A. (1938). Sampling errors in timber survey. Journal of Agricultural Research, 57(10), 713–736. Hasel, A. A. (1942). Estimation of volume timber stands by strip sampling. Annal of Mathematics and Statistics, 13(2), 179–206. Hatfield, J. L., Sauer, T. J., & Cruse, R. M. (2017). Soil: The forgotten piece of the water, food, energy nexus. Advances in Agronomy, 143, 1–46. https://doi.org/10.1016/bs.agron.2017.02.001 Hessburg, P. F., Miller, C. L., Parks, S. A., Povak, N. A., Taylor, A. H., Higuera, P. E., Prichard, S. J., North, M. P., Collins, B. M., Hurteau, M. D., Larson, A. J., Allen, C. D., & Stephens, S. L. (2019). Climate, environment, and disturbance history govern resilience of western north american forests. Forinteirs in Ecology and Evolution 7, 1–27. https://doi.org/10.3389/fevo.2019.00239 Hessburg, P. F., MItchell, R. G., & Filip, G. M. (1994). Historical and current roles of insect and pathogens in eastern Oregon and Washington forestsed landscapes. U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station, 72p. Hessburg, P. F., Smith, B. G., & Salter, R. B. (1999). Detecting Change in Forest Spatial Patterns from refrence conditions. Ecological Applications, 9(4), 1232–1252. https://doi.org/10.1890/1051-0761(1999)009[1232:DCIFSP]2.0.CO;2 Holling, C. S. (1992). Cross-Scale Morphology , Geometry , and Dynamics of Ecosystems. Ecological Monographs, 62(4), 447–502. Hood, S. M., Smith, S. L., & Cluck, D. R. (2010). Predicting mortality for five California conifers following wildfire. Forest Ecology and Management, 260 (2010), 750-762. 129 Hosmer, D. W., & Lemeshow, S. (1989). Applied logistic regression. 2nd ed. John Wiley & Sons, New York. Hotelling, H. (1933). Analysis of a complex of statistical variables into Principal Components. Jour. Educ. Psych., 24, 417-441, 498-520. The Journal of Educational Psychology, 24, 417–441. https://doi.org/10.1037/h0071325. Hurteau, M. D., Liang, S., Martin, K. L., North, M. P., Koch, G. W., & Hungate, B. A. (2016). Restoring forest structure and process stabilizes forest carbon in wildfire-prone southwestern ponderosa pine forests. Ecological Applications, 26(2), 382–391. https://doi.org/10.1890/15-03371.1/ Isaaks, E.H., & Srivastava, R.M. (1989). An Introduction to Applied Geostatistics. Oxford University Press, New York, 413. Jackson, D. A. (1993). Stopping rules in principal components analysis : a comparison of heuristical and statistical approaches. Journal of Ecology, 74(8), 2204–2214. https://doi.org/10.2307/1939574. Jain, A.K., & Dubes, R.C., 1988. Algorithms for clustering data, in: Jain, A.K., Dubes, R. C. (Eds.), Clustering methods and algorithms. Prentice Hall, New Jersey, pp. 55-142. Jaquette, M., Sánchez Meador, A. J., Huffman, D. W., & Bowker, M. A. (2021). Mid-scale drivers of variability in dry mixed-conifer forests of the mogollon rim, arizona. Forests, 12(5). https://doi.org/10.3390/f12050622 Jenkins, J. C., Chojnacky, D. C., Heath, L. S., & Birdsey, R. A. (2003). Comprehensive database of diameter-based biomass resgressions for north american tree species. Department of Agriculture, Forest Service, Rocky Mountain Research Station. Jia, G. J., Burke, I. C., Goetz, A. F. H., Kaufmann, M. R., & Kindel, B. C. (2006). Assessing 130 spatial patterns of forest fuel using AVIRIS data. Remote sensing Environement,102(3-4), 318–327. https://doi.org/10.1016/j.rse.2006.02.025 Kaiser, H. F. (1960). The application of electronic computers to factor analysis. Educational and Psychological Measurement, 20(1), 141–151. https://doi.org/10.1177/001316446002000116 Kashian, D. M., Romme, W. H., Tinker, D. B., Turner, M. G., & Ryan, M. G. (2006). Carbon storage on landscapes with stand-replacing fires. BioScience 56:598-606. https://doi.org/10.1641/0006-3568(2006)56[598:CSOLWS]2.0.CO;2 Keane, R. E., (2001). Mapping wildland fuels for fire management across multiple scales: Integrating remote sensing, GIS, and biophysical modeling. International Journal of Wildland Fire, 10(5), 301-319. https://doi.org/10.1071/WF01028 Keane, R. E., Gray, K., Bacciu, V., & Leirfallom, S. (2012). Spatial scaling of wildland fuels for six forest and rangeland ecosystems of the northern Rocky Mountains, USA. Landscape Ecology, 27(8), 1213–1234. https://doi.org/10.1007/s10980-012-9773-9 Keane, R. E., Hessburg, P. F., Landres, P. B., & Swanson, F. J. (2009). The use of historical range and variability (HRV) in landscape management. Forest Ecology and Management, 258(7), 1025–1037. https://doi.org/10.1016/j.foreco.2009.05.035 Keen, F. P. (1936). Relative Suceptibility of Ponderosa Pines to Bark-Beetle Attack. Journal of Fore, 34, 919–927. https://doi.org/10.1093/jof/34.10.919. Keeton, W. S., & Franklin, J. F. (2005). Do remnant old-growth trees accelerate rates of succession in mature Douglas-fir forests? Ecological Monographs, 75(1), 103–118. https://doi.org/10.1890/03-0626 Kerns, B. K., Thies, W. G., & Niwa, C. G. (2006). Season and severity of prescribed burn in ponderosa pine forests: Implications for understory native and exotic plants. Ecoscience, 131 13(1), 44–55. https://doi.org/10.2980/1195-6860(2006)13[44:SASOPB]2.0.CO;2 Kerry, R., & Oliver, M. A. (2008). Determining nugget:sill ratios of standardized variograms from aerial photographs to krige sparse soil data. Precision Agriculture, 9(1–2), 33–56. https://doi.org/10.1007/s11119-008-9058-0 Keyes, C. R., & Manso González, R. (2015). Climate-influenced ponderosa pine (Pinus ponderosa) seed masting trends in western Montana, USA. Forest Systems, 24(1), 1–7. https://doi.org/10.5424/fs/2015241-05606 Keyser, T. L., Lentile, L. B., Smith, F. W., & Shepperd, W. D. (2008). Changes in forest structure after a large, mixed-severity wildfire in ponderosa pine forests of the black hills, South Dakota, USA. Forest Science, 54(3), 328–338. https://doi.org/10.1093/forestscience/54.3.328 Kilgore, B. M. (1973). The ecological role of fire in Sierran conifer forests. Its application to National Park management. Quaternary Research, 3(3), 496–513. https://doi.org/10.1016/0033-5894(73)90010-0 Knapp, E. E., Lydersen, J. M., North, M. P., & Collins, B. M. (2017). Efficacy of variable density thinning and prescribed fire for restoring forest heterogeneity to mixed-conifer forest in the central Sierra Nevada, CA. Forest Ecology and Management, 406(August), 228–241. https://doi.org/10.1016/j.foreco.2017.08.028 Knapp, E. E., Skinner, C. N., North, M. P., & Estes, B. L. (2013). Long-term overstory and understory change following logging and fire exclusion in a Sierra Nevada mixed-conifer forest. Forest Ecology and Management, 310, 903–919. https://doi.org/10.1016/j.foreco.2013.09.041 Knight, C. A., Cogbill, C. V., Potts, M. D., Wanket, J. A., & Battles, J. J. (2020). Settlement-era 132 Forest structure and composition in the Klamath Mountains: reconstructing a historical baseline. Ecosphere 11(9):e03250. 10.1002/ecs2.3250 Kocher, D. S., & Harris, R. (2007). Tree Growth and competition. In ANR Publication 8235 (Issue 1). https://doi.org/10.1093/oxfordjournals.oep.a040940 Kolb, T. E., Agee, J. K., Fulé, P. Z., McDowell, N. G., Pearson, K., Sala, A., & Waring, R. H. (2007). Perpetuating old ponderosa pine. Forest Ecology and Management, 249(3), 141–157. https://doi.org/10.1016/j.foreco.2007.06.002 Kuehne, C., Weiskittel, A. R., Wagner, R. G., & Roth, B. E. (2016). Development and evaluation of individual tree-and stand-level approaches for predicting spruce-fir response to commercial thinning in Maine, USA. Forest Ecology and Management. 376, 84–95. https://doi.org/10.1016/j.foreco.2016.06.013 Kuhn, M. (2008). Building predictive models in R using the caret package. Journal of Statistical Software, 28(5), 1–26. https://doi.org/10.18637/jss.v028.i05. Kruskal, W. H., & Wallis, W. A. (1952). Use of ranks in one-criterion variance analysis. Journal of the American Statistical Association, 47(260), 583–621. Larson, A. J., & Churchill, D. (2012). Tree spatial patterns in fire-frequent forests of western North America, including mechanisms of pattern formation and implications for designing fuel reduction and restoration treatments. Forest Ecology and Management, 267, 74–92. https://doi.org/10.1016/j.foreco.2011.11.038 Larson, A. J., Stover, K. C., & Keyes, C. R. (2012). Effects of restoration thinning on spatial heterogeneity in mixed-conifer forest. Canadian Journal of Forest Research, 42, 1505–1517. https://doi.org/10.1139/X2012-100 Larsson, S., Oren, R., Waring, R. H., & Barrett, J. W. (1983). Attacks of mountain pine beetle as 133 related to tree vigor of ponderosa pine. Forest Science, 29(2), 395–402. https://doi.org/10.1093/forestscience/29.2.395 Lauvaux, A., Skinner, C. N., & Taylor, A. H. (2016). High severity fire and mixed conifer forest-chaparral dynamics in the southern Cascade Range, USA. Forest Ecology and Management, 363, 74–85. https://doi.org/10.1016/j.foreco.2015.12.016 League, K., & Veblen, T. (2006). Climatic variability and episodic Pinus ponderosa establishment along the forest-grassland ecotones of Colorado. Forest Ecology and Management, 228(1–3), 98–107. https://doi.org/10.1016/j.foreco.2006.02.030 Lefsky, M. A., Cohen, W. B., Parker, G. G., & Harding, D. J. (1999). Lidar Remote Sensing for Ecosystem Studies: Lidar, an emerging remote sensing technology that directly measures the three-dimensional distribution of plant canopies, can accurately estimate vegetation structural attributes and should be of particular interest to forest, landscape, and global ecologists, BioScience, 52(1), 9–30, https://doi.org/10.1641/0006 3568(2002)052[0019:LRSFES]2.0.CO;2 Legendre, P., & Fortin, M. (1989). Spatial pattern and ecological analysis. Vegetatio, 80, 107–138. Lentz, G. H. (1931). The forest survey in the bottomland hardwoods of the Mississippi delta. Journal of Forestry, 29(7), 1046–1059. https://doi.org/10.1093/jof/29.7.1046 Dolph, L., K., Mori, S. R., & Oliver, W. W. (1995). Height-diameter relationship for conifer species on the Blacks Mountain Experimental Forest. Research Paper - US Department of Agriculture, Forest Service, 5. Levin, S. A. (1992). The problem of pattern and scale in ecology. Ecology, 73(6), 1943–1967. https://doi.org/10.2307/1941447 134 Long, J. N. (1985). A practical approach to density management. The Forestry Chronicle, 61(1), 23–27. https://doi.org/10.5558/tfc61023-1 Long, J. N., & Daniel, T. W. (1990). Assessment of growing stock in uneven-aged stands, Western Journal of Applied Forestry, 5(3), 93–96. https://doi.org/10.1093/wjaf/5.3.93. Long, J. N., & Shaw, J. D. (2005). A density management diagram for even-aged sierra nevada mixed-conifer stands. Western Journal of Applied Forestry, 27(4), 187–195. https://doi.org/10.5849/wjaf.11-036 Long, J. N., Windmuller-Campione, M., & De Rose, R. J. (2018). Building resistance and resilience: Regeneration should not be left to chance. Forests, 9(5), 270. https://doi.org/10.3390/f9050270 Lowdermilk, W. C. (1927). A method for rapid surveys of vegetation. Journal of Forestry, 25(2), 181–185. https://doi.org/10.1093/jof/25.2.181 Lydersen, J. M., North, M. P., Knapp, E. E., & Collins, B. M. (2013). Quantifying spatial patterns of tree groups and gaps in mixed-conifer forests: Reference conditions and long-term changes following fire suppression and logging. Forest Ecology and Management, 304, 370–382. https://doi.org/10.1016/j.foreco.2013.05.023 Lydersen, J., & North, M. (2012). Topographic variation in structure of mixed-conifer forests under an active-fire regime. Ecosystems, 15, 1134–1146. https://doi.org/10.1007/s10021-012-9573-8 Lynch, D. W., & Schumacher, F. X. (1941). Concerning the dispersion of natural regeneration. Journal of Forestry, 39(1), 49–51. https://doi.org/10.1093/jof/39.1.49 MacArthur, H. (1957). On the relative abundance of bird species. PNAS, 293–295. 135 https://doi.org/10.1073 Mauro, F., Ritchie, M., Wing, B., Frank, B., Monleon, V., Temesgen, H., & Hudak, A. (2019). Estimation of changes of forest structural attributes at three different spatial aggregation levels in northern california using multitemporal LiDAR. Remote Sensing, 11, 1–31. https://doi.org/10.3390/rs11080923 McCullough, G. D., & Wagner, R. (1987). Evaluation of four techniques to assess vvigor of water-stressed ponderosa pine. Canadian Journal of Forest Research, 17, 138–145. https://doi.org/10.1139/x87-025 McDonald, P. M. (1983). Clearcutting and natural regeneration: management implications for the northern Sierra Nevada. United States Department of Agriculture, Forest Service PSW-GTR-70, 11p. McKelvey, K.S., Skinner, C. N., Chang, C., Erman, D. C., Husari, S. J., Parsons, D. J., Van Wagtendonk, J. W., & Weatherspoon, C. P. (1996). An overview of fire in the Sierra Nevada. Sierra Nevada Ecosystem Project Final Report to Congress, Water Resources Center Report No. 37. Davis, CA: Centers for Water and Wildland Resources, University of California, 8p. https://doi.org/10.1016/S0273-1223(96)00755-X McKelvey, Kevin S, & Johnston, J. (1992). Historical perspectives on forests of the Sierra Nevada and the Transverse Ranges of southern California: Forest conditions at the turn of the century. US Department of Agriculture, Forest Service, PSW-GTR-133, January, 225–246. McRoberts, R. E., Tomppo, E. O., & Næsset, E. (2010). Advances and emerging issues in national forest inventories. Scandinavian Journal of Forest Research, 25(4), 368–381. https://doi.org/10.1080/02827581.2010.496739 136 Merschel, A. G., Beedlow, P. A., Shaw, D. C., Woodruff, D. R., Lee, E. H., Cline, S. P., Comeleo, R. L., Hagmann, R. K., & Reilly, M. J. (2021). An ecological perspective on living with fire in ponderosa pine forests of Oregon and Washington: Resistance, gone but not forgotten. Trees, Forests and People, 4(18), 100074. https://doi.org/10.1016/j.tfp.2021.100074 Millar, C. I., Stephenson, N. L., & Stephens, S. L. (2007). Climate change and forests of the future: Managing in the face of uncertainity. Ecological Applications, 17(8), 2145–2151. https://doi.org/10.1890/06-1715.1 Moore, D. A., & Carpenter, T. E. (1999). Spatial analytical methods and geographic information systems: Use in health research and epidemiology. Epidemiologic Reviews, 21(2), 143–161. https://doi.org/10.1093/oxfordjournals.epirev.a017993 Moral, R. A., Hinde, J., & Demétrio, C. G. B. (2017). Half-normal plots and overdispersed models in R: The hnp package. Journal of Statistical Software, 81(10), 1–23. https://doi.org/10.18637/jss.v081.i10 Moran, P. (1950). A test for the serial independence of residuals. Biometrika, 37, 178–181. http://dx.doi.org/10.1093/biomet/37.1-2.178 Norman, S. P., & Taylor, A. H. (2003). Tropical and north Pacific teleconnections influence fire regimes in pine-dominated forests of north-eastern California, USA. Journal of Biogeography, 30: 1081-1092. https://doi.org/10.1046/j.1365-2699.2003.00889.x Norman, S. P., & Taylor, A. H. (2005). Pine forest expansion along a forest-meadow ecotone in northeastern California, USA. Forest Ecology and Management, 215(1–3), 51–68. https://doi.org/10.1016/j.foreco.2005.05.003. North, M., Chen, J., Oakley, B., Song, B., Rudnicki, M., Gray, A., & Innes, J. (2004). Forest 137 stand structure and pattern of old-growth western hemlock/Douglas-fir and mixed-conifer forests. Forest Science, 50(3), 299–311. North, M., Oakley, B., Chen, J., Erickson, H., Gray, A., Izzo, A., Johnson, D., Ma, S., Marra, J., Meyer, M., Purcell, K., Rambo, T., Rizzo, D., Roath, B., & Schowalter, T. (2002). Vegetation and ecological characteristics of mixed conifer and red fir forests at the Teakettle Experimental Forest. General Technical Report PSW-GTR-186, 52. http://ir.library.oregonstate.edu/xmlui/handle/1957/17581. North, M. P., Tompkins, R. E., Bernal, A. A., Collins, B. M., Stephens, S. L., & York, R. A. (2022). Operational resilience in western US frequent-fire forests. Forest Ecology and Management, 507(December 2021). https://doi.org/https://doi.org/10.1016/j.foreco.2021.120004. North, M., Stine, P., O’Hara, K., Zielinski, W., & Stephens, S. (2009). An Ecosystem Management Strategy for Sierran Mixed- Conifer Forests. Forest Service U.S. Department of Agriculture, 49 p. Oester, P. T., Shaw, D. C., & Filip, G. M. (2018). Managing insects and disease of Oregon conifers. Oregon State University Press, Corvallis Oregon State University Extension Publication EM8980, 134p. Oksanen, A. J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., Hara, R. B. O., Simpson, G. L., Solymos, P., Stevens, M. H. H., & Szoecs, E. (2022).Vegan: Community ecology package. R Package Version 2.6-2. http://cran.r-project.org/web/packages/vegan/index.html. Oliver, C. D., & Larson, B. C. (1996). Forest stand dynamics: updated edition. In Forest stand dynamics. Formerly published by John Wiley & Sons; Copyright now held by C.D.Oliver & 138 B.C. Larson. Oliver, W. (1995). Is self-thinning in ponderosa pine ruled by Dendroctonus bark beetles? General Technical Report - US Department of Agriculture Forest Service, RM-GTR-267, 213–218. Oliver, W. (2000). Ecological research at the Blacks Mountain Experimental Forest in northeastern California. United Stated Forest Service, General Technical Report , PSW-GTR-179, 66p. Oliver, W. W., & Dolph, L. K. (1992). Mixed-conifer seedling growth varies in response to overstory release. Forest Ecology and Management, 48(1–2), 179–183. https://doi.org/10.1016/0378-1127(92)90130-2 Oliver, W. W., & Power, R. F. (1971). Early height growth of ponderosa pine: Forecast dominance in plantations. USDA Forest Service Research Note PSW-250 ., 0–3. Oliver, W. (1967). Ponderosa pine can stagnate on a good site. Journal of Forestry, 65(11), 814–816. https://doi.org/10.1093/jof/65.11.814 Olsen, W. K., Schmid, J. M., & Mata, S. A. (1995). Stand characteristics associated with mountain pine beetle infestations in ponderosa pine. Forest Science, 42(3), 310–327. https://doi.org/10.1093/forestscience/42.3.310 Oswalt, S. N., Smith, W. B., Miles, P. D., & Pugh, S. A. (2014). Forest resources of the united states , 2012 : A Technical Document Supporting the Forest Service Update of the 2010 RPA Assessment. Forest Service U.S. Department of Agriculture, October, 218 p. Ouzts, J., Kolb, T., Huffman, D., & Sánchez Meador, A. (2015). Post-fire ponderosa pine regeneration with and without planting in Arizona and New Mexico. Forest Ecology and Management, 354(June), 281–290. https://doi.org/10.1016/j.foreco.2015.06.001 139 Patton, S. C., Grulke, N., & Bienz, C. (2021). Assessment of ponderosa pine vigor using four-band aerial imagery in south central oregon: Crown objects to landscapes. Forests, 12(5). https://doi.org/10.3390/f12050612 Pawlikowski, N. C., Coppoletta, M., Knapp, E., & Taylor, A. H. (2019). Spatial dynamics of tree group and gap structure in an old-growth ponderosa pine-California black oak forest burned by repeated wildfires. Forest Ecology and Management, 434(October 2018), 289–302. https://doi.org/10.1016/j.foreco.2018.12.016 Pearson, G. A. (1950). Management of ponderosa pine in the Southwest: As developed by research and experimental practice. US Department of Agriculture, Forest Service, Agriculture Monograph No.6, 218p. https://doi.org/10.4324/9781351250726-14 Pebesma, E. J. (2004). Multivariable geostatistics in S: the gstat package. Computers & Geosciences, 30, 683–691. https://doi.org/10.1016/j.cageo.2004.03.012 Perrakis, D. D. B., Agee, J. K., & Eglitis, A. (2011). Effects of prescribed burning on mortality and resin defenses in old growth ponderosa pine (Crater Lake, Oregon): Four years of post-fire monitoring. Natural Areas Journal, 31(1), 14–25. https://doi.org/10.3375/043.031.0103 Perry, G. L. W., Enright, N. J., Miller, B. P., & Lamont, B. B. (2009). Nearest-neighbour interactions in species-rich shrublands : the roles of abundance , spatial patterns and resources. Oikos, 118(2), 161-174. https://doi.org/10.1111/j.1600-0706.2008.16947.x Pielou, E. C. (1962). The use of plant-to-neighbour distances for the detection of competition. Journal of Ecology, 50(2), 357–367. https://www.jstor.org/stable/2257448 Plant, R. E. (2012). Spatial data analysis in ecology and agriculture using R. Taylor & Francis Group. Boca Raton, New York. PSW Research Station. (1934). Instruction for type map and inventory of experimental forest; 140 Technical Report; United States Dep. Agric. For. Serv.: Redding, CA, USA, p. 9. Puhlick, J. J., Laughlin, D. C., & Moore, M. M. (2012). Factors influencing ponderosa pine regeneration in the southwestern USA. Forest Ecology and Management, 264, 10–19. https://doi.org/10.1016/j.foreco.2011.10.002 Radeloff, V. C., Miller, T. F., He, H. S., & Mladenoff, D. J. (2000). Periodicity in spatial data and geostatistical models: Autocorrelation between patches. Ecography, 23(1), 81–91. https://doi.org/10.1111/j.1600-0587.2000.tb00263.x Reed, T. W., Gulland, E., West, G., Mcmeekin, D. A., & Moncrieff, S. (2016). Geographic metadata searching with semantic and spatial filtering methods. GEOProcessing 2016: The Eighth International Conference on Advanced Geographic Information Systems, Applications, and Services. Reynolds, R. T., Sánchez Meador, A. J., Youtz, J. A., Nicolet, T., Matonis, M. S., Jackson, P. L., Delorenzo, D. G., Graves, A. D., Richard, T., Meador, S., Andrew, J., James, A., Megan, S., Patrick, L., Donald, G., & Andrew, D. (2013). Restoring composition and structure in Southwestern frequent-fire forests: A science-based framework for improving ecosystem resiliency. General Technical Report RMRS-GTR-310, 76. Riegel, G. M., Miller, R. F., & Krueger, W. (1992). Competition for resources between understory vegetation and overstory pinus ponderosa in northeastern Oregon. Ecological Applications, 2 (1), 71–85. https://doi.org/10.2307/1941890 Ritchie, M. W. (2016). Multi-scale reference conditions in an interior pine-dominated landscape in northeastern California. Forest Ecology and Management, 378, 233–243. https://doi.org/10.1016/j.foreco.2016.07.017 Ritchie, M. W. (2020). Relating stocking and density for natural regeneration of conifers in 141 northern California. Forest Science, 66(5), 521–525. https://doi.org/10.1093/forsci/fxaa005 Ritchie, M. W., Zhang, J., & Hamilton, T. A. (2013). Aboveground tree biomass for Pinus ponderosa in northeastern California. Forests, 4(1), 179–196. https://doi.org/10.3390/f4010179 Roel, A., & Plant, R. E. (2004). Site specific analysis. Factors underlying yield variability in two California rice fields. Agronomy Journal, 96(5), 1481–1494. https://dopi.org/10.2134/agronj2004.1481. Rollins, M. G., Keane, R. E., & Parson, R. A. (2004). Mapping fuels and fire regimes using remote sensing, ecosystem simulation, and gradient modeling. Ecological Applications, 14(1), 75–95. https://www.jstor.org/stable/4493523 Rossi, R. E., Mulla, D. J., Journel, A. G., & Franz, E. H. (1992). Geostatistical tools for modeling and interpreting ecological spatial dependence. Ecological Monographs, 62(2), 277–314. https://doi.org/10.2307/2937096 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. (2020). lidR: An R package for analysis of Airborne Laser Scanning (ALS) data. Remote Sensing of Environment, 251(February), 112061. https://doi.org/10.1016/j.rse.2020.112061 Sackett, S. S., Haase, S. M., & Harrington, M. G. (1996). Prescribed burning in Southwestern ponderosa pine. Effects Of Fire On Madrean Province Ecosystems - A Symposium Proceedings, 289, 178–186. Salman, K. A., & Bongberg, J. W. (1942). Logging high risk trees to control insects in the pine stands of northern California. Journal of Forestry, 40, 533–549. Sanchez Meador, A. J., Parysow, P. F., & Moore, M. M. (2010). Historical stem-mapped 142 permanent plots increase precision of reconstructed reference data in ponderosa pine Forests of Northern Arizona. Restoration Ecology, 18(2), 224–234. https://doi.org/10.1111/j.1526-100X.2008.00442.x Savage, M., Brown, P. M., & Feddema, J. (1996). The role of climate in a pine forest regeneration pulse in the southwestern United States. Ecoscience, 3(3), 310–318. https://doi.org/10.1080/11956860.1996.11682348 Savage, M., & Mast, J. N. (2005). How resilient are southwestern ponderosa pine forests after crown fires? Canadian Journal of Forest Research, 35(4), 967–977. https://doi.org/10.1139/x05-028 Scholl, A. E., & Taylor, A. H. (2010). Fire regimes, forest change, and self-organization in an old-growth mixed-conifer forest, Yosemite National Park, USA. Ecological Applications, 20(2), 362–380. https://doi.org/10.1890/08-2324.1 Schubert, G. H. (1974). Silviculture of southwestern ponderosa pine : the status of our knowledge. United Stated Forest Service, Research,Paper RM-123, pp.71. https://doi.org/10.5962/bhl.title.98816 Seidl, R., Spies, T. A., Peterson, D. L., Stephens, S. L., & Hickle, J. A. (2016). Searching for resilience: addressing the impacts of changing disturbance regimes on forest ecosystem services. Journal of Applied Ecology, 53(1), 120–29. https://doi.org/10.1111/1365-2664.12511 Shan, Y., Bettinger, P., Cieszewski, C. J., & Li, R. T. (2009). Trends in spatial forest planning. International Journal of Mathematical and Computational Forestry & Natural Rsources, 1(2), 86-112. Shepperd, W. D., Edminster, C. B., & Mata, S. A. (2006). Long-term seedfall, establishment, 143 survival, and growth of natural and planted ponderosa pine in the Colorado Front Range. Western Journal of Applied Forestry, 21(1), 19–26. https://doi.org/10.1093/wjaf/21.1.19 Sheridan, R. D., Popescu, S. C., Gatziolis, D., Morgan, C. L. S., & Ku, N. W. (2015). Modeling forest aboveground biomass and volume using airborne LiDAR metrics and forest inventory and analysis data in the Pacific Northwest. Remote Sensing, 7, 229–255. https://doi.org/10.3390/rs70100229 Sheskin, D. J. (2000). Handbook of parametric and nonparametric statistical procedures. In Chapman & Hall/CRC. Shive, K. L., Preisler, H. K., Welch, K. R., Safford, H. D., Butz, R. J., O’Hara, K. L., & Stephens, S. L. (2018). From the stand scale to the landscape scale: predicting the spatial patterns of forest regeneration after disturbance. Ecological Applications, 28(6), 1626–1639. https://doi.org/10.1002/eap.1756 Show, S. B., & Kotok, E. I. (1924). The role of fire in the California pine forests. United States Department of Agriculture, Technical BUlletin 1294, 80p. https://doi.org/10.5962/bhl.title.107935 Shoot, C., Andersen, H. E., Monika Moskal, L., Babcock, C., Cook, B. D., & Morton, D. C. (2021). Classifying forest type in the national forest inventory context with airborne hyperspectral and lidar data. Remote Sensing, 13(10), 1–18. https://doi.org/10.3390/rs13101863 Skinner, C. N., & Taylor, A. H. (2006). Southern cascades bioregion. In Sugihara, N. G.; van Wagtendonk, J. W.; Fites-Kaufmann, J.; Shaffer, K. E.; Thode, A. E. (Ed.), Fire in California’s ecosystems. (p. pp.195-224). University of California Press, Berkeley. Skinner, C.N., & Taylor, A.H. (2018). Chapter twelve. Southern Cascades Bioregion. In Fire in 144 California’s ecosystems, second edition, ed. Jan van Wagtendonk, 195–218. Berkeley: University of California Press. https://doi.org/10.1525/ 9780520961913-015. Skowronski, N. S., & Lister, A. J. (2012). The utility of lidar for large area forest inventory applications. Forest Inventory and Analysis Symposium GTR-NRS-P-105, April. Smith, T. F., Rizzo, D. M., & North, M. P. (2005). Patterns of mortality in an old-growth mixed-conifer forest of the Southern Sierra Nevada, California. Forest Science, 51(3), 266–275. https://doi.org/10.1093/forestscience/51.3.266 Soil Survey Staff. (1951). Soil survey manual. USDA Handbook No. 18, 503 p. Stage, A. (1976). An expression for the effect of aspect, slope, and habitat type on tree growth. Forest Science, 22(4), 457–460. https://doi.org/10.1093/forestscience/22.4.457 Stein, S. J., & Kimberling, D. N. (2003). Germination,establishment, and mortality of naturally seeded southwestern ponderosa pine. Western Journal of Applied Forestry, 18(2), 109–114. https://search.proquest.com/docview/220297257?accountid=12834 Stein, W. I. (1992). Regeneration surveys and evaluation. In S.D. Hobbs (Ed.), Reforestation practices in southwestern Oregon and northern California (pp. 346–382). Corvallis, Or. : Forest Research Laboratory, Oregon State University. Stephens, S. L. (2000). Mixed conifer and red fir forest structure and uses in 1899 from the central and northern Sierra Nevada, California. Madroño, 47(1), 43–52. https://doi.org/10.1890/ES14-00046.1 Stephens, S. L., Fry, D. L., & Franco-Vizcaíno, E. (2008). Wildfire and spatial patterns in forests in Northwestern Mexico: United States wishes it had similar fire problems. Ecology and Society, 13(2). https://doi.org/10.5751/ES-02380-130210 Stephens, S. L., Lydersen, J. M., Collins, B. M., Fry, D. L., & Meyer, M. D. (2015). Historical 145 and current landscape-scale ponderosa pine and mixed conifer forest structure in the southern Sierra Nevada. Ecosphere, 6(5), 1–63. https://doi.org/10.1890/ES14-00379.1 Stephens, S. L., Skinner, C. N., & Gill, S. J. (2003). Dendrochronology-based fire history of Jeffrey pine - Mixed conifer forests in the Sierra San Pedro Martir, Mexico. Canadian Journal of Forest Research, 33(6), 1090–1101. https://doi.org/10.1139/x03-031 Stephenson, N. L. (1998). Actual evapotranspiration and deficit: Biologically meaningful correlates of vegetation distribution across spatial scales. Journal of Biogeography, 25(5), 855–870. https://doi.org/10.1046/j.1365-2699.1998.00233.x Stoddard, J. L., Larsen, D. P., Hawkins, C. P., Johnson, R. K., & Norris, R. H. (2006). Setting expectations for the ecological condition of streams: The concept of reference condition. Ecological Applications, 16(4), 1267–1276. https://doi.org/10.1890/1051-0761(2006)016[1267:SEFTEC]2.0.CO;2 Swetnam, T. W. (1993). Fire history and climate chage in Giant Sequoia groves. Science, 262, https://885–889. doi/10.1126/science.262.5135.885 Taylor, A. H. (2000). Fire regimes and forest changes in mid and upper montane forests of the southern Cascades , Lassen. Journal of Biogeography, 27, 87–104. http://dx.doi.org/10.1046/j.1365-2699.2000.00353.x Taylor, A. H. (2010). Fire disturbance and forest structure in an old-growth Pinus ponderosa forest, southern Cascades, USA. Journal of Vegetation Science, 21(3), 561–572. https://doi.org/10.1111/j.1654-1103.2009.01164.x Taylor, A. H., & Skinner, C. N. (2003). Spatial patterns and controls on historical fire regimes and forest structure in the Klamath Mountains. Ecological Applications, 13(3), 704–719. https://doi.org/10.1890/1051-0761(2003)013[0704:SPACOH]2.0.CO;2 146 Trotsiuk, V., Pederson, N., Druckenbrod, D. L., Orwig, D. A., Bishop, D. A., Barker-Plotkin, A., Fraver, S., & Martin-Benito, D. (2018). Testing the efficacy of tree-ring methods for detecting past disturbances. Forest Ecology and Management, 425(March), 59–67. https://doi.org/10.1016/j.foreco.2018.05.045 van Mantgem, P. J., & Stephenson, N. L. (2007). Apparent climatically induced increase of tree mortality rates in a temperate forest. Ecology Letters, 10, 909–916. https://doi.org/10.1111/j.1461-0248.2007.01080.x Van Pelt, R., & Franklin, J. F. (1999). Response of understory trees to experimental gaps in old-growth Douglas-fir forests. Ecological Applications, 9(2), 504–512. https://doi.org/10.1890/1051-0761(1999)009[0504:ROUTTE]2.0.CO;2 Van Wagtendonk, J. W., Sugihara, N. G., Stephens, S. L., Thode, A. E., & Fites-Kaufman, J. (2006). Sierra Nevada bioregion. In N. G. Sugihara, J. W. Van Wagtendonk, J. Fites-Kaufman, K. E. Shaffer, & A. E. Thode (Eds.), Fire in California’s Ecosystems (pp. 264–294). University of California Press. https://doi.org/10.1525/california/9780520246058.003.0012 Vankat, J. L., & Major, J. (1978). Vegetation changes in sequoia national park, California. Journal of Biogeography, 5(4), 377–402. https://www.jstor.org/stable/3038030. Vorster, A. G., Evangelista, P. H., Stovall, A. E. L., & Ex, S. (2020). Variability and uncertainty in forest biomass estimates from the tree to landscape scale : the role of allometric equations. Carbon Balance and Management, 1–20. https://doi.org/10.1186/s13021-020-00143-6 Weaver, H. (1943). Fire as an ecological and silvicultural factor in the ponderosa pine region of the pacific slope. Fire Ecology, 41(1), 7–15. https://doi.org/10.4996/fireecology.1001001 147 Wilcoxon, F. (1945). Individual comparisons by ranking methods. Biometrics Bulletin, 1(6), 80–83. https://doi.org/10.2307/3001968 Wilson, J. P. (2012). GIScience research at the thirty-second annual Esri International User Conference. Transactions in GIS, 16(3), 267–269. https://doi.org/10.1111/j.1467-9671.2012.01345.x Wing, B. M., Ritchie, M. W., Boston, K., Cohen, W. B., Gitelman, A., & Olsen, M. J. (2012). Prediction of understory vegetation cover with airborne lidar in an interior ponderosa pine forest. Remote Sensing of Environment, 124, 730–741. https://doi.org/10.1016/j.rse.2012.06.024 Wing, M. W., Boston, K., & Ritchie, M. W. (2019). A Technique for Implementing Group Selection Treatments with Multiple Objectives Using an Airborne Lidar-Derived Stem Map in a Heuristic Environment, Forest Science, 65(2), 211- 222, https://doi.org/10.1093/forsci/fxy050 Woodall, C. W., Miles, P. D., & Vissage, J. S. (2005). Determining maximum stand density index in mixed species stands for strategic-scale stocking assessments. Forest Ecology and Management, 216(1–3), 367–377. https://doi.org/10.1016/j.foreco.2005.05.050 Youngblood, A., Max, T., & Coe, K. (2004). Stand structure in eastside old-growth ponderosa pine forests of Oregon and northern California. Forest Ecology and Management, 199(2–3), 191–217. https://doi.org/10.1016/j.foreco.2004.05.056 Zald, H. S. J., Gray, A. N., North, M., & Kern, R. A. (2008). Initial tree regeneration responses to fire and thinning treatments in a Sierra Nevada mixed-conifer forest, USA. Forest Ecology and Management, 256(1–2), 168–179. https://doi.org/10.1016/j.foreco.2008.04.022 Zawadzki, J., Cieszewski, C. J., Zasada, M., & Lowe, R. C. (2005). Applying geostatistics for 148 investigations of forest ecosystems using remote sensing imagery. Silva Fennica Monographs, 39(4), 599–617. https://doi.org/10.14214/sf.369 Zhang, J., Finley, K. A., Johnson, N. G., & Ritchie, M. W. (2019). Lowering stand density enhances resiliency of ponderosa pine forests to disturbances and climate change. Forest Science, 65(4), 496–507. https://doi.org/10.1093/forsci/fxz006 Zhang, J., Oliver, W. W., Ritchie, M. W., & Neal, D. L. (2013). Overstory and understory dynamics in a ponderosa pine plantation vary with stand density in the Sierra Nevada: 40-year results. Forest Science, 59(6), 670–680. https://doi.org/10.5849/forsci.10-033 Zhang, J., Ritchie, M. W., & Oliver, W. W. (2008). Vegetation responses to stand structure and prescribed fire in an interior ponderosa pine ecosystem. Canadian Journal of Forest Research, 38, 909-918. https://doi.org/10.1139/X07-230 Ziegler, J. P., Hoffman, C. M., Fornwalt, P. J., Sieg, C. H., Battaglia, M. A., Chambers, M. E., & Iniguez, J. M. (2017). Tree regeneration spatial patterns in ponderosa pine forests following stand-replacing fire : Influence of topography and neighbors. Forests, 8(391),1–15. https://doi.org/10.3390/f8100391 Ziegler, J. P., Hoffman, C., Battaglia, M., & Mell, W. (2017). Spatially explicit measurements of forest structure and fire behavior following restoration treatments in dry forests. Forest Ecology and Management, 386, 1–12. https://doi.org/10.1016/j.foreco.2016.12.002 Zou, B. C., Breshears, D. D., Newman, D, B., Wilcox, P. B., Gard, O. M., & Rich, M. P. (2008). Soil water dynamics under low-versus high-ponderosa pine tree density: ecohydrological functioning and restoration implications. Ecohydrology, 1, 309-315. https://doi.org/10.1002/eco Zouhar, K. (2001). Abies concolor. In: Fire effects information system, [Online]. U.S. 149 Department of Agriculture, Forest Service, Rocky Mountain Research Station, Fire Sciences Laboratory (Producer). https://www.fs.fed.us/database/feis/plants/tree/abicon/all.html Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A., & Smith, G. M. (2009). Mixed effects models and extensions in ecology with R. Springer; https://doi.org/10.1007/978-0-387-87458-6_9 150 Appendices Appendix A Additional tables for Chapter 2 and 3 Table A.1 Summary of the applied combination of structural diversity treatments and prescribed burn in various units of BMERP. Unit Structural Treatments Treatment Year Prescribed Burn Burn Year # of 1 ha Plots RNA-A No treatment - Unburned - 55 RNA-B No treatment - Burned 1997 63 RNA-C No treatment - Burned 1997 76 RNA-D No treatment - Unburned - 35 UNIT-39 LOD 1996 Burned 1997 89 UNIT-39 LOD 1996 Unburned - 85 UNIT-40 LOD 1998 Burned 2000 87 UNIT40 LOD 1998 Unburned - 79 UNIT-43 LOD 1996 Burned 1997 68 UNIT-43 LOD 1996 Unburned - 90 UNIT-44 LOD 1997 Burned 1999 56 UNIT-44 LOD 1997 Unburned - 49 UNIT-45 LOD 1997 Burned 1999 57 UNIT-45 LOD 1997 Unburned - 75 UNIT-38 HID 1996 Burned 1997 87 UNIT-38 HID 1996 Unburned - 93 UNIT-41 HID 1996 Burned 1997 82 UNIT-41 HID 1996 Unburned - 79 UNIT-42 HID 1997 Burned 1999 83 UNIT-42 HID 1997 Unburned - 101 UNIT-47 HID 1997 Burned 1999 59 UNIT-47 HID 1997 Unburned - 48 UNIT-48 HID 1998 Burned 2000 81 UNIT-48 HID 1998 Unburned - 76 Total 1720 HID: high structural diversity treatment; LOD: low structural diversity treatment; RNA: research natural area 151 Table A.2 Brief description of each Dunning’s vigor classes (Dunning 1928). Class Age Crown position Crown size Vigor Crown height 1 Young or Thrifty (50-150) Isolated or dominant Average or wider good 60 % or more of total height 2 Young or Thrifty (50-150) Codominant Average or narrower Good or moderate Less than 60 % of total height 3 Mature (150-300) Isolated or dominant Average or wider Moderate 65 % or more of total height 4 Mature (150-300) Codominant Average or narrower poor Less than 65 % of total height 5 Over mature (>300) Isolated or dominant Any size Poor Less than 65 % of total height 6 Young or Thrifty (50-150) Intermediate or suppressed Any size Good Less than 65 % of total height 7 Mature or Over mature Intermediate or suppressed Any size Poor Less than 65 % of total height 152 Table A.3 Summary of observed and broken-stick model eigenvalues along with variance explained by each of the principal components. Principal components (PC) Eigenvalue Variance percentage Cumulative variance Observed Broken Stick Model PC1 3.45 3.25 24.66 31.51 PC2 2.25 2.24 16.04 40.70 PC3 1.74 1.75 12.45 53.16 PC4 1.07 1.41 7.71 60.87 PC5 0.96 1.16 8.87 67.74 PC6 0.83 0.96 5.93 73.67 PC7 0.75 0.80 5.38 79.06 PC8 0.58 0.65 4.19 83.26 PC9 0.52 0.53 3.77 87.03 PC10 0.46 0.42 3.71 90.75 PC=principal components; PCA; principal component analysis. 153 Table A.4 Summary of gamma statistics (γ) and Moran’s I for cluster sizes k = 2 to 10 Cluster size k Gamma statistic (γ) Moran’s I 2 0.99 0.59 3 0.59 0.51 4 0.63 0.31 5 1.00 0.10 6 0.54 0.14 7 0.47 0.12 8 0.41 0.25 9 0.51 0.10 10 0.29 0.09 154 Table A.5 Stand density index (SDI) by seven vigor classes for two species group: ponderosa pine (PP) and cedar-fir (CF) with standard deviation in parentheses. Ponderosa pine Clusters SDIPP1 SDIPP2 SDIPP3 SDIPP4 SDIPP5 SDIPP6 SDIPP7 Number (%) Mixed-High 35.53 (31.33) 17.13 (15.91) 32.67 (24.60) 10.70 (11.37) 100.70 (43.29) 8.61 (8.26) 1.42 (2.11) 304 (8) PP-Low 30.94 (16.61) 13.66 (10.17) 71.17 (41.77) 38.02 (20.32) 147.29 (59.58) 3.91 (4.22) 6.28 (4.27) 723 (19) PP-High 62.55 (31.85) 22.53 (15.88) 28.82 (26.06) 8.18 (10.06) 90.21 (50.81) 5.34 (6.12) 0.95 (1.74) 1821 (47) Mixed-Low 29.58 (17.30) 12.04 (8.15) 39.14 (30.18) 11.61 (10.19) 128.87 (52.95) 3.16 (3.01) 1.53 (1.34) 982 (26) Total 158.60 65.36 171.80 68.51 467.06 21.02 10.17 3830 Cedar-fir SDICF1 SDICF2 SDICF3 SDICF4 SDICF5 SDICF6 SDICF7 Number (%) Mixed-High 40.80 (18.82) 28.75 (14.56) 19.13 (13.87) 7.47 (6.8) 37.56 (26.88) 13.54 (11.26) 0.57 (0.37) 304 (8) PP-Low 6.40 (7.19) 3.38 (4.94) 5.86 (9.26) 7.16 (12.73) 20.98 (18.56) 0.73 (0.68) 0.23 (0.45) 723 (19) PP-High 1.36 (4.11) 0.66 (2.17) 0.51 (2.19) 0.20 (1.23) 2.13 (7.14) 0.24 (1.06) 0.00 (0.00) 1821 (47) Mixed-Low 19.46 (11.43) 10.82 (8.49) 11.11 (9.27) 4.15 (5.98) 42.70 (33.03) 2.49 (2.45) 0.18 (0.14) 982 (26) Total 68.02 43.61 36.62 18.98 103.37 17.00 0.98 3830 SDI: stand density index; PP = ponderosa pine; CF = incense-cedar and white fir, number in the subscript to the seven Dunning vigor classes, where 1, 2, 3, and 4 are high and 4, 5, and 7 are low vigor classes 155 A6. Process of Cluster size determination The cluster size k determines how many clusters are identified as part of the cluster analysis. For example, the analysis where cluster size is set to k = 4 resulted in four clusters of plots, where each plot belonged to one of the four clusters. We then evaluated the results for each cluster size k with regards to internal consistency, spatial structure, and hierarchical arrangement of the plots within each cluster following Plant (2012) and Roel and Plant (2004). The results from the cluster analysis for each cluster size k were translated into two-dimensional PCA-biplots of the first two principal components (Hancock & Zvelebil, 2004). Gamma statistic to assess the internal consistency of the clusters—Internal consistency indicates the similarity of plots within a cluster for a given cluster size k in terms of their Euclidean distance and the values of their overstory variables (Plant, 2012; Roel and Plant, 2004). To calculate the gamma statistic, for each cluster size k, we replicated the cluster analysis with ten different initial seeds (Plant, 2012). For a given cluster size k, internal consistency is achieved when the plots are partitioned into the same cluster in each replication (Plant, 2012; Roel and Plant, 2004). Internal consistency was measured using the gamma statistic (γ) for each cluster size k = 2 to 10 (Plant, 2012; Roel and Plant, 2004). Values of γ closer to 1 indicate that a given cluster size k resulted in internally consistent clusters. Values closer to 0 indicate a lack of internal consistency (Plant, 2012; Roel and Plant, 2004). Moran’s I to assess spatial autocorrelation of the clusters—We calculated Moran’s I (Moran, 1950) for each of the cluster sizes k = 2 to 10 to detect spatial autocorrelation as an indicator of spatial aggregation of plots belonging to a cluster. Moran’s I range from -1 to 1. Values between 0 and 1 indicate positive spatial autocorrelation, meaning that plots similar in 156 SDI by species group and vigor classes are aggregated together. A value of 0 indicates perfect spatial randomness of plots belonging to a cluster (Moran, 1950). Visual assessment for indication of hierarchical arrangement—Moran’s I and γ may be high for the cluster size k+1 compared to cluster size k, yet clusters in cluster size k+1 may be hierarchically arranged (Plant, 2012; Roel and Plant, 2004). When the cluster size switched from k to k+1, we looked at whether the additional cluster in cluster size k+1 was a subset of the existing cluster in the cluster size k as an indication of hierarchical arrangement (Roel and Plant, 2004). Therefore, we selected a cluster size for which the clusters were not hierarchically arranged. 157 Appendix B LiDAR point cloud processing and model building process B1. LiDAR processing The LiDAR point cloud was processed using the ‘lidR’ package in R (Roussel et al., 2020). I filtered the duplicates for the overlapping region of the point clouds and then classified the ground returns (Roussel et al., 2020). All the ground returns were used to obtain a digital terrain model of 1 x 1 m grid cell size (Roussel et al., 2020). The point cloud was normalized by subtracting the digital terrain model from all returns to remove the influence of terrain on above-ground returns (Roussel et al., 2020). From the normalized point cloud, I used the area-based approach (Roussel et al., 2020) to calculate the LiDAR metrics (Table C.1)—hereafter referred to as auxiliary variables—using the ‘cloudmetrics’ function related to the 154 circular ground plots. I also extracted average elevation (Table C.1) 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. (Wilson, 2012). B2. Model building and selection process To avoid over fitting of the model and address the need for 10-15 observations per auxiliary variable, I divided the auxiliary variables into four groups (Mauro et al. 2019; Harrell et al., 1984). Within each variable group, following the top-down model-building strategy for glm models described by Zuur et al. (2009), the LiDAR-derived variables (Table C.1) 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 (Mauro et al., 2019). The variables that were not statistically significant (p-value > 0.05) and had a variance inflation factor (VIF) greater 10 were dropped sequentially (R. Sheridan et al., 2015; Sheskin, 158 2000; Zuur et al., 2009b). 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 (Sheskin, 2000; Zuur et al., 2009). 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 D.3) and graphs of the residual deviance (Zuur et al., 2009). I selected the biomass model with the lowest AIC and lowest residual deviance (Table C.2) as our final model. Our final model consisted of the mean height for all first returns above 2 m (zmean, p=1.4 x10-7), the height for the first return in the 5th percentile (zq5, p = 0.05) and crossed treatments as fixed effects (p=0.05) and unit level random effects (Table C.2). 159 Appendix C Additional tables for Chapter 5 Table C.1 LiDAR-derived and topography auxiliary variables were used for biomass modeling. Group LiDAR Variables Description Group A zmean Mean elevation for all first returns above 2m zq5 Elevation for first return in 5th percentile zq10 Elevation for first return in 10th percentile zq15 Elevation for first return in 15th percentile zq20 Elevation for first return in 20th percentile zq25 Elevation for first return in 25th percentile Group B zq30 Elevation for first return in 30th percentile zq35 Elevation for first return in 35th percentile zq40 Elevation for first return in 40th percentile zq45 Elevation for first return in 45th percentile zq50 Elevation for first return in 50th percentile zq55 Elevation for first return in 55th percentile zq60 Elevation for first return in 60th percentile zq65 Elevation for first return in 65th percentile zq70 Elevation for first return in 70th percentile zq75 Elevation for first return in 75th percentile Group D zq80 Elevation for first return in 80th percentile zq85 Elevation for first return in 85th percentile zq90 Elevation for first return in 90th percentile zq95 Elevation for first return in 95th percentile pzmean Percentage of first return above mean Elev (m) Elevation from digital elevation model @1 m Treatment X prescribe burn 160 Table C.2 Summary of the best-selected models from each auxiliary variable group. Highlighted in bold is the final model selected for above-ground biomass (Mgha-1) prediction. Variable group Model AIC Residual deviance p-value Group A Fixed effects = zmean+zq5+Treatment x prescribe fire Random effect = units 1577.9 1557.9 zmean (p=1.4 x 10-7) zq5 (p=0.00087) Group B Fixed effects = zq50+ Treatment x prescribe fire Random effect = units 1590.01 1572.0 zq50 (p=0.00034) Group C Fixed effects = zq75+ Treatment x prescribe fire Random effect = units 1582.9 1564.1 zq75 (p=4.7X10-6) Group D Fixed effects = zq95+ Treatment x prescribe fire Random effect = units 1581.9 1561.9 zq95 (p=8.1X10^-6) 161 C3. Summary of the prediction using biomass model The range of predicted values of overstory AGB for the cells was within the range of the observed values for circular ground plots (Table C.3, Figure D.1). The mean and median predicted values were very close to the mean and median observed values (Table C.3, Figure D.4). Overstory AGB predicted within the low structural diversity (LOD) units was low and consistent with observed overstory AGB (Figure D.5). Predicted overstory AGB in high structural diversity (HID) units and research natural areas (RNAs) was high compared to LOD and within the range of observed overstory AGB for HID and RNA units (Figure D.2 vs. Figure D.5). Table C.3 Summary statistics of the observed and predicted above-ground biomass (Mgha-1) values for circular plots (n=153) and cells (n=12,647). Above-ground biomass value (Mgha-1) Min Max Mean Median Observed (Circular plots) 4.33 327.53 101.71 83.48 Predicted (32x32 m cells) 10.33 381.65 101.02 89.03 162 Table C.4 Summary of RMSE and values for the spherical, exponential, and Gaussian models selected for Burned and Unburned halves of HID and LOD units. Unit Fire Model RMSE Sill LOD UNIT-39 Burned Gaussian 0.54 0.86 UNIT-39 Unburned Gaussian 0.71 1.04 UNIT-40 Burned Gaussian 0.72 0.93 UNIT-40 Unburned Gaussian 0.73 1.05 UNIT-43 Burned Spherical 0.75 1.03 UNIT-43 Unburned Spherical 0.74 1.00 UNIT-44 Burned Gaussian 0.64 1.15 UNIT-44 Unburned Gaussian 0.49 0.92 UNIT-45 Burned Gaussian 0.58 1.14 UNIT-45 Unburned Gaussian 0.65 1.07 HID UNIT-38 Burned Gaussian 0.65 0.92 UNIT-38 Unburned Gaussian 0.59 0.87 UNIT-41 Burned Spherical 0.74 0.84 UNIT-41 Unburned Spherical 0.76 0.98 UNIT-42 Burned Exponential 0.74 1.13 UNIT-42 Unburned Exponential 0.73 1.32 UNIT-47 Burned Gaussian 0.63 0.90 UNIT-47 Unburned Gaussian 0.68 2.81 UNIT-48 Burned Gaussian 0.70 0.92 UNIT-48 Unburned Gaussian 0.74 0.88 HID: high structural diversity, LOD: low structural diversity 163 Appendix D Additional figures for Chapter 5 Figure D.1 Histogram for above-ground biomass (Mgha-1) for 2016 ground data. The red line indicates mean, and blue indicates median. The red line indicates mean while blue line median. 164 Figure D.2 Predicted above-ground biomass (Mgha-1) using the final model and LiDAR derived metrics in 2016 at 32 m X 32 m pixel level. HID: high structural diversity treatment; LOD: low structural diversity treatment; RNA: research natural area; BMEF: Blacks Mountain Experimental Forest. 165 Figure D.3 Observed vs. predicted above-ground biomass (Mgha-1) values for the best model from each of the variable groups along with the final model. 166 Figure D.4 Predicted overstory above-ground biomass (AGB, Mgha-1) for 1 A) predicted overstory AGB at cell level, and B) observed overstory AGB for circular ground truthing plots. The blue line represents the median and red line represents the mean for AGB. 167 Figure D.5 Observed overstory above-ground biomass (AGB, Mgha-1) for circular ground truthing plots in 2016. 168 Figure D.6 Variogram models for research natural areas (RNA) in 1934 and 2016. 169 Figure D.7 Variogram models for low structural diversity units (LOD) in 1934 and 2016. 170 Figure D.8 Variogram models for high structural diversity units (HID) in 1934 and 2016. 171 Figure D.9 Differences in Moran’s I correlogram between 1934 vs. 2016 within different units under study in Blacks Mountain Experimental Forest. Dotted lines represent the lower and upper limit of the no significant spatial autocorrelation using 95% confidence envelopes for 1934 (red) and 2016 (blue) from Monte Carlo simulations. Points connected with solid lines indicate Moran’s I at a given lag-distance. Points above the dotted lines in upper part of envelope indicate a positive Moran’s I that is significantly different from 0, while points below the dotted lines in lower part of the envelope show negative Moran’s I significantly different from 0. Points within the dotted envelope shows no significant Moran’s I. 172 Figure D.10 Variogram models for Burned and Unburned halves of low structural diversity units (LOD) in 2016. 173 Figure D.11 Variogram models for Burned and Unburned halves of high structural diversity units (HID) in 2016. 174 Figure D.12 Differences in Moran’s I correlogram between 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 no significant spatial autocorrelation using 95% confidence envelopes for Burned (red) and Unburned (blue) from Monte Carlo simulations. Points connected with solid lines indicate Moran’s I at a given lag-distance (m). Points above the dotted lines in upper part of envelope indicate a positive Moran’s I that is significantly different from 0, while points below the dotted lines in lower part of the envelope show negative Moran’s I significantly different from 0. Points within the dotted envelope shows no significant Moran’s I.
Thesis/Dissertation
eng
Forestry
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Eskelson, Bianca N. I.
Spatial patterns of stand structure in Pinus ponderosa dominated forests within the Blacks Mountain Experimental Forest, California, USA, in 1930s and contemporary conditions
Text