Woody Surface Area Measurements with Terrestrial Laser Scanning Relate to the Anatomical and Structural Complexity of Urban Trees

: Urban forests are part of the global forest network, providing important beneﬁts to human societies. Advances in remote-sensing technology can create detailed 3D images of trees, giving novel insights into tree structure and function. We used terrestrial laser scanning and quantitative structural models to provide comprehensive characterizations of the woody surface area allometry of urban trees and relate them to urban tree anatomy, physiology, and structural complexity. Fifty-six trees of three species ( Gleditsia triacanthos L., Quercus macrocarpa Michx., Metasequoia glyptostroboides Hu & W.C. Cheng) were sampled on the Michigan State University campus. Variations in surface area allocation to non-photosynthesizing components (main stem, branches) are related to the fractal dimension of tree architecture, in terms of structural complexity (box-dimension metric) and the distribution of “path” lengths from the tree base to every branch tip. The total woody surface area increased with the box-dimension metric, but it was most strongly correlated with the 25th percentile of path lengths. These urban trees mainly allocated the woody surface area to branches, which changed with branch order, branch-base diameter, and branch-base height. The branch-to-stem area ratio differed among species and increased with the box-dimension metric. Finally, the woody surface area increased with the crown surface area of the study trees across all species combined and within each species. The results of this study provide novel data and new insights into the surface area properties of urban tree species and the links with structural complexity and constraints on tree morphology.


Introduction
Urban trees provide a wide range of important ecosystem services, including temperature regulation, carbon storage, water purification, air pollutants uptake, biodiversity, etc. [1][2][3][4][5][6][7][8]. Trees in urban areas are growing under very different conditions than forestgrown counterparts and it is important to study their architecture and their physiological performance [9] to increase our understanding of their adaptability to urban environments [10].
In this study, we used TLS technology to produce novel woody surface area data for urban trees. The objectives of the study are the following: (i) to measure the total above-ground woody surface area of urban trees of different species; (ii) to examine the above-ground woody surface area allocation into stems and branches of different sizes and positions in tree crowns (branch orders, branch-base diameter, and branch-base height classes); (iii) to quantify the relationship between tree woody surface area and tree fractalstructural complexity (i.e., the box-dimension metric) and different Euclidean measures of tree architecture (i.e., metrics that account for the length of all paths from the tree base to each branch tip, and crown surface area).

Urban Tree Data
We chose open-grown urban trees as our object of study, since we expected to get better TLS-based measurements of tree attributes, without occlusion from neighboring trees, but also since we expected that the low or complete absence of competition from tree neighbors in cities should allow trees to better express their inherent fractal character in terms of structural complexity [25][26][27][28].
Fifty-six trees, of three species, that represent different tree functional types were sampled on the Michigan State University campus: 18 Gleditsia triacanthos L. trees (compoundleaved, deciduous angiosperms), 15 Quercus macrocarpa Michx. trees (entire-leaved, deciduous angiosperms), and 23 Metasequoia glyptostroboides Hu & W.C. Cheng trees (needleleaved, deciduous gymnosperms). The trees were selected to cover a large range of sizes (Table 1, from 10.9 to 122.2 cm stem diameter at breast height (DBH)). Since all species are deciduous, we were able to study their total woody surface area during the leaf-off period ( Figure 1). Of particular interest was M. glyptostroboides, which was selected since we wanted to study the structure of an urban-grown, needle-leaved gymnosperm, that could be scanned (alive) in a leaf-off condition; it is difficult to get complete, non-occluded scans of the stems and branches of needle-leaved evergreen species [68].  Laser scanning of the G. triacanthos and Q. macrocarpa trees was accomplished, with leaves-on, in July and August 2019. The M. glyptostroboides trees were scanned, with leaves-on, in August 2020. The leaf-on tree point clouds allowed the estimation of the crown surface area of the trees (see specific methods below). The same trees were also laserscanned in leaf-off condition, in January, February, and March 2020, in order to estimate their woody surface area (see specific methods below). Before re-scanning the study trees, we confirmed that they were not pruned by Michigan State University arborists between July 2019 and August 2020. Therefore, pruning did not impose any bias in the estimation of the woody surface area and crown surface area of the study trees, during the study period. Of course, past pruning events could have affected the tree structure observed.

Terrestrial Laser Scanning and Point Cloud Processing
All trees were scanned with the FARO Focus 3D X 330 terrestrial laser scanner (FARO Technologies Inc., Lake Mary, FL, USA), which operates with a laser light of 1550 nm wavelength, typical beam divergence 0.19 mrad, and range of 0.6-330 m. Each individual tree was scanned with high resolution from a minimum of four different directions and distances in order to minimize occlusion effects, and five reference target-spheres were placed around each focal tree to spatially reference all scans and create a single point cloud for each tree following the field scanning protocols suggested by Wilkes et al. [70]. The first two scans of each tree were conducted in opposite directions from distances that allowed the top of the tree to be clearly visible. The other two scans were conducted in opposite directions (perpendicularly to the first two scans) from a closer distance to the focal tree to better capture its branching architecture and get closer views of the main stem. For large trees with complex crowns two or three additional scans were conducted below the tree crown to capture more dense point clouds of the branches. All laser scans were conducted when there was little or no wind. Spatial co-registration and noise-filtering of all scans was automatically performed using the software SCENE 2019.2 (FARO Technologies Inc., Lake Mary, FL, USA, 2019.2). Using the same software, each tree was manually separated from the point cloud of the urban site background. This process was judged as an accurate alternative to an automatic segmentation process [48].

Tree Reconstruction from Quantitative Structure Models
Quantitative Structure Models (QSMs) [71][72][73][74] describe the three-dimensional architecture of trees by fitting cylinders to a tree's point cloud. QSMs preserve branch and stem topology and provide information about the size, the location, the hierarchy, and the orientation of the branching network. They are currently considered to be the most robust method to model tree volume and tree-architecture [75].
Quantitative structure models were generated from the leaf-off point clouds of the trees, with the algorithm TreeQSM v.2.3.0 [76] (see example in Figure 2). There are two main steps for the tree reconstruction from a single point cloud based on this algorithm. The first step is the point cloud segmentation into stem and branches based on cover sets, and the second step includes the reconstruction of the volume and the surface area of the segments with cylinders [45,77]. TreeQSM generated multiple QSMs for each tree point cloud based on different values for the minimum and maximum size of the cover sets and it finally determined the optimal QSM for each study tree [71]. Based on the optimal QSM parameters the algorithm produced 30 additional QSMs for each study tree, which allowed for quantification of the variation of the modeled tree variables, due to the inherent stochastic component of TreeQSM [71]. The definition of the main stem of a tree according to TreeQSM is based on three criteria: (1) the main stem extends near the top of a tree, (2) it goes almost straight up, and (3) it is not too curved which means that the ratio of the stem length to the stem base-tip distance, must be the minimum among all candidate main stems (Raumonen, P. personal communication, 2 June 2020). After the main stem has been determined, the first-order branches (i.e., branches attached to the main stem) are defined based on the following criteria: they are the farthest-reaching candidates, with the ratio of the branch length over the branch base-tip distance to be less than 1.2, and the branch base-tip distance to be over 75% of the maximum. The branch length/distance ratio will iteratively increase if no candidates with a ratio equal to 1.2 exist. The second-order branches are attached to the first-order branches, and they include the candidates with the longest branch base-tip distance. Branches of higher order are defined following the same rules attached to the second, third, fourth order, etc. (Raumonen, P. personal communication, 2 June 2020). Figure 2B,C shows a color-coded example of a result of this procedure.

Tree Woody Surface Area Computation
From the optimal QSMs of the leaf-off scans of the urban trees, their total woody surface area (WSA, m 2 ) was computed as the total surface area of the cylinders that were fitted to the point cloud of each tree ( Figure 2). Note: total "woody" surface area in this study is technically the surface area outside of the bark tissues. Next, the total WSA of each tree was separated between the main stem woody component and the branches woody component. The branch WSA was further analyzed by branch order (there was a maximum of eleven branching orders across all species combined), by branch-base diameter classes of 1 cm size from the diameter of the cylinder at the base of a branch (there were 48 classes across all species combined), and by branch-base height classes of 1 m size, based on the height from the base of a tree to the base cylinder of a branch (there were 25 classes across all species combined; described in detail in the Results Section).

Computation of Other Tree Structural Metrics
We wanted to study how WSA related to other published metrics of tree architectural complexity. Smith et al. [78] examined a metric called the "path fraction", which is the mean length of all "paths" through the branching network, from the stem base to all branch tips, divided by the maximum path length. Path lengths are Euclidean metrics of tree structure and they can be calculated from a QSM based on the lengths of the interconnected cylinders whose topological hierarchy is preserved in a QSM. However, we did not simply use the QSM for each tree to compute Smith et al.'s path fraction, but instead we looked at various statistics from the distribution of path lengths, to take advantage of the rich data provided. This included the quantiles of the path lengths (25th, 50th, and 75th percentiles), as well as minimum, maximum, mean, and standard deviation of path lengths.
We computed the box-dimension (D b ) metric [27], as a direct measure of above-ground structural complexity, calculated directly from the leaf-off point cloud of each tree. D b has the advantage of not having to apply a QSM to the data, it uses only the original tree point cloud. The D b metric takes into account the number of boxes that are needed to encapsulate all points of a laser-scanned tree, and how the number of boxes varies with the ratio of the box size to the original box size, which is defined as the smallest box that encapsulates the whole tree [15]. The smallest box encapsulating the entire tree point cloud is the so-called "upper cut-off", as it represents the largest box applied to the tree for counting the number of consecutive boxes needed. Consecutive boxes always have half the edge length of the previous box so that eight of them fit exactly in the initial box. The smallest box-size among all boxes is the so-called "lower cut-off", and it was defined to be 10 cm in this study ( Figure 3A). It is a very liberal estimate of the maximum distance between two neighboring laser points at any given location in the tree. The "lower cut-off" must ensure that no virtual box is considered empty only because it fits in the "unsampled" space that was not reached by any laser beam of the laser scans. This "unsampled" space may be the result of the diverging beams emitted from the scanner leaving unscanned areas at greater distances to the scanner or simply due to occlusion effects in the tree. (B) The log-log plot for the quantification of the box-dimension metric for the same tree. The regression line slope is the box-dimension of the tree, i.e., D b = 2.05. N is the number of boxes required to capture all points of the tree point cloud, s is the size of the length of each box, and s_initial is the size of the length of the initial box that encapsulates the whole tree. The 95% confidence interval has been plotted around the regression line. D b is equal to the slope of the least-squares line when the logarithm of the number of boxes is regressed against the logarithm of the inverse of the size of a box relative to the size of the initial box [15,27] (Figure 3B). D b , which is unitless, takes values between one and three. Values smaller than one are only possible if the "lower cut-off" has not been properly chosen (i.e., mean distance between neighboring points is greater than the edge-length of the smallest box). Values of three (or greater) are not possible in reality, since it would imply that a tree is a solid cube. D b values close to but smaller than three imply trees with greater crown complexity and "space-filling character", whereas, a perfectly cylindrical stem without branches would have D b equal to one [27]. Both the path fraction [78] and the D b [27] metrics are meant to capture the fractal-like nature of trees [42][43][44], which should explain a portion of the variation in their WSA.
Finally, the crown surface area of the study trees was computed as the convex hull from the leaf-on laser points of a tree's crown using Heron's formula to quantify the triangles that create the surface of this hull [24]. In this study, it refers to the photosynthetically-active surface area of a tree [15,16].

Statistical Analyses
All statistical analyses for this study were done with custom coding and available packages written in the R software language [79]. Regression analysis was used to relate the total WSA, and the branch to stem WSA, with the metrics of crown complexity and tree architecture (see Section 2.5). Correlation strengths were quantified with the Pearson correlation coefficient (r) and the statistical significance of the relationships was assessed at α = 5% level.
The total WSA of the trees was modeled as a power function of the metrics described in Section 2.5 in order to explore the relationship between WSA and these metrics. The power function form was selected since it had a better fit to the data compared to the linear model form and since power functions better describe the multiplicative processes of tree allometry (e.g., WSA allometry), and they are scale-invariant [80]. The species was added in the candidate models as a random, grouping variable that influences the exponents of the predictor variables.
The mixed-effects model is of the form: where WSA is the total woody surface area (m 2 ) of the trees, b is the normalization constant, D b is the box-dimension (unitless), L is one of the path length metrics in meters that were described previously, c is the scaling exponent parameter of the box-dimension (fixed effect), d is the scaling exponent parameter of the path length metrics (fixed effect), and S c and S d are the species random effects added to the candidate models to modify the c and d parameters, respectively, and they have three levels (i.e., G. triacanthos, Q. macrocarpa, and M. glyptostroboides). The error term (ε) has a multiplicative structure, which is additive on a log-log scale. Assumptions of variance homoscedasticity and error normality were checked by plotting the model residuals against the fitted values, and the Q-Q plots and the histograms of the model residuals. The "nlme" function of the linear and nonlinear mixed effects models ("nlme") R package [81] was used to fit the models. The best models were selected considering both the coefficient of determination (adjusted R 2 ) and the Akaike information criterion (AIC). A one-way analysis of variance (ANOVA) test, with unequal variances, was used to evaluate differences in the mean value of the branch to stem WSA ratio across the three species combined (i.e., G. triacanthos, Q. macrocarpa, and M. glyptostroboides). A one-way ANOVA test was also used to evaluate differences in mean WSA of branches per branch order, per branch-base diameter class, and per branch-base height class, across and within the above-mentioned species. For these tests, the WSA of all branches belonging to different classes for every study tree was considered. In all ANOVA tests, the normality of the data in each group was checked with Q-Q plots, and significant differences in group means were assessed at α = 5%. Finally, the coefficient of variation was used to quantify the uncertainty in estimating total WSA from the consecutive QSM reconstructions of the same point cloud of a tree.

Estimated Total and Component Woody Surface Areas
Basic tree measurements and surface areas computed for the study trees are shown in Table 1, along with other tree statistics (discussed later). The data show that the study trees varied widely in their surface areas and other metrics of complexity.
Branches comprised the greatest portion of WSA of the urban trees studied. The branch to the main stem woody surface area (BMS) ratio ranged between 4.3 and 38.6 with the mean value of 16.3 across all study trees combined. ANOVA showed that the mean BMS differed significantly among the three species, i.e., M. glyptostroboides (MEGL), G. triacanthos (GLTR), and Q. macrocarpa (QUMA) (p < 0.001). G. triacanthos had the highest mean BMS value compared to the trees of the other two species: Mean BMS GLTR = 24.1, mean BMS QUMA = 13.4, and mean BMS MEGL = 12.1. Furthermore, a strong positive relationship was found between the BMS ratio and the D b metric of the trees (Pearson's r = 0.6, p < 0.001).
The median branch order was five across all study tree species combined (range 1 to 11, Table 1), with M. glyptostroboides showing fewer branch orders than the two angiosperm species (median 4, range 1 to 9). ANOVA showed that the branch woody surface area (BWSA) significantly differed among the different branch orders across all study tree species combined and within each species (p < 0.001). BWSA was mainly accumulated in lower branch orders and the distribution of surface area was positively skewed (Figure 4). Second-and third-order branches supplied the greatest amount of BWSA across all study tree species combined, and for Q. macrocarpa trees ( Figure 4A,C). BWSA came mainly from second-, third-and fourth-branch orders in G. triacanthos trees ( Figure 4B), and from lower order (first, second and third) in M. glyptostroboides trees ( Figure 4D). Examination of BWSA by branch basal diameter ( Figure 5) indicated that the BWSA followed a positive skewness, but with a somewhat bimodal distribution (except for the M. glyptostroboides trees). Medium-sized branches (between 4 and 11 cm base diameter approximately) and large branches (more than 35 cm base diameter approximately) accumulated much of the BWSA, while small branches and twigs (less than 4 cm base diameter), though numerous, accumulated a relatively small portion of the BWSA ( Figure 5). ANOVA confirmed that the BWSA differed statistically among the different branch-base diameter classes across all species combined and within each species (p < 0.001). Some large trees showed very large branches, with a basal diameter greater than 35 cm ( Figure 5). These "branches" were actually large forks in the stem, common to large, open-grown, urban trees (See Figure 1), which the TreeQSM algorithm defined as branches. At a major fork, the QSM determines the longest, straightest stem to the top of the tree as main stem (see details in the Methods Section above), and calls the others branches.
Since the topology of the trees is captured by the QSM, we were also able to examine how the surface area was distributed vertically in the trees ( Figure 6). The BWSA differed statistically (as assessed with ANOVA) among the different branch-base height classes across all study tree species combined and within each species (p < 0.001). Graphical analysis ( Figure 6) shows a parabolic distribution of relative BWSA for all study tree species combined peaking near the midpoint of the crown (0.5 on the y axis in Figure 6A). Relative BWSA peaked higher up in the tree for G. triacanthos trees ( Figure 6B), about the midpoint for Q. macrocarpa trees ( Figure 6C), and below the midpoint for M. glyptostroboides trees ( Figure 6D).

Uncertainty Analysis of the Estimated Woody Surface Areas
The coefficient of variation of the WSA of the trees indicated that the uncertainty due to the consecutive QSM reconstructions of the same point cloud of a tree was on average 2.4% of the mean WSA per tree across all study tree species combined, and G. triacanthos trees had the highest uncertainty (on average 2.7% of the mean WSA per tree, Table 1). The distribution of the coefficient of variation of the WSA of the trees was positively skewed across all study tree species combined and within each species and bimodal for the G.

Relationships between Woody Surface Area and Metrics of Tree Architecture and Structural Complexity
Significant, positive relationships were found between the WSA of the urban trees, and the D b metric, and the different metrics that account for the length of all paths from the tree base to each branch tip (Table 2, Figure 8A-H). The strongest positive relationship was found between the WSA of the trees and the 25th percentile of path lengths (Pearson's r = 0.87, p < 0.001, Figure 8F). However, the relationships between the WSA and the 25th percentile of path lengths, the mean path length, and the 50th percentile of path lengths, were not very different ( Figure 8). The best and most parsimonious predictors of WSA (Equation (1)) were the combination of the D b metric and the 25th percentile of path lengths with species effects ( Table 2). The correlation between the D b metric and the other predictor variables in each model of WSA (Table 2) was not statistically significant (i.e., p > 5%). Table 2. Woody surface area models with the highest adjusted R 2 and lowest AIC values among all candidate models fitted to the data. Tree woody surface area (WSA) was modeled as a power function of different predictor-combinations (Equation (1)), including box-dimension (D b ) and various statistics of path length (L), i.e., mean and the 25th, 50th, and 75th percentiles of path lengths. The character "|spp" denotes that species was added as a random effect, modifying the exponent of each predictor variable in the model. The best model by each statistic is highlighted in bold.

Model
Adjusted

Advances in Urban Tree Surface Area Measurement
In this study, we used active remote sensing (TLS) to produce detailed WSA data for urban trees. Measuring the total surface area of the woody parts of trees has been challenging in the past with the only direct method via destructive sampling, which is particularly challenging and undesirable for large trees in urban areas. This study provided the first comprehensive measurements of the total above-ground WSA of urban trees with TLS, including the relative surface area of branch versus stem WSA and complete vertical characterizations of BWSAs for branches of different size and order. TLS has become an important tool used to quantify the three-dimensional structure of trees [35,75] and more accurate measurements of tree surface area may be the most important new advance in tree measurements associated with this technology. Data from mobile laser scanning (MLS) could also be used to study the WSA of trees covering larger spatial scales if occlusion effects in the point clouds are not significant. According to Dorji et al. [82], MLS data can be used to study the structural complexity of trees based on fractal analysis and quantified by the D b metric.
With any new measurement system come new sources of uncertainty. Our field procedure was designed to minimize occlusion effects in the tree point clouds by scanning the study trees from multiple directions and distances (see Section 2.2). This reduced the estimation uncertainty due to cylinder size and cylinder fitting errors in the generated QSMs [45]. The uncertainty in the estimates due to the consecutive QSM reconstructions of the same point cloud of a tree, comprised only a small portion of the estimated WSA across all study trees combined. This was on average 2.4% of the mean WSA per tree across all species combined (Table 1), while very few trees had a coefficient of variation of their WSA larger than 5% ( Figure 7A-D). Therefore, the consecutive QSM reconstructions of a tree provide precise WSA estimates. This does not mean that the QSMs do not introduce bias, such as systematically over-or under-estimating surface areas of different parts of the trees, when identifying them from the point clouds.
Some large study trees had large branches with the max QSM base-diameter greater than 35 cm. Diameter overestimation of large branches (i.e., larger than 40 cm) is usually quite small in the QSMs generated by the TreeQSM algorithm (Raumonen, P. personal communication, 4 March 2021), but parts of forked stems can also be interpreted as branches. In other studies where branches were destructively sampled, an underestimation of 6% in QSM base-diameter was found for branches, with the actual base-diameter greater than 60 cm, while a diameter underestimation of 8% was observed for branches with diameters between 20 and 60 cm [83,84]. Therefore, we are confident in the general accuracy of the BWSA values produced in this study.

Relationships of the Woody Surface Area of Trees Explained by Major Theories of Tree Structure (WBE Model and Pipe Model Theory)
It has been suggested that the variation in branch area is related to the diameter of a branch and its position in the crown [29,32]. It was found here that medium and largesized branches (based on their basal-diameter), of lower branching orders, accumulated the largest portion of the total BWSA. This pattern can be interpreted in the light of the pipe model theory [41] and the WBE model [42], which connect the tree structure with tree physiology. Both theories assume a fractal branching architecture whose vascular structure is an assemblage of tubes that taper from base to tip. Therefore, larger, lowerorder branches accumulate more conducting and non-conducting tubes over their length, resulting in greater cumulative volume that scales with WSA [85]. Similarly, Weiskittel and McGuire [29], found that on average 82% of the total branch surface area in Douglas-fir (Pseudotsuga menziesii) trees was allocated into primary branches (those attached to the main stem). However, Meir et al. [34] found that small branches significantly contributed to the WSA of trees growing in a tropical rainforest.
According to the models produced in this study, much of the variation in the nonphotosynthetic surface areas of urban trees can be explained by a combination of fractalstructural complexity (quantified by the D b metric), and "hydraulic" size (quantified by the Euclidean metric of the 25th percentile of the path lengths; see Smith et al. [78]), which are constrained by the genetics of tree species.
Smith et al. [78] defined the path fraction metric as the ratio of the mean path length to the maximum path length from the tree base to each branch tip, in order to quantify to what extent a real branch network differs from an ideally fractal branch network, such as that described by the WBE theory. In this study, the path fraction was not significantly related to the WSA of the trees, but as expected, significant relationships were found between the WSA and various statistics from the distribution of path lengths. This suggests that the absolute mean, variation, and distribution of path lengths may better help characterize the surface area complexity than the mean relative path length (a.k.a. the path fraction of Smith et al. [78]). Weiskittel and Maguire [29] found that the WSA of Douglas-fir (Pseudotsuga menziesii) trees increased with crown length, which agrees with the positive relationship that was found in this analysis between the total WSA of the urban trees and the 25th percentile of the path lengths, which is the frequency of the short-path lengths that affects crown length.
The WSA of the study urban trees was found to increase with their fractal-like architecture, as quantified by the D b metric [15]. According to the WBE model [42], this pattern could imply efficient respiration rates and sufficient supply for energy demanding units, e.g., leaves, chloroplasts [86], since the inherent fractal character of trees allows them to maximize the scaling of their external surface areas for gas exchange with the atmosphere, while minimizing the internal vascular distances for transferring and allocating the available resources to different organs and tissues [85][86][87].
An important issue to consider, when analyzing the relationships between surface areas and metrics of crown fractal complexity, is whether the observed patterns are confounded with tree size (e.g., DBH, total tree height, etc.). The D b metric is reported to be scale and tree-size independent [15,48], and therefore, we can use it to compare trees of different sizes [15]. Further analysis showed that the relationship between the total woody volume of the study urban trees computed from QSMs and the D b metric was not statistically significant (p > 5%), suggesting that both smaller-and larger-volume trees can be structurally complex. This could mean that architectural changes that occur through the ontogeny of trees, e.g., development of higher order branches and altered stem to branch relationships [48], might explain more complex structures in larger trees, more than their size, per se.

Anatomical and Physiological Implications of Surface Area Allocation Patterns
The surface area distribution found for these urban trees and the theories described in the previous sections have implications for understanding the anatomical structure and physiological function of urban trees, and trees in general. This study enabled not only the computation of the total WSA of trees, but also the analysis of its distribution into different components (stem and branches).
As expected, the branch-to-stem surface area ratio of the trees was found to increase with their structural complexity (as captured by the D b metric), underscoring the contribution of branching to crown complexity [15]. This ratio was found to significantly differ among the studied species, so, in this sense, it describes the resource allocation "decision" of different species to invest in increasing branch versus stem surface areas, as a functional response to urban environments. The squat form of urban trees (i.e., a wide tree crown with a short trunk) gives them mechanical stability against wind loads in cities [88], and reflects the tendency of trees to allocate less resources to growing a taller main stem as the crowding conditions decrease [88,89]. Mäkelä [90] found that "branchiness" of Scots pine (Pinus sylvestris) trees (described as the ratio of total branch cross-sectional area to stem surface area), increased as stand density decreased. Therefore, this pattern appears to hold for trees growing in both rural and urban areas. As such, the branch to stem surface area ratio could be an important component of the plant "structural economics spectrum", which explains species-structural diversity in terms of tree architectural traits along a spectrum balancing light interception, carbon allocation, and mechanical stability [91].
The WSA of trees relates to their respiration rates and captures broad maintenance costs [17][18][19][20][21] and crown surface area refers to the photosynthetically active surface of trees and their energy income [15,16]. Therefore, the strong and positive relationship that was observed between WSA and CSA (across all study tree species combined and within each species) implies that as the respiration rate of a tree increases, its production efficiency should also increase in order to maintain sustainable growth. Otherwise, trees should lose vigor.
For trees in natural forests and plantations, the distribution of branches and foliage is heavily influenced by shading or sheltering from neighboring trees [88], but particularly the need to maintain a positive carbon balance in the leaves. Weiskittel and McGuire [29] showed that branch surface area peaked a bit below the middle of the crowns of Douglas-fir (Pseudotsuga menziesii) trees, since smaller branches near the top had considerably less surface area, while large, lower branches, with greater surface area, tended to dieback over time, due to a reduction in sustaining leaves on branches near the base of the crown. Xu and Harrington [92] found that sub-dominant trees in a plantation of Loblolly pine (Pinus taeda) distributed most of their foliage at the top third of their crowns. These patterns might be expected for trees of species which are less tolerant of shade. Niinemets and Valladares [93] produced numerical tolerance indices ranging from 1 to 5 (1 = very intolerant to 5 = very tolerant) ranking Pinus taeda = 1.99 and Pseudotsuga menziesii = 2.78.
The three urban tree species studied here showed a branch surface area vertical distribution inversely corresponding to Niinemets and Valladare's [93] shade tolerance indices, with Gleditsia triacanthos = 1.61, Quercus macrocarpa = 2.71, and Metasequoia glyptostroboides = 3.00, showing patterns of branch area peaking in the upper-mid, middle, and lower-mid crowns, respectively ( Figure 6). This result was somewhat surprising, since these trees were open-grown and not shaded by other trees. This suggests that self-shading of leaves and branches could be an important element of the branching architecture of open-grown trees [94] along with inherent shade tolerance [95], but also suggests that a mechanism other than maintaining positive carbon balance in leaves might be at play.
Another physiological explanation of this pattern of branch woody surface distribution could be a need to counterbalance optimizing light energy capture with the need to minimize the surface area for heat gain due to incoming solar radiation, and water loss through transpiration. It is well known that warm temperatures and heat islands in cities [96][97][98] cause increased rates of leaf transpiration [12]. Niinemets and Valladare's [93] also published drought tolerance indices for the study species (drought tolerance for Gleditsia triacanthos = 4.98, Quercus macrocarpa = 3.85, and Metasequoia glyptostroboides = 2.38), showing that the more drought-tolerant the species, the more it concentrated the branch surface area towards the upper crown (see Figure 6).
Niinemets and Valladare [93] showed, over many species, that inherent traits of shade and drought tolerance of species were often negatively correlated, so the above pattern was expected. However, this study suggests that, at least for the case of open-grown, urban trees, building a branching architecture that optimizes drought tolerance may be a good explanation for branch surface area distribution. It may be that hydraulic limitations are not only an important force limiting the size and hydraulic architecture of tall trees [99,100], but also for the branching architecture of any tree whose physiological water stress exceeds their photosynthetic capacity. These findings have important implications for the management of urban forests, particularly the selection of species for urban plantings, given expected, continued increases in global temperatures and urbanization.
While surface area data have long been available for leaves [101][102][103][104][105][106][107][108], TLS, in combination with QSMs, can now be used to quantify the surface area of the "woody skeleton" of trees, which plays a vital role in gas exchange with the atmosphere. Tree respiration rates are closely related to their WSA [18][19][20] since respiration of non-photosynthetic tissues mainly occurs in the cambial sheath and the living annual growing rings around the dead heartwood [17]. Nonetheless, Sprugel [109] suggested that a forest stand with high bole WSA does not necessarily have high rates of respiration, so there is still a need to scale up from tree to stand to forest-level process modeling. Respiration still contributes a significant portion of uncertainty related to the carbon budget and offset potential of urban forests [8]. The type of data that was produced in this study could form the basis to develop new process models that describe the carbon balance and growth of urban forests under a changing climate, which has been an important focus of forest process modeling for decades [90,110].

Conclusions
In this study, we demonstrated the use of TLS technology to produce detailed data that quantify the total above-ground WSA of urban trees (first study objective) and we found that the study trees varied widely in their WSA. Furthermore, based on TLS data we studied the allocation patterns of WSA to different components of the woody skeleton of trees, i.e., stem and branches of different order, base-diameter, and base-height classes (second study objective), and we found that the urban trees allocated their WSA mainly to branches, while branch order, branch-base diameter, and branch-base height influenced the observed allocation pattern. The vertical allocation pattern of branch surface area was inversely related to the shade tolerance of the species.
Measuring the WSA of trees with TLS is a non-destructive method that allows explicitly accounting for the above-ground structural-fractal complexity of trees, and it does not rely on any biological assumptions for tree architecture [35], in comparison with previous methods that approximated branch and stem geometry or estimated WSAs from allometric equations [13,19,29,31,32,36,37].
This study showed that WSA is a function of the above-ground fractal-structural complexity of trees, and their "hydraulic" size quantified by different Euclidean metrics of path lengths from the tree base to each branch tip (third study objective). The observed positive relationship between the crown surface area and the woody surface area of the trees (third study objective) implies a physiological mechanism for maintaining a positive carbon balance at tree scale. In general, the type of data produced in this study describes tree surface allometry, and it can be used to develop new or inform existing process models that quantify the growth and productivity of urban forests.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.