Simulation and prediction of changes in tree species composition in subtropical forests of China using a nonlinear difference equation system model

Changes in tree species composition are one of the key aspects of forest succession. In recent decades, significant changes have occurred in the tree species composition of subtropical forests in China, with a decrease in coniferous trees and an increase in broad-leaved trees. This study focuses on Zhejiang Province, located in the subtropical region of China, and utilizes seven inventories from the National Continuous Forest Inventory (NCFI) System spanning 30 years (1989-2019) for modeling and analysis. We categorized tree species into three groups: pine, fir, and broadleaf. We used the proportion of biomass in a sample plot as a measure of the relative abundance of each tree species group. A novel nonlinear difference equation system (NDES) model was proposed. A NDES model was established based on two consecutive survey datasets. A total of six models were established in this study. The results indicated that during the first two re-examination periods (1989-1994, 1994-1999), there was significant fluctuation in the trend of tree species abundance, with no consistent pattern of change. During the latter four re-examination periods (1999-2004, 2004-2009, 2009-2014, 2014-2019), a consistent trend was observed, whereby the abundance of the pine group and the fir group decreased while the abundance of the broad-leaved group increased. Moreover, over time, this pattern became increasingly stable. Although the abundances of the pine group and the fir group have been steadily declining, neither group is expected to become extinct. The NDES model not only facilitates short-term, medium-term, and even long-term predictions but also employs limit analysis to reveal currently obscure changing trends in tree species composition.

1 Introduction

As an essential part of the global ecosystem, the forest ecosystem is a complex heterogeneous hierarchical system characterized by dynamic species replacement (Bruelheide et al., 2011). Forest succession can reveal changes in the whole forest, affect the wood supply of the forest (Lantz et al., 2022), and influence ecological functions, such as the supply of ecosystem services (Cortés-Calderón et al., 2021), carbon sequestration (Dulamsuren, 2021), soil carbon accumulation (Zhu et al., 2021), and animal and plant habitation (González et al., 2022). Species composition is a vital expression of forest succession (Mi et al., 2016). Species composition can reflect the composition of tree species in forest ecosystems and the reasonable proportion of each tree species (He et al., 2015; Pelissari et al., 2017). It is necessary to characterize the changes in species composition and predict the forest succession process (Tian et al., 2022), which helps us better understand the effects of human disturbance and climate change on species composition and provide sustainable forest management (Rutishauser et al., 2015).

Mathematical models have been widely used to simulate and analyze the process of forest succession. There is a category known as forest gap models, which primarily focus on small-scale forest succession. The patterns of forest succession, such as the optimal growth, replacement, and death of a single tree or the gap phase dynamics at the plot scale, can be well demonstrated using the succession gap model (Larocque and Bell, 2021). Examples of these models include JABOWA (Botkin et al., 1972), FORET(Shugart and West, 1977; Bormann and Likens, 1979), SORTIE (Pacala et al., 1993; Bugmann, 2001), and FAREAST (Yan and Shugart, 2005). Some of researchers prefer highly parametric models to simulate the photosynthesis and transpiration of individual leaves, as these models can well explain the process of dynamic forest succession (Medlyn et al., 2005; Ibrom et al., 2006). In addition, climate-sensitive gap models were applied to capture the effects of climate and ecological processes on future forest dynamics, such as ForClim (Mina et al., 2015; Martin-Benito et al., 2022), TreeMig (Thuiller et al., 2008), FireBGCv2 (Keane et al., 2011), PnET-Succession(Scheller et al., 2011), and FAREAST (Hu et al., 2018).

With increasing simulation scale, stand-level and landscape-level models are being proposed in forest succession research. For example, LANDIS-II (De Jager et al., 2019), FATES (Koven et al., 2020), FORMIND (Rödig et al., 2017) and SIMREG (Perin et al., 2021) have been applied to simulate the formation, development, disturbance, degradation, and other processes of forests on a wide range of spatial and temporal scales. Some of these researchers in this area have paid more attention to the basal area and compositional changes during tropical forest succession, and demographic forest models have been developed to analyze the growth-survival and stature-recruitment trade-offs (Ouyang et al., 2016; Hu et al., 2018; Ruger et al., 2020). Other researchers have focused on the dynamic disturbances involved, such as the survival rate of tree species in the process of forest succession (Peng et al., 2010), the regeneration rate after disturbance (Blanco et al., 2021), and the occurrence of natural disasters (Yin et al., 2009). The Markov chain model is used to describe the state of the forest ecosystem in forest growth scenarios, in which dynamic disturbances are treated as special events of forest growth scenarios that are constantly transferred according to a certain probability (Risch et al., 2009). Ecological process models are commonly employed for predicting changes in forest tree species composition, such as forecasting the loss and replacement of 111 tree species in China under climate change (Li et al., 2020).

In terms of research content, some researchers have emphatically explored the vertical structure changes of subtropical tree species (Yang et al., 2014), community structure changes (Xu et al., 2020), phylogenetic turnover (Chen et al., 2017; Michalski et al., 2017; Zhang et al., 2021), and ecological memory (Sun et al., 2013). Some researchers are concerned about how interacting forest components affect forest succession (Norden et al., 2015).

Most of the current research on forest succession is based on time series data (Clara et al., 2018; Mathieu et al., 2022). Studies have shown that difference equation theory can provide novel insights into forestry (Tomé et al., 2006; Chen, 2021). For example, Gilbertson and Kot (2021) proposed applying the Pioneer-Climax model based on difference equations to conduct a nine-case steady-state analysis of the long-term behavior of pioneer and climax species. Fort and Grigera (2021) developed a logistic difference model based on linear Lotka-Volterra equations to model and predict the tree species composition of the Barro Colorado tropical forest. Utilizing national forest resource continuous inventory data provides numerous conveniences for simulating and predicting forest tree species composition. For instance, Nikolova et al. (2019) employed generalized linear models based on forest inventory data from the National Continuous Forest Inventory (NCFI) to simulate and predict tree species composition in Swiss forests.

The tree species in Zhejiang Province effectively represent the forest ecosystems of South China and East China (Wu, 1980; Zhang et al., 2022). The distributed forests of Masson pine (Pinus massoniana), mixed conifer-broadleaf forests, and broad-leaved forests within the region represent different stages of forest succession, reflecting the temporal sequence of tree species composition evolution in the subtropical region of China (Sun et al., 2013; Li et al., 2021; Zhang et al., 2011; Piao et al., 2015).

Over the past few decades, due to sustained economic development and increased awareness of ecological benefits, as well as the impact of pine wilt disease (Hao et al., 2022) and a decline in timber prices, the forest species composition in Zhejiang Province has undergone significant changes. The proportion of pine and fir tree species has consistently decreased, while the broadleaf tree species have notably increased (Lin et al., 2012). This shift is not unique to Zhejiang Province but is observed across much of the subtropical regions in China. Simulating and analyzing this change and forecasting its trends is evidently crucial and meaningful. Where will this transformation ultimately lead? Will pine and fir species become extinct? To date, there is limited reporting on macro-temporal studies of tree species composition changes in the subtropical regions of China. This study utilized NCFI data from Zhejiang Province to investigate the relationship and dynamics among three major tree species groups: pines, firs, and broadleaves. The tree species group abundance (referred to as species abundance in some cases) was quantified using the proportion of biomass in each plot. Based on the plot-level species abundance, a nonlinear difference equation system (NDES) model comprising three equations was developed to simulate the temporal changes in tree species composition over the past 30 years (1989-2019) in Zhejiang Province. Furthermore, long-term trends and ultimate scenarios (limits) of future changes were analyzed.

2 Materials and methods

2.1 Study area

Zhejiang Province is located along the southeastern coast of mainland China (118°01’ to 123°10’ E, 27°02’ to 31°11’ N) (Figure 1). The region belongs to the subtropical monsoon humid climate zone. The annual average temperature ranges from 16.1 to 18.6°C, and the annual average precipitation ranges from 1109.1 to 2132.3 millimeters. The southwestern part of the region is characterized by higher elevation, while the northeastern part is relatively lower. The total land area is approximately 105,500 square kilometers, with approximately 70% consisting of mountains and hills. Plains are predominantly distributed in the northern region. In 2020, the forest area in the study area was 60,800 square kilometers, accounting for 61.17% of the total land area in the province. The vegetation in the region is characterized by zonal vegetation, predominantly consisting of subtropical evergreen broad-leaved forests. The main broad-leaved tree species include Quercus L, Cinnamomum camphora, and Schima superba Gardn. et Champ. The main coniferous tree species include Pinus massoniana, Pinus elliottii, Pinus taiwanensis Hayata, and Cunninghamia lanceolate (Tao et al., 2022).

FIGURE 1

Figure 1 The distribution of forest resources in Zhejiang Province, China. (Note:Image cited from:http://bzdt.ch.mnr.gov.cn/browse.html?picId=%224o28b0625501ad13015501ad2bfc0290%22.).

2.2 Experimental design and data preprocessing

The data were collected from forest plots within the permanent plots of the NCFI in Zhejiang Province, China (Figure 2). A total of 4250 permanent plots were established across the entire province using a grid of 4 km × 6 km. The plots were square with an area of 800 square meters. Individual measurements were conducted for all trees within the plots with a diameter at breast height (DBH) equal to or greater than 5 centimeters (NFGA, 2020). The survey interval for NCFI was 5 years. The NCFI data selected for analysis were obtained from surveys conducted in 1989, 1994, 1999, 2004, 2009, 2014, and 2019, spanning a total period of 30 years. Subsequently, the entire set of seven datasets was divided into six periods: 1989-1994, 1994-1999, 1999-2004, 2004-2009, 2009-2014, and 2014-2019. The plots within the same period were paired. The plots were required to be forested at the beginning of the reexamination period (early stage).

FIGURE 2

Figure 2 Distribution of NCFI plots in Zhejiang Province, China.

In forest management practices in subtropical China, tree species are often categorized into three groups: pine, fir, and broadleaf. Such classification meets the requirements of general forest management practices. Classifying tree species is crucial for determining dominant tree species within a community (Kazmierczak et al., 2014). The pine group includes Pinus massoniana Lamb, Pinus thunbergii Parl, Pinus taiwanensis Hayata, Pinus elliottii Engelm, and others. Among them, Pinus massoniana Lamb (Masson pine) accounts for the vast majority. The fir group includes Cunninghamia lanceolata (Lamb.) Hook, Cryptomeria fortune Hooibrenk ex Otto et Dietr, Metasequoia glyptostroboides Hu & W. C. Cheng, Taxodium ascendens Brongn, and Taxus wallichiana var. chinensis (Pilg.) Florin, and others. Among them, Cunninghamia lanceolata (Lamb.) Hook (Chinese fir) constitutes the vast majority. The broad-leaved group includes Quercus, Castanopsis, Liquidambar, Schima superba Gardn. et Champ., Lauraceae (Sun et al., 2020), and others.

Table 1 displays the number of remeasured forest plots in each period based on dominant tree species group statistics. The tree species group that exhibited the highest proportion of biomass within a plot was considered the dominant tree species group. The forest plots in the first three periods were relatively fewer compared to the latter three periods because in the surveys conducted in 1994, 1999, and 2004, one-third of the plots were relocated and could not be paired. The tree biomass was calculated based on plot stem measurement data (Ji et al., 2012). The biomass proportion of each tree species group in each plot was calculated as an indicator of species group abundance, where the sum of abundance indices for the three tree species groups within each plot was equal to 1. The variation in actual tree species abundance across years is illustrated in Figure 3, revealing significant trends. The pine group exhibits a continuous decline, transitioning from absolute dominance to absolute inferiority, whereas the broad-leaved group exhibits the opposite pattern.

TABLE 1

Table 1 The numbers of plots in each reexamination period.

FIGURE 3

Figure 3 Changes in the average abundance of each species group (pine, fir, and broadleaf) during 1989-2019 based on the CFI data.

2.3 Nonlinear difference equation system model

In this study, a nonlinear difference equation system (NDES) model was developed, consisting of three difference equations (the model construction and derivation process are described in Appendix A). The NDES model is presented as follows:

$\begin{array}{ll}\{\begin{array}{c}{y}_{11}=\frac{1}{1+t_1{e}^{t_3{y}_{10}+t_4{y}_{20}}+t_2{e}^{t_5{y}_{10}+t_6{y}_{20}}}\equiv f_1\\ y_{21}=\frac{1}{(1/t_1)e^{-t_3{y}_{10}-t_4{y}_{20}}+1+(t_2/t_1)e^{(t_5-t_3)y_{10}+(t_6-t_4)y_{20}}}\equiv f_2\\ y_{31}=\frac{1}{(1/t_2)e^{-t_5{y}_{10}-t_6{y}_{20}}+(t_1/t_2)e^{(t_3-t_5)y_{10}+(t_4-t_6)y_{20}}+1}\equiv f_3\end{array}& (1)\end{array}$

where y_ij represents the species abundance, i denotes the species group, i=1 (pine), 2 (fir), 3 (broadleaf), and j represents the temporal sequence of the remeasured plots, j=0 (former), 1 (latter). Auto conform 0≤y_ij ≤ 1, y₁₁+y₂₁+y_31 = 1. t₁, t₂, t₃, t₄, t₅, and t₆ are 6 parameters of the NDES model.

Rewriting Eq. (1) in vector form:

$\begin{array}{ll}{Y}_1=F(Y_0)& (2)\end{array}$

where $Y_1={(y_{11},y_{21},y_{31})}^{\text{'}}$, $Y_0={(y_{10},y_{20},y_{30})}^{\text{'}}$, and $F={(f_1,f_2,f_3)}^{\text{'}}$. By substituting Y₀ into the right side of Eq. (2), the estimated value Y₁ for 5 years later (the remeasurement interval in this study is 5 years) can be obtained. Similarly, by substituting Y₁ into the right side of Eq. (2), the estimated value Y₂ for 10 years later can be obtained. This process continues iteratively, such that if Y_k(k=0,1,2,…) is already known, Y_k+1 can be calculated as Y_k+1=F(Y_k).

The estimation of the model parameters is based on plot-level data. For each plot, the species abundance (y_ij) is calculated, including both the early and late stages. First, the total biomass of the plot is computed, and the proportion of biomass for each species group to the total biomass represents the species abundance of that particular species group.

2.4 Prediction of species abundance limits, interval analysis and hypothesis

If $\underset{k\to \infty }{\lim }F(Y_k)=Y^{*}$, i.e., $Y^{*}=F(Y^{*})$, where $Y^{*}={(y_1^{*},y_2^{*},y_3^{*})}^{\text{'}}$is a constant vector, then $Y^{*}$ is referred to as the equilibrium point of the difference equation system, and such an equilibrium point is inherently stable (Elaydi, 2005). A model with a stable equilibrium point possesses significant biological significance, as it facilitates theoretical analysis and long-term predictions. In this study, the equilibrium point is referred to as the limit, representing the theoretical ultimate state of tree species abundance variations.

Reliability analysis is necessary for limits. This paper simulates the variance of limit $Y^{*}$ through Monte Carlo methods. Let $T={(t_1,\cdots ,t_6)}^{\text{'}}$ be the estimated value of the model parameter, with its covariance matrix denoted as ${\sum }^{}$. Let $T^0\sim N(T,{}^{\sum })$, and generate a set of random numbers $T^0={(t_1^0,\cdots ,t_6^0)}^{\text{'}}$for each plot. Taking the latter observed abundance of the plot as the initial value, we calculated the limit of the model by considering T⁰ as the model parameter. We estimate the covariance matrix $D(Y^{*})$by computing the covariance matrix ${\Sigma }_{3\times 3}^{*}$ based on the simulated limits of all n plots and consider it as an estimate for $D(Y^{*})$ . Since $y_1^{*}+y_2^{*}+y_3^{*}=1$ , the sum of the elements in matrix ${\Sigma }_{3\times 3}^{*}$ is zero.

Interval estimation of the abundance limit for a single tree species group. Assume that the abundance limit follows a normal distribution. Because the abundance limit is bounded between 0 and 1, it is theoretically appropriate to use a truncated normal distribution. Let $y^{*}$ be the point estimate of the abundance limit for a certain tree species group and ${\sigma }^2$ be the variance of $y^{*}$, represented by the diagonal elements in matrix

${\Sigma }_{3\times 3}^{*}$ . Let $p(y,y^{*},{\sigma }^2)$ denote the probability density function of the normal distribution. The probability density function of the truncated normal distribution is given by:

$\begin{array}{ll}{p}_{\text{tr}}(y,y^{*},{\sigma }^2)=p(y,y^{*},{\sigma }^2)/\tau ,\tau ={\displaystyle {\int }_0^{1}p(y,y^{*},{\sigma }^2)dy}& (3)\end{array}$

Overall hypothesis testing of the limit vector. The overall hypothesis testing of the limit vector for tree species abundance was conducted using a multivariate normal distribution. Let $Y^{*}={(y_1^{*},y_2^{*},y_3^{*})}^{\text{'}}$ be the point estimate of the limit vector and $Y^{0*}={(y_1^{0*},y_2^{0*},y_3^{0*})}^{\text{'}}$ be another vector. The aim is to assess the possibility of $Y^{*}=Y^{0*}$ . The sum of the three components of the abundance limit is 1, implying that it has 2 degrees of freedom (3-1 = 2). Consequently, the rank of matrix ${\Sigma }_{3\times 3}^{*}$ is also 2, indicating that the matrix is not full rank and thus cannot be inverted. Considering the abundance limits of the first two tree species groups, namely, pines and firs, the determination of the abundance limit for broadleaves was subsequently established. Therefore, we only need to examine the limits of the first two tree species groups.

Let $Y^{*(2)}={(y_1^{*},y_2^{*})}^{\text{'}}$, $Y^{0*(2)}={(y_1^{0*},y_2^{0*})}^{\text{'}}$ and ${\Sigma }_{2\times 2}^{*(2)}=\left(\begin{array}{cc}{\sigma }_1^2& {\sigma }_{1·2}\\ {\sigma }_{2·1}& {\sigma }_2^2\end{array}\right)$, and compute the F-statistic:

$\begin{array}{ll}F=\frac{(n-p)}{(n-1)p}{(Y^{*(2)}-Y^{0*(2)})}^{\text{'}}{({\Sigma }_{2\times 2}^{*(2)})}^{-1}(Y^{*(2)}-Y^{0*(2)})\sim F(p,n-p)& (4)\end{array}$

The F distribution has the first degree of freedom equal to p and the second degree of freedom equal to n-p, where p is 2. To test whether both pines and firs may converge to zero in the future, a hypothesis test can be conducted by setting $Y^{0*(2)}={(0,0)}^{\text{'}}$.

2.5 Long-term prediction of tree species abundance

To uniformly forecast tree species abundance until 2119, a comparison of the prediction curves from models of different periods was enabled. This analysis can then be used to examine the changes in tree species composition over the past 30 years within the modeling dataset. Based on the plot data and using the initial species abundance from the latter observations of the modeling dataset, iterative calculations can be performed to obtain predicted values of species abundance for every 5-year interval in the plots. The average of the predicted species abundance values for the plots was calculated as the predicted average value of species abundance for the entire province. A predicted value for a plot is dependent only on the previous adjacent state and is not influenced by other preceding states (Scherstjanoi et al., 2014).

3 Results

3.1 Model parameter estimation results and the actual prediction accuracy

The parameters of the nonlinear differential equation model in each reexamination period are shown in Table 2. The parameters had a certain regularity. t1 ranged from 0.722 to 1.5111, t2 ranged from 20.4263 to 42.8872, t3 ranged from -3.3836 to -2.3495, t4 ranged from 2.5992 to 3.5865, t5 ranged from -6.9519 to -5.4936, and t6 ranged from -4.4867 to -2.8888 (Table 2). The same parameter in different reexamination periods had the same trend label, i.e., positive or negative. The coefficients of determination were greater than 0.8.

TABLE 2

Table 2 The parameters of the nonlinear differential equation model in each reexamination period.

The predicted data are shown in Table 2. The abundance data in the table represent predicted values. The relative error is defined as the difference between the predicted values and actual values divided by the actual values. The actual values are shown in Figure 3. In general, as the prediction horizon increases, the error also increases. However, for the same prediction horizon, the models in the first two periods (1989-1994 and 1994-1999) exhibit significantly larger prediction errors compared to the models in the subsequent four periods.

3.2 Prediction of species abundance limit

To analyze the ultimate trend of the models in each period, the models were subjected to limit prediction. Figure 4 displays the predicted limit values of species abundance for each period model. It is evident that the timeline can be divided into two distinct phases, with 1999 as the dividing point: two periods prior to 1999 (1989-1994 and 1994-1999) and four periods following 1999 (1999-2004, 2004-2009, 2009-2014, and 2014-2019). The former phase exhibits significant fluctuations, whereas the latter phase shows minimal variations.

FIGURE 4

Figure 4 Estimation results of the limit value of the abundance of the pine, fir, and broad-leaved groups in each reexamination period (y-axis: the limit value of the abundance of the species groups).

3.3 Long-term prediction of species abundance

Based on the preceding limit analysis, it is evident that the trends of the models in the latter four periods are similar. In other words, over the past 20 years, the changes in species abundance have exhibited a consistent pattern. Now, a comparative analysis of long-term predictions is conducted based on the models from these four periods. Each period model predicts until the year 2119, with the shortest prediction spanning 100 years (2019-2119) and the longest prediction spanning 115 years (2004-2119) (Figure 5). Based on the predicted trends, there is a significant decrease in the abundance of pines, a decrease in firs, and a noticeable increase in broadleaves. This trend is consistent with the field CFI data. The trend lines of species abundance share a common characteristic: lines ① and ② are the furthest apart, while lines ③ and ④ are closest to each other and sandwiched between lines ① and ②. This indicates that the two models from the first 10 years exhibit greater fluctuations compared to the two models from the latter 10 years. The patterns of species abundance have become increasingly stable over the past 20 years.

FIGURE 5

Figure 5 The prediction of abundance of species groups in different reexamination periods (model 1: model based on the reexamination period, 1999-2004; model 2: model based on the reexamination period, 2004-2009; model 3: model based on the reexamination period, 2009-2014; model 4: model based on the reexamination period, 2014-2019).

3.4 Interval estimation and hypothesis testing

The interval estimation and hypothesis testing procedures for the limit values of each period model are identical. Take the most recent period (2014-2019) as an example. In our study, the predicted values of the species abundance of broad-leaved trees fluctuated significantly with the use of model one and model two, while the predicted curves fluctuated less between model three and model four. In addition, the difference in the predicted species abundance of broad-leaved trees was the lowest (gray line one) when model one was used. The number of plots, n, is 1712. The model parameter T can be found in the last column of Table 3. Table 4 presents the interval estimation results of species abundance for pines, firs, and broadleaves based on the 2014-2019 model obtained using the truncated normal distribution method.

TABLE 3

Table 3 Comparison of relative prediction error of abundance in different reexamination periods.

TABLE 4

Table 4 The interval estimation results of the abundance for a single tree species group.

The limit vector hypothesis testing results obtained using Eq. (4) are presented in Table 5. Hypothesis 1: $Y^{0*}$represents the trend of the abundance of pines approaching zero, while there is no change in the abundance of firs; Hypothesis 2: $Y^{0*}$represents the trend of the abundance of firs approaching zero, while there is no change in the abundance of pines; Hypothesis 3: $Y^{0*}$represents the trend of the species abundance of both pines and firs approaching zero. The results indicate that none of these three hypotheses can be established. This indicates that despite both pines and firs tending toward absolute disadvantage, their species abundance cannot approach zero (i.e., extinction).

TABLE 5

Table 5 Results of hypothesis testing of the abundance in three tree species groups.

4 Discussion

4.1 Through limit analysis, hidden patterns of variation can be revealed

The NDES model is developed based on the data of repeated measurements of permanent plots. Therefore, the variation trend in tree species abundance during the period of data collection determines the overall trajectory of the model. Under the assumption that the forest continues to develop according to the current trend, the model is expected to provide stable and reliable predictive performance. When the forest deviates from the current trend of development, the model’s predictions will deviate from the actual observations. This behavior is observed in nearly all models. A good model is capable of consistently forecasting this trend forward during predictions without experiencing significant biases solely due to extrapolation. From Table 2, it can be observed that for the same prediction timeframe, the models during the first two periods (1989-1994 and 1994-1999) exhibit significantly larger prediction errors compared to the models during the subsequent four periods. This indicates that there may be distinct differences in the patterns of species abundance variation between the first two periods and the latter four periods. Consequently, the models developed for the first two periods are not suitable for predicting the latter four periods. During the latter four periods, the forest exhibits a relatively consistent pattern of variation. As a result, the models developed for this phase all yield favorable predictive outcomes (refer to Table 2, Figure 5).

If the forest continues to develop steadily in a certain trend and a model can consistently reflect this developmental trend, then the theoretical limits of the model should reasonably reflect this trend. The following discussion addresses this issue from the perspective of limits. From Figure 3, it can be observed that the actual species abundance exhibits gradual changes over time. The abundance of pines shows an absolute dominance starting in 1989, reaching its peak value (0.5810) in 1994 and subsequently declining to gradually become an absolute inferiority. For firs, the abundance fluctuated upward from 1989 to 2004, reaching its peak (0.3203) in 2004 and then decreasing, transitioning from moderate dominance to absolute inferiority. The abundance of broadleaves maintained a stable growth trend and surpassed that of both pines and firs in 2009, becoming dominant and expanding its dominance thereafter. Despite significant changes over the 30-year period, there are no abrupt shifts observed in between.

The trend of abundance limits depicted in Figure 4 is distinct from the actual value trend shown in Figure 3. The limits represent the ultimate state of abundance, which differs from the current actual values and signifies the long-term trend of abundance. From Figure 4, it can be observed that the limit distribution can be divided into two distinct phases, with 1999 serving as the boundary. Before 1999, the limits of the models during the periods of 1989-1994 and 1994-1999 exhibit fluctuations. The limit for pines shows a significant decline but remains in clear dominance along with firs. The limit for broadleaves shows an upward trend, while the limit for firs demonstrates a relatively rapid increase. This indicates that during these two periods, there are distinct differences in the trends of species abundance, and a stable pattern of variation has not yet been established. After 1999, the limits of species abundance for the four periods of 1999-2004, 2004-2009, 2009-2014, and 2014-2019 demonstrated remarkable stability. Throughout these periods, broadleaves consistently maintain absolute dominance, while pines and firs consistently remain in a state of absolute inferiority. This indicates that over the course of the subsequent 20 years, although there have been significant changes in the actual species abundance, the trend of these changes has become stable and relatively consistent, remaining within the same overall pattern. It can be observed that through limit analysis using the NDES model, certain patterns of variation can be revealed that may not be readily apparent. The NDES model is capable of adequately simulating the trend of variation and maintaining this trend consistently during long-term predictions, as long as there are no significant changes in the trend of forest tree species abundance. Therefore, the risks associated with long-term prediction (extrapolation) are minimal.

The limit analysis reveals the ultimate state of tree species abundance variation, while the interval analysis of the limits indicates the potential range of variation for this ultimate state. Although Figure 3 indicates a decreasing trend in the actual abundance of pines and firs, the predictive figures in Figure 5 also demonstrate the same pattern. However, the interval analysis of the limits suggests that the probability of their abundance approaching zero is nearly zero. In other words, while pines and firs may perpetually remain in a state of absolute inferiority, they are highly unlikely to go extinct.

4.2 Comparison with existing methods

In terms of model classification, the NDES model proposed in this study is an empirical model, which fundamentally differs from the category of gap model and ecological process model. In the empirical modeling approach, the model proposed in this study possesses its own unique structure, which is entirely distinct from the models employed by Nikolova et al. (2019); Fort and Grigera (2021), and others.

In terms of the number of tree species (groups) that the model can simulate, the model proposed in this study is capable of simulating only three tree species (groups), which is significantly fewer compared to the models employed by Nikolova et al. (2019); Fort and Grigera (2021), and others. However, the model proposed in this study is designed to meet the practical needs of forest management in the subtropical forests of China. It is capable of fulfilling the application requirements in this specific context. From an ecological research perspective, the model’s results can provide a macroscopic understanding of the changes in tree species composition within a large-scale forest region. The number of difference equations in the model proposed in this study is equal to the number of tree species (groups) under investigation. It is evident that the model can be expanded according to the methods provided in Appendix A to meet the research requirements as needed.

The predicted tree species abundances from the NDES model are constrained to be nonnegative and cannot exceed 1. The sum of the abundances is guaranteed to be equal to 1. Such properties enhance the stability of the model. The limits of the NDES model effectively capture the trend of variation in the modeling data. Therefore, a key application of the NDES model is to predict potential future trends and ultimate states through limit analysis. It is a big characteristic of the proposed model. Currently, there is scarce literature reporting on such research.

4.3 Analysis of the driving forces behind changes in forest tree species abundance in Zhejiang Province

The primary forest plant community in Zhejiang Province is characterized by evergreen broad-leaved forests. However, prolonged anthropogenic activities have led to significant degradation of these forests in most areas, resulting in the emergence of substantial secondary forests, primarily dominated by Masson’s pine, and mixed coniferous and broad-leaved secondary forests. In recent decades, the province has experienced economic development, a reduced dependence on forest-based economies, and an increased awareness of forest conservation, which has facilitated rapid forest recovery. The decrease in coniferous tree populations is attributed to the prevalence of pine wilt disease and fluctuations in timber prices. This reduction in coniferous trees has created ecological opportunities for the growth of broad-leaved species. Broad-leaved trees typically exhibit greater shade tolerance compared to conifers, leading to an increase in forest density. Rising labor costs have resulted in reduced forest management intensity, further contributing to increased forest density. When forest density reaches a certain threshold, strong shade-intolerant tree species such as pines and firs are outcompeted by broad-leaved species, leading to increased mortality. Under these circumstances, even in the absence of other influencing factors, the proportion of coniferous trees is expected to decline. These factors have resulted in the forest ecosystem of Zhejiang Province undergoing a forest succession process that closely approximates a natural state. Therefore, over the past 20 years, the species composition of forests in Zhejiang Province has maintained a fundamentally similar trend and has become increasingly stable. Among the four predicted curves in Figure 5 for the most recent four periods, curve ① from 1999-2004 and curve ② from 2004-2009 are far apart, while curve ③ from 2009-2014 and curve ④ from 2014-2019 are closer to each other, indicating a diminishing fluctuation between the curves. These circumstances indicate that over the past few decades, the changes in Zhejiang Province’s forests have primarily been driven by economic development and an increased awareness of conservation. The spread of pine wilt disease and fluctuations in timber prices have accelerated these changes. Once these changes reach a certain threshold, natural succession becomes the dominant driving force.

5 Conclusions

Due to various reasons, pine and fir tree species in Zhejiang Province have been consistently decreasing, and this trend is expected to persist over the long term. The fir group shifted from a relative advantage to an absolute disadvantage. They will perpetually remain at an absolute disadvantage but will never go extinct. As the Climax community in subtropical forests is characterized by evergreen broad-leaved forests, the observed changes in forest tree species composition over the past few decades conform to the natural succession patterns of subtropical forests.

The nonlinear difference equation system (NDES) model proposed in this paper is suitable for simulating and predicting the forest tree species composition in the macroscopic region of subtropical China. It not only enables short-term, medium-term, and even long-term predictions but also allows in-depth exploration of currently obscure trends and their stability in forest tree species composition through limit analysis of the model. The application of the model can provide a basis for the development of sustainable forest management measures.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

Author contributions

BJ: Funding acquisition, Writing – original draft, Data curation, Investigation, Methodology. KY: Conceptualization, Formal Analysis, Writing – original draft. FW: Methodology, Validation, Writing – original draft. HG: Methodology, Validation, Writing – original draft. LJ: Funding acquisition, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was funded by The National Nature Science Foundation of China (32271876), the “Pioneer” and “Leading Goose” R&D Program of Zhejiang (2022C02053), and the Province-Academy Cooperative Forestry Science and Technology Project of Zhejiang Province and Chinese Academy of Forestry (2021SY05).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Blanco, J. A., Lo, Y.-H., Kimmins, J. P., Weber, A. (2021). Highlighting complex long-term succession pathways in mixed forests of the Pacific Northwest: A Markov chain modelling approach. Forests 12. doi: 10.3390/f12121770

CrossRef Full Text | Google Scholar

Bormann, F. H., Likens, G. E. (1979). Pattern and process in a forested ecosystem (New York, NY: Springer). doi: 10.1007/978-1-4612-6232-9

CrossRef Full Text | Google Scholar

Botkin, D. B., Janak, J. F., Wallis, J. R. (1972). Some ecological consequences of a computer model of forest growth. J. Ecol. 60, 849–872. doi: 10.2307/2258570

CrossRef Full Text | Google Scholar

Bruelheide, H., Böhnke, M., Both, S., Fang, T., Assmann, T., Baruffol, M., et al. (2011). Community assembly during secondary forest succession in a Chinese subtropical forest. Ecol. Monogr. 81, 25–41. doi: 10.1890/09-2172.1

CrossRef Full Text | Google Scholar

Bugmann, H. (2001). Review of forest gap models. Clim. Change 51, 259–305. doi: 10.1023/A:1012525626267

CrossRef Full Text | Google Scholar

Chen, J., Zhou, N., Li, C., Xu, Q. (2017). Application of meta-population competition mechasism in forest succession simulation of Mount Lushan. Chin. J. Ecol. 36, 862–868. doi: 10.13292/j.1000-4890.201703.003

CrossRef Full Text | Google Scholar

Chen, Y. (2021). Development of stand density logistic growth model based on difference equation. M.S. thesis. (In Chinese: Zhejiang A&F University).

Google Scholar

Greig, C., Robertson, C., Lacerda, A. E. B. (2018). Spectral-temporal modelling of bamboo-dominated forest succession in the Atlantic Forest of Southern Brazil. Ecol. Model. 384, 316–332. doi: 10.1016/j.ecolmodel.2018.06.028

CrossRef Full Text | Google Scholar

Cortés-Calderón, S., Mora, F., Arreola-Villa, F., Balvanera, P. (2021). Ecosystem services supply and interactions along secondary tropical dry forests succession. For. Ecol. Manage 482, 118858. doi: 10.1016/j.foreco.2020.118858

CrossRef Full Text | Google Scholar

De Jager, N. R., Van Appledorn, M., Fox, T. J., Rohweder, J. J., Guyon, L. J., Meier, A. R., et al. (2019). Spatially explicit modelling of floodplain forest succession: interactions among flood inundation, forest successional processes, and other disturbances in the Upper Mississippi River floodplain, USA. Ecol. Modell. 405, 15–32. doi: 10.1016/j.ecolmodel.2019.05.002

CrossRef Full Text | Google Scholar

Dulamsuren, C. (2021). Organic carbon stock losses by disturbance: Comparing broadleaved pioneer and late-successional conifer forests in Mongolia’s boreal forest. For. Ecol. Manage 499, 119636. doi: 10.1016/j.foreco.2021.119636

CrossRef Full Text | Google Scholar

Elaydi, S. (2005). An introduction to difference equations (New York, NY: Springer).

Google Scholar

Fort, H., Grigera, T. S. (2021). A method for predicting species trajectories tested with trees in barro colorado tropical forest. Ecol. Modell. 446, 109504. doi: 10.1016/j.ecolmodel.2021.109504

CrossRef Full Text | Google Scholar

Gilbertson, N. M., Kot, M. (2021). Dynamics of a discrete-time pioneer–climax model. Theor. Ecol. 14, 501–523. doi: 10.1007/s12080-021-00511-z

CrossRef Full Text | Google Scholar

González, C., Macip-Ríos, R., Suazo-Ortuño, I. (2022). Phylogenetic structure and diversity among herpetofaunal communities along a successional gradient of a tropical dry forest in Mexico. Perspect. Ecol. Conserv. 20, 249–255. doi: 10.1016/j.pecon.2022.05.004

CrossRef Full Text | Google Scholar

Hao, Z., Huang, J., Li, X., Sun, H., Fang, G. (2022). A multi-point aggregation trend of the outbreak of pine wilt disease in China over the past 20 years. For. Ecol. Manage. 505, 119890. doi: 10.1016/j.foreco.2021.119890

CrossRef Full Text | Google Scholar

He, Y., Zhang, Q., Zhang, M., Xu, G., Yang, Y. (2015). Spatial distribution patterns of Polar-birch secondary forest in different succession process stages in Changbai Mountain. J. Northwest For Univ 30, 8–13.

Google Scholar

Hu, S., Ma, J., Shugart, H. H., Yan, X. (2018). Evaluating the impacts of slope aspect on forest dynamic succession in Northwest China based on FAREAST model. Environ. Res. Lett. 13. doi: 10.1088/1748-9326/aaa7bd

CrossRef Full Text | Google Scholar

Ibrom, A., Jarvis, P. G., Clement, R., Morgenstern, K., Oltchev, A., Medlyn, B. E., et al. (2006). A comparative analysis of simulated and observed photosynthetic CO2 uptake in two coniferous forest canopies. Tree Physiol. 26, 845–864. doi: 10.1093/treephys/26.7.845

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, B., Tao, J., Zhang, G., Du, Q., Yao, H., Xu, J. (2012). Zhejiang province’s forest vegetation biomass assessment for guaranteed accuracy. J. Zhejiang A F Univ(In Chinese) 29, 328–334.

Google Scholar

Kazmierczak, M., Wiegand, T., Huth, A. (2014). A neutral vs. non-neutral parametrizations of a physiological forest gap model. Ecol. Modell. 288, 94–102. doi: 10.1016/j.ecolmodel.2014.05.002

CrossRef Full Text | Google Scholar

Keane, R. E., Loehman, R. A., Holsinger, L. M. (2011). “The FireBGCv2 landscape fire and succession model: a research simulation platform for exploring fire and vegetation dynamics,” in Fort collins, CO (U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station), 137. doi: 10.2737/RMRS-GTR-255

CrossRef Full Text | Google Scholar

Koven, C. D., Knox, R. G., Fisher, R. A., Chambers, J. Q., Christoffersen, B. O., Davies, S. J., et al. (2020). Benchmarking and parameter sensitivity of physiological and vegetation dynamics using the Functionally Assembled Terrestrial Ecosystem Simulator (FATES) at Barro Colorado Island, Panama. Biogeosciences 17, 3017–3044. doi: 10.5194/bg-17-3017-2020

CrossRef Full Text | Google Scholar

Lantz, V., McMonagle, G., Hennigar, C., Sharma, C., Withey, P., Ochuodho, T. (2022). Forest succession, management and the economy under a changing climate: Coupling economic and forest management models to assess impacts and adaptation options. For Policy Econ 142, 102781. doi: 10.1016/j.forpol.2022.102781

CrossRef Full Text | Google Scholar

Larocque, G. R., Bell, F. W. (2021). Evaluating the performance of a forest succession model to predict the long-term dynamics of tree species in mixed boreal forests using historical data in Northern Ontario, Canada. Forests 12 (9), 1181. doi: 10.3390/f12091181

CrossRef Full Text | Google Scholar

Li, W.-Q., Huang, Y.-X., Chen, F.-S., Liu, Y.-Q., Lin, X.-F., Zong, Y.-Y., et al. (2021). Mixing with broad-leaved trees shapes the rhizosphere soil fungal communities of coniferous tree species in subtropical forests. For. Ecol. Manage 480, 118664. doi: 10.1016/j.foreco.2020.118664

CrossRef Full Text | Google Scholar

Li, G., Huang, J., Guo, H., Du, S. (2020). Projecting species loss and turnover under climate change for 111 Chinese tree species. For. Ecol. Manage. 477, 118488. doi: 10.1016/j.foreco.2020.118488

CrossRef Full Text | Google Scholar

Lin, D., Lai, J., Muller-Landau, H. C., Mi, X., Ma, K. (2012). Topographic variation in aboveground biomass in a subtropical evergreen broad-leaved forest in China. PloS One 7, e48244. doi: 10.1371/journal.pone.0048244

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin-Benito, D., Molina-Valero, J. A., Pérez-Cruzado, C., Bigler, C., Bugmann, H. (2022). Development and long-term dynamics of old-growth beech-fir forests in the Pyrenees: Evidence from dendroecology and dynamic vegetation modelling. For. Ecol. Manage. 524, 120541. doi: 10.1016/j.foreco.2022.120541

CrossRef Full Text | Google Scholar

Mathieu, Decuyper., Roberto, O., Chávez, Madelon., Lohbeck, et al. (2022). Continuous monitoring of forest change dynamics with satellite time series. Remote Sens. Environ. 269, 112829. doi: 10.1016/j.rse.2021.112829

CrossRef Full Text | Google Scholar

Medlyn, B. E., Berbigier, P., Clement, R., Grelle, A., Loustau, D., Linder, S., et al. (2005). Carbon balance of coniferous forests growing in contrasting climates: Model-based analysis. Agric. For Meteorol 131, 97–124. doi: 10.1016/j.agrformet.2005.05.004

CrossRef Full Text | Google Scholar

Mi, X., Swenson, N. G., Jia, Q., Rao, M., Feng, G., Ren, H., et al. (2016). Stochastic assembly in a subtropical forest chronosequence: evidence from contrasting changes of species, phylogenetic and functional dissimilarity over succession. Sci. Rep. 6, 32596. doi: 10.1038/srep32596

PubMed Abstract | CrossRef Full Text | Google Scholar

Michalski, S. G., Bruelheide, H., Durka, W. (2017). Phylogenetic turnover during subtropical forest succession across environmental and phylogenetic scales. Ecol. Evol. 7, 11079–11091. doi: 10.1002/ece3.3564

PubMed Abstract | CrossRef Full Text | Google Scholar

Mina, M., Bugmann, H., Klopcic, M., Cailleret, M. (2015). Accurate modeling of harvesting is key for projecting future forest dynamics: a case study in the Slovenian mountains. Reg. Environ. Change 17, 49–64. doi: 10.1007/s10113-015-0902-2

CrossRef Full Text | Google Scholar

NFGA (2020). Technical regulations for continuous forest inventory (Beijing: Standards press of China).

Google Scholar

Nikolova, P. S., Rohner, B., Zell, J., Brang, P. (2019). Tree species dynamics in Swiss forests as affected by site, stand and management: A retrospective analysis. For. Ecol. Manage. 448, 278–293. doi: 10.1016/j.foreco.2019.06.012

CrossRef Full Text | Google Scholar

Norden, N., Angarita, H. A., Bongers, F., Martinez-Ramos, M., Granzow-de la Cerda, I., van Breugel, M., et al. (2015). Successional dynamics in Neotropical forests are as uncertain as they are predictable. Proc. Natl. Acad. Sci. U.S.A. 112, 8013–8018. doi: 10.1073/pnas.150040311

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouyang, S., Xiang, W., Wang, X., Zeng, Y., Lei, P., Deng, X., et al. (2016). Significant effects of biodiversity on forest biomass during the succession of subtropical forest in south China. For. Ecol. Manage. 372, 291–302. doi: 10.1016/j.foreco.2016.04.020

CrossRef Full Text | Google Scholar

Pacala, S. W., Canham, C. D., Silander, J. A., Jr. (1993). Forest models defined by field measurements: I. design northeastern For. simulator. Can. J. For. Res. 23, 1980–1988. doi: 10.1139/x93-249

CrossRef Full Text | Google Scholar

Pelissari, A. L., Figueiredo Filho, A., Péllico Netto, S., Ebling, A. A., Roveda, M., Sanquetta, C. R. (2017). Geostatistical modeling applied to spatiotemporal dynamics of successional tree species groups in a natural Mixed Tropical Forest. Ecol. Indic. 78, 1–7. doi: 10.1016/j.ecolind.2017.02.044

CrossRef Full Text | Google Scholar

Peng, S.-L., Hou, Y.-P., Chen, B.-M. (2010). Establishment of Markov successional model and its application for forest restoration reference in Southern China. Ecol. Modell. 221, 1317–1324. doi: 10.1016/j.ecolmodel.2010.01.016

CrossRef Full Text | Google Scholar

Perin, J., Pitchugin, M., Hébert, J., Brostaux, Y., Lejeune, P., Ligot, G. (2021). SIMREG, a tree-level distance-independent model to simulate forest dynamics and management from national forest inventory (NFI) data. Ecol. Modell. 440, 109382. doi: 10.1016/j.ecolmodel.2020.109382

CrossRef Full Text | Google Scholar

Piao, S., Yin, G., Tan, J., Cheng, L., Huang, M., Li, Y., et al. (2015). Detection and attribution of vegetation greening trend in China over the last 30 years. Glob Chang Biol. 21, 1601–1609. doi: 10.1111/gcb.12795

PubMed Abstract | CrossRef Full Text | Google Scholar

Risch, A. C., Schütz, M., Bugmann, H. (2009). Predicting long-term development of abandoned subalpine conifer forests in the Swiss National Park. Ecol. Modell. 220, 1578–1585. doi: 10.1016/j.ecolmodel.2009.04.007

CrossRef Full Text | Google Scholar

Rödig, E., Cuntz, M., Heinke, J., Rammig, A., Huth, A. (2017). Spatial heterogeneity of biomass and forest structure of the Amazon rain forest: Linking remote sensing, forest modelling and field inventory. Glob. Ecol. Biogeogr. 26, 1292–1302. doi: 10.1111/geb.12639

CrossRef Full Text | Google Scholar

Ruger, N., Condit, R., Dent, D. H., DeWalt, S. J., Hubbell, S. P., Lichstein, J. W., et al. (2020). Demographic trade-offs predict tropical forest dynamics. Science 368, 165–168. doi: 10.1126/science.aaz4797

PubMed Abstract | CrossRef Full Text | Google Scholar

Rutishauser, E., Herault, B., Baraloto, C., Blanc, L., Descroix, L., Sotta, E. D., et al. (2015). Rapid tree carbon stock recovery in managed Amazonian forests. Curr. Biol. 25, R787–R788. doi: 10.1016/j.cub.2015.07.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheller, R. M., Hua, D., Bolstad, P. V., Birdsey, R. A., Mladenoff, D. J. (2011). The effects of forest harvest intensity in combination with wind disturbance on carbon dynamics in Lake States Mesic Forests. Ecol. Modell. 222, 144–153. doi: 10.1016/j.ecolmodel.2010.09.009

CrossRef Full Text | Google Scholar

Scherstjanoi, M., Kaplan, J. O., Poulter, B., Lischke, H. (2014). Challenges in developing a computationally efficient plant physiological height-class-structured forest model. Ecol. Complex. 19, 96–110. doi: 10.1016/j.ecocom.2014.05.009

CrossRef Full Text | Google Scholar

Shugart, H., West, D. (1977). Development of an Appalachian deciduous forest succession model and its application to assessment of the impact of the chestnut blight. J. Environ. Manage. 5, 161–179.

Google Scholar

Sun, J., Guo, J., Shen, A., Xu, X., Feng, H., Zhang, S., et al. (2020). Composition and environmental interpretation of the communities of Sassafras tzumu, a protected species, at Zhejiang province in eastern China. Glob. Ecol. Conserv. 24, e01218. doi: 10.1016/j.gecco.2020.e01218

CrossRef Full Text | Google Scholar

Sun, Z., Ren, H., Schaefer, V., Lu, H., Wang, J., Li, L., et al. (2013). Quantifying ecological memory during forest succession: A case study from lower subtropical forest ecosystems in South China. Ecol. Indic. 34, 192–203. doi: 10.1016/j.ecolind.2013.05.010

CrossRef Full Text | Google Scholar

Tao, J., Zhang, G., Ji, B. (2022). Forest resources monitoring in Zhejiang province (Beijing, China: China Forestry Publishing House).

Google Scholar

Thuiller, W., Albert, C., Araújo, M. B., Berry, P. M., Cabeza, M., Guisan, A., et al. (2008). Predicting global change impacts on plant species’ distributions: Future challenges. Perspect. Plant Ecol. Evol. Syst. 9, 137–152. doi: 10.1016/j.ppees.2007.09.004

CrossRef Full Text | Google Scholar

Tian, H., Zhu, J., He, X., Chen, X., Jian, Z., Li, C., et al. (2022). Using machine learning algorithms to estimate stand volume growth of Larix and Quercus forests based on national-scale Forest Inventory data in China. For. Ecosyst. 9, 100037. doi: 10.1016/j.fecs.2022.100037

CrossRef Full Text | Google Scholar

Tomé, J., Tomé, M., Barreiro, S., Paulo, J. A. (2006). Age-independent difference equations for modelling tree and stand growth. Can. J. For. Res. 36, 1621–1630. doi: 10.1139/x06-065

CrossRef Full Text | Google Scholar

Wu, Z. (1980). Vegetation of China (Beijing, China: Science Press).

Google Scholar

Xu, M., Liu, T., Xie, P., Chen, H., Su, Z. (2020). Temporal changes in community structure over a 5-year successional stage in a subtropical forest. Forests 11 (4), 438. doi: 10.3390/f11040438

CrossRef Full Text | Google Scholar

Yan, X., Shugart, H. (2005). FAREAST: a forest gap model to simulate dynamics and patterns of eastern Eurasian forests. J. Biogeograpgy 32, 1641–1658. doi: 10.1111/j.1365-2699.2005.01293.x

CrossRef Full Text | Google Scholar

Yang, X.-D., Yan, E.-R., Chang, S. X., Da, L.-J., Wang, X.-H. (2014). Tree architecture varies with forest succession in evergreen broad-leaved forests in Eastern China. Trees 29, 43–57. doi: 10.1007/s00468-014-1054-6

CrossRef Full Text | Google Scholar

Yin, Y., Wu, Y., Bartell, S. M. (2009). A spatial simulation model for forest succession in the Upper Mississippi River floodplain. Ecol. Complex. 6, 494–502. doi: 10.1016/j.ecocom.2009.09.003

CrossRef Full Text | Google Scholar

Zhang, Q., Han, R., Zou, F. (2011). Effects of artificial afforestation and successional stage on a lowland forest bird community in southern China. For. Ecol. Manage. 261, 1738–1749. doi: 10.1016/j.foreco.2011.01.025

CrossRef Full Text | Google Scholar

Zhang, Y., Li, X., Kong, Z., Du, N., Wu, M. (2022). Subtropical forest vegetation development and climate change in Baishanzu area of Zhejiang Province, China, since the Holocene. Palaeogeogr. Palaeoclimatol. Palaeoecol. 608, 111293. doi: 10.1016/j.palaeo.2022.111293

CrossRef Full Text | Google Scholar

Zhang, H., Yang, Q., Zhou, D., Xu, W., Gao, J., Wang, Z. (2021). How evergreen and deciduous trees coexist during secondary forest succession: Insights into forest restoration mechanisms in Chinese subtropical forest. Glob. Ecol. Conserv. 25, e01418. doi: 10.1016/j.gecco.2020.e01418

CrossRef Full Text | Google Scholar

Zhu, X. F., Fang, X., Wang, L. F., Xiang, W. H., Alharbi, H. A., Lei, P. F., et al. (2021). Regulation of soil phosphorus availability and composition during forest succession in subtropics. Forest Ecol. Manag. 502, 119706. doi: 10.1016/j.foreco.2021.119706

CrossRef Full Text | Google Scholar

Appendix A. Derivation process of the nonlinear difference equation model for the abundance of the tree species group

The base model of the abundance of the tree species groups is set as the following exponential model:

$\begin{array}{ll}\{\begin{array}{c}{y}_{11}=a_{10}{e}^{a_{11}{y}_{10}+a_{12}{y}_{20}+a_{13}{y}_{30}}\\ y_{21}=a_{20}{e}^{a_{21}{y}_{10}+a_{22}{y}_{20}+a_{23}{y}_{30}}\\ y_{31}=a_{30}{e}^{a_{31}{y}_{10}+a_{32}{y}_{20}+a_{33}{y}_{30}}\end{array}& \text{A-Eq. (1)}\end{array}$

where y₁₀, y₂₀, and y₃₀ represent the early abundances of the pine, fir and broad-leaved groups, respectively; y₁₁, y₂₁, and y₃₁ represent the late stage; $y_{11}+y_{21}+y_{31}=1$; $y_{10}+y_{20}+y_{30}=1$; and 0≤y₁₁, y₂₁, y₃₁, y₁₀, y₂₀, y₃₀ ≤ 1 must be satisfied. $a_1{}_0,a_{11},a_{12},a_{13}$, $a_2{}_0,a_{21},a_{22},a_{23}$, $a_{30},a_{31},a_{32},and a_{33}$are the 12 parameters of the model. Since $y_{11}+y_{21}+y_{31}=1$, that is:

$a_{10}{e}^{a_{11}{y}_{10}+a_{12}{y}_{20}+a_{13}{y}_{30}}+a_{20}{e}^{a_{21}{y}_{10}+a_{22}{y}_{20}+a_{23}{y}_{30}}+a_{30}{e}^{a_{31}{y}_{10}+a_{32}{y}_{20}+a_{33}{y}_{30}}=1$ must be satisfied, and Equation A-Eq. (1) is adjusted to a proportional model:

$\begin{array}{ll}\{\begin{array}{l}{y}_{11}=\frac{a_{10}{e}^{a_{11}{y}_{10}+a_{12}{y}_{20}+a_{13}{y}_{30}}}{a_{10}{e}^{a_{11}{y}_{10}+a_{12}{y}_{20}+a_{13}{y}_{30}}+a_{20}{e}^{a_{21}{y}_{10}+a_{22}{y}_{20}+a_{23}{y}_{30}}+a_{30}{e}^{a_{31}{y}_{10}+a_{32}{y}_{20}+a_{33}{y}_{30}}}\hfill \\ y_{12}=\frac{a_{20}{e}^{a_{21}{y}_{10}+a_{22}{y}_{20}+a_{23}{y}_{30}}}{a_{10}{e}^{a_{11}{y}_{10}+a_{12}{y}_{20}+a_{13}{y}_{30}}+a_{20}{e}^{a_{21}{y}_{10}+a_{22}{y}_{20}+a_{23}{y}_{30}}+a_{30}{e}^{a_{31}{y}_{10}+a_{32}{y}_{20}+a_{33}{y}_{30}}}\hfill \\ y_{13}=\frac{a_{30}{e}^{a_{31}{y}_{10}+a_{32}{y}_{20}+a_{33}{y}_{30}}}{a_{10}{e}^{a_{11}{y}_{10}+a_{12}{y}_{20}+a_{13}{y}_{30}}+a_{20}{e}^{a_{21}{y}_{10}+a_{22}{y}_{20}+a_{23}{y}_{30}}+a_{30}{e}^{a_{31}{y}_{10}+a_{32}{y}_{20}+a_{33}{y}_{30}}}\hfill \end{array}& \text{A-Eq. (2)}\end{array}$

A-Eq. (2) is reduced to:

$\begin{array}{ll}\{\begin{array}{l}{y}_{11}=\frac{1}{1+(a_{20}/a_{10})e^{(a_{21}-a_{11})y_{10}+(a_{22}-a_{12})y_{20}+(a_{23}-a_{13})y_{30}}+(a_{30}/a_{10})e^{(a_{31}-a_{11})y_{10}+(a_{32}-a_{12})y_{20}+(a_{33}-a_{13})y_{30}}}\hfill \\ y_{21}=\frac{1}{(a_{10}/a_{20})e^{(a_{11}-a_{21})y_{10}+(a_{12}-a_{22})y_{20}+(a_{13}-a_{23})y_{30}}+1+(a_{30}/a_{20})e^{(a_{31}-a_{21})y_{10}+(a_{32}-a_{22})y_{20}+(a_{33}-a_{23})y_{30}}}\hfill \\ y_{31}=\frac{1}{(a_{10}/a_{30})e^{(a_{11}-a_{31})y_{10}+(a_{12}-a_{32})y_{20}+(a_{13}-a_{33})y_{30}}+(a_{20}/a_{30})e^{(a_{21}-a_{31})y_{10}+(a_{22}-a_{32})y_{20}+(a_{23}-a_{33})y_{30}}+1}\hfill \end{array}& \text{A-Eq. (3)}\end{array}$

Let

$\begin{array}{ll}\{\begin{array}{l}{k}_1=a_{20}/a_{10},k_2=a_{30}/a_{10},a_{30}/a_{20}=k_2/k_1\hfill \\ k_3=a_{21}-a_{11},k_4=a_{22}-a_{12},k_5=a_{23}-a_{13},\hfill \\ k_6=a_{31}-a_{11},k_7=a_{32}-a_{12},k_8=a_{33}-a_{13}\end{array}& \text{A-Eq. (4)}\end{array}$

According to Eq. (4), we can obtain:

$\begin{array}{ll}\{\begin{array}{l}{a}_{11}-a_{21}=-k_3,a_{12}-a_{22}=-k_4,a_{13}-a_{23}=-k_5\hfill \\ a_{31}-a_{21}=(a_{31}-a_{11})-(a_{21}-a_{11})=k_6-k_3\hfill \\ a_{32}-a_{22}=(a_{32}-a_{12})-(a_{22}-a_{12})=k_7-k_4\hfill \\ a_{33}-a_{23}=(a_{33}-a_{13})-(a_{23}-a_{13})=k_8-k_5\hfill \\ a_{11}-a_{31}=-k_6,a_{12}-a_{32}=-k_7,a_{13}-a_{33}=-k_8\hfill \\ a_{21}-a_{31}=k_3-k_6,a_{22}-a_{32}=k_4-k_7,a_{23}-a_{33}=k_5-k_8\hfill \end{array}& \text{A-Eq. (5)}\end{array}$

Then, Eq. (3) can be rewritten as:

$\begin{array}{ll}\{\begin{array}{l}{y}_{11}=\frac{1}{1+k_1{e}^{k_3{y}_{10}+k_4{y}_{20}+k_5{y}_{30}}+k_2{e}^{k_6{y}_{10}+k_7{y}_{20}+k_8{y}_{30}}}\hfill \\ y_{21}=\frac{1}{(1/k_1)e^{-k_3{y}_{10}-k_4{y}_{20}-k_5{y}_{30}}+1+(k_2/k_1)e^{(k_6-k_3)y_{10}+(k_7-k_4)y_{20}+(k_8-k_5)y_{30}}}\hfill \\ y_{31}=\frac{1}{(1/k_2)e^{-k_6{y}_{10}-k_7{y}_{20}-k_8{y}_{30}}+(k_1/k_2)e^{(k_3-k_6)y_{10}+(k_4-k_7)y_{20}+(k_5-k_8)y_{30}}+1}\hfill \end{array}& \text{A-Eq. (6)}\end{array}$

According to $y_{10}+y_{20}+y_{30}=1$, the elimination of $y_{30}$gives:

$\begin{array}{ll}\{\begin{array}{l}{y}_{11}=\frac{1}{1+k_1{e}^{k_5+(k_3-k_5)y_{10}+(k_4-k_5)y_{20}}+k_2{e}^{k_8+(k_6-k_8)y_{10}+(k_7-k_8)y_{20}}}\hfill \\ y_{21}=\frac{1}{(1/k_1)e^{-k_5+(k_5-k_3)y_{10}+(k_5-k_4)y_{20}}+1+(k_2/k_1)e^{(k_8-k_5)+(k_6-k_3-k_8+k_5)y_{10}+(k_7-k_4-k_8+k_5)y_{20}}}\hfill \\ y_{31}=\frac{1}{(1/k_2)e^{-k_8+(k_8-k_6)y_{10}+(k_8-k_7)y_{20}}+(k_1/k_2)e^{(k_5-k_8)+(k_3-k_6-k_5+k_8)y_{10}+(k_4-k_7-k_5+k_8)y_{20}}+1}\hfill \end{array}& \text{A-Eq. (7)}\end{array}$

Let $k_1{e}^{k_5}=t_1,k_2{e}^{k_8}=t_2,k_3-k_5=t_3,k_4-k_5=t_4,k_6-k_8=t_5,\text{ and }{k}_7-k_8=t_6$; then, the final result is:

$\begin{array}{ll}\{\begin{array}{l}{y}_{11}=\frac{1}{1+t_1{e}^{t_3{y}_{10}+t_4{y}_{20}}+t_2{e}^{t_5{y}_{10}+t_6{y}_{20}}}\hfill \\ y_{21}=\frac{1}{(1/t_1)e^{-t_3{y}_{10}-t_4{y}_{20}}+1+(t_2/t_1)e^{(t_5-t_3)y_{10}+(t_6-t_4)y_{20}}}\hfill \\ y_{31}=\frac{1}{(1/t_2)e^{-t_5{y}_{10}-t_6{y}_{20}}+(t_1/t_2)e^{(t_3-t_5)y_{10}+(t_4-t_6)y_{20}}+1}\hfill \end{array}& \text{A-Eq. (8)}\end{array}$

A-Eq. (8) is Eq. (1). The derivation shows that t₁ and t₂ must be greater than zero.