A New Model for Size-Dependent Tree Growth in Forests

Tree growth, especially diameter growth of tree stems, is an important issue for understanding the productivity and dynamics of forest stands. Metabolic scaling theory predicted that the 2/3 power of stem diameter at a certain time is a linear function of the 2/3 power of the initial diameter and that the diameter growth rate scales to the 1/3 power of the initial diameter. We tested these predictions of the metabolic scaling theory for 11 Japanese secondary forests at various growth stages. The predictions were not supported by the data, especially in younger stands. Alternatively, we proposed a new theoretical model for stem diameter growth on the basis of six assumptions. All these assumptions were supported by the data. The model produced a nearly linear to curvilinear relationship between the 2/3 power of stem diameters at two different times. It also fitted well to the curvilinear relationship between diameter growth rate and the initial diameter. Our model fitted better than the metabolic scaling theory, suggesting the importance of asymmetric competition among trees, which has not been incorporated in the metabolic scaling theory.


Introduction
It is important to understand and predict tree growth because it is related to the productivity and dynamics of forest stands. Furthermore, tree growth rate is an important demographic parameter, together with mortality and recruitment rates [1]. Among various measurements of tree growth, diameter of tree stem has been most widely monitored because of the easiness and accuracy of the measurement compared to tree height.
There are several theoretical models to express stem diameter growth of individual tree. One of the simplest is a model based on the metabolic scaling theory. Enquist et al. [2] predicted that stem diameter growth rate scales to the 1/3 power of the diameter. This prediction was deduced from the assumption that growth rate of tree biomass is proportional to its gross photosynthetic rate, which in turn is determined by leaf biomass. This assumption means that both small and large trees have the same resource availability per unit leaf biomass.
However, resource availability may depend on tree size. In a closed forest stand, light is a limited resource and tree growth is governed mainly by light conditions [3][4][5]. Availability and competition for light is size-dependent and asymmetric [6][7][8]. Larger trees have an advantage over smaller trees, since the former shade the latter but the latter rarely do so [9]. In some extreme cases, competition is one-sided: larger tree suppress the growth of smaller trees but not vice versa [7,10,11]. Recent studies [6,12] extended the metabolic scaling theory to incorporate asymmetric competition. Their model fitted better to the data of natural forests than that of Enquist et al. [2], suggesting the importance of asymmetric competition.
Here, we propose a new theoretical model for tree diameter growth that implicitly considers asymmetric competition by extending the model of Kikuzawa [10]. Our model is based on a tree size distribution function. Frequency distribution of the tree biomass is often inverse Jshaped with many small trees and few larger trees. Small size differences among trees at the seedling or sapling stage due to genetic factors or slight differences in environmental conditions at the micro site will be amplified by competition among individuals. The initial normal frequency distribution in biomass becomes skewed to the right under asymmetric competition [7,11,13,14]. Such skewed size distribution can be expressed in Y-N curve: Y = N/(AN + B) [15] where Y is the cumulative tree biomass in the stand summed for the largest N trees and N is the rank of a tree ordered from the largest to the smallest tree. This relationship fitted to both even-aged plantation and uneven-aged natural forests [15,16]. Our model is based on this Y-N curve, which was later derived theoretically from a growth model assuming asymmetric competition [17]. We tested assumptions of our model and the fit of our model against tree growth data of multi-species secondary deciduous broadleaf forests in Hokkaido, northern Japan. We compared our model to the model by Enquist et al. [2].

Model of Enquist et al. (1999)
Enquist et al. [2] proposed a model for stem diameter growth on the basis of the metabolic scaling theory. They started from the assumption that whole plant biomass growth rate (dM/dt) is proportional to leaf biomass (L B ) where M is whole plant biomass, t is time, C G is a constant [18,19]. Here, leaf biomass L B is proportional to the 3/4 power of individual tree volume V (= M/ρ; ρ being the density of tree body). On the other hand, stem diameter of a tree is proportional to the 3/8 power of the tree's volume. By using these allometric correlations, Enquist et al. [2] obtained the following equation for stem diameter growth: where C is the proportionality constant. By substitution and integration from t = 0 to t = T, they derived the following equation: The 2/3 power of stem diameter at time T (D T ) can be expressed as a linear function of the 2/3 power of the initial stem diameter (D 0 ), having a positive intercept and a slope of unity. Enquist et al. [2] argued that this equation was held in a secondary tropical forest.

Alternative new model
We developed an alternative tree growth model by extending that of Kikuzawa [10]. This model is based on six assumptions that will be described in the following parts: Assumption 1: Y-N curve fits to the size frequency distribution of trees. Hozumi et al. [15] observed in natural forests and a plantation that the size frequency distribution of trees is inverse J-shaped, with many small trees and few larger trees. In many cases [15,16], such size hierarchy can be expressed by the following equation (assumption 1): where Y is the cumulative biomass in the stand summed from the individual with largest biomass to the N-th largest tree, N is the rank of a tree ordered from the largest individual, and A and B are parameters that are constants but change with time. Later, this equation was derived theoretically from the tree growth model that assumed asymmetrical competition [17]. Therefore, this relationship stands on a firm ground. Under Eq (4), the biomass of an individual tree (M) that comprises the stand is expressed by the following equation [15]: This equation means that individual tree's biomass is expressed as the function of N, the tree's rank. Individual tree biomass of the same rank but different time (t = 0, T) will be expressed by the following equations: [10] suggested that under one-sided competition, individual tree's rank changes little. Assuming that an individual tree's rank does not change between two different times (assumption 2), we can eliminate N from Eqs (5') and (5") and obtain the following relationship between M 0 and M T : Parameters A and B have biological meanings, and there is a relationship between parameters A and B [10]. Parameter A is a reciprocal of stand biomass (Y max ) [15]: Y max was assumed to be expressed as follows where d is stand biomass divided by the height of the highest individual in the stand (H max ). Kira and Shidei [20] define d as dry matter density and stated that d takes a constant value of 1-1.5 kg/m 3 , irrespective of stand height. By generalizing their findings, Kikuzawa [10] found that d is also a function of H max : where a is a parameter and K 1 is a constant. When a = 0, Eq (9) will be consistent with the statement of Kira and Shidei [20]. Substituting Eq (8) by Eq (9) will yield the following power function between stand biomass and the height of the highest individual in the stand (assumption 3): Substituting Eq (7) by Eq (10), we obtain Parameter B is the reciprocal of the maximum individual biomass of the stand (M max ) [15]: Kikuzawa [10] assumed that tree biomass M has an allometric relationship with tree height (H) (assumption 4) [21]: where b is a parameter and K 2 is a constant. Applying this allometric relation to the maximum individual biomass (M max ), we obtain: By applying Eqs (14) to (12), we obtain: From Eqs (11) and (15), the scaling relationship between A and B is obtained (assumption 5) for full-density stands as follows: where K 3 and c are constants. Substitution of Eq (16)) in Eq (6) will yield the following: This equation shows the relationship between the biomass of the same ranked tree at a different time under the condition that Eq (4) is equally held at different times.
West et al. [22] and Enquist et al. [2] showed that the biomass of an individual tree is proportional to the 8/3 power of stem diameter D, M = K 4 D 8/3 (assumption 6). Eq (17) can be rewritten considering this allometric relationship as follows: This equation expresses the relationship between the 2/3 power of stem diameter at t = T to the 2/3 power of the initial stem diameter (t = 0) and may be an alternative to the prediction of Enquist et al. (Eq 3).
Moreover, the diameter growth rate, dD/dt % (D T -D 0 )/Δt, can be obtained using Eq (18): By using this equation, we can examine the relationship between diameter growth rate and initial size (D 0 ).

Data set
To test the six assumptions of our model and the fit of our model (Eqs 18 and 19) and the model of Enquist et al. [2], we used data sets of 11 secondary deciduous broadleaf forests in Hokkaido, northern Japan (Table 1). Mean annual temperature and annual precipitation during past 30 years at the nearest meteorological stations of 11 forests were 5.5-7.3°C and 887-1415 mm, respectively. Permanent plots between 25 m 2 and 2500 m 2 in area were set. The number of trees per hectare ranged from 2,000 to 200,000 (Table 1). Stand age and maximum stem diameter in each plot were highly correlated (r = 0.978, P<0.001). Therefore, we used stand age and the maximum stem diameter as a surrogate of stand maturity and named the plots in order of maturity: plot 1 the youngest and plot 11 the most mature.
Plots 1 to 5 were young forests that naturally regenerated after scarifications, dominated by birch trees (Betula ermanii and B. maximowiczii) [10,23,24] (S1 Table). Plot 5 was artificially thinned. Plots 6, 7, 8, 9, and 10 were relatively mature mixed forests of birches and some other tree species such as oak (Quercus crispula), maple (Acer mono and A. japonicum), and others that naturally regenerated after forest fires. Plot 11 was the most mature secondary forest that had been regenerated after fuel wood production and composed of oak, Carpinus laxiflora, Prunus sargentii, and other tree species. A few conifer trees (Abies sachalinensis) grew in plot 10 and 11.
Stem diameters at breast height (1.3 m aboveground) were measured using a measuring tape or a caliper for a 1-or 2-year interval for 7 to 15 years. The position of diameter measurement was marked with paint. Tree height was measured only in plots 1-5 for all censused years. All trees in each plot were analyzed together, regardless of species, but trees that died during the census were excluded from the analyses. When D T < D 0 , diameter growth (dD) was set as 0.

Analysis
To test the equation of Enquist et al. [2] (Eq 3), we conducted reduced major axis (RMA) regression of the 2/3 power of the initial stem diameter (D 0 ) to the 2/3 power of the diameter at final census t = T (D T ). RMA regression, one of model II regression methods, was used because the relationship between D 0 and D T does not imply a direct causal relationship [25]. The lmo-del2 function of R 2.15.3 was used to calculate the slope and intercept for RMA regression. We tested the six assumptions used to derive our model. To test assumption 1 and 3-6, individual tree aboveground biomass was estimated by using an allometric equation. Although site-and species-specific allometric equations may give most accurate biomass estimates, we did not have an allometric equation for the studied forests and for each species. Recent studies [26,27] showed that generic multi-species equation created from compiled data is the best available choice when site-specific equation is not available. Thus, we decided to use a generic allometric equation created for natural forests in Japan [27]: where M is aboveground biomass of a tree, D is stem diameter at breast height, δ is species-specific wood density, and CF is the correction factor, which is 1.029 in this case. This equation was created from the compiled data of 1203 trees belonging to 102 species (60 deciduous angiosperm, 32 evergreen angiosperm, and 10 evergreen gymnosperm species) harvested from 70 natural forests all over Japan. This equation provided a similar biomass estimate as an allometric equation that had tree height as an additional predictive variable [27].
To test the assumption 1 (Y-N curve fits to the size-frequency distribution of trees), Eq (4) was fitted to the data by non-linear regression for each plot and each year via nls function of the R software.
To test the assumption 2 (a tree's rank does not change with time), Spearman's ρ was calculated between the rank at the initial census and that at the last census.
Assumptions 3 and 4 were tested in plots 1-5 where tree height was measured. In plots 6-11, tree height was too high to be measured accurately and only assumption 5 that was derived from the assumption 3 and 4 was tested. To test the assumption 3 (stand biomass follows a power function of the height of the highest tree in the stand), we calculated stand biomass by summing aboveground biomass for all trees in each plot for each census year. Logtransformed stand biomass was regressed to the log-transformed height of the highest individuals by ordinary least square (OLS) regression. To test the assumption 4 (allometric relationship between tree biomass and tree height), log-transformed aboveground biomass was regressed to log-transformed tree height by OLS regression for each plot and census year. Note that the above generic allometric equation used to estimate aboveground biomass [27] did not include tree height as a predictive variable, so the relationship between tree biomass and tree height, if it does exist, is not an artifact.
We tested the assumption 5 about the relationship between A and B by fitting ln (B) = ln (K 3 ) + c × ln (A) to the parameter estimates of A and B obtained after fitting Y-N curve (assumption 1) by ordinary least square (OLS) regression and parameters K 3 and c were obtained.
To test the assumption 6 (individual tree biomass scales with the 8/3 power of stem diameter), ln (M) = ln (K 4 ) + e × ln (D) was fitted to log-transformed tree aboveground biomass (M) and log-transformed stem diameter (D) by RMA and OLS regressions. Because M was estimated from the allometric equation, which is the polynomial function of stem diameter, the scaling of M to D may be an artifact. Therefore, we also tested this assumption by using another dataset that had been used to create the generic allometric equation in Japan [27,28]. The dataset included measured aboveground biomass and stem diameter of 1203 trees destructively harvested from 70 boreal, temperate, and subtropical forests in Japan (hereafter called harvested dataset). This dataset was also used to additionally test assumption 4.
Using the estimates of parameters A 0 , A T , K 3 , K 4 and c, we tested the fit of Eqs (18) and (19) to the 11 plots. We did not fit these equations directly to the data because these equations are highly flexible and can generate various curves depending on the five parameters. Thus, we first obtained the plot-specific estimates of A 0 , A T , K 3 , and c through tests of assumptions 1 and 5 (S2 Table). For K 4 , we did not use plot-specific estimate because it seemed to be affected by artifact (see previous paragraph). Instead, we used a common value that was estimated from harvested dataset by fitting ln (M) = ln (K 4 ) + 8/3 × ln (D). Then, we fitted Eqs (18) and (19) using these fixed parameter values.

Evaluation of the model of Enquist et al. (1999)
Eq (3) proposed by Enquist et al. [2] fitted well only to some mature stands. Assumption 1: Y-N curve fits to the size frequency distribution of trees Y-N curve between cumulative biomass in the stand summed from the largest tree to the N-th largest tree (Y) and tree's rank (N) fitted well to the actual size structure in all stands (Fig 2).

Assumption 2: Tree's rank does not change with time
Individual tree's rank at the initial census and that at the final census were strongly correlated (Spearman's ρ = 0.903-0.996, P<0.001). Therefore, we can assume no or few changes in a tree's rank with time.

Assumption 3: Stand biomass follows a power function of the height of the highest individual in the stand
In plots 1-5 where tree height was measured, stand total biomass was well expressed by Y max = K 1 H max a+1 (Eq 8, Fig 3, P < 0.002, adjusted r 2 = 0.74-0.95). The 95% CI for a + 1 was larger than 1 in all five plots. Therefore, assumption 3 was supported for all the five plots.

Assumption 4: Allometric relationship between tree biomass and tree height
Although the generic allometric equation used to estimate aboveground biomass did not have tree height as a predictive variable, estimated aboveground biomass showed an allometric relationship to tree height in plots 1-5 in each year, i.e., M = K 2 H b , with b = 1.29-3.68 significantly  different from 0 (Fig 4, P < 0.0001, adjusted r 2 = 0.56-0.95, except r 2 = 0.24 for 1 year in plot 1). Measured tree biomass of harvested dataset [27] also showed an allometric relationship to tree height (P < 0.0001, adjusted r 2 = 0.89, b = 3.28). Therefore, assumption 4 was supported for all five plots and for the compiled harvested dataset.   Assumption 6: Individual tree biomass scales with the 8/3 power of stem diameter Estimated tree aboveground biomass in the 11 forest plots showed a scaling relationship to stem diameter (Fig 6). However, the 95% CI of scaling exponent e of M = K 4 D e did not include 8/3, and the estimates were smaller than 8/3 (i.e., %2.67) in all years and plots (e = 1.71-2.54 by OLS regression, 1.74-2.55 by RMA regression). In particular, younger stands had a smaller scaling exponent (Fig 6; S1 Fig). Measured biomass of harvested dataset [27] also showed a scaling relationship to stem diameter. Again, the 95% CI of e was smaller than 8/3 (e = 2.33-2.36 by OLS regression, 2.35-2.38 by RMA regression). Therefore, assumption 6 was supported partially.  Table). Eq (18) produced a curvilinear relationship between D 0 2/3 and D T 2/3 from the obvious curve relationship in younger stands to the almost linear relationship in mature stands. The adjusted coefficient of determination was 0.64-0.99. In plots 1-7, Eq (18) explained 2-14% more variation in D T 2/3 than Eq (3) by

Relationship between initial and final stem diameters
Enquist et al. [2] with a slope of unity. For plots [8][9][10][11], there was no difference in explained variation between two equations (0.0-0.8%). Relationship between stem diameter growth rate and initial diameter Eq (19) fitted to the observed convex relationship between the initial stem diameter and diameter growth rate in all plots (Fig 7, black curve, adjusted r 2 = 0.05-0.67, parameter values are given in S2 Table). The model of Enquist et al. [2] for tree growth rate (Eq 2) did not fit to the data (dashed line).

Discussion
Our Eq (19) fitted to the observed convex relationship between growth rate of stem diameter and the initial stem diameter better than that of Enquist et al. [2] for all forests (Fig 7). Convex relationship between growth rate and size was shown to suggest asymmetric competition, while the linear-like relationship suggesting no competition or symmetric competition [7,29]. Indeed, asymmetric competition is the theoretical basis underpinning our model. Eq (18) was derived from Eq (4), which was mathematically derived from the growth model assuming asymmetric competition [17]. On the other hand, the model of Enquist et al. [2] is based on the assumption that tree growth is linear to leaf biomass. This means that both small and large trees have the same resource availability per unit leaf biomass. Such assumption may not be always true especially in closed forests where asymmetric competition for light is prominent than symmetric competition [6][7][8]. Smaller trees obtain less light resource per unit leaf biomass due to shading by larger trees resulting in reduced diameter growth than expected by the model of Enquist et al. (Fig 7). Thus, our results suggest that asymmetric competition is an important process in tree diameter growth.
The effect of asymmetric competition on tree growth can be seen in the difference between plots 4 and 5. Plot 5, which was adjacent to plot 4, was heavily thinned (95% of individual trees were thinned). The curve by Eq (19) was less convex in plot 5 than plot 4 (Fig 7). This suggests that smaller trees are less suppressed in their growth in plot 5 than those trees in plot 4, which may be the result of weakened asymmetric competition due to heavy thinning.
The model by Enquist et al. [2] and our model both fitted to the linear-like relationship beteween D T 2/3 and D 0 2/3 in mature forests (Fig 1) even though two models are based on different ecological assumptions. This is because the linear-like relationship between D T 2/3 and D 0 2/3 is a spurious trend resulting from autocorrelation. Since D T = D 0 + TÁdD/dt, the correlation between D T 2/3 and D 0 2/3 should include an autocorrelation of D 0 [30]. The degree of autocorrelation may be greater when the contribution of D 0 to D T is greater. In other words, the autocorrelation is stronger in mature forests, as shown in Fig 1. To avoid artifact due to autocorrelation, it may be better to analyze the relationship between diameter growth rate (dD/dt) and initial diameter (D 0 ) rather than between D 0 and D T . The theoretical relationship between ΔD/Δt and D 0 was easily obtained as Eq (19). Our model is a non-spatial model that considers asymmetric competition. Various growth models that incorporated asymmetric competition have been proposed [6,8,12,31,32]. Some studies are based on sigmoid growth curve such as logistic model and Birch's model [31,32]. However, there has been debate about the mathematical form of plant growth or even about sigmoid growth in trees [33][34][35][36]. Our model is not affected by the form of growth curve because we start from size frequency distribution (Y-N curve). Similarly, the growth model by Yokozawa & Hara [8] assumed no a priori growth function. Their model can evaluate the strength of symmetric and asymmetric competition from the field data. However, to apply to field data, spatial data of tree distribution was needed [37,38]. The models that extended the metabolic scaling theory [6,12] also required spatial data as well as data of light availability and crown area to parameterize the effect of asymmetric competition. Our model does not require these spatial, environmental, and morphological data and therefore may be applicable to various forest inventory data, which usually lack these detailed information.
Another characteristic of our model is that it deals with averaged growth rate for a given size. It cannot express the variation in individual tree growth due to local crowdedness as recent spatially explicit individual-based modeling can (eg. [39][40][41]). However, "macroscopic equation" like our model may be more useful for the future prediction of forest dynamics using dynamic global vegetation model [3].
Eq (19) described the curvelinear relationship of diameter growth rate to initial diameter (Fig 7). However, in some stands, the r 2 value was not high. There may be some reasons to dilute the size dependency of diameter growth. We used the assumption that tree biomass scales to the 8/3 of diameter, but this relationship was not fully supported as shown by previous studies [42,43]. The growth rate may be different among species or growth form even if their size is same [12,30,44]. Because our model is based on asymmetric competition, our model may not perform well in forests where symmetric competition predominates or where interaction among individuals is weak. Examples of such forests may be the forests that have low tree density and crowdedness [38] or forests developed at low-productivity or high-disturbance sites [6,45]. Future researches are needed to improve our model or test the generality of our model.