Modeling recent bark stripping by red deer (Cervus elaphus) in South Belgium coniferous stands

Abstract

• Context

Over the past few decades, the impact of large herbivorous ungulates on forest vegetation has been clearly highlighted. Among those impacts, bark stripping of coniferous trees is one of the most damaging. Bark stripping leads to rot development, inducing serious loss of timber value.

• Aims

The present study aimed firstly at evidencing the factors explaining the variations observed in fresh bark peeling rate for spruce and Douglas-fir in southern Belgium and secondly at identifying the key factors to consider when setting up a deer management plan.

• Method

Fresh bark peeling rate was recorded with a systematic sampling survey from 2004 to 2007. The covered territory was then divided into 63 distinct hunting zones of area ranging from 1,000 to 25,000 ha. About 5,000 plots were monitored annually. Each zone was characterized with a large number of explanatory variables. The explanatory variables were integrated firstly into fixed linear models using a stepwise procedure, and then into a mixed model.

• Results

The significant variables included in the model (R ² = 44 %) are (by decreasing order of importance) red deer densities, proportion of coniferous stands and agricultural areas, snow cover, distance to urban habitats, and species diversity in the understory.

• Conclusion

The models revealed the impacts of several factors on bark peeling: deer density, deer-carrying capacity of the territory, landscape structure, and severity of winter conditions. The adjusted model allowed subtracting the impact of winter conditions in order to produce a relevant indicator for hunting management. In addition, the model was used to assess the sensitivity of a forested area to bark peeling based on its environmental characteristics.

1 Introduction

The impact of large herbivorous ungulates on forest vegetation has been clearly evidenced (Reimoser et al. 1999). Over the past few decades, the increasing level of damage caused by these animal populations has become a source of concern for forest managers in Europe (Jerina et al. 2008; Kiffner et al. 2008; Klopcic et al. 2010; Vasaitis et al. 2012; Verheyden et al. 2006) as well as in other parts of the world (Takatsuki 2009; Hall and Gill 2005; Dolman and Wäber 2008). The causes of this damage can be classified into two types: (1) grazing, affecting regeneration mostly in its installation or growing phases; and (2) bark stripping, plaguing already established stands (Kiffner et al. 2008). Coniferous trees are more affected by bark stripping than broad-leaved trees (Vospernik 2006).

The consequences of bark stripping are closely linked to the tree’s specific sensitivity to pathogens once the wood has been deprived of its protective bark (Kiffner et al. 2008). Norway spruce represents one of the most sensitive species to pathogens (Jerina et al. 2008; Verheyden et al. 2006; Vospernik 2006). Consequently, the proportion of stripped stems that develop rot can reach 60–100 % (Girompaire and Ballon 1992; Vasiliauskas 2001). Once rotting has set in, the lower part of the trunk, the most valuable part, cannot be used for the most profitable uses (Heikkilä and Härkönen 1996; Simon and Kolá 2001).

This problem of tree damage caused by bark stripping calls for the setting up of systems that can monitor regularly the pressure exerted by herbivores, and more specifically by red deer (Cervus elaphus), on the forest ecosystem. Tools also need to be developed that are capable of analyzing the environmental factors interacting with the damage that occurs (Bugmann and Weisberg 2003). Indeed, although the increasing deer population may explain the rise in the level of bark stripping damage, the evolving characteristics of the ungulate habitat cannot be overlooked as a contributory factor (Reimoser et al. 2003).

To date, quantitative studies aiming at analyzing the variations in bark stripping activity have been rare or limited to rather small study areas (Honda et al. 2008; Jerina et al. 2008; Kiffner et al. 2008; Vospernik 2006). Models in those studies are functions of deer density (Jerina et al. 2008; Kiffner et al. 2008); stand characteristics such as age, density, or composition (Honda et al. 2008; Kiffner et al. 2008; Vospernik 2006); the surrounding environment, characterized, for example, by proximity to forest edges (Honda et al. 2008; Kiffner et al. 2008); topographic position (Jerina et al. 2008; Kiffner et al. 2008; Vospernik 2006); and elevation or the mean annual snow depth (Honda et al. 2008; Kiffner et al. 2008; Vospernik 2006).

The present study aims at evidencing the factors that might explain the variations observed in the rate of fresh bark peeling damage found in spruce and Douglas-fir stands located in Wallonia (southern Belgium). Whereas previous studies have established models of bark stripping at tree or plot level, we study instead the damage caused by deer on a broader scale. Our aims justify this approach since the scale of the inventory zones studied corresponds to the scale of hunting management planning for the area (for example, culling plans). In addition, key contributory factors are identified for the total fresh bark peeling damage and also separately for fresh damage caused during the winter and during the summer. Finally, we illustrate how those findings can be utilized in drawing up a hunting management plan.

2 Materials and methods

2.1 Study area

The bark stripping damage survey took place in the Belgian Ardenne (southeast Belgium). This region features a high forest cover (52 %), which corresponds to 62 % of the total Walloon forest (Rondeux and Lecomte 2002). The region is largely covered with coniferous stands (63 %), the majority of which are even-aged spruce stands (78 %). The existing Walloon forest inventory reveals that a large proportion (29 %) of coniferous stands exhibit bark peeling damage.

Inside the Ardenne region, we defined 63 survey zones to assess the level of fresh bark stripping damage (Gheysen et al. 2011). These survey zones were drawn up in collaboration with the Walloon Forest Services (DNF) so that they coincided with relatively homogenous hunting zones; the zones were based on knowledge of current red deer distribution (Fig. 1). The areas of the survey zones ranged from 1,000 to 25,000 ha, with an average of 9,000 ha. Altitude varied between 228 and 643 m above sea level (429 m on average) and forest cover reached 17–94 % (60 % on average). Coniferous stands represented 17–91 % of the forested area (63 % on average). The estimated deer density ranged from 0.2 to 8 individuals per square kilometer, with an average of 2.8 individuals per square kilometer.

Fig. 1

Delineated inventory zone for fresh bark peeling damage in the Walloon Region (upper left corner) and forest areas defined in this study (relatively homogeneous hunting zones) (right hand side)

2.2 Bark stripping inventory

Observations of recent bark stripping damage were recorded annually. We used two-phase systematic sampling with a square sampling grid of 200 × 200 m. Firstly, we selected the grid intersections that corresponded to sensitive stands. The sensitive stands were spruce stands of between 8 and 30 years old, and Douglas fir stands of between 8 and 36 years old (Gheysen et al. 2011). Secondly, three circular sampling subunits with six trees were delineated following a north–south transect (Fig. 2). The inventories were carried out between mid-April and mid-May.

Fig. 2

Sampling design for inventorying fresh bark peeling damage in coniferous plantations in Wallonia since 2004. a First step: square sampling grid with studied plots. b Second step: coniferous stand sampling. c Example of plot with its three subplots

For the purpose of this study, we took into account four annual data sets recorded between 2004 and 2007. Those data sets encompassed between 4,286 and 5,214 plots, and the number of observed trees varied between 64,395 and 75,424. This allowed us to compute the estimates of total fresh bark peeling rate (BPR), summer fresh bark peeling rate (BPRS), and winter fresh bark peeling rate (BPRW) for each zone and year (Table 1). BPR was the sum of BPRW and BPRS. More detailed information regarding these three models can be found in Gheysen et al. (2011).

Table 1 Summary of results from the inventory of fresh bark peeling damage for the inventory campaigns of 2004, 2005, 2006, and 2007

2.3 Explanatory variables

Bark peeling intensity, expressed for each inventory zone by BPR calculations, was analyzed with reference to various environmental variables. These variables were selected based on not only the literature regarding bark stripping by deer (Gill 1992; Kiffner et al. 2008; Verheyden et al. 2006; Vospernik 2006) but also on the relevance and precision of the available variables for the study area. Finally, 28 explanatory variables were selected from various sources (Table 2): Mapping data (such as topographic situation and land use) were drawn from processing topographic maps at 1:50,000 (Institut Géographique National de Belgique, 1994–2001), the digital land model at 1:50,000 (Institut Géographique National de Belgique, 1994), and the land use map (Service public de Wallonie—Direction générale opérationnelle: Aménagement du territoire, Logement, Patrimoine et Energie, 1993). The Navstreet databank (http://www.navmart.com/) was then used to determine the proximity of the road infrastructure. These various layers of data were processed using the ArcGIS® software, and the Patch Analyst software was particularly useful for describing the forest patchwork (Wheatley 2010). Weather data were provided by the Belgian Institute of Meteorology (RMI). The level of human disturbance was also evaluated with an ordinal level of visitation divided into four classes (Colson et al. 2010).

Table 2 Number of explanatory variables and sources of information

The selected explanatory variables were thus predominantly environmental ones (Table 3). They referred mostly to spatial variability and were considered as constant for the duration of the bark peeling observations. The only variable involving temporal variation concerned the severity of the winter, described as the number of days when the ground snow cover equalled or exceeded 2 cm at 8:00 am. This variable was defined by the RMI in order to take into account only the days when snow completely covered the whole area, and not just certain parts. The average number of snow days varied between 7 (in 2007) and 60 (in 2006).

Table 3 List of the main explanatory variables whose impact on the severity of bark peeling rate was tested in the present study

Deer density had previously been estimated by the DNF using kilometric indexes, spotlight counts, and counts near supplemental feedings. Deer density was assumed to be constant during the present study because of the low level of precision of those estimations. The relevance of this variable was mainly to be found in its relative value among the survey zones. Deer density was also evaluated on the basis of hunting bag statistics (number of shot animals per thousand hectares). Indeed, every animal that is shot or found dead has to be reported to the DNF (Fichant 2003).

3 Statistical method

The level of BPR could be explained by deer density and the environmental carrying capacity for the species. We modeled these phenomena using two-part models (Zuur et al. 2009): on the one hand, the absence–presence of bark peeling and, on the other hand, the bark peeling intensity. Firstly, we modeled the absence or presence of damage using a logistic regression based on deer density in the 63 zones. The chi-square test was then used to compare the model obtained with the null model. Secondly, we used multivariate regression analysis to study BPR, BPRS, and BPRW. These three models were adjusted with only the zones where bark peeling damage had been observed. These models were thus used to analyze observations of fresh bark peeling damage for 51 zones during the 2004–2007 inventory campaigns (n = 204).

Next, as the dependent variables each corresponded to a proportion of the zones, we applied an angular transformation (Eq. 1) in order to guarantee appropriate application conditions in the case of linear regressions (Dagnelie 2006):

$$ Y\prime = 2\arcsin \sqrt {Y} $$

(1)

where Y′ represents the transformed variable and Y the original one. Moreover, the Box and Tidwell transformation (Palm et al. 2011) was used to test whether it was necessary to transform the explanatory variables.

After these preliminary treatments, a forward and backward stepwise regression procedure was applied in order to select the most relevant explanatory variables and to design the corresponding model (Kowalski and Montgomery 2011). Although logistic regression is used more frequently to model the occurrence of deer damage at tree level (Honda et al. 2008; Jerina et al. 2008; Kiffner et al. 2008; Vospernik 2006), we chose to use linear regression as we were modeling deer damage at the inventory zone scale.

As the data recording had been repeated over time, our observations were not independent. In order to obtain unbiased estimations, mixed models were adjusted to take into account the random effect of the inventory zones (Bates and Maechler 2009).

Finally, the contribution of each variable was evaluated by calculating the type II sum of squares (hierarchical or partially sequential sum of squares). These sums of squares measured the decreased residual sum of squares due to an added variable once all the other variables have been introduced into the model (Langsrud 2003). Here, this indicator made it possible to classify the variables according to their importance in the model.

All of these analyses were performed with the R software, version 2.9.2 (R Development Core Team 2008).

4 Results

Table 3 below gives an overview of the tested explanatory variables, the value of the correlation coefficients between these variables, and the bark peeling rate (BPR), as well as other abbreviations used thereafter.

The results of the logistic regression (Fig. 3) suggested that deer density had a highly significant effect (χ ² = 37.8, p < 0.001, n = 62) on the presence of BPR.

Fig. 3

Logistic regression of the absence–presence of damage per forest area in terms of deer density. It should be noted that the probability of finding fresh bark peeling in forest areas with fewer than 10 deer/10,000 ha is below 20 %

The stepwise procedure for BPR gave the following model:

$$ \mathrm{ASinBPR}={{{-0.204-0.00315\,\mathrm{AgriProp}+0.00260\,\mathrm{ConifProp}+0.0606\,Sqrt\left( {DeerDensity} \right)+0.0000283\,\mathrm{SnowCove}{{\mathrm{r}}^2}+0.280}} \left/ {{SaplingDiversity-0.0000500\,\mathrm{UrbanDist}}} \right.} $$

The model adjusted determination coefficient ($R_{\mathrm{adj}}^2$) was 44 %. This value represented the variability part of the transformed dependent variable (AsinBPR) explained by the model. Estimating this indicator in the case of the initial explanatory variable (BPR) was the result of Eq. 2 (Kowalski and Montgomery 2011).

$$ R_{\mathrm{adj}}^2=1-\left[ {\frac{{\sum {{{{\left( {{y_i}-y_i^{\prime }} \right)}}^2}} }}{{\sum {{{{\left( {{y_i}-{y_{\mathrm{moy}}}} \right)}}^2}} }}} \right]\left( {\frac{n-1 }{n-p-1 }} \right) $$

(2)

where y _i represents the ith answer observed; $y_i^{\prime }$, the ith answer estimated; y _moy, the average answer; n, the number of observations; and p, the number of variables in the model. As a result, $R_{\mathrm{adj}}^2$ equaled 42 % for the nontransformed BPR.

We then verified the absence of heteroscedasticity in the residuals by examining the graphs of residuals in relation to estimated values. The normality of the residuals was also tested using the Ryan and Joiner test (r _j  = 0.996, p > 0.1) (Kowalski and Montgomery 2011).

When we subsequently looked at the residuals of the 51 studied zones (Fig. 4), we noticed immediately that the average residuals of some zones were different from zero. This meant that we had to adjust a new model, which was a mixed model, by adding a random factor that grouped the data by zone. The estimated coefficients remained unchanged, but the procedure allowed us to correct our estimation of their variance (Table 4). The likelihood ratio test indicated that the addition of random effects was very significant (L = 51.570, p < 0.0001). The residual standard deviation was 0.0859, and the standard deviation of the zone grouping factor was 0.0838, confirming the importance of this random effect.

Fig. 4

Distribution of residuals per forest area. It should be noted that a large part of the variability remains unexplained at the forest area level

Table 4 Estimated parameters of the BPR model for the 51 study areas

The BPRW and BPRS models are introduced in Tables 5 and 6. Their adjusted determination coefficients ($R_{\mathrm{adj}}^2$) were, respectively, 46 and 32 %. Winter fresh damage accounted for approximately 81 % of the total fresh damage. Those models were close to the BPR model. Nevertheless, we observed some differences between the selected explanatory variables. The BPRS model understandably did not include the number of snow days, as this variable was not significant (df = 197, F = 0.455, P = 0.501). For the BPRW model, altitude was a better predictor of the bark peeling rate than the proportion of coniferous stands.

Table 5 Estimated parameters of the BPRW model for the 51 study areas

Table 6 Estimated parameters of the BPRS model for the 51 study areas

5 Discussion

The model enabled us to identify the significant impact on bark peeling rate of six of the 28 tested variables. As a result, the estimated deer density, the average number of snow days, and the local proportion of coniferous stands were all factors that contributed to an increase in the bark peeling rate. Conversely, the proportion of farming areas, the number of species found in the herbaceous stratum observed in neighboring stands, as well the average distance to urban areas led to a decrease in the bark peeling rate. The partial sequential sum of squares allowed us to classify each variable by degree of importance (Table 4). The most important variables were, in decreasing order of importance, DeerDensity, ConifProp, AgriProp, SnowCover, UrbanDist, and SaplingDiversity.

In order to better understand and interpret the model, we may group the explanatory variables into four categories based on their significance at the environmental level:

• Level of deer populations (DeerDensity)

• Richness of the environment and available food (SaplingDiversity)

• Spatial structure of the biotope (AgriProp, ConifProp, and UrbanDist)

• Seasonal food scarcity (SnowCover).

The type II sums of squares for these groups were, respectively, 39, 9, 41, and 11 %. These results show, first of all, that estimated deer population level is an undeniably important factor in accounting for the level of bark peeling damage. Nevertheless, other factors need to be taken into account: food availability, mentioned repeatedly in the literature (Gill 1992; Verheyden et al. 2006), also stands out in our model, where it is characterized by an abundance of conifers (ConifProp) and by the vegetal diversity found in the herbaceous stratum (SaplingDiversity). In the Ardenne region, coniferous stands—spruce stands mostly—are usually dense and rather poor in terms of food resources for ungulates (Lecomte et al. 2003). Moreover, coniferous stands constitute areas of cover that are appreciated by game. Indeed, several studies conducted in the Ardenne region (Fichant 2003; Licoppe 2006) have shown that these areas, and in particular young coniferous stands, are occupied preferentially by deer during the winter or in the day time. This refuge-providing role combined with a scarcity of food resources would thus partly explain the severity of bark peeling damage in coniferous stands.

Moreover, the presence of the UrbanDist variable confirms that human activity has a significant impact on the level of bark peeling (Petrak 1998). Based on the correlation coefficients obtained in Table 3, the disturbance would be linked mostly to the proximity of human activities and, to a lesser extent, to road infrastructure.

Contrary to our expectations, the “altitude” factor did not appear in the BPR model, although several studies have evidenced the importance of this factor (Debeljak et al. 2001; Vospernik 2006). However, in the study area, altitude was found to correlate strongly with other spatial variables ConifProp (r = 0.779, P < 0.001) and UrbanDist (r = 0.292, P < 0.001). The impact of altitude was thus concealed by the presence of other variables already included in the model. Our results showed that bark peeling damage consisted mainly of winter damage (80 %). Winter and summer damages depended on the level of the deer population, the richness of the environment, and the structure of the biotope. In contrast to summer damage, winter damage increased with the number of snow days and hence the rigors of winter and altitude. Summer damage appeared to be somehow linked to the proportion of coniferous stands. Those stands (mainly composed of Norway spruce) provided lower level of mineral salts, a factor that may contribute to a rise in the occurrence of bark peeling (Gill 1992).

Inventorying of fresh bark peeling damage constitutes an invaluable tool to assist forest managers in drawing up hunting management plans. Indeed, it allows them to follow the evolution of bark peeling damage over several years (Fichant 2003). The evolution of the bark peeling rate is, however, not only an indicator of the evolution of deer populations; it is also an indicator of forest-carrying capacity, mostly linked to the yearly availability of food, and thus to the severity of a particular winter (Gill 1992; Putman and Moore 1998), and also to the presence of mast. The model we have developed makes it possible to “correct” the observed bark peeling rate in order to moderate the effect of snow cover on fluctuations in bark peeling rate. This correction consists of estimating the total fresh bark peeling rate that would have been observed over a set number of snow days. Figure 5 illustrates this correction, where the number of snow days has been set at 30; this corresponds to the average for southern Belgium according to the Belgian Institute of Meteorology (RMI). The advantage of the “corrected” rate is that it more accurately reflects the evolution of deer populations, since fluctuations in those populations due to the thickness of snow cover have been suppressed or at least mitigated. For example, the corrected rate seems to indicate that red deer populations would have decreased in Elsenborn (zone ID = 15), whereas the numbers would have increased in Trois-Pont (zone ID = 22).

Fig. 5

A comparison of the evolution of the total fresh bark peeling rate observed and corrected in the 51 zones. In 2007, the winter was mild in contrast with the other 3 years, with 2006 showing the maximum severity. The corrected bark peeling rate is an indicator of the evolution of the deer population

Another advantage of our model lies in the possibility of attributing to each zone a level of sensitivity to bark peeling damage based on the zone’s intrinsic characteristics (proportion of conifers and plains, and average number of species involved in regeneration). Consequently, areas that are very sensitive to bark peeling can be identified, so that the bark peeling damage rate can be better interpreted. For example, the model can explain why some areas with a low red deer density show higher bark peeling rates than some other areas with a high red deer density. Hence, we were able to create a map showing the studied area as being subdivided according to the three levels of bark peeling sensitivity determined here (Fig. 6). This map gives forest managers a tool for gaining a better understanding of the balance between game density and the territory’s carrying capacity for deer. Indeed, in more sensitive areas, the carrying capacity will run out more quickly. Deer abundance in these areas therefore needs to be more strictly controlled.

Fig. 6

Classification of the inventory zone by degree of sensitivity to bark peeling damage: black highly sensitive, dark gray moderately sensitive, and light gray slightly sensitive

6 Conclusion

The results of this study, as well as the perspectives it has generated, show how useful inventorying bark peeling damage can be when setting up a deer hunting plan. When deer populations are estimated at over ten animals per thousand hectares, the present model indicates the impact on the level of bark peeling damage caused by deer density, by the zone’s carrying capacity for deer, by the landscape structure, and by the severity of winter conditions.

Moreover, this model makes it possible to quantify the variability of bark peeling due to weather changes, and in particular due to the severity of winter. Thus, the estimated impact of winter conditions can be subtracted from the estimated level of bark peeling damage, providing a relevant indicator for use in hunting management and, specifically, culling plans.

In addition, an index of sensitivity to bark peeling can be calculated for each zone, which is based on the zone’s intrinsic characteristics. With the help of this index, the bark peeling rate, computed as a result of the many data collection campaigns, can be interpreted more judiciously and with reference to a specific geographic area.