Comparing strategies for representing individual-tree secondary growth in mixed-species stands in the Acadian Forest region

Forest Ecology and Management 459 (2020) 117823

Contents lists available at ScienceDirect

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

Comparing strategies for representing individual-tree secondary growth in mixed-species stands in the Acadian Forest region

T

⁎

C. Kuehnea, , M.B. Russellb, A.R. Weiskittela, J.A. Kershaw Jr.c a

Center for Research on Sustainable Forests, University of Maine, Orono, ME 04469, USA Department of Forest Resources, University of Minnesota, St. Paul, MN 55108, USA c Faculty of Forestry and Environmental Management, University of New Brunswick, Fredericton, NB E3B 5A3, Canada b

A R T I C LE I N FO

A B S T R A C T

Keywords: Stem diameter growth Basal area increment Two-stage potential-modifier approach Realized diameter increment Mixed effect modeling Mixed species forests Forest growth and yield modeling

Tree diameter increment (ΔDBH) is a key component of a forest growth and yield model as predictions are passed to other submodels and tree-level estimates are scaled up to represent plot- and stand-level measures. A common problem faced in mixed-species stands is that ΔDBH needs to be characterized for numerous species, each with varying growth rates, shade tolerances, and competitive abilities. In addition, a variety of approaches have been used to model ΔDBH with unclear implications for general suitability for each species and overall prediction accuracy. This analysis used remeasurement data comprising 2,656,354 observations from 16,204 permanent sample plots across the Acadian Forest region of North America to develop and compare alternative approaches to estimating ΔDBH as well as stem basal area increment (ΔBA). Sixty-one species or genera including several with N < 10 were represented where observed mean growth rates ranged from 0 to 1.08 cm yr−1, depending on species. Analyzing several modeling approaches to project DBH of the 15 most abundant species using various evaluation statistics to quantify prediction performance, this study showed that i) modeling ΔDBH was generally superior compared to approaches that estimated ΔBA, ii) a two-stage modeling procedure predicting potential growth and a corresponding multiplicative modifier to derive ultimate increment was mostly inferior compared to strategies predicting realized ΔDBH or ΔBA in a unified model form, and iii) species-specific, realized increment models exhibited similar behavior and accuracy compared to models fitted with modeling species as random effect. These key findings became even more evident when projection lengths increased (≥ 30 years). Our study thus showed the efficiency and flexibility of diameter predictions by including tree species as a random effect to account for ΔDBH differences of trees in mixed-species stands, including infrequent species. However, curves for rare species derived with the mixed effects modeling approach still need to be evaluated for biological plausibility as unbalanced or biased data can lead to uncharacteristic and potentially illogical behavior. Overall, the study highlights the challenges of accurately predicting ΔDBH across a range of species and conditions, but offers a general framework for future analyses in mixed species forests.

1. Introduction Forests across Maine and the Canadian Maritime provinces are composed of tree species that occur in various assemblages and differ in terms of their growth rate, shade tolerance, and competitive ability. As a transition zone between the boreal forest to the north and the temperate deciduous forest to the south, this mixed-species forest type, termed the Acadian Forest, contains many species growing at the limits of their natural range (Braun, 1950; Rowe, 1972). Given the complexities of mixed-species stands, modeling individual-tree secondary growth often requires a novel approach. Furthermore, interpreting tree diameter measurements from permanent sample plots requires

⁎

statistical control for tree diameter growth, making growth dynamics studies complex (Bowman et al., 2013). As most stands of the Acadian Forest are regenerated by natural means, spatial information for trees within these stands is often lacking, which currently limits a model developer to constructing distance-independent growth models (Kuehne et al. 2019). In addition, the multi-cohort structure of many stands throughout the Acadian Forest dismisses the possibility of using stand or tree age as a potential covariate in explaining the growth of individual trees. The ability to accurately forecast tree diameter for each of the various tree species found within mixed-species, multi-cohort stands is essential to determining future forest structure and composition and for evaluating effective management regimes,

Corresponding author. E-mail address: [email protected] (C. Kuehne).

https://doi.org/10.1016/j.foreco.2019.117823 Received 7 November 2019; Received in revised form 12 December 2019; Accepted 13 December 2019 Available online 09 January 2020 0378-1127/ © 2019 Elsevier B.V. All rights reserved.

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

in mixed-species stands (Phillips et al., 2002; Picard et al., 2010; Vanclay, 1991; Zhao et al., 2006). However, several problems of grouping species remain. First, the species groupings that are derived are often specific to the dataset used in model parameterization (Vanclay, 1995; Zhao et al., 2006) and can result in non-logical groupings (e.g. Russell et al., 2014). Second, the question of the optimal number of groups to employ remains to be answered (Picard et al., 2010). Species groupings according to taxonomy may seem fitting in a biological sense, but may be inappropriate in terms of quantifying growth (Zhao et al., 2006). Lastly, while the method of grouping may successfully categorize the most common species into unique groups, determining how rare and/or infrequent species are accounted for remains an important and unanswered question (Picard et al., 2010). These infrequent species may be placed into groups subjectively (Phillips et al., 2002), while some argue that the characteristics of these rare species should not go into defining groupings (Picard et al., 2010). Given these important uncertainties, additional methods for developing species-specific equations, particularly for rare species, are needed. Statistical techniques have advanced rapidly in recent years for fitting and evaluating individual-tree increment models. For example, nonlinear mixed-effects models (NLME) have been shown to outperform annual ΔDBH predictions fitted with generalized nonlinear leastsquares (Weiskittel et al., 2007). In naturally-regenerated, mixed-species stands, growth equations have primarily been developed on a species-by-species basis (e.g. Weiskittel et al., 2016). A quantitative strategy that advances and extends the NLME approach and thus eliminates the need to obtain a set of parameters for each species (or species group) would be to consider each species as a random element of the stand mixture. The use of species as random effect has been used for a variety of growth and yield models including individual tree taper/volume (Weiskittel et al., 2015), biomass (Colmanetti et al., 2018), height (Lam et al., 2016), and height increment (Russell et al., 2014). The extension of this approach to ΔDBH is logical. In addition, future efforts thus could consider model formulations for ΔDBH or ΔBA which generate a universal growth equation that can be applied to trees of all species in mixed-species stand structures. The overall objective of this analysis was to evaluate several methods for incorporating tree species in individual-tree secondary growth equations for use in the mixed-species Acadian Forest of northeastern North America. Specific objectives were to (1) develop and compare ΔDBH equations with models predicting ΔBA, (2) develop and compare increment equations derived with the POTMOD approach with models based on the REAL approach, (3) develop and compare species-specific increment equations with models incorporating tree species as a random effect, and (4) evaluate and validate all model forms to assess their performance.

particularly with respect to carbon sequestration potential. Projecting individual tree diameter at breast height (DBH) can be accomplished by a variety of means with the most common approaches predicting diameter (ΔDBH) or stem basal area increment (ΔBA). The optimal dependent variable has often been debated (e.g. Cole and Stage, 1972; MacFarlane and Kobe, 2006; Russell et al., 2011). Given that most studies measure stem circumference and convert it to a diameter under the assumption of a circle, the use of ΔBA as a dependent variable further extends the assumption of a circle, resulting in more mensuration error in the response variable and more unexplainable variation about the equation (Kmenta, 1997; Weiskittel et al., 2011b). For example, Hann and Larsen (1991) found that transforming BA to DBH resulted in unreasonable predictions for trees with small diameters. Similarly, Russell et al. (2011) indicated that prediction of ΔDBH was superior to ΔBA, particularly as projection lengths increased. This is likely because basal area assumes that trees are circular and noneccentric with the importance of this assumption potentially increasing with time (Russell et al., 2011). However, these findings were based on a limited dataset on a few species so it would be important to evaluate it on a larger dataset with extensive long-term measurements and a greater number of available species. Two primary strategies for model formulations have been employed for estimating secondary growth of individual trees (Weiskittel et al., 2011b). The first strategy predicts increment in two stages: first by estimating an annual maximum or potential increment for an individual tree and second by adjusting those estimates using multiplicative modifiers. This potential-modifier approach (POTMOD) is advantageous since it constrains the ultimate prediction of ΔDBH or ΔBA not to exceed maximum increment and has been widely used in a variety of studies (Larocque et al., 2011, Huang et al. 2013, Strimbu et al., 2017). Estimating potential ΔDBH in eastern North America has traditionally been accomplished by parameterizing models using i) a subset of the available data, i.e. dominant/co-dominant individuals only or a certain percentile of the largest or fastest growing trees of a stand (Wensel et al.,1987; Teck and Hilt, 1991), or ii) open-grown trees (Moore, 1989). Recently, nonlinear quantile regression methodologies (Koenker and Hallock, 2001) have become available in statistical software packages that enable one to quantitatively estimate potential ΔDBH. This method negates the need to subjectively select a subset of trees for model parameterization as demonstrated by Pretzsch and Biber (2010). A second strategy predicts secondary tree growth in a unified approach where increment is estimated using tree size and vigor, competition, and site quality in a single equation. This realized diameter increment approach (REAL) has seen limited use in the forests which comprise northeastern North America; however, the methods have performed well in other regions (Wykoff, 1990; Monserud and Sterba, 1996; Cao, 2000; Hann et al., 2003; Weiskittel et al., 2007). Regardless of whether the POTMOD or REAL approach is used, prediction of future growth has tremendous implications in terms of estimating future forest structure and composition as DBH is generally the greatest source of prediction uncertainty in a growth model (Wilson et al., 2019), particularly in the mixed species forest of the northeastern US (Russell et al., 2013). In estimating ΔDBH or ΔBA in two stages, Lessard et al. (2001) suggests that quantifying model prediction uncertainty cannot be accomplished due to the difficulties in estimating the covariance between the maximum and modifier equations. Assessing the performance of using both approaches would need to be determined for use in stands with various species assemblages in even- and uneven-aged stand structures. Increment models are regularly developed for individual species or species groups. As an example, the Forest Vegetation SimulatorNortheast Variant (FVS-NE) is a widely-used distance-independent growth model in the region that includes sets of model coefficients predicting ΔDBH for 26 different species or species groups (Dixon et al., 2018), but its general performance can be highly variable spatially for a given species in this region (Russell et al., 2013). Grouping species with similar attributes has historically been a method for quantifying growth

2. Methods 2.1. Study area The Acadian Forest region resides in the transition zone between the softwood-dominant boreal forests to the north and the broadleaf forests to the south (Braun, 1950; Rowe, 1972). Three Canadian Maritime Provinces (New Brunswick, Nova Scotia, and Prince Edward Island) are found in the region, along with southern portions of Québec, and much of the US state of Maine. Across the region, climate estimates indicate average annual precipitation is 113 cm with a range of 87 to 175 cm, while mean growing degree days (sum of temperature > 5 °C) is 1,625 with a range of 726 to 2,292 days (Rehfeldt, 2006). Glacial till is the principal soil parent material with soil types ranging from well-drained loams and sandy loams on glacial till ridges to poorly and very poorly drained loams on flat areas between low-profile ridges. The majority of the Acadian Forest is dominated by naturally-regenerated, mixed-species stands that display either even- or unevenaged stand structures, while some portions of New Brunswick contain 2

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

intensively-managed, single-species plantations. Over 60 species dominate the region, and common forest types in the region include black spruce (Picea mariana (Miller) B.S.P) - northern white-cedar swamps (Thuja occidentalis L.), red spruce (Picea rubens Sarg.) - balsam fir (Abies balsamea L.) flats, red spruce - balsam fir - red maple (Acer rubrum L.) yellow birch (Betula alleghaniensis Britton), sugar maple (Acer saccharum Marsh.) - American beech (Fagus grandifolia Ehrh.) – yellow birch, aspen (Populus spp.) - paper birch (Betula papyrifera Marsh.), eastern white pine (Pinus strobus L.) - mixed softwood forest, and northern red oak (Quercus rubra L.) – eastern white pine (Eyre, 1980). Other species present include white spruce (Picea glauca (Moench) Voss.), Norway spruce (Picea abies (L.) Karst.), red pine (Pinus resinosa Ait.), Scots pine (Pinus sylvestris L.), eastern hemlock (Tsuga canadensis (L.) Carr.), paper birch (Betula papyrifera Marsh.), quaking aspen (Populus tremuloides Michx.), bigtooth aspen (Populus grandidentata Michx), black cherry (Prunus serotina Ehrh.), and white ash (Fraxinus americana L.).

Experimental Forest (PEF), located in the towns of Bradley and Eddington, Maine (44° 52′ N, 68° 38′ W). The PEF is characterized by a mixture of northern softwood and hardwood species that dominate its forest cover. As part of a replicated long-term silvicultural study, ten contrasting silvicultural systems were installed by USFS between 1952 and 1957. Beginning in 1974, individual trees were numbered in a network of PSPs established along transects traversing each stand. PSPs consisted of a nested design with 0.081- and 0.020-ha, circular plots sharing the same plot center. PSPs were initially measured then remeasured at 5-year intervals, provided no harvesting was scheduled within the stand. Beginning in 2000, a 0.008-ha plot was incorporated into the nested plot design and remeasurements were recorded on a 10year interval. A complete description of the measurements collected and the data collection procedures can be found in Kenefic et al. (2015). A total of 226 PEF PSPs and measurements through 2011 were included in this analysis.

2.2. Data

2.2.3. Commercial Thinning Research Network (CTRN) Individual tree growth records were obtained from 16 research locations across Maine as part of the Commercial Thinning Research Network (CTRN) established by the Cooperative Forestry Research Unit (Wagner and Seymour, 2006). Three controlled experiments, established in mostly red spruce-balsam fir dominated stands comprise the CTRN (Bataineh et al., 2013; Kuehne et al., 2018c). A total of 140 0.081 ha measurement plots installed across the research locations were used in this study. A thinning prescription was applied in a research area surrounding each measurement plot that differed in terms of type, intensity, and timing of thinning (Kuehne et al., 2016), which also resulted in varied residual tree spatial patterns (Kuehne et al., 2018c). Tree DBH measurements were initially made in 2001 and remeasured annually or every two years through 2015.

Individual tree DBH measurements for this study were obtained from an extensive database of fixed-area permanent sample plots (PSPs) compiled from a variety of data sources throughout the region (Weiskittel et al., 2010; Li et al., 2011). A summary of plot-level metrics is provided in Table 1. A more specific description of each dataset is described below. 2.2.1. Maine Forest Inventory and Analysis (FIA) Tree measurements were obtained from a network of PSPs at various locations across Maine as part of the US Forest Service (USFS) Forest Inventory and Analysis (FIA) Program (Bechtold and Patterson, 2005). Approximately one PSP location was established for each 2,400 ha in a hexagonal grid. Four 7.32 m radius, fixed area subplots were established at each location and trees ≥ 12.7 cm DBH were measured. Trees ≥ 2.5 cm DBH < 12.7 cm were measured on a smaller 2.07-m radius microplot nested within the larger one. Measurement of plots began in 1999 when the FIA annual inventory design was implemented across the US (McRoberts et al., 2005). Given the 5year measurement interval, data were recorded up to three times. Data were obtained from the online FIA database at https://apps.fs.usda. gov/fia/datamart/CSV/datamart_csv.html and downloaded in October 2018 with a total of 2,941 and 10,557 PSPs and subplots, respectively, used for this study. Stand metrics were calculated at the subplot-level.

2.2.4. Maine Ecological Reserves Network The Maine Natural Areas Program’s Ecological Reserves Network currently contains 50 long-term permanent reserves across the state of Maine (Kuehne et al., 2018a). Starting in 2002, mostly five to six PSPs were established along randomly established transects with PSPs spaced approximately 241 m apart along these transects. The nested plot design consisted of a 17.95 m radius plot for trees ≥ 51 cm DBH and a 7.32 m radius plot for trees ≥ 12.7 cm DBH sharing the same plot center (Kuehne et al., 2018b). As of 2017, a subset of PSPs in 20 reserves were remeasured 10 years after initial inventory and 323 of these PSPs were used here.

2.2.2. Penobscot Experimental Forest (PEF) Long-term individual tree growth records were obtained from various research installations occurring at the USFS Penobscot

2.2.5. New Brunswick PSP Individual DBH measurements were obtained from 1,398 PSPs that were 0.04 ha in size established across New Brunswick in a mix of multi-aged, mixed-species stand types (Province of New Brunswick, 2005; McGarrigle et al., 2011). Data collection was initiated in 1985 with most measurement intervals lasting three to five years.

Table 1 Plot-level (N = 16,204) summaries including standard deviation (SD) for diameter increment data gathered from mixed-species stands across the Acadian Forest region of North America. Attribute

Mean

SD

Min

Max

Plot size (m2) Interval length (years) Longitude (degrees) Latitude (degrees) Elevation (m) Climate site index (m) Stem density (trees ha−1) Stand density index (trees ha−1) Relative density Basal area (m2 ha−1) Percent basal area in hardwoods (%) Quadratic mean diameter (cm) Species diversity (# plot−1) Shannon diversity index for species

253.62 10.67 −68.78 45.80 255.94 13.90 2,450.71 399.39 0.44 22.98 41.18

135.86 7.53 2.16 1.15 188.43 2.42 3,208.58 261.94 0.27 13.58 36.74

168.11 1 −73.25 43.11 0 4.77 9.88 0 0 0.02 0

810.11 40 −59.81 49.22 1,095 31.03 35,685.67 1,858.76 1.79 118.56 100

14.49 3.79 0.89

6.55 1.71 0.46

2.04 1 0

76.5 12 2.11

2.2.6. Québec PSP Tree DBH measurements were obtained from a variety of PSPs located across Québec from various datasets (Weiskittel et al., 2010; Li et al., 2011). Datasets included: the first plot network that was established between 1970 and 1977 (PQ-BAS1); the second plot network that was established in 1989 (PQ-BAS2); Fédération des producteurs de bois du Québec (PQ-FEDE); Parks Canada (PQ-PACA); Service de la Comptablité Forestière (PQ-SCOF); Service de la protection des insectes et des maladies (PQ-SPIM); and University of Laval (UNLA). Plot size across the 2,337 PSPs used here was 0.04 ha and mean measurement intervals ranged from approximately five to 11 years. 2.2.7. Nova Scotia PSP Individual DBH measurements were obtained from 1,227 PSPs that were 0.04 ha in size located across Nova Scotia (Weiskittel et al., 2010; 3

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

ΔDBH = ΔDBHPOT × ΔDBHMOD

Li et al., 2011). Tree DBH measurements began in 1965, and measurement intervals averaged five years.

where ΔDBHPOT and ΔDBH are species-specific potential and ultimate predicted annual diameter increment (cm yr−1), respectively, and ΔDBHMOD is the multiplicative modifier. The ΔDBHPOT equation is often characterized by using measures of initial tree size and vigor as well as of site quality, while ΔDBHMOD is often predicted using measures of one- (e.g., basal area in larger trees, BAL) and two-sided competition (e.g., stand-level basal area or crown competition factor, see e.g. Kuehne et al., 2019). An accurate estimate of site productivity proved problematic because height-age information was extremely limited across the study area and the mixed-species uneven-aged stand structures prevented the use of traditional site index (SI) equations. Because SI has been a significant explanatory variable of ΔDBHPOT in several studies (e.g., Dixon et al., 2018; Pretzsch and Biber, 2010), a climate-derived SI (CSI, m) was estimated. CSI is a base–age 50 years estimate of site index based on downscaled climate variables (http://forest.moscowfsl.wsu.edu/ climate/) and a random forest nonparametric regression model built on field–observed values. (Weiskittel et al., 2011a; Weiskittel et al., 2011c). A model form modified after that of Pretzsch and Biber (2010) resulting in a sigmoidal model curve was carried out in the first modeling stage to estimate ΔDBHPOT:

2.3. Data preparation All available measurements from the various data sets were merged into a common format and converted to metric units. Individual tree total height and height to crown base values missing were imputed by using species-specific equations published in Russell et al. (2014). Instead of modeling ΔDBH (cm yr−1) and ΔBA (m2, yr−1) based on successive measurements only, similar to Chen et al. (2017) measurements of all possible intervals were used to increase number of observations. More precisely, growth data were not just derived from consecutive inventories (e.g. year1-year2, year2-year3, year3-year4, and so on) but all potential combinations, i.e. including year1-year3, year1-year4, year2year4, and so on. A total of 621 observations or 0.0002% of the available data where the diameter of an individual tree in subsequent inventories was less than that from a previous census were removed from the analysis. In addition, measurement periods including a cleaning or harvest operation were also not included in the analysis. A total of 2,656,354 observations were available for this study. Overall, 54 species recorded to the species level were present, including 39 hardwoods and 15 softwoods (Table 2). In addition, 0.1% or 2,727 of all observations were recorded to the genus level, including Alnus spp., Amelanchier spp., Cornus spp., Crataegus spp., Malus spp., Salix spp., and Sorbus spp. (all hardwoods). The top 15 most frequent species (8 softwoods and 7 hardwoods) comprised 96.8% or 2,572,361 of all observations and were in decreasing order of abundance: balsam fir, red spruce, red maple, black spruce, paper birch, white spruce, sugar maple, northern white-cedar, yellow birch, eastern hemlock, eastern white pine, quaking aspen, American beech, tamarack (Larix laricina (Du Roi) K. Koch), and northern red oak.

ΔDBHPOT = α 01 exp(−a11 DBH ) DBH a21 exp(α31 CR + 1) exp(α 41 ln(CSI )) (2) where DBH is initial tree diameter (cm), CR is crown ratio (ratio of crown length and total tree height, %), and αijs are species-specific parameters to be estimated. Nonlinear quantile regression was used to estimate ΔDBHPOT with the nlrq function found in the quantreg package (Koenker, 2013). Specifying the 99th percentile, predictions were assumed to represent the species maximum diameter growth potential (Pretzsch and Biber, 2010). The second stage of SS-POTMOD was to estimate a multiplicative diameter increment modifier (ΔDBHMOD) for each species. This was characterized by one- and two-sided competition:

2.4. Model development For modeling ΔDBH and ΔBA, single-species equations (SS) were first developed using the POTMOD (SS-POTMOD) and REAL (SS-REAL) approaches based on observations of the 15 most abundant species. Then, tree observations from all species (AS) were used to develop a single equation that incorporated species as a random effect (AS-REAL). With the exception of tree size (DBH vs BA), all increment equations were fitted with the same set of predefined explanatory variables to maximize comparability across the evaluated modeling approaches. To overcome problems of varying measurement intervals (1–40 years) observed in the data and in order to provide a finer resolution of treeand stand-level dynamics, parameters were annualized using an iterative technique of Weiskittel et al. (2007). Based on Cao (2000) the right side of the equation was a function which summed the estimated annual ΔDBH or ΔBA estimates, respectively, over the number of growing seasons during the observed growth period using the updated parameter estimates from the optimization algorithms. For each growing season during the growth period, DBH or BA was subsequently updated, while all explanatory variables were linearly interpolated between their beginning values and ending values except CSI which was assumed to be constant over time. All models were fitted using the programming software R (R Development Core Team, 2018). Each modeling approach is further described below.

ΔDBHMOD = exp[β01 + β11 BAL SW + β21 BALHW + β31 RD]

(1a)

ΔDBHMOD = f (one - , two - sided competition)

(1b)

(3)

where BALSW, and BALHW are the basal area in larger trees in softwood and hardwood trees, respectively, RD is relative density, and ßij are species specific parameters estimated with the gnls function found in the nlme package (Pinheiro et al., 2012). Separating BAL into softwood and hardwood components allowed for species-type differences that has shown to work well for stands with mixed species assemblages (Nunifu, 2009). For this as well as all following model forms that include metrics of site occupancy, plot-level measures of competition were tested that are applicable for these stand types with various species assemblages. Stand-level basal area (m2 ha−1), crown competition factor (%), and RD were examined with RD performing best among the three competition measures in preliminary analyses. RD is defined as the ratio of stand density index (SDI) and maximum stand density index (SDIMAX). SDI was computed based on the summation method (Weiskittel et al. 2011b) while SDIMAX was quantified using the equation provided in Weiskittel and Kuehne (In Press) for mixed-species stands of the region. As for all subsequent modeling approaches, explanatory variables including RD were removed from a specific model if they were insignificant or showed illogical behavior, i.e. a wrong sign on a parameter estimate was observed. Parameters of Eq. (3) of stage 2 were derived by fitting Eq. (1c) after incorporation of ΔDBHPOT estimates from stage 1. Consequently, final ΔDBH predictions for the SSPOTMOD approach were computed as ΔDBHPOT × ΔDBHMOD. Similarly, individual tree annual ΔBA was modeled using the SSPOTMOD approach after DBH measurements were converted to basal area (BA, m2) to compute stem basal area increment (ΔBA, m2 yr−1). However, due to model stability and convergence issues model form for

2.4.1. Single-Species, Potential-Modifier (SS-POTMOD) The process of modeling annual ΔDBH of trees using the SSPOTMOD approach included two separate modeling stages and took the following general form:

ΔDBHPOT = f (tree size, tree vigor, site productivity)

(1c)

4

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

Table 2 Species acronyms and names with total number of measurements (N) and statistics including standard deviation (SD) for initial diameter at breast (DBH, cm) and mean annual DBH increment (ΔDBH, cm yr−1) gathered from mixed-species stands across the Acadian Forest region of North America. Acronym

AB AE AH AI AL AP AW BA BC BE BF BH BN BO BP BR BS BT BW CC DW EC EH GA GB HH HT JP MA MM NM NS PB PP PR PY QA RM RN RO RS SB SC SE SH SM SO ST SV SW SY TA TM WA WC WI WO WP WS YB YP Overall

Scientific name

Fagus grandifolia Ulmus americana Carpinus caroliniana Ailanthus altissima Alnus spp. Malus spp. Chamaecyparis thyoides Fraxinus nigra Prunus serotina Acer negundo Abies balsamea Carya cordiformis Juglans cinerea Quercus velutina Populus balsamifera Quercus macrocarpa Picea mariana Populus grandidentata Tilia americana Prunus virginiana Cornus spp. Populus deltoides Tsuga canadensis Fraxinus pennsylvanica Betula populifolia Ostrya virginiana Crataegus spp. Pinus banksiana Sorbus spp. Acer spicatum Acer platanoides Picea abies Betula papyrifera Pinus rigida Prunus pensylvanica Populus heterophylla Populus tremuloides Acer rubrum Pinus resinosa Quercus rubra Picea rubens Betula lenta Pinus sylvestris Amelanchier spp. Carya ovata Acer saccharum Quercus coccinea Acer pensylvanicum Acer saccharinum Quercus bicolor Platanus occidentalis Larix laricina Pinus pungens Fraxinus americana Thuja occidentalis Salix spp. Quercus alba Pinus strobus Picea glauca Betula alleghaniensis Liriodendron tulipifera -

N

43,726 689 89 1 304 203 10 3842 1724 7 853,205 7 18 251 2675 1 224,734 10,675 340 90 3 15 62,619 298 14,356 3629 37 11,651 1502 1964 7 68 142,947 185 3994 3 51,628 295,414 3804 14,586 482,112 138 30 303 13 85,809 3 9116 211 4 1 18,969 2 11,004 73,828 375 349 54,476 97,347 70,961 2 2,656,354

ΔDBH

DBH Mean

SD

Min

Max

Mean

SD

Min

Max

14.55 14.98 5.97 3.56 6.31 15.37 21.08 12.17 14.55 14.59 12.38 17.23 22.47 23.74 16.83 12.7 12.36 17.73 18.84 4.94 2.28 8.97 18.46 15.5 8.18 11.45 5.68 17.48 11.18 5.35 6.42 16.07 13.19 24.43 9.6 11.6 17.09 14.64 22.41 18.7 15.61 19.62 12.14 9.3 14.15 17.13 14.99 6.92 17.92 31.05 13.97 15.75 5.33 16.45 17.62 11.94 18.09 22.1 17.12 18.63 16.76 14.46

8.16 7.2 5.16 – 1.2 7.01 5.75 7.06 6.77 4.42 5.91 4.13 8.73 8.19 10.39 – 5.72 8.68 7.61 6.33 0.44 9.23 11.14 9.21 4.53 6.02 2.64 5.28 5.17 2.06 4.98 10.1 6.65 7.73 4.49 7.19 8.59 7.53 10.54 8.57 7.16 7.5 2.23 4.61 5.14 9.67 0 4.07 9.81 0.48 – 6.31 0 8.52 8.53 4.94 7.27 12.64 7.52 10.71 0 7.67

1.27 2.54 2.54 – 5.1 3.3 12.7 1 2.54 10.41 1 13.5 9.4 5.08 1.27 – 1 1.78 1.27 2.54 2.03 2.79 1.1 2.54 1.27 1.27 2.54 3.2 1.1 1.6 3.81 4.4 1 12.95 1.27 3.3 1.27 1 2.54 1 1 5.08 9.1 1.27 5.33 1 14.99 1 9.1 30.73 – 1 5.33 1.52 1 3.05 2.54 1.27 1 1 16.76 1

63.6 54.86 48.51 – 11.3 42.42 30.48 52 42.42 22.1 48.2 24.8 44.4 45.97 63.75 – 70 69.6 48.3 44.2 2.79 32 88.6 45.21 32.51 35.2 13.97 38 49.1 20.32 17.53 46 64.52 56.39 34.54 15.75 64 77.98 70.87 84.07 65.2 42.16 16 22.35 18.8 85.85 14.99 25.8 63.25 31.75 – 65 5.33 93 98.7 74 39.62 105 68.9 82 16.76 105

0.19 0.37 0.09 0.20 0.07 0.17 0.29 0.15 0.2 1.08 0.26 0.3 0.66 0.42 0.26 0.3 0.13 0.33 0.24 0.12 0.07 0.56 0.28 0.23 0.15 0.11 0.16 0.14 0.18 0.1 0.38 0.58 0.14 0.25 0.17 0.22 0.29 0.18 0.32 0.25 0.2 0.19 0.39 0.1 0.1 0.19 0.53 0.17 0.42 0.28 0.51 0.21 0 0.24 0.17 0.15 0.19 0.45 0.25 0.22 0.27 0.22

0.15 0.26 0.08 – 0.06 0.16 0.1 0.12 0.19 0.12 0.2 0.1 0.39 0.22 0.2 – 0.1 0.21 0.19 0.12 0.06 0.72 0.19 0.18 0.14 0.1 0.16 0.11 0.16 0.1 0.33 0.29 0.12 0.21 0.15 0.16 0.2 0.14 0.23 0.19 0.16 0.14 0.19 0.08 0.06 0.15 0.07 0.13 0.25 0.15 – 0.16 0 0.19 0.12 0.13 0.14 0.32 0.19 0.17 0.09 0.18

0 0 0 – 0 0 0.16 0 0 0.86 0 0.17 0.14 0 0 – 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.2 0.03 0 0 0 0.05 0 0 0 0 0 0 0.06 0 0 0 0.46 0 0 0.2 – 0 0 0 0 0 0 0 0 0 0.20 0

1.47 1.32 0.41 – 0.3 0.76 0.52 1.83 1.17 1.22 2.72 0.42 1.32 1.12 1.73 – 1.83 1.83 0.89 0.56 0.13 2.47 1.88 1.12 2.54 0.97 0.66 1.03 1.33 0.82 1.12 1.21 2.54 0.98 1.06 0.36 2.79 2.48 1.57 1.42 2.2 0.71 0.93 0.41 0.2 2.48 0.58 1.29 1.1 0.51 – 1.47 0 1.32 1.78 0.76 0.61 2.38 2.04 2.34 0.33 2.79

estimating potential ΔBA was changed and ultimate ΔBA was thus derived as follows:

ΔBAPOT = exp(α 02 + α12 ln(BA) + α22 BA + α32 CR + α42ln(CSI))

(4a)

ΔBAMOD = exp(β02 + β12 BAL SW + β22 BALHW + β32 RD)

(4b)

ΔBA = ΔBAPOT × ΔBAMOD

(4c)

where ΔBAPOT and ΔBA are species-specific potential and ultimate predicted annual basal area increment (m2 yr−1), respectively, ΔBAMOD is the multiplicative modifier, and BA is stem basal area (m2) and all other variables and parameters are defined above.

5

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

2.4.2. Single-Species, Realized (SS-REAL) For the SS-REAL approach, a unified equation was fitted to measurements of each species which combined tree size and vigor, competition, and site productivity variables in a single equation. The equation was modified after Weiskittel et al. (2016):

MB% =

MAB =

ΔDBH

(yi − yi )̂ ∗ 100⎞⎟ yi ⎝ ⎠

n

∑i=1 ⎛

⎜

n

∑i=1 (|yi − yi|)̂

MAB% =

(5a)

RMSE =

= exp(γ02 + γ12ln(BA) + γ22 BA + γ32 (5b)

where all variables are as previously defined and parameters γ ij were estimated with the gnls function found in the nlme package (Pinheiro et al., 2012). 2.4.3. All Species, Realized, species as a random effect (AS-REAL) Incorporating tree species as a random effect within a ΔDBH or ΔBA equation would seemingly account for variation in growth associated with each species. Here, fixed-effect parameters could account for population-level factors influencing increment, while random effects could account for subject-specific (i.e., species) influences that lead to differences in growth. Careful consideration and evaluation would be needed to determine which parameters to consider as random in a mixed-modeling framework. Incorporating species as a random effect would be advantageous in that it would reduce i) the amount of coefficients to be used for numerous species growing in mixed-species stands as well as ii) the complexity of fitting species-specific equations. Modeling species as random effect could also be applied to predict growth for rare and/or infrequent species where having enough observations to generate a biologically-consistent and plausible growth curve is often a limitation. Here, we optimized the parameters to be random based on preliminary analysis by testing various combinations of random effects with the best approach selected based on Akaike’s information criterion (AIC). For the AS-REAL approach, all tree observations regardless of species were used as the data source:

EI =

̂ (|yi − yi|) ∗ 100⎞⎟ y i ⎝ ⎠

n

⎜

(6a)

(6b) where parameters δij and dij,SP are fixed and species-specific random effects, respectively, estimated with the nlme function found in the nlme package (Pinheiro et al., 2012) and all other variables are as previously defined. 2.5. Model evaluation Uncertainty in model predictions includes both systematic and random variation (Kangas 1999). Based on the 15 most abundant species, we thus calculated various measures of model prediction accuracy to evaluate both types, namely mean bias (MB), relative MB, mean absolute bias (MAB), relative MAB, and root mean square error (RMSE): n

n

n (7d)

n

(7e)

l

( (∑

nk i=1

)

|yik − y^ik| nk DBHk

)

(7f)

(m − 1)(Si − Smin ) Smax − Smin

(8)

where rRi is the relative rank of modeling approach i (i = 1,2,…,m), Si is the evaluation metric derived for approach i, and Smin and Smax are minimum and maximum values of the evaluation metric across all approaches, respectively. Absolute values had to be used in order to derive meaningful rR for MB and MB%. Because the magnitude and not only the order of S is taken into consideration, rR provides more information than traditional ordinal ranks (Poudel and Cao, 2013), especially when averaged across a number of statistics for a final appraisal as done in this study. We further evaluated model fit of the various approaches by calculating the above described evaluation metrics for i) trees ≥ 50 cm DBH, and ii) observations with a measurement period length ≥ 30 years as well as by comparing prediction accuracy for iii) hard- and softwoods, and iv) species of varying shade tolerance. Shade tolerance was derived from the shade tolerance scale by Niinemets and Valladares (2006) with species-specific values ranging from 0.98 to 2.32, 2.33 to 3.67, and 3.68 to 5.01 indicating low (shade-intolerant), moderate (intermediate or moderately shade-tolerant), and high shade tolerance (shade-tolerant), respectively.

CR+δ42 BAL SW + (δ52 + d52,SP)BALHW + δ62 RD + δ72ln(CSI))

∑i=1 (yi − yi)̂

n

∑i=1 (yi − yi)̂ 2

∑k=1

rRi = 1 +

ΔBA = exp(δ02 + d 02,SP + (δ12+d12,SP)ln(BA) + δ22 BA + (δ32 + d32,SP)

MB =

(7c)

where yik and yik̂ are the ith observed and predicted DBH in 5 cmdiameter class k, nk is the total number of observations in diameter class k, l is the total number of 5 cm-diameter classes, and DBHk is the midpoint DBH of diameter class k. MB and MB% closer to zero signal greater prediction accuracy while the lower MAB, MAB%, RMSE, and EI the better. Finally, we employed Poudel and Cao’s (2013) relative rank (rR) to rate each modeling approach based on the studied evaluation metrics. rR is defined as:

ΔDBH = exp(δ01 + d 01,SP + (δ11+d11,SP)ln(DBH) + δ21 DBH + (δ31 + d31,SP) CR+δ41 BAL SW + (δ51 + d51,SP)BALHW + δ61 RD + δ71ln(CSI))

n

where yi is the observed DBH, yi ̂ is the predicted DBH, and n is the number of observations. Predicted DBH were calculated by applying the derived annualized ΔDBH or ΔBA equations of each modeling approach in a stepwise, i.e. annual manner while linearly interpolating explanatory variables in the same stepwise fashion. Given its application in other subsequent individual-tree submodels of forest growth and yield simulators (e.g. total height and mortality) as well as its importance in the quantification of stand-level summary statistics (e.g. basal area, total stem volume), DBH and not ΔDBH and ΔBA, respectively, was used as the evaluation metric. Based on the notion that bias in larger, i.e. more valuable trees is of greater significance (Pogoda et al., 2019) we also derived the combined error index (EI) after Reynolds et al. (1988) which is the sum of weighted MAB by diameter class and calculated as:

ΔBA CR + γ42 BAL SW + γ52 BALHW + γ62 RD + γ72ln(CSI))

(7b)

∑i=1 ⎛

= exp(γ01 + γ11ln(DBH) + γ21 DBH + γ31 CR + γ41 BAL SW + γ51 BALHW + γ61 RD + γ71ln(CSI))

n

3. Results For the data used in model development, DBH of all observations across all species/genera averaged 14.46 ± 7.67 cm (mean ± SD, Table 2). Minimum DBH was 2.5 cm or lower for most species/genera, while maximum DBH varied substantially with species/genera averaging 44.59 ± 25.82 cm. Annual ΔDBH for all observations across all

(7a) 6

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

species/genera averaged 0.22 ± 0.18 cm yr−1. For species or genera with > 10 observations, mean ΔDBH was a low as 0.07 ± 0.06 cm yr−1 for alder (Alnus spp.) and as high as 0.66 ± 0.39 cm yr−1 for butternut (Juglans cinera L.) (Table 2). Among the 15 most abundant species, eastern white pine exhibited the greatest (0.45 ± 0.32 cm yr−1) and black spruce the lowest mean annual ΔDBH (0.13 ± 0.10 cm yr−1). ΔDBHPOT curves developed for the 15 most abundant species in the SS-POTMOD approach showed a sigmoidal behavior and indicated a strong relationship with initial DBH, CR, and CSI, while most respective species-specific ΔDBHMOD equations exhibited strong relationships with BALSW, BALHW, and RD. In contrast, a second initial BA term, i.e. in addition to ln(BA), could not be incorporated into any ΔBAPOT equation derived with the SS-POTMOD approach despite a change in model formulation as result of convergence issues (see equations (2) and (4a)). Consequently, ΔBAPOT curves showed a unimodal concave instead of a sigmoidal behavior. In addition, CSI and RD were often not influential in the respective ΔBAPOT or ΔBAMOD functions, respectively. In general, there was great similarity with regard to influential explanatory variables in ΔDBH and ΔBA equations derived in the SS-REAL approach. Irrespective of the examined response, RD and CSI were insignificant or biologically not plausible for some of the 15 most abundant species and thus omitted from the corresponding SS-REAL models. Except for the tamarack and red spruce ΔBA equations, a second significant DBH or BA term, respectively, was part of all SS-REAL models (see equations (5a) and (5b)). The fixed effect parameter estimate for RD exhibited a positive sign in the ΔDBH and ΔBA equation derived in the NLME approach with species put as random. RD was thus removed from the final AS-REAL ΔDBH and ΔBA models (Table S1). All other fixed parameter estimates (see equations (6a) and (6b)) were significant and performed as expected and therefore were retained in the respective ΔDBH and ΔBA equation. Putting CSI as random resulted in overall negative parameter signs (combined fixed and random effect parameters) for about half of the 15 common as well as many other more infrequent species. CSI was thus not allowed to vary random in the final AS-REAL models. Incorporating extracted parameter estimates of the random effects into ΔDBH and ΔBA predictions substantially improved performance of the final AS-REAL models. MAB and RMSE decreased by 10 and 12% for the 15 most abundant and all other less common species, respectively, when comparing fixed effects only predictions with predictions considering fixed and random effects (data not shown). ΔDBH curves of the SS-REAL and AS-REAL approach (fixed and random effects) often displayed similar behavior, while the respective SS-POTMOD curves somewhat deviated at least for some of the 15 most abundant species (Fig. 1). Based on the derived evaluation metrics, comparing the performance of ΔDBH and ΔBA equations of the same approach indicated that modeling ΔDBH often resulted in greater prediction accuracy than modeling ΔBA (Fig. 2). In addition, the SS-REAL and specifically the ASREAL approaches often exhibited greater model fit than their SSPOTMOD counterparts. These patterns changed only slightly when evaluating trees ≥ 50 cm DBH (Table S2) or comparing soft- and hardwoods (Fig. 3) as well as species of varying shade tolerance (Fig. S1), with species of intermediate tolerance to shade exhibiting the lowest prediction accuracy. We also found no substantial differences between approaches across the various databases and thus corresponding DBH measurement thresholds used (Fig. S2). In contrast, differences in model fit between the two response variables (ΔDBH vs ΔBA) as well as across the various modeling approaches appeared to intensify when analyzing observations with a measurement period length ≥ 30 years (Fig. 4, Table S3). Consequently, ΔDBH equations for the 15 most abundant species received mostly higher relative ranks with the exception of MAB (SS-REAL, AS-REAL) and EI (SS-POTMOD) (Table 3). Averaging a value of 1.90, mean relative rank across all evaluation metrics was lowest for ΔDBH models derived in the AS-REAL

approach followed by ΔDBH equations of the SS-REAL and SS-POTMOD approach with average relative ranks of 2.07 and 2.10, respectively (Table 3). The AS-REAL approach performed best among the ΔBA equations with a mean relative rank of 2.52 closely followed by the SSREAL approach with 2.58. The striking difference in ranking among approaches for evaluation metric EI compared to all other metrics was a result of the superior prediction performance of the two SS-POTMOD approaches in trees of very large DBH (Table 3). Summing DBH classspecific MAB over the five largest 5 cm-DBH classes found in our data (80–105 cm) resulted in the same ranking of approaches as for EI, while the sum of MAB of all other smaller DBH classes exhibited the inverse ranking of approaches (data not shown). Consequently, the difference in ranking among approaches for evaluation metric EI was even more evident when evaluating trees ≥ 50 cm DBH only (Table S3). In contrast, there was no such apparent ranking difference in EI when only considering observations with a measurement period length ≥ 30 years because individuals > 70 cm DBH were missing in this evaluation (Table S2). Using the extracted species-specific parameter estimates of the random effects, the ΔDBH AS-REAL approach was robust to account for the diameter growth patterns of some species with a lower number of observations. Swamp white oak (Quercus bicolor Muenchh., N = 4) and pitch pine (Pinus rigida Mill., N = 185), for example, displayed ΔDBH curves characteristic for such medium-sized, shade-intolerant species by reaching a comparatively early peak (Fig. 5). In contrast, the latesuccessional, more shade-tolerant shagbark hickory (Carya ovata (Mill.) K. Koch, N = 13) and Norway spruce (Picea abies (L.) H.Karst, N = 68) showed a later peak in ΔDBH and/or a more moderate decline with increasing DBH (Fig. 5). However, not all rare species equations including butternut (N = 18) and Scots pine (Pinus sylvestris L., N = 30) were biologically plausible depicting either a generally uncharacteristic behavior or resulting in unrealistically large ΔDBH (Fig. 5). In addition, AS-REAL ΔDBH curves of a few species including black cherry, black spruce, and red pine displayed a unimodal convex behavior (Fig. 5), while combining fixed and random effect parameter estimates resulted in a negative parameter sign for explanatory variable CR in black cherry, northern red oak, and white ash (Table S2). 4. Discussion Stem diameter increment (ΔDBH) is a key component within a forest growth and yield model as predictions are passed through to other submodels and tree-level estimates are scaled up to represent plot- and stand-level measures. Evaluating various approaches to project DBH through time, this study showed that modeling realized ΔDBH with a unified model form was mostly superior when compared to approaches that estimate basal area increment (ΔBA). This was specifically evident for a two-stage procedure predicting maximum or potential growth (ΔBAPOT) and a corresponding modifier (ΔBAMOD) to derive the increment. Our study also showed the efficiency and flexibility of diameter predictions by including tree species as a random effect to account for ΔDBH differences of trees in mixed-species stands. In contrast to this study, previous analyses directly comparing performance of different equations for various species in other regions found no appreciable differences in the prediction of future DBH using ΔDBH or ΔBA models (West, 1980; Shifley, 1987). As stated by Russell et al. (2011), however, the works by West (1980) and Shifley (1987) were limited in terms of the projection lengths by only examining measurement intervals of 3 to 6 and 10 years, respectively. In contrast, Russell et al. (2011) found that ΔDBH models outperformed ΔBA equations in longer-term projections of > 25 years. Using measurement periods up to 40 years, this study also found a superior performance of ΔDBH equations compared to comparable ΔBA models. Moreover, the observed differences in prediction accuracy became more evident when evaluating observations with measurement intervals ≥ 30 years only. As stated previously, the use of ΔBA might have 7

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

Fig. 1. Annual diameter increment (ΔDBH, cm yr−1) predictions for six common Acadian tree species of varying shade tolerance over tree diameter (DBH, cm) for the average tree and stand conditions. Curves represent equations fitted with i) nonlinear mixed-effects models (NLME) that employed fixed plus species-specific random effects (AS-REAL), ii) generalized nonlinear least squares (GNLS) using a species-specific equation with a unified model form (SS-REAL), as well as iii) quantile regression (potential or maximum increment) and GNLS (modifier) for species-specific potential-modifier models (SS-POTMOD).

better model fit statistics when compared to ΔDBH, but the approach relies on a critical assumption of trees maintaining a non-eccentric and circular form through time. This assumption has a large influence as trees become larger, particularly hardwoods, which can become more eccentric with time. Results showed that the SS-REAL approach outperformed that of the SS-POTMOD for the 15 most common species. This was also shown by Subedi and Sharma (2011) when examining jack pine (Pinus banksiana Lamb.) and black spruce in Ontario. It is hypothesized that the potential increment accounted for a large amount of the variability in the performance of the SS-POTMOD approach. This likely was especially true for ΔBAPOT equations of this study as a second significant BA term that allows for a sigmoidal curve, i.e. a decline in ΔBA for larger DBH could not be incorporated. ΔDBHMOD and ΔBAMOD thus may not have reduced ΔDBHPOT and ΔBAPOT, respectively, enough to account for the true growth patterns. However and as revealed in the EI ranking for all data (Table 3), the SS-POTMOD approach worked particularly well for ΔBA predictions in trees of very large DBH where the general performance of the ΔBAMOD equation likely is of lower importance because of the superior competitive status of these dominant individuals. Moreover, procedures for estimating the uncertainty of ΔDBH predictions with the two SS-POTMOD stages remain unknown (Lessard et al., 2001), especially when techniques such as annualization are employed in model fitting. Quantile regression generally worked well and was able to smoothly derive ΔDBHPOT models (Pretzsch and Biber, 2010). We

believe that model convergence and stability issues encountered when fitting ΔBAPOT equations did not stem from the applied statistical procedure, but were a result of the structure of the underlying ΔBA data and the initially targeted model form. Despite the biological appeal and constrained behavior, the POTMOD approach likely constrains model behavior unnecessarily and complicates model fitting. Similar findings as this analysis were also identified by Russell et al. (2014) for height increment modeling. Difficulties associated with the POTMOD approach and how to potentially overcome them were also discussed by others (e.g. Lessard et al. (2001) and Pokharel and Froese (2008)). The AS-REAL approach was flexible in generating curves that were expected for the general growth pattern of most species, including some infrequent species. For example, 185 growth measurements were available for pitch pine, an exceedingly variable tree in terms of form and growth, often putting on diameter growth in installments, and at its northern range limit in Maine (Hardin et al., 2001). The NLME procedure provided a corresponding ΔDBH curve that was expected for the species, by reaching a peak at a small DBH and diminishing thereafter. To construct an equation depicting ΔDBH for these kinds of species, modelers would be faced with the difficulties of dealing with small sample sizes to fit an accurate growth curve for a target species. Similarly, data that are available for these species likely come from a limited number of plots that may not represent the variability in terms of climate and stand structures that the species is subject to. Similarly to what Weiskittel et al. (2007) observed after fitting NLME equations 8

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

Fig. 2. Mean bias (MB), relative MB (MB%), mean absolute bias (MAB), relative MAB (MAB%), root mean square error (RMSE), and error index (EI) for the evaluated annual diameter (ΔDBH, cm yr−1) and basal area increment (ΔBA, m2 yr−1) modeling approaches including single species potential-modifier (SS-POTMOD), single species realized (SS-REAL), and all species realized (AS-REAL).

estimates. In species-specific models, such outcomes can be avoided by excluding the specific explanatory variable from the equation. However, only very few species exhibited such a reversal of sign for the combined fixed plus random effect parameter estimate of CR. In contrast, species-specific random effects for the site productivity metric CSI resulted in overall unexpected negative parameter signs for over half of all common as well as infrequent species as shown by preliminary analysis. Since the fixed parameter estimate indicated a general and significantly positive effect on diameter growth, predictor CSI was retained in the ΔDBH and ΔBA AS-REAL models, respectively, but not allowed to vary random. Moreover, not all AS-REAL equations of infrequent species displayed the expected behavior often producing uncharacteristic increment values for specific DBH. We would suggest that sufficient data quality (i.e. measurements that reliably reflect the course of an increment curve to a minimum) appears to be still necessary to derive plausible models for rare species from the AS-REAL approach analyzed here. As initially assumed, however, data quantity (i.e. number of observations) seems of lower significance as compared to approaches aiming at species-specific growth equations. In contrast, species-specific model equations for rare species often do not converge as a result of the lack of sufficient data.

using installation and plot as random effects, improvement in model performance was noted when extracted random effects were incorporated into ΔDBH predictions, as indicated by lower MAB and RMSE for the 15 most abundant species investigated as well as for all less common species/genera. Results highlight the ability that population-level fixed-effects can have in accounting for plot-and tree-level differences in growth, but also how extracted random effects for species can provide robust predictions of ΔDBH for mixed-species stands. However, in comparison to the two approaches that derived speciesspecific equations, the AS-REAL approach also held specific shortcomings that should be noted and addressed in future applications. Due to the biologically implausible behavior of its estimated fixed effect parameter, RD had to be removed from the AS-REAL model. In contrast, various SS-POTMOD and SS-REAL equations included significantly negative parameters for RD in addition to BALSW and BALHW. Including both one- (BAL) and two-sided metrics of competition (RD) more effectively represents competition processes in growth equations and thus offers greater potential to improve prediction accuracy (Weiskittel et al., 2011b). Similarly, incorporating species-specific random effects can substantially alter predictor effects when plausible parameter signs are reversed as a result of summing fixed and random effects parameter 9

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

Fig. 3. Mean bias (MB), relative MB (MB%), mean absolute bias (MAB), relative MAB (MAB%), root mean square error (RMSE), and error index (EI) for the evaluated annual diameter (ΔDBH, cm yr−1) and basal area increment (ΔBA, m2 yr−1) modeling approaches including single species potential-modifier (SS-POTMOD), single species realized (SS-REAL), and all species realized (AS-REAL) by species type, namely hard- (HW, left column) and softwoods (SW, right column). Only trees < 90 cm DBH were analyzed for comparative purposes.

With the exception of RD, signs and magnitude of the fixed-effects parameters for the AS-REAL approach were generally what was expected for the growth of the species in the region (Table S1). The BALSW and BALHW terms were both significantly negative, indicating that competition from both softwood and hardwood species diminishes diameter growth. This was considered to be biologically appropriate, yet differed somewhat when compared to other studies. For example, Nunifu (2009) observed that an increase in BALHW favored diameter growth for white spruce and lodgepole pine (Pinus contorta Dougl. ex Loud. var. latifolia Engelm.) in Alberta. Results observed here are likely due to the fact that a greater number species occurred in these stand types and stands observed by Nunifu (2009) contained a smaller proportion of hardwood dominance compared to this study. ΔDBH curves of shade-intolerant and shade-tolerant species often display substantial differences. While the former generally exhibit an early peak at smaller DBH and a swift decline thereafter, the latter peak at greater DBH often at a lower magnitude followed by a more gradual, less pronounced decrease (e.g. Castle et al. 2018). Such distinctive

patterns are less evident in species of intermediate shade tolerance. Variation in the DBH-ΔDBH relationship across as well as within intermediate species thus might be a reason why secondary growth predictions of this study were least accurate for these species. An important future consideration will be in testing the ability that the extracted parameter estimates of the random effects for species (Table S1) can represent ΔDBH for new observations. Models of ΔDBH have commonly been designed that employ stand and/or plot as the random effect (e.g., Weiskittel et al., 2007; Fortin et al., 2008; Sudebi and Sharma, 2011). Calibration approaches using such stand and plot as well as time and/or tree random effect structures have been evaluated and appraised previously (e.g. Calama and Montero, 2005; Lhotka and Loewenstein, 2011). However, an assessment of how best to predict random effects for a new stand offers further investigation. This might include selecting a random tree of each species within a plot to improve the predictive ability for the ΔDBH equation when applied to these new stands. In addition, complications arise when using periodic measurements for calibration of annualized predictions or using tree ring

10

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

Fig. 4. Mean bias (MB), relative MB (MB%), mean absolute bias (MAB), relative MAB (MAB%), root mean square error (RMSE), and error index (EI) for the evaluated annual diameter (ΔDBH, cm yr−1) and basal area increment (ΔBA, m2 yr−1) modeling approaches including single species potential-modifier (SS-POTMOD), single species realized (SS-REAL), and all species realized (AS-REAL) for observations with a measurement period length ≥ 30 years. Table 3 Relative ranks after Poudel and Cao (2013) for the evaluated annual diameter (ΔDBH, cm yr−1) and basal area increment (ΔBA, m2 yr−1) modeling approaches including single species potential-modifier (SS-POTMOD), single species realized (SS-REAL), and all species realized (AS-REAL) by evaluation metric, namely mean bias (MB), relative MB (MB%), mean absolute bias (MAB), relative MAB (MAB%), root mean square error (RMSE), and error index (EI) after Reynolds et al. (1988). Approach

MB

MB%

MAB

MAB%

RMSE

EI

Overall

ΔDBH: SS-POTMOD ΔDBH: SS-REAL ΔDBH: AS-REAL ΔBA: SS-POTMOD ΔBA: SS-REAL ΔBA: AS-REAL

2.69 2.69 2.97 6 1 1.33

1 1.72 1.66 6 2.9 2.87

1.85 1 1.28 6 1.96 2.31

1 1.12 1.17 6 2.25 2.42

3.32 1 1.49 6 1.39 1.97

2.58 5.07 2.81 1 6 4.23

2.07 2.10 1.90 5.17 2.58 2.52

Fig. 5. Annual diameter increment (ΔDBH, cm yr−1) predictions for four infrequent Acadian tree species over tree diameter (DBH, cm) for the average tree and stand conditions. Curves were derived from the AS-REAL approach using nonlinear mixed-effects modeling (NLME) that employed fixed plus speciesspecific random effects.

observations from one or even two tree cores (Weiskittel et al., 2011b). 5. Conclusions For estimating species-specific growth patterns, this analysis found that modeling ΔDBH and ΔBA using a single, unified equation outperformed one which estimated potential increment and a multiplicative modifier in two stages. Incorporating species as a random effect using a nonlinear mixed-effect modeling (NLME) approach resulted in predictions that were equivalent or superior in terms of prediction accuracy when compared to species-specific unified equations. The ASREAL strategy shown here therefore can be useful for modelers seeking to capture the growth dynamics of individual trees found in complex stand structures with various species assemblages. However, NLME curves for rare species still need to be evaluated for biological plausibility as unbalanced or biased data can lead to uncharacteristic and likely illogical behavior. Although the AS-REAL offers a simple and

effective approach for quickly determining genus- and/or species-level predictions, there is an important cost and it may be a less than perfect solution, so widespread application without through assessment is not suggested. Despite the general findings, the derived increment equations of this analysis are likely not the most accurate ones possible given the fixed set of predefined explanatory variables, which were used to maximize comparability across all studied modeling approaches. Given the importance of DBH in growth and yield models, an accurate approach to predict future DBH is required. Despite the importance, a variety of approaches have been used in the past and longterm assessments are relatively rare. Clear differences of these 11

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

contrasting approaches were evident across species, tree size, and projection length. In particular, divergence across methods was observed as projection length increased with 40 years being the maximum length available for assessment in this analysis, while most operational growth and yield projections are 50–100 years or even longer in length (Weiskittel et al., 2011b). Of course, continued evaluation of these approaches across various species and forest types is recommended, particularly the use of species as a random effect. The general balance and number of observations across species as well as inclusion of additional hierarchies such as genus (e.g. Lam et al., 2016), forest type, or ecoregion will likely influence the overall effectiveness of this method. In addition, there are additional methodologies for predicting future DBH such as self-referencing (e.g. Sharma et al., 2017) or multi-dimensional system approaches (Garcia, 2018) that should also be examined. Overall, the strategies employed in this analysis offer a general framework for predicting future DBH for mixed species stands and should be further examined in additional regions with this forest type.

and carbon of the highly diverse Atlantic Forest in Brazil: comparison of alternative individual tree modeling and prediction strategies. Carbon Manage. 9 (4), 383-397. Dixon, G.E., Keyser, C.E., 2018. Northeast (NE) variant overview-Forest Vegetation Simulator (revised March 16, 2018). Internal Rep. Fort Collins, CO: USDA For. Serv., For. Manage. Serv. Cen. Fort Collins, CO. 40 p. Fortin, M., Bédard, S., DeBlois, J., Meunier, S., 2008. Accounting for error correlations in diameter increment modelling: a case study applied to northern hardwood stands in Quebec, Canada. Can. J. For. Res. 38 (8), 2274–2286. Garcia, O., 2018. Reverse causality in size-dependent growth. Math. Comput. For. Nat. Resour. Sci. 10 (1), 1–5. Hann, D.W., Larsen, D.R., 1991. Diameter growth equations for fourteen tree species in southwest Oregon. Oregon State Univ. For. Res. Lab. Res. Bull. 69. Hann, D.W., Marshall, D.D., Hanus, M.L., 2003. Equations for predicting height-to-crownbase, 5-year diameter-growth rate, 5-year height-growth rate, 5-year mortality rate, and maximum size-density trajectory for Douglas-fir and western hemlock in the coastal region of the Pacific Northwest. Oregon State Univ. For. Res. Lab. Res. Cont. 40. Hardin, J.W., Leopold, D.J., White, F.M., 2001. Harlow & Harrar's textbook of dendrology, 9th ed. McGraw Hill. 534 pp. Huang, J.-G., Stadt, K.J., Dawson, A., Comeau, P.G., 2013. Modelling growth-competition relationships in trembling aspen and white spruce mixed boreal forests of Western Canada. PLoS One 8, e77607. Kangas, A.S., 1999. Methods for assessing uncertainty of growth and yield predictions. Can. J. For. Res. 29 (9), 1357–1364. Kenefic, L.S., Rogers, N.S., Puhlick, J.J., Waskiewicz, J.D., Brissette, J.C., 2015. Overstory tree and regeneration data from the “Silvicultural Effects on Composition, Structure, and Growth” study at Penobscot Experimental Forest. 2nd Edition. Fort Collins, CO: Forest Service Research Data Archive. https://doi.org/10.2737/RDS-2012-0008-2 (Accessed 23 September 2019). Kmenta, J., 1997. Elements of econometrics. Ann Arbor, MI: University of Michigan Press. Koenker, R., 2013. quantreg: Quantile regression. R package version 4.98. Available from http://www.r-project.org. Koenker, R., Hallock, K.F., 2001. Quantile regression. J. Econom. Persp. 15 (4), 143–156. Kuehne, C., Puhlick, J.J., Weiskittel, A.R., 2018a. Ecological reserves in Maine: initial results of long-term monitoring. Orono (ME): University of Maine, Center for Research on Sustainable Forests. [accessed 2018 December 6]. http://www. nefismembers.org/documents/ecological-reserves-in-maine-initial-results-of-longterm-monitoring. Kuehne, C., Puhlick, J., Weiskittel, A., Cutko, A., Cameron, D., Sferra, N., Schlawin, J., 2018b. Metrics for comparing stand structure and dynamics between ecological reserves and managed forest of Maine, USA. Ecology 99 (12), 2876. Kuehne, C., Weiskittel, A., Pommerening, A., Wagner, R.G., 2018c. Evaluation of 10-year temporal and spatial variability in structure and growth across contrasting commercial thinning treatments in spruce-fir forests of northern Maine, USA. Ann. For. Sci. 75, 20. Kuehne, C., Weiskittel, A.R., Wagner, R.G., Roth, B.E., 2016. Development and evaluation of individual tree-and stand-level approaches for predicting spruce-fir response to commercial thinning in Maine. USA. For. Ecol. Manage. 376, 84–95. Kuehne, C., Weiskittel, A.R., Waskiewicz, J., 2019. Comparing performance of contrasting distance-independent and distance-dependent competition metrics in predicting individual tree diameter increment and survival within structurally-heterogeneous, mixed-species forests of Northeastern United States. For. Ecol. Manage. 433, 205–216. Lam, T.Y., Kershaw, J.A., Hajar, Z.S.N., Rahman, K.A., Weiskittel, A.R., Potts, M.D., 2016. Evaluating and modelling genus and species variation in height-to-diameter relationships for Tropical Hill Forests in Peninsular Malaysia. Forestry 90 (2), 268–278. Larocque, G.R., Archambault, L., Delisle, C., 2011. Development of the gap model ZELIGCFS to predict the dynamics of North American mixed forest types with complex structures. Ecol. Model. 222, 2570–2583. Lessard, V.C., McRoberts, R.E., Holdaway, M.R., 2001. Diameter growth models using Minnesota forest inventory and analysis data. For. Sci. 47 (3), 301–310. Li, R., Weiskittel, A.R., Kershaw Jr., J.A., 2011. Modeling annualized occurrence, frequency, and composition of ingrowth using mixed-effects zero-inflated models and permanent plots in the Acadian Forest region of North America. Can. J. For. Res. 41, 2077–2089. Lhotka, J.M., Loewenstein, E.F., 2011. An individual-tree diameter growth model for managed uneven-aged oak-shortleaf pine stands in the Ozark Highlands of Missouri, USA. For. Ecol. Manage. 261 (3), 770–778. MacFarlane, D.W., Kobe, R.K., 2006. Selecting models for capturing tree-size effects on growth resource relationships. Can. J. For. Res. 36 (7), 1695–1704. McGarrigle, E., Kerhsaw Jr., J.A., Lavigne, M.B., Weiskittel, A.R., Ducey, M., 2011. Predicting the number of trees in small diameter classes using predictions from a twoparameter Weibull distribution. Forestry 84 (4), 431–439. McRoberts, R.E., Bechtold, W.A., Patterson, P.L., Scott, C.T., Reams, G.A., 2005. The enhanced forest inventory and analysis program of the USDA forest service: historical perspective and announcement of statistical documentation. J. For. 103 (6), 304–308. Monserud, R.A., Sterba, H., 1996. A basal area increment model for individual trees growing in even- and uneven-aged forest stands in Austria. For. Ecol. Manage. 80 (1–3), 57–80. Moore, A.D., 1989. On the maximum growth equation used in forest gap simulation models. Ecol. Model. 45, 63–67. Niinemets, Ü., Valladares, F., 2006. Tolerance to shade, drought, and waterlogging of temperate Northern Hemisphere trees and shrubs. Ecol. Monogr. 76, 521–547. Nunifu, T.K., 2009. Compatible diameter and height increment models for lodgepole pine, trembling aspen, and white spruce. Can. J. For. Res. 39 (1), 180–192. Phillips, P.D., Yasman, I., Brash, T.E., van Gardingen, P.R., 2002. Grouping tree species

CRediT authorship contribution statement C. Kuehne: Methodology, Data curation, Formal analysis, Writing original draft. M.B. Russell: Conceptualization, Methodology, Writing original draft, Writing - review & editing. A.R. Weiskittel: Conceptualization, Funding acquisition, Writing - review & editing. J.A. Kershaw: Data curation, Writing - review & editing. Acknowledgments Data for this analysis was provided by numerous organizations including Maine Department of Agriculture Conservation and Forestry, New Brunswick Department of Natural Resources, Nova Scotia Department of Natural Resources, Penobscot Experimental Forest, Quebec Ministry of Forests, The Nature Conservancy Maine, and US Forest Service. Funding was provided by National Science Foundation Center for Advanced Forestry Systems (CAFS) and USDA National Institute of Food and Agriculture, McIntire-Stennis Project Number ME041516 through the Maine Agricultural and Forest Experimental Station. Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References Bataineh, M., Wagner, R.G., Weiskittel, A.R., 2013. Long-term response of spruce-fir stands to herbicide and precommercial thinning: Observed and projected growth, yield, and financial returns in central Maine, USA. Can. J. For. Res 43, 385–395. Bechtold, W.A., Patterson, P.L. (eds)., 2005. Forest Inventory and Analysis national sample design and estimation procedures. USDA For. Serv. Gen. Tech. Rep. SRSGTR-80. Bowman, D.M.J.S., Brienen, R.J.W., Gloor, E., Phillips, O.L., Prior, L.D., 2013. Detecting trends in tree growth: not so simple. Trends Plant Sci. 18 (1), 11–17. Braun, E.L., 1950. Deciduous forests of eastern North America. Hafner, New York, NY. 596 pp. Calama, R., Montero, G., 2005. Multilevel linear mixed model for tree diameter increment in stone pine (Pinus pinea): a calibrating approach. Silva Fenn 39 (1), 37–54. Cao, Q.V., 2000. Prediction of annual diameter growth and survival for individual trees from periodic measurements. For. Sci. 46 (1), 127–131. Castle, M., Weiskittel, A., Wagner, R., Ducey, M., Frank, J., Pelletier, G., 2018. Evaluating the influence of stem form and damage on individual-tree diameter increment and survival in the Acadian Region: implications for predicting future value of northern commercial hardwood stands. Can. J. For. Res. 48, 1007–1019. Chen, C., Weiskittel, A., Bataineh, M., MacLean, D.A., 2017. Evaluating the influence of varying levels of spruce budworm defoliation on annualized individual tree growth and mortality in Maine, USA and New Brunswick, Canada. For. Ecol. Manage. 396, 184–194. Cole, D.M., Stage, A.R., 1972. Estimating future diameter of lodgepole pine trees. Research Paper INT-131. U.S. Department of Agriculture, Forest Service. Intermountain Forest and Range Experiment Station, Ogden, UT, 20 pp. Colmanetti, M.A.A., Weiskittel, A., Barbosa, L.M., Shirasuna, R.T., de Lima, F.C., Ortiz, P. R.T., Catharino, E.L.M., Barbosa, T.C., do Couto, H.T.Z., 2018. Aboveground biomass

12

Forest Ecology and Management 459 (2020) 117823

C. Kuehne, et al.

tropical rain-forests. For. Ecol. Manage. 42 (3–4), 143–168. Vanclay, J.K., 1995. Growth models for tropical forests - a synthesis of models and methods. For. Sci. 41 (1), 7–42. Wagner, R.G., Seymour, R.S., 2006. Commercial thinning research network. In: 20042005 Cooperative Forestry Research Unit Annual Report. Univ. of Maine, Orono, ME. pp. 34–41. Weiskittel, A.R., Crookston, N.L., Radtke, P.J., 2011a. Linking climate, gross primary productivity, and site index across forests of the western United States. Can. J. For. Res. 41, 1710–1721. Weiskittel, A., Frank, J., Walker, D., Radtke, P., Macfarlane, D. Westfall, J., 2015. Advancing individual tree biomass prediction: assessment and alternatives to the component ratio method. In: Stanton, Sharon M.; Christensen, Glenn A., comps. 2015. Pushing boundaries: new directions in inventory techniques and applications: Forest Inventory and Analysis (FIA) symposium 2015. 2015 December 8–10; Portland, Oregon. Gen. Tech. Rep. PNW-GTR-931. US Department of Agriculture, Forest Service, Pacific Northwest Research Station, Portland, OR, 125–132. Weiskittel, A.R., Garber, S.M., Johnson, G.P., Maguire, D.A., Monserud, R.A., 2007. Annualized diameter and height growth equations for pacific northwest plantationgrown Douglas-fir, western hemlock, and red alder. For. Ecol. Manage. 250 (3), 266–278. Weiskittel, A. R., Hann, D. W., Kershaw Jr, J. A., & Vanclay, J. K. 2011b. Forest growth and yield modeling. John Wiley & Sons. Weiskittel AR, Kuehne C. Evaluating and modeling variation in site-level maximum carrying capacity of mixed-species forest stands in the Acadian Region of northeastern North America. For Chron. In press. Weiskittel, A., Kuehne, C., McTague, J.P., Oppenheimer, M., 2016. Development and evaluation of an individual tree growth and yield model for the mixed species forest of the Adirondacks Region of New York, USA. For. Ecosyst. 3 (1), 26. Weiskittel, A.R., Wagner, R.G., Seymour, R.S., 2010. Refinement of the Forest Vegetation Simulator, northeastern variant growth and yield model: Phase 1. In: 2009 Cooperative Forestry Research Unit Annual Report (S.R. Meyer, ed.). Univ. of Maine, Orono, ME. pp. 44–48. Weiskittel, A.R., Wagner, R.G., Seymour, R.S., 2011c. Refinement of the forest vegetation simulator, northeastern Variant growth and yield model: phase 2. In: Mercier W.J. and Nelson A.S. (eds.) 2010 Cooperative Forestry Research Unit Annual Report. University of Maine, School of Forest Resources, Orono, pp. 23–30. Wensel, L.C., Meerschaert, W.J., Biging, G.S., 1987. Tree height and diameter growth model for northern California coast. Hilgardia 55, 1–20. West, P.W., 1980. Use of diameter increment and basal area increment in tree growthstudies. Can. J. For. Res. 10 (1), 71–77. Wilson, D., Monleon, V., Weiskittel, A., 2019. Quantification and incorporation of uncertainty in forest growth and yield projections using a Bayesian probabilistic framework: a demonstration for plantation coastal Douglas-fir in the Pacific Northwest, USA. Math. Comput. For. Nat.-Resource Sci. 11 (2), 264–285. Wykoff, W.R., 1990. A basal area increment model for individual conifers in the northern Rocky Mountains. For. Sci. 36, 1077–1104. Zhao, D., Borders, B., Wilson, M., Rathbun, S.L., 2006. Modeling neighborhood effects on the growth and survival of individual trees in a natural temperate species-rich forest. Ecol. Model. 196, 90–102.

for analysis of forest data in Kalimantan (Indonesian Borneo). For. Ecol. Manage. 157, 205–216. Picard, N., Mortier, F., Rossi, V., Gourlet-Fleury, S., 2010. Clustering species using a model of population dynamics and aggregation theory. Ecol. Model. 221, 152–160. Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., 2012. nlme: Linear and nonlinear mixed effects models. R package version 3.1-103. Available from http://www.r-project.org. Pogoda, P., Ocha, W., Orze, S., 2019. Modeling diameter distribution of black alder (Alnus glutinosa (L.) Gaertn.) stands in Poland. Forests 10, 412. Pokharel, B., Froese E., R., 2008. Evaluating alternative implementations of the Lake States FVS diameter increment model. For. Ecol. Manage. 255 (5–6), 1759–1771. Poudel, K.P., Cao, Q.V., 2013. Evaluation of methods to predict Weibull parameters for characterizing diameter distributions. For. Sci. 59 (2), 243–252. Pretzsch, H., Biber, P., 2010. Size-symmetric versus size-asymmetric competition and growth partitioning among trees in forest stands along an ecological gradient in central Europe. Can. J. For. Res. 40 (2), 370–384. Province of New Brunswick. 2005. Partial harvest permanent sample plot establishment manual. Internal report. 27 pp. R Development Core Team. 2018. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Available from http://www.r-project.org. Rehfeldt, G.E., 2006. A spline model of climate for the western United States. USDA For. Serv. Gen. Tech. Rep. RMRS-165. Reynolds, M.R., Burk, T.E., Huang, W.C., 1988. Goodness-of-fit tests and model selection procedures for diameter distribution models. For. Sci. 34 (2), 373–399. Rowe, J.S., 1972. Forest regions of Canada. Publ. 1300. Dept. of the Env., Can. For. Serv., 251 Ottawa. 172 pp. Russell, M.B., Weiskittel, A.R., Kershaw Jr., J.A., 2011. Assessing model performance in forecasting long-term individual tree diameter versus basal area increment for the primary Acadian tree species. Can. J. For. Res. 41, 2267–2275. Russell, M., Weiskittel, A., Kershaw Jr., J., 2013. Benchmarking and calibration of Forest Vegetation Simulator individual tree attribute predictions across the northeastern United States. Northern J. Appl. For. 30 (2), 75–84. https://doi.org/10.5849/njaf.12034. Russell, M.B., Weiskittel, A.R., Kershaw, J.A., 2014. Comparing strategies for modeling individual-tree height and height-to-crown base increment in mixed-species Acadian Forests of northeastern North America. Eur. J. For. Res. 133, 1121–1135. Sharma, R.P., Vacek, Z., Vacek, S., Jansa, V., 2017. Modelling individual tree diameter growth for Norway spruce in the Czech Republic using a generalized algebraic difference approach. J. For. Sci. 63 (5), 227–238. Shifley, S.R., 1987. A generalized system of models forecasting Central States tree growth. USDA For. Serv. Res. Pap. NC-279. Strimbu, V., Bokalo, M., Comeau, P., 2017. Deterministic models of growth and mortality for jack pine in boreal forests of western Canada. Forests 8 (11), 410. Subedi, N., Sharma, M., 2011. Individual-tree diameter growth models for black spruce and jack pine plantations in northern Ontario. For. Ecol. Manage. 261 (11), 2140–2148. Teck, R.M., Hilt, D.E., 1991. Individual-tree diameter growth model for the northeastern United States. USDA For. Serv. Res. Pap. NE-649. 11 pp. Vanclay, J.K., 1991. Aggregating tree species to develop diameter increment equations for

13