Development and Evaluation of Black Spruce ( Picea mariana (Miller) B.S.P.) Diameter Increment Models across Silvicultural Treatments in Northern Minnesota, USA

: The black spruce cover type occupies roughly 10% of Minnesota’s 7 million hectares of forestland, and is an important species, both ecologically and economically. A clearcut regeneration harvest is the main silvicultural system in black spruce in this region. The effects of managing black spruce with alternative silvicultural methods in the Lake States remains largely understudied. Here, we examine a silviculture study in lowland black spruce to assess the performance of two diameter growth models ﬁt to this data compared to a widely-used model. Six silvicultural treatments (clearcut strips, clearcut patches, thinning, group selection, single-tree selection, and shelterwood) and a control were treated and measured around 1950, with a follow-up measurement occurring 10 years later. Fixed- and mixed-effects growth-models were adapted from the previous work, and ﬁt to 10,231 observations and compared to a recently released diameter growth model. The mixed-effects model using treatment, compartment, and plot as nested random effects outperformed the ﬁxed-effects model, and outperformed a model proposed for use in the Lake States variant of the Forest Vegetation Simulator that was ﬁt to this data. This modeling approach of localized growth models across a wide-range of diameters (9.1–32.1 cm) more accurately predicted the diameter growth in lowland black spruce than the conventional approach of using separate models for large (>12.7 cm) and small ( ≤ 12.7 cm) diameter trees.


Introduction
Black spruce (Picea mariana (Miller) B.S.P.) is distributed broadly across North American boreal forests, and is an important species, both economically and ecologically [1]. Black spruce is the most important pulpwood species in Canada, and it is a major commercial species in the Lake States region (Michigan, Minnesota, and Wisconsin) of the United States [1]. In Minnesota, black spruce comprises 648,000 hectares of the 7.04 million forested hectares [2] and is the second-most harvested pulpwood species by volume [3].
Across the US Lakes States, black spruce is near its southern and western range limit. Black spruce can occur under a variety of environmental conditions in both pure and mixed composition stands. Black spruce is the dominant overstory species in the black spruce/feather moss-lichen forest type, and is a co-dominant species in five other forest types [4]. In Minnesota, black spruce can occur on

Study Area
The Compartment Study is a long-term experiment located on the Big Falls Experimental Forest (BFEF) near Big Falls, MN, USA (48 • 21 N, 93 • 46 W, 371 m above sea level.). Minnesota is one of three states in the Lake States region, an area near the western Great Lakes and south of the border with Canada. The study occupies the lowland, or peatland, black spruce cover type, black spruce group (12) [4], and Native Plant Community type APn81 [26]. The climate is continental with short, warm summers and long, cold winters. Maximum summer temperatures can exceed 32 • C, with high humidity (80 percent), and minimum winter temperatures can reach −35 • C. Precipitation ranges from 500-640 mm, with average snowfall around 150 cm [27].
The BFEF was established in 1948 by an agreement between the Minnesota Department of Conservation (now the Minnesota Department of Natural Resources) and the United States Forest Service (USFS) Lake States Forest Experiment Station (now the USFS Northern Research Station). The agreement reserved 828 hectares of state-owned land for the purpose of silviculture, harvesting, and utilization, and economic studies in lowland conifers. The cover type at the time was being exploited and severely overcut [28]. The original goals of the Compartment Study were to examine the extension of the life and improvement in the productivity of these forests, and also included an investigation of regeneration after harvest, gathering information on site conditions and site determination, and studying the economic viability of various management practices in lowland black spruce [29,30].

Field Methods
Six silvicultural treatments and a control were implemented in 1948: Three even-aged treatments (clearcut strips, clearcut patches, shelterwood), two uneven-aged treatments (group selection, single-tree selection), and one intermediate treatment (light thinning). Treatments were carried out in a single winter between 1949 and 1951. The clearcut strips treatment consisted of 20-to 40 m wide strips oriented north-south to favor seed dispersal from westerly prevailing winds. Clearcut patches were 0.1 to 0.2 ha patches of any shape ( Figure 1). Three or four entries over 17 years were utilized to cut the entire stand into both clearcut treatments. The shelterwood treatment was implemented with a shelterwood cut to an average basal area of 11.5 m 2 ha −1 , leaving healthy and vigorous trees. The overstory removal cut was implemented 10 years after the shelterwood cut, with care taken to avoid damage to the understory. The single-tree and group selection cuts targeted individual trees or groups of four to eight trees, with openings not larger than nine to 12 m in diameter, with the goal of develop all-age conditions through variable cutting cycles. The thinning treatment was a free thinning to capture mortality through a variable cutting cycle of every eight to 10 years. The ideal scenario was that the light thinning treatment would be followed by a regeneration harvest; however, ultimately only the initial treatment was implemented.
Each treatment was replicated three times, totaling 21 compartments. Depending on the treatment, compartments contained eight to 16 permanent 0.4 ha plots, and ranged in size from 2.3 to 3.8 ha. Compartments were all located within two kilometers of each other (Figure 2). At the time of establishment, the 21 compartments ranged from 66 to 176 years old, stand structure was even-or uneven-aged, and site index (total height of dominant and co-dominant trees at 100 years) ranged from 9.1 to 18.3 m. determined by coring two or three black spruce nearest to the plot center at breast height. Individual tree variables included species, diameter to 0.1 cm (all trees over 9.1 cm in DBH), crown class, crown vigor, and form or health descriptions of each tree. Individual trees within each plot were stemmapped and numbered, and variables were recorded for specific trees with each measurement. Repeated measurements on tagged trees allowed for the diameter and other variables to be monitored over time.

Analysis
Ten-year DBH growth was calculated as the difference between DBH at measurement one (around 1950) and DBH at measurement three (around 1960). Some compartments were re-measured nine or 11 years after the initial measurement. In such instances, 10-year DBH change was standardized by calculating the difference in DBH from measurement one to measurement three, dividing by the number of years between measurements, and multiplying by 10. Plot-level basal area was calculated by expanding to a per hectare basis.

Ten-year Diameter Growth Model
A series of diameter growth models was formulated to predict 10-year diameter increment from the initial to the 10-year measurement. A model from Hann et al. [22] was adapted, as it has proven to be successful in modeling the diameter growth in other species [31], including black spruce [23]. The initial measurement was taken prior to the initial treatment in all compartments between 1948 and 1950. Two additional measurements were taken-one about five years after the initial treatment, and another about 10 years after the initial treatment (referred to as the "10-year re-measurement"), just prior to the second treatment. All compartments were treated both initially and at 10 years, except for the control (which was left untreated) and the light thinning, which only received the initial treatment.
In the clearcut treatments, only plots that were not cut in the original treatment were re-measured, due to the immature size class of treated plots at the time of re-measurement. These uncut portions were removed in later treatments after the 10-year re-measurement, when strips and patches were expanded. Residual plot centers were as close as 15 m from areas that were removed with the initial treatment, meaning that trees on the outer edge of plots may have been within several meters from cut portions of the stand ( Figure 1). As such, they are considered to be different from the controls, which had no portions of any plot come within the same proximity to harvesting.
Plot-level variables recorded for each plot were site index and age for each plot. Age was determined by coring two or three black spruce nearest to the plot center at breast height. Individual tree variables included species, diameter to 0.1 cm (all trees over 9.1 cm in DBH), crown class, crown vigor, and form or health descriptions of each tree. Individual trees within each plot were stem-mapped and numbered, and variables were recorded for specific trees with each measurement. Repeated measurements on tagged trees allowed for the diameter and other variables to be monitored over time. A link to the data used can be found in the Supplementary Materials.

Analysis
Ten-year DBH growth was calculated as the difference between DBH at measurement one (around 1950) and DBH at measurement three (around 1960). Some compartments were re-measured nine or 11 years after the initial measurement. In such instances, 10-year DBH change was standardized by calculating the difference in DBH from measurement one to measurement three, dividing by the number of years between measurements, and multiplying by 10. Plot-level basal area was calculated by expanding to a per hectare basis.

Ten-Year Diameter Growth Model
A series of diameter growth models was formulated to predict 10-year diameter increment from the initial to the 10-year measurement. A model from Hann et al. [22] was adapted, as it has proven to be successful in modeling the diameter growth in other species [31], including black spruce [23].
Crown position was originally recorded for each tree on a scale of seven classes: head dominants, strong dominants, conditional dominants and codominants, weak dominants and codominants, intermediates, suppressed trees, and open grown trees. To align with the more commonly used indicators, these were converted to four crown classes that are often used, by using descriptions of the original classes in the original working plan [30]: dominant, codominant, intermediate, and suppressed. The revised four-class crown position variable was then fit with a dummy variable where dominant and codominant trees were given a value of 1, and where intermediate and suppressed trees were 0.
The Hann et al. [22] model includes the crown ratio as a predictor variable; however, crown ratio was not measured in the Compartment Study. Instead, the revised four-class crown position variable mentioned above was fit with the dummy variable for dominant/codominant and intermediate/suppressed trees. Predicting the crown ratio was considered instead of using a crown position indicator variable, but this was ultimately rejected in favor of fitting a model using only variables that were measured in the study. Additionally, age was taken at the plot level in the Compartment Study, and this was assessed as a possible predictor variable. The possible full equation based on all assessed variables was: where ∆DBH is 10-year diameter increment in cm, DBH is the diameter of the tree at initial measurement (cm), SI is the site index in m at 100 years, BAL is an individual tree asymmetric-based variable equaling the basal area of trees larger than the subject tree in a plot, BA is the stand basal area of the initial stand after all initial treatments, CP is crown class dummy variable equal to 1 for dominant and codominant trees, and 0 for all other crown classes., and ε a is the within-tree error.

Model Fitting
All model fitting was completed using R [32] and the NLME package [33]. Model fitting began by using only fixed effects as a generalized nonlinear least squares (GNLS) model (Equation (1)). Variables were added one at a time from the form suggested by Weiskittel et al. [31], and they were kept in the model by assessing how the variable affected the Akaike information criterion (AIC) value for the model, and whether the coefficient for the variable was significant (p < 0.05). AIC is defined as AIC = −2 ln(L) + 2k where ln is the natural logarithm, L is the likelihood function, and k is the number of parameters in the model. AIC favors a small residual error, but penalizes additional predictors that over-fit a model. In our case, a variable would stay in the model if the coefficient was significantly different than zero, and if the AIC of that model was lowered by more than 10 points [34]. Model fitting was complete once no additional variables remained, or if additional variables did not meaningfully reduce the AIC. Root mean square error ( and assessed to compare model performance. After the model with only fixed effects was established, random effects using treatment, compartment nested within treatment, and plot nested within compartment nested within treatment were considered as a non-linear mixed effects model (NLME). The GNLS and NLME models were compared by assessing the difference in the AIC of adding random effects, R 2 , mean bias, and the RMSE of predicted diameter growth. Treatment is added as a random effect, due to it being a discrete variable with seven levels (six treatments and a control).
As Weiskittel et al. [31] pointed out, when the assumption of constant variance is violated, a power variance function of the initial diameter can be included, which gives less influence to trees with larger initial diameters. This was used in both the GNLS and NLME models. The power variance function is defined as s 2 (v) = |v| 2t , where s 2 v is the variance function evaluated at v, v is the variance covariate, and t is the variance function coefficient [31,32]. An assessment of the residual plot of fitted values confirmed no trends in the residuals.
The NLME model was also compared to the 10-year diameter growth model for black spruce over 12.7 cm DBH, developed by Deo & Froese [25] as a part of their redevelopment of the equations for large trees used in the Lake States Variant of FVS. Discrepancies currently exist over proper coefficients used in the work of Deo & Froese, and those implemented into FVS [35]. Therefore, the entire set of possible predictor variables used by Deo & Froese to model the diameter increment for all Lake States species were considered, and the resultant re-specified model was fit to the data from the Compartment Study. Using forward stepwise regression and the 'leaps' package [36], variables were selected that resulted in the best fitting model using the entire Compartment Study dataset, which included tree diameters as small as 9.1 cm, instead of the 12.7 cm minimum used in the original modeling by Deo & Froese. To overcome the uncompacted crown ratio (total height of crown from crown base to total tree height) not being measured in the Compartment Study, it was predicted using the equation in FVS described by Dixon and Keyser [24].
The full range of diameters in the dataset (12.7 cm and greater) were used in comparing the model performances of the GNLS, NLME, and refit Deo & Froese models. Fit statistics (RMSE and mean bias) were compared to assess the better fitting model. While the equations for other species have been implemented into FVS, the black spruce equation has yet to be implemented. Therefore, comparisons will be made to the refit Deo & Froese model, and referred to as such, rather than the FVS model.

Results
A large range of stand conditions at the plot level-both before (Table 1) and after ( Table 2) the initial treatment-and individual tree level (Table 3) existed in each treatment. However, stand conditions prior to initial treatment were not significantly different among treatments (degree of freedom, df = 6, p > 0.05).   Model fitting for the GNLS and NLME equations resulted in six fixed effects, as shown in Equation (2).
After fitting the GNLS model, nested random effects were added in fitting the NLME model. Plot nested within the compartment nested within the treatment resulted in the best model fit (AIC: 19,258) (Equation (3)). The final NLME model improved the coefficient of determination (R 2 ) by 51% from 0.2056 to 0.3114 over the GNLS model. The RMSE decreased 10% from 0.6772 to 0.6104 cm. Parameter estimates and their standard errors are in Table 4. Fit statistics for the final 10-year diameter growth models are in Tables 5 and 6. The final mixed effects model form was: where b ijk is a unique intercept term added to β 0 for each plot nested within the compartment and treatment. Refitting the model of Deo & Froese resulted in all but two of their independent variables being used to predict diameter increment (Equation (4)). The variables DBH QMD and DBH 2 QMD were not significant to the model: The NLME model outperformed the refit model of Deo & Froese, resulting in an R 2 value that was 230% higher (0.3113 vs. 0.1352) and a lower RMSE (0.6104 vs. 0.6849). The GNLS model saw mixed results, with an adjusted R 2 value of 0.2047 and a lower RMSE than the refit Deo & Froese model (0.6773).
The NLME model was much better in predicting the diameter growth over the refit Deo & Froese model across all treatments ( Table 5). The overall mean bias of the NLME model was 0.0054 cm, compared with −0.1469 cm for the Deo & Froese refit model. The Deo & Froese refit model underestimated diameter growth consistently across all treatments. The mean of the NLME model-predicted diameter growth was nearly identical to the observed growth, at 1.536 cm, compared to the observed growth of 1.539 cm over 10 years, while the Deo & Froese refit model averaged 1.357 cm ( Table 5).
For individual treatments, RMSE for the NLME model was the smallest in the two clearcut (0.5932 and 0.5490 cm) and thinning treatments (0.5697), and largest in the shelterwood treatment (0.8345 cm). The mean bias was the smallest in the group selection method (−0.0018 cm) and largest in the thinning treatment. The model overpredicted diameter growth slightly in the shelterwood and clearcut methods, and underpredicted diameter growth in the selection methods and the thinning ( Table 5).

Discussion
Ten-year diameter growth varied between treatments. The shelterwood treatment had the greatest diameter growth, likely a factor of it benefiting from increased available resources, due to having the lowest residual BA and TPA of all treatments after the initial treatment ( Table 3). The higher intensity of cutting resulted in greater growing space for residual trees in this treatment. Both selection treatments had the next highest average diameter growth, and they correspondingly had the next lowest residual BA. This pattern continued for the thinning treatment. Notably, plots that were re-measured in 1960 in the clearcut treatments were those that were not harvested in the initial 1950 treatment. Plots in patches or strips that were harvested were not re-measured, due to the immature size class of regeneration. The harvest area during the initial treatment in the clearcut strips or patches may have been as close as a few meters from trees on the outer portions of unharvested plots, based on the visual evidence of treatments and the cutting cycles of the experiment (Figure 1). Girona et al. [8] documented that trees at the edge of skidding trails showed significantly increased growth compared to trees in the interior of a stand. In future studies, location relative to harvest strips and patches could be assessed for their impact on diameter growth.
We adapted the Hann [22] model in as simple a matter as possible using two frameworks. First, we manually added basic variables (e.g., initial diameter), and progressed by adding more complex variables (Table 4). Second, we used AIC to balance the model performance with the model complexity.
Fitting the diameter growth model by adapting the model of Hann [22] and manually adding variables resulted in fit statistics (R 2 and RMSE) that were comparable to other modeling studies of diameter growth [22,23,25,32]. Using the power variance structure by weighting the variance for trees with larger diameters was effective in overcoming the non-constant variance to reduce the bias introduced by large trees. Stand age was tested as a possible independent variable, but this did not meaningfully reduce AIC. Age was not a factor in previous fittings using this model [22,32], or in Subedi and Sharma's [23] model fitting for black spruce in Ontario, Canada. Additionally, including age as a predictor variable may limit the model's applicability to stands where age is not measured or are uneven-aged.
Adding random effects (i.e., the NLME model) after fitting the model with all fixed effects (i.e., the GNLS model) was successful in improving the model performance, likely due to the variety of stand conditions in lowland black spruce. The method of first adding the treatment, and then the compartment nested within treatment, and finally the plot nested within the compartment nested within the treatment, was logical in terms of sequentially increasing the model complexity. Adding the compartment as a nested random effect with treatment was particularly effective in improving the model performance through reducing AIC. Compartments ranged in size from 2.3 to 3.8 hectares, and they were located within about two kilometers of one another. However, black spruce can grow on a variety of sites with characteristics that can change even over short distances, including on this study. Site and environmental factors that were not measured, such as peat depth, water table level, soil moisture, or pH may help explain the additional variability that is not captured in the model [19,37]. This could explain the success of adding random effects at any level (i.e., treatment, compartment, and/or plot).
The NLME model outperformed the refit Deo & Froese model across treatments, and the model performance varied slightly by treatment. The higher RMSE in the shelterwood and selection treatments likely points to the difficulty in predicting the diameter across treatments that experience patchy residual areas. The shelterwood treatment aimed to leave healthy and vigorous trees, while the group and individual tree selection methods aimed to result in uneven-aged stands. The subjective nature of cutting rather than a systematic treatment, as in the clearcut treatments, likely resulted in the patchy distribution of residual trees, and variable resources.
This modeling approach fits a single growth model to all diameters in the dataset (9.1-32.5 cm), which may be more appropriate than fitting two models for large and small trees. The work by Deo & Froese [25] refit the large tree diameter growth equations (trees greater than 12.7 cm) used in the Lake States Variant of FVS. The FVS framework models diameter growth differently for large and small trees (those below 12.7 cm DBH, and those equal to or greater than 12.7 cm). Fitting a single growth model for black spruce and other slow-growing species that attain smaller diameters may be advantageous over the current approach in FVS. The diameter growth curve ( Figure 3) under average stand conditions shows that diameter growth in peatland black spruce peaks around 12.7 cm. Since black spruce grows under a variety of site conditions, and maximum sizes vary largely, depending on those conditions, having separate equations for small and large trees with a transition near the peak of diameter growth for a species may not be ideal. Additionally, FVS diameter increment models for large and small trees were originally developed for the Intermountain West, where average and maximum diameters for tree species are generally greater than that of Lake States species. Fitting separate equations with a breakpoint diameter of 12.7 cm in those geographic regions may be more appropriate than in the Lake States, particularly for species with smaller maximum diameters, like black spruce. study. Site and environmental factors that were not measured, such as peat depth, water table level, soil moisture, or pH may help explain the additional variability that is not captured in the model [19,38]. This could explain the success of adding random effects at any level (i.e., treatment, compartment, and/or plot). The NLME model outperformed the refit Deo & Froese model across treatments, and the model performance varied slightly by treatment. The higher RMSE in the shelterwood and selection treatments likely points to the difficulty in predicting the diameter across treatments that experience patchy residual areas. The shelterwood treatment aimed to leave healthy and vigorous trees, while the group and individual tree selection methods aimed to result in uneven-aged stands. The subjective nature of cutting rather than a systematic treatment, as in the clearcut treatments, likely resulted in the patchy distribution of residual trees, and variable resources.
This modeling approach fits a single growth model to all diameters in the dataset (9.1-32.5 cm), which may be more appropriate than fitting two models for large and small trees. The work by Deo & Froese [25] refit the large tree diameter growth equations (trees greater than 12.7 cm) used in the Lake States Variant of FVS. The FVS framework models diameter growth differently for large and small trees (those below 12.7 cm DBH, and those equal to or greater than 12.7 cm). Fitting a single growth model for black spruce and other slow-growing species that attain smaller diameters may be advantageous over the current approach in FVS. The diameter growth curve ( Figure 3) under average stand conditions shows that diameter growth in peatland black spruce peaks around 12.7 cm. Since black spruce grows under a variety of site conditions, and maximum sizes vary largely, depending on those conditions, having separate equations for small and large trees with a transition near the peak of diameter growth for a species may not be ideal. Additionally, FVS diameter increment models for large and small trees were originally developed for the Intermountain West, where average and maximum diameters for tree species are generally greater than that of Lake States species. Fitting separate equations with a breakpoint diameter of 12.7 cm in those geographic regions may be more appropriate than in the Lake States, particularly for species with smaller maximum diameters, like black spruce. It is important that diameter growth models used by managers are applicable to the conditions where they are applied. This is of particular important in lowland black spruce in the Lake States, It is important that diameter growth models used by managers are applicable to the conditions where they are applied. This is of particular important in lowland black spruce in the Lake States, where the cover type tends to grow in pure stands, and growing conditions vary widely depending on its location and proximity to its range edges [1]. Understanding how silvicultural systems and standand tree-level attributes affect diameter growth and its modeling in varying geographical locations and physiological conditions is critical for managers to be able to accurately model future stand conditions and volumes.
This model fitting resulted in the same variables predicting the diameter growth, as was found with the same model for other species by Hann [22] and Weiskittel et al. [31], with the notable substitution of a crown position indicator variable for crown ratios. This substitution seems adequate, given that the crown ratio was not measured in the original study. While crown ratio as an absolute measure would be a more accurate portrayal of a tree's crown than the subjective measure of crown class, crown class may provide an alternative measurement to the crown ratio. This is of particular relevance in black spruce peatlands, where visibility in the canopy can be low, due to dense crowns and high stand densities.
The signs and values for both the GNLS and NLME models matched biological expectations for diameter growth. Positive coefficients occurred with log(DBH + 1), log(SI − 1.37), and the crown position dummy variable, indicating diameter growth, was positively correlated with these variables. Diameter growth was negatively correlated with BA 1/2 , DBH 2 , and BAL 2 log(DBH+5) . Notably, DBH growth is strongly positively correlated with log(DBH + 1), and negatively correlated with DBH 2 , showing that black spruce individuals with particularly large initial diameters will grow more slowly than trees with smaller diameters. This is evidenced when examining diameter growth in average stand conditions using the NLME parameters ( Figure 3).
One important factor in black spruce growth that was not considered in this study is eastern spruce dwarf mistletoe. While it was not present on the study site, dwarf mistletoe is a prevalent health concern in black spruce forests in Minnesota, and more broadly, in the Lake States, and it considerably slows growth and reduces stand volumes [15,16]. Future studies in the diameter growth of black spruce could focus on the effects of mistletoe.

Conclusions
Individual-tree non-linear fixed effects (GNLS) and mixed effects (NLME) 10-year diameter growth models were adapted and fit to lowland black spruce across six silvicultural treatments. Adapting the model by Hann [22] was successful in predicting 10-year diameter growth in lowland black spruce across the breadth of treatments, resulting in similar fit statistics to other diameter growth models (R 2 and RMSE). Substituting crown class for crown ratio as a predictor variable due to variable availability proved to be an effective substitution. Adding a random effect of plot nested within compartment nested within treatment in the NLME model improved its fit over the GNLS model. The GNLS and NLME modeling approach more accurately predicted diameter growth compared to the conventional approach used in FVS of using separate small and large diameter models. We recommend this approach as an improved alternative in the lowland black spruce cover type to the conventional procedures in FVS, where separate diameter growth models are used for small-and large-diameter trees.