Shifts in Foliage Biomass and Its Vertical Distribution in Response to Operational Nitrogen Fertilization of Douglas-Fir in Western Oregon

: Nitrogen (N) fertilization is a commonly applied silvicultural treatment in intensively managed coast Douglas-ﬁr ( Pseudotsuga menziesii (Mirb.) Franco var. menziesii ) plantations. Field trials were established in a randomized complete block design by Stimson Lumber Company (Gaston, Oregon), to test the economic viability of N fertilization on their ownership and to better understand Douglas-ﬁr growth responses. The 23 stands comprising the trials were Douglas-ﬁr dominated, had a total age of 16–24 years, had been precommercially thinned, and had a density of 386–1021 trees ha − 1 . Fertilizer was applied aerially at a rate of 224 kg N ha − 1 as urea during the 2009–2010 dormant season. In the dormant season of 2016–2017, seven growing seasons following application, 40 trees were felled and measured with the objective of assessing crown attributes and aboveground allometrics. Branch-level foliage mass equations were developed from 267 subsampled branches and were applied to the 40 felled sample trees on which the basal diameter and height of all live branches were measured, allowing estimation of both the total amount of foliage and its vertical distribution. A right-truncated Weibull distribution was ﬁtted to data, with the truncation point speciﬁed as the base of live tree crown. The resulting tree-level parameter estimates were modeled as functions of tree-level variables. Stand-level factors not explicitly measured were captured through the use of linear and nonlinear mixed-e ﬀ ects models with random stand e ﬀ ects. Fertilization resulted in more total crown foliage mass in the middle crown-third and caused a downward shift in the vertical distribution of foliage, with implications for feedback responses in crown development and photosynthetic capacity. Deﬁning the morphological responses of Douglas-ﬁr crowns to nitrogen fertilization provides a framework for studying inﬂuences on stand dynamics and should ultimately facilitate improved site-speciﬁc predictions of stem-volume growth. relative distribution


Introduction
Net primary production (NPP) of forest stands can be estimated through quantification of ecophysiological processes, including mechanisms that influence the quantity and photosynthetic efficiency of foliage (e.g., [1]). NPP is simply the difference between gains accrued through net photosynthesis and losses to construction and maintenance respiration, with the net difference measurable as dry plant matter [2]. When attempting to quantify growth responses to spatially varying environmental conditions, changing climate, or alternative silvicultural regimes, identifying the effects on fundamental ecophysiological mechanisms can provide unique insights, particularly if these mechanisms can be integrated with empirical relationships into "hybrid" growth models [3].
If environmental variables relevant to fundamental ecophysiological processes can be measured, estimated, or forecast in a cost-effective manner for operational stands, prediction of tree growth under a range of alternative management activities and environmental conditions should be enhanced [4,5]. Under this scenario, any anthropogenic manipulation of resource availability (e.g., through thinning or fertilization) could be accounted for at a mechanistic level to predict growth responses, facilitating more accurate predictions where spatial or temporal variation in environmental conditions can be adequately characterized or predicted. Reliable quantification of foliage amount and its vertical distribution are key components in the hybrid modeling approach because they are driving factors for the amount of intercepted photosynthetically active radiation and are closely linked to growth distribution and changes in allometric relationships [6,7].
Silvicultural treatments such as thinning or fertilization have been shown to influence foliage production [8][9][10], foliage distribution [11][12][13], and related gross crown dimensions [9][10][11][12][13][14]. Quantifying effects of silvicultural treatments on foliage quantity and distribution therefore incorporate important feedback mechanisms between silvicultural treatments, growth responses, crown and other morphological changes, and associated ecophysiological processes. Recognition of these mechanisms has motived efforts to model the vertical distribution of foliage by fitting continuous probability distributions to foliage mass (or area) binned by vertical crown segment [11]. Several distributions have been explored, including the Weibull [12,13,[15][16][17][18][19][20][21], normal [22], generalized logistic [23], and beta (β) [11,[24][25][26]. The β-distribution offers the advantage of extreme flexibility over a domain formed by the closed interval [0, 1] [11]. Foliage distributions can vary among trees by level of species shade-tolerance [27,28] and factors such as social position and stand density, particularly in even-aged stands [11,13,20]. Incorporating the influence of silvicultural treatments on foliage distribution into growth models quantifies the recognized link between foliage distribution and the vertical distribution of stem increment on individual trees [29,30]. If quantified with sufficient accuracy, representation of this mechanism in growth models can refine predicted responses of stem form and total stem volume. Models of vertical foliage distribution on individual trees can also facilitate simulation of canopy processes, including net photosynthesis over various time scales and net primary production (NPP) on annual cycles. These production measures can supplement measured or predicted site indices in hybrid growth models (e.g., [31,32]). Empirical models have traditionally been limited to conventional inventory data, making them highly dependent on a static index of site productivity [7]. In theory, potential response to silvicultural treatments such as fertilization should be dependent on resource availability as determined by soil and climatic variables. Availability of resources other than nutrients, in particular water, often dominates long-term forest productivity [33], controls interannual variability in forest productivity [34], and even affects the ability of trees to respond to fertilization [35,36]. The importance of water limits the efficacy of nutrient availability alone to predict the growth response to fertilization. Improvements in characterizing site quality through soil attributes and seasonal weather patterns, and the ecophysiological responses to these conditions, are therefore potentially advantageous for predicting growth responses to silvicultural manipulations.
Commercial application of nitrogen (N) fertilizer is a common silvicultural tool for increasing volume production of coast Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco var. menziesii) plantations. Regional fertilization trials have demonstrated that N is the most limiting nutrient to growth in the Pacific Northwest [37][38][39][40][41], but growth response is extremely variable among different stands, sites [40], and geographic regions. This reality has compelled researchers to identify physiologically relevant site and stand attributes that may indicate site capacity for potential growth response. These site and stand attributes have included slope, elevation, forest floor carbon to nitrogen ratio (C:N), and relative stand density [42]. Water availability has also been shown to influence the magnitude of growth response to fertilization (e.g., [35,43,44]), probably due to coupling of water and nutrient uptake, as well as the constraint imposed by water availability on carrying capacity of the site for leaf area [45,46]. Transpiration rates and carbon assimilation by forests are influenced by foliage amount and its vertical distribution [47][48][49]. The following crown responses to fertilization are therefore relevant: (1) an Forests 2020, 11, 511 3 of 27 increase in total foliage quantity [8,9]; (2) a shift in vertical distribution of foliage [24,26,50,51]; and (3) an increase in water use efficiency defined as stem growth per unit of transpired water [44]. Accurate quantification of the relationship between foliage dynamics and growth responses to fertilization should ultimately lead to reliable identification of responding sites. The core of this predictive capacity should be the magnitude and duration of growth responses in Douglas-fir plantations [26].
The goal of this study was to examine the influence of nitrogen fertilization on foliage quantity and its vertical distribution on individual Douglas-fir trees, in part to facilitate future hypothesis tests regarding effects of these foliage attributes on corresponding changes in stem form on fertilized and unfertilized control trees [52]. The specific objectives of this analysis were to: (1) Develop branch-level prediction equations for total branch foliage mass and test for differences between fertilized and unfertilized control branches; (2) Develop tree-level prediction equations for total foliage mass and test for differences between fertilized and unfertilized control trees; and (3) Develop models for describing relative vertical distribution of foliage within individual Douglas-fir tree crowns and test for fertilization effects on vertical foliage distribution. This study did not allow examination of differences among crown classes or site types, but tree-level covariates and tree or stand random effects were introduced to account for these sources of variation and isolate the effects of nitrogen fertilization. More broadly, this study aimed to augment the extensive, historical fertilization studies in the Pacific Northwest implemented by the Nutrition Project of the Stand Management Cooperative (SMC) and its predecessor, the Regional Forest Nutrition Research Project (RFNRP), specifically to further the understanding of fundamental physiological mechanisms driving growth response of intensively managed Douglas-fir plantations to nitrogen fertilization.

Study Area
In 2009 Stimson Lumber Company (SLC) installed a fertilization field trial to test the economic viability of accelerating growth by operational application of nitrogen (N) fertilizer. SLC distributed the installations across their ownership in the northern Coast Range of western Oregon during the dormant season of 2009-2010 ( Figure 1). Implementation of the field trial conformed to a randomized complete block design with twenty-three operational stands serving as blocks. Stands were selected for this study based on relative uniformity (in terms of stand composition and structure) to minimize experimental error, and on geographic distribution to represent the geographic range of SLC ownership. Stands selected for this study met the following criteria: (1) Douglas-fir comprised ≥90% of the stand basal area; (2) All stands were precommercially thinned to a residual density of 565-865 trees per hectare; (3) Birth age from seed germination ranged from 15 to 25 years; (4) Top height growth indicated a 50-year site index from 34 to 37 m [53]; and (5) Access was sufficient for expedient remeasurement of permanent plots. Each stand was divided into two parts, or experimental units, with fertilization randomly assigned to one (hereafter referred to "fertilized" treatment) and no fertilization assigned to the other (hereafter referred to as "control" treatment). Operational terrain largely dictated partitioning of stands (blocks) into two experimental units. The fertilization experimental units were designed sufficiently large to facilitate further splitting into yet smaller experimental units to receive a later second application of fertilizer. The broader objective of testing multiple fertilizer applications was outside of the scope of the analysis presented here. Nitrogen was applied aerially as pelletized urea (46% N) at a rate of 224 kg N ha −1 to all experimental units designated for fertilization during the 2009-2010 dormant season. Application was completed when weather conditions were cool (i.e., <21 • C), and when precipitation was expected within a day or two to reduce the risk of volatilization. The expected precipitation events also maximized the probability that urea prills caught in tree crowns would reach the ground via wind and rain. A total of 562 hectares were fertilized in this study, with fertilized and control experimental units ranging in size from 4 to 44 hectares and from 2 to 16 hectares, respectively.

Data Collection
Permanent plots were established within each experimental unit in the 2009-2010 dormant season (just before fertilization), including three plots in the control experimental units and three or six plots in the fertilized experimental units. Plots were circular and covered 0.04 ha (11.35 m radius). All trees with diameter at breast height (D) > 5.0 cm within the plots were numbered and measured for D (nearest 0.25 cm), total height (H; nearest 0.03 m), and height-to-crown base (HCB; height to the lowest live branch; nearest 0.03 m). Any visible damage or deformity was recorded along with the measurements. Trees within the plots were remeasured during the 2016-2017 dormant season, providing growth responses for the first seven-year period following application. During the 2016-2017 dormant season, two trees from each experimental unit were felled in each of ten stands randomly selected from the twenty-three constituting the fertilization trial, yielding a 40-tree sample of felled trees. One felled sample tree had D equal to the quadratic mean diameter (QMD) for the pooled plot data in a given experimental unit, and one felled sample tree had D equal to the 90th percentile of the diameter distribution for the pooled plot data in a given experimental unit (DBH90). The latter tree was selected to represent a combination of the top height or site tree component of the stand and future crop trees at rotation age. All felled sample trees were free from visually obvious damage (i.e., bear or porcupine scarring, broken top, fork), and were located in an area with local

Data Collection
Permanent plots were established within each experimental unit in the 2009-2010 dormant season (just before fertilization), including three plots in the control experimental units and three or six plots in the fertilized experimental units. Plots were circular and covered 0.04 ha (11.35 m radius). All trees with diameter at breast height (D) > 5.0 cm within the plots were numbered and measured for D (nearest 0.25 cm), total height (H; nearest 0.03 m), and height-to-crown base (HCB; height to the lowest live branch; nearest 0.03 m). Any visible damage or deformity was recorded along with the measurements. Trees within the plots were remeasured during the 2016-2017 dormant season, providing growth responses for the first seven-year period following application. During the 2016-2017 dormant season, two trees from each experimental unit were felled in each of ten stands randomly selected from the twenty-three constituting the fertilization trial, yielding a 40-tree sample of felled trees. One felled sample tree had D equal to the quadratic mean diameter (QMD) for the pooled plot data in a given experimental unit, and one felled sample tree had D equal to the 90th percentile of the diameter distribution for the pooled plot data in a given experimental unit (DBH 90 ). The latter tree was selected to represent a combination of the top height or site tree component of the stand and future crop trees at rotation age. All felled sample trees were free from visually obvious damage (i.e., bear or porcupine scarring, broken top, fork), and were located in an area with local stand structure (e.g., density) consistent with that on the permanent plots in the same experimental unit. Selected trees were visible from permanent plot boundaries (i.e., within approximately 15 m). Felled sample tree measurements were consistent with the standing tree measurements in the permanent plots, so included D (nearest 0.25 cm), H (nearest 0.03 m), and HCB (nearest 0.03 m) ( Table 1). Live crown length (CL; nearest 0.03 m) was calculated as the distance between tree tip and crown base (H-HCB), and crown ratio (CR; %) was calculated as the ratio of CL to H, expressed as a percentage (i.e., 100 × (CL/H)). Height to point-of-insertion (BH; nearest 0.01 m) and basal diameter (BD; nearest 0.1 mm) of every live branch from the base of the live crown to the tree tip were measured. We defined basal branch diameter as the diameter at a distance from bole approximately equal to one branch diameter to avoidbasal bulges that vary in length proportional to branch diameter. Basal diameter of some branches that were damaged from tree felling were estimated based on the approximate taper of intact branches. Each tree crown was divided vertically into thirds of equal length. Branches were numbered consecutively in each crown third, starting with 1 at the base of the third. Three branches were collected randomly from each of the crown thirds, using a random number generator to identify each candidate branch. The sampling protocol required that two branches have a basal diameter >15 mm and one branch a basal diameter >5 mm but <15 mm. Additionally, the total branch length (BL; nearest 1 cm) and distance from tree bole to first live foliage (LLF; nearest 1 cm) on the branch were recorded. Branches were transported to the lab, foliage was separated from branchwood by age class (maximum of six annual cohorts), and each age class was placed in a separate container for drying. After drying for at least 48 hours at 70 • C, each foliage age class was weighed (nearest 0.1 g), providing an oven-dry foliage mass for each sample branch by needle age class. To assess the seven-year responses of foliage production to fertilization, total branch foliage mass (BFM) was expressed as the summed mass of all age classes ( Figure 2). Branch height ranged from 2.88 to 24.42 m, basal diameter ranged from 5.0 to 45.6 mm, total branch length ranged from 18 to 430 cm, distance from bole to first live foliage on a branch ranged from 0 to 175 cm, and foliage mass ranged from 3.2 to 750.2 g ( Table 2).

Branch-Level Foliage Mass
Based on model forms used in previous studies, numerous linear, log-transformed [55], and weighted and unweighted nonlinear [11,20,25] models were tested to develop the following form of a branch-level equation for predicting total foliage mass (g or kg) from branch level variables: where relHACB was the relative height above crown base (1.1 − (DINC/CL)), I Fert was an indicator variable for branches on fertilized trees (1 if fertilized; 0 otherwise), BFM was total foliage mass, BD was basal branch diameter, LLF was distance from tree bole to closest live foliage on the branch, BL was total branch length, DINC was absolute depth into crown (tree height − branch height), and RDINC was relative depth into crown ([tree height − branch height]/crown length). A value of 1.1 was used as a surrogate for maximum relHACB to ensure a nonzero value for the lowest live branch (DINC = CL), thereby avoiding computational problems in the log-transformed and nonlinear regressions [55]. Model errors were assumed additive, random, normal, and independent with variances proportional to a power of branch basal diameter (BD) to be determined from the data. To account for site/stand effects previously observed in the relationship described by Equation (1) (e.g., [11]), stand-level random effects were included in preliminary modeling and assessed by trial and error [12,13,25,27]. Preliminary models also tested for differences in foliage mass between treatments by including an indicator variable for fertilized tree branches (i.e., I Fert ); however, no significant differences were found. Models were fitted in R using the base lm function for linear models and using the lme or nlme functions within the nlme package for linear and nonlinear mixed effects models [56,57]. All model parameter estimates were tested for significant difference from zero at α = 0.05. Final model selection was based on distribution of the residuals, biological relationships of the various predictors, and the following goodness-of-fit criteria considering the alternative weights (BD −m ; m = 0, 1, . . . , 4): likelihood ratio tests, Akaike's Information Criterion (AIC; [58]), generalized R g 2 [59], and unweighted and weighted root mean square error (RMSE and wRMSE, respectively). The latter three criteria were defined as follows: where Y i was the measured foliage mass (g) of a given branch i (i = 1, 2, . . . , N; N = 267);Ŷ i was the model predicted foliage mass for branch i; Y i was the average mass of all foliage sample branches; and w i were the normalized weights. Normalized weights were estimated as the weight for each observation estimated from the variance function divided by the sum of weights for all observations. For example, the normalized weights for a power variance function using branch basal diameter (BD) as the covariate was computed as: was the weight estimated for branch i from the power variance function of basal diameter (BD i ); and t was the power variance coefficient optimized from the data.

Total Crown Foliage Mass
Using the best performing branch-level foliage mass equation, foliage mass was predicted for all live branches on the 40 felled sampled trees and summed for an estimate of total crown foliage mass (TFM; kg). Equation (5) was developed for estimating total crown foliage mass from tree-level predictors by testing various weighted and unweighted linear and nonlinear models. Model forms were based on previous studies [11][12][13]60,61]. To test the hypothesis that, seven years after application, fertilization increased total crown foliage mass on trees with otherwise identical diameter, height, and crown dimensions, an indicator variable (i.e., I Fert ) was included for fertilized trees. The general formulation of the model with all potential predictor variables was as follows: where B was tree basal area (m 2 ); R was crown ratio above breast height (CL/(H − 1.37)); HMC was the height to the middle of the crown ((m); HCB + (CL/2)); RHMC was relative height to middle crown (HMC/H); I Fert was an indicator variable for fertilized trees (1 if fertilized; 0 otherwise); and D, H, CL, and CR were as described above. The parameters on which random stand effects or tree size (i.e., QMD or DBH 90 ) effects were tested were selected by an iterative process between residual plots and biological expectation. Models were fit in R using the base lm function for linear models and using the lme or nlme functions within the nlme package for linear and nonlinear mixed effects models [56,57]. All model parameter estimates were considered significantly different from zero if p ≤ 0.05. Significance of the indicator variable (I Fert ) would suggest a significant fertilization effect assuming least squares parameter estimates are minimum variance and unbiased [62,63].

Vertical Distribution of Foliage
The vertical distribution of foliage mass was estimated empirically for each of the 40 fully measured felled trees by applying the selected branch-level foliage mass equation to each live branch on each felled tree. All foliage on a given branch was assigned to the height of branch attachment. Several alternative probability density functions were fitted in two different ways to the resulting distribution of foliage between the tip of each tree to the live crown base (HCB). In the first approach branch foliage mass estimates were assigned to the height of branch insertion into the tree bole and to the corresponding depth into live crown (DINC). In the second approach, branch foliage mass estimates were grouped into DINC bins of equal length. The probability density functions (PDFs) fitted to the data included the following: (1) a right-truncated Weibull distribution (RTW); (2) an untruncated Weibull distribution (UTW); (3) a Johnson's S B distribution (JSB); and (4) a beta distribution (BET) (see Equations (6)-(8) below). The minima of the distributions were assumed zero (i.e., tree tip), and the maximum values (Johnson and beta) or truncation values (Weibull) were assumed equal to HCB. Preliminary analyses of the ability of different combinations of PDF and resolution of vertical binning indicated that 20 bins were sufficient if not preferable. Summing foliage mass within a fixed segment of crown length reduced the noise introduced by inter-whorl branches and variations in branch size within annual shoots of the main stem [12,13,25,64].
The foliage height assumption described for estimating empirical distributions does not account for the branch angle of origin [11] or branch angle of termination. Significant curvature in the primary branch axis, particularly in older branches, often causes these two angles to differ substantially. In addition, young secondary branches off the primary branch are arranged in generally circular distribution around the branch axis and likewise for higher order branches, all of which can also be somewhat pendant. In short, it would be very difficult to correct simultaneously for branch angle of origin and branch angle of termination to improve the accuracy of empirical estimates of foliage height; similarly, correcting for the variation in spatial arrangement of higher order branches would be even more difficult. Maguire and Bennett [11] and Xu and Harrington [20] argued that adjusting for branch angle would offer only a slight benefit because estimating the effects of branch angle, branch curvature, and the orientation of higher order branches simultaneously would be very complicated. The net improvement in accuracy is likely to remain low relative to assuming that all foliage on a given branch is held at the level of branch attachment to the main stem. Any bias from this assumption would most likely have the effect of underestimating the height of some foliage, with bias increasing with height on the stem height and perhaps reversing in larger trees near the live crown base.
Procedures for obtaining maximum likelihood estimates were as described in Weiskittel et al. [12] and Nelson et al. [27], using an expectation/maximization (EM) algorithm modified from Robinson [65], obtaining initial values for the algorithm using moment-based estimators. The forms of the PDFs fitted to foliage distribution were as follows: Forests 2020, 11, 511 10 of 27 where X represented absolute or relative depth into crown (DINC or RDINC); f 3 (X), f 4 (X), and f 5 (X) were relative density of foliage mass per m of crown length for RTW and per unit relative crown length for the JSB and BET distributions; β, η and Ψ were the Weibull shape, scale, and truncation parameters (Equation (6)); τ and ω were the two JSB shape parameters (Equation (7)); a and b were the two beta shape parameters (Equation (8)); and Γ(x) was the gamma function. Performance of both DINC and RDINC were compared as alternative variables for representing vertical position of foliage in preliminary fitting of the RTW and JSB distributions. Scaling to relative depth into crown (RDINC) facilitated comparison across varying crown lengths [12,13,25,27], and fitting of BET, whose domain is [0,1]. However, the absolute scale performed better and was more biologically interpretable, so was retained in final model selection (with the exception of BET). After predicting foliage mass (g) from each PDF fitted to each tree, RMSE and mean absolute bias (MAB) were computed as the criteria for selecting the best-fitting PDF. After identifying the best-fitting PDF, equations were developed to predict maximum likelihood estimates of the tree-level PDF parameters from tree-level variables. To test for a fertilization effect on vertical distribution of foliage (seven years after application), a fertilization indicator variable was added to the models. A series of linear and nonlinear fixed-effects models and linear and nonlinear mixed-effects models were evaluated starting with model forms from previous studies [11][12][13]25,27]. Inclusion of random effects and procedures for final model selection followed those for the branch-level and total crown foliage mass equations. All distributions and parameter prediction models were fitted in R [57]. Maximum likelihood estimates of distribution parameters and estimated distributions of foliage density and foliage mass were computed with R functions developed by John Kershaw (University of New Brunswick, Fredericton; [66]) and previously applied by Weiskittel et al. [12]. Tree-level parameter prediction models were fitted using the base lm function for linear models and using the lme or nlme functions within the nlme package for linear and nonlinear mixed effects models [56,57]. All model parameter estimates were tested for statistical significance at α = 0.05. Models for estimating each parameter of the fitted distributions were selected and fitted individually. However, because the set of parameters for a given PDF are contemporaneously correlated, estimates from individual equation fits are inconsistent and inefficient. Therefore, after final models were selected for predicting each of the parameters for each PDF, the models were refitted as a system-of-equations using iterative three-stage least squares (3SLS; [67]). In the presence of contemporaneously correlations among parameter estimates for a PDF, three-stage least squares lead to consistent and asymptotically more efficient estimates. The 3SLS systems were fitted using the systemfit package in R [57,68], following procedures from Schmidt [69] to obtain 3SLS estimates. Once the three-stage least squares estimates were obtained, performance was evaluated by comparing the RMSE from the system fit (3SLS) to RMSE from ordinary least squares (OLS). Goodness of fit for the system was evaluated by McElroy's R M 2 [70]. Statistical significance (i.e., α = 0.05) of the fertilization indicator variable in one or more of the equations for the system would indicate a significant difference in one or more PDF parameter estimates between fertilized and control trees, suggesting a significant fertilization effect on relative vertical distribution of foliage.

Branch-Level Foliage Mass
The best performing model was based on the form presented by Garber and Maguire [25] for three conifer species in central Oregon and by Maguire and Bennett [11] for coastal Douglas-fir. This model form had the highest log likelihood score and R g 2 , and lowest AIC and RMSE relative to alternative models. The final model (Equation (9)) included a random stand effect on the first relative depth into crown (RDINC) term and was weighted using a power variance function of BD to correct for heteroscedasticity. The final model form was: where FM ijk was the observed foliage mass (g) for branch i (i = 1, 2, . . . , n j ) on tree j (j = 1, 2, 3, 4) in stand k (k = 1, 2, . . . , 10); BD ijk was basal diameter (mm) for the sample branch; RDINC ijk was relative depth into the tree crown of the sample branch; α 0 , α 1 , α 2 , and α 3 were parameters to be estimated from the data; δ 3,k was a random effect of the kth stand; and ε 1,ijk was the random error term for the ith branch on the jth tree in the kth stand. Random effects δ 3,k and ε 1,ijk were assumed to have a multivariate normal distribution specified as follows: where δ 3 was the 10x1 vector of random stand effects; ε 1 was the 267x1 vector of random branch errors; 0 δ3 was the mean vector of random stand effects; 0 ε1 was the mean vector of random errors; σ δ3 2 I was the variance-covariance matrix for the random stand effects; and σ ε1 2 Σ was the block-diagonal variance-covariance matrix for the random errors, with the block diagonals allowing for potential correlations between branch mass observations within a tree. Zero covariance was assumed between random stand effects and random branch errors. Branch diameter and relative depth into crown imposed highly significant fixed effects on branch foliage mass, and the random stand effect was also significant (Table 3). Fertilization was not detected to have any marginal effect on the amount of foliage held on a branch of given basal diameter at a given relative depth into crown, apparently because any increase in foliage mass was proportional to branch diameter. Table 3. Parameters, estimates, standard errors, and p-values for the final branch-level foliage mass model (Equation (9)). The final model was applied to estimate total branch foliage mass from measured RDINC and BD of all live branches on the 40 felled sample trees. The model indicated that for a given BD, branch foliage mass peaked approximately halfway between crown base and tree tip (Figure 3).

Parameter
The final model was applied to estimate total branch foliage mass from measured RDINC and BD of all live branches on the 40 felled sample trees. The model indicated that for a given BD, branch foliage mass peaked approximately halfway between crown base and tree tip (Figure 3).  (9)) for all live branches measured from felled sampled trees. RDINC = 1 at crown base and RDINC = 0 at tree tip.

Total Crown Foliage Mass
The best performing model for predicting total tree foliage mass was based on forms tested by Maguire and Bennett [11] and Williams et al. [13]. The final model (Equation (10)) was a nonlinear mixed-effects model that included a nested random effect of experimental unit within a given stand. The final model took the following form: where TFMjkl was the total foliage mass (kg) for tree j in experimental unit l in stand k (j = 1, 2; l = 1, 2; k = 1, 2, ..., 10); Bjkl was tree basal area (m 2 ); Rjkl was crown ratio above breast height (CL/(H − 1.37)); HMCjkl was height to middle of crown length (m); Hjkl was tree height (m); IFert was an fertilization indicator variable (1 if fertilized; 0 otherwise); β0, β1, β2, β3, and β4 were parameters to be estimated from the data; δ1,kl was the random effect of the l th experimental unit in the k th stand; and ε2,jkl was the random tree error for the j th tree in the l th experimental unit in the k th stand. Random effects δ1,kl and ε2,jkl were assumed to have a multivariate normal distribution specified as follows:  (9)) for all live branches measured from felled sampled trees. RDINC = 1 at crown base and RDINC = 0 at tree tip.

Total Crown Foliage Mass
The best performing model for predicting total tree foliage mass was based on forms tested by Maguire and Bennett [11] and Williams et al. [13]. The final model (Equation (10)) was a nonlinear mixed-effects model that included a nested random effect of experimental unit within a given stand. The final model took the following form: where TFM jkl was the total foliage mass (kg) for tree j in experimental unit l in stand k (j = 1, 2; l = 1, 2; k = 1, 2, . . . , 10); B jkl was tree basal area (m 2 ); R jkl was crown ratio above breast height (CL/(H − 1.37)); HMC jkl was height to middle of crown length (m); H jkl was tree height (m); I Fert was an fertilization indicator variable (1 if fertilized; 0 otherwise); β 0 , β 1 , β 2 , β 3 , and β 4 were parameters to be estimated from the data; δ 1,kl was the random effect of the lth experimental unit in the kth stand; and ε 2,jkl was the random tree error for the jth tree in the lth experimental unit in the kth stand. Random effects δ 1,kl and ε 2,jkl were assumed to have a multivariate normal distribution specified as follows: where δ 1 was a 20x1 vector of random effects for the lth experimental unit in the kth stand; ε 2 was a 39x1 vector of random tree errors; 0 δ1 was the mean vector of random experimental unit effects; 0 ε2 was the mean vector of random tree errors; σ δ1 2 Λ was the variance-covariance matrix for the random experimental unit effects; and σ ε2 2 Σ was the block-diagonal variance-covariance matrix for random errors, with the block diagonals allowing for potential correlations between foliage mass of trees within an experimental unit. Zero covariance was assumed between experimental unit random effects and tree random errors. Equation (10) had the highest R g 2 and log likelihood score and lowest AIC and RMSE of all candidate models. One tree was excluded from model fitting because it imposed unusually strong influence on model behavior and appeared as an unrepresentative outlier, perhaps due to an asymmetric crown or measurement errors. All tree-level covariates imposed significant effects on total tree foliage mass (all p < 0.05; Table 4). Table 4. Parameters, estimates, standard error, and p-values for the final tree-level foliage mass model (Equation (10)).

Vertical Distribution of Foliage
Based on the root mean square error (RMSE) and mean absolute bias (MAB), performance was similar between the truncated and untruncated Weibull distributions, regardless of whether parameters were estimated from individual branch data or data aggregated into vertical segments of constant length. In contrast, the beta and Johnson distributions performed better with data aggregated by live crown segment, but the higher variances on parameter estimates from all distributions using aggregated data suggested that aggregation could have masked patterns in finer-scale tree-to-tree variability that have been observed in other studies (e.g., [26]). Overall, the best performing PDF was a right-truncated Weibull fitted to unaggregated data, as measured by RMSE, compared to beta and Johnson distributions fitted to the same data (Table 5). Therefore, the right-truncated Weibull was selected as the final model for characterizing vertical distribution of foliage on individual Douglas-fir trees. To test for potential differences in vertical foliage distribution between fertilized and unfertilized Douglas-fir trees in the SLC study, equations were developed for predicting tree-level parameter estimates obtained from the right-truncated Weibull distributions fitted to the data. Tree-level dimensions were screened in both a transformed and untransformed state, and an indicator variable for fertilization was also tested in various components. Preliminary models suggested the most influential variables included diameter at breast height (D), a surrogate for stem taper and/or relative crown size (e.g., D/H), total crown foliage mass (TFM), and a direct measure of crown size (e.g., crown length, height-to-crown base, crown ratio). Because total crown foliage mass is predicted (i.e., Equation (10)), and subsequently used as a predictor, the possibility of correlation among errors has to be accounted for. This becomes particularly important when applying this system to an independent tree or tree-list. Therefore, predicted total crown foliage mass (TFM) from the final model (Equation (10)) was used instead of observed TFM, in which case, TFM can be viewed as a transformation of exogenous tree-level variables. The final equations took the following form: whereη j was the expectation/maximization (EM) prediction of the scale parameter in Equation (6) and β j was the EM-prediction of the shape parameter in Equation (6) for tree j;TFM j was the predicted total crown foliage mass (kg) from Equation (10); ε 3 was the residual error for prediction of estimatedη j ; ε 4 was the residual error around the predictions of estimatedβ j ; and all other variables were described previously. Prior to refitting as a system-of-equations, the residual plots of individual equations were examined and no evidence of increase or decreasing variance across the range of either right-truncated Weibull distribution parameters was found. The iterative three-stage simultaneous least squares (3SLS) estimates of all tree-level parameters were significantly different from zero at α = 0.05 ( Table 6). The system explained 91.6% and 41.7% of the original variation in the scale and shape parameter estimates, respectively, and the McElroy R M 2 [70] indicated that the system accounted for 85.6% of the combined variation in estimates of the two right-truncated Weibull parameters. The significance of the parameter estimate on the fertilization indicator (κ 3 : p-value = 0.024) suggested a significant effect of fertilization on the vertical distribution of foliage. Table 6. Parameters, estimates, standard errors, p-values, and fit statistics for the models predicting EM-estimates (see text) of tree-level, right-truncated Weibull parameters η (Equation (11)) and β (Equation (12)) fitted by iterative three-stage least squares (3SLS). To compare the foliage distributions between the control and fertilized treatments, total crown foliage mass and Weibull distribution parameters were estimated based on the average size of the 40 felled sample trees in the study. Relative foliage mass density was then estimated from the right-truncated Weibull PDF (Equation (6)) and plotted on relative depth into crown (RDINC; Figure 4). The distribution of relative foliage mass density peaked at approximately mid-crown, and the peak location was nearly identical between the control and fertilized treatments (0.55 and 0.56, respectively). Fertilization induced a slight decrease in relative foliage mass density in the top compared to unfertilized control trees. The final model for total crown foliage (Equation (10)) predicted approximately 10.35% more foliar mass (23.30 kg on fertilized trees versus 21.11 kg on unfertilized controls). This additional foliage mass apparently accumulated at mid-crown. The fertilization effect on the shape (skew) parameter β appeared to drive the decreases in relative foliage mass density observed near tree tip and near crown base and the slight downward shift in the mode of the foliage distribution. The increase in total crown foliage mass induced by fertilization drove the relative foliage mass density increases observed at mid-crown. These results confirm that nitrogen fertilization increases total crown mass and shifts the relative vertical distribution of foliage mass on individual Douglas-fir trees. To examine further the vertical distribution of foliage, the cumulative form of the right-truncated Weibull distribution (CDF) was estimated by dividing the crown into 100 segments of depth into crown, determining the cumulative proportion of foliage mass at a given depth, and multiplying the respective cumulative proportion by the total crown foliage mass. As with relative foliage mass density per unit vertical segment (Figure 4), a direct treatment effect was evident in the cumulative foliage mass for the fertilized and unfertilized tree of average size ( Figure 5). Based on the CDF, foliage mass first increased with increasing RDINC just above mid-crown and continued to increase to within approximately 15% of crown base. Very slight decreases apparent in the top-third of the crown were probably not statistically or biologically significant. To examine further the vertical distribution of foliage, the cumulative form of the right-truncated Weibull distribution (CDF) was estimated by dividing the crown into 100 segments of depth into crown, determining the cumulative proportion of foliage mass at a given depth, and multiplying the respective cumulative proportion by the total crown foliage mass. As with relative foliage mass density per unit vertical segment (Figure 4), a direct treatment effect was evident in the cumulative foliage mass for the fertilized and unfertilized tree of average size ( Figure 5). Based on the CDF, foliage mass first increased with increasing RDINC just above mid-crown and continued to increase to within approximately 15% of crown base. Very slight decreases apparent in the top-third of the crown were probably not statistically or biologically significant.

Branch-Level Foliage Mass
The first objective of this study was addressed by developing a robust model for estimating branch foliage mass in the target population and testing for fertilization effects on total branch foliage mass. The weighted, nonlinear mixed-effects model utilizing branch-and tree-level variables proved most effective among alternatives for predicting foliage mass at the branch-level. This model form has previously been found the best among alternatives explored for Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) [11], as well as for other conifer species (e.g., grand fir (Abies grandis (Dougl. ex D. Don) Lindle.), lodgepole pine (Pinus contorta Dougl. ex Loud.), and ponderosa pine (Pinus ponderosa Dougl. ex Laws.) [25]; balsam fir (Abies balsamea (L.) Mill.), northern white-cedar (Thuja occidentalis (L.)), eastern hemlock (Tsuga canadenis (L.) Carr.), eastern white pine (Pinus strobus (L.)), and red spruce (Picea rubens (Sarg.)) [12]). As with this and previous models, branch foliage mass increased with increasing branch diameter; however, as branches near crown base start losing foliage due to shading, relative depth into crown becomes an important variable for modeling the decline in foliage mass as branches of large diameter become relegated to a lower canopy position near crown base. This decline in foliage mass is slower in more shade-intolerant species [13,25,71], but photosynthetically active radiation (PAR) eventually falls below the light compensation point and branches can no longer survive. As a result, branch foliage mass increases with branch diameter while

Branch-Level Foliage Mass
The first objective of this study was addressed by developing a robust model for estimating branch foliage mass in the target population and testing for fertilization effects on total branch foliage mass. The weighted, nonlinear mixed-effects model utilizing branch-and tree-level variables proved most effective among alternatives for predicting foliage mass at the branch-level. This model form has previously been found the best among alternatives explored for Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) [11], as well as for other conifer species (e.g., grand fir (Abies grandis (Dougl. ex D. Don) Lindle.), lodgepole pine (Pinus contorta Dougl. ex Loud.), and ponderosa pine (Pinus ponderosa Dougl. ex Laws.) [25]; balsam fir (Abies balsamea (L.) Mill.), northern white-cedar (Thuja occidentalis (L.)), eastern hemlock (Tsuga canadenis (L.) Carr.), eastern white pine (Pinus strobus (L.)), and red spruce (Picea rubens (Sarg.)) [12]). As with this and previous models, branch foliage mass increased with increasing branch diameter; however, as branches near crown base start losing foliage due to shading, relative depth into crown becomes an important variable for modeling the decline in foliage mass as branches of large diameter become relegated to a lower canopy position near crown base. This decline in foliage mass is slower in more shade-intolerant species [13,25,71], but photosynthetically active radiation (PAR) eventually falls below the light compensation point and branches can no longer survive. As a result, branch foliage mass increases with branch diameter while simultaneously decreasing with depth into crown. The net effect is an increase in branch foliage mass with increasing distance from tree tip until reaching a peak mid-crown, below which branch foliage mass decreases toward crown base. This pattern has been consistently observed in other studies and demonstrates the importance of including both branch diameter and some measure of position in crown as predictors [11,12,18,25,31,55,72,73]. The random stand effect interacted with relative depth into crown but was consistent with observed differences among sites in the foliage mass supported by branches of a given diameter and crown depth (e.g., [11]). Although the stands within this study were selected in part based on their similar structure and silvicultural history, the degree to which small differences in age, site quality, stand density and other potential factors affect the foliage mass carried by branches of similar size and crown position is poorly understood.
No significant fertilization effects on foliage mass of branches of a given size were apparent based on either the final model or extensive preliminary model fitting. The second phase of the first objective therefore resulted in failure to reject the implied null hypothesis that branches on fertilized trees in fertilized stands do not hold greater amounts of foliage biomass after accounting for their size and crown position. Fertilization has been demonstrated to increase branch size [9], probably to an extent that is commensurate with any increase in foliage mass and increase in tree size, and to have no effect on number of whorl or inter-whorl branches per unit length of main stem [74]; hence, the model accounted for the increase in foliage mass of fertilized trees predominantly through the increase in branch diameters. The height within the crown where fertilization increased foliage mass was evident by plotting branch foliage mass on relative depth into crown (Figure 2). It is unclear from the analysis presented whether increases in branch size were correlated with foliated branch length, branch width, or branch length, any of which may indicate finer scale mechanisms explaining the greater foliage mass on branches with larger diameters. Additional branch-level terms within models have the potential to improve model performance, and subsequently account for these differences; however, including some of these variables may also exacerbate potential adverse consequences of multicollinearity [55,63]. Increases in branch diameter were most likely correlated with increases in branch length, and depth into crown should account for the decline in foliated branch length with increasing shading toward the crown interior and at crown base. Because stand density and water relations typically affect the latter, the significant stand-level random effect presumably accounted for some of these relationships.

Total Crown Foliage Mass
The nonlinear mixed-effects model for estimating total tree foliage mass on both fertilized and unfertilized trees addressed the second objective of this analysis. At the tree level, a significant effect of fertilization proved significant in a final model that performed best among many alternatives for estimating the total foliage mass of the Douglas-fir trees in the study.
Fertilization increased total tree foliage mass through the increases in foliage production and branch size observed in previous studies [8,9]. The model incorporated various transformations of tree-level variables previously used by Maguire and Bennett [11] and Williams et al. [13]. Among the strongest predictors were a combination of tree basal area (B), crown ratio above breast height (R), height to mid-crown (HMC), and the ratio of tree basal area to relative height to mid-crown. The product of crown ratio and tree basal area can be viewed as a surrogate for basal area at crown base, with CR serving as a crude taper model for basal area (B). The ratio of tree basal area to relative height to mid-crown served as a similar measure, and the combination agreed with previous observations that crown mass increased with crown length and decreased with increasing height [11][12][13]61,75]. Crown size for a given DBH is highly variable among stands of varying density and top height, suggesting that these supplemental predictors were highly beneficial [11]. Previous research has shown that diameter at crown base [76], sapwood area at breast height [77], sapwood area at crown base [78], and gross crown dimensions [79] can all serve as strong predictors of total crown foliage mass [11]. The combination of transformed DBH and crown variables used in the final tree-level model in this study likely represented the functional effect of these other direct measures, and perhaps with less measurement error.
The second objective of this study implied testing the null hypothesis that fertilization did not affect the amount of foliage held by individual trees, after accounting for diameter, height, and crown size of the tree. The final tree foliage model firmly rejected this null hypothesis. For trees of equal size initial size, fertilization resulted in approximately 2.2 kg (~10%) more foliar mass, seven years following application. Growth responses to fertilization treatment further indicated that indirect growth responses were driven by positive feedbacks through associated increase in foliage mass. Although crown lengths on fertilized tree were slightly longer than those on control trees in this study (average crown length 13.34 and 13.64 m, respectively), this was offset by slightly taller height-to-crown base (average height to crown base: 7.14 and 7.41 m, respectively). The fertilization effect through the height-to-crown base term may be attributable to differences in diameter at crown base (e.g., [75]), as the largest increases in inside-bark stem diameter as a result of fertilization have been shown to be at or near this point [52]. Several studies have concluded that short-term (2-5-year) growth response to fertilization was attributable to increased photosynthetic efficiency [9] and/or increased water use efficiency [44]. In contrast, the primary factor in the sustained long-term volume growth responses to fertilization was greater foliar area [8,9,35,44,80]. This increase in foliar area has also been found contingent upon initial foliage amount and light limitations at time of fertilization [8,9,81]. Similar conditions prevailed in this study, where the largest increases in foliage mass were observed mid-crown, where favorable light conditions may still exist among the largest, most competitive branches. The individual-tree foliage and corresponding growth responses to fertilization summed to significant growth responses at the stand level as well, due to increases in stand-level foliage mass and faster accumulation of growing stock [52].
The significant interaction between fertilization and height-to-crown base (HCB) may have been of particular importance in this study given that the positive effect of fertilization on foliage mass increased with increasing initial height-to-crown base. This positive interaction suggested a greater response by increasingly dominant trees, a seven-year direct response to fertilization not captured by tree growth responses in diameter, height, or crown length. The growth response linked to crown base was also consistent with the influence of fertilization on stem form [52] as mentioned previously. The increase in total crown foliage mass following fertilization aligns with observations from previous studies, as fertilization has been demonstrated to increase production of foliage, increase needle and branch size, increase longevity of lower branches, and subsequently increase live crown length [8,9].

Vertical Distribution of Foliage
The third objective of this analysis was addressed by first estimating the empirical distribution of foliage mass over relative depth into crown using the branch-level foliage equations. Different probability density functions were then fitted to these empirical distributions for each tree, and prediction equations for the parameters estimated in the selected probability density function (PDF) were developed and tested for any fertilization effects. As has been found in previous studies, the Weibull distribution was the best performing distribution for characterizing the vertical distribution of foliage [12,13,[15][16][17][18][19][20][21]. Right-truncation at crown base provided a more biologically appropriate model than the unmodified Weibull distribution with a domain of [0,∞). That fact that at least some live foliage remains at crown base probably gave an edge to the right-truncated Weibull distribution over the beta distribution. In the final model vertical distribution of foliage depended on DBH, total crown foliage mass, a surrogate for stem form (i.e., H/D), and crown-related variables. As with other studies, crown size and foliar mass have consistently emerged as strong predictors of the parameters controlling the shape of the distribution [11][12][13]. The combination of tree basal area and the relative position of mid-crown was an effective index of foliage distribution.
Meeting the third objective also required testing the null hypothesis that fertilization had no effect on the vertical distribution of foliage mass. In fact, fertilization imposed a significant effect on vertical foliage distribution. The distributional patterns in response to fertilization revealed that the largest increases in foliage mass occurred near mid-crown, where large established branches were still receiving significant amounts of light but could benefit from more. The decrease in foliage mass compared to the control near the upper portion of the stem (Figure 3), could potentially be the result of increased height growth in response to fertilization. Also, the pattern of diminished foliage mass near crown base is driven by increased shading in response to foliage production at higher levels in the crown as the tree grows in height. The acceleration in foliage production stimulated by fertilization almost certainly exacerbated this foliage loss. Also, differences in foliage distribution between control and fertilized trees entailed not only increases in total crown foliage mass, but also increases in the shape (skew) parameter β, suggesting that the mode of the distribution is shifted slightly downward with a longer tail toward the tip of the tree. There was no significant effect of fertilization on the scale (kurtosis) parameter. Based on the distribution of branch-level foliage (Figure 2), the relative height of the peak in the distribution is similar between treatments, but the differences in maxima indicate the increase in total crown foliage mass after fertilization. These patterns are best observed by examining the vertical distribution of foliage on an absolute scale for the average size of the twenty felled fertilized trees and the average size of the twenty felled unfertilized trees ( Figure 6; Table 1).
Forests 2020, 11, x FOR PEER REVIEW 20 of 27 compared to the control near the upper portion of the stem (Figure 3), could potentially be the result of increased height growth in response to fertilization. Also, the pattern of diminished foliage mass near crown base is driven by increased shading in response to foliage production at higher levels in the crown as the tree grows in height. The acceleration in foliage production stimulated by fertilization almost certainly exacerbated this foliage loss. Also, differences in foliage distribution between control and fertilized trees entailed not only increases in total crown foliage mass, but also increases in the shape (skew) parameter β, suggesting that the mode of the distribution is shifted slightly downward with a longer tail toward the tip of the tree. There was no significant effect of fertilization on the scale (kurtosis) parameter. Based on the distribution of branch-level foliage (Figure 2), the relative height of the peak in the distribution is similar between treatments, but the differences in maxima indicate the increase in total crown foliage mass after fertilization. These patterns are best observed by examining the vertical distribution of foliage on an absolute scale for the average size of the twenty felled fertilized trees and the average size of the twenty felled unfertilized trees ( Figure 6; Table 1).  Table 1). RDINC = 1 at crown base and RDINC = 0 at tree tip.
Water availability, vapor pressure deficit, and other crown microclimatic factors may be contributing to the responses observed in foliage distributions on these fertilized trees. Gravitational water potentials are more negative [82] and vapor pressure deficits are likely to be higher in the upper canopy, which may limit foliage production. A higher severity of this effect is expected for the more Figure 6. Effect of nitrogen fertilization on vertical foliage distribution over depth into crown (DINC; m) as modeled with a right-truncated Weibull distribution (control = solid line; fertilized = dashed line). Vertical foliage distribution is represented as absolute foliage mass (kg). Distributions were standardized to the average tree size by treatment in the study (See Table 1). RDINC = 1 at crown base and RDINC = 0 at tree tip.
Water availability, vapor pressure deficit, and other crown microclimatic factors may be contributing to the responses observed in foliage distributions on these fertilized trees. Gravitational water potentials are more negative [82] and vapor pressure deficits are likely to be higher in the upper canopy, which may limit foliage production. A higher severity of this effect is expected for the more dominant trees and may partially explain the lower foliage mass per unit relative crown length in the upper crown of fertilized trees, despite more than adequate light availability. The largest increases in foliage mass observed at mid-crown may represent a balance between minimizing water loss and maximizing light capture, contingent upon the degree of crown shading [13]. Factors forcing this balance may be particularly strong in the middle and eastern parts of the Oregon Coast Range where summer water availability is routinely limiting. Determining soil water availability may therefore be critical in determining seasonal and diurnal patterns in stomatal conductance and associated capacity for response to fertilization [35]. Douglas-fir foliage distributions have exhibited an upward shift with increasing crown competition [11]. The results of this study suggest that fertilization may shift the relative distribution of foliage mass in the opposite direction toward the midpoint of the live crown length.
Trees sampled in this study ranged from the middle (i.e., QMD) to upper end (i.e., DBH 90 ) of the diameter distributions, limiting the opportunity to assess effects across the full range of current tree-crown social position. This sampling strategy may limit short-term stand-level inference but does focus on the most likely crop trees at rotation age. In many other studies of foliage distribution, sampling often targets the most vigorous trees in the stand [12]. Factors such as stand age [83], stand density [20], and tree relative height [11], influence vertical foliage distribution [12]. Future studies could include trees from the full mid-rotation range of diameters and crown positions to clarify silvicultural options and strengthen inferences on early stand-level responses to silvicultural treatments. However, the major gains in wood production and value recovery will depend on tree near and above mid-rotation quadratic mean diameter.
The described data will facilitate improvements in models for simulating combination of fertilization with other silvicultural treatments such as thinning. Fertilization effects on the entire stand canopy and internal stand dynamics will advance our understanding of both direct and indirect effects of this silvicultural practice, particularly if better quantification of growing stock accumulation emerges and leads to improvements in estimating indirect effects of increased basal area, cambial surface area, crown surface area, and foliage area or mass. Changes to many measures of crown size influence dynamic factors such as potential height growth, crown recession, and probability of mortality, as well as ecophysiological mechanisms. Crown responses therefore play a central role in refining growth and yield models (e.g., ORGANON [84]). Understanding the implications of shifts in foliage distribution in response to fertilization is just one step toward obtaining more precise and reliable predictions of tree growth and stand development over time.

Conclusions
Light capture and the subsequent process of converting fixed carbon to stem-wood requires accurate quantification of complex canopy structures, including spatial distribution of foliage mass and area [85][86][87][88]. Our understanding of growth efficiency, commonly defined as stem volume growth per unit leaf area, has improved in recent years [71] and can be used to assess the relationship between tree growth, stand structure, and stand productivity [88][89][90]. In even-aged stands with a relatively simple stand structure, growth efficiency decreases with increasing total leaf area [71,88,91]. This measure of tree productivity has the potential to become more useful and insightful if some of its variability within and between stands can be predicted from the quantity and distribution of foliage area on subject trees. Pressler's hypothesis states that stem increment along the stem is proportional to the quantity of foliage above that stem point (as cited in [92]). Given the predicted cumulative distribution of foliage (Figure 4), the cross-sectional increment would be expected to be larger in the lower two-thirds of a given tree crown, and smaller from that point to tree tip. Putney [52] developed a variable exponent taper model based on the same destructively sampled trees in this study, finding that the largest growth response of diameter inside-bark to fertilization occurred at mid-stem, and that a slight decrease in stem diameter growth response occurred near tree tip. The results presented in this study provide evidence for the "cause" of this "effect." Given successful quantification of cumulative foliage distribution predictable from silvicultural treatments and their effects on tree-level variables, a unique opportunity presents itself for refining periodic and cumulative growth predictions by supplementing empirical predictors with more physiologically direct mechanisms. Furthermore, because the foliage analyses presented in this study are directly relevant to photosynthetic capacity and net primary productivity (NPP), an opportunity is open to explore methods of incorporating these responses and relationships into growth and yield models. Hybrid models with empirical and mechanistic elements that combine statistical-and ecophysiological-based approaches offer a potential advantage of increased flexibility and simplicity over purely statistical or process-based models. A working hypothesis is that hybrid models can mechanistically represent silvicultural influences through relatively few, wisely selected physiological principles [93]. These models typically derive stand-level productivity predictions directly from physiological processes to supplement or replace covariates in statistical models (e.g., estimating net photosynthesis to supplement site index; [31,32]). Hybrid modeling has been coined the future of forest growth modeling [94] and has been successfully utilized to reduce mean square error (MSE) of specific output in numerous growth models [7,32,[95][96][97][98]. However, improvements to these models are only achievable through successful identification and quantification of only key mechanistic processes, so that the quantity and quality of data required is operationally feasible. Developing the framework and methodology for quantification of mechanisms is therefore crucial to furthering the development of these models.
In this study, fertilization resulted in increased crown foliar mass of individual trees, concentrated at mid-crown ( Figure 6). The decrease in foliage near the upper crown and near crown base resulted from growth reallocation to the middle third of the crown, and increased shading of the lower crown by the increased foliar mass above. The model presented in this study demonstrated that light availability will be crucial for the maximization of foliage production, particularly in the first 3-4 years following fertilization [9] and suggests that available soil water and concurrent stand density regulation are essential. While not novel, the stand selection criterion in this study that required precommercially thinned stands corroborates the notion that stand density must be optimal prior to fertilization (i.e., well below thresholds where density-induced competition may occur, where fertilized trees do not get adequate light, or where the stand has already reached maximum foliage mass carrying capacity). Prior to fertilization, tree crowns must be maintained in a vigorous and well-formed condition, despite inevitable natural disturbances, to ensure maximal capacity for response. Prime candidates for fertilization should be free as possible from disease and storm damage, exhibit optimal or near optimal pretreatment density, and possess healthy crowns.