IntroductionTop

Windthrows are natural hazards and a significant body of literature has been published on this particular issue (Ulanova, 2000; Klaus et al., 2011; Nilsson et al., 2004; Schliemann & Bockheim, 2011; Boon, 2012). Nowadays, windthrows are included into the large concept of natural disturbances (Popa, 2008), along with wildfires, droughts, and insects’ attacks.

One important outcome and symptom of the climate change is the frequency and geographical distribution of endemic windthrows (EW), which are affecting larger and larger areas (Res et al., 2006; Senf & Seidl, 2017). EW were defined by Lanquaye-Opoku & Mitchell, (2005) as disturbances produced by peak winds, with return intervals shorter than 5 years. The amount of wood blown down may vary from tens of cubic meters to hundreds and this wide interval of variation is caused by the ‘domino effect’, quite difficult to take into consideration (Nolet, 2012; Nolet & Béland, 2017). Lohmander & Helles, (1987) summarized a series of theories on the causes of windthrows while (Mitchell, 1998) came up with a methodological framework based on the so-called windthrow triangle, which makes use of data concerning the stand condition, soil, and exposure.

Due to the large variety in geographic distribution and economic loss, studying natural disturbances under climate change is very provocative in statistical terms. Yousefpour et al., (2012) ranked the methodological approaches and found out that the most used models are geometric Brownian motion and probability distribution. The authors also produced some interesting statistics of the articles published on risk-related issues on forestry and concluded that the trend of articles published each year is exponential. Markov Chains and Bayesian models, we have used in our exercise, were among the least used approaches.

Long-term sustained yield has long been modeled as MC (Kouba, 2002) and Buongiorno & Zhou, (2011) demonstrated that even the Faustmann’s formula is a particular case of MC, assuming that all transitions from one phase to the next one are certain, meaning that all stands will reach the prescribed maturity age.

Holecy & Hanewinkel (2006) produced a forest insurance model based on probabilities given by Weibull functions plugged into a MC model while Strigul et al., (2012) developed a MC model to study the dynamic mosaic of forest stands; they concluded that age-class structures of the forest are skewed to the right, which is consistent with a MC stationary distribution.

The new paradigm of adaptive forest management stems from a mixture of extensive use of recent knowledge, risk management and a series of measures able to maintain the ecosystem resilience (Bolte et al., 2009). In line with this new paradigm, the forest management planning shall come up with a more flexible approach of forest sustainability and, in particular, a deeper understanding of the sustained yield principle.

An extensive literature review on assessing the natural hazard risks in forestry has been published by Hanewinkel et al. (2010) who came up with a flowchart about the risks pertaining to different geographical scales of the decision-making process. However, the first step of the procedure, that is the framework analysis, takes into consideration just the regional climate models and the general circulation model, which are relevant for extreme events only (Popa, 2008).

The sustained yield principle (SYP) has long been adopted by Romanian forest policy makers and, it has been implemented into the algorithm of calculating the AAC too. As anyone would expect, SYP was an important pillar of forest management under the command and control economy, and it made sense for many reasons. During the transition to market economy SYP has been justified by the social concern that freeing too much the forest planning will encourage the forest owners to use any legal loophole to make the most of their own forests (Bouriand, 2005). Obviously, many illegal loggings reported by the media and the Court of Accounts in 2013 (Anonymous, 2013) gave to politicians a convincing reason for keeping a tight control on the annual harvest, despite that, after 2008, the allowable cut adopted by the forest management planning had been rigorously checked by the public authority.

However, on a long run, pursuing only SYP could be ineffective and somehow dangerous for the forest economy and forest ecosystems, for two reasons at least: 1) SYP cannot go along with any ecological leeway simply because the planning methods based on SYP do not account for any uncertainty; and 2) pursuing SYP per se has created a positive feed-back loop because more salvage fellings in stands older than 60 years old have replaced mature stands initially included into the cutting budget and the harvesting plan.

Having salvage fellings instead of regular harvesting operations is even more convenient from a sheer economic standing point: according to the Forest Act, the regeneration fund, which is a sort of insurance deposit for having enough money to regenerate a given stand after final felling, is fed by 10% of the income provided by the main yield only. Such a ‘tax relief’ for salvage cuttings makes sense for catastrophic events, like the windthrow that took place in 1995, when more than 8 million m3 have fallen in one night (Popa, 2005). At that time, the broken trees were left un-hauled in remote areas for more than three years, due to extremely difficult harvesting conditions. Noteworthy, the Forest Act was adopted in 1996 when the consequences of that catastrophic event were apparent. On the contrary, when it comes to EW such a legal provision turns into a loophole, because more salvage fellings instead of regular harvest allow forest owners to avoid feeding this regeneration fund.

The aforementioned feedback loop has been imple­mented into an official regulation which states that every year the annual harvest is made of 80% of the current AAC, plus 20% of AAC of the previous year unless salvage products are to be harvested from stands older than 60 years.

If salvage fellings occur in stands older than 60 ye­­ars, due to whatever reason, their total volume is deducted from the previous year’s left AACR. If the salvage volume is higher than the previous year reserve, the difference will be deducted from the next year 20% AAC, and so forth. Eventually, having salvage cuttings each year, by the end of the planning period (every ten years) at least 20% of the harvestable volume will have been left uncut. This 20% AACR has neither economic nor technical support, being based on a rough assessment of the contribution of the catastrophic windthrows to the national cutting budget since 1995.

From the standpoint of mid and long-term planning, sooner or later the mature stands, postponed from harvesting, will eventually pile up into the last age class and the AAC stops being a means to normalizing the forest age structure. This particular situation shall also be addressed, in order to find out how the process of getting the mature stands older and older can be stopped.

Often, due to logging settings, salvage cuttings brought about by the windthrow are not applied only to the broken trees but also to some healthy trees that have to be removed in order to make broken trees accessible, which makes sense from the technological point of view (Waldron et al., 2013). Therefore, the scope of any perturbation is always larger than the forest area covered by the damaged trees.

When it comes to windthrows, the gap between case-studies based on different methodologies and effective headways towards a better forest planning is still large. Therefore, this study aims to producing an algorithm that utilizes the data referring to EW occurred in the past to improve the decisions regarding the annual cutting budget every next year. We developed a cascading Bayes model (Keith et al., 2003) to assess the EW likelihoods considering all possible permutations of three EW contributing factors (proximity of a prior EW, terrain conditions and slenderness). The EW likelihoods are further plugged into a MC model that shows the tendency of having a slant age class structure, regardless the normal age pursued by the forest planning. Finally, the share of each age class produced by MC is plugged into an intuitive formula that helps the forester decide the share of each annual allowable cut (AAC) that must be put off for the coming year in order to compensate for the EW volume that might occur in the current year.

Material and methodsTop

Pilot area

The proposed algorithm in this study uses data collected from Brosteni FD located in Eastern Carpathians, in Romania (Figure 1). The training dataset refers to the EW occurred between 1999 and 2009, while the testing data refer to the EW occurred after 2009. The main features of the two data sets are summarized in Table 1.

Figure 1. Location of Brosteni FD.

According to the current management plan Brosteni FD covers 21,013 hectares of full-stocked forest, with elevations ranging between 137 and 1650 m a.s.l. Part of the forest (7341.9 hectares) is strictly protected for different reasons (steep slopes with outcrops, patches of old growth forest, biodiversity and conservation). The main species is Norway spruce (Picea abies L.,) encompassing 88% of the area, followed by silver fir (Abies alba) with 4%, beech (Fagus sylvatica, 3%, and some other less important broadleaf species such as birch (Betula alba). The main yield allowable cut for 1999 forest management plan was 55,081 m3 per year; the total length of the forest roads network is 137 km (an average 6.5 m/ha). The forest area restituted to individual forest owners is labeled on the map as ‘no data available’ (Figure 2). Prior to forest restitution, which started in 2000, the total forest area was 30,957 ha.

Figure 2. Age class distribution of Brosteni FD.

Eight years later, when the management plan was updated, the area shrank to 21,032.70 ha and the allowable cut decreased to 36,626 m3 per year. The total current growth of commercial stands, which may be harvested if the age structure is normalized, is 60,352 m3, meaning that the age structure is unbalanced. The total area of commercial forests is 12,941 ha (61% of the total forest area).

Bayesian assessment of EW likelihoods

We actually adopted the windthrow triangle defined by Mitchell (1995), without hypothesizing the relative importance of any enabling factor. Based on literature (Scott & Mitchell 2005, Panferov & Sogachev 2008) and data-screening, we considered that EW depends on the following three factors: the stem taper, or slenderness (Bošela et al., 2014) (symbolized by ‘S’), ‘pits and mounds’ terrain condition (symbolized by ‘T’), and the vicinity to a gap in tree canopy (‘V’), which might be a previous EW. The existence of a canopy gap makes an EW more likely, due to the margin effect (Gadow, 2000) and fragmentation (Holeksa et al., 2017). We calculated stem tapper as ratio between the height, in meters, and the diameter at breast height in centimeters, according to Olofsson & Blennow (2005).

All three factors were conveyed in binomial variables for all stands either affected by EW or not affected at all. Having the contingency table of EW occurred between 1999 and 2008, we calculated the likelihood of having an EW in any compartment, after 2008, in six different ways, considering all permutations of the three enabling factors as follows: STV, SVT, TSV, TVS, VTS, and VST, where each capital letter signifies one of the three enabling factors.

The Bayesian cascade works as follows: the posterior likelihood given by the interaction between wind and the first factor turns into the prior of the interaction of wind with the second factor, and with the third factor respectively. For example, equations 1-3 show how we estimated the EW posterior probability, given the terrain conditions, gap vicinity and slenderness. More precisely, equation (1) gives the posterior of having an EW on pits and mounds (terrain conditions), equation (2) gives the posterior of having EW on pits and mounds and near an existing canopy gap, having the prior given by equation (2), and equation (3) gives the final EW posterior, given terrain condition, vicinity, and slenderness (stem tapper higher than one). The posterior likelihood given by the interaction between wind and the first factor turns into the prior of the interaction of wind with the second factor, and with the third factor respectively.

(1)

(2)

(3)

The symbols have the following meanings: P(W|T) – the probability of EW, given the stand location on pits and mounds (location sensitivity); P(T|W) – conditional probability the stand is located on pits and mounds, given EW; P(W) - windthrow prevalence, or the share of the stands blown down in the last 10 years irrespective to their location, stem tapper or vicinity; P(T) – likelihood that a stand is located on pits and mounds, P(V) – likelihood of being adjacent to a canopy gap (vicinity), P(S) – probability of stem taper higher than one.

The affected and not affected areas were broken down in a contingency table (Table 2), given the EW recorded by FD staff after 1999. The EW taken into consideration occurred on commercial or noncommercial forests (i.e. forests managed for wood production, or for conservation purposes, respectively).

The locations of stands with EW produced in the last decade are presented in Figure 3. The vicinity of a new stand with a prior EW was established in a buffer of 30 m, as it approximates the mean height of a mature Norway spruce stand. The GIS buffers were applied to both 1998 and 2008 GIS maps. Further, the two buffers were intersected, resulting an area that contains the attributes of the two adjacent sub-compartments, with new EW next to prior EW within a maximum distance of 60 m. A screenshot of the GIS buffers is presented in Figure 4. The location of stands with stem tapper greater than one is shown in Figure 5 while the stands located on pits and mounds are shown in Figure 6.

Figure 3. Spatial distribution of EW against age groups.

Figure 4. Buffers used to tag EW vicinity.

Figure 5. Spatial distribution of stands with stem tapper higher than one.

Figure 6. Pits and mounds distribution.

For each sub-compartment, the volume affected by EW between 1999 and 2008 was summed up. The corresponding area was estimated by diving the EW volume to the average growing stock per hectare of each sub-compartment. Even though the patches of EW couldn’t be precisely located, their areas were better estimated. If two or more consecutive windthrows are reported in the same compartment, the area of that compartment goes entirely into the corresponding cell of the contingency table shown in Table 2.

This sketchy manner to operate with data is precise enough to address EW, without any account on catastrophic events. After having calculated six arrays of EW likelihoods, based on the training dataset of decade 1999-2008, we tested the likelihoods on the next decade dataset (Figure 7).

Figure 7. Flowchart of calculating AACR.

For both periods we split the data into three 40-year subsets to take into consideration the age of each stand. Having the contingency table of EW occurred in each age group, and considering all six permutations of enabling factors, we calculated 48 likelihoods of having EW, considering all possible combination of the binary values taken by the enabling factors. In so doing we have computed six different likelihoods of EW for any stand, (one likelihood for each permutation of factors) according to the values taken by the three enabling factors.

Markov Chains

We used the stochastic age class model produced for the first time by Kouba (2002) described as follows:

(4)

where ai is the share of the ith age class into the whole forest, and pi,i+1 is the probability of having the ith age class passed into the following age class (transition probability).

When all transition probabilities are equal to one the transition process is fully deterministic and ai=1/n, for all i=1,n where n is the number of age classes. If all or some transition likelihoods are smaller than one, then the age structure is skewed and the sum of the ai is smaller than one, meaning that at any moment a certain portion of the forest is waiting for being replanted after a disturbance (in our case, an EW).

Adjusting the transition probabilities of age classes to forecast accuracy

Appraising the transition likelihood for each age class was a real challenge, for the many reasons that have been already discussed in the introductory section. A Bayesian method was more tempting because it allows a step-by-step assessment of the windthrow likelihood, based on the prior events produced in the last decade. However, the forecast accuracy shall be considered, in order to avoid any overestimation of EW likelihoods. Therefore, for each permutation of enabling factors, within each age group, we calculated the Receiver Operating Characteristic (ROC) and the corresponding Area Under Curve (AUC), to adjust the likelihoods produced by Bayesian inferences, as equation (5) shows:

(5)

The ROC curve is the geometrical position of the cumulative frequencies of true positives (TP) and false positives (FP), after having the posterior likelihoods of EW ranked in descending order. Based on the coordinates of ROC points, we further calculated the AUC having true positives (TP) and false positives (FP) estimated by equations (5) and (6).

(6)

(7)

Where Si is the area of the ith stand, pi is the poste­rior likelihood of EW, and fi is the real state of that stand (1 =affected by EW in the last decade; 0 –not affected by EW).

The annual AACR vs. the cutting budget to EW occurrence

Having an AAC given by the forest plan and the expected returns to the first age class from older stands we can estimate the AACR lagged for the next year (R) by equation (8).

(8)

where si is the area of the ith age class (i=1,..,n), and gi is the current growth per year and hectare for the ith age class. When all ai equal 1/n, meaning that no EW had occurred, the AACR is zero, and the whole AAC can be entirely harvested in the current year.

As already mentioned, the SYP has been adopted long ago and it means more than equal yields year after year. Therefore, two different methods are used to calculate AAC and both methods take into account any likely short supply of wood that may occur in the next 10, 20, 40 and 60 years. By “short supply” we mean the lowest ratio between the total volume that will be harvestable in a given period of time and the length, in years, of that period (next 10, 20, 40 or 60 years).

The whole algorithm has been described at length in Dragoi (2010). The indicative growth method we have referred to, makes use of a ratio labeled with q, which makes the difference between two different situations: too much yield in the next 20 years, due to numerous mature and over-mature stands (q>1) or too few mature stands (q<1). For each situation different algorithms are being used. For Brosteni FD the second situation applies (deficit of mature stands) and there is no risk of triggering the positive feed-back we highlighted in the introduction.

If the age class distribution were imbalanced and both first and third age groups exceed the second one (the age-class structure of commercial forests resembles a “U” shape), it means that many mature stands were put off over and over after each salvage (i.e. EW) cutting. In such circumstances, indicating a positive feed-back loop, the AAC shall not be diminished with the AACR, because the AAC has long stopped steering the forest to a normal structure, where all age classes are, more or less, even.

When a higher AACR is the norm, not the exception, the “deterministic” age structure will never be balanced but skewed, and it makes no sense to pursue a normal age-class structure, without any account of EW occurrence.

Adjusting the demand for seedlings

Assuming a deterministic model, and no shelterwood system (i.e. no natural regeneration), every ten years the demand of seedlings shall cover the corresponding area of mature stands clear-felled before. In case of EW, the total area which shall be regenerated within each decade is larger, because not only mature stands were harvested but also the damaged trees and the damaged portions of the younger stands. If all pi,i+1 are smaller than one, the sum of all ai given by equation (4) is also smaller than one, and the additional area that shall be replanted is given by equation (9):

(9)

Where P is the additional area to be planted in the next 20 years, F is the total area of commercial forests.

ResultsTop

Adapting cutting budget to EW occurrence

The contingency table used to calculate priors for every permutation of factors is presented in Table 2. Table 3 shows the posterior likelihoods considering all permutations of enabling factors (slenderness, gap vicinity, and terrain conditions) and all possible combinations of values taken by the three enabling factors. Having these probabilities stored in a matrix, we assigned to each sub-compartment included into the testing dataset six probabilities, one for each permutation of enabling factors, considering the age group wherein each sub-compartment falls. The AUC for each permutation and age groups are presented in Table 4 while the EW likelihoods, weighted on area for each age class are shown in Table 5. The AUC for all permutations within the third age group are zero, meaning that posterior likelihoods of having EW are not reliable at all for further calculations. All input variable plugged into equation 8 are shown in Table 6 where the AACR is shown in the last row.

As aforementioned the total AAC of Brosteni FD is 36,626 m3 and according to the current regulations, the formal reserve for every next year is 7325 m3. According to the Markovian model, an AACR of 2,696 m3 is enough to compensate for the effect of EW on sustained yields.

Assessing the extra-demand of seedlings

The gaps left into the forest after EW must be regenerated within a two year period after having hauled all the affected trees. The actual management planning system has no allowance for replanting the gaps produced by EW. For our pilot areas, having the EW likelihoods, the additional area that shall be replanted every ten years is 336 ha, which means an average 33.6 ha/year. Having a norm of 5000 seedlings/ha, the additional number of seedlings per year reaches 168500 plants.

DiscussionTop

The key issue discussed here is whether or not the sheer Bayes’ rule is a good option for making statistical inferences. The Bayes’ rule assumes interchangeable roles between hypotheses and evidences (Nöth & Weber, 2003). This assumption is counterintuitive because most of the people are inclined rather to think in terms of cause and effect; once this epistemological difficulty is overcome things become clearer if we associate the hypothesis with a causal chain strengthen by more evidence, which turns into a new hypothesis, and so forth (Wintle & Lindenmayer, 2008).

Based on the old concept of windthrow triangle (Mitchell, 1995) we considered just three enabling factors, in all possible orders conveyed into six permutations. Having one more attribute, say the position on the slope, given that in many cases the stands located on foothills are blown down by the wind coming from the opposite direction, the number of all possible permutations (i.e. hierarchies of factors to be tested) rises to 24. From the computational point of view this is not big challenge because all calculations can be plugged into spreadsheet and the input data are produced by a simple cross-tabulation, provided that all important attributes had been included into the GIS table of attributes (deFillipi et al., 2004).

A 30 m wide buffer may indicate as adjacent two stands separated by a 59.99 m wide strip of trees which were not yet blown down, if the two stands affected by EW are separated by a stand of trees narrower than 60 m. Such situations are represented by very thin buffers, some of them being screen-captured in Figure 4 and marked with yellow circles. Hence, prior EW vicinities might have been overestimated in some cases, but such situations can instantly be spotted.

For the first age group of stands the best predictability corresponds to a seemingly unexpected permutation of factors: vicinity-terrain-slenderness. It is quite surprising because most of the young stands are not actually thinned, and slender. Even more surprisingly, slenderness plays the most important role for the EW produced in middle-aged stands (Konopka et al., 1987). For the last age group none of the permutations of factors proved to be a serious threat, considering the combinations of EW and values taken by the three enabling factors.

We tested the accuracy of predictions with the ROC-AUC procedures. The main disadvantage of gauging any prediction made by discrete-output procedures – and Bayes’ rule turns discrete likelihoods – is the ROC shape (Schindler et al., 2009). All discrete values taken by EW likelihoods are shown in Table 3 but, associated to the stands, they appear in lesser combinations, depending on how numerous are the stands that fall into one of the three age groups. Yet the most complex situations are characterized by eight likelihoods, only the ones with likelihoods grated than 0.5 are being considered, as shown in Figure 8. The first one (a) was the most complex (with four discrete values) but most of other ROC look like (b), with three discrete probabilities and (c), with only two discrete values.

Figure 8. Four typical ROC curve for discrete EW likelihoods.

The fourth situation (d) typifies all ROC curves for the last age group. In Figure 8 d) we have only two true positives situations, without any false positives. Because the cumulative frequencies are plotted against the two axis, summed up from the lowest to the highest likelihood, the two pints are vertically aligned to the right, having no area below. However, the loss of information is not important because the EW likelihood of the last age class must be forced to one, because all mature stands are supposed to be harvested within the next two decades.

If the young stands had higher growths, the AACR would have been higher: an average of 3.7 m3/year/ha, as Table 6 shows, is a hint of having unstocked stands in the first two age classes.

By cascading the Bayes tree probabilities for each development stage, we assumed that it must be an order in which all three enabling factors weaken most of the trees, although the order is yet unknown. Hence, we had to conceive all possible orders of factors and see which ones behave better in terms of accuracy of prediction.

Our case study showed that what really matters is not the predefined set of enabling factors, but the number of cases analyzed in order to train the model. Having more cases in the first age group, four ROC, out of six, are meaningful (greater than 0.5), while the second age group has just two significant ROC (for STV, and VST respectively, see Table 4).

We have confined our study to two periods of ten years each to match the data required by the model with the data supplied at no additional cost by the Romanian forest planning system. Worth noting, for Brosteni FD, as well as for any public FDs in Romania, the last and the last but one forest management plans have been produced in GIS technology, allowing cheap data processing; consequently, we couldn’t go back too much in the past. In other conditions, having a greater amount of data and more GIS maps for the same area.

ConclusionsTop

This study is the first attempt in calculating AACR, which shall be put off for the next year, given a certain pattern of EW (or whatever disturbance) for Romanian forest planning. The study demonstrated that 10% AACR is enough in mountainous area, where EW are quite frequent. The 10% or 20% of the AAC left unharvested from the previous year, and the same amount of wood left uncut for the coming year is not, allegedly, a substantial improvement of forest management. However, a slimmer AACR is advantageous for logging companies who are bidding for harvesting contracts every year, because they can make more combinations with the tracts worth bidding for, simply because the share of the main yield auctioned every year by the forest managers is 10 % greater.

If the SYP fails to normalize the age-class structure of the forest due to the positive feed-back loop previously described, the model allows a better control on the total yield, and eventually can be used to improve the calculation of AAC, considering a slant age-class structure, not an even one. Having better predictions of EW, knowing their future locations (elevation, slope and aspect) and the hydrological effects of a clear-felling, the forest planners will be able to integrate this information into more complex decision-making models, able to better balance the wood production with water runoff management.

We have considered only three enabling factors, but the list is open: besides slenderness, terrain conditions and EW vicinity, other factors can be added, such as vertical structure of the stands, presence of wind-resisting species of trees, or their location (uphill or downhill). The hints for taking into consideration a new factor are the values stored into the contingency table; and the more permutation of enabling factors are at hand, the more curved ROC will be, allowing better predictions of further EW.

AcknowledgmentsTop

We are deeply grateful to the two peer reviewers for their supportive suggestions. We also acknowledge the staff of Brosteni FD who provided us the input data.