Evaluating individual-based tree mortality modeling with temporal observation data collected from a large forest plot

Evaluating individual-based tree mortality modeling with temporal observation data collected from a large forest plot

Forest Ecology and Management 450 (2019) 117496 Contents lists available at ScienceDirect Forest Ecology and Management journal homepage: www.elsevi...

2MB Sizes 1 Downloads 248 Views

Forest Ecology and Management 450 (2019) 117496

Contents lists available at ScienceDirect

Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco

Evaluating individual-based tree mortality modeling with temporal observation data collected from a large forest plot Yu Zhua, Zhili Liua,b, Guangze Jina,b, a b

T



Center for Ecological Research, Northeast Forestry University, Harbin 150040, China Key Laboratory of Sustainable Forest Ecosystem Management-Ministry of Education, Northeast Forestry University, Harbin 150040, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Forest dynamics plot (FDP) Individual-based model (IBM) Tree mortality Model evaluation Temperate forest

In forests, tree mortality strongly determines forest dynamics, creates gaps for recruitment and contributes to the coexistence of tree species. Individual-based models (IBMs) of tree mortality are key submodels in forest gap models and have been shown to drive simulated long-term forest dynamics. However, tree mortality IBMs remain poorly evaluated at the stand scale, particularly quantitatively, due to the lack of adequate tree mortality demographic data. Tree mortality dynamic data have often been absent from previous mortality IBMs, resulting in difficulty during model evaluation at the stand scale. The goals of this study were (1) to develop a spatiallyexplicit tree mortality IBM derived from the FORSKA forest gap model and (2) to evaluate its performance by both qualitative and quantitative fits based on a 10-year interval of tree mortality demographic data under a 9ha forest dynamics plot (FDP) in an old-growth temperate forest in Northeast China. The results showed that, for model qualitative evaluation, the observed dead trees were mainly distributed in the southwest and northeast parts, but the predicted dead trees were centrally clustered only in the southwest part of the FDP. For model quantitative evaluation, the overall model error (Error (%)) for all trees was −10.5%. However, the absolute values of model error for some species (i.e., Acer spp., Betula spp., Tilia spp., Fraxinus mandschurica, and Ulmus spp.) and some diameter at breast height (DBH) classes (i.e., 30–40, 40–50, 50–60, 60–70, and ≥70 cm) were greater than or equal to 20%. Our study highlights that the tree mortality IBM at the stand scale still requires further improvement to enhance the model performance and that the FDPs could provide crucial data support for initialization, parameterization and evaluation in the tree mortality IBM at the stand scale in forest dynamic modeling.

1. Introduction Tree mortality is a key driver of forest dynamics, and it shapes forest structure and successional trajectories, creates gaps for recruitment and contributes to the coexistence of tree species and biodiversity in forests (Franklin et al., 1987; Canham et al., 2001; Uriarte et al., 2012). Characterizing spatiotemporal patterns of tree mortality is critical for understanding forest community dynamics (Gonzalez-Akre et al., 2016). Tree mortality is also one of the most uncertain processes in regard to modeling forest response to climate change (Bugmann et al., 2019). One of the earliest works that used gap models to simulate forest responses to the climate change was that of Shao (1996). In traditional forest modeling, three dimensions (functional, structural and spatial complexity) have generally been separately addressed by different families of forest models (physiological, forest gap and forest landscape models) (Seidl et al., 2012). An individual-based model



(IBM) of tree mortality is a key submodel in dynamic global vegetation models (Sato et al., 2007), forest gap models (Bugmann, 2001) and forest landscape models (Seidl et al., 2012). Gap models (Bugmann, 2001) have been developed to study the structural and compositional dynamics of forest ecosystems as mediated by the environment, and thus mainly deal with structural complexity. Currently, gap models have grown from stand scale to continental scale and even global scale applications to help assess the potential consequences of climate change on natural forests (Shugart et al., 2018). Therefore, tree mortality is a critical process in forest gap formulation and further impacts forest structures and their successional trajectories. As a key module in forest gap models, tree mortality submodels have been shown to drive simulated long-term forest dynamics (Bugmann et al., 2019). Gap models commonly considers the maximum age as an independent variable of tree mortality (Shugart, 1984). Shao et al., (1994) considered tree mortality in two situations in their KOPIDE model: (1) senescence

Corresponding author at: 26 Hexing Road, Xiangfang District, Harbin 150040, China. E-mail address: [email protected] (G. Jin).

https://doi.org/10.1016/j.foreco.2019.117496 Received 29 April 2019; Received in revised form 21 July 2019; Accepted 22 July 2019 0378-1127/ © 2019 Elsevier B.V. All rights reserved.

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

under a 9-ha FDP in an old-growth mixed broadleaved-Korean pine (Pinus koraiensis) temperate forest in Northeast China.

(assumes that 1% of trees will live to maximum age) and (2) when respiration nearly equals assimilation due to critically reduced carbon gain. These two situations represent the first modification of tree mortality modeling with gap models. Bugmann (1996) considered tree mortality in two situations in their FORCLIM model: (1) base mortality (life table with constant mortality) and (2) stress-induced mortality (combination of absolute and relative thresholds of “slow growth”). Köhler and Huth (1998) considered tree mortality in three situations in their FORMIND model: (1) background mortality rate, (2) stem diameter-dependent mortality rate, and (3) diameter increment-dependent mortality rate. Moorcroft et al., (2001) considered tree mortality based on three process-based mortality algorithms in their Ecosystem Demography (ED) model: (1) carbon starvation, (2) hydraulic failure, and (3) phloem failure. However, although many studies have evaluated tree mortality IBMs at global and regional scales (Bond-Lamberty et al., 2015), the individual-based tree mortality models remain poorly evaluated at the stand scale, particularly quantitatively, due to the lack of adequate tree mortality demographic data. Tree coordinates, habitat heterogeneity information and mortality dynamics data were often absent or not sufficient in most of the previous tree mortality IBMs of gap models, resulting in difficulties during model initialization and evaluation at the stand scale. Model simulations have to begin at the bare forest floor or from a set of small sample plots (Botkin et al., 1972; Leemans and Prentice, 1987; Yan and Shugart, 2005); see Hülsmann et al., (2017). In recent decades, large forest dynamics plots (FDPs) have been increasingly used to study biodiversity-maintaining mechanisms and to understand the pattern of tree growth, mortality and recruitment as well as community dynamics (Harms et al., 2000; Uriarte et al., 2004; Wright et al., 2010; Wang et al., 2012). Core censuses of FDPs often have the following advantages (Anderson-Teixeira et al., 2015): (1) the plot size is very large, ranging from 2 to 120 ha, which makes model calibration and validation possible. (2) Every freestanding woody individual tree ≥1 cm diameter at breast height (DBH) is measured, and individual tree biomass is derived using an allometric relationship; consequently, individual tree level biomass can be estimated. (3) One could carry out spatially-explicit individual-based modeling and analyze neighborhood interactions that benefit from the mapping of all individual trees and fine scale topography (e.g., elevation, slope, aspect, and convexity) in the census. (4) The census typically is repeated every 5 years (i.e., survival, growth and recruitment), making it possible to characterize demographic dynamics. Therefore, we can combine the core census data and supplementary measurement (e.g., plant functional traits, climate and edaphic dataset) from large FDPs to parameterize and evaluate ecosystem and earth system models (e.g., an ecosystem demography model derived from forest gap models) (Moorcroft et al., 2001; Medvigy et al., 2009; Anderson-Teixeira et al., 2015). However, no evaluation attempt for tree mortality IBMs at the stand scale using the census data and measurement from large FDPs dataset have been reported. In addition, recent research findings of studies on species coexistence mechanisms (e.g., variation of neighborhood radii among species, habitat heterogeneity, and negative density dependence, NDD) (Janzen, 1970; Connell, 1971; Piao et al., 2013; Zhu et al., 2017) in forest communities can help improve and correct the structure of tree mortality IBMs. For example, some studies have proposed that neighbors belonging to species that are phylogenetically more closely related to a focal individual have a significantly negative effects on the survival of that individual (Webb et al., 2006; Paine et al., 2012; Zhu et al., 2015). These studies on biodiversity-maintaining mechanisms in forests are often carried out in FDPs. In this paper, we take advantage of such a dataset from a 9-ha FDP in China. The objectives in this study were (1) to develop a spatiallyexplicit tree mortality IBM derived from the FORSKA forest gap model and (2) to evaluate its performance by both qualitative and quantitative fits based on a 10-year interval of tree mortality demographic data

2. Materials and methods 2.1. A large forest dynamics plot 2.1.1. Study site The study site is located in the Heilongjiang Liangshui National Natural Reserve (47°10′50″N, 128°53′20″E), Xiaoxing’an Mountains, Northeast China (Appendix S1). The area has been spared from logging and other major human disturbances since 1952. The reserve was established in 1980 and became part of the China’s Man and the Biosphere Reserve Network in September 1997. It was promoted to a national nature reserve with the approval of the Chinese State Council in December 1997 to protect its old-growth mixed broadleaved-Korean pine forest ecosystem. The reserve is approximately 12, 133 ha and is characterized by rolling mountainous terrain with elevations ranging from 280 m to 707 m a.s.l. The soil is classified as dark brown forest soil. The mean annual precipitation is 676 mm, with 78% relative humidity and an annual evaporation of 805 mm. The mean annual temperature is −0.3 °C, with a minimum mean of −6.6 °C and a maximum mean of 7.5 °C. The study site is located in a temperate continental monsoon climate area with high winds in the spring. The core zone of this reserve has never been logged. This forest is typical in structure and composition of the mature forests in this region. 2.1.2. Data collection A 9-ha (300 m × 300 m) FDP was established within the core zone of the reserve where human influence is minimal. The Liangshui FDP is part of the Chinese Forest Biodiversity Monitoring Network (CForBio). The plot was divided into standard quadrats (10 m × 10 m, 900 total). In 2005, all the freestanding live trees in the plot ≥2 cm DBH (at 1.3 m) were mapped, measured, identified by species and tagged. Shortly afterwards, the canopy width, tree height and bole height of the trees ≥10 m tree height were also measured. Tree height and bole height were measured by a Vertex III instrument (Haglöf Sweden AB, Långsele, Sweden). The status of each tree, alive or dead, and DBH were recorded in the 2015 census. A phylogenetic tree was established using the Phylomatic software based on the angiosperm phylogeny group (APG) III backbone phylogeny (http://phylodiversity.net/phylomatic/) (Zanne et al., 2014). We created a 5-m scale raster digital elevation model (DEM) based on 1-m elevation contours in ArcGIS 10.1 (ESRI, Redlands, CA, USA). Then, we calculated hillshade at 5-m scale from the DEM. We used a time-domain reflectometry (TDR) probe (IMKO, Ettlingen, Germany) to measure soil moistures (i.e., volumetric moisture content, %) in 2011 (Shi et al., 2015). 2.2. A spatially-explicit individual-based model of tree mortality 2.2.1. Model structure We shaped the model description of a spatially-explicit individualbased model (IBM) of tree mortality derived from FORSKA gap model in the context of the ODD (overview, design concepts, details) Protocol for IBMs to enhance understanding and comparability (Grimm et al., 2006; Grimm et al., 2010). The FORSKA gap model was selected for shaping the tree mortality IBM because it has been widely parameterized and applied in the mixed broadleaved-Korean pine temperate forest, Northeast China (Shao, 1991; Ge, 1996; Sang and Li, 1998). (1) Purpose: Develop a spatially-explicit tree mortality IBM at the stand scale, mainly derived from the mortality module of classical forest gap models (i.e., FORSKA) and modified based on species coexistence mechanisms. We developed this IBM to simulate the tree mortality dynamics pattern at the stand scale. (2) Entities, state variables, and scales: The entities included individual tree and grid cell (i.e., subquadrat, 5-m scale, one quadrat can 2

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

Table 1 Explanation of major parameters of the spatially-explicit individual-based model (IBM) of tree mortality (Botkin et al., 1972; Leemans and Prentice, 1987; Leemans, 1991). Explanation Environmental parameters Growing season-average incident light intensity Light extinction coefficient Growing degree-days Drought index Species parameters Half saturation point Light compensation point Growth scaling factor Sapwood maintenance cost factor Intrinsic mortality rate Growth-dependent mortality rate Minimum growing degree-days Maximum growing degree-days Maximum aridity index Threshold value for index of vigor Neighborhood radius Maximum possible DBH

Symbol

Dimension

I0 k DEGD DI

μmol m−2 s−1

α c γ δ U0 U1 DEGDmin DEGDmax DImax θ R DBHmax

μmol m−1 s−1 μmol m−1 s−1 cm2 m−1 a−1 cm2 m−2 a−1

Table 2 Biological parameters of the tree mortality IBM.

°C d a−1

°C d a−1 °C d a−1

m cm

be divided into four subquadrats) in the tree mortality IBM. Each individual tree is characterized by the state variables or attributes: species, spatial location (x and y coordinates), DBH, tree height, bole height, canopy width, leaf area, and biomass. Trees are grouped in different plant functional types (PFTs, i.e., shallow or deep root species). Each grid cell is characterized by the state variables or attributes of hillshade and soil moisture content. The model simulates a forest mortality in annual time steps, and the simulation runs usually comprise many years (10 years in this study). One grid cell represents 25 m2 (i.e., grain) and the model stand comprises 9 ha (i.e., extent). (3) Parameterization: We used 19 focal tree species from the Liangshui FDP in the tree mortality IBM. We adopted the three biological parameter sets from other studies (Shao, 1991; Ge, 1996; Sang and Li, 1998) that were all conducted in the mixed broadleaved-Korean pine forests in China (Tables 1 and 2), and we hypothesized that the biological parameters are phylogenetically conserved; we thus approximated the parameters of other species that are not included in the above parameter sets by applying phylogenetic pairwise distance information. We adopted the site environmental parameter from Sang and Li (1998) and Sang et al., (1999), in which studies were also conducted in the mixed broadleaved-Korean pine forests (Tables 1 and 3). We hypothesized that all the freestanding live individual trees were healthy in 2005. We began the simulation in 2005 in the plot since the freestanding live individual trees were mapped in 2005. Please see Appendix S2 for details of the parameter estimation methods in the IBM of tree mortality. (4) Process overview and scheduling: The model includes the following main processes: growth, photosynthetic production, respiration, and tree mortality. The annual growth of each tree is calculated on the basis of the main phenomenological photosynthesis assimilation and respiration. Growth process equations are partly taken from the FORSKA model. Allometric functions relate tree height, biomass and leaf area. Tree mortality can occur either by intrinsic mortality, growthdependent mortality, or windthrows, which are frequent in the study area. Processes are scheduled in the following order: calculation of light climate in the forest, photosynthetic production of leaves, respiration, and tree mortality (i.e., intrinsic, growth-dependent, and extrinsic mortality). Fig. 1 shows the flowchart of the tree mortality IBM. (5) Design concepts: Emergence – All analyzed forest attributes, such as forest biomass, directly emerge from the interactions of the individual trees. Interaction – The processes growth and mortality of a tree directly depend on other individuals. Interactions among neighboring trees occur via competition for light and space (self-thinning). Sensing –

Species

c

α

γ

δ

DBHmax

U0

Pinus koraiensis Picea koraiensis Picea jezoensis Abies nephrolepis Tilia amurensis Tilia mandshurica Fraxinus mandschurica Betula costata Betula platyphylla Ulmus japonica Ulmus laciniata Juglans mandshurica Phellodendron amurense Quercus mongolica Populus davidiana Populus ussuriensis Acer mono Acer tegmentosum Acer ukurunduense

14.78 5.69 5.69 8.53 12.09 12 24 30 50 25.4 24 20 19 41.4 37.5 36 28 28 28

250 175 175 180 260 250 400 420 650 300 300 400 400 600 550 550 300 300 300

14.56 5.81 5.81 9.11 9.08 12.2 31 50.24 78.39 21.06 21.78 22.07 19.81 133.7 30.41 50.04 15.71 15.71 15.71

0.045 0.016 0.016 0.024 0.024 0.026 0.094 0.131 0.136 0.055 0.051 0.038 0.038 0.324 0.066 0.13 0.027 0.027 0.027

200 100 100 60 100 100 100 82 60 90 70 80 100 60 60 150 60 30 30

0.009 0.015 0.015 0.023 0.015 0.015 0.018 0.023 0.046 0.023 0.023 0.018 0.018 0.013 0.031 0.023 0.031 0.031 0.031

Species

U1

DEGDmin

DEGDmax

DImax

θ

Pinus koraiensis Picea koraiensis Picea jezoensis Abies nephrolepis Tilia amurensis Tilia mandshurica Fraxinus mandschurica Betula costata Betula platyphylla Ulmus japonica Ulmus laciniata Juglans mandshurica Phellodendron amurense Quercus mongolica Populus davidiana Populus ussuriensis Acer mono Acer tegmentosum Acer ukurunduense

0.04 0.04 0.04 0.04 0.06 0.06 0.06 0.12 0.12 0.06 0.06 0.06 0.06 0.06 0.12 0.12 0.06 0.06 0.06

1350 800 800 650 1500 1500 1600 1350 1350 1250 1250 1600 1650 1450 1400 1250 1450 1450 1450

3250 1800 1800 1750 3600 3600 3400 4000 4000 3600 3600 3800 3670 3800 4000 3600 3800 3800 3800

1.05 0.89 0.89 0.96 1.4 1.4 1.25 1.7 1.7 1.55 1.55 1.32 1.31 1.64 1.75 1.75 1.39 1.39 1.39

0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.03 0.04 0.025 0.025 0.025 0.025 0.025 0.04 0.04 0.025 0.025 0.025

Table 3 Environmental parameters of the tree mortality IBM. Symbol

I0

k

DEGD

DI

Parameter

450

0.4

1700

0.74

Growth of a tree is strongly influenced by the light climate, which is determined by the crown distribution of the local tree population. Mortality also depends on the density of the local tree population and the falling down of dying neighboring trees. Trees sense other individuals only indirectly via the local light climate or canopy crowding. Stochasticity – Mortality is modeled as a stochastic process. We conducted and reported the results of the replicated model runs for 50 simulations to calculate the mean predicted tree mortality rate. All facets of mortality were interpreted on the basis of probabilities. (6) Initialization: We began the simulation in 2005 in the 9-ha FDP since the freestanding live trees were mapped in 2005. The initial values of state variables of entities (i.e., tree and grid) were not arbitrarily chosen but were based on data from the FDP. (7) Time series input data: The model does not use input data to represent time-varying processes. Site conditions were assumed to be homogeneous, and there was no interannual variability of environmental conditions. (8) Submodels: The spatially-explicit tree mortality IBM is composed of corrected intrinsic mortality (U0′), corrected growth-dependent mortality (U1′), and extrinsic mortality (U2) as followed: 3

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

Fig. 1. Flowchart of the tree mortality IBM. The IBM is composed of intrinsic mortality, growth-dependent mortality, and extrinsic mortality. Growth-dependent mortality depends on tree growth efficiency. Extrinsic mortality depends on whether the species is shallow root or not.

X = U0′ + U1′ + U2

(1)

Q = 1 − Exp (−X )

(2)

Pn = 1 − (1 − Q1 )(1 − Q2 )(1 − Q3). ..(1 − Qi ). ..(1 − Qn )

(3)

species in some previous gap models (e.g., JABOWA, FORSKA, and FAREAST) (Botkin et al., 1972; Leemans and Prentice, 1987; Yan and Shugart, 2005). However, in previous studies there is a trade-off between growth rate and longevity, i.e., slow-growing trees (e.g., latesuccessional species) may become older than fast-growing trees (e.g., pioneer species) (Bigler and Veblen, 2009; Rötheli et al., 2011). Here, “slow/fast-growing” is a species-level characteristic. Therefore, we assumed that for any fast-growing trees of pioneer species (e.g., Betula and Populus spp.) with excessively slow growth under the growing pressure (e.g., competition effect from neighbors) in successive years, only 10% are expected to survive for 20 years. Therefore, the uncorrected annual growth-dependent mortality is 0.12 (Eqn 5 ). We also used growth efficiency (i.e., Erel) to determine if a tree was under the growing pressure at year i in the tree mortality IBM, just as in the FORSKA gap model. Growing pressure is related to growth efficiency (i.e., Erel), and the explanation of Erel is below (Eqn 6). For those trees of late-successional species (e.g., Pinus spp.) with excessively slow growth under the continuous growth pressure, only 10% are expected to survive for 60 years, and the uncorrected annual probability of growthdependent death is 0.04. For those trees of other species with excessively slow growth under the continuous growth pressure, only 10% are expected to survive for 40 years, and the uncorrected growth-dependent mortality rate is 0.06 in the IBM.

where X represents the mortality probability, Q represents the annual mortality probability, Qi represents the mortality probability at year i, and Pn represents the cumulative mortality probability within n years. (1) Intrinsic mortality (size-dependent mortality): Most previous models apply a constant annual intrinsic mortality probability throughout a tree’s lifespan, which leads to a higher cumulative mortality probability for older trees. However, the intrinsic mortality probabilities (U0) of saplings and older trees are often higher than those of adult trees (Liang, 2010). In addition, Manusch et al. (2012) suggested that it was preferable to simulate the intrinsic mortality in dynamic vegetation models (DVMs) by size rather than age (‘size-dependent’ intrinsic mortality rate). Therefore, we adopted a sine function (Liang, 2010) and used a size approach instead of an age approach to construct a new intrinsic mortality function (U0′) as follows:

DBH U0′ = U0·⎡1 − μ·sin ⎛ ·π ⎞⎤ ⎢ ⎝ DBHmax ⎠⎥ ⎣ ⎦ ⎜



(4)

where U0 represents the uncorrected annual intrinsic mortality probability from the FORSKA model, DBH represents the current DBH of the individual tree, DBHmax represents the maximum possible DBH of the species, and μ represents a corrected coefficient used to prevent U0′ from being zero, because even adult trees also have a certain degree of intrinsic mortality probability. U0′ represents the corrected annual intrinsic mortality probability. The coefficient μ is set to 0.25. (2) Growth-dependent mortality: The annual probabilities of growth-dependent mortality were assigned the same value for all

U1 = 1 − Exp ⎡ ⎣

ln (SurvProb) ⎤ SurvYrs ⎦

(5)

where U1 represents the uncorrected annual growth-dependent mortality probability, SurvProb represents survival probability, and SurvYrs represents survival years.

4

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

U1′ =

U1 1 + (Erel/ θ) ρ

comparison on the level of individual trees is overly rigorous. The model error was calculated by:

(6)

where U1′ represents the corrected annual growth-dependent mortality probability, Erel represents the growth efficiency (a detailed explanation is in the next paragraph), and ρ represents the steepness parameter, usually set to 999 according to the FORSKA model (Leemans, 1991), θ represents a species-specific threshold value for the index of vigor. The FORSKA model used growth efficiency (i.e., Erel) as a stress indicator to correct the ‘growth-dependent’ mortality rate. The growth efficiency of a tree is defined as the ratio of the realized growth rate and the maximum possible growth rate (Leemans and Prentice, 1987; Leemans, 1991). A tree that fails to maintain a certain minimum growth rate cannot survive for a long time, and the probability of tree mortality will increase. Leemans and Prentice (1987) regarded the species-specific certain minimum growth rate as a threshold value for the index of vigor (θ). We also adopted the concept of growth efficiency in the tree mortality IBM. The trees with excessively slow growth have higher mortality; and differences in predicted growth arise from the effects of topography and neighboring crowns on incident light, and competition from neighboring trees (which increases with neighbor biomass, proximity, and phylogenetic relatedness). A detailed description of the growth-dependent mortality was enclosed in the supplementary material (Appendix S3). (3) Extrinsic disturbance mortality: Our study site is located in a temperate continental monsoon climate region with high winds in the spring. The Heilongjiang Liangshui National Nature Reserve was heavily damaged by severe windstorms in 2008 and 2009, causing abundant windthrows (Zhen et al., 2013). Abies nephrolepis and Picea spp. are evergreen needle and late-successional species with shallow roots in the mixed broadleaved-Korean pine forest (Zhang, 2008). A. nephrolepis is also a hygrophilous species and grows in valleys (e.g., the northeast part of the plot) due to greater water availability (Zhang, 2008). Therefore, these trees are more vulnerable, are easily influenced by monsoons and tend to be uprooted and blown-down windthrown trees, especially when those trees are suppressed in the plot. The windthrow mortality may increase with soil moisture. A recent study conducted in a spruce fir valley forest (A. nephrolepis and Picea spp. are dominant species) in the same reserve (i.e., the Liangshui reserve) showed that soil moisture had a negative effect on sapling survival (Pu and Jin, 2018). Therefore, the tree mortality IBM considered extrinsic disturbance mortality (i.e., windthrow). We assumed that the trees will suffer annual windthrow mortality if the tree species belongs to the shallow root type (e.g., A. nephrolepis and Picea spp.).

−3·SM ⎞⎤ U2 = WMP0·⎡1 − Exp⎛ ⎢ ⎝ SM0 ⎠⎥ ⎣ ⎦ ⎜

Error (%) =

PredMort − ObsMort ·100% ObsMort

(8)

where PredMort and ObsMort represent the mean predicted (50 simulation) and observed tree mortality rates, respectively. In addition, we also evaluated the simulated tree growth based on DBH from 2005 to 2015 because the growth-dependent mortality (U1′) strongly depends on the simulated level of competition in the canopy; this mortality is a function of tree growth. We used the mean error (ME) to evaluate the performance of simulated tree growth: N

ME =

∑ i=1

PredGrowthi − ObsGrowthi N

(9)

where PredGrowthi and ObsGrowthi represent the predicted and observed tree growth for the ith evaluation, respectively, and N represents the number of trees. 2.2.3. Sensitivity analysis To assess the sensitivity of model parameters in the tree mortality IBM, the sensitivity of output variable (Y) to input parameters (X) (or the effect of parameter X on the variable Y), ΔY/ΔX, was calculated (Eqn 10) as a ratio of output variable change to parameter change (both in %) (Tatarinov and Cienciala, 2006). For the absolute quantity (|ΔY/ ΔX|), the parameters were ranked in terms of their effect on the modeled variable as (1) parameters with a high sensitivity (|ΔY/ΔX|≥0.5), (2) parameters with a medium sensitivity (0.3≤|ΔY/ΔX| < 0.5) and (3) parameters with a low sensitivity (|ΔY/ΔX| < 0.3). In the tree mortality IBM, the output variable (Y) was the predicted tree mortality rate and the input parameters (X) were selected as PD0, PDW0, SurvYrs (for U1′, i.e., growth-dependent mortality), WMP0, and SM0 (for U2′, i.e., windthrow mortality).

Sensitivity =

ΔY ΔX

(10)

where ΔY represents the output variable change (in %) and ΔX represents the parameter change (in %). We assigned six change ratios to the ΔX: 0.05, 0.10, 0.20, −0.05, −0.10, and −0.20. Similarly, we also calculated and compared the model errors for different settings (i.e., phylogenetic weighing vs. no weighting, fixed vs. variable neighbor radius, spatial homogeneity vs. heterogeneity, and cylinder crown vs. cone/ellipsoid of rotation architectures) to determine if the effects of these settings improved the model predictions. For example, to evaluate if the level of wind disturbance in the time window (i.e., 2005–2015) is realistic, we calculated and compared the model errors (i.e., Error (%)) for tree mortality rates with and without windthrow mortality in the tree mortality IBM.



(7)

where U2 represents annual extrinsic mortality probability, WMP0 represents maximum annual windthrow mortality probability (assume 0.03), SM represents soil moisture (i.e., volumetric moisture content, %), and SM0 represents the maximum soil moisture in the plot (assume 50). When SM is greater than this value, U2 will level off and be asymptotic to WMP0. Other algorithms involved DBH and the tree height function, the leaf area function and the growth efficiency in the tree mortality IBM, which adopted methods from the FORSKA model.

3. Results 3.1. Model qualitative evaluation We found that the rates of dead trees of the pioneer and companion species (i.e., Betula spp., Populus spp., Acer spp.) were relatively high based on the observed dead trees for main species (Table 4). Although the mean DBH of dead P. koraiensis trees was the largest among the species (Table 5), the proportion of dead trees was relatively low (Table 4). We found that the observed dead trees were mainly distributed in the southwest and northeast parts; however, the predicted dead trees were centrally clustered only in the southwest part of the FDP (Fig. 2). The dead trees were less distributed in the southeast part of the FDP, based on both observed and predicted mortality patterns (Fig. 2).

2.2.2. Model evaluation Model evaluation included qualitative and quantitative fits. For the qualitative fit, we compared the spatial patterns of predicted vs. observed tree mortality in 2015 in the Liangshui 9-ha FDP to determine if the simulated tree mortality appeared to match the spatial patterns from the FDP demographic dataset. For the quantitative fit, we compared the predicted vs. observed tree mortality rates of different species and DBH classes because a 5

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

3.3. Sensitivity analysis

Table 4 Observed vs. predicted tree mortality rates for the main species in the Liangshui 9-ha temperate forest dynamics plot, Northeast China. We combined the species with the same genus into a species group (i.e., Acer spp., Betula spp., Picea spp., Populus spp., Tilia spp., and Ulmus spp.) that had similar parameterization in the tree mortality IBM. We also dropped Juglans mandshurica, Phellodendron amurense, and Quercus mongolica because the abundances of these species are relatively small (N < 30). Species

Observed mortality rate (%)

Predicted mortality rate (%) (Mean)*

Error (%)**

Abies nephrolepis Acer spp. Betula spp. Fraxinus mandschurica Picea spp. Pinus koraiensis Populus spp. Tilia spp. Ulmus spp. All species

6.8 6.9 8.9 1.9 11.4 3.5 19.6 3.7 4.1 5.7

6.8 5.2 14.0 3.2 10.2 2.8 23.2 5.2 5.6 5.1

0.0 −24.6 +57.3 +68.4 −10.5 −20.0 +18.4 +40.5 +36.6 −10.5

For the parameter sensitivity analysis, there were two (2/6) high sensitivity grades found for parameters PD0 and SM0, and there were three (3/6) and four (4/6) high sensitivity grades found for parameters WMP0 and PDW0, respectively, for six change ratios (Table 7). However, there were six (6/6) high sensitivity grades found for parameter SurvYrs (Table 7). The absolute values of model errors for spatial homogeneity (−7%) and cylinder crown (−7%) settings in the tree mortality IBM were smaller than that for spatial heterogeneity and the cone/ellipsoid of rotation architectures used in the default model settings (-10.5%) (Table 8). However, the absolute values of model errors for no phylogenetic weighting (−12.2%) and fixed neighbor radius (+19.3%) settings in the tree mortality IBM were larger than that for phylogenetic weighting and variable neighbor radius used in the default model settings (−10.5%) (Table 8). The tree mortality IBM would underestimate the tree mortality rates when any mortality formulation (i.e., intrinsic, growth-dependent, and windthrow) was dropped, and the absolute values of model errors would be very large (Table 8).

* 50 simulations. Error (%) = (mean predicted mortality rate − observed mortality rate)/observed mortality rate·100%. ** +: The model is overestimating; −: underestimating.

4. Discussion 4.1. Simulated vs. observed tree mortality

Table 5 Mean and range of DBH for observed dead trees of the main species in the Liangshui 9-ha temperate forest dynamics plot, Northeast China. We combined the species with the same genus into a species group (i.e., Acer spp., Betula spp., Picea spp., Populus spp., Tilia spp., and Ulmus spp.) that had similar parameterization in the tree mortality IBM. We also dropped Juglans mandshurica, Phellodendron amurense, and Quercus mongolica because the abundances of these species are relatively small (N < 30). Species

Abies nephrolepis Acer spp. Betula spp. Fraxinus mandschurica Picea spp. Pinus koraiensis Populus spp. Tilia spp. Ulmus spp. All species

The dead trees were more distributed in the southwest part in the FDP based on both observed and predicted mortality patterns (Fig. 2), probably because this area had previously been mildly disturbed (Xu and Jin, 2012). Betula costata, Betula platyphylla and Populus davidiana are light-demanding pioneer species with the highest light compensation points (LCPs) among focal species in the mixed broadleavedKorean pine forest (Zhang, 2008). The intraspecific and interspecific competition for light and other resources among closely related species is often strong due to shared or similar resource requirements (Burns and Strauss, 2011), resulting in low tree growth efficiency and an intense negative effect on the survival of a focal tree (i.e., Conspecific NDD and Phylogenetic NDD mortality). Light-demanding trees often have shorter canopy heights than those of shade-tolerant trees, resulting in strong self-thinning. In addition, pioneer trees are often fastgrowing species with short lifespans and high intrinsic mortality probabilities (Manusch et al., 2012). Indeed, the rate of observed dead trees for pioneer species (i.e., Betula spp. and Populus spp.) were relatively high (Table 4). Most of the pioneer trees in 2005 were still young, with smaller DBH, but a small part of these pioneers were also already old and of larger size (Fig. 3). We found the observed dead trees were also more distributed in the northeast part according to the observed mortality pattern; however, this was not true for the predicted mortality pattern in the FDP (Fig. 2). The overall model error (i.e., Error (%)) for all trees was −10.5% based on the observed and predicted mortality rates (Tables 4 and 6). Acer tegmentosum and Acer ukurunduense are the dominant companion species with higher LCPs than P. koraiensis but lower LCPs than pioneer species in the mixed forest (Zhang, 2008). These two Acer species are also understory species and are distributed in clumps in the flat and mild northeast part of the plot (Xu and Jin, 2012). There are no trees that belong to these two species that live in the canopy or subcanopy positions, with relatively weak light conditions due to shading effects from overstory trees (Hao et al., 2009; Liu and Jin, 2010). In forests, competition for light highly depends on the spatial arrangement of trees and their canopies. A recent study conducted in American temperate forests showed that light is the most crucial resource for tree growth (McMahon et al., 2011). Similarly to B. costata, B. platyphylla and P. davidiana, intraspecific and interspecific competition for light and other resources among closely related species (i.e., Acer spp.) is often strong, which results in low tree growth efficiency and NDD mortality. Neighbors often have two kinds of NDD effects on a focal tree’s

DBH (cm) Mean

Min

Max

25.2 16.1 20.9 30.3 34.2 48.3 22.1 24.8 23.7 27.2

9.8 5.0 4.5 6.7 17.1 12.0 4.5 10.2 7.8 4.2

45.0 30.5 69.0 47.8 53.3 77.0 133.0 53.0 49.3 133.0

3.2. Model quantitative evaluation The overall model error (i.e., Error (%)) for all trees was −10.5% based on the observed and predicted mortality rates (Tables 4 and 6). The absolute values of model errors for some species (i.e., A. nephrolepis, Picea spp., P. koraiensis and Populus spp.) were less than or equal to 20%, but much larger errors were also found for the other species (i.e., Acer spp., Betula spp., Tilia spp., Fraxinus mandschurica, and Ulmus spp.) (Table 4). The absolute values of model errors for some smaller DBH classes (i.e., < 10, 10–20, and 20–30 cm) were less than 20% but much larger errors were also found for the larger DBH classes (i.e., 30–40, 40–50, 50–60, 60–70, and ≥70 cm) (Table 6). For the tree growth evaluation, the ME for all trees was −0.19 cm yr−1 (Appendixes S4-5). The absolute values of the MEs for some species (i.e., Acer spp., Betula spp., P. koraiensis, and Tilia spp.) were less than 0.20 cm yr−1, but much larger errors were also found for the others species (Appendix S4). The absolute values of the MEs for some smaller DBH classes (i.e., < 10 cm, 10–20 cm, 60–70 cm, and ≥70 cm) were less than 0.20 cm yr−1, but much larger MEs were also found for the medium and larger DBH classes (Appendix S5).

6

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

Fig. 2. Spatial pattern of observed vs. predicted (one simulation) tree mortality in 2015 in the Liangshui 9-ha temperate forest dynamics plot, Northeast China. (A) Observed mortality; (B) Predicted mortality.

probability will quickly increase when trees are close or even reach the maximum possible species-specific DBH, possibly due to the relatively small maximum possible DBH (30 cm) of A. tegmentosum and A. ukurunduense.

performance: the most crucial resource competition, and the light shading and extinction effect from the canopy of higher neighbors, emphasized by the Beer-Lambert law (Monsi and Saeki, 1953; Hirose, 2005; Monsi and Saeki, 2005). Moreover, the intrinsic mortality 7

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

Table 6 Observed vs. predicted tree mortality rates for different DBH classes in the Liangshui 9-ha temperate forest dynamics plot, Northeast China. DBH class (cm)

Observed mortality rate (%)

Predicted mortality rate (%) (Mean)*

Error (%)**

< 10 10–20 20–30 30–40 40–50 50–60 60–70 ≥70 All DBH classes

10.1 5.6 4.8 5.3 4.8 4.2 5.3 5.9 5.7

8.8 6.6 4.2 3.8 2.5 2.7 4.1 9.3 5.1

−12.9 +17.9 −12.5 −28.3 −47.9 −35.7 –22.6 +57.6 −10.5

Table 8 Observed vs. predicted tree mortality rates for different settings in the Liangshui 9-ha temperate forest dynamics plot, Northeast China.

* 50 simulations. Error (%) = (mean predicted mortality rate − observed mortality rate)/observed mortality rate·100%. ** +: The model is overestimating; −: underestimating.

PD0

PDW0

WMP0

SM0

SurvYrs

Default value

0.6 0.6 0.6 0.6 0.6 0.6 600 600 600 600 600 600 0.03 0.03 0.03 0.03 0.03 0.03 50 50 50 50 50 50 20** 60*** 40**** 20** 60*** 40**** 20** 60*** 40**** 20** 60*** 40**** 20** 60*** 40**** 20** 60*** 40****

ΔX

ΔY (Mean)*

Sensitivity

Sensitivity grade

0.05 0.10 0.20 −0.05 −0.10 −0.20 0.05 0.10 0.20 −0.05 −0.10 −0.20 0.05 0.10 0.20 −0.05 −0.10 −0.20 0.05 0.10 0.20 −0.05 −0.10 −0.20 0.05

0.0784 −0.1765 −0.0196 0.0196 0.0392 0.0392 0.0980 0.0588 0.0000 0.1373 0.0588 −0.0784 0.0392 0.0588 0. 0000 −0.0784 −0.0196 0.0392 −0.0392 0.0196 −0.0784 0.0980 −0.0392 −0.0784 −0.0980

1.5686 −1.7647 −0.0980 −0.3922 −0.3922 −0.1961 1.9608 0.5882 0. 0000 −2.7451 −0.5882 0.3922 0.7843 0.5882 0. 0000 1.5686 0.1961 −0.1961 −0.7843 0.1961 −0.3922 −1.9608 0.3922 0.3922 −1.9608

High High Low Medium Medium Low High High Low High High Medium High High Low High Low Low High Low Medium High Medium Medium High

0.10

−0.3137

−3.1373

High

0.20

−0.2745

−1.3725

High

−0.05

−0.0392

0.7843

High

−0.10

0.0784

−0.7843

High

−0.20

0.3725

−1.8627

High

Observed mortality rate (%)

Predicted mortality rate (%) (Mean)*

Error (%)**

Default setting No phylogenetic weighting Fixed neighbor radius Spatial homogeneity Cylinder crown No windthrow mortality No intrinsic mortality No growth dependent mortality

5.7 5.7

5.1 5.0

−10.5 −12.2

5.7 5.7 5.7 5.7 5.7 5.7

6.8 5.3 5.3 4.4 2.9 1.3

+19.3 −7.0 −7.0 –22.8 −49.1 −77.2

* 50 simulations. Default setting means the phylogenetic weighting, variable neighbor radius, spatial heterogeneity, cone/ellipsoid of rotation crown architectures, and three mortality formulations were used simultaneously in the IBM of tree mortality. Error (%) = (mean predicted mortality rate − observed mortality rate)/observed mortality rate·100%. ** +: The model is overestimating; −: underestimating.

Table 7 Sensitivity (expressed as ratio (% of variable change)/(% of parameter change)) of model output to quantitative parameters in the tree mortality IBM. Parameter

Model setting

koraiensis trees have shallow and developed lateral roots and prefer to live on hillsides or ridges with moist but well-drained soil (e.g., steep area). Adult trees of P. koraiensis are distributed in clumps in the steep habitat of the plot (Xu and Jin, 2012). Habitat heterogeneity and species-specific habitat preferences can result in spatial clustering and locally restricted distributions (Piao et al., 2013). Environmental filtering (Chesson, 2000; HilleRisLambers et al., 2012), arises from speciesspecific habitat preferences, or the inability of a species to persist in all local habitats, often occurring during earlier life stages (Hutchinson, 1957; Baldeck et al., 2017). Small trees tend to be more vulnerable to abiotic habitat impacts, and large trees perform well in their preferred habitats due to abiotic filtering that occurs for smaller size classes (Russo et al., 2005). A previous study in this plot showed that most of the focal species exhibited habitat preferences caused by large-scale habitat heterogeneity (Piao et al., 2013). This was consistent with the idea that environmental filtering at earlier life stages results in adult tree habitat associations (Comita and Engelbrecht, 2009). Adult trees of P. koraiensis (the proportion of DBH > 30 cm is greater than 80%) have experienced habitat filtering during earlier life-history stages and lived in their preferred habitat (i.e., steep areas). In addition, the intraspecific competition among these large canopy trees for light is relatively weak, and NDD mortality is unlikely. In contrast, physiological senescence (i.e., intrinsic mortality) may be the main reason for the mortality of adult trees of P. koraiensis in the FDP (Liu and Jin, 2010).

4.2. Structural, parametric and input data uncertainties in tree mortality IBMs For model quantitative evaluation, the overall model error (i.e., Error (%)) for all trees was −10.5%. However, the absolute values of model errors for some species (i.e., Acer spp., Betula spp., Tilia spp., F. mandschurica, and Ulmus spp.) or some DBH classes (i.e., 30–40, 40–50, 50–60, 60–70, and ≥70 cm) were greater than or equal to 20% (Tables 4 and 6). For the tree growth evaluation, the ME for all trees was −0.19 cm yr−1 (Appendixes S4-5). Consequently, the tree mortality IBM still requires further improvement to enhance the model performance. Tree mortality is one of the most uncertain processes, and growth-dependent and disturbance-induced tree mortality often have some uncertainties in the mortality submodel of forest models (Bugmann et al., 2019). In this study, the low performance obtained for the comparison between predicted and observed tree mortality is likely to result from uncertainties in structure, parameters, and input of the tree mortality IBM.

* 50 simulations. ** Pioneer species. *** Late-successsional species. **** Other species.

The dead trees were less distributed in the southeast part of the FDP from both observed and predicted mortality patterns (Fig. 2), probably because many P. koraiensis trees with large size grow in these areas. There are two steeper areas across the FDP: the east and the southeast parts of the FDP. Many P. koraiensis trees live in these areas and have large diameters, long lifespans and relatively low mortality probabilities (Xu and Jin, 2012). In contrast, many pioneer species with short lifespans and high light demands distributed in the southwest area (Xu and Jin, 2012) have relatively high mortality probabilities. P. 8

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

Fig. 3. Smooth density plot showing the DBH distribution for the main species in 2005 in the Liangshui 9-ha temperate forest dynamics plot, Northeast China. We combined the species with the same genus into a species group (i.e., Acer spp., Betula spp., Picea spp., Populus spp., Tilia spp., and Ulmus spp.) that had similar parameterization in the tree mortality IBM. We also dropped Juglans mandshurica, Phellodendron amurense, and Quercus mongolica because the abundances of these species are relatively small (N < 30). The vertical dashed lines represent the mean value of DBH.

9

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

Fig. 3. (continued)

weighting and variable neighbor radius in the competition index in the tree mortality IBM improved the model performance) (Table 8); however others showed that model structure improvements (e.g., spatial heterogeneity and the cone/ellipsoid of rotation crown architectures) failed to enhance the model performance (Table 8). Nonetheless, any tree mortality submodel (i.e., intrinsic, growth-dependent, and windthrow formulations) should be retained in the tree mortality IBM; otherwise, the tree mortality rates would be underestimated and the model errors would increase (Table 8). For example, the model error will rise to –22.8% if we neglected windthrow mortality in the FDP (Table 8). Parametric uncertainties: Parametric uncertainty is usually unavoidable due to the complexity and limited knowledge of the ecosystem structure and function (Song et al., 2013). Several parameter assumptions were used in the study. In the growth-dependent mortality formation, we assumed that the uncorrected annual growth-dependent mortality of any fast-growing trees of pioneer species, late-successional species and those trees of other species were 0.12, 0.04 and 0.06, respectively (i.e., SurvProb parameter), based on the assumptions that there is a trade-off between growth rate and longevity among species (Bigler and Veblen, 2009; Rötheli et al., 2011) and that these three species groups were expected to survive 20, 60, and 40 years, respectively (i.e., SurvYrs parameter), when the trees experienced excessively slow growth under continuous growth pressure (e.g., competition effect from neighbors). In addition, we also modified the neighbor

Structural uncertainties: The growth-related tree mortality formulations in forest models can be splitted into three broad categories: theoretical, empirical, and highly mechanistic formulations (Bugmann et al., 2019). The growth-related tree mortality formulations of the tree mortality IBM in this study is derived from the FORSKA gap model, whose growth-dependent mortality formulation are often based on theoretical reasoning (Bugmann et al., 2019). The simulation on forest structure and dynamics using theoretical mortality formulations have been shown to be highly sensitivity to the exact assumptions being made (Bugmann, 2001). It is almost impossible to determine a priori which theoretical mortality formulation in the forest models is correct (Bugmann et al., 2019). The empirically-based mortality formulations are often based on either tree-ring width (consequently, early-warning signals metrics) or forest inventory data (Mamet et al., 2015; Bugmann et al., 2019; Cailleret et al., 2019). However, the shortcoming is that different empirically-based tree mortality formulations have been shown to lead to widely different forest trajectories in long-term simulations (Bircher et al., 2015). The identification of a suitable mortality formulation in tree mortality submodel in forest models is often impossible when based on past data that typically cover only several decades (Bugmann et al., 2019). In this study, we found some structural uncertainties in the tree mortality IBM. We conducted comparisons between observed vs. predicted mortality rates for different model settings. Some results among these comparisons met our expectations (e.g., applying phylogenetic 10

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

mortality to improve the robustness of simulated future forest dynamics (Bugmann et al., 2019). The future development of tree mortality IBM may rely more on integrated approaches that fuse various data sources instead of formulations that are from spatiotemporally-limited, single datasets (Bugmann et al., 2019). Large FDPs could be an important data set source to support the development of the tree mortality IBM at the stand scale: (1) Model initialization: Although there are many gap models that start from bare ground (i.e., spin-up), the development of forest surveys (e.g., large FDPs) to provide the starting points (initialization) for simulations still needs improvement for forest IBMs (Shugart et al., 2018). In this study, the initial values of the state variables of entities (i.e., individual tree and grid) were not arbitrarily chosen but were based on measured data from the FDP. (2) Model parameterization: The combination of core census data (e.g., tree mortality dynamics and mapping of all stems and fine-scale topography) and supplementary measurements (e.g., plant functional traits, ecosystem measurements, soils data, and weather data) from large FDPs are invaluable for parameterizing ecosystem and earth system models (Anderson-Teixeira et al., 2015), and most of these kind of models are individual-based (e.g., the ED/ED2 model is derived from forest gap models) (Moorcroft et al., 2001; Medvigy et al., 2009). (3) Model evaluation: Previous studies have evaluated tree mortality IBMs, but mostly at the global and regional scales, and no evaluation attempts for this type of IBM at the stand scale have been reported using large FDP datasets. For example, the evaluation of tree mortality in the ED model is an indirect evaluation at regional/global scales via the way biomass declines (Frasson et al., 2015). In this study, we used tree mortality dynamics data from the FDP to directly evaluate tree mortality through the combination of qualitative (i.e., simulated vs. observed spatial patterns of mortality) and quantitative (i.e., Error (%)) methods at the stand scale. Recent study has shown that the forest gap model has been applied to the community assembly in forests (large FDPs are ideal platforms for the study of community assembly mechanisms) (Chauvet et al., 2017); however, the gap model applied to the study of the FDPs is rarely reported. In addition, accounting for biodiversity in ecosystem models will be important for predicting forest responses to climate change (AndersonTeixeira et al., 2015), and the next-generation IBMs must integrate biodiversity and ecosystems (Grimm et al., 2016). Large FDPs are exactly crosscutting research platform for biodiversity science (Ma, 2017). In the Barro Colorado Island (BCI) 50-ha FDP in Panama, censuses have been carried out for eight time periods (1981–1983, 1985, 1990, 1995, 2000, 2005, 2010, and 2015) (http://ctfs.si.edu/webatlas/datasets/ bci/), and the long-term time-series recensus dataset has been used to quantitatively evaluate an individual-based forest dynamics model (Ngugi and Botkin, 2011). Future work may extend the evaluations of tree mortality IBMs by including long term simulations and by assessing the effects of parameters and/or process formulations. We can expect that large FDPs will play a more critical role in research on IBMs in the future.

competition equation (Appendix S3 and S6) by applying the phylogenetic weighting function between a focal individual and its neighbors to distinguish the difference in competition intensity between intraspecific and interspecific neighbors (e.g., PNDD) (Piao et al., 2013; Zhu et al., 2015). We assumed that the minimum phylogenetic weight (i.e., PD0 parameter) and the maximum phylogenetic distance (i.e., PDW0 parameter) were 0.6 and 600 Ma (based on a phylogenetic tree used in the study), respectively. Moreover, we considered the windthrow mortality (Eqn 7, Appendix S6) as an extrinsic disturbance mortality probability in the tree mortality IBM due to high winds in the spring at the study site. We assumed that the maximum annual windthrow mortality probability (i.e., WMP0 parameter) and the maximum soil moisture (i.e., SM0 parameter) were 0.03 and 50 (based on soil volumetric moisture content data used in the study), respectively. In this study, we found that the parameter SurvYrs (survival years when trees with excessively slow growth under growth pressure, e.g., light competition in the canopy in successive years), was likely to be the most sensitive parameter, with six (6/6) high sensitivity grades, among five parameters (i.e., PD0, PDW0, WMP0, SM0, and SurvYrs) in the tree mortality IBM according to the sensitivity analysis (Table 7). Input data uncertainties: One of the uncertainty sources that has not been extensively explored was model inputs derived from empirical data (Fuentes et al., 2006). We found that the simulated tree growth in the tree mortality IBM performed poorly (the ME for all trees was −0.19 cm yr−1, Appendices S4-5). The simulated tree growth (and consequent simulated growth-dependent mortality) seems to not be realistic. The growth-dependent mortality strongly depends on the simulated level of competition in the canopy. Consequently, the source of input uncertainties in the tree mortality IBM considered is likely to be the canopy-related tree variables (e.g., tree height, bole height, crown width, and leaf area). The tree growth module in the tree mortality IBM was sensitive to these canopy-related tree variables (Eqn 12 in Appendix S3) (Leemans and Prentice, 1987). However, the measurement instrument (Vertex III) used in the traditional field survey in the study inevitably had certain measurement errors for these canopy-related tree variables. Recent studies showed that the unmanned aerial vehicle-light detecting and ranging (UAV-LiDAR) system was able to accurately measure the individual tree metrics (e.g., tree height and crown width) (Wallace et al., 2012), and we expect that the UAV-LiDAR system will help improve the measurement accuracy and consequently reduce the input uncertainties in the tree mortality IBM in the future work. In addition, the future development of the tree mortality IBM should consider the following limitations in the current model: (1) The use of an annual drought index in the tree mortality IBM can be further improved, especially for deciduous species; drought during the growing season is biologically more meaningful. In addition, the tree mortality IBM only uses climate averages but drought is more related to variation in water supply than average levels of precipitation, and interannual variation (i.e., a climate time series) among climate data is thus even more important. (2) Setting maximum and minimum DEGD based on species range map assumes that the range boundary is determined by the levels where survival is zero (Botkin et al., 1972); this assumption may be unlikely because the range boundary may be determined by any number of other limiting factors (e.g., competition, recruitment, growth, etc.), or it may be determined by a decreasing but nonzero survival rate.

5. Conclusions We developed a spatially-explicit tree mortality IBM derived from the FORSKA forest gap model and evaluate its performance by both qualitative and quantitative fits based on a 10-year interval of tree mortality demographic data under a 9-ha FDP in an old-growth temperate forest in Northeast China. The results showed that, for model qualitative evaluation, the observed dead trees were mainly distributed in the southwest and northeast parts, however the predicted dead trees were centrally clustered only in the southwest part in the FDP. For model quantitative evaluation, the overall model error (i.e., Error (%)) for all trees was −10.5%. However, the absolute values of model errors for some species (i.e., Acer spp., Betula spp., Tilia spp., F. mandschurica, and Ulmus spp.) and some DBH classes (i.e., 30–40, 40–50, 50–60, 60–70, and ≥70 cm) were greater than or equal to 20%. Our study

4.3. Large forest dynamics plots provide data set support for individualbased mortality modeling Large FDPs may provide important data set support for initialization, parameterization and evaluation in the spatially-explicit individual-based modeling of tree mortality at the stand scale in forest dynamic modeling. The development of the tree mortality IBM (e.g., model structure, initialization, parameterization and model evaluation, etc.) needs more data and a better process understanding of tree 11

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

highlights that the tree mortality IBM at the stand scale still requires further improvement to enhance model performance and that FDPs could provide crucial data support for initialization, parameterization and evaluation in the tree mortality IBM at the stand scale in forest dynamics modeling.

species. Can. J. For. Res. 31, 1–10. Chauvet, M., Kunstler, G., Roy, J., Morin, X., Russo, S., 2017. Using a forest dynamics model to link community assembly processes and traits structure. Funct. Ecol. 31, 1452–1461. Chesson, P., 2000. Mechanisms of maintenance of species diversity. Annu. Rev. Ecol. Syst. 31, 343–366. Comita, L.S., Engelbrecht, B.M., 2009. Seasonal and spatial variation in water availability drive habitat associations in a tropical forest. Ecology 90, 2755–2765. Connell, J.H., 1971. On the role of natural enemies in preventing competitive exclusion in some marine animals and in rain forest trees. Centre for Agricultural Publication and Documentation, Wageningen, Netherlands. Franklin, J.F., Shugart, H., Harmon, M.E., 1987. Tree death as an ecological process. Bioscience 37, 550–556. Frasson, R.P.D.M., Bohrer, G., Medvigy, D., Matheny, A.M., Morin, T.H., Vogel, C.S., Gough, C.M., Maurer, K.D., Curtis, P.S., 2015. Modeling forest carbon cycle response to tree mortality: Effects of plant functional type and disturbance intensity. J. Geophys. Res. Biogeosci. 120, 2178–2193. Fuentes, M., Kittel, T.G.F., Nychka, D., 2006. Sensitivity of ecological models to their climate drivers: statistical ensembles for forcing. Ecol. Appl. 16 (1), 99–116. Ge, J., 1996. Forest ecological modeling and simulation. Northeast Forestry University Press, Harbin, China. Gonzalez-Akre, E., Meakem, V., Eng, C.Y., Tepley, A.J., Bourg, N.A., McShea, W., Davies, S.J., Anderson-Teixeira, K., 2016. Patterns of tree mortality in a temperate deciduous forest derived from a large forest dynamics plot. Ecosphere 7, 17. Grimm, V., Ayllón, D., Railsback, S.F., 2016. Next-generation individual-based models integrate biodiversity and ecosystems: yes we can, and yes we must. Ecosystems 20, 229–236. Grimm, V., Berger, U., Bastiansen, F., Eliassen, S., Ginot, V., Giske, J., Goss-Custard, J., Grand, T., Heinz, S.K., Huse, G., Huth, A., Jepsen, J.U., Jørgensen, C., Mooij, W.M., Müller, B., Pe’er, G., Piou, C., Railsback, S.F., Robbins, A.M., Robbins, M.M., Rossmanith, E., Rüger, N., Strand, E., Souissi, S., Stillman, R.A., Vabø, R., Visser, U., DeAngelis, D.L., 2006. A standard protocol for describing individual-based and agentbased models. Ecol. Model. 198, 115–126. Grimm, V., Berger, U., DeAngelis, D.L., Polhill, J.G., Giske, J., Railsback, S.F., 2010. The ODD protocol: A review and first update. Ecol. Model. 221, 2760–2768. Hao, Z.Q., Salame, R., Lascoux, N., Salmon, E., Maioli, P., Kasparian, J., Wolf, J.P., 2009. Multiple filamentation of non-uniformly focused ultrashort laser pulses. Appl. Phys. B 94, 243–247. Harms, K.E., Wright, S.J., Calderón, O., Hernández, A., Herre, E.A., 2000. Pervasive density-dependent recruitment enhances seedling diversity in a tropical forest. Nature 404, 493–495. HilleRisLambers, J., Adler, P., Harpole, W., Levine, J., Mayfield, M., 2012. Rethinking community assembly through the lens of coexistence theory. Annu. Rev. Ecol. Evol. Syst. 43, 227–248. Hirose, T., 2005. Development of the Monsi-Saeki theory on canopy structure and function. Ann. Bot. 95, 483–494. Hülsmann, L., Bugmann, H., Brang, P., 2017. How to predict tree death from inventory data-lessons from a systematic assessment of European tree mortality models. Can. J. For. Res. 47, 890–900. Hutchinson, G.E., 1957. Concluding remarks. Cold Spring Harb. Symp. Quant. Biol. 22, 415–427. Janzen, D.H., 1970. Herbivores and the number of tree species in tropical forests. Am. Nat. 104, 501–528. Köhler, P., Huth, A., 1998. The effects of tree species grouping in tropical rainforest modelling: simulations with the individual-based model FORMIND. Ecol. Model. 109, 301–321. Leemans, R., 1991. Sensitivity analysis of a forest succession model. Ecol. Model. 53, 247–262. Leemans, R., Prentice, I.C., 1987. Description and simulation of tree-layer composition and size distributions in a primaeval Picea-Pinus forest. Vegetatio 69, 147–156. Liang, Y., 2010. A simulation system of disturbance effects on structure and process of secondary forest in Northeast China. In. Northeast Forestry University, Harbin, China. Liu, Y., Jin, G., 2010. Character of coarse woody debris in a mixed broad-leaved Korean pine forest in Xiaoxing' an Mountains, China. Scientia Silvae Sinicae 46, 8–14. Ma, K., 2017. Forest dynamics plot is a crosscutting research platform for biodiversity science. Biodivers. Sci. 25, 227–228. Mamet, S.D., Chun, K.P., Metsaranta, J.M., Barr, A.G., Johnstone, J.F., 2015. Tree rings provide early warning signals of jack pine mortality across a moisture gradient in the southern boreal forest. Environ. Res. Lett. 10, 084021. Manusch, C., Bugmann, H., Heiri, C., Wolf, A., 2012. Tree mortality in dynamic vegetation models-A key feature for accurately simulating forest properties. Ecol. Model. 243, 101–111. McMahon, S.M., Metcalf, C.J.E., Woodall, C.W., 2011. High-dimensional coexistence of temperate tree species: functional traits, demographic rates, life history stages, and their physical context. PLoS ONE 6, e16253. Medvigy, D., Wofsy, S.C., Munger, J.W., Hollinger, D.Y., Moorcroft, P.R., 2009. Mechanistic scaling of ecosystem function and dynamics in space and time: Ecosystem Demography model version 2. J. Geophys. Res. 114, G01002. Monsi, M., Saeki, T., 1953. Über den Lichtfaktor in den Pflanzengesellschaften und seine Bedeutung für die Stoffproduktion (On the factor light in plant communities and its importance for matter production). Jpn. J. Bot. 14, 22–52. Monsi, M., Saeki, T., 2005. On the factor light in plant communities and its importance for matter production. 1953. Ann. Bot. 95, 549–567. Moorcroft, P.R., Hurtt, G.C., Pacala, S.W., 2001. A method for scaling vegetation dynamics: the ecosystem demography model (ED). Ecol. Monogr. 71, 557–586. Ngugi, M.R., Botkin, D.B., 2011. Validation of a multispecies forest dynamics model using

Acknowledgements This study was financially supported by the National Natural Science Foundation of China (31870399), and the Fundamental Research Funds for the Central Universities (2572019CP14). We thank Dr. Harald Bugmann for helpful suggestions. We thank the editor and two anonymous reviewers for their constructive comments. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.foreco.2019.117496. References Anderson-Teixeira, K.J., Davies, S.J., Bennett, A.C., Gonzalez-Akre, E.B., Muller-Landau, H.C., Wright, S.J., Abu Salim, K., Almeyda Zambrano, A.M., Alonso, A., Baltzer, J.L., Basset, Y., Bourg, N.A., Broadbent, E.N., Brockelman, W.Y., Bunyavejchewin, S., Burslem, D.F., Butt, N., Cao, M., Cardenas, D., Chuyong, G.B., Clay, K., Cordell, S., Dattaraja, H.S., Deng, X., Detto, M., Du, X., Duque, A., Erikson, D.L., Ewango, C.E., Fischer, G.A., Fletcher, C., Foster, R.B., Giardina, C.P., Gilbert, G.S., Gunatilleke, N., Gunatilleke, S., Hao, Z., Hargrove, W.W., Hart, T.B., Hau, B.C., He, F., Hoffman, F.M., Howe, R.W., Hubbell, S.P., Inman-Narahari, F.M., Jansen, P.A., Jiang, M., Johnson, D.J., Kanzaki, M., Kassim, A.R., Kenfack, D., Kibet, S., Kinnaird, M.F., Korte, L., Kral, K., Kumar, J., Larson, A.J., Li, Y., Li, X., Liu, S., Lum, S.K., Lutz, J.A., Ma, K., Maddalena, D.M., Makana, J.R., Malhi, Y., Marthews, T., Mat Serudin, R., McMahon, S.M., McShea, W.J., Memiaghe, H.R., Mi, X., Mizuno, T., Morecroft, M., Myers, J.A., Novotny, V., de Oliveira, A.A., Ong, P.S., Orwig, D.A., Ostertag, R., den Ouden, J., Parker, G.G., Phillips, R.P., Sack, L., Sainge, M.N., Sang, W., Sri-Ngernyuang, K., Sukumar, R., Sun, I.F., Sungpalee, W., Suresh, H.S., Tan, S., Thomas, S.C., Thomas, D.W., Thompson, J., Turner, B.L., Uriarte, M., Valencia, R., Vallejo, M.I., Vicentini, A., Vrska, T., Wang, X., Wang, X., Weiblen, G., Wolf, A., Xu, H., Yap, S., Zimmerman, J., 2015. CTFS-ForestGEO: a worldwide network monitoring forests in an era of global change. Glob. Change Biol. 21, 528–549. Baldeck, C.A., Harms, K.E., Yavitt, J.B., John, R., Turner, B.L., Valencia, R., Navarrete, H., Bunyavejchewin, S., Kiratiprayoon, S., Yaacob, A., Supardi, M.N., Davies, S.J., Hubbell, S.P., Chuyong, G.B., Kenfack, D., Thomas, D.W., Dalling, J.W., 2017. Habitat filtering across tree life stages in tropical forest communities. Proc. Biol. Sci. 280, 20130548. Bigler, C., Veblen, T.T., 2009. Increased early growth rates decrease longevities of conifers in subalpine forests. Oikos 118, 1130–1138. Bircher, N., Cailleret, M., Bugmann, H., 2015. The agony of choice: different empirical mortality models lead to sharply different future forest dynamics. Ecol. Appl. 25, 1303–1318. Bond-Lamberty, B., Fisk, J.P., Holm, J.A., Bailey, V., Bohrer, G., Gough, C.M., 2015. Moderate forest disturbance as a stringent test for gap and big-leaf models. Biogeosciences 12, 513–526. 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. Bugmann, H., 2001. A review of forest gap models. Clim. Change 51, 259–305. Bugmann, H., 1996. A simplified forest model to study composition along climate gradients. Ecology 77, 2055–2074. Bugmann, H., Seidl, R., Hartig, F., Bohn, F., Brůna, J., Cailleret, M., François, L., Heinke, J., Henrot, A.-J., Hickler, T., Hülsmann, L., Huth, A., Jacquemin, I., Kollas, C., LaschBorn, P., Lexer, M.J., Merganič, J., Merganičová, K., Mette, T., Miranda, B.R., NadalSala, D., Rammer, W., Rammig, A., Reineking, B., Roedig, E., Sabaté, S., Steinkamp, J., Suckow, F., Vacchiano, G., Wild, J., Xu, C., Reyer, C.P.O., 2019. Tree mortality submodels drive simulated long-term forest dynamics: assessing 15 models from the stand to global scale. Ecosphere 10, e02616. Burns, J.H., Strauss, S.Y., 2011. More closely related species are more ecologically similar in an experimental test. PNAS 108, 5302–5307. Cailleret, M., Dakos, V., Jansen, S., Robert, E.M.R., Aakala, T., Amoroso, M.M., Antos, J.A., Bigler, C., Bugmann, H., Caccianaga, M., Camarero, J.-J., Cherubini, P., Coyea, M.R., Čufar, K., Das, A.J., Davi, H., Gea-Izquierdo, G., Gillner, S., Haavik, L.J., Hartmann, H., Hereş, A.-M., Hultine, K.R., Janda, P., Kane, J.M., Kharuk, V.I., Kitzberger, T., Klein, T., Levanic, T., Linares, J.-C., Lombardi, F., Mäkinen, H., Mészáros, I., Metsaranta, J.M., Oberhuber, W., Papadopoulos, A., Petritan, A.M., Rohner, B., Sangüesa-Barreda, G., Smith, J.M., Stan, A.B., Stojanovic, D.B., Suarez, M.-L., Svoboda, M., Trotsiuk, V., Villalba, R., Westwood, A.R., Wyckoff, P.H., Martínez-Vilalta, J., 2019. Early-warning signals of individual tree mortality based on annual radial growth. Front. Plant Sci. 9, 1964. Canham, C.D., Papaik, M.J., Latty, E.F., 2001. Interspecific variation in susceptibility to windthrow as a function of tree size and storm severity for northern temperate tree

12

Forest Ecology and Management 450 (2019) 117496

Y. Zhu, et al.

50-year growth from Eucalyptus forests in eastern Australia. Ecol. Model. 222, 3261–3270. Paine, C.E., Norden, N., Chave, J., Forget, P.M., Fortunel, C., Dexter, K.G., Baraloto, C., 2012. Phylogenetic density dependence and environmental filtering predict seedling mortality in a tropical forest. Ecol. Lett. 15, 34–41. Piao, T., Comita, L.S., Jin, G., Kim, J.H., 2013. Density dependence across multiple life stages in a temperate old-growth forest of northeast China. Oecologia 172, 207–217. Pu, X., Jin, G., 2018. Conspecific and phylogenetic density-dependent survival differs across life stages in two temperate old-growth forests in Northeast China. For. Ecol. Manage. 424, 95–104. Rötheli, E., Heiri, C., Bigler, C., 2011. Effects of growth rates, tree morphology and site conditions on longevity of Norway spruce in the northern Swiss Alps. Eur. J. For. Res. 131, 1117–1125. Russo, S.E., Davies, S.J., King, D.A., Tan, S., 2005. Soil-related performance variation and distributions of tree species in a Bornean rain forest. J. Ecol. 93, 879–889. Sang, W., Chen, L., Ma, K., 1999. Research on succession model FOROAK of Mongolian oak-Korean pine (Quercus mongolica-Pinus koraiensis) forest. Acta Botanica Sinica 41, 658–668. Sang, W., Li, J., 1998. Dynamics modeling of Korean pine forest in southern Lesser Xiangan Mountains of China. Acta Ecologica Sinica 18, 38–47. Sato, H., Itoh, A., Kohyama, T., 2007. SEIB–DGVM: A new Dynamic Global Vegetation Model using a spatially explicit individual-based approach. Ecol. Model. 200, 279–307. Seidl, R., Rammer, W., Scheller, R.M., Spies, T.A., 2012. An individual-based process model to simulate landscape-scale forest ecosystem dynamics. Ecol. Model. 231, 87–100. Shao, G., 1991. Moisture-therm indices and optimum-growth modeling for the main species of Korean pine/deciduous mixed forests. Scientia Silvae Sinicae 27, 21–27. Shao, G., Schall, P., Weishampel, J.F., 1994. Dynamic simulations of mixed broadleavedPinus koraiensis forests in the Changbaishan biosphere reserve of China. For. Ecol. Manage. 70 (1–3), 169–181. Shao, G., 1996. Potential impacts of climate change on a mixed broadleaved-Korean pine forest stand: A gap model approach. Clim. Change 34, 263–268. Shi, B., Gao, W., Jin, G., 2015. Effects on rhizospheric and heterotrophic respiration of conversion from primary forest to secondary forest and plantations in northeast China. Eur. J. Soil Biol. 66, 11–18. Shugart, H.H., 1984. A theory of forest dynamics. Springer, New York, pp. 278. Shugart, H.H., Wang, B., Fischer, R., Ma, J., Fang, J., Yan, X., Huth, A., Armstrong, A.H., 2018. Gap models and their individual-based relatives in the assessment of the consequences of global change. Environ. Res. Lett. 13, 033001. Song, X., Bryan, B.A., Almeida, A.C., Paul, K.I., Zhao, G., Ren, Y., 2013. Time-dependent

sensitivity of a process-based ecological model. Ecol. Model. 265, 114–123. Tatarinov, F.A., Cienciala, E., 2006. Application of BIOME-BGC model to managed forests: 1. Sensitivity analysis. For. Ecol. Manage. 237 (1–3), 267–279. Uriarte, M., Canham, C.D., Thompson, J., Zimmerman, J.K., 2004. A neighborhood analysis of tree growth and survival in a hurricane-driven tropical forest. Ecol. Monogr. 74, 591–614. Uriarte, M., Clark, J.S., Zimmerman, J.K., Comita, L.S., Forero-Montana, J., Thompson, J., 2012. Multidimensional trade-offs in species responses to disturbance: implications for diversity in a subtropical forest. Ecology 93, 191–205. Wallace, L., Lucieer, A., Watson, C., Turner, D., 2012. Development of a UAV-LiDAR system with application to forest inventory. Remote Sens. 4, 1519–1543. Wang, X., Comita, L.S., Hao, Z., Davies, S.J., Ye, J., Lin, F., Yuan, Z., 2012. Local-scale drivers of tree survival in a temperate forest. PLoS ONE 7, e29469. Webb, C.O., Gilbert, G.S., Donoghue, M.J., 2006. Phylodiversity-dependent seedling mortality, size structure, and disease in a Bornean rain forest. Ecology 87, S123–S131. Wright, S.J., Kitajima, K., Kraft, N.J.B., Reich, P.B., Wright, I.J., Bunker, D.E., Condit, R., Dalling, J.W., Davies, S.J., Díaz, S., Engelbrecht, B.M.J., Harms, K.E., Hubbell, S.P., Marks, C.O., Ruiz-Jaen, M.C., Salvador, C.M., Zanne, A.E., 2010. Functional traits and the growth-mortality trade-off in tropical trees. Ecology 91, 3664–3674. Xu, L., Jin, G., 2012. Species composition and community structure of a typical mixed broadleaved-Korean pine (Pinus koraiensis) forest plot in Liangshui Nature Reserve, Northeast China. Biodivers. Sci. 20, 470–481. Yan, X., Shugart, H.H., 2005. FAREAST: a forest gap model to simulate dynamics and patterns of eastern Eurasian forests. J. Biogeogr. 32, 1641–1658. Zanne, A.E., Tank, D.C., Cornwell, W.K., Eastman, J.M., Smith, S.A., FitzJohn, R.G., McGlinn, D.J., O'Meara, B.C., Moles, A.T., Reich, P.B., Royer, D.L., Soltis, D.E., Stevens, P.F., Westoby, M., Wright, I.J., Aarssen, L., Bertin, R.I., Calaminus, A., Govaerts, R., Hemmings, F., Leishman, M.R., Oleksyn, J., Soltis, P.S., Swenson, N.G., Warman, L., Beaulieu, J.M., 2014. Three keys to the radiation of angiosperms into freezing environments. Nature 506, 89–92. Zhang, Z., 2008. Dendrology (North China Edition). China Forestry Publishing House, Beijing, China. Zhen, Z., Li, F., Liu, Z., Liu, C., Zhao, Y., Ma, Z., Zhang, L., 2013. Geographically local modeling of occurrence, count, and volume of downwood in Northeast China. Appl. Geogr. 37, 114–126. Zhu, Y., Cai, H., Jiang, F., Jin, G., 2017. Variation of biotic neighbourhood and topographic effects on tree survival in an old-growth temperate forest. J. Veg. Sci. 28, 1166–1177. Zhu, Y., Comita, L.S., Hubbell, S.P., Ma, K., 2015. Conspecific and phylogenetic densitydependent survival differs across life stages in a tropical forest. J. Ecol. 103, 957–966.

13