Effect of strip shelterwood-cuts on the crown morphology plasticity of natural regenerated Pinus tabuliformis saplings in northeastern China

Study area

The research was conducted at Liaoning ecological and experimental farms located in the semiarid areas of the western part of Liaoning Province, northeastern China (120°15′–121°18′E, 41°23′–42°17′N). The area has a continental monsoon climate in the northern temperate zone. The mean air temperature of this area is 8.3 °C, with a maximum temperature of 40.7 °C and minimal temperature of −26.4 °C with in a year. The mean annual precipitation ranges from 450 mm to 550 mm, with approximately 70% distributed within July and August. The elevation ranges from 200 m to 1074 m, and the frost-free period is usually 145–150 days. Brown forest soil is the main soil type in this area.

Our study area is located in the mature Pinus tabuliformis plantation that is 45 years old that has a density of 1,035 trees per hectare in the Linghe experimental region of the Liaoning ecological and experimental farm. In the spring of 2014, the shelterwood-cut strip experiment was conducted. Part of the framework of the research is shown in Fig. 1. The widths of the shelterwood-cut strips varied from 13 to 25 m, and the widths of the uncut stand strips were 5–20 m. Timber and slash were removed immediately after the strip-cut process. The surface soil and the regenerated saplings were not destroyed during the strip-cut work. To reduce the effect of strip width on our results, all the measurements were taken from shelterwood-cut strips with widths of 22 m and 25 m and from uncut strips with widths of 17 m and 22 m. The orientation of the shelterwood-cut strips were in the south–north direction. The average tree height of the mature pine was 15.9 m. All the saplings in this area naturally regenerated within the stands.

Data collection

At the end of the growing season in November 2017, a total of 50 saplings from the shelterwood-cut strip and 30 saplings from the uncut strip were selected. All the saplings were selected from the 2-m wide belts from the centre of the shelterwood-cut and uncut strips to ensure homogeneous conditions regarding light availability for the saplings from shelterwood-cut and uncut strips. For all saplings, the total tree height (HT, 0.1 m) was greater than 1.3 m, and the diameter at breast height (0.1 cm) was less than 5.0 cm. HT and DBH were measured before felling for all the sampled saplings (Table 1).

For the sampled saplings, the whorl containing at least one branch that was continuous to the former whorl was defined as the crown base, and the length from the sapling tip to the crown base was the crown length (CL). A straight line from the ground level to breast height was marked on the trunk to indicate the north side, which facilitated the determination of the azimuth of the branches, and the line was extended to the sapling tip after the saplings were carefully felled. After felling, the HTs of the saplings were remeasured. To facilitate the measurement, the saplings were divided into 1-m sections, while any section less than 1 m was considered a sapling tip. Each section was placed upright on the ground. For each branch from each whorl from every sapling, the vertical distance along the trunk from the branch base to the sapling tip (L) was measured first. Then, the branch diameter (BD), branch length (BL), branch chord length (BC), branch angle (VA) and azimuth, which denotes the angle deviated from the north side in the clockwise direction, were measured (Fig. 2). To detect the effect of shelterwood-cut strips on the annual growth of branch length, one 5-year-old sample branch with from the fifth whorl was selected from each sapling. The annual growth of the branch length for all sampled branches was measured.

To detect possible recording or measurement errors in the field work, a visual inspection of the crown profiles was carefully implemented for each sapling. One specific sapling from the shelterwood-cut strip with a DBH of 0.73 cm and a tree height of 1.4 m showed a strange profile and was considered biologically unreasonable. Thus, this sapling was removed from the data set, leaving 49 saplings from the shelterwood-cut strips with a total of 507 whorls and 1,398 branches for the analysis and modelling. For the 30 saplings from the uncut strips, a total of 363 whorls and 753 branches were harvested. The absolute depth into the crown radius of interest (DINC) was calculated as the difference between the distances from the sapling tip to the branch base and the projection of the branch onto the trunk. The relative depth into the crown radius of interest (RDINC) was calculated as RDINC = DINC/CL. The outermost crown radius (OCR) was calculated as the branch chord length multiplied by the sine value of the branch angle. The characteristics of the saplings and branches from the shelterwood-cut strips and uncut strips are shown in Table 1.

Parametric regression approach for the crown profile model

For the parametric regression approach, the models that perform well with other species when targeting crown profile, largest branch radius distribution and taper description of adult trees were selected as the candidate models. Due to the irregularity of the crown profiles for the saplings, the models used in the description of the crown profiles of the adult trees were reconstructed in the present study (Roeh & Maguire, 1997; Gao et al., 2019).

The power-exponential equation performs well in modelling the Pinus sylvestris var. Mongolica plantation and was used in the present study (Sun, Gao & Li, 2017). The basic form of the power-exponential equation is shown in Eq. (1).

(1) $$\mathrm{O}\mathrm{C}\mathrm{R}=a_1\cdot \mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}{\mathrm{C}}^{a_2}{\mathrm{e}}^{a_3\cdot \mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{C}}$$where OCR is the outermost crown radius within the crown, RDINC is the relative depth into the crown of interest, and a₁–a₃ are the model parameters to be estimated.

The Kozak equation is a continuous variable-exponential equation with a changing exponent from the ground to the top of a tree that describes the different shapes within the entire bole (Kozak, 1988). Due to its flexibility, it has been modified to be applied to the trend in branch diameter within acrown and the crown profile description (Maguire, Johnston & Cahill, 1999; Weiskittel et al., 2010; Gao, Bi & Li, 2017). The basic model form is shown in Eq. (2).

(2) $$\mathrm{O}\mathrm{C}\mathrm{R}=c_1\cdot X^{c_2}=c_1\cdot {\left[{\displaystyle \frac{1-(1-\mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{C}{)}^{0.5}}{1-p^{0.5}}}\right]}^{c_2}$$where the base X is a function of RDINC, the point of inflection p and the two other parameters c₁ and c₂ are variables to be estimated. The OCR is zero at the tip of the tree when RDINC equals zero and reaches the maximum, that is the largest crown radius at the point of inflection when RDINC equals p.

Both basic models (Eqs. (1) and (2)) were further modified or tested in terms of their performance in crown profile modelling for the saplings. Ordinary least square estimation (OLS) was used to estimate the model parameters on the R platform (R Core Team, 2018). The dataset including the branch with the largest radius from each whorl was used to fit these candidate models. The goodness-of-fit statistics of adjusted R² (R_adj²) and root of mean square error (RMSE) (Eqs. (3) and (4)) were used to compare the models. The model with the largest R_adj² and the smallest RMSE was selected as the best model to describe the outermost crown profiles of the Pinus tabuliformis saplings.

(3) $${R_{\mathrm{a}\mathrm{d}\mathrm{j}}}^2=1-\left(1-R^2\right)\left({\displaystyle \frac{n-1}{n-p}}\right)=1-\mathrm{S}\mathrm{S}\mathrm{E}\left({\displaystyle \frac{n-1}{\mathrm{S}\mathrm{S}\mathrm{T}}}\right)$$ (4) $$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{\sum _{i=1}^M\sum _{j=1}^{n_i}{(\mathrm{O}\mathrm{C}{\mathrm{R}}_{ij}-\mathrm{O}\hat{\mathrm{C}}{\mathrm{R}}_{ij})}^2}{n-1}}}$$where $R^2=1-\sum _{i=1}^M\sum _{j=1}^{n_i}{(\mathrm{O}\mathrm{C}{\mathrm{R}}_{ij}-\mathrm{O}\hat{\mathrm{C}}{\mathrm{R}}_{ij})}^2/\sum _{i=1}^M\sum _{j=1}^{n_i}{(\mathrm{O}\mathrm{C}{\mathrm{R}}_{ij}-\mathrm{O}\overline{\mathrm{C}}{\mathrm{R}}_{ij})}^2$, OCR_ij is the outer crown radius of the jth branch from the ith sapling, $\mathrm{O}\hat{\mathrm{C}}{\mathrm{R}}_{ij}$ is the predicted outermost crown radius of the jth branch from the ith sapling, M is the number of saplings, and n_i is the number of branches from the ith sapling.

Nonlinear quantile regression with an interior point algorithm proposed by Koenker & Park (1996) was used to model the outermost crown profiles. The general form of the nonlinear quantile regression model is shown in Eq. (5), and the nonlinear quantile regression estimator was defined as in Eq. (6) (Koenker, 2015).

The estimation for the parameters of the specific quantile was implemented by the quantreg package of R software (Koenker, 2018). In the present study, the outermost crown profile model using a nonlinear quantile regression approach with the same deterministic nonlinear function is shown in Eq. (7).

(7) $$\mathrm{O}\mathrm{C}{\mathrm{R}}_{ij}=f(X,\theta )+{\epsilon }_{ij}$$where OCR_ij is the outer crown radius of the jth branch from the ith sapling, f(·) is the nonlinear equation form, X represents the variables, θ represents the parameters, and ε_ij is the error term.

Outermost crown profile description for the saplings

A nonlinear quantile regression model was used to describe the outermost crown profiles of the Pinus tabuliformis saplings from the shelterwood-cut and uncut strips. Nonparametric boundary regression with a polynomial frontier estimator was used to determine the specific quantile for the nonlinear quantile regression to model the outermost crown profiles. For the nonparametric boundary regression approach, the data pairs (RDINC_i, OCR_i) where i indicates the number of branches per tree were used to characterise the branch tips. A single polynomial equation with a known degree parameter p was used to address the full observations as follows (Abdelaati, Thibault & Hohsuk, 2017): (8) $$\mathrm{O}\mathrm{C}\mathrm{R}={\phi }_{\theta }(\mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{C})={\theta }_0+{\theta }_1\cdot \mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{C}+\cdots +{\theta }_p\cdot {\mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{C}}^p\mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{C}\in [a,b]$$where a and b are the lower and upper endpoints, respectively, of the design points RDINC₁,…, RDINC_n, where n indicates the number of divisions for the interval of [a, b]. In the present study, a equalled 0, and b equalled 1 according to the biological characteristics of the crown. The optimal polynomial degree was obtained by minimising AIC and BIC. The nonparametric boundary estimator and the optimal polynomial degree were estimated by the new R software package ‘npbr’ (Abdelaati, Thibault & Hohsuk, 2017; Koenker, 2018).

A set of quantile values ranging from 0.5 to 0.99 with an even interval of 0.01 was specified for the nonlinear quantile regression model. Nonparametric boundary regression was used to describe the outermost crown profiles for all the saplings from the shelterwood-cut strips and uncut strips. RDINC was divided into 10,000 intervals by an even interval of 0.0001. The crown radius of all the specified RDINCs was calculated for each sapling by a series of nonlinear quantile regression models and nonparametric boundary regression. The absolute distance between the two curves modelled by nonlinear quantile regression and nonparametric boundary regression for each sapling was calculated (Gao, Bi & Li, 2017). The specific nonlinear quantile regression model with the smallest distance to the nonparametric boundary was selected to describe the outermost crown profile for each sapling (Gao, Bi & Li, 2017). Finally, the quantiles used to describe the outermost crown profiles of the saplings for the shelterwood-cut strips and uncut strips were calculated by averaging all the outermost quantiles for the two types of strips.

Effect of shelterwood-cut strips on the crown morphology plasticity of saplings

The outermost crown profile model developed by the nonlinear quantile regression approach was used to depict the crown morphology plasticity of the saplings from shelterwood-cut and uncut strips. The inflection points were calculated for each sapling. The difference in the inflection points of the saplings from shelterwood-cut strips and uncut strips was compared. The crown volume for the upper crown was calculated by a rotational integral using the predicted outermost crown profile curves. The crown morphology plasticity, including the largest crown radius, inflection points, and annual growth of branch length of the two types of saplings, was also compared.

Basic model selection for the outermost crown profile

All parameters of the power-exponential equation and Kozak equation (Eqs. (1) and (2)) were reparameterized to the sapling variables. R_adj² and RMSE were used to compare different models, and the best model was selected for each basic equation, as shown in Eqs. (9) and (10).

For the power-exponential equation, a₁, a₂, and a₃ were reparameterized by introducing easily measured sapling variables. The final model form is shown in Eq. (9).

(9) $$\mathrm{O}\mathrm{C}\mathrm{R}=a_1\cdot \mathrm{D}\mathrm{B}{\mathrm{H}}^{a_2}\mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}{\mathrm{C}}^{(a_3+a_4\cdot \mathrm{C}\mathrm{R})}\mathrm{e}\mathrm{x}\mathrm{p}(a_5\cdot \mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{C})$$where OCR is the outermost crown radius within the crown, DBH is the diameter at breast height, CR is the crown ratio, a₁–a₅ are the parameters to be estimated, and the other variables are defined as those above.

For the Kozak equation, we further modified the basic equation to fit the crown profiles of the saplings, and it is shown as follows (Eq. (10)): (10) $$\mathrm{O}\mathrm{C}\mathrm{R}=c_1\cdot {\mathrm{D}\mathrm{B}\mathrm{H}}^{c_2}\cdot {\left(\frac{1-(1-\mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{C}{)}^{0.5}}{1-(c_3\cdot {\mathrm{C}\mathrm{R}}^{c_4}{)}^{0.5}}\right)}^{(c_5\cdot (1-\mathrm{R}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{C}))}$$where c₁–c₅ are the parameters to be estimated, and the other variables were defined as those above.

Parameter estimates and goodness-of-fit statistics for Eqs. (9) and (10) are shown in Table 2. The R_adj² for the power-exponential equation and the modified Kozak equation with only DBH and CR as independent variables was 0.62. The power-exponential equation was much more flexible in describing the outermost crown radius and was thus selected as the best model to describe the outermost crown profiles for the Pinus tabuliformis saplings.

Outermost crown profile model development using nonlinear quantile regression

The power-exponential equation was used as the basic equation to model the outermost crown profiles for the Pinus tabuliformis saplings. The nonlinear quantile regression models were developed for a set of quantile values from 0.5 to 0.99 with an even interval of 0.01 for the saplings from the shelterwood-cut strip and uncut strip. All sapling branches from the shelterwood-cut strip and uncut strip were used to estimate the model parameters. The estimates of the parameters for the 0.75, 0.85, 0.95, and 0.99 quantiles of the shelterwood-cut and uncut strips are shown in Table 3.

The results indicated that the crown curves based on the nonlinear quantile regression with the same quantile value for the saplings from shelterwood-cut and uncut strips were different in terms of parameter estimates. The distances between the predicted crown profiles from the nonlinear quantile regression with crown quantiles from 0.50 to 0.99 with aneven interval of 0.01 and the fitted outermost crown profiles by nonparametric boundary regression were calculated for all saplings. For the shelterwood-cut strip saplings, the nonlinear quantile regression crown profile model with the quantile of 0.88 had the smallest distance with the nonparametric boundary regression approach. For the uncut strip saplings, the smallest distance was in the 0.92 quantile. The estimates for the two nonlinear quantile regression models for the outermost crown profilesare shown in Table 3 and were used to compare crown morphology plasticity.

Comparison of crown morphology for saplings between shelterwood-cut and uncut strips

With the nonparametric boundary regression, the nonlinear quantile regression was powerful and suitable in describing and predicting the outermost crown profile for the individual saplings from both the shelterwood-cut (Fig. 3) and uncut strips (Fig. 4). Thus, nonlinear quantile regression was used in determining the crown morphology plasticity for the two types of saplings.

The relationship between the characteristics of the saplings and crown regularity was compared between the shelterwood-cut strip and uncut strip saplings. The CR was fixed at 0.80, approximately the mean CR value of our sampled saplings, and the relationship between the crown radius and DBH for the two strip types was analysed (Figs. 5A and 5B). Similarly, the relationship between the CR and crown radius for the saplings from the two strip types is shown in Figs. 5C and 5D, when DBH was fixed at the approximate mean value of 2.0 cm. Under the condition of fixing CR, the crown radius increased with increasing DBH for all saplings from the shelterwood-cut and uncut strips. The crown radius values of the shelterwood-cut saplings werelarger than those of the uncut saplings within the range of 0–0.2 for the RDINC. The crown radius values in which the RDINC were larger than 0.20 for the saplings from the shelterwood-cut strips were smaller than the crown radius values from the uncut strip, and the difference decreased with increasing DBH. The crown radius values increased with increasing CR when the DBH was fixed, while the CR had a much more significant effect on the shelterwood-cut strip sapling crowns than on the uncut strip sapling crowns (Figs. 5C and 5D). With increasing CR, the largest crown radius for both shelterwood-cut strip and uncut strip saplings only slightly increased.

The crown characteristics of the saplings from shelterwood-cut strips and uncut strips were analysed. For the largest crown radius, saplings from the uncut strips were larger than those of the saplings from the shelterwood-cut strips when the saplings from the two strips were the same size (Figs. 5A and 5B). The inflection points of the shelterwood-cut strip sapling crowns were larger than those of the uncut strip sapling crowns. The mean value of the inflection points was 0.81 for the shelterwood-cut strip saplings and 0.65 for the uncut strip saplings. The crown volume of the small uncut strip saplings was larger than that of the shelterwood-cut strip saplings of the same size, and the difference decreased with increasing sapling size.

The annual growth of branch lengths for the saplings from the shelterwood-cut strips and uncut strips were analysed and compared, as shown in Fig. 6. In 2013, before the implementation of the shelterwood-cut strip experiment, the mean annual growth of branch length for shelterwood-cut strip saplings was 11.2 cm, and that of the uncut strip saplings was 11.5 cm. Under similar growing conditions, there was almost no difference in the annual growth of branch length between the two types of saplings. In 2014, the mean annual growth of branch length was 12.9 cm and 14.0 cm for the shelterwood-cut strip and uncut strip saplings, respectively, and the annual growth of branch length for the uncut strip saplings was significantly larger than that for the shelterwood-cut strip saplings. The differences in annual growth between the two types of saplings continued to increase until 2015. However, there was no significant difference in the annual growth of branch length before 2015 (p > 0.05). After 2015, the growth rate of the shelterwood-cut strip saplings increased and was significantly larger than that of the uncut strip saplings (p < 0.05).