^{1}

^{*}

^{2}

^{1}

Natural regeneration is crucial for silvicultural approaches based on the continuous presence of a forest cover, or Continuous Cover Forestry (CCF). Sitka spruce (

Continuous cover forestry (CCF), a range of silvicultural approaches involving uninterrupted maintenance of forest cover and avoidance of clearcutting (

Sitka spruce (

Seed availability is undoubtedly the first crucial factor: Sitka spruce seeds, like those of most temperate forest tree species, have a low survival rate in the forest soil and do not produce a viable seed bank. Sitka spruce in the UK starts to have a good seed crop at 25-35 years, after which the seed production increases with age and can reach high levels already at 35-40 years, depending also on the stand density (

Seed germination is highly dependent on the seedbed characteristics.

Stand structure can affect regeneration through different mechanisms. A certain level of overstorey cover was found to be beneficial for Sitka spruce regeneration both in its natural range (

Thinning interventions have been shown to have a positive effect by creating a favourable light environment. Studies carried out both in the UK and in North America generally found more Sitka spruce seedlings in the stands with lower densities, which were either more recently or more heavily thinned (

The UK has been defined “data-poor” regarding natural regeneration (

In the absence of studies following the development of regeneration over time, the dataset generated by

Age of the plantation in years (from now on simply Age), Soil Nutrient Regime (SNR), time after last thinning, and Deer Impact Index (DII) were recorded at stand level. We calculated from the original inventory the plot level values for basal area (BA), stems per ha (SPH), and the maximum DBH (maxDBH). From those values we calculated at plot level the quadratic mean diameter (QMD, the diameter of a tree considered as having the average basal area); and the Global Site Factor (GSF), an indication of the canopy light transmittance, using the relationship established from

As an indication of seed availability, we investigated the use of Age and two possible alternatives.

The SNR was estimated by the original field surveyor from analysis of the ground vegetation following the Ecological Site Classification criteria (

We considered the plot-level stand density measures of BA, SPH and maxDBH as a negative proxy for the light regime under the forest cover (higher stand density, lower light level) and so expected to be negatively correlated with regeneration presence. On the other hand, GSF is a direct indication the light regime under the forest cover, expected to be positively correlated with regeneration presence. The time since the last thinning was estimated for each stand using both historical records and evidence on the ground; the expected effect was a negative correlation between the time since the intervention and the likelihood of regeneration. We divided the stands in the present study into three different Thinning Classes (TC) as in

The Deer Impact Index (DII) was visually estimated as low (no browsing observed), moderate (browsing damage on up to 25% of the regeneration) and high (browsing damage on more than 25% of the regeneration). Because of the unbalanced distribution (

Additionally, we retrieved stand-level geographical variables from topographic maps, namely northing, easting, elevation and aspect, and stand-level climatic variables from the Forestry Commission’s decision support system ESC-DSS (

The density of Sitka spruce seedlings per plot was very different between the Thinning Classes (

For independent validation, we assessed in 2016 four Sitka spruce-dominated stands in Clocaenog forest, Denbighshire, Wales (53° 04′ N, 3° 25′ W; altitude: 390-430 m a.s.l.), and four in Kielder forest, Northumberland, England (55° 10′ N, 2° 29′ W; altitude: 200-250 m a.s.l.). Both forests were originally artificial plantations that have been managed in recent years according to different CCF principles, using silvicultural systems ranging from irregular shelterwood to group selection. All stands belonged to Thinning Class 2, but most of them were thinned more frequently or with higher intensity in the past than stands in the calibration dataset. The situation in all stands was generally a lower tree density than under the traditional management (as defined by

We carried out all the analyses using R statistical software (

We assessed the accuracy of the best model with a cross-validation technique (

In the first method, we defined a cut-off likelihood value using the Receiver Operator Characteristics (ROC) curve method with the package pROC (

We used two approaches. First, we investigated GLMMs using the same random and fixed effects as described above, using the sub-dataset for plots with presence of regeneration (n = 138), and a Gamma distribution with log-link to approximate the seedling distribution. No preliminary model based on all plots with presence of regeneration (n = 138) could converge (results not shown). The importance of the Thinning Class was evident from the sharp difference in seedling distribution amongst the classes, so we decided to calibrate separate models for TC1 and TC2 & TC3 (pooled together due to the lower number of observations). For those two subsets of data, we prepared GLMMs using the same random and fixed effects as described above (excluding Thinning Class). Then we removed non-significant parameters using a step-wise approach aimed at reducing the AIC to select the best model. We evaluated its accuracy through comparing predicted and observed values at plot and stand level. For the second approach, we simulated the seedling density by generating random numbers that approximated the observed density distribution for each Thinning Class (

We calculated the likelihood of regeneration presence in the independent validation plots using the best model above selected (calibrated on the full dataset). Then, we used the same two methods as before to assign the presence of regeneration. First, we considered the same cut-off likelihood value previously determined with the ROC method, assigning the status of presence of regeneration to all plots above that threshold. Second, we used the stochastic method to randomly determine the presence or absence of regeneration. For both methods, we built contingency tables at plot level and examined the performance at stand level by comparing the total numbers of simulated and observed plots with regeneration, with the same procedures described above for the cross-validation. Then, we used both seedling density modelling methods prepared with the calibration dataset to simulate the density in the plots of the independent datasets with observed presence of regeneration. The simulated seedling density was compared with the observed values. A summary of the workflow for the statistical analysis is shown in

The model structure after the step-wise AIC reduction process is shown in Model 1 (

where

The model did not converge when the forest-level random effect was included, so we maintained only the stand-level effect. The effect of bare ground was not significant, and it had a weak negative relationship with regeneration, contrary to the hypothesis. Only the quadratic term for BA remained in the best model structure amongst the stand density indicators. Note that values of BA were divided by 100 since they were on a different scale from the other variables.

^{2} ha^{-1}. For TC2, only in old stands (more than 70 years old) was the probability of regeneration above 0.5, while for TC3 the likelihood was always low. The effect of mosses on regeneration likelihood was more linear.

After the cross-validation analysis, with the ROC method, the cut-off likelihood value for the regeneration presence probability was 0.3.

When the results were aggregated at stand level for the ROC method, 21 stands out of 34 had a difference between total observed and predicted regeneration plots equal to or lower than 20% (11 with no difference), while five had a difference equal to or larger than 50% (worse than chance). For the stochastic method, very similar results were obtained: 22 stands out of 34 had a difference between total observed and predicted regeneration plots equal to or lower than 20% (10 with no difference), while five had a difference equal to or larger than 50%.

The worst simulated stands were almost the same stands in both methods. The field notes provided additional insights about them, showing that they were generally the ones subjected to heterogeneous thinning interventions within the same stand, suggesting that the TC class was inaccurate. In stand with fewer simulated regenerating plots than observed, it was also observed that windblow events had opened gaps comparable to a thinning, or that there was precocious cone production in young stands. In stands with more simulated regenerating plots than observed, it was noted that in stands favourable for regeneration according to all the model variables, the limiting factors were likely to be: (i) competing ground vegetation; (ii) presence of deer browsing; and (iii) lack of cone production. In the two worst over-simulated stands with both methods, everything seemed suitable for regeneration based on the field notes, and its total absence was inexplicable for the surveyor too.

In the GLMMs calibrated for TC1 and TC2 & TC3, only the effect of BA was significant, but with a positive relationship with seedling density in the former class (TC1) and a negative relationship for the latter group (TC2 & TC3). However, both models showed a very poor fit between the simulated and observed density values and they were discarded (results not shown).

The Weibull distributions fitted to seedling density distribution in each TC are described by the parameters in ^{-1}, but with extremes of -177,500 and 59,000 seedlings ha^{-1}.

We used

Regeneration density was then estimated in the plots with observed regeneration presence (n=62). Only the Weibull distribution approach was used, with the function previously calibrated for TC2. The GLMM approach was already deemed too inaccurate. After averaging the results, there was no good correspondence between the simulated and observed values, and no stands had a simulated density ± 20% of the observed value. On average, the difference between simulated and observed values was -34,570 seedlings ha^{-1}, with extremes of -155,800 and 4,500 seedlings ha^{-1}.

The model predicting regeneration presence was based on the established knowledge of the biological and ecological characteristics of Sitka spruce. The effect of time since the last thinning showed the strongest significance in the model, and the largest coefficient after standardization. Consistently with

The age of the plantation emerged as the second most important factor. Such a positive effect in the artificial plantations of the present study can be explained by the larger seed production of older trees, and possibly also by the higher number of gaps that can naturally occur in a mature canopy past the self-thinning stage. We tested the use of maximum DBH (at plot level) and quadratic mean diameter (at stand level) as possible alternatives to age, but in this research they were both negatively correlated with regeneration presence. For maximum DBH, it is likely that large trees present in the small study plots (5.6 m radius) were shading the ground and dispersing their seed outside the plots.

Mosses showed a positive effect on regeneration consistent with previous findings. A thin layer of mosses cover is favourable for germination due to their water retention capacity, but heavy mosses can prevent roots from reaching the mineral soil (

Increasing competition from the overstorey, expressed here as the quadratic term of basal area, influenced the regeneration negatively. However, for Thinning Class 1, at low overstorey levels the effect was relatively low and almost flat, likely confirming the benefit of a certain amount of shading. The same levels of basal area can be obtained with different numbers of trees, resulting in different canopy structures and thus light availability on the ground. When the number of trees is lower for a given basal area, there are likely to be more gaps between crowns and significantly more light at ground level (

None of the topographic and climatic variables tested showed significance. The climatic data were interpolations for 10 km grid squares of average climatic data collected during 1960-1990. They had already been found in another study to lack the precision needed for stand analysis (

The cross-validation process with the use of the Response Operator Characteristics curve showed satisfactory statistical results at plot level: 73% of plots were correctly simulated, with similar values for specificity and sensitivity. The stochastic method, such as is employed by various models, showed worse results: 66% of total plots were correctly simulated, with a larger difference between sensitivity and specificity. However, when aggregating the results at stand level and considering the difference between the total simulated and total observed plots with regeneration, the results were similar between methods: around two-thirds of the stands showed an acceptable error (simulated values within ± 20% of observed values). In a non-spatial forest growth simulator (

The results of the independent validation with the Response Operator Characteristics curve method were not satisfactory since the model predicted regeneration in almost all plots, even if the total accuracy was 82%. Using the stochastic method, the total accuracy was worst (64%), although there was a slightly better balance between sensitivity and specificity. It is evident that the independent dataset is describing a situation largely different from the calibration dataset, noting the differences both in the stand variables (

The models tested here for regeneration density did not give results of acceptable accuracy. Generating random numbers from Weibull distributions was, in the present study, the only option found and still produced inadequate results both during the auto-validation and independent validation. Nonetheless, even if the models were deemed too inaccurate, it is interesting to note that the effect of basal area was significant and positive in the seedling density model based only on plots belonging to Thinning Class 1, suggesting a possible mother-tree positive effect. In the model for Thinning Class 2 & 3, basal area had a negative effect, maybe because the already-lower light availability is aggravated by bigger tree size and the overstorey competition effect becomes predominant. Similar results were observed by

A very important limitation of both models was the lack of data on the regeneration size or age. Both the regeneration presence and density model did not consider the possibility of other tree species germinating and competing with Sitka spruce, likely another crucial limitation of the use of these models in mixed forest stands resulting from continuous cover forestry practices. Presence of deer browsing, although not statistically significant in this analysis, was found in the field notes as a possible cause of limiting factor for regeneration in some sites where all the model variables were at a beneficial level for regeneration.

Concluding, the tools here described can be used to simulate regeneration presence in traditional Sitka spruce plantations in the UK. Then, the growth of the regeneration can be predicted with the light-growth models presented by

We wish to thank: Dr. Gary Kerr, Dr. Victoria Stokes and Barnaby Wilder (Forest Research) for collecting and providing the original dataset used in this research; Dr. Catia Arcangeli (Forest Research) and Dr. Caterine Cahalan (Bangor University) for their comments on the study; and the Forestry Commission and the Scottish Forestry Trust for funding the study.

Frequency of seedlings per hectare in different Thinning Classes (TC1, TC2 and TC3), only plots with presence of regeneration, for the calibration dataset (left) and the independent validation dataset (right).

A summary of the workflow for the statistical analysis.

Regeneration presence likelihood (_{regen}) as a function of the model variables. In each graph, the likelihood was estimated with only one variable varying across all its range (plotted on the

Coefficient values after standardization of the model variables (BA = Basal area, TC = Thinning class). The dot corresponds to the mean values, the blue line to the 90% confidence interval.

Receiver Operator Characteristics (ROC) curve for the cross-validation method. The dot represents the point with the highest sum of the specificity and sensitivity values (presented between parentheses) and shows the corresponding cut-off likelihood value.

Some of the crucial factors influencing Sitka spruce natural regeneration, the general conclusions drawn in the literature about them, and the evidence quality of such conclusions. Adapted from

Factor | Conclusions | Evidence Quality |
---|---|---|

Seed availability | Mast years very important, in British Sitka spruce stands happening every 4-5 years | Good-Moderate |

Germination conditions | Favourable seedbed conditions: moist soils with needle litter or light moss cover | Moderate-Poor |

Vegetation competition | Avoid fertile sites or competition from ericaceous vegetation | Moderate |

Understorey microclimate | Retain some canopy cover to limit frost damage but provide adequate light | Moderate |

Light requirements for growth | At least 20% of full light, plus an overstorey with basal area of 30 m^{2 }ha^{-1} and reduced tree density |
Good |

Browsing pressure | Keep deer population below 5 animals per 100 ha | Moderate |

Details of the calibration dataset. Values at stand (Age, Quadratic Mean Diameter, Soil nutrient regime, Time after last thinning and Deer Impact Index) and plot level (the remaining parameters).

Variable | Min. | 1^{st} Qu. |
Mean | 3^{rd} Qu. |
Max. |
---|---|---|---|---|---|

Age (years) | 32 | 39 | 54.5 | 64 | 85 |

Basal area (m^{2 }ha^{-1}) |
1.6 | 43.6 | 58 | 70 | 196 |

Stems per hectare (n ha^{-1}) |
0 | 400 | 700 | 900 | 2.200 |

Quadratic mean diameter (cm) | 0 | 27 | 36.3 | 43.1 | 83 |

Maximum diameter breast height (cm) | 0 | 36 | 45.3 | 52 | 90 |

Global Site Factor | 0.02 | 0.16 | 0.21 | 0.26 | 0.55 |

Bare ground (%) | 0 | 0 | 1.2 | 0 | 85 |

Mosses (%) | 0 | 5 | 41.6 | 80 | 95 |

Seedling density (ha^{-1}) |
0 | 0 | 20.780 | 10.000 | 450.000 |

Soil Nutrient Regime (no. plots) | Very Rich (10) | Rich(10) | Medium(130) | Poor(180) | Very poor(10) |

Time after last thinning (no. plots) | Class 11-5 yrs(170) | - | Class 26-10 yrs(90) | - | Class 310+ yrs (80) |

Deer Impact Index(no. plots) | Low(10) | - | Moderate(290) | - | High(40) |

Details of the validation dataset. Values at stand (Age and Time after last thinning) and plot level (the remaining parameters).

Variable | Min. | 1^{st} Qu. |
Mean | 3^{rd} Qu. |
Max. |
---|---|---|---|---|---|

Age (years) | 60 | 65 | 69 | 77 | 80 |

Basal area (m^{2 }ha^{-1}) |
7.6 | 28.8 | 41.4 | 53.8 | 107.2 |

Stems per hectare (n ha^{-1}) |
50 | 200 | 284 | 400 | 1.100 |

Maximum diameter at breast height (cm) | 35 | 44 | 50.6 | 55.8 | 85 |

Global site factor | 0.08 | 0.22 | 0.28 | 0.34 | 0.49 |

Bare ground (%) | 0 | 0 | 0.1 | 0 | 4 |

Mosses (%) | 0 | 70.6 | 82.3 | 99.7 | 100 |

Seedling density (ha^{-1}) |
0 | 2.500 | 48.460 | 52.500 | 417.500 |

Time after last thinning(no. plots) | Class 11-5 yrs(0) | - | Class 26-10 yrs(78) | - | Class 310+ yrs(0) |

Details of coefficients of Model 1, estimating the likelihood of Sitka spruce regeneration presence at plot-level.

Effects | Variables/Stats | Value | St. Err | p-value |
---|---|---|---|---|

Fixed | (Intercept) | -2.693 | 1.703 | 0.114 |

Thinning class | -1.864 | 0.556 | 0.001 | |

Age | 0.087 | 0.03 | 0.003 | |

Mosses | 0.02 | 0.007 | 0.005 | |

(Basal area/100)^{2} |
-1.569 | 0.886 | 0.076 | |

Random Effects (intercept) | Variance | 3.726 | - | - |

Standard Deviation | 1.93 | - | - |

Contingency tables for both the cross-validation and the independent validation results, using both the Response Operator Curve (ROC) method and stochastic method. “Yes” indicates the presence of regeneration, “No” the absence.

Validation | Method | Observed | Predicted | PartialAccuracy | Overallaccuracy | ||
---|---|---|---|---|---|---|---|

Yes | No | Total | |||||

Cross-validation | ROC | Yes | 105 | 33 | 138 | 0.76 | 0.73 |

No | 58 | 144 | 202 | 0.71 | |||

Total | 163 | 177 | 340 | - | |||

Stochastic | Yes | 76 | 62 | 138 | 0.55 | 0.66 | |

No | 52 | 150 | 202 | 0.74 | |||

Total | 128 | 212 | 340 | - | |||

IndependentValidation | ROC | Yes | 62 | 0 | 62 | 1 | 0.82 |

No | 14 | 2 | 16 | 0.12 | |||

Total | 76 | 2 | 78 | - | |||

Stochastic | Yes | 44 | 18 | 62 | 0.71 | 0.64 | |

No | 10 | 6 | 16 | 0.39 | |||

Total | 54 | 24 | 78 | - |

Parameters for the Weibull distributions fitted to seedling density per ha (× 1000).

ThinningClass | Shape | Rate |
---|---|---|

TC1 | 0.696 | 52.555 |

TC2 | 0.871 | 14.134 |

TC3 | 1.834 | 4.651 |