DRIVERS OF 3-D CANOPY FUNCTIONAL DIVERSITY IN A TEMPERATE FOREST: MAPPING FUNCTION WITH HYPERSPECTRAL AND LIDAR REMOTE SENSING By Anthony Winston Bowman A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Geography—Master of Science 2024 ABSTRACT Leaf-level functional traits important to photosynthetic carbon (C) assimilation vary widely across biomes, between plant functional types, and coinciding with vertical structural variation in the canopy, which necessitates the use of tools capable of mapping this diversity over large scales. Advances in the spectral and spatial resolution in multi- and hyperspectral imagery (HSI) over the past three decades have allowed for increasingly fine resolution mapping of top- of-canopy (TOC) functional traits at the landscape scale from airborne spectrometers (Wang et al. 2020), while lidar data have been shown to effectively quantify structural attributes of forest canopies important to C assimilation (Dahlin et al. 2013; Asner 2015; Kamoske et al. 2019). However, these datasets are seldom spatiotemporally contiguous, thus limited studies have combined these datasets to capture a complete picture of canopy functional diversity. Here, I integrate field and remote sensing data from a mixed deciduous forest in Virginia, USA, to generate estimates of TOC percent and total-canopy (g/mGround 2) nitrogen content before comparing the spatial patterns and drivers of these two metrics. Conversely to trends in TOC percent N, we found that evergreen needleleaf canopies had on average higher total nitrogen (gN/m2) and foliar biomass (gLeaf/m2) than broadleaf deciduous canopies. Upon a closer analysis of the spatial trends and environmental drivers of TOC percent and total N between these plant functional types, we found that the change in plant health (i.e., NDVI) during an outbreak of the spring cankerworm (Palecrita vernata) was a strong predictor of total canopy N, suggesting this disturbance perturbs a complete understanding of whole-canopy function in this system. Overall, this thesis highlights the importance of local and regional-scale ecological understanding of plant form and function, as my results reinforce similar studies (e.g., Kamoske et al. 2021) but contradict long held views of differences in productivity between plant functional types. This thesis is dedicated to my partner, Brianna. I am forever grateful for your unwavering support, belief, and encouragement. iii ACKNOWLEDGEMENTS I would like to first extend thanks to the National Science Foundation’s Macrosystems Biology program which funded the latter half of my degree. I next want to thank my advisor, Dr. Kyla Dahlin, for her constant encouragement and belief through all adversity along this journey. I also thank my committee members, mentors, and those that helped to lay the foundation for this research, Drs. Lifeng Luo, Lars Brudvig, Aaron Kamoske, Ethan Theuerkauf, Randall Schaetzl, Scott Stark, and Shawn Serbin for their contributions, support, and guidance that made this thesis possible. Further, I would like to thank my fellow ERSAM Lab member, Dr. Adrianna Uscanga, and the 2024 Spectral Ecology (SPEC) School participants for their aid in data collection during the final stages of this research. Finally, I extend immense gratitude to my family and friends for their continuous motivation during the best and worst of times – your support made this dream a reality. iv TABLE OF CONTENTS Introduction ..................................................................................................................................... 1 Methods........................................................................................................................................... 6 Results ........................................................................................................................................... 17 Discussion ..................................................................................................................................... 26 BIBLIOGRAPHY ......................................................................................................................... 30 APPENDIX A ............................................................................................................................... 36 v Introduction Forest canopy functional and structural diversity drive resource use efficiency (RUE), production, and subsequently carbon (C) assimilation (Smith et al. 2002; Hardiman et al. 2011; LaRue et al. 2016; Gough et al. 2019). Research in this area has been an ongoing frontier in ecosystem ecology at local to global scales (Shugart et al. 2010; Niinemets et al. 2012; Stark et al. 2012; Atkins et al. 2018; Kamoske et al. 2021). However, understanding how functional and structural diversity interact to effect whole-canopy function (and thus C sequestration) has been limited due to the disconnect in spatiotemporal consistency between functional and structural measurements across spatial scales (Asner et al. 2015). As forests globally decline in species and functional diversity, ensuring effective modeling and management of functional and structural diversity has become even more crucial (Wilfarht et al. 2014). Leaf-level functional or response traits reflect the tradeoff between growth, production, and survival, causing predictable variation in leaf form and function along the leaf economic spectrum under given environmental conditions (e.g., resource availability and competition) (Violle et al. 2007; Diaz et al. 2016). While individual plant functional traits (e.g., leaf area) influence light acquisition and C sequestration at the leaf-level, canopy scale structural heterogeneity has been shown to drive large scale patterns of C assimilation (Hardiman et al. 2011). It has been hypothesized that structural heterogeneity creates diverse light environments allowing radiation to be absorbed throughout the canopy volume, maintaining production in vertical strata with little direct incoming light (Hardiman et al. 2013; Atkins et al. 2018). Forest processes in the form of net primary production and C sequestration are inherently tied to canopy light availability, which shows a strong relationship with vertical variation of leaf area and arrangement in traditional field studies (Brown & Parker 1994; Shugart et al. 2010; 1 Hardiman et al. 2011; Hardiman et al. 2013; Stark et al. 2015). This 3-D variation in leaf distribution at the canopy scale is caused by an array of historical and current forces not limited to succession, tree mortality, and impacts from pests (Hardiman 2011; Wilbur et al. 2021; Choi et al. 2023). At the stand and landscape scales, these factors create a diversity of light environments which are known to heavily influence variation in leaf-level functional traits (i.e., leaf mass per area, chemical composition) as plants maximize their photosynthetic capacity throughout the vertical extent of the canopy, with cascading effects on production and C assimilation (Field & Mooney 1986; Lefsky et al. 1999; Smith et al. 2002; Niinemets et al. 2012; Niinemets et al. 2015). Thus, it may be critical to use an integration of measurements on canopy-scale structure and leaf-level traits to truly describe forest function at the landscape scale and beyond. Though much work has been done to address the question of how forest structure impacts function at the traditional plot and landscape scales (Hardiman et al. 2013; LaRue et al. 2016), limited efforts have had the opportunity to use spatiotemporally comparable remotely sensed structural and functional information at the landscape and regional extent (Townsend 2003; Asner et al. 2015; Kamoske et al. 2021). Further, these studies have been localized to distinct ecosystem types, none of which are applicable to understanding the 3-dimensional functional variation in the temperate mixed forests of Appalachia. While leaf functional traits (particularly foliar N concentration) studied in forests undergoing successional change (Fowler et al. 2015; Wang et al. 2020) reflect latitudinal gradients of N deposition that may positively impact C assimilation, it is unclear if these landscape scale trends will persist if whole-canopy structural information from lidar is considered (Townsend et al. 2003; Thomas et al. 2010; Knyazikhin et al. 2013; Wang et al. 2020). Nitrogen (N) concentration has been effectively modeled using HSI in a variety of forest ecosystems, though no current studies have investigated the role of canopy 2 structure on large scale functional traits prescriptive of C sequestration in the Central Appalachians (Townsend 2003; Wang et al. 2020; Kamoske 2020). Many studies over the past three decades have reliably modeled and scaled foliar functional traits beyond the leaf-level using passive multispectral and more recently hyperspectral imagery (HSI) with in situ (or field-based) measurements to generate 2- dimensional estimates of canopy traits (Townsend et al. 2003; Dahlin et al. 2013; Wang et al. 2022). Recent advances in the spectral and spatial resolution in multispectral and subsequently HSI have allowed for increasingly fine scale mapping of top-of-canopy (TOC) foliar functional traits at the landscape scale from air and space-borne spectrometers (Lepine et al. 2016; Wang et al. 2022). Active remote sensing systems in the form of waveform and discrete-return light detection and ranging (lidar) have been effectively used to model the vertical structure of forest canopies in recent years, (Lefsky et al. 1999; Harding et al. 2001; Hardiman et al. 2011; Kamoske et al. 2019) and can yield important information about the distribution of leaves in the canopy and subsequent light availability through vertical canopy strata (Kane et al. 2010; Hardiman 2013). This allows for interpretation of canopy characteristics beyond TOC reflectance in the form of vertical structure, and thus can become complementary to HSI data for modeling the functional and structural traits of a canopy at larger scales than previously attainable with traditional in situ measurements alone (Asner et al. 2015; Kamoske et al. 2021). However, lidar data and HSI are seldomly recorded simultaneously and therefore a limited number of studies have allowed for integration of these measurements to understand the functional and structural characteristics of a given forest ecosystem (Asner 2015; Kamoske et al. 2021). The National Ecological Observatory Network (NEON; Kampe et al. 2010) works to circumvent this hinderance to understanding a complete picture of canopy function by operating 3 a hyperspectral imager measuring hundreds of wavelengths in conjunction with waveform and discrete return lidar sensors. These datasets offer the opportunity to study forest structure and function with studies leveraging this data to answer an array of questions about forest ecosystems, what drives their functional diversity, and how this may be impacted by disturbance (Kamoske et al. 2021; Choi et al. 2023). In the mixed deciduous forests of the Appalachian Mountains, disturbances from historical overlogging, development, species invasions, and degradation of the understory by animals (particularly white-tail deer (Odocoileus virginianus)) diversify light environments throughout the canopy volume by influencing stand density, canopy clumping, among other structural characteristics which shape functional and structural diversity at various spatial scales (Tremblay et al. 2004; Small & Chamberlain 2015; Wilbur et al. 2021; Choi et al. 2023). The Central Appalachians in southwestern Viriginia were deforested throughout the 19th century due to logging for timber production and clear-cutting for cattle pastures, which coincided with the with the sharp decline of the locally dominant American Chestnut (Castanea dentata) and emergence of red oak (Quercus rubra) as the dominant species in this region (Redmond et al. 2012; Small & Chamberlain 2015; Wilbur et al. 2021). The recovery of these forests has resulted in a range of broadleaf deciduous and evergreen needleleaf stands assorted along elevational gradients typical of this region, where pine species and eastern hemlock (Tsuga canadensis) dominate valley floors while successional oak-hickory stands persist along more well-drained upland sites. Through this functional and structural diversity secondary forests have been shown to increase or maintain levels of production and C sequestration as they age (Hardiman et al. 2013), though the long-term effects of species diversity decline due to increased N deposition, 4 species invasions, and development over the last century remain unclear (Butler et al. 2015; Thomas et al. 2010; Fowler et al. 2015). Temperate forests in this region of the world are ideal model systems to advance a methodology for measuring whole-forest function using a fusion of tools as they are often complex in terms of species and plant functional diversity, which allows for comparison of forest function across plant communities dominated by different functional types and environmental gradients at finer scales than current satellite products offer (Zhang et al. 2017). Due to the role of temperate Appalachian forests in regional C dynamics and global C storage and productivity, the impact of vertical structural heterogeneity on forest function should be quantified to inform C cycling parameterizations in models generated from remotely sensed data (Liu et al. 2006; Pan et al. 2011; Fowler et al. 2015; Zhang et al. 2017; Canham et al. 2024). In this study, we investigate whether the incorporation of vertical structural information can elicit more robust models of canopy functional traits than top-of-canopy information alone by leveraging a fusion of field, lidar, and hyperspectral data. Following Kamoske et al. (2021), the questions I seek to address in this research are: 1) How can a fusion of lidar, hyperspectral, and in situ data be used to accurately model 3D canopy functional traits in an Appalachian forest diverse in species, plant functional types (PFTs, e.g., needleleaf vs deciduous broadleaf trees), and topography? 2) How do the spatial trends of these traits differ from the TOC to whole canopy levels? and 3) What are the environmental drivers of top of canopy and total N in a heterogeneous central Appalachian forest? 5 Methods Site description Mountain Lake Biological Station (MLBS) is nested within the Appalachian Mountain chain atop Salt Pond Mountain in southwestern Virginia, USA. MLBS has a cool, continental climate with high seasonality and low average temperature at 8.8°C due to harsh winters typical at this elevation (Wilbur et al. 2021; NEON Site Report). Elevation at this site ranges from 800 to 1,329 meters, with a mean canopy height around 18m. Mean annual precipitation is roughly 1277mm, with sporadic dry and wet spells characteristic of the central Appalachians. MLBS is surrounded by a mixed oak-hickory forest consisting of northern red oak, white oak (Quercus alba), sweet birch (Betula lenta), eastern hemlock (Tsuga canadensis), and shagbark hickory (Carya ovata), with stands of other gymnosperms such as red cedar (Juniperus virginiana) also present. At higher elevations, canopies mainly consist of red maple (Acer rubrum) and white oak while at lower elevations, more fertile soils support species such as blackgum (Nyssa sylvatica) and eastern white pine (Pinus strobus). The forest understory is composed of multiple species of ferns, Rhododendron species, and mountain laurel (Kalmia latifolia). The study area used here (Fig. 1) is a 2 x 2 km square centered on the NEON eddy covariance flux tower and surrounding MLBS itself, which has been owned and managed by the University of Virginia since its inception in 1936 as a field research site. MLBS has served as a terrestrial NEON field site in Domain 07 since 2017, with three separate plots adjacent to existing MLBS research sites and actively managed U.S. Forest Service land within the Jefferson National Forest (Sarvis et al. 1993; Wilbur et al. 2021). The following year in 2018, an outbreak of the Spring cankerworm (Palecrita vernata; hereafter cankerworm), occurred throughout the study area (NEON 2021b; Choi et al. 2023). The larva of this generalist pest defoliates deciduous 6 vegetation throughout the canopy volume and understory by rapidly consuming leaves, which can drastically change the light environment and short-term structure of canopies in affected stands (Fig. 2.; Darr & Coyle 2021). A recent study by Choi et al. (2023) found that in all 23 NEON plots (40 x 40 m) examined at MLBS, roughly 60% of oak species showed moderate-high levels of defoliation with the anticipated impacts of gap fraction increase and canopy density (or leaf area density (LAD; mLeaf 2/mGround 2)) reduction which could impact accurate estimations of canopy structural traits otherwise interpretable in absence of this defoliation event. Figure 1. Map of the study area extent around MLBS used in this project centered around the NEON flux tower, denoted as a blue diamond north of the research station, and MLBS-owned land as a dashed line excluding roads. 7 Figure 2. Example of defoliated broadleaf stand in August during the 2018 cankerworm outbreak. (photo credit: Aaron Kamoske) Field Sampling & Lab Methods To generate models of within-canopy functional traits from both remotely sensed and field data, leaf samples were collected between early May and mid-June 2018, proximal and coinciding with the NEON Airborne Observation Platform (AOP) flight at MLBS on June 15th, 2018, and, coincidentally, during the cankerworm outbreak (NEON Site Report 2021b). Throughout the vertical extent of the canopy, foliar samples (n = 111) were collected from the dominant species within the canopy using a pole pruner or BigShot Line Launcher (SherrillTree, Greensboro, North Carolina, USA). Each foliar sample height was estimated using a laser range finder. 8 During collection, each sample location was documented with a Trimble GEO7x GPS (Trimble, Sunnyvale, California, USA) and differentially corrected using Trimble’s GPS Pathfinder Office software afterwards. Upon collection, each leaf sample was wrapped in a damp paper towel and placed in a cooler with ice packs prior to collecting reflectance measurements the same day in our mobile field laboratory. Three reflectance measurements were taken from different leaves in each sample using an SVC HR-1024i Spectroradiometer equipped with an LC-RP-Pro leaf clip foreoptic (Spectra Vista Corporation, Poughkeepsie, New York, USA). The SVC HR-1024i Spectroradiometer collects data in wavelength range of 340 – 2,500nm, with roughly 1050 bands binned at approximately 2nm. Between each sample, the instrument was recalibrated using a white Spectralon panel. For broadleaf samples, individual leaves were placed flat directly into the leaf clip, whereas for needleleaf samples mats were created by aligning needles vertically aside one another and taping their ends together prior to being placed in the leaf clip. At least 500 mg of leaf material was collected from each sample using scissors sterilized between collections, placed on a flatbed scanner and analyzed for area using the ImageJ software. After imaging, the leaf material collected from each sample was placed in a paper coin envelope and dried at 70°C for at least 48hr. Leaf mass per area (LMA; g/mL 2, where mL is meters of leaf area) was subsequently calculated using area measurements from the ImageJ software and dry weight (g) for each sample. A subset of samples (n=34, ~30%) was then collected, redried, and ground using a ball mill (2000 Geno Grinder; SpexSample Prep, Cridersville, Ohio, USA) into a fine powder. This subset was sent to Brookhaven National Laboratory (Upton, New York, USA) for determination of elemental N content (g N/gleaf, or percent N) and sample C:N ratio. From each sample, 1.50-2.50 mg were weighed using 0.1 mil tin foil vials (AX26DR; Mettler Toledo, Columbus, Ohio, USA) and analyzed using a 2400 Series II CHNS/0 Analyzer (Perkin Elmer, 9 Waltham, MA, USA) following Kamoske et al. (2021). To retrieve a field-based measure of vertical leaf density and arrangement and support future comparisons of whole-canopy traits across ecoregions, we followed the protocol described in Kamoske et al. (2019) and collected 44 upward facing hemispherical photographs at this site. Hyperspectral Data The hyperspectral imagery used to generate the top-of-canopy (TOC) functional trait data employed here was collected by the NEON AOP on August 5, 2017, the summer prior to the cankerworm outbreak. The NEON AOP is equipped with a hyperspectral sensor that measures 426 distinct bands in the visible-SWIR wavelength range (380-2500nm) in 5nm increments at 1m spatial resolution (Kampe et al. 2010). Using hyperspectral reflectance and leaf samples collected in the sun-lit portions of the canopy crowns, a partial least squares regression (PLSR) approach was used by the Townsend Lab at the University of Wisconsin-Madison (Wang et al. 2021) to model and generate maps of top-of-canopy LMA and N content (mg N/gLeaf; converted to a percentage here for comparability with other studies) for the AOP flightline extent around MLBS (Data publicly available at https://uwmadison.app.box.com/s/xll3ljh5ar8rh0geztywd5w2qkfgzcib). These maps were leveraged here in lieu of deriving our own maps of TOC LMA and N content for MLBS since 82 samples were collected to model TOC functional traits at MLBS by Wang et al. (2021), and only 34 TOC samples were collected during our field campaign. Moreover, during the cankerworm outbreak leafless patches of the canopy would allow understory plants to influence the airborne reflectance measurements. In order to match the lidar data discussed in the next section, the 1 x 1 m functional trait maps were aggregated to 25m2 by taking the average of all pixels in each 5 x 5 m extent. All processing was done in the R programming language. 10 Lidar data processing Discrete return lidar point cloud data were acquired for 2018 NEON AOP flightlines over MLBS and rasterized into 5 x 5 m voxels for each 1 m vertical slice of the canopy. After leveling the lidar point clouds in each cell so that the elevation of all ground returns was equal to 0, LAD (mLeaf 2/mG 3 , where mG 3 is a volume of airspace above the ground) was calculated for all voxels within the AOP flightlines (Sabol et al. 2014; Kamoske et al. 2019). LAD here was calculated relying on theory from MacArthur and Horn (1969), where LAD is equal to the total number of lidar pulses that enter and exit each voxel that have at minimum of one ground and canopy return using the canopyLazR package (Stark et al 2012; Kamoske et al 2019; Kamoske et al 2021). To avoid noise from understory vegetation and sharp topographic variations on the ground (Fig. 3) affecting interpretation of canopy structural traits, the bottom 5 meters of the canopy were removed before further processing. To calibrate lidar LAD estimates to field-based measurements, hemispherical upward facing fisheye photographs collected in the field were processed for LAI. All fisheye photographs were taken with the Canon Rebel T-6 and Bower 8mm F3.5 Ultra-Fast Lens and processed using the hemispeR package, which utilizes an angular gap fraction approach to calculate LAI, effective leaf area index (LAIe), and other canopy structural attributes (Canon Inc. (Ota, Tokyo, Japan), Bower Inc. (New York, New York, USA); Macfarlane et al. 2012; Chiannuci & Cutini 2013; Chiannuci & Macek 2023). The hemispheR package generates estimates of LAI by relying on the blue region of the visible spectrum in RGB photographs and binarizes each into black (canopy) or white (sky) based on an image-specific calculated reflectance threshold in the blue band. This method takes advantage of the stark difference in 11 diffuse and direct radiation in the blue band between sky and canopy pixels in upward facing photos. Figure 3. Rock cavern near Bear Cliff (VA). Each fisheye image was binarized into sky or canopy and divided into five zenith rings and eight azimuth segments at a maximum of 27° to mimic estimates from an LAI 2000/2200 Plant Canopy Analyzer, and replicate the extent of the AOP swath width, respectively. (LI-COR study; Chiannuci et al. 2015; Chiannuci & Macek 2023). From the binarized and segmented fisheye image, gap fraction estimates were then derived to calculate LAI and effective leaf area index (LAIe). For consistency due to the inability to discern woody material from leaves in the lidar returns, LAIe values (hereafter LAI) extracted from the hemispherical photos were used and 12 compared against lidar-derived LAI (Fig. A.1.) to derive a Beer-Lambert extinction coefficient (k = .53; Fig. A.2.). The Beer-Lambert coefficient in this context is the slope of the relationship between the lidar- and digital hemispherical photography (DHP) derived LAI estimates, which is then multiplied by the voxel height to account for the exponential decrease in light penetration as it pulses penetrate further into the canopy (Parker et al. 2001, Stark et al. 2012; Kamoske et al. 2019; Kamoske et al. 2021). From these calibrated LAD estimates, we calculated 27 structural attribute maps (Table A.1.) representing different aspects of the canopy architecture for each 1m slice of the canopy. Within-Canopy Modeling To build a model for within-canopy functional traits at this site, we first extracted data from the 26 lidar-derived structural attribute maps and the TOC LMA and N (%) maps at pixels encompassing our 111 leaf sample locations. To keep consistency with the hyperspectral and lidar data, we removed all samples below 5 m in the canopy and similarly those at the top of the canopy, as similar samples were used to train the PLSR model for the TOC traits. This resulted in a final set of 41 samples, which were split into ~80% training (n=33) and testing 20% (n=8) datasets. To prevent multicollinearity among predictors in our model, we removed all variables that were too correlated (R > 0.5) which resulted in a retention of 11 predictor variables. These final predictors were used to train an ordinary least squares regression (OLS) model to predict within-canopy LMA, from which the best set of predictor variables was selected using a backwards stepwise model selection approach using the Akaike Information Criterion (AIC). Preliminary OLS regression results showed that only TOC LMA was a statistically significant predictor of within-canopy LMA when included in models along with the 26 structural attribute rasters (Appendix A). Thus, in lieu of using modeled within-canopy LMA values from the OLS 13 regression equation which would have likely added uncertainty, we used TOC LMA values and calculated total-canopy N (g/mG 2, where mG 2 is meters of ground) using the below equation: 𝑁𝑡𝑜𝑡 = 𝑁𝑇𝑂𝐶 × 𝐿𝑀𝐴𝑇𝑂𝐶 × 𝐿𝐴𝐼 where h is the maximum height of each column of voxels, NTOC is the TOC N value at a given pixel (converted to percent N to be comparable with Kamoske et al. 2021) from Wang et al. (2020) and i represents individual voxels in the canopy column beginning at 5 m, as lidar pulse returns diminish or become unreliable due to noise from understory plants and topographic variation (Kamoske et al. 2019; 2021). In addition to Ntot, we used the same equation to calculate foliar biomass (g/m2) for the entire canopy volume by withholding the TOC N values and summing the product of LMA and LAD for each 1m slice of the canopy. To understand whether the spatial patterns of TOC percent N and total canopy N correlate, potentially allowing for extrapolation of whole canopy traits from remotely sensed TOC traits, we extracted 10,000 points from each raster at the 5 x 5 m resolution. These samples were used to fit variograms for comparison of the two metrics, and gauged their respective spatial autocorrelation using the range, nugget, and sill of the variograms. Planet data processing To map deciduous vs evergreen forests and the impact of the cankerworm outbreak, we acquired Planet Labs (PlanetScope; Planet Team 2020) 3m mosaicked 4-band (red, blue, green, and near- IR) surface reflectance scenes at MLBS for 2016, 2018, and 2019 (cloud-free images were not available for 2017). Planet operates a collection of 430+ small cube satellites (called Cubesat Doves) which work in unison to offer high spatiotemporal resolution at 3m spatial and on- average daily revisit time, respectively. To delineate the PFTs within the study area, we derived NDVI (red – NIR / red + NIR) from leaf-off reflectance data (early March 2016) and set a 14 threshold of 0.6 to classify each pixel into deciduous (< 0.6) or evergreen (> 0.6). Additionally, we calculated NDVI for leaf-on data for pre- and post-cankerworm outbreak (August 2016 and 2019) to assess the extent of defoliation due to this event and see how it may affect estimates of TOC and total N in broadleaf deciduous stands. Environmental driver analysis In order to address our second and third research questions of how top-of-canopy and total- canopy functional traits differ spatially and in their respective drivers, we used Moran’s I to estimate spatial autocorrelation and a Monte Carlo simulation approach to understand the respective influences of the underlying abiotic (e.g., slope) and biotic conditions (i.e., cankerworm outbreak) as drivers of total (g/m2) and top-of-canopy percent N. To generate a dataset representative of the abiotic variables that may be driving spatial patterns of canopy N at this site, we generated 11 raster layers representing topographic characteristics of MLBS. Topographic characteristics (e.g., topographic position index, flow accumulation) were derived using a digital terrain model generated from the 5 x 5 m lidar data (Table A.2.). Due to the lack of variation in geologic substrate within the study area (i.e., a continuous substrate type throughout) no geologic variables were used. A Monte Carlo test with 1,000 simulations was then conducted to calculate an average R2 and Moran’s I of model residuals, and a distribution of model coefficients. During each simulation, 10,000 points were extracted from the predictor variable raster layers to build the models. After standardizing all non-binary variables (Gelman 2006; mean = 0, standard deviation = 0.5) to allow for comparison across predictors, we generated regression models for top-of-canopy percent N and total canopy N individually. To avoid multicollinearity among the predictor variables, we assessed the correlation (Pearson’s R) between each and removed the variable least correlated with canopy N when two predictors had 15 an R ≥ 0.5 for each regression model and simulation. After filtering out highly correlated variables with little impact in explaining variation in total or TOC percent N, an initial OLS regression equation was developed using the remaining predictors. Following Kamoske et al. (2021), we used backward stepwise AIC variable selection to reduce our set of predictors to only those found to be statistically significant (P < 0.05). This final subset of predictors was then used in a final OLS regression model for each simulation. We assessed spatial autocorrelation using Moran’s I on the residuals of the individual final models. All processing was completed with ArcGIS Pro and the R programming language (ESRI 2023; R Core Team 2021). 16 Results Modeling within-canopy functional traits with lidar and hyperspectral data The final OLS regression model for predicting within-canopy LMA with the lowest AIC score consisted of three components: top-of-canopy (TOC) LMA, within-canopy rugosity, and sample height. Due to TOC LMA being the only statistically significant predictor variable in final model (p < 0.01), within-canopy rugosity and sample height were removed, and TOC LMA alone was used to model total canopy N (g/mG 2 2; Fig. 5). Similarly, we calculated LAI (mL /mG 2; Fig. 4) using the summation of LAD for each 1m slice in every pixel. LAD estimates here were calibrated with LAI estimates from in situ hemispherical photos to derive a Beer-Lambert extinction coefficient (k = .53) to account for diminishing lidar returns as they penetrate further into the canopy. Using the same calibrated LAD rasters for each 1m slice of the canopy and LAI for the study area, we masked all pixels from the structural attribute rasters that were found to be outliers in the distributions of LAI and canopy height values outside the published ranges for this ecoregion (Albaugh et al. 2020). This caused 5,837 pixels (~27%) to be removed that either had an LAI > 6 or canopy height > 40, which exceed the respective ranges published from past studies at this site (Wilbur et al. 2021; Prior et al. 2022). 17 Figure 4. Map of leaf area index (LAI; mLeaf 2/mGround 2). Figure 5. Maps of top-of-canopy N (%; Wang et al. 2020) and total canopy N (g/mG 2) for MLBS. 18 Using this LAI raster and TOC LMA, we also calculated foliar biomass (g/mG 2; Figure 6). Foliar biomass was calculated in this way to keep consistency with how we calculated total canopy N as only TOC LMA was a strong predictor of in situ samples regardless of location in the canopy. We then removed extreme outliers from our foliar biomass estimate using Tukey’s test (k =3), which resulted in 898 cells (~0.6%) being removed from the final raster. Figure 6. Map of foliar biomass (g/mG 2). 19 We assessed the relationship between top-of-canopy percent N and total canopy N after normalizing both distributions (mean = 0, SD = 1) and found a weak but significant relationship between the two metrics (R2 = 0.01, P < 0.05). To further understand the comparative spatial patterns between top-of-canopy percent and total N at MLBS, we calculated variograms for each dataset (Fig 7.). For the top-of-canopy percent N map, we found that this metric is slightly more autocorrelated (Moran’s I = .68) than the map of total canopy N (Moran’s I = .55). Both metrics had similar autocorrelation among residuals, with Moran’s I of the top-of-canopy percent N residuals equal to 0.014 with a standard deviation of 0.001, and total-canopy N with a Moran’s I of 0.015 and standard deviation of < 0.001. These values indicate minor autocorrelation in top- of-canopy and total-canopy functional traits that our models may not capture with the current suite of predictor variables, though this study was focused on the comparative influence of abiotic and biotic variables of these two metrics rather than using this regression in a predictive manner. A closer assessment of the variogram components for TOC (%) and total N reveals that TOC percent N samples were spatially autocorrelated up to about 97 meters, whereas total N is spatially autocorrelated at nearly over twice this distance (Fig. 7.; range = 196 meters). Unlike the distinct difference in range of spatial autocorrelation among the two metrics, partial sill (Psill) and nugget values for top-of-canopy percent N (Psill = .34, nugget = .64) are comparable with total canopy N (Psill = .25, nugget = .70) indicating similar clustering of values between the two datasets. 20 Figure 7. Variograms for normalized datasets (mean = 0, SD = 1, 10,000 samples each) of top- of-canopy percent N and total canopy N (g/mG 2). Psill = partial sill. Environmental driver analysis High and low elevation regions of the study area exhibited distinct patterns between top-of- canopy percent and total canopy N (Fig. 5), as oak-hickory upland stands give way to hemlocks and eastern white pine in more nutrient rich lowlands. As such, plant functional type (PFT) and elevation seemed to influence the contrasting spatial patterns between leaf-level and whole- canopy traits. Elevation was loosely related to both top-of-canopy (TOC) percent N (R2 = 0.04, P < 0.01) and total canopy N (R2 = 0.06, P < 0.01). Understandably, the relationships between PFT and metrics of canopy N were also minor, as functional types and elevation both follow trends from deciduous to more evergreen stands along an elevational gradient in this system (aside from a planted pine stand directly east of the station itself). PFT was related to TOC percent N (R2 = 0.06, P < 0.05), but showed a weaker relationship to total canopy N (R = 0.004, 21 P < 0.05). Additionally, the difference in NDVI between 2016 and 2018 in upland areas with deciduous canopies showed a clear visual decline likely due to the fall cankerworm outbreak (NEON 2021b; Choi et al. 2023), which may have influenced these relationships. To investigate these relationships more thoroughly, we assessed the relative influences of topography, disturbance, and functional type on top-of-canopy percent and total canopy N (predictors listed in Table 1). Models for TOC percent N had a mean R2 = 0.10 and standard deviation < 0.05, with seven variables appearing in more than 20% of models (Table 1.), three variables appearing in no models, and four variables appearing in all models. Elevation (DTM) was the only major predictor (coefficient = .11) of top-of-canopy N with a positive coefficient. No other predictors had strong positive or negative influence from the compiled results (coefficient > 0.1 or < -0.1). The models for total canopy N performed negligibly better, with an R2 = .11 and standard deviation also less than 0.05. For total canopy N (g/mG 2), eight variables appeared in more than 20% of models, two variables appeared in no models, and three variables appeared in all models (Table 1.). Easting (distance east from collection boundary) was a strong negative topographic predictor (coefficient = -0.99) of total canopy N, along with topographic position index (-0.32), topographic wetness index (-0.33), plant functional type (-0.31) and solar radiation at the summer solstice (-0.24), whereas flow accumulation, slope, and aspect were major positive predictors (coefficient > 0.1). The difference in NDVI from 2018 to 2016 was the stronger positive predictor of total canopy N (coefficient = 1.0). 22 Figure 8. Monte Carlo simulation coefficients from standardized predictor variables (mean = 0, sd = 0.5) that appeared in at minimum 20% of models and are statistically significant (P < 0.05) for top-of-canopy percent N or total N (g/mG 2). 23 Figure 9. Maps of total canopy N (g/mG 2) and within-canopy nitrogen (g/mG 3) profiles for two dominant canopy species at Mountain Lake Biological Station: (a) white oak total foliar N: 9.49 g/mG 2) and (b) Eastern hemlock (total foliar N: 10.01 g/mG 2). 24 Table 1. Mean standardized coefficients (mean = 0; SD = 0.5), standard deviation of coefficients, and percentage of models’ predictors were present in Monte Carlo simulations Variables Aspect DTM (Elevation) Easting Northing Top-of-Canopy Total Canopy Mean coefficient Standard deviation Models present (%) Mean coefficient Standard deviation Models present (%) 0.0459 0.006 100 0.0629 0.110 4.3 0.1174 0.014 100 N/A N/A 0 -0.0334 0.006 99.7 -0.9987 0.056 100 N/A N/A 0 -0.1519 0.032 46.6 65.3 Flow accumulation -0.1391 0.007 14.3 Slope 0.0839 0.006 58.6 0.1734 0.2658 0.041 N/A 0 Topographic Position Index (TPI) Topographic Wetness Index (TWI) Topographic Roughness index (TRI) Solar radiation at summer solstice Solar radiation at winter solstice Difference in NDVI (2018-2016) Plant functional type 0.0549 0.006 100 -0.3205 0.063 99.8 -0.0133 0.013 35.1 -0.3354 0.069 100 -0.0074 0.014 4.5 0.1361 0.024 7.5 N/A N/A 0 -0.2424 0.058 94.7 -0.0377 0.005 41.4 N/A N/A 0 N/A N/A 0 1.001 0.073 100 -0.2184 0.017 100 -0.318 0.102 38.3 25 Discussion We leveraged a combination of field, hyperspectral, and discrete return lidar data to understand how a fusion of these tools can be applied to understanding whole-canopy function in a temperate, high-elevation Central Appalachian forest. Compared to foundational work fusing hyperspectral and lidar data to model and compare whole canopy traits in Southern Appalachia (Kamoske et al. 2021), we found that lidar-derived structural attributes such as canopy rugosity and leaf area density we not reliable predictors of within-canopy functional traits (i.e., LMA) as our field samples only correlated well with top-of-canopy estimates out of all predictor variables. Thus, we adapted our methodology accordingly to address our first research question of how an integration of data sources can be used to model whole-canopy functional traits in this system by simplifying our equation for total canopy N (g/m2), in contrast to previous work in taller, more vertically complex canopies (Kamoske et al. 2021). In addressing our second research question of how top-of-canopy and whole-forest functional traits compare to each other spatially, our results suggest maps of top-of-canopy percent and total canopy N (g/m2) are divergent in spatial patterns at this site but similarly spatially autocorrelated, with areas of low top-of-canopy nitrogen coinciding directly with high total canopy nitrogen between the two maps (Fig. 5). Similarly, top-of-canopy percent and total canopy N had distinct environmental drivers, though comparatively similar strengths for the overall regression models in predicting the two metrics. Because of the divergence of top-of-canopy and total canopy N between broadleaf and needleleaf evergreen stands at this site, plant functional type was presumably driving top-of- canopy and whole-canopy functional traits but had a minor influence on either metric from our compiled Monte Carlo simulation results. Therefore, it was clear from our environmental driver 26 analysis that due to the strong relationship between NDVI difference pre- and post-cankerworm outbreak with total canopy N, disturbance was driving the trends in whole-canopy function. The 2018 cankerworm outbreak likely caused underestimation in total canopy N (g/mG 2) values in broadleaf stands where defoliation was mild to severe (Choi et al. 2023) due to an underestimate of LAI. Because the NDVI difference between pre-defoliation and post-defoliation was a strong predictor of whole-canopy nitrogen but did not appear in any models of top-of- canopy N, it appears that this underestimation of LAI may have reduced total canopy N predictions particularly in broadleaf stands where defoliation was most severe. As disturbances (e.g., severe wind events, insect outbreaks, among others) that may confound estimations of whole-canopy traits become more common along with warming temperatures and shifts in community diversity and function across latitudes (Wilfarht et al. 2014), periodic modeling of whole-canopy traits will be essential for efficient forest management practices due to the divergence from spatial trends in top-of-canopy estimates of foliar N found in this and previous studies (Kamoske et al. 2021). While regular modeling of whole-canopy traits using this methodology may be inefficient in terms of time and resources, our results here demonstrate that modeling 3-D functional traits for a single year, or over two years with the integration of lidar and hyperspectral data from consecutive leaf-on NEON AOP flights, may not fully capture whole-canopy traits when climatic and biotic disturbances have occurred. At the leaf scale, needleleaf evergreen leaves are associated with low N concentrations and low productivity, while broadleaf deciduous leaves are associated with higher N concentrations and higher productivity (Diaz et al. 2016). Similarly, at continental scales we tend to associate broadleaf deciduous and mixed forests with higher productivity, relative to dryland and high latitude needleleaf evergreen forests (Zhang et al. 2016). Kamoske et al. 2021 show that 27 in a single landscape, stands that contain or are dominated by needleleaf evergreen trees tend to have higher LAI, higher leaf biomass, and higher total canopy N concentration (Fig. 9), suggesting they may be functionally similar compared to broadleaf deciduous dominated stands. Whether this observation translates to differences in productivity is difficult to measure, as eddy covariance tower locations rarely vary only in PFT distribution and no other confounding factors. While the results presented here may reflect these trends found in Southern Appalachia (i.e., Kamoske et al. 2021), the defoliation event coinciding with lidar data collection heavily impacted estimates of canopy structural traits (e.g., LAI) in broadleaf stands, which prevents a complete comparison between these two sites and clouds any concrete conclusions about how whole-canopy function across forests at MLBS is influenced by abiotic or biotic (i.e., plant functional type) conditions. Decreasing plant biodiversity due to shifts in climatic regimes (Wilfahrt et al. 2014), species invasions (Simberloff & von Holle 1999; Mollot et al. 2017), and disturbance (Dale et al. 2001; Tremblay et al. 2004) has been shown to negatively affect forest functional and structural diversity within the canopy, characteristics which are known to positively impact net primary production and carbon assimilation (Hardiman et al. 2011; Hardiman et al. 2013). The disconnect between current understanding of productivity between PFTs (i.e., evergreen vs deciduous) at varying scales shown by Kamoske et al. (2021) and reflected here despite the clear impact of using non-contemporaneous data affected by disturbance, demonstrates that more work is needed within the range of temperate forest ecosystems to fully disentangle how leaf-level and canopy- scale functional traits influence overall productivity more broadly. In this study, we show that while disturbance in distinct communities (i.e., broadleaf deciduous stands) can perturb accurate estimations of total-canopy function in a single year, 28 divergent patterns exist between canonical estimates of canopy function (i.e., TOC percent N) and total-canopy function when structure is incorporated in undisturbed canopies. Thus, it is important to understand how forests are changing not only out of concern for biodiversity, but how broad scale whole-forest function may be altered in the future due to disturbance along with structural, land use, and species composition changes in these ecosystems with cascading effects on carbon storage (LaRue et al. 2016; Lawrence et al. 2019). 29 BIBLIOGRAPHY Albaugh, T.J., C.A. Maier, O.C. Campe, M.A. Yáñez, E.D. Carbaugh, D.R. Carter, R.L. Cook, R.A. Rubilar, and T.R. Fox. 2020. Crown architecture, crown leaf area distribution, and individual tree growth efficiency vary across site, genetic entry, and planting density. Trees 34, 73–88. Asner, G.P., R.E. Martin, C.B. Anderson, and D. E. Knapp. 2015. Quantifying forest canopy traits: Imaging spectroscopy versus field survey. Remote Sensing of Environment 158: 15- 27. Atkins, J.W., R.T. Fahey, B.H. Hardiman, and C.M. Gough. 2018. Forest canopy structural complexity and light absorption relationships at the subcontinental scale. Journal of Geo- physical Research: Biogeosciences 123: 1387–1405. Butler, P.R., L. Iverson, F.R. Thompson III, L. Brandt, S. Handler, M. Janowiak, P. D. Shannon, C. Swanston, K. Karriker, J. Bartig, S, Connolly, W. Dijak, S. Bearer, S. Blatt, A. Brandon, E. Byers, C. Coon, T. Culbreth, J. Daly, W. Dorsey, D. Ede, C. Euler, N. Gillies, D.M. Hix, C. Johnson, L. Lyte, S. Matthews, D. McCarthy, D. Minney, D. Murphy, C. O’Dea, R. Orwan, M. Peters, A. Prasad, R. Cotton, J. Reed, C. Sandeno, T. Schuler, L. Sneddon, B. Stanley, A. Steele, S. Stout, R. Swaty, J. Teets, T. Tomon, J. Vanderhorst, J. Whatley, and N. Zegre. 2015. Central Appalachians forest ecosystem vulnerability assessment and synthesis: a report from the Central Appalachians Climate Change Response Framework project. 2015. Gen. Tech. Rep. NRS-146. Newtown Square, PA: U.S. Department of Agriculture, Forest Service, Northern Research Station. 310 p. Brown, M.J., and G.G. Parker. 1994. Canopy light transmittance in a chronosequence of mixed- species deciduous forests. Canadian Journal of Forest Research 24: 1694–1703. Canham, C.D., Murphy, L., & Hansen, W.D. 2024. Successional dynamics of carbon sequestration in forests of the eastern United States. Ecosphere 15(4): e4838. Chianucci, F., C. Macfarlane, J. Pisek, A. Cutini and R. Casa. 2015. Estimation of foliage clumping from the LAI-2000 Plant Canopy Analyzer: effect of view caps. Trees 29: 355– 366 Chianucci, F., J. Zou, P. Leng, Y. Zhuang, and C. Ferrara. 2019. A New Method to Estimate Clumping Index Integrating Gap Fraction Averaging with the Analysis of Gap Size Distribution. Canadian Journal of Forest Research 49(5): 471–79. Chiannuci, F. and M. Macek. 2023. hemispheR: an R package for fisheye canopy image analysis. Agricultural and Forest Meterology 336; 109470. Choi, D.H., LaRue, E.A., J.W. Atkins, J.R. Foster, J.H. Matthes, R.T. Fahey, B. Thapa, S. Fei, and B.S. Hardiman. 2023. Short-term effects of moderate severity disturbances on forest canopy structure. Journal of Ecology 111(9): 1866-1881. 30 Dahlin, K.M., G.P. Asner and C.B. Field. 2013. Environmental and community controls on plant canopy chemistry in a Mediterranean-type ecosystem. Proceedings of the National Academy of the Sciences 110: 6895-6900. Dale, V.H., L.A. Joyce, S. McNulty, R.P. Neilson, M.P. Ayres, M.D. Flannigan, P.J. Hanson, L.C. Irland, A.E. Lugo, C.J. Peterson, D. Simberloff, F.J. Swanson, B. J. Stocks, B.M. Wotton. 2001. Climate Change and Forest Disturbances: Climate change can affect forests by altering the frequency, intensity, duration, and timing of fire, drought, introduced species, insect and pathogen outbreaks, hurricanes, windstorms, ice storms, or landslides. BioScience 51(9): 723–734. Darr, M.N., and D.R. Coyle 2021. Fall cankerworm (Lepidoptera: Geometridae), a native defoliator of broadleaved trees and shrubs in North America. Journal of Integrated Pest Management 12(1), 1–9. Dıaz, S., J. Kattge, J.H.C. Cornelissen, I.J. Wright, S. Lavorel, S. Dray, B. Reu, M. Kleyer, C. Wirth, I. C. Prentice, E. Garnier, G. Bönisch, M. Westoby, H. Poorter, P.B. Reich, A.T. Moles, J. Dickie, A.N. Gillison, A.E. Zanne, J. Chave, S.J. Wright, S.N. Sheremet’ev, H. Jactel, C. Baraloto, B. Cerabolini, S. Pierce, B. Shipley, D. Kirkup, F. Casanoves, J.S. Joswig, A. Günther, V. Falczuk, N. Rüger, M.D. Mahecha, and L.D. Gorné. 2016. The global spectrum of plant form and function. Nature 529: 167–171. ESRI 2023. ArcGIS Pro. Release: 3.1.2. Redlands, CA: Environmental Systems Research Institute. Field, C.B., and H.A. Mooney. 1986. The photosynthesis-nitrogen relationship in wild plants. Pages 25–56 in T.J. Givnish, editor. On the Economy of Plant Form and Function. Cambridge University Press, Cambridge. Fowler, Z.K., M.B. Adams, and W.T. Peterjohn. 2015. Will more nitrogen enhance carbon storage in young forest stands in central Appalachia? Forest Ecology & Management 337: 144-152. Gara, T.W., R. Darvishzadeh, A.K. Skidmore, T. Wang, and M. Heurich. 2019. Accurate modelling of canopy traits from seasonal Sentinel-2 imagery based on the vertical distribution of leaf traits. ISPRS Journal of Photogrammetry and Remote Sensing 157: 108-123. Gelman, A. 2007. Scaling regression inputs by dividing by two standard deviations. Statistics in Medicine 27: 2865–2873. Gough, C.M., J.W. Atkins, R.T. Fahey, & B.S. Hardiman. 2019. High rates of primary productivity in structurally complex forests. Ecology 100(10): e02864. Hardiman, B.S., G. Bohrer, C.M. Gough, C.S. Vogel, P.S. Curtis, S. Vogel, S. Curtis, and B.S. Hardiman. 2011. The role of canopy structural complexity in wood net primary production of a maturing northern deciduous forest. Ecology 92: 1818–1827. 31 Hardiman, B.S., C.M. Gough, A. Halperin, K.L. Hofmeister, L.E. Nave, G. Bohrer, and P.S. Curtis. 2013. Maintaining high rates of carbon storage in old forests: a mechanism linking canopy structure to forest function. Forest Ecology and Management 298: 111– 119. Harding, D. J., M. A. Lefsky, G. G. Parker, and J. B. Blair. 2001. Laser altimeter canopy profiles methods and validation for closed-canopy, broadleaf forests. Remote Sensing of Environment 76: 283–297. Hijmans, R.J. 2022. Terra: Spatial Data Analysis. R package version 1.7-46. https://CRAN.R- project.org/package=terra Kamoske, A.G., K. M. Dahlin, S. C. Stark, and S. P. Serbin. 2019. Leaf area density from airborne LiDAR: Comparing sensors and resolutions in a temperate broadleaf forest ecosystem. Forest Ecology and Management 433: 364–375. Kamoske, A.G., K.M. Dahlin, S.P. Serbin and S.C. Stark. 2021. Leaf traits and canopy structure together explain canopy functional diversity: an airborne remote sensing approach. Ecological Applications 31: e02230. Kamoske, A. G., K. M. Dahlin, Q. Read, S. Record, S. C. Stark, S. P. Serbin, and P. L. Zarnetske. 2022. Towards mapping biodiversity from above: Can fusing lidar and hyperspectral remote sensing predict taxonomic, functional, and phylogenetic tree diversity in temperate forests? Global Ecology and Biogeography 31(7): 1440-1460. Kampe, T. U., B.R. Johnson, M. Kuester, and M. Keller. 2010. NEON: the first continental-scale ecological observatory with airborne remote sensing of vegetation canopy biochemistry and structure. Journal of Applied Remote Sensing 4: e043510. Kane, V.R., J.D. Bakker, R.J. McGaughey, J.A. Lutz, R.F. Gersonde and J.F. Franklin. 2010. Examining conifer canopy structural complexity across forest ages and elevations with LiDAR data. Canadian Journal of Forest Research 40(4): 774-787. Knyazikhin, Y., M. A. Schull, P. Stenberg, M. Mõttus, M. Rautiainen, Y. Yang, A. Marshak, P. L. Carmona, R. K. Kaufmann, P. Lewis, M. I. Disney, V. Vanderbilt, A. B. Davis, F. Baret, S. Jacquemoud, A. Lyapustin, and R. B. Myneni. 2013. Hyperspectral remote sensing of foliar nitrogen. PNAS 110(3): 185-192. LaRue, E.A., B.S. Hardiman, J.M. Elliot, and Songlin Fei. 2019. Structural diversity as a predictor of ecosystem function. Environmental Research Letters 14(11): 114011. Lawrence, D. M., R. A. Fisher, C. D. Koven, K. W. Oleson, S. C. Swenson, G. Bonan, N. Collier, B. Ghimire, L. van Kampenhout, D. Kennedy, E. Kluzek, P. J. Lawrence, F. Li, H. Li, D. Lombardozzi, W. J. Riley, W. J. Sacks, M. Shi, M. Vertenstein, W. R. Wieder, C. Xu, A.A. Ali, A.M. Badger, G. Bisht, M. van den Broeke, M.A. Brunke, S.P. Burns, J. Buzan, M. Clark, A. Craig, K. Dahlin, B., Drewniak, J.B., Fisher, M., Flanner, A.M., Fox, P. Gentine, F. Hoffman, G. Keppel-Aleks, R. Knox, S. Kumar, J. Lenaerts, L.R. Leung, W.H. Lipscomb, Y. Lu, A. Pandey, J.D. Pelletier, J. Perket, J.T. Randerson, D.M. 32 Ricciuto, B.M. Sanderson, A. Slater, Z.M. Subin, J. Tang, R.Q. Thomas, M. Val Martin, and X. Zeng. 2019. The Community Land Model Version 5: Description of new features, benchmarking, and impact of forcing uncertainty. Journal of Advances in Modeling Earth Systems 11, 4245–4287. Lefsky, M.A., W.B. Cohen, S.A. Acker, G.G. Parker, T.A. Spies, and D. Harding. 1999. Lidar remote sensing of the canopy structure and biophysical properties of Douglas-fir Western Hemlock forests. Remote Sensing of Environment 70:339–361. Liu, J., S. Liu, and T.R. Loveland. 2006. Temporal evolution of carbon budgets of the Appalachian forests in the U.S. from 1972 to 2000. Forest Ecology & Management 222(1-3): 191-201. Macfarlane, C., A. Grigg, and C. Evangelista. 2007. Estimating forest leaf area using cover and fullframe fisheye photography: Thinking inside the circle. Agricultural and Forest Meteorology 146(1-2): 1–12. Mollot, G., J.H. Pantel, and T.N. Romanuk. 2017. Chapter two – the effects of invasive species on the decline in species richness: a global meta-analysis. Advances in Ecological Research 56: 61- 83. NEON. 2021b. Site management and event reporting (DP1.10111.001). National Ecological Observatory Network (NEON). Niinemets, U. 2007. Photosynthesis and resource distribution through plant canopies. Plant, Cell & Environment 30(9): 1052-1071. Niinemets, U. 2012. Optimization of foliage photosynthetic capacity in tree canopies: towards identifying missing constraints. Tree Physiology 32: 505–509. Niinemets, U., T. F. Keenan, and L. Hallik. 2015. A worldwide analysis of within-canopy variations in leaf structural, chemical and physiological traits across plant functional types. New Phytologist 205: 973–993. Pan, Y., R.A. Birdsey, J. Fang, R. Houghton, P.E. Kauppi, W.A. Kurz, O.L. Phillips, A. Shviden ko, S.L. Lewis, J.G. Canadell, P. Ciais, R.B. Jackson, S.W. Pacala, A.D. McGuire, S. Pia o, A. Rautiainen, S. Sitch, and D. Hayes. 2011. A large and persistent carbon sink in the world’s forests. Science 333: 988-993. Parker, G.G., M.A. Lefsky, and D.J. Harding. 2001. Light transmittance in forest canopies determined using airborne laser altimetry and in-canopy quantum measurements. Remote Sensing of Environment 76(3): 298-309. Planet Team. 2020. “Planet Surface Reflectance Product v2.” Planet Labs, Inc, Accessed 18.08.2020. https://assets.planet.com/marketing/PDF/Planet_Surface_Reflectance_ Technical_White_Paper.pdf. Prior, E., V. Thomas, and R. Wynne. 2022. Estimation of mean dominant height using NAIP digital aerial photogrammetry and lidar over mixed deciduous forest in the southeastern USA. International Journal of Applied Earth Observation and Geoinformation 110 (e102813). 33 Redmond, M.D., R.B. Wilbur, and H.M. Wilbur. 2012. Recruitment and Dominance of Quercus rubra and Quercus alba in a previous Oak-Chestnut Forest from the 1980s to 2008. The American Midland Naturalist 168(2): 427-442. R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/. Sabol, J., Z. Patocka, and T. Mikita. 2014. Usage of Lidar data for leaf area index estimation. GeoScience Engineering 60: 10–18. Schimel, D.S. 1995. Terrestrial biogeochemical cycles: global estimates with remote sensing. Remote Sensing of Environment 51: 49-56. Serbin, S.P., A. Singh, B.E. McNeil, C.C. Kingdon, and P. Townsend. 2014. Spectroscopic determination of leaf morphological and biochemical traits for northern temperate and boreal tree species. Ecological Applications 24: 1651-1669. Shugart, H.H., S. Saatchi, and F.G. Hall. 2010. Importance of structure and its measurement in quantifying function of forest ecosystems. Journal of Geophysical Research, 115: G00E13. Simberloff, D., and B. von Holle. 1999. Positive interaction of nonindigenous species: invasional meltdown? Biological Invasions 1: 21-32. Small, C.J. and J.L. Chamberlain. 2015. Forest diversity and disturbance: changing influences and the future of Virginia’s forests. Virginia Journal of Science 66(3): Article 7. Smith, M., S.V. Ollinger, M. E. Martin, J.D. Aber, R. A. Hallett, and C.L. Goodale. 2002. Direct estimation of aboveground forest productivity through hyperspectral remote sensing of canopy nitrogen. Ecological Applications 12(5): 1286-1302. Stark, S.C., V. Leitold, J.L. Wu, M.O. Hunter, C.V. de Castilho, F.R.C. Costa, S.M. McMahon, G.G. Parker, Mônica, T. Shimabukuro, M.A. Lefsky, M. Keller, L.F. Alves, J. Schietti, Y.E. Shimabukuro, D.O. Brandão, T.K. Woodcock, N. Higuchi, P.B. de Camargo, R.C. de Oliveira, Scott R. Saleska. 2012. Amazon forest carbon dynamics predicted by profiles of canopy leaf area and light environment. Ecology Letters 15(12): 1406-1414. Stark, S.C., B.J. Enquist, S.R. Saleska, V. Leitold, J. Schietti, M. Longo, L.F. Alves, P.B. Camargo, and R.C. Oliveira. 2015. Linking canopy leaf area and light environments with tree size distributions to explain Amazon forest demography. Ecology Letters 18: 636– 645. Thomas, R.Q., C.D. Canham, K.C. Weathers, and C.L. Goodale. 2010. Increased tree carbon storage in response to nitrogen deposition in the US. Nature Geoscience 3: 13–17. Townsend, P. A., J.R. Foster, R.A. Chastian Jr., and W.S. Currie. 2003. Canopy nitrogen in the forests of the Central Appalachian Mountains using Hyperion and AVIRIS. IEEE Transactions on Geoscience and Remote Sensing 41: 1347–1354. 34 Tremblay, J. P., S.D. Côté, T.P. Rooney, and C. Dussault, and D. M. Waller. 2004. Ecological impacts of deer over-abundance. Annual Review of Ecology, Evolution, and Systematics 35:113–147. Violle, C., M-L. Navas, D. Vile, E. Kazakou, C. Fortunel, I. Hummel, and E. Garnier. 2007. Let the concept of trait be functional! OIKOS 116(5): 882-892. Wang, Z., A. Chlus, R. Geygan, Z. Ye, T. Zheng, A. Singh, J.J. Couture, J. Cavender-Bares, E.L. Kruger, and P.A. Townsend. 2020. Foliar functional traits from imaging spectroscopy across biomes in eastern North America. New Phytologist 228: 494-511. Wilbur, R.B. H. M. Wilbur, and N.J. Peterson. 2021. Flora of the Forest at Mountain Lake Biological Station, Giles County, Virginia. Castanea 85(2): 376-403. Wilfahrt, P.A., B. Collins, and P.S. White. 2014. Shifts in functional traits among tree communities across succession in eastern deciduous forests. Forest Ecology and Management 324: 179-185. Zhang, Y., X. Xiao, C. Jin, J. Dong, S. Zhou, P. Wagle, J. Joiner, L. Guanter, Y. Zhang, G. Zhang, Y. Qin, J. Wang, and B. Moore III. 2016. Consistency between sun-induced chlorophyll fluorescence and gross primary production of vegetation in North America. Remote Sensing of Environment 183: 154-169. Zhang, Y., X. Xiao, X. Wu, S. Zhou, G. Zhang, Y. Qin, and J. Dong. 2017. A global moderate resolution dataset of gross primary production of vegetation for 2000–2016. Sci Data 4: 170165. 35 APPENDIX A Figure A.1. Scatterplot of Lidar-derived LAI and hemispherical photography derived LAI pre calibration with Beer-Lambert extinction coefficient (k = NULL). 36 Figure A.2. Scatterplot of Lidar-derived LAI and hemispherical photography derived LAI post calibration with Beer-Lambert extinction coefficient (k = .53). 37 Table A.1. lidar derived metrics Metric Source Empty voxel ratio Filled voxel ratio Canopy rugosity Top-of-canopy rugosity Euphotic zone volume Euphotic zone total leaf area Oligophotic zone volume Oligophotic zone total leaf area NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language Canopy oligophotic zone total leaf area (3x3 window) NEON AOP lidar; R programming Canopy oligophotic zone volume (3x3 window) Canopy empty volume (3x3 window) Euphotic depth Canopy Euphotic zone volume Canopy Euphotic zone total leaf area Max lad height Porosity ratio LAD height quantiles (x4) Canopy height model Leaf area index (LAI) LAD column standard deviation language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language 38 Table A.1. (cont’d) Digital terrain model (DTM) Digital surface model (DSM) Table A.2. Environmental drivers NEON AOP lidar; R programming language NEON AOP lidar; R programming language Metric Source Aspect DTM (Elevation) Easting Northing Flow accumulation Slope Topographic Position Index (TPI) Topographic Wetness Index (TWI) Topographic Roughness index (TRI) Solar radiation at summer solstice Solar radiation at winter solstice Difference in NDVI (2018-2016) Plant functional type NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language NEON AOP lidar; R programming language PlanetScope Surface Reflectance V2; R programming language PlanetScope Surface Reflectance V2; R programming language 39