https://orcid.org/0000-0001-9666-5805
Abstract
Pinus edulis Engelm. is a short-stature, drought-tolerant tree species that is abundant in piñon-juniper woodlands throughout semiarid ecosystems of the American Southwest. P. edulis is a model species among ecophysiological disciplines, with considerable research focus given to hydraulic functioning and carbon partitioning relating to mechanisms of tree mortality. Many ecological studies require robust estimates of tree structural traits such as biomass, active sapwood area, and leaf area. We harvested twenty trees from Central New Mexico ranging in size from 1.3 to 22.7 cm root crown diameter (RCD) to derive allometric relationships from measurements of RCD, maximum height, canopy area (CA), aboveground biomass (AGB), sapwood area (AS), and leaf area (AL). Total foliar mass was measured from a subset of individuals and scaled to AL from estimates of leaf mass per area. We report a strong nonlinear relationship to AGB as a function of both RCD and height, whereas CA scaled linearly. Total AS expressed a power relationship with RCD. Both AS and CA exhibited strong linear relationships with AL (R2 = 0.99), whereas RCD increased nonlinearly with AL. We improve on current models by expanding the size range of sampled trees and supplement the existing literature for this species.
Study Implications: Land managers need to better understand carbon and water dynamics in changing ecosystems to understand how those ecosystems can be sustainably used now and in the future. This study of two-needle pinon (Pinus edulis Engelm.) trees in New Mexico, USA, uses observations from unoccupied aerial vehicles, field measurements, and harvesting followed by laboratory analysis to develop allometric models for this widespread species. These models can be used to understand plant traits such biomass partitioning and sap flow, which in turn will help scientists and land managers better understand the ecosystem services provided by pinon pine across North America.
Two-needle piñon pine (Pinus edulis Engelm.) is among the most abundant and widely distributed tree species throughout the southwestern US (Anderson 2002), with a range spanning Utah, Colorado, Arizona, and New Mexico and minor populations in California, Wyoming, Oklahoma, and Texas. Several Juniperus species (J. monosperma Engelm., J. osteosperma [Torr.] Little, J. erythrocarpa [Martinez] Gaussen ex R.P. Adams, J. deppeana Steud., and J. scopulorum Sarg.) are commonly associated with P. edulis and three closely related piñon pines (P. cembroides Zucc., P. discolor D.K. Bailey & Hawksw., and P. monophylla Torr. & Frem) in mixed woodlands, with the piñon-juniper cover type covering 19–22 million hectares throughout the intermountain West (Shaw 2006, Miller et al. 2008, 2019). P. edulis provides essential ecosystem services in this region including regulatory services such as carbon sequestration (Huang et al. 2009), soil water and temperature regulation and partitioning of evapotranspiration (Morillas et al. 2017), and provisioning services, such as furnishing wildlife habitat (Finch and Ruggiero 1993, Bender et al. 2007, Bombaci and Pejchar 2016). The pine nut industry is associated with important cultural and tribal activities throughout the southwestern US (Tarancón et al. 2020). Piñon-juniper woodlands have a heterogeneous canopy cover throughout the landscape, with high bare ground connectivity that drive gradients of soil water potential and surface evaporation (Newman et al. 2010, Puttock et al. 2013). The physiology of P. edulis is well studied and it is regarded as a model species for plant hydraulic research because it exercises tight stomatal control in response to moisture deficits (McDowell et al. 2008, Plaut et al. 2012, Sevanto et al. 2014, Hudson et al. 2018).
Despite its status as among the most drought-tolerant conifer species, P. edulis was abundant during a time (ca. 23,400 years before present) when regional precipitation was considerably higher than at present (Cole et al. 2013). Extensive woodland expansion throughout the last century is attributed to anthropogenic practices such as fire suppression, livestock grazing of competing vegetation, and elevated levels of carbon dioxide (Romme et al. 2009). The expansion of the piñon-juniper woodland since European settlement in the 19th century was primarily into sagebrush ecosystems (Reinhardt et al. 2020). However, climate change within the Anthropocene has increased the frequency and severity of droughts throughout the range of P. edulis with the potential to alter the current species distribution (Garfin et al. 2013, Hatten et al. 2016). Drought is a common inciting factor of diseases and can compromise a tree’s ability to synthesize and induce chemical defenses (e.g., resin and terpene compounds) due to a lack of available carbohydrate, especially over prolonged dry periods. Consequently, P. edulis throughout the southwestern US has become susceptible to invasion by native bark beetles (Ips confusus Leconte), resulting in rapid and often widespread mortality events (Meddens et al. 2015). Based on the expected continuation of increased temperature and shifts in precipitation in the southwestern US in the coming decades, the range of P. edulis is projected to contract due to the lethal combination of hydraulic and biotic stressors (Greenwood and Weisberg 2008, Adams et al. 2009, Breshears et al. 2018, Flake and Weisberg 2019).
Due to the ever-changing nature of the piñon-juniper cover type across the landscape, the surface fuel load it represents (Wozniak et al. 2020), and important ecosystem services provided, it is necessary to monitor physical traits in an efficient manner that can be scaled easily to the landscape level. This requires a diverse toolbox of allometric models. With the burgeoning adoption of aerial photography techniques using unoccupied aerial vehicles (UAV or drones), it is now possible to survey areas up to 100 ha rapidly and efficiently (Anderson and Gaston 2013, Fu et al. 2021). Specifically, the ability to use top-down photography to estimate tree height and canopy cover has become common practice for both research and management objectives (She et al. 2014, Cunliffe et al. 2016, Jiménez-Brenes et al. 2017, Panagiotidis et al. 2017). Previous work by Darling (1967), Grier et al. (1992), and Sprinkle and Klepac (2015) derived diameter to aboveground biomass relationships for P. edulis throughout Arizona and Utah, and later work by Jenkins et al. (2003) and Chojnacky et al. (2014) incorporates these models within generalized allometric equations for softwood and woodland pine species, respectively, intended for continental-scale application. Furthermore, there are only two reports estimating sapwood area from diameter for this species (West et al. 2008, Pangle et al. 2015), and a single report for estimating total leaf area from diameter through scaling of branch segments (Pangle et al. 2015).
Here, we investigated several allometric relationships of P. edulis that can be used for predicting meaningful physical attributes for this important tree species. We focus on deriving relationships between stem diameter, maximum height, aboveground biomass (AGB), canopy area (CA), total leaf area (AL), and cross-sectional sapwood area (AS). To accomplish this, we conducted aerial surveys coupled with traditional ground measurement of tree attributes, followed by harvesting and dissection of selected individual trees.
Methods
Site Description
We conducted this study in a ~6 ha area of piñon-juniper woodland, dominated by P. edulis (two-needle piñon) and J. monosperma (oneseed juniper), with some C4 grasses (esp. Bouteloua gracilis) and cacti (esp. Cylindropuntia imbricata). The site is located within the Deer Canyon Preserve, near Mountainair in central New Mexico, USA (34.439°N, −106.234°E) at an elevation of 2,150 m above sea level. The mean annual precipitation is 359 mm yr−1 and mean annual temperature is 10.3°C (for 1895–2019; PRISM Climate Group). Soils are classified as Turkey Springs stony loam (NRCS 2018) with the presence of a caliche petrocalcic horizon between the depths of 30 to 80 cm (Morillas et al. 2017). Trees were sampled approximately 0.3 km northeast and downwind of the US-Mpj flux tower, part of the AmeriFlux network (Anderson-Teixeira et al. 2011). Widespread bark beetle (Pinon ips)-driven mortality of piñon occurred at this site from 2013 to 2016, preceded by 2 years of severe drought (Liebrecht 2018).
Tree Selection and Field Measurement
We selected twenty P. edulis trees for our study, sampling evenly across a distribution of sizes from small saplings to the largest individuals observed in the study area. Individuals were selected using a stratified design based on the natural scale of observed size classes (Roxburgh et al. 2015). Our focus was on relatively isolated “open-grown” trees instead of more coalesced clusters whose canopy properties might be more influenced by adjacent trees (Scher et al. 2019). Although piñon pine can often bear multiple primary stems branching near ground level, we selected trees in which canopies originated from a single main stem (bole), consistent with a phenotypic resistance to the shoot-boring moth (Dioryctria albovittella) described by Sthultz et al. (2009). Measured parameters for each tree are reported in Table S1. Following severe droughts in 2011–2012, standing dead piñon were common within the stand, with many live trees exhibiting advanced crown dieback, which was in part exacerbated by the presence of dwarf mistletoe (Arceuthobium divaricatum Hulst) (see Table S2, P-12 and P-15). We visually assessed crown dieback of potential target trees using the average estimate of two trained operators (Schomaker et al. 2007), and selecting individuals with observed crown dieback <30%. Photographs of harvested individuals are included in Table S2. For each tree, we measured the total height and height to the base of the live crown using a telescoping leveling rod prior to harvest. Bole diameter at the root crown (RCD) was measured at ground level to best capture the total basal area of each tree.
Aerial Survey of Canopy Area and Height
We surveyed the selected individuals using a lightweight consumer UAV to acquire aerial images on October 11, 2019. The survey protocol is described in detail in Cunliffe and Anderson (2019). Briefly, we flew surveys using a DJI Phantom 4 Advanced multi-rotor equipped with a built-in CMOS sensor (FC6310 camera) at 8.8 mm focal length capturing 20M effective pixels. We flew two sets of survey flights, the first obtaining nadir imagery at ~26 m above ground level (a.g.l.), and secondly, obtaining oblique (~20° from nadir) images at ca. 31 m a.g.l. Both survey flights obtained 75% forward and side overlap, together capturing at least thirty-four images for each part of the study area. Camera aperture was f5, sensitivity (ISO) was 200, shutter speed was faster than 1/800th of a second, and images were underexposed by 1/3rd of a stop to optimize image quality. To constrain the photogrammetric reconstruction spatially, we deployed thirteen 20 × 20 cm ground control markers across the survey area and geolocated these using an real time kinetic-global navigation satellite system (GNSS).
We processed aerial images using structure-from-motion photogrammetry and a workflow based on Cunliffe et al. (2020). Geotagged image data and marker coordinates were imported into Agisoft PhotoScan (v1.4.3) and converted to a common coordinate reference system (EPSG:: 26913). Image sharpness was assessed using PhotoScan’s image quality tool, and all images had a sharpness score ≥0.84. We matched photos and aligned cameras using the highest quality setting, key point limit of 40,000, tie point limit of 8,000, generic and reference pair preselection enabled, and adaptive camera model fitting disabled. We used the following reference settings: camera location accuracy = XY ± 20 m, Z ± 50 m; marker location accuracy = XY ± 0.02 m, Z ± 0.05 m; marker projection accuracy = two pixels; tie point accuracy = maximum (mean root mean square reprojection error or one). We filtered the sparse cloud, excluded points with reprojection error above 0.45 from further analysis, and estimated camera positions to verify their plausibility. All images were aligned and used for further processing. We placed geolocated markers on ten projected images for each of the thirteen ground control points and deselected three markers used for independent accuracy assessment. We optimized the bundle adjustment using the filtered point cloud with the following lens parameters: focal length (f), principal point (cx, cy), radial distortion (k1, k2), tangential distortion (p1, p2), aspect ratio, and skew coefficient (b1, b2). Multi-view stereopsis (dense point cloud generation) was undertaken using the ultrahigh quality setting, mild depth filtering, and calculation of point colors enabled. We built a triangular irregular network from the dense point cloud via Delaunay triangulation, using a high face count (1/5 the number of points in the dense cloud) and height field approach with interpolation enabled. We built an orthomosaic at the native spatial resolution projected onto the triangular irregular network surface (“model data”) using the mosaic blending mode with fill holes enabled. The orthomosaic was exported as uncompressed GeoTIFF at a spatial grain of 0.01 m. The orthomosaic was loaded into a GIS (ESRI ArcPro v2.1.3), and the edge of each individual’s canopy was manually digitized by the same operator at a scale of 1:40. The area (m2) of each canopy polygon was extracted (CA1), and the length (m) of the widest (a) and perpendicular (b) canopy diameters were extracted using the minimum bounding geometry utility in ArcPro. To test the similarity between canopy area derived directly from polygons (CA1) against canopy area derived from typical field observations (CA2), canopy area was also calculated as follows:
To extract maximum canopy height from the photogrammetric point clouds, we used PDAL (v2.1.0) (PDAL Contributors 2020). We used Delaunay triangulation to interpolate a terrain surface between the GNSS-observed corners of a plot around each individual tree (Cunliffe et al. 2021). The point cloud representing each harvest plot was subset using the corner coordinates. Within each plot, we extracted the maximum canopy height above the terrain model.
Biomass Measurement
We harvested individual trees over 4 days (see Table S1) and separated biomass from each tree into two components: (1) wood tissue (>3 cm stem diameter, including any bole wood) and (2) leaf and twig tissue (<3 cm diameter, hereafter termed “leaf/twig”, Cunliffe et al. 2020). We harvested trees to the ground line and measured the wet weights of both components within 1 hour of each tree being felled. Subsamples of the wood and leaf/twig components for each component of each tree were collected, weighed to determine wet weight, and oven-dried at 80°C for ≥150 hours to a constant weight (defined as <0.1% change in mass in 24 hours) to determine water content. These dried subsamples were up to 14 kg, accounting for at least 3% and up to 100% of each tree. We calculated water contents on a green weight basis for each tree component; that is, (wet mass−dry mass)/wet mass, to convert wet mass to dry mass. We calculated component proportions on a dry-weight basis and AGB as the sum of wood tissue and leaf/twig for each individual tree.
Sapwood Area and Bark Thickness
We defined the RCD as the lowest primary branching point for a single stem, where a stem either intersected with a main bole or with the ground (Grier et al. 1992, Chojnacky and Rogers 1999, Krofcheck et al. 2016). All cross sections (hereafter “disks”) were measured in the field using a thin-line diameter tape for samples >5 cm, or using the mean of two measurements about the major and minor axis with a digital caliper for samples < 5cm to determine the fresh RCD (RCDWet). We obtained a disk from all trees at RCD. After air-drying the disks over ≥30 days to minimize cracking, we remeasured them to determine the dried RCD (RCDDry). We polished dried disks using a combination of belt and orbital sanders at progressively finer grits to clearly distinguish the sapwood-heartwood interface. In some instances, applying a small amount of water evenly across the disk surface with a paper towel prior to imaging provided additional clarity for distinguishing the heartwood-sapwood interface. We captured fine resolution images (45.7 Megapixels) for each disk using a DSLR camera (Nikon D850) and measured the total sapwood area (AS) of each disk using the image processing software ImageJ (Schneider et al. 2012). It should be noted that total sapwood area can differ from the conductive xylem area, which is typically measured via dye perfusion or tomography, and was not assessed within this study.
Bark thickness was also estimated from imaged disks as:
where BA is the outside basal area and IBA is the inside-bark basal area. This method for computing bark thickness provides an integrated average of bark thickness around the disk, which can be variable (Stängle and Dormann 2018), as opposed to conducting multiple measurements about the outer radii of each sample.
Needle Mass and Foliar Area
Total needle dry mass was measured for five trees spanning in RCD from 1.2 to 19.3 cm. All needles were dried attached to stem segments using a combination of oven and air drying. Once dry, needles were separated from stems by hand and weighed with an electronic balance to 0.1 g. A subset of fresh needles was imaged on a flatbed scanner and the projected area measured using ImageJ software (Glozer 2008). The leaf mass per area (LMA, g m−2) was calculated by dividing the total dry weight for a combination of the needle age classes represented from the subsamples by the total projected scanned area. Total AL for the five target trees was calculated by scaling the total weight of dry foliage using the mean LMA.
Statistical Analysis
Statistical analyses were conducted in R (v4.0.2, R Core Team 2020). To evaluate relationships between maximum height measured on the ground versus that extracted from the photogrammetry and for comparing canopy area estimates of CA1 and CA2, we applied total least square regression to best account for uncertainty in both dependent and independent variables. For allometry of tree height as a function of RCD, bark thickness as a function of fresh disk diameter, and predictions of biomass as a function of canopy areas, we fitted linear models with ordinary least squares regression. Power models were applied for predictions of basal sapwood area as a function of fresh disk diameter, biomass from maximum tree height, and biomass as a function of RCD.
To evaluate total least square linear regression fits, we used Lin’s concordance correlation coefficient (CCC) (Lin 1989) using the DescTools package (Signorell et al. 2021), and adjusted coefficient of determination (R2) was used for total ordinary least square models. Nonlinear models were evaluated by the root mean square error (RMSE) using the Metrics package (Hamner et al. 2018). To test for potential influence of canopy dieback on AGB estimates, we performed ANOVA on model residuals using three binned categories of observed dieback (0%–10%, 10%–20%, and 20%–30%).
Results
Tree Height and Canopy Areas
Key tree metrics are provided in Table S1. For the twenty trees harvested in this study, RCD ranged from 1.3 to 22.7 cm (mean = 12.0 cm, SD = 7.5 cm) with maximum tree height ranging from 0.3 to 6.7 m (mean = 3.3 m, SD = 2.2 m). Live crown ratios ranged from 63% to 100% (mean = 86%, SD = 14%). Even trees with observed crown dieback of 20%–30% exhibited live crown ratios >75% on average. We measured only one tree with a dieback rating of 30% but a live crown that extended to ground level, and several dead branches within the midcanopy (P-15, Table S2). We found that RCD scaled linearly with total tree height measured in the field (R2 = 0.93, RMSE = 0.56, P < 0.001). Field measurements of tree height were very strongly correlated with maximum height derived from UAV photogrammetry (CCC = 0.979; Figure 1A). Estimates of projected canopy area derived from digitized polygons (CA1) and from orthogonal dimensions of the canopy polygons (CA2) were highly correlated with a linear relationship (CCC = 0.996; Figure 1B).
(A) Maximum canopy height estimated from UAV as a function of measured height using a telescoping leveling rod. (B) Comparison of UAV derived canopy area using polygons (CA1) and canopy area estimated from perpendicular crown diameters (CA2). Fitted lines (solid) show the total least squares regression, dashed lines show 1:1.
Aboveground Biomass
A total of fifty-five subsamples totaling 104.9 kg of fresh biomass were used to determine tree water content and estimate dry biomass. Water content for biomass subsamples in the <3 cm size class ranged from 36% to 63% (mean = 46%, SD = 7%) and in the larger >3 cm size class from 34% to 48% (mean = 39%, SD = 4%). The higher water content within the <3 cm size class was expected due to the inclusion of foliage. The fraction of dry biomass <3 cm had a significant negative linear relationship to RCD (R2 = 0.83, RMSE = 0.098, P < 0.001). As trees grow larger, the fine branches and foliage account for a smaller proportion of total mass. Total dry biomass ranged from 0.03 to 171.01 kg (mean = 43.9, SD = 53.2) across all trees. RCD was a strong predictor of total tree dry biomass using a power model, with RMSE = 18.18 (Figure 2A, Table 1). The largest trees in our study (RCD > 20 cm, n = 5) accounted for most of the variance within our model. These large trees accounted for a mean absolute deviation of 21.60 kg, whereas all other trees (n = 15) accounted for just 8.74 kg. For the largest trees, our model both over- and underpredicted total biomass, although we found no effect of crown dieback on model residuals that may have influenced the total standing live biomass (p = 0.233). Tree height was also a strong predictor of total dry biomass (Figure 2B, Table 1) using a power model, although we noted one outlier among sampled trees that considerably influenced the fit of this model. This outlying individual displayed a wide open-grown crown that was atypical of other P. edulis within the stand (tree P-14, Table S2). By removing this outlier, the RMSE of this model was reduced by 66%, from RMSE = 28.11 kg to RMSE = 9.47 kg, providing an overall better predictive model than the biomass relationship to RCD. Finally, metrics of canopy area provided reliable predictions of AGB, although neither canopy estimate (CA1 or CA2) performed statistically better than the other (R2 = 0.924; Figure 2C, Table 1).
Parameters and equations for statistical models.
| Dependent variable | Predictor variable | Model form | a (SE) | b (SE) | P value | Adj R2 | RMSE | n |
|---|---|---|---|---|---|---|---|---|
| AGB (kg) | RCD (cm) | Y = a × Xb | 0.018 (0.031) | 2.851 (0.552) | — | — | 18.119 | 20 |
| AGB (kg) | Height (m) | Y = a × X b | 0.971 (0.478) | 2.563 (0.272) | — | — | 9.456 | 20 |
| AGB (kg) | CA1 (m2) | Y = a + b × X | −6.099 (4.635) | 10.144 (0.623) | <0.001 | 0.924 | 13.901 | 20 |
| AGB (kg) | CA2 (m2) | Y = a + b × X | −4.685 (4.573) | 9.643 (0.633) | <0.001 | 0.924 | 13.911 | 20 |
| AS (cm2) | RCD (cm) | Y = a × X b | 0.878 (0.561) | 1.695 (0.216) | — | — | 17.969 | 20 |
| AS (cm2) | CA1 (m2) | Y = a + b × X | 12.595 (9.631) | 11.196 (1.382) | <0.001 | 0.773 | 28.886 | 20 |
| AS (cm2) | CA2 (m2) | Y = a + b × X | 14.506 (9.720) | 10.574 (1.345) | <0.001 | 0.762 | 29.568 | 20 |
| AL (m2) | AS (cm2) | Y = a + b × X | −0.917 (1.054) | 0.413 (0.017) | <0.001 | 0.994 | 1.444 | 5 |
| AL (m2) | RCD (cm) | Y = a × X b | 0.004 (0.002) | 3.216 (0.203) | — | — | 0.969 | 5 |
| AL (m2) | CA1 (m2) | Y = a + b × X | −1.877 (1.168) | 7.585 (0.331) | <0.001 | 0.992 | 1.564 | 5 |
| AL (m2) | CA2 (m2) | Y = a + b × X | −1.903 (1.221) | 7.741 (0.352) | <0.001 | 0.992 | 1.634 | 5 |
| Height (m) | RCD (cm) | Y = a + b × X | −0.042 (0.250) | 0.276 (0.018) | <0.001 | 0.926 | 0.556 | 20 |
| Bark thickness (cm) | RCD (cm) | Y = a + b × X | 0.129 (0.064) | 0.055 (0.005) | <0.001 | 0.879 | 0.141 | 20 |
| Dependent variable | Predictor variable | Model form | a (SE) | b (SE) | P value | Adj R2 | RMSE | n |
|---|---|---|---|---|---|---|---|---|
| AGB (kg) | RCD (cm) | Y = a × Xb | 0.018 (0.031) | 2.851 (0.552) | — | — | 18.119 | 20 |
| AGB (kg) | Height (m) | Y = a × X b | 0.971 (0.478) | 2.563 (0.272) | — | — | 9.456 | 20 |
| AGB (kg) | CA1 (m2) | Y = a + b × X | −6.099 (4.635) | 10.144 (0.623) | <0.001 | 0.924 | 13.901 | 20 |
| AGB (kg) | CA2 (m2) | Y = a + b × X | −4.685 (4.573) | 9.643 (0.633) | <0.001 | 0.924 | 13.911 | 20 |
| AS (cm2) | RCD (cm) | Y = a × X b | 0.878 (0.561) | 1.695 (0.216) | — | — | 17.969 | 20 |
| AS (cm2) | CA1 (m2) | Y = a + b × X | 12.595 (9.631) | 11.196 (1.382) | <0.001 | 0.773 | 28.886 | 20 |
| AS (cm2) | CA2 (m2) | Y = a + b × X | 14.506 (9.720) | 10.574 (1.345) | <0.001 | 0.762 | 29.568 | 20 |
| AL (m2) | AS (cm2) | Y = a + b × X | −0.917 (1.054) | 0.413 (0.017) | <0.001 | 0.994 | 1.444 | 5 |
| AL (m2) | RCD (cm) | Y = a × X b | 0.004 (0.002) | 3.216 (0.203) | — | — | 0.969 | 5 |
| AL (m2) | CA1 (m2) | Y = a + b × X | −1.877 (1.168) | 7.585 (0.331) | <0.001 | 0.992 | 1.564 | 5 |
| AL (m2) | CA2 (m2) | Y = a + b × X | −1.903 (1.221) | 7.741 (0.352) | <0.001 | 0.992 | 1.634 | 5 |
| Height (m) | RCD (cm) | Y = a + b × X | −0.042 (0.250) | 0.276 (0.018) | <0.001 | 0.926 | 0.556 | 20 |
| Bark thickness (cm) | RCD (cm) | Y = a + b × X | 0.129 (0.064) | 0.055 (0.005) | <0.001 | 0.879 | 0.141 | 20 |
AGB, aboveground biomass; RCD, root crown diameter; CA1, canopy area extracted from polygons; CA2, canopy area calculated from two perpendicular measurements of canopy diameter; AS, sapwood area; SE, standard error; Adj R2, adjusted coefficient of determination; RMSE, root mean square error; n, number of individuals used to fit the model.
Parameters and equations for statistical models.
| Dependent variable | Predictor variable | Model form | a (SE) | b (SE) | P value | Adj R2 | RMSE | n |
|---|---|---|---|---|---|---|---|---|
| AGB (kg) | RCD (cm) | Y = a × Xb | 0.018 (0.031) | 2.851 (0.552) | — | — | 18.119 | 20 |
| AGB (kg) | Height (m) | Y = a × X b | 0.971 (0.478) | 2.563 (0.272) | — | — | 9.456 | 20 |
| AGB (kg) | CA1 (m2) | Y = a + b × X | −6.099 (4.635) | 10.144 (0.623) | <0.001 | 0.924 | 13.901 | 20 |
| AGB (kg) | CA2 (m2) | Y = a + b × X | −4.685 (4.573) | 9.643 (0.633) | <0.001 | 0.924 | 13.911 | 20 |
| AS (cm2) | RCD (cm) | Y = a × X b | 0.878 (0.561) | 1.695 (0.216) | — | — | 17.969 | 20 |
| AS (cm2) | CA1 (m2) | Y = a + b × X | 12.595 (9.631) | 11.196 (1.382) | <0.001 | 0.773 | 28.886 | 20 |
| AS (cm2) | CA2 (m2) | Y = a + b × X | 14.506 (9.720) | 10.574 (1.345) | <0.001 | 0.762 | 29.568 | 20 |
| AL (m2) | AS (cm2) | Y = a + b × X | −0.917 (1.054) | 0.413 (0.017) | <0.001 | 0.994 | 1.444 | 5 |
| AL (m2) | RCD (cm) | Y = a × X b | 0.004 (0.002) | 3.216 (0.203) | — | — | 0.969 | 5 |
| AL (m2) | CA1 (m2) | Y = a + b × X | −1.877 (1.168) | 7.585 (0.331) | <0.001 | 0.992 | 1.564 | 5 |
| AL (m2) | CA2 (m2) | Y = a + b × X | −1.903 (1.221) | 7.741 (0.352) | <0.001 | 0.992 | 1.634 | 5 |
| Height (m) | RCD (cm) | Y = a + b × X | −0.042 (0.250) | 0.276 (0.018) | <0.001 | 0.926 | 0.556 | 20 |
| Bark thickness (cm) | RCD (cm) | Y = a + b × X | 0.129 (0.064) | 0.055 (0.005) | <0.001 | 0.879 | 0.141 | 20 |
| Dependent variable | Predictor variable | Model form | a (SE) | b (SE) | P value | Adj R2 | RMSE | n |
|---|---|---|---|---|---|---|---|---|
| AGB (kg) | RCD (cm) | Y = a × Xb | 0.018 (0.031) | 2.851 (0.552) | — | — | 18.119 | 20 |
| AGB (kg) | Height (m) | Y = a × X b | 0.971 (0.478) | 2.563 (0.272) | — | — | 9.456 | 20 |
| AGB (kg) | CA1 (m2) | Y = a + b × X | −6.099 (4.635) | 10.144 (0.623) | <0.001 | 0.924 | 13.901 | 20 |
| AGB (kg) | CA2 (m2) | Y = a + b × X | −4.685 (4.573) | 9.643 (0.633) | <0.001 | 0.924 | 13.911 | 20 |
| AS (cm2) | RCD (cm) | Y = a × X b | 0.878 (0.561) | 1.695 (0.216) | — | — | 17.969 | 20 |
| AS (cm2) | CA1 (m2) | Y = a + b × X | 12.595 (9.631) | 11.196 (1.382) | <0.001 | 0.773 | 28.886 | 20 |
| AS (cm2) | CA2 (m2) | Y = a + b × X | 14.506 (9.720) | 10.574 (1.345) | <0.001 | 0.762 | 29.568 | 20 |
| AL (m2) | AS (cm2) | Y = a + b × X | −0.917 (1.054) | 0.413 (0.017) | <0.001 | 0.994 | 1.444 | 5 |
| AL (m2) | RCD (cm) | Y = a × X b | 0.004 (0.002) | 3.216 (0.203) | — | — | 0.969 | 5 |
| AL (m2) | CA1 (m2) | Y = a + b × X | −1.877 (1.168) | 7.585 (0.331) | <0.001 | 0.992 | 1.564 | 5 |
| AL (m2) | CA2 (m2) | Y = a + b × X | −1.903 (1.221) | 7.741 (0.352) | <0.001 | 0.992 | 1.634 | 5 |
| Height (m) | RCD (cm) | Y = a + b × X | −0.042 (0.250) | 0.276 (0.018) | <0.001 | 0.926 | 0.556 | 20 |
| Bark thickness (cm) | RCD (cm) | Y = a + b × X | 0.129 (0.064) | 0.055 (0.005) | <0.001 | 0.879 | 0.141 | 20 |
AGB, aboveground biomass; RCD, root crown diameter; CA1, canopy area extracted from polygons; CA2, canopy area calculated from two perpendicular measurements of canopy diameter; AS, sapwood area; SE, standard error; Adj R2, adjusted coefficient of determination; RMSE, root mean square error; n, number of individuals used to fit the model.
(A) Total dry biomass as a function of RCD. Power relationships from previously published works are shown in colored dashed lines. (B) Total dry biomass as a function of measured tree height. The red dashed line is the model fitted to all observations while the solid black line is the model fitted excluding one outlier (ID = P14, red open circle) which exhibited atypical growth habit, being of shorter stature with an open-grown spreading crown. (C) Total dry biomass as a function of estimated canopy areas derived from polygons (CA1, solid line and filled circles) and from canopy diameters (CA2, dashed line and open circles). Shaded areas denote the 95% confidence interval of linear regressions.
Sapwood Area and Bark Thickness
For disks sampled for the assessment of sapwood area and bark thickness, the fresh diameter was on average 2.8% lower than RCD measured in the field. Compared with wet disks, the dry disk diameter was reduced by 5.9% across all samples. Both wet and dry disk diameters were strongly correlated with RCD measured in the field (P < 0.001, R2 > 0.99; Fig. S2). Stem diameter was a strong predictor of sapwood area using a power relationship (Figure 3A, Table 1). Bark thickness scaled linearly with wet disk diameter (R2 = 0.88, RMSE = 0.14, P < 0.001), with a maximum back thickness of 1.5 cm among trees sampled in this study (Fig. S3, Table 1).
(A) Total sapwood area as a function of fresh stem diameter. Fitted models from this study, Pangel et al. (2015), and West et al. (2008) are shown with solid black, red dashed, blue dashed lines respectively. (B) Total sapwood area as a function of CA1 (solid line and filled circles) and CA2 (dashed line and open circles). Shaded areas denote the 95% confidence interval of linear regression.
Leaf Area
Mean LMA was 263.0 g m−2 (SD ± 28.9) for foliar subsamples in this study. The total projected AL for our five subsampled trees ranged from 0.1 to 54.0 m2, with total foliar dry mass ranging from 0.02 to 14.45 kg. The relationship of AL as a function of RCD indicated a strong fit using a power model (RMSE = 0.969; Figure 4A). We found a highly significant linear relationship between AS and AL (R2 = 0.990, P < 0.001; Figure 4B). Canopy area strongly predicted AL following a linear relationship (R2 = 0.992, P < 0.001; Figure 4C). The ratio of estimated canopy area (CA1 and CA2) to AL was between 0.13 and 0.70 m2 m−2, indicating that observed crown spread represents only a fraction of total projected leaf area. Among the five trees sampled for AL we found a mean Huber value (AS:AL), of 4.04 (SD ± 1.84) cm2 m−2.
(A) Total projected leaf area as a function of RCD. Solid black line and dashed red line indicate model fits from our study and Pangle et al. (2015) respectively. (B) Relationship between total sapwood area and projected leaf area. (C) Total projected leaf area as a function of the canopy areas CA1 (filled symbol, solid line) and CA2 (open symbol, dashed line). Shaded areas denote the 95% confidence interval of the linear models.
Discussion
Aboveground Biomass
This short-stature tree typically ranges from 3.0 to 15.5 m in height with stems up to 76 cm in diameter depending on site conditions, although most individuals are considerably smaller, standing <10.7 m tall and <46 cm in diameter when mature (Burns and Honkala 1990). The largest tree we sampled was 6.7 m tall and 22.7 cm in diameter. We found that several attributes were strong predictors of P. edulis AGB, described using both linear and power models. Our study supplements the existing diameter-dependent allometric literature for this taxon: the P. edulis model from (Darling 1967) (n = 10 trees), the P. edulis model from Grier et al. (1992) for predicting total biomass in combined young and mature stands (n = 15 trees), the Jenkins et al. (2003) generalized model for softwood pine (n = 34 species), and the woodland Pinaceae model from Chojnacky et al. (2014) aggregated from five equations (n = 2 for P. edulis, n = 1 for P. monophylla, and n = 1 for P. cembroides). Finally, the Sprinkle and Klepac (2015) fresh biomass model (n= 20 trees) was included while adjusting parameters to metric and converting to dry mass by accounting for a bulk moisture content of 38% reported by the authors. As the Jenkins et al. (2003) model scales from measurements taken at diameter at breast height, we converted from our measurements of RCD following (Chojnacky and Rogers 1999), which excludes any trees <1.3 m tall.
Compared with these previously published AGB models for this species (see Table S3 for details of comparable AGB models), our model estimates using RCD as a predictor variable yields higher estimates for AGB among larger stems. For stems <14 cm, our model is in close agreement with that of Darling (1967), Grier et al. (1992), and Chojnacky et al. (2014). The generalized equation by (Jenkins et al. 2003) for Pinaceae underpredicted P. edulis biomass of this study, although the Jenkins model converges with our own at RCD = 35.2 cm. It is important to note that the two generalized equations derived in Jenkins et al. (2003) and Chojnacky et al. (2014) incorporate those of Darling (1967) and Grier et al. (1992) in addition to allometric equations from other related species. Although we demonstrate the applicability of deriving structural vegetation traits using UAV-based methodology developed for surveys of woodlands (Cunliffe et al. 2016, Kachamba et al. 2016, Krofcheck et al. 2019, Karl et al. 2020), traditional survey techniques may be potentially as efficient in these open woodland settings depending on the scale of application.
Maximum tree height was the strongest overall predictor of AGB (Figure 2B) based on the lowest RMSE among dependent variables we considered. This is an encouraging finding for aerial survey applications, as maximum height can be extracted from point cloud photography supported by suitable terrain models, and we found a near 1:1 correlation with maximum height estimates from photogrammetry and traditional field measurement (Figure 1A), comparable to other recent reports (Birdal et al. 2017, Krause et al. 2019). Similarly, the linear relationship of AGB as a function of canopy area expressed a lower RMSE than that of RCD. As both polygon-derived canopy area (CA1) and lengthwise perpendicular canopy area (CA2) yielded similar values for each of the target trees (Figure 1B), neither estimate of canopy area performs significantly better than the other for predicting AGB (Figure 2C). Note that CA2 can be readily estimated from ground-based measurements using the sighting method (Pretzsch et al. 2015). Unlike RCD, both maximum height and canopy area can be measured using aerial photogrammetry, foregoing the need for expansive ground-level surveys to produce rigorous estimates of P. edulis AGB on a landscape extent of approximately 1–100 ha (Cunliffe et al. 2016, 2021).
Sapwood Area and Bark Thickness
Fresh and dried cross-sectional disks demonstrated tight agreement with our RCD measures of the standing live trees (Fig. S2). We report estimated bark thickness of P. edulis (Fig. S1, Table 1), as this is an important parameter for estimating the inside bark volume for timber merchandising and the resulting bark residue (Marshall et al. 2006, Thomas and Bennett 2014), as well as a functional trait indicator of fire resistance (Lawes et al. 2013, Pausas 2015) and a constitutive defense against wood-boring pests (Franceschi et al. 2005).
For estimates of total AS, we found strong agreement between our power model and a previously published model for P. edulis (Figure 3A; see Table S4 for details of comparable As models). One previous model was derived from n = 16 stem samples collected from a site 28 km further west in New Mexico, ranging in RCD from 1.8 to 15.9 cm (Pangle et al. 2015); the other from n = 22 stems ranging in diameter (~1 m above ground level) from 1.8 to 41.6 cm collected 550 km northwest in southern Utah (West et al. 2008; although this equation has been corrected to account for a typographical error in the original publication, West, personal communication, 2021).
The high level of reproducibility for this allometry increases confidence in the use of either model for nondestructive predictions of AS for future studies. This agreement is notable, although our study includes seven trees above the maximum diameter evaluated in the earlier study (Pangle et al. 2015). All three AS models yield comparable values for diameters <10 cm, although the linear model by West et al. (2008) suggests an underprediction of AS for trees above this threshold. However, the field sites for this study and Pangle et al. (2015) share similar qualities; thus, it is conceivable that AS allometry could exhibit ecophenotypic variation elsewhere throughout the extensive range of P. edulis. The model reported by West et al. (2008) resulted in lower estimates of AS across a greater range of diameters, highlighting the potential geographic variability of this species. Reliable measurements of AS were aided by the presence of blue-stain fungus (Endoconidiophora spp.) throughout the active xylem of several samples, helping make the distinction of the heartwood-sapwood interface (see also Fig. S4A), and differentiation of the heartwood-sapwood interface was aided in other samples by evenly wetting the surface with water and allowing it to air dry for a matter of minutes (Fig. S4B). With regard to measuring sap fluxes, the use of AS is often desirable to gauge the depth of stem sapwood prior to installation of sensor probes, as it has been well established that the rate of sap flux can vary along the radial profile (Ford et al. 2004, Alvarado-Barrientos et al. 2013). To this end, the linear relationship of sapwood depth to RCD can be found in Fig. S5 and accounts for the bark thickness along the stem. Among our single-stemmed trees, the distribution of sapwood around the bole was visually quite uniform, unlike what has been observed for multi-stemmed trees such as J. monosperma or some single-stemmed hardwood species (Aparecido et al. 2019, Cunliffe et al. 2020).
Although canopy area (CA1, CA2) was not found to be the most robust independent variable for estimating AS (R2 = 0.76 – 0.77), we considered it valuable to include here (Figure 2C, Table 1) for studies that may desire a method to assess sapwood areas using aerial photography. As basal sapwood area is a physical trait necessary for upscaling point measurements of sap flux to derive stand-level transpiration (Cermák et al. 1998, Mitchell et al. 2009, Hernandez-Santana et al. 2015), this allometric equation could be used in lieu of RCD for aerial datasets, which may quantify canopy areas.
Leaf Area
Our measured LMA of 263.0 g m−2 in this study is in very close agreement with a previously published LMA of 261.7 g m−2 reported for P. edulis growing under ambient conditions in northern New Mexico over a 3-year period (Grossiord et al. 2017). As direct measurements of total tree leaf area can be exceptionally time consuming, most studies subsample branches of large trees and use aggregated branch diameters to scale to whole-tree AL (Jonckheere et al. 2005, Pangle et al. 2015, Lovynska et al. 2018). Despite the relatively short height of P. edulis, total collection of needles from individual trees was labor intensive, with the two largest trees (RCD = 13.4, 19.3 cm) requiring >20 person hours of effort to process. Although conducting a whole-tree harvest has the potential to provide a more accurate assessment of tree-level total AL, we acknowledge a comparatively small sample (n = 5 trees) does not provide a robust statistical model for widespread application. Regardless of the scaling methodology used, AL is an important functional trait linked to gross carbon assimilation (Murthy et al. 2005, Falster et al. 2011, Stephenson et al. 2014) and transpiration (Vertessy et al. 1995, Ryan et al. 2000). We found that the total projected leaf area of P. edulis was highly correlated with total sapwood area (Figure 4B, Table 1), a persistent allometric relationship that has been identified across conifer and broadleaf species alike (Coyea and Margolis 1992, Meadows and Hodges 2002, Stancioiu and O’Hara 2005), as AL typically dictates functional xylem (Huber 1928, Waring et al. 1982). Interestingly, the xylem anatomy of P. edulis is sensitive to precipitation trends, where drought conditions can alter the AS:AL relationship due to isohydric functioning (Guérin et al. 2018, Hudson et al. 2018). Notably, as trees limit transpiration under water stress, the amount of functional AS required per unit AL is reduced to buffer against cavitation and subsequent hydraulic failure. Therefore, while we report a species-specific value for AS:AL from this study, it is important to acknowledge that this reflects ecophenotypic variation at the time and location this research was conducted. It is also worth noting that our equation for estimating AL from RCD corresponds well with the previously published model (Pangle et al. 2015) within the size range sampled by that study (up to 16 cm RCD), but the models diverge considerably for larger P. edulis trees. Our model intersects with that of Pangle et al. (2015) at an RCD of 10.8 cm, such that our equations produce lower AL estimates below this threshold and substantially higher estimates above it (Figure 4A). The largest stem diameter sampled by Pangle et al. (2015) was 15.9 cm; thus, it is unlikely to be representative for whole trees greater than that size, whereas the same limitation applies to our study for trees >19.3 cm.
Conclusion
Accurate estimates of size dependent tree attributes are fundamental to studies of ecosystem structure and function. Allometric relationships provide a robust means to scale cost and labor prohibitive measurements such as AGB across stand and regional scales, while also lending themselves to applications using existing forest inventories. This is of particular importance among trees such as piñon pine, which are abundant on the landscape and serve as carbon sinks. In this study, we destructively harvested twenty trees with a range of dry mass from 0.03 to 171 kg. Compared with previously published studies, our model for estimating AGB from diameter indicates a greater amount of biomass among stems in the larger size classes, although is in tight agreement with existing models for RCD < 14 cm. However, we present models derived from maximum height and projected canopy areas that are indicated as being more robust than the frequently used RCD metric. We also produced a model for estimating AS from RCD, which has validated the utility of an existing model for this species and may serve to upscale estimates of transpiration. Among five trees subsampled throughout our size range, we completely harvested all foliage, deriving several relationships proportional to total leaf area. Altogether, this study has demonstrated the applicability of using aerial photography to derive physical attributes of this widespread and ecological important tree species, while improving on, and in some cases validating, existing models pertaining to allometry of P. edulis across the American Southwest.
Acknowledgments
We thank Tom Carrol and the Deer Canyon Preserve for their permission to conduct research on their land, Stephanie Baker and Joe Bradley for essential logistical support, and Nick Smith, Karen Anderson, Dominic Fawcett, and Nathan Robertson for assistance with data collection.
Funding
This research was funded by the US National Science Foundation (DEB #1557262) awarded to M.L. and Natural Environment Research Council (NE/R00062X/1) awarded to A.M.C.
Author Contributions
C.D.M and A.M.C. conceived the research idea and developed the experimental design. M.L and A.M.C. acquired funding. C.D.M., A.M.C., and F.B. undertook the investigation. A.M.C. and F.B. collected and analyzed the UAV data. C.D.M. processed disk samples and conducted needle harvesting. C.D.M. and A.M.C performed the analysis and wrote the manuscript with input from all authors.
Conflicts of Interest
The authors declare no conflicts of interest.
Data Accessibility
Data and code are available at https://zenodo.org/badge/latestdoi/307378812.
Literature Cited
Author notes
These authors contributed equally to this work.
