Self-thinning tree mortality models that account for vertical stand structure, species mixing and climate Forest Ecology and Management

Self-thinning dynamics are often considered when managing stand density in forests and are used to constrain forest growth models. However, self-thinning relationships are often quantified using only data at a con-ceptualised self-thinning line, even though self-thinning can begin before the stand actually reaches a self- thinning line. Also, few self-thinning relationships account for the effects of species composition in mixed-species forests, and stand structure such as relative height of species (in mixtures), and/or size or age cohorts in uneven-aged forests. Such considerations may be important given the effects of global climate change and interest in mixed-species and uneven-aged forests. The objective of this study was to develop self-thinning relationships based on changes in the tree density relative to mean tree diameter, instead of focusing only on data for state variables (e.g. tree density) at the self- thinning line. This was done while also considering how the change in tree density is influenced by site quality and stand structure (species composition and relative height). The relationships were modelled using data from temperate Australian Eucalyptus plantations (436 plots), subtropical forests in China (88 plots), and temperate forests in Switzerland (1055 plots). Zero-inflated and hurdle generalized linear models with Poisson and negative binomial distributions were fit for several species, as well as for all-species equations. The intercepts and slopes of the self-thinning lines were higher than many published studies which may have resulted from both the less restrictive equation form and data selection. The rates of self-thinning often decreased as the proportion of the object species increased, as relative height increased (species or size cohort became more dominant), and as site (quality) index increased. The effects of aridity varied between species, with self-thinning increasing with aridity index for Abies alba, Pinus sylvestris, Quercus petraea and Quercus robur , but decreasing with aridity index for Eucalyptus nitens, Fagus sylvatica and Picea abies as sites became wetter and cooler. Self- thinning model parameters were not correlated with species traits, including specific leaf area, wood basic density or crown diameter – stem diameter allometry. All-species self-thinning relationships based on all data could be adjusted using a correction factor for rarer species where there were insufficient data to develop species-specific equations. The approach and equations developed could be used in forest growth models to calculate how the tree density declines as mean tree size increases, as height changes relative to other cohorts or species, as species proportions change, and as climatic and edaphic conditions change.


Introduction
Many decades after the development of indices of the effects of density on self-thinning and of maximum sizedensity relationships in even-aged monospecific forests (Reineke, 1933;Yoda et al., 1963) there remains debate about data selection and statistical methods for fitting the relationship (Vanclay and Sands, 2009;Burkhart and Tomé, 2012;Trouvé et al., 2017). Moreover, self-thinning relationships for mixedspecies or uneven-aged forests have received much less attention than even-aged monocultures.
Self-thinning relationships are often described as a straight-line relationship between logarithmically transformed tree density versus mean tree size: where N is the tree density (i.e. number of trees per unit area), d q is the quadratic mean tree-stem diameter, v is the mean tree-stem volume, m is the mean tree mass, and h is mean tree height. When these relationships are fitted to data on the "frontier" (outer edge) of a data set, they provide a self-thinning line that is assumed to represent a barrier such that any further increase in mean tree size is only possible if the density of trees declines along the self-thinning line. This theoretical frontier or barrier has been given different names, including the maximum density (Reineke, 1933; Eq. (2)), self-thinning line (Yoda et al., 1963;Eq. (3)), maximum population density (Enquist et al., 1998;Eq. (4)) or natural basal area (Pienaar and Turnbull, 1973;Vanclay, 1994). Additional variations have also been proposed (García, 2009;Eq. (5)), including the competition density rule (Sterba, 1987). These equations are useful simplifications of stand-level allometry, and Eqs. (2) and (3) have been used in forest growth models to predict density-dependent mortality (e.g. Landsberg and Waring, 1997;Burkhart and Tomé, 2012;Härkönen et al., 2019). A limitation is that contrary to the implications of Eqs. (1)-(5), self-thinning can begin before the stand actually reaches the self-thinning line, not "at" the line (e.g. Drew and Flewelling, 1979), and many methods used for the calibration of the relationships require a selection of data that are assumed to be at the theoretical self-thinning frontier (Vanclay and Sands, 2009;Burkhart and Tomé, 2012). Self-thinning relationships can also vary depending on the methods used to calibrate them, e.g. ordinary least squares, non-linear least squares, reduced major axis, stochastic frontier analysis, quantile regression, segmented regression, generalized linear models with different error distributions (e.g. Poisson, negative binomial) (Pretzsch and Biber, 2005;Zhang et al., 2005;Burkhart and Tomé, 2012;Possato et al., 2016;Trouvé et al., 2017), and they tend to be very data demanding to calibrate (Vanclay and Sands, 2009). Moreover, managed stands do not often lie close to posited self-thinning frontiers, somewhat limiting their practical significance (García, 2011).
To overcome these limitations, several studies have suggested modelling changes in tree density relative to changes in tree size instead of the state variables (Eqs. (1)-(5)) (García, 2009;Vanclay and Sands, 2009;García, 2011;Trouvé et al., 2017): where N is tree density, S is the mean size of the trees, and β 0 , β 1 and β 2 are constants. This enables the simulation of the conditions that cause tree density to decline, even before the stand reaches the theoretical frontier. It also makes use of all data, not only data close to the theoretical frontier, and therefore reduces potential subjectivity of data selection. Eq. (6) can be solved, with the separation of variables method, to obtain N as a function S, from which the following self-thinning line can be derived (Supplementary material S1.2 of Trouvé et al., 2017): where From the solution of Eq. (6), a self-thinning trajectory can be derived that predicts tree density at the end of a period N t2 as a function of the density at the start of the period N t1 and the mean tree size at the start S t1 and end S t2 of the period (Trouvé et al., 2017): The slopes and intercepts of Eqs.
Self-thinning equations calibrated for a specific combination of species (e.g. Binkley, 1984;Puettmann et al., 1992;Pretzsch and Forrester, 2017) are limited to that species combination and often do not quantify the effect of species proportion, which may change between stands and with stand age. Stand-level and cohort-level forest growth simulation models require self-thinning equations that can predict changes in stand density for each cohort (species, size class, age class) and any combination of species, size classes or age classes that it is competing with.
The species or tree size cohort composition can be quantified as the proportion of the stand density (for example as tree density, or basal area) contributed by that species or size cohort (García, 2013). However, for a given tree size (diameter, height, etc.), different species have different competitive abilities, so self-thinning patterns for a given species will probably differ depending on whether it is mixed with a strong competitor, or a weaker competitor. Based on an analysis of 140,000 plots across the world Kunstler et al. (2016) found that as stemwood density increased, there was a decrease in the growth response to neighbourhood basal area from other species and the competitive effect on other species increased. Although stem-wood density can vary for a given species (e.g. with growth rate, age, resource availability, climate and silviculture (Beadle et al., 2011;Pretzsch et al., 2018;Rocha et al., 2020)), it may still quantify some of the effects of species composition on self-thinning. Species composition could be quantified as a species proportion (p) in terms of the contribution the given species makes to the stand basal area weighted by the species competitive ability as expressed by wood basic density. Similarly, the mean stand wood basic density was used to quantify stand density indices in mixed forests of North America (Woodall et al., 2005;Ducey and Knapp, 2010).
Vertical stand structure may also need to be accounted for because the mortality rates for a given species or size cohort could depend on canopy position (dominant, subdominant, etc). Vertical structure can be quantified as the relative height (rh), which is the mean height of the object species (or cohort) divided by the mean height of all trees in the stand.
The use of species or size cohort proportions (p) and relative height (rh) would extend Eq. (6) to Like Eq. (6), Eq. (11) can be solved to obtain N as a function of S (and rh and p), and from this solution the intercept parameter for the selfthinning line of Eq. (7) can be derived (the slope parameter is as shown by Eq. (9)): The self-thinning trajectory function becomes Self-thinning equations also need to be appropriate for stands with different size or age distributions, such as from even-aged forests with unimodal-shaped distributions to single-tree selection forests with reverse-J-shaped distributions. That is, except where there is no significant relationship between tree size and mortality rate, the shape of the size distribution will influence the stand mortality rate because the size distribution defines how much of the stand density is composed of a given tree size, and hence the proportion of the stand composed of tree sizes susceptible to mortality (Forrester, 2019). That size distributions influence self-thinning patterns has been known for several decades (Sterba and Monserud, 1993), and the effect can be at least partly accounted for by applying Eq. (13) to several cohorts assumed to have roughly unimodal-shaped distributions, instead of treating the whole stand as a single cohort. Thus, the response of each cohort can differ depending on the exponents β 3 and β 4 (Eqs. (11) and (13)), and the response of the entire stand will be the sum of all the cohorts.
Dividing a stand into cohorts has also been used to avoid the issue of Jensen's inequality (Jensen, 1906) when calculating stand density indices based on Eq. (2), i.e. the mean of a function is not necessarily the same as the function of the mean. For example, several studies have calculated stand density indices by summing the stand density index for each size class within a stand (Stage, 1968;Long and Daniel, 1990). However, this approach can potentially produce biased estimates by overestimating the relative density of small size classes (Woodall et al., 2003). Summation approaches have also been applied to mixed-species forests (Woodall et al., 2005;Ducey and Knapp, 2010;Ducey et al., 2017).
Eq. (11) may need to account for the effects of climatic and edaphic factors on self-thinning relationships. A potential site effect is consistent with suggestions to use dominant height as the size variable (García, 2009). However, Trouvé et al. (2017) suggested that mean diameter would be appropriate because it is more closely related to crown width and space occupancy.
The fitting of self-thinning equations (e.g. Eqs. (1)-(5)) requires relatively large data sets that include most situations where the equations are to be applied, such as a wide range of stand basal areas, mean diameters, tree densities, species compositions, and vertical structures, while not having any recruitment or applied thinning during the periods examined. Such data sets are not often available for many species or species mixtures (Woodall et al., 2005;Ducey and Knapp, 2010). A possible solution to this problem is to develop all-species equations, where data from many species are combined. This requires that species identity can be quantified using a variable that is correlated with inter-specific differences in self-thinning patterns.
Several studies have found that the basic density of the wood of a species is negatively correlated with maximum stand density indices and hence self-thinning patterns (Dean and Baldwin, 1996;Woodall et al., 2005;Ducey and Knapp, 2010;Ducey et al., 2017;Andrews et al., 2018). The reasoning is that for a given diameter, species with higher wood densities can support larger branches and leaf areas and therefore fewer trees per unit area are required to close the canopy; species with higher wood densities should have lower maximum stand density indices (Dean and Baldwin, 1996). Shade tolerance has also been considered, however, it is less consistently correlated with maximum stand density indices (Lonsdale, 1990;Jack and Long, 1996;Woodall et al., 2005;Weiskittel et al., 2009;Charru et al., 2012;Ducey et al., 2017;Andrews et al., 2018). Shade tolerance is also a complex variable that is unavailable for many species, thereby reducing its practicality. While these traits have been used to develop all-species stand density index equations, we are not aware of any self-thinning equations that have been developed using such traits.
The first objective of this study was to examine how species proportion, relative height, site quality and climate affect self-thinning. It was hypothesized that the intercept of self-thinning lines will be greater when (1) the species proportion of the object species decreases (assuming generally positive mixing effects on stand density), (2) a species is more dominant, in terms of relative height, and (3) as climatic or edaphic conditions improve site quality. The second objective was to examine whether species traits could be used to develop all-species selfthinning equations by adding trait variables to Eqs. (11), (12) and (13).

Experimental sites
Data were obtained from permanent sample plots and long-term silvicultural experiments in natural and planted forests in Australia, China and Switzerland. The Australian data (436 plots from 29 sites) were from Eucalyptus plantation experiments in temperate and Mediterranean climatic regions, examining planting spacing, thinning, fertiliser application or species mixing (Table 1). The Chinese data were from tree-species rich subtropical forests (88 Chinese national forest inventory plots). The Swiss data were from temperate forests that were either unmanaged, including 71 plots from 21 sites in the Swiss Forest Reserve Research (FRR) network (Hobi et al., 2020), or experiments in managed forests, including 984 plots from 594 sites in the Experimental Forest Management (EFM) network . We only included data from plots where there were at least 5 trees of the object species, and it existed as a single cohort (determined by visual inspection of size distributions). This still resulted in negative exponential shaped size distributions in the data set, but was done to avoid introducing potential subjectivity from dividing an object species into multiple cohorts.

Climate data
The effect of climate on self-thinning was quantified using an aridity index (de Martonne, 1926) that was calculated for each site as (mean annual precipitation in mm) ÷ (mean annual temperature in • C + 10). The mean temperatures and precipitation used to calculate the aridity index were obtained: (i) for the Australian sites from interpolated national climate surfaces ANUCLIM MTHCLM (Xu and Hutchinson, 2011), (ii) for the Swiss sites by interpolation (after Thornton et al. (1997)) from data collected by the Federal Office of Meteorology and Climatology MeteoSwiss. For the Chinese data, monthly precipitation and temperature were calculated by spatial interpolation of climate data collected at 369 meteorological stations in Hunan Province (Zhao et al., 2013).

Tree measurements
In the Australian plots, tree diameter at 1.3 m height (d) was measured for all trees when they were at least 1.3 m tall, and tree height (h) was measured on either all trees or a sample spanning their size range. For some sites, crown diameters were also measured for a sample of trees. In the Chinese plots, d was measured on all trees with d ≥ 5 cm, but no heights were measured. In the Swiss EFM plots and FRR plots, d was measured on all trees with d ≥ 8 cm or d ≥ 4 cm, respectively, and tree heights and four crown radii were measured for a sample of trees (the 100 largest-diameter trees ha − 1 and 20% of the rest). Unmeasured tree heights in the Australian and Swiss plots were estimated using plot-, year-and species-specific regressions after Forrester et al. (2019). The measurement intervals for the Australian, Chinese and Swiss plots varied between about one and twelve years, depending on silvicultural treatments and growth rates. In the data we only included growth periods without any recruitment.

Calculations
For each species in each plot, quadratic mean d (d q ), tree density (N, trees ha − 1 ), species proportion (p), and relative height (rh) were calculated for alive trees at each measurement age. Changes in N (ΔN) and d q (Δd q ) between successive pairs of measurements were annualised by dividing by the length of the observation period. Because of the highly diverse species composition of the Chinese plots, N for each species was too small for species-level analyses and therefore the species were not differentiated and p and rh were not calculated. For all other mixed-species plots, the N of a given species was divided by its proportional contribution to stand basal area to provide a monoculture equivalent. This was required because the self-thinning equations are for individual species or cohorts, so regardless of whether they are growing in a mixture or monocultures, the N of an object species needs to represent that N that would exist in a monoculture that has the same total stand basal area as the mixture. The species proportion was quantified in terms of the contribution the object species made to stand basal area (G, m 2 ha − 1 ) weighted by the wood basic density (ρ; oven dry mass per fresh volume; g cm − 3 ): This weighting was to represent inter-specific differences in growth responses to neighbourhood basal area (competitive ability) for a tree of a given d, which has been found to be correlated with ρ (Kunstler et al., 2016). For a given proportional contribution to basal area, p will be lower when a species is mixed with a high-wood-density species than when it is mixed with a low-wood-density species. A p of 1 indicates that the species is in a monoculture, or that all species have the same wood density. The ρ for each species was obtained from the literature (Tables A1 and A2).
Because the aridity index will only partly account for all site effects (climatic and edaphic), a site (quality) index was also calculated. While this will be influenced by the climatic factors that influence the aridity index, the aridity index and site index were not correlated for any of the species in this study and both indices were used to make comparisons with previous studies where these indices have been used to examine self-thinning relationships. Site index (SI, m) for each species at each site was quantified using relationships between dominant height and age. In monocultures, the dominant height was the mean height of the largestd 100 trees ha − 1 . In mixtures, to avoid underestimating the site index, the dominant height was calculated from fewer trees: [100 × species proportion in terms of basal area] trees ha − 1 . The relationship between dominant height and age for a given species was estimated by fitting all the data for the species to a von Bertalanffy equation (Pienaar and Turnbull, 1973) where H d is the dominant height (m), A is age (years) and a H , b H and c H are parameters. Eq. (15) was fit as a non-linear model using the nlxb function of the nlmrt package (Nash, 2016) in R (R Core Team, 2019).
The SI was then calculated from where A r is the reference age for SI: 10 years for the Australian plots, which were relatively faster growing, and 100 years for Swiss plots, which were slower growing. SI could not be calculated for the Chinese plots because height data were not available.

All-species equations
Three species traits were used to examine whether they could be used to develop all-species self-thinning equations: (i) wood basic density (ρ), (ii) specific leaf area (SLA, m 2 kg − 1 ), and (iii) the relationship between tree crown diameter (c d , m) and stem diameter (d). These values were obtained from literature review (Table A1). The c d -d trait was determined as the two parameters a C and b C in the equation These parameters were used because for a given d, species with longer branches may require fewer trees per unit area to close the canopy, and thus may manifest inter-specific differences in self-thinning. The c d -d parameters were previously calculated from crown diameter measurements in the EFM and FRR plots for the Swiss species (Forrester et al., 2021), and for E. nitens by Forrester et al. (2012a) (Table A1). E. globulus c d -d parameters were calculated using data from the experiment described in Forrester et al. (2011), which is also included in this study.
To determine any inter-specific differences in self-thinning, the slope and intercept as in Eqs. (1) and (2) were calculated from the parameters of the species-specific equations (e.g. Eqs. (13)) using Eqs. (9) and (12). The slope and intercepts were then regressed against the traits.
The all-species equations were based on Eq. (11), but each explanatory variable was replaced with an interaction between the given variable and an individual trait or trait parameter. In order to include each species in the all-species equation, SI and aridity index were constrained to a maximum of 1. That is, the index for a given species was divided by the maximum value for that species. This normalisation was required because some species can, on average, occur on sites with higher SI or aridity indices than others and thus confound species effects with the SI or aridity effects.

Model fitting
The self-thinning lines (Eqs. (9) and (12)) and self-thinning trajectory functions (Eq. (13)) were calibrated by fitting the rate of change equation (Eq. (11)) as a generalized linear model (GLM) with random effects. Eq. (11) was transformed into by approximating the differentials (dN,dS = dd q ) with differences ( − ΔN, ΔS = Δd q ) and by taking the logarithm on both sides. By considering the mortality count ΔN ijk of observation period k on plot j on site i as a random variable with a discrete probability distribution, Eq. (18) could be described as a GLM with a logarithmic link function and log ( Δd q ) as an offset term: is the mean of the distribution of ΔN ijk , assumed to depend on the observed values of d q,ijk , N ijk , rh ijk , and p ijk in the beginning of the period k and of Δd q,ijk during the period k, and b i and b ij are random effects for the site and for the plot within each site, respectively, assumed to be independent and normally distributed with means 0 and separate constant variances. The random effects account for the dependence between the repeated measurements on a plot, and between multiple plots at the same site. Several models following Eq. (19) but with different probability distributions for ΔN ijk were compared. First, a model with Poisson distribution (Poisson) was fitted (cf. Trouvé et al., 2017). The problem with the Poisson distribution is that it involves only one parameter (mean μ ijk ) and assumes the variance to be equal to the mean, whereas in practice count data are more dispersed. To deal with the problem of overdispersion, a model with a negative binomial distribution (NB) was also fitted (cf. Affleck, 2006;Trouvé et al., 2017). This distribution involves an additional dispersion parameter, which can depend on the mean linearly or quadratically; in this study, the parameterisation with quadratic dependence was used (Var In the data, the number of zero mortality counts was likely to be larger than what would be predicted by a Poisson or negative binomial model. To deal with the problem of zero-inflation, models with zeroinflated, or hurdle, Poisson (ziP, hP) and negative binomial (zNB, hNB) distributions were fitted. These distributions involve a separate process that generates extra zeros with probability π ijk .
A zero-inflated distribution (Lambert, 1992;Greene, 1994) assumes two types of observation periods: those where mortality can never occur (ΔN is always zero; source of extra zeros), and those where mortality can occur but does not always (normal mortality process). Zero mortality counts can thus come from both types of periods. The zero-inflated probability distribution function is where f is the probability distribution function of Poisson distribution with parameter Θ ijk = μ ijk or negative binomial distribution with pa- A hurdle distribution (Cragg, 1971;Mullahy, 1986) assumes that the observation periods are either where mortality can never occur, or where mortality always occurs (at least one dead tree is observed). Here, zero mortality counts can thus come only from periods of the first type. The hurdle probability distribution function is where the notation is as above. When using zero-inflated and hurdle distributions, the probability π ijk generating extra zeros was modelled with a GLM involving Bernoulli distribution, logit link function and log or both as the explanatory variables. The models were fitted using the glmmTMB function in the R package glmmTMB (Brooks et al., 2017). An example R script to fit these models is provided as supplementary material.
To evaluate the goodness-of-fit for each model, we calculated the mean of prediction errors the root mean square error and the coefficient of determination where y k and ŷ k are the observed and the fitted values of mortality count ΔN k for period k. All analyses were performed using R software (R Core Team, 2019). The best models for each species were generally considered those with the lowest AIC and BIC, but visual inspection of predicted and observed ΔN, and the other statistical information in Table 2 were also used.

Results
Self-thinning trajectories were fit for species with at least 100 pairs of successive measurements (Table 2). Fits using smaller sample sizes were less reliable, based on the statistics used to judge the quality of fit and a visual inspection of the predicted self-thinning trajectories.
While there were often differences in parameter estimates depending on the models used (Table 2), there were generally only minor differences in the predicted self-thinning trajectories (results not shown), and no fitting type was consistently better or worse than the others. The sign of the parameters, indicating the direction of the relationship, was the same for a given variable and species.
Relative height (rh) was significant for A. alba, Q. petraea and Q. robur, such that mortality declined as the relative height (i.e. dominance) increased (Fig. 1). This could not be tested for the Eucalyptus species or P. sylvestris because there was not enough range in rh for these shade intolerant species. Relative height had no significant effect on the self-thinning of F. sylvatica or P. abies.
Species proportion (p) was significant for 4 of the 5 species where it could be tested (was not significant for P. sylvestris) (Figs. 1 and 2). In all cases the mortality rates declined as the proportions of the object species increased. The aridity index was significant for all species except Eucalyptus (Figs. 1 and 2; Table 2). The mortality rates declined with increasing moisture for P. abies and F. sylvatica, but increased for A. alba, P. sylvestris and Quercus. Site index (SI) was significant for P. sylvestris and E. globulus, and in both cases the rates of self-thinning decreased as Summaries of plot and site characteristics for the 11 object species, as well as the species-diverse subtropical forest in China. d q is quadratic mean of tree diameters at 1.3 m height, p is the proportion of stand basal area contributed by the object species, and rh is the relative height of the object species (mean height of trees of object species divided by mean height of all trees in the plot).   * not calculated due to the high species diversity that resulted in lower sample sizes for any given object species. ** January and July were considered the warmest and the coldest month in the southern hemisphere, respectively, and the coldest and the warmest month in the northern hemisphere, respectively.

Table 2
Parameter estimates and fit statistics for self-thinning models based on Eqs. (11) and (13) ΔN)). Note that for individual species equations the site index and aridity index were not standardised, as opposed to the all-species equation.  site index increased ( Fig. 3; Table 2). None of the species traits or trait parameters examined were related to the species-specific slopes or intercepts (calculated using Eqs. (9) and (12)) (Fig. 4). Therefore, the all-species models that included the traits are not presented. Instead all-species models were fit to all data without differentiating between species (Fig. 5; Table 2). To apply this to species for which there were not enough data, a species-specific correction factor (CF) was calculated as the ratio of mean dN of measured data divided by the mean dN of the predicted data for the given species. This was applied as shown by Eq. (25) for the subtropical Chinese plots, L. decidua, P. strobus, and P. menziesii ( Fig. 6; Table 2). Since the Chinese data did not contain the stand structural variables (rh and p), the allspecies equation used also did not contain these variables, and is presented in Table 2. For the three conifers, the CF was negative indicating that before applying the CF, the all-species model overestimated rates of self-thinning, whereas for the Chinese data, the all-species model underestimated rates of self-thinning.   Table 2. The sd is 1 standard deviation of the indicated variables.  Table 2. The sd is 1 standard deviation of the indicated variables.  Table 2. The sd is 1 standard deviation of the indicated variables.  Inset for (d) shows non-transformed data. The equation used is described in Table 2. The sd is 1 standard deviation.

Fig. 4.
Scatter plots indicating the lack of relationships between species traits: specific leaf area (SLA), wood basic density (ρ), and c d -d parameters (a C and b C from Eq. (17)) and the intercept and slope (e.g. of Eqs. (1) or (2)) that were calculated using Eqs. (9) and (12). The round symbols indicate published values (Table A1). The star symbols were the fitted slope and intercept from the present study, which were calculated using rh = 1, p = 1, aridity index = 56 (mean for data set) and SI = 30 (mean for data set).

Discussion
Using rates of change to model self-thinning provided similar results to previous studies that used state variables (Eqs. (1) and (2)) and focused on data close to the self-thinning frontier (Fig. 7). However, the slopes and intercepts in this study were often greater than those reported in the literature (Fig. 8). There is no reason that the equations in this study are less accurate than those from the literature because they were based on large data sets, including many plots that had never been thinned (e.g. E. globulus and E. nitens plots) or only lightly thinned, nor do the statistics suggest a low quality of fit. Instead the greater intercepts and slopes in this study may have resulted from the less restrictive equation form and data selection. Similarly, the non-linear self-thinning frontiers calculated by Charru et al. (2012) appear to be more consistent with those in this study than in the studies that fit linear relationships, and several studies have suggested that self-thinning frontiers are not linear (Zeide, 1987;Cao et al., 2000). The greatest discrepancy was for large tree sizes, such that the linear relationships often predict that stands with d q = 100 cm will have at least 100 trees ha − 1 , whereas the equations in this study predict lower stand densities.

Species mixing
All stand structural and site variables influenced self-thinning relationships. Self-thinning trajectories shifted upwards with increasing species proportion p, indicating that monospecific stands supported higher densities than mixed-species stands. This was consistent, occurring for each species where p effects were significant; A. alba, F. sylvatica, P. abies, and Q. robur or Q. petraea. This p effect was also found for several North American species (Weiskittel et al., 2009). Therefore, while potential maximum stand density is sometimes higher in mixtures for certain complementary combinations of species (Binkley, 1984;Binkley et al., 2003;Pretzsch and Forrester, 2017) the results of this study show that averaged across many species combinations, mixing generally reduced potential maximum stand densities for a given species. The use of p assumes that the effects of species interactions on selfthinning are correlated with the species proportion in terms of basal area weighted by the wood basic density of each species. If necessary, more precise self-thinning equations could be developed by using data from a specific species mixture.

Effects of relative height and uneven-aged stands
Mortality rates declined as rh increased for each species where rh was significant (A. alba, Q. robur and Q. petraea). The higher mortality rates with lower relative heights is consistent with findings of higher survival rates with increasing light availability (Osunkjoy et al., 1992;Inman-Narahari et al., 2016). Relative height is rarely considered in mortality studies, but these results indicate that the dominance of a species or cohort can be as influential as the effects of species proportions and larger, as well as more consistent, than climatic effects such as those quantified using the aridity index. The effect of rh on P. sylvestris or the Eucalyptus species could not be examined due to the small range in rh in the data set for these species. The small ranges resulted because these species were often in monocultures or only occurred in the upper canopy.
By accounting for the effects of species or size cohort composition (p) and the vertical position of each species or size cohort (rh), the equations can be used for different cohorts within uneven-aged forests. Similarly, the influence of size distributions on the self-thinning of a whole stand can be accounted for by calculating, and then summing, the changes in the tree density of several cohorts. This could account for the effects of size distributions on self-thinning (Sterba and Monserud, 1993) caused by size-related differences in physiology, morphology and phenology that influence interactions between different size classes (Forrester, 2019). It can also reduce the bias caused by Jensen's inequality in stand density index calculations by summing the stand density indices of each Fig. 6. The all-species model applied to specific species or the Chinese data after applying a speciesspecific correction factor (CF). The blue lines are the uncorrected all-species model, while the red lines are the corrected model, the grey lines show the raw data for the object species. The parameters for the all-species model (containing only d q and N) are shown in Table 2. The bias reduction (bias red.) is the mean of the differences between observed dN and predicted dN, and was reduced to 0 for all species after applying the CF. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 7. Comparison of the fitted equations in this study (coloured lines, Table 2) with published equations where ln(N) is either a linear or quadratic function of ln(d q ) (Hynynen, 1993;del Río et al., 2001;Pretzsch and Biber, 2005;Charru et al., 2012;Guzmán et al., 2012;Rivoire and Moguedec, 2012;Vospernik and Sterba, 2015;Trouvé et al., 2017;de Prado et al., 2020;Pretzsch and del Río, 2020).
cohort within a stand (e.g. Long and Daniel, 1990). The rh and p exponents (β 3 and β 4 in Eqs. (11) and (13)), additionally allow for a greater range in species or size cohort interactions (i.e. symmetric and asymmetric competition) than approaches where the same exponent (e.g. 1.6 in Eq. (2)) is assumed for all size classes.

Site index and aridity index
Self-thinning occurred at higher stand densities as site index increased for P. sylvestris and for E. globulus, as found in many other studies (Zeide, 1987;Hynynen, 1993;Bi, 2001;Weiskittel et al., 2009;Reyes-Hernandez et al., 2013;Zhang et al., 2013). The aridity index was significant for all species except E. globulus. Climatic variables have been shown to influence self-thinning lines or maximum stand density indices previously (Condés et al., 2017;Aguirre et al., 2018;Andrews et al., 2018;de Prado et al., 2020). The size of the effect was similar to the stand structural variables (p, rh) but inconsistent, such that for P. abies and F. sylvatica the mortality rates were higher at moister and cooler sites, but for the other species they were higher at dryer and hotter sites. If the effects of a given climatic factor on growth depend on other climatic or soil factors (Etzold et al., 2019), the interacting effects of different variables could reduce the effect size of the aridity index (by averaging contrasting effects) or distort the effect. This makes it more difficult to account for climatic effects and requires data from a wider range of sites that include multiple gradients in climate and stand structures. Another complication of using climatic variables could be that temporal climate effects need to be quantified using a temporal resolution that matches the resolution at which climate influences selfthinning (e.g. monthly resolutions for droughts) rather than the approximately 8-year resolution of some of the inventory data used in this study (Landsberg and Sands, 2011). If the resolution is not fine enough, the effects of climate could be averaged out and underestimated.
An argument against including climate in self-thinning equations that are to be used in process-based forest growth models is that the applicability of such self-thinning equations for future climatic conditions is unknown. The use of such equations for predictions under new climatic conditions that were not included in the empirical data, requires an assumption that space-for-time substitution is appropriate. However, this is not necessarily a valid assumption (Trotsiuk et al., 2020). Therefore, although excluding climate may reduce precision for predicting self-thinning in past and current conditions, it may also prevent the addition of bias in simulations of future conditions.

All-species equations
All-species equations based on species traits could not be developed because no traits were correlated with the slope and intercept of the selfthinning lines. This contrasts with studies reporting correlations between maximum stand density indices and traits such as wood basic density, shade tolerance or specific leaf area (Lonsdale, 1990;Dean and Baldwin, 1996;Jack and Long, 1996;Woodall et al., 2005;Weiskittel et al., 2009;Ducey and Knapp, 2010;Charru et al., 2012;Ducey et al., 2017;Andrews et al., 2018). Although the traits could not be incorporated into all-species equations, the all-species equation could be applied to rarer species by calculating a correction factor. This could be used for species where there is not yet enough data to fit species specific equations, but there is some data available to calculate the correction factor (Eq. (25)).

Conclusions
The use of rates of change, such as changes in tree density relative to mean tree diameter, in self-thinning equations provided similar results to previous studies that focused on state variables at the self-thinning line (e.g. Eq. (2)), but allowed the quantification of self-thinning that occurs before the self-thinning lines defined by Eqs. (1), (2) and (7). The intercepts and slopes of these lines were greater than many previous studies, which may have resulted from the less restrictive equation form and data selection. The mortality rates declined with increasing species proportion and relative height within a stand, and site index.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
See Tables A1 and A2.   (2)) that were calculated using Eqs. (12) and (9) and the parameters for the all species equation in Table 2. The round symbols indicate published values (Table A1). The star symbols were the fitted slope and intercept from the present study, which were calculated using rh = 1, p = 1, aridity index = 56 (mean for data set) and SI = 30 (mean for data set).