Predicting Maximum Tree Heights and Other Traits from Allometric Scaling and Resource Limitations

Terrestrial vegetation plays a central role in regulating the carbon and water cycles, and adjusting planetary albedo. As such, a clear understanding and accurate characterization of vegetation dynamics is critical to understanding and modeling the broader climate system. Maximum tree height is an important feature of forest vegetation because it is directly related to the overall scale of many ecological and environmental quantities and is an important indicator for understanding several properties of plant communities, including total standing biomass and resource use. We present a model that predicts local maximal tree height across the entire continental United States, in good agreement with data. The model combines scaling laws, which encode the average, base-line behavior of many tree characteristics, with energy budgets constrained by local resource limitations, such as precipitation, temperature and solar radiation. In addition to predicting maximum tree height in an environment, our framework can be extended to predict how other tree traits, such as stomatal density, depend on these resource constraints. Furthermore, it offers predictions for the relationship between height and whole canopy albedo, which is important for understanding the Earth's radiative budget, a critical component of the climate system. Because our model focuses on dominant features, which are represented by a small set of mechanisms, it can be easily integrated into more complicated ecological or climate models.

Our framework for predicting height is based on a tree's ability to collect sufficient water and sunlight to meet its basal metabolic needs without exceeding the availability of these resources. These requirements can be summarized by inequalities relating the required flow rate Q 0 , the potential evaporative flow rate Q e , and the available flow rate Q p : Q p − Q 0 ≥ 0 (S1) Eq. S1 ensures that the tree must receive enough water to maintain its basal metabolic flow. It is through this constraint that we couple size to water resources. Eq. S2 states that precipitation must meet or exceed evaporative flow. Eq. S3 ensures that the energy a tree receives from its environment, which is translated into evaporative flow, meets its basal metabolic needs. These statements are summarized by a single condition: Thus, Q 0 and Q p set the boundaries of acceptable flow in an environment. Maximum tree height can then be predicted by finding the largest tree for which this relationship holds. That is, our strategy searches for trees that meet their metabolic needs without exceeding their water or solar resources.
branches at each branching generation, k, whose radius and length are given by r k and l k , respectively, it can be shown that h ≈ l 0 1 − n −1/3 , r k r N = n (N −k)(a/2) , and where, typically, n = 2, and a ≈ 1 is the branching exponent [1]. For terminal units, r N ≈ 0.4 mm, and l N ≈ 4 cm [1]. Note that plant height can also be expressed as Many observed allometric relationships are given in terms of the basal stem diameter, D. In the relationships above D = 2r 0 , and thus D ∝ h 3/2 , which agrees with the findings of reference [3] (h ∝ D 2/3 for large trees). This allows us to re-write scaling relationships for stem diameter in units of tree height. For example, the basal flow rate of fluid through a tree given by [4] becomes with η 2 = 3η 1 a/2 ≈ 2.7 and β 2 = β 1 2r N 1−n −1/3 l N 3a/2 η1 ≈ 9.2 × 10 −7 (liter day −1 cm −η2 ; for h in cm), given the parameters shown in Table 1. It should be noted that these values are found from measurements of trees in situ which may not be operating at the absolute basal metabolism. Future work should focus on determining the minimal requirement of a tree and this will likely adjust the value of β 1 but not the exponent η 1 . This branching architecture also defines several geometric characteristics of the canopy and root systems required for calculating Q p and Q e .

Relationships governing the available flow Q p
The available flow rate is determined by the incoming rate of precipitation, p inc , along with the capture area and efficiency, γ, of the roots. We define the capture area in terms of the radial extent of the root system, r root , so that We find that on average Q p matches Q 0 for the tallest trees taking γ = 1/3, that is, solving for Q p (h max ) = Q 0 (h max ) yields a mean value of γ = 1/3, and this is the value that we use in our model. Calibrating γ in this way allows us to later solve for height using only the single limitation of Q p . Adding variations in soil type and hydrology to calculate γ locally is an important area of future research. The root vascular system is assumed to obey similar constraints as the canopy except that gravitational biomechanical limitations, important for canopy branches, are not important for roots. Roots may exhibit a different architecture, or overall shape, from that of the canopy, but, nevertheless, obey similar hydrodynamic and minimization constraints as the above-ground branching network. Branches and stems grow in diameter through a process of older vascular tubes becoming heartwood and new vascular tubes forming new sapwood. Since vascular tubes run the entire length of the tree, the above-and below-ground networks are not independent. As such, branch tube formation and mortality (above-ground sapwood and heartwood formation) impose constraints on root structure. This is supported by data which show that the branching exponent a for the root system has a mean value close to unity at each branching generation, similar to that for above-ground branches [5]; (our reanalysis of data for the two species presented in Fig. 1 of ref. [5] gives a = 1.03 ± .33 and a = 1.12 ± .45 for n = 2). Combining all of these constraints results in the scaling and topology of the root system being essentially identical to that of the canopy, although the detailed architectures may be different. Thus, the total mass of the roots M R should scale isometrically with the mass of the stems M S , a prediction supported by the data: [3,6].
For optimal access to moisture the roots maximize the volume of soil occupied. The radial extent of the roots should therefore be equal to the longest tubes found in the root system. For stems, the length of a tube running from the base of the trunk to the leaves is given by L S ≡ h = β 4 M 1/4 S [1] so, given the above argument for roots, it follows that L R ≡ r root = β 4 M 1/4 R . Thus the radial extent of the roots is given by Later we relax these assumptions regarding the scaling of the root system with overall tree height and show how optimizing the root scaling exponent can significantly reduce the error between observed maximum tree height and our predictions.

Relationships governing the evaporative flow Q e
Our relationship for the evaporative flow rate of a tree is determined by the interaction of the canopy with the local meteorological conditions. Here we will analyze the energy budget of the canopy, and in doing so, determine how the allometric properties of the canopy govern various heat fluxes.

Energy Budget
Our analysis of the energy budget follows ref. [7]. Conservation of energy requires: where R abs is the total rate at which radiation is absorbed by the canopy (W) and each of the other terms is an energy flux (W m −2 ) away from the canopy multiplied by an effective area, a i=g,j,f , over which that flux occurs (the scaling of these areas are described below). L is the emitted thermal radiation, H is the sensible heat loss, λE is the latent heat loss with λ being the latent heat of evaporation (J mol −1 ) and E the evaporative molar flux (mol m −2 s −1 ). In the definitions below it can be seen that each of the flux terms depends on the leaf temperature, T L (in degrees Celsius; temperatures in degrees Kelvin will be denoted using boldface T). It is generally desirable to eliminate T L from the system of equations in order to focus on common measurements such as the local air temperature. The Penman-Monteith method gives an approximation of T L by linearizing each component of the energy budget in terms of ∆T = T A − T L and solving for T L , where T A is the air temperature. Assuming ∆T << T L and ∆T << T A we can approximate the emitted thermal radiation as where l is the emmisivity of the leaf, σ (W m −2 K −4 ) is the Stefan-Boltzmann constant, c p (J mol −1 C −1 ) is the specific heat of air, and g r = 4 l σT 3 a /c p (mol m −2 s −1 ) is the radiative conductance (we will discuss various conductances below).

The sensible heat loss is similarly defined in terms of ∆T and is given by
where g Ha (mol m −2 s −1 ) is the heat conductance of air [7]. For the latent heat loss we make the approximation is the Tetans formula for saturation vapor pressure (kPa) with b 1 = .611 (kPa), b 2 = 17.502 (dimensionless), and b 3 = 240.97 ( • C) [7].
Solving this system of equations we obtain which, using the approximation in Eq. S21, leads to an equation for the evaporative flux The total volume flow rate, Q e , is obtained by multiplying E can by the area a f , and then converting from a molar rate to a volume flow rate. Thus Q e = a f Eµ w /ρ w given the molar mass, µ w (kg mol −1 ), and density, ρ w (kg m −3 ), of water. This can be written in full as

Allometric canopy geometry and solar radiation
In this section we derive the rate of absorbed solar radiation, R abs , a necessary component of Q e . This absorption rate represents the relationship between available solar resources, plant structure, and, ultimately, tree height. The collection of solar radiation by a plant can be described by coefficients for canopy transmission, τ , canopy reflection, ξ c , and soil reflection ξ s all of which can be encapsulated as the canopy absorption coefficient [8]. A commonly used model [7,8], which we find connects easily to the branching architecture discussed above, assumes that leaves are uniformly distributed over a spheroidal canopy of surface area, S can , defined by the semi-axes h can = h can /2 and r can . In that case, the transmission coefficient is given by a classic absorption formula τ = e − √ αK(ψ)LAI (S27) [7,8], where α is leaf absorptivity and LAI is the leaf area index, a measure of the total leaf area per unit ground below the canopy, which is given by The extinction coefficient, K (ψ) = 2P can /S can , describes the ratio of canopy projected area, P can , to canopy surface area, S can , and is a function of the solar zenith angle ψ. Given a spheroidal canopy, it follows that P can = πr can (r can cos θ + h can sin θ tan −ψ), and [8,9] where θ = arctan (−h can /r can tan ψ) is the angle defining the point of ray tangency on the canopy surface and h can = h can /2. The transmission coefficient, τ , can be calculated by combining the above equations. The total solar radiation captured by the canopy is then given by where R inc (W m −2 ) is the incoming radiation per unit area (normal to the ground), and ξ c and ξ s are described below.

Heat flux areas
Our canopy energy budget is derived by scaling up from the budget for a single leaf whose properties are assumed invariant across trees of different size. In this picture, changes in energy fluxes are predominantly due to allometric relationships for the canopy shape and number of leaves. Leaves are assumed to act in parallel, and thus the total area for sensible heat loss is a j = 2a L , where a L is the total one-sided area of all the leaves on the tree, given by a L = a l n N ∝ h 3 ∝ M 3/4 [1]. The evaporative flux, on the other hand, occurs over the total area of the stomatal openings in the leaves, a f = 2a L δ s a s , where δ s is the density of stomata on a leaf, and a s is the area of a single stoma. If L is the thermal radiative flux emanating from the tree, then a good approximation for the effective radiative surface area is the total one-sided leaf area, so a g ≈ a L . Since S can ∼ h 2 , whereas a L ∼ h 3 , leaves will overlap on the canopy surface for sufficiently large trees so that the above approximation for a g becomes less accurate. We tested other schemes for a g , such as using S can , or the minimum of S can and a L , and found similar results to the ones presented here.

Scaling of the canopy radius and height
In previous work, it has been shown that the canopy radius scales as r 3 can ∝ r 2 0 [10], or which follow from Eqs. S5 and S6. The scaling with height is consistent with empirical data where it has been shown that r can = β 5 h η3 , with η 3 = 1.14 , and β 5 = 35.24 (cm m −η3 ; for h in meters) (see SI of ref. [11]). Unlike r can , we were unable to locate data characterizing the scaling of h can with total tree height, h. Following ref. [12], we assume that the canopy scales isometrically so that h can = β 6 h, where β 6 is a constant. A lower bound for the canopy height can be obtained from a configuration in which all of the branches of the tree emanate directly upwards from the top of the trunk leading to h can ≥ h − l 0 ; clearly, an upper bound is given by the total height of the tree: for n = 2. Note, further, that, from data the diameter of the canopy, d can = 0.7h 1.1 , indicating that, generally speaking, canopy height exceeds canopy diameter, h can > d can , suggesting that the canopy forms an approximate prolate spheroid. For the discussion presented in the text, we took β 6 = 1. Varying β 6 between 0.8 and 1 slightly shifted the mean of the distribution in Fig. 1C of the main text.

Canopy Radiation Coefficients
Above, we discussed the derivation of τ , which is one of the radiation coefficients involved in calculating the total radiation absorbed by a canopy. The other two coefficients correspond to the canopy and soil reflection, respectively, ξ c and ξ s . The soil reflectivity can be taken as a constant value which is given in Table 1, but the canopy reflectivity is expected to change based on the size and scattering properties of the canopy. It has been shown that the canopy reflection, as a function of the leaf area index, LAI and extinction coefficient K (ψ) is well approximated by and ξ * c is the deep canopy reflection coefficient, which is constant [7,8]. For small LAI, Eq. S35 is approximately the ground reflectivity ξ s , and for large LAI, it is ξ * c . Combining equations S26, S27, and S35 it is possible to determine the total absorption coefficient of the canopy and this is the calculation used to produce the curve in Fig. 4 of the main text.

Canopy conductances
Each of the above terms depends on the conductances of heat, vapor, and radiation which are described below. Assuming that the canopy surface is exposed to forced laminar convection, the heat conductance in air is given by where Re a = ud/ν a is the dimensionless Reynolds number for air, P r a = ν a /D H is the dimensionless Prandtl number for air, D H (m 2 s −1 ) is the heat diffusivity, ν a (T a ) (m 2 s −1 ) is the kinematic viscosity of air, u (m s −1 ) is the wind speed, d = 0.81 2 (a l /π) 1/2 (m) is the characteristic dimension of a circular leaf of area a l (m), andρ a is the molar density of air given by the Boyle-Charles law aŝ for Kelvin temperature T A and air pressure p a = 101.3e −A/8200 (kPa) at altitude A (m) [7].
We assume that the canopy is in series with the atmospheric boundary layer, so g v = 1/ 1 g vL + 1 gva , where g vL and g va are, respectively, the vapor conductances for the leaf and atmosphere. Assuming forced laminar convection, as before, the boundary layer conductance is where Sc a = νa Dv is the dimensionless Schmidt number, and D v (m 2 s −1 ) is the diffusivity for water vapor in air. The leaf conductance is calculated by considering the stomata as acting in parallel. We model a single stoma as a cylinder of depth z s (m) in which case Fick's law gives for the conductance of a single stoma. Summing over the entire leaf area, the leaf conductance is given by where a s (m 2 ) is the area of the opening of a single stoma, and δ s (stomata m −2 ) is the stomatal density on a leaf.

Trait values
For each tree trait we examined the literature in order to gain a sense of the variation of each trait across many species, plant sizes, and environments and we picked values that were representative of that variation. For several of these values we were able to check for agreement with means from the TRY database [13] which is a collection of multiple databases containing global values for various plant traits across numerous environments and species. For each trait there are several ways to analyze the distribution of trait values produced by the TRY database. The first is to take the straight mean, which we present in Table 1. The second is to recognize that most of these traits appear to be lognormally distributed in which case it is possible to logarithmically transform the data, find the transformed mean, and then take the exponential of that value. Formally this is the median assuming a lognormal distribution and these values are also given in Table 1. These two analyses typically provide a range of possible representative values for each trait, and the values that we used fall within this range.
Relationships between Q 0 , Q e , and Q p for observed tallest trees We tested the model by looking at observed maximum tree heights in different environments to check whether they do indeed adhere to the given resource constraints governing the required, evaporative, and available flow rates Q p (h obs ) ≥ Q e (h obs ) ≥ Q 0 (h obs ). We calculate each of these quantities using observed maximum heights and measured environmental conditions. For trees which have grown to reach the upper bound on height in a given environment we would expect that Q p (h max ) = Q e (h max ) = Q 0 (h max ). Figure S1A-C shows the resulting pair-wise relationship between each of these quantities.
The derived data generally fall along the green line which represents equality between the two quantities being considered. In each plot the scaling of the two flows is independent of the choice of γ, but the overall intercept will depend on this value.

Detailed examination of predicted vs. observed tree heights
The error between predictions and observations, as shown in Fig. 1C of the main text, is slightly bimodal with the secondary mode corresponding to situations where we under -predict tree height. Most trees in the under-prediction mode are found in mountainous regions (elevations greater than 1000m) where there is often high spatial variability in precipitation. Thus the interpolated meteorology may not be representative of the environment experienced at the tree location. We tested this hypothesis by comparing precipitation estimates from the Parameter-elevation Regression on Independent Slopes Model (PRISM) [14,15] to estimates from the North American Regional Reanalysis (NARR) [16,17]. Large discrepancies between the two databases may indicate either differences in the modeling approach used to assimilate individual station data or high spatial variability in the meteorology. Fig. S2A shows the distribution of the discrepancies between the two analyses. The two data-sets are typically in good agreement, with the PRISM data generally being a little larger than the NARR data. When we removed locations where the error between the two precipitation estimates was more than one standard deviation away from the mean then the under-prediction mode was diminished (Fig. S2B) in the distribution of error between observations and predictions of maximum tree height.

Temperature Shifts
A central question in ecology is how ecosystems respond to environmental shifts. Temperature plays a dominant role in controlling evaporative rates in our framework and here we investigate how our predictions for maximum tree height respond to changes in mean annual temperature. We consider a fairly large change in mean annual temperature of ±2 • C. The maps in Figures S3A and B illustrate the resulting percentage change in predicted maximum height. Averaging over the entire continent the maximum tree height decreases by 11% for a +2 • C shift while maximum height increases by 13% for a −2 • C shift. For an increase in temperature nearly all predicted maximum heights will decrease, although the percentage change varies regionally, and for a decrease in temperature nearly all predicted heights will increase. Looking across different environments we find that warmer environments are less sensitive to a ±2 • C shift than colder environments, hence the strong gradient in percentage change with latitude and/or elevation. The northern and mountainous regions are more sensitive to an absolute shift in temperature. However, it should be noted that a ±2 • C shift in a colder environment represents a much larger percentage change in mean annual temperature compared with a warmer environment. Another analysis which accounts for this effect is to change temperature by ±10% in every environment. Doing this we find that warmer environments are more sensitive to a constant percentage change than colder environments (Figs. S3C and D). For the +10% change in mean annual temperature we find that the average change in maximum tree height is −8% across the entire continental United States. For the −10% change maximum heights change by an average of +9%.

Parameter Sensitivity
It should be noted that the empirical exponents used in this study are fits to data and are accompanied by some amount of error. In addition there are several assumptions or derivations that we have made regarding the scaling of plant features for which there is no data. An important question is how robust our model is to small deviations in the values of these empirical or analytic exponents. Yet understanding the covariation of all of the exponents such that they all are in agreement with data is an interesting and complicated problem (see refs. [18][19][20]). For the suite of exponents involved in this paper the mechanisms of covariation are simply not known, and we do not have data for each property. Thus a full-scale, empirically-realistic sensitivity analysis to parameter combinations is not possible. Given these limitations we can however examine the effects of perturbing each exponent independently of the others. Figure S4 shows changes in the median relative error, , of our predictions given a percent change in the dominant exponents of our model (including analytic assumptions or results). The dominant exponents are those which control water acquisition (the root radius scaling, r root ), and those which regulate the canopy energy budget via its geometry and heat flux areas (canopy height, h can , and radius r can , and the thermal, sensible, and latent heat flux areas, a g , a j , a f respectively). Our original assumption is that all heat flux areas scale according to the scaling of leaf area (see above). In Figure S4 we test changes to this leaf area scaling as well as changes to the scaling of each heat flux area individually. In each plot zero represents the parameter value used in our model, and it can be seen that for each parameter there is a minimum in the median relative error (Fig. S4). This minimum often does not occur at the value used for our original predictions, revealing how our model could be optimized to more accurately predict data. In the future it may be possible to optimize all exponents simultaneously in order to achieve greater predictive power. However, our goal here is to test the zeroth order theory based on the available observations. For the canopy height and radius and the scaling of the total leaf area, small changes to the exponents result in small changes to the median relative error. However, the median relative error is very sensitive to the scaling of the root radius which directly affects moisture gathering, and to differences in the scaling of the heat flux areas compared to the overall leaf area. The thermal and sensible heat flux areas contain multiple minima in the relative error and these two parameters have competing effects with one another illustrating the complicated connection between exponents and the importance of parameter covariation for a true assessment of sensitivity. Each heat flux exponent represents a physiological response to the environment, and given the complicated connection of these exponents it is of future interest to test whether our framework can be used to predict changes in exponents across environments as observed by [18].

Improving predictions using parameter optimization
From Figure 1C of the main text, it is apparent that, despite overall good agreement with the model predictions, there are significant outliers representing deviations between observations and our predictions. This is hardly surprising, given the many potential sources of error such as a possible mismatch between meteorological measurements and the actual conditions experienced by the observed trees (see discussion regarding Fig. S2), and processes such as logging, mortality, and fire which are not included in the scope of our analysis. Additionally, the measured tallest trees may not have reached the environmentallydetermined upper bound. As we have seen in Figure 5 of the main text, our model anticipates the observed environmental dependence of optimal plant traits and this implies that using only a single set of plant traits across all environments may also be a source of error. Finally, some of our assumptions regarding unmeasured scaling relationships may be inaccurate and contribute to error within the predictions. Our goal was to capture the central tendencies, or average behavior, given a set of known (both observed and theoretically motivated) scaling relationships and average plant traits. However, here we show how our predictions can be improved by finding additional scaling relationships which optimize the model.
In the main text we solved for the stomatal density which maximized height at a given location and set of environmental conditions (Fig. 5). Looking across all of the tree sites we found a relationship between stomatal density and temperature. Similarly, we can solve for the stomatal density, δ s , which gives Q P (h obs ) = Q E (h obs , δ s , {m}) for every observed tallest tree (given the local meteorological conditions {m}). Doing this we find that stomatal density should scale with tree height as δ s ∝ h −.71 (assuming a constant stomatal area). Incorporating this scaling into our model then improves our predictions as illustrated in Figure S5B. In our model we also use a constant value of the root absorption efficiency γ and we can also optimize the scaling of this trait alone where we find that γ ∝ h −.78 also reduces the error in the predictions (Fig. S5C). This demonstrates the potential utility of the model in anticipating scaling laws for which measurements do not already exist. However, it should be noted that these two analyses optimize scaling with respect to only a single tree feature one at a time, and, for comparison with future data, it is likely necessary to co-optimize all exponents of relevance. This analysis is an example to illustrate that adding additional scalings to the model can improve predictability and reduce the variance of the deviations; however finding the realistic set of exponents subject to multiple limitations is a subject of ongoing research and requires a more advanced optimization study --- * It should be noted that our model relies on the bulk property of "Stomatal Area per Leaf Area" for calculations requiring an evaporative flux area a f . Thus the sub-parameters as and δs need not agree with the TRY database as long as the bulk property δsas does.