Introduction

A crucial goal for contemporary forestry is to develop management methods capable of combining economic sustainability, maintenance of biodiversity and adaptability of forest ecosystems to on-going environmental changes (Messier et al. 2014; Bauhus et al. 2017). An increasing understanding of these challenges has led to growing interest in the functioning of complex natural ecosystems. Such ecosystems are usually characterized by higher structural heterogeneity, which is crucial for habitat diversity as well resilience against disturbances and recovery potential (Röhrig et al. 2006; Brang et al. 2014). Nevertheless, after the long history of human colonization in Europe, only relicts of the former pristine forests remain. Among these remnants, one of the best represented is montane mixed-species forests consisting of European beech (Fagus sylvatica L.), silver fir (Abies alba Mill.) and Norway spruce (Picea abies (L.) H. Karst.), which can still be found in fragmented patches in the Carpathians, Balkans and Alps (Mayer 1986; Korpeľ 1993; Veen et al. 2010).

Until the 1990s the leading view on the dynamics of natural montane mixed-species forests was the forest cycle concept, which assumed that their structural heterogeneity is an effect of the coexistence of integral and spatially distinguishable stand patches of a size corresponding to areas occupied by several canopy trees and representing different phases of a linear developmental cycle (regeneration–growth–maturity–senescence–breakdown and back to regeneration) (Korpeľ 1993; Leibundgut 1993). Because research was conventionally conducted on single study plots selected subjectively to characterize a given type of stand structure (e.g., developmental phase), and therefore not representative for the stand as a whole (Peck et al. 2015), or on small widely spaced inventory plots that did not allow for multi-scale inference about structural heterogeneity, effective testing of the forest cycle and inherent patch–mosaic hypotheses was not possible. However, recently, large-scale surveys and sophisticated analytic tools have been employed to gain insight into these issues. The view emerging from these studies is that the dynamics of natural montane mixed-species forests is driven by disturbances of minor intensity that kill single trees or lead to the death of small groups of trees (Standovár and Kenderes 2003; Splechtna et al. 2005; Zeibig et al. 2005; Firm et al. 2009; Trotsiuk et al. 2012; Motta et al. 2015; Petritan et al. 2015; Orman and Dobrowolska 2017). Severe disturbances occurring in the form of insect outbreaks, wind and/or ice storms are infrequent (Nagel and Diaci 2006; Jaloviar et al. 2017; Janda et al. 2017), and the pattern of live canopy tree distribution and their mortality is close-to-random (Szwagrzyk and Czerwczak 1993; von Oheimb et al. 2005; Paluch 2007; Šebková et al. 2011; Paluch et al. 2015). As a result, canopy heterogeneity is usually characterized by a reverse J-shaped gap size frequency distribution, with larger gaps usually containing retained trees and gap expansion playing only a minor role (Kenderes et al. 2009; Kucbel et al. 2010; Nagel and Svoboda 2008; Standovár and Kenderes 2003; Zeibig et al. 2005; Nagel and Diaci 2006; Parobeková et al. 2018). Multiple releases followed by periods of suppression have frequently been found in the life histories of canopy trees (Splechtna et al. 2005; Zeibig et al. 2005; Nagel and Svoboda 2008; Firm et al. 2009; Trotsiuk et al. 2012), which underscores the complexity of the gap-filling process. Studies by Nagel and Diaci (2006) and Nagel et al. (2014) also suggest also that a portion of temperate forest landscapes is in some stage of recovery from intermediate-severity wind damage and thus represents rather non-equilibrium conditions. A consequence of recurrent light or intermediate disturbances is the widespread intermixture of small and large trees. Although homogeneously structured patches occur, they cover a smaller proportion of the total forest area than those that are heterogeneously structured (Commarmot et al. 2005; Motta et al. 2011; Šebková et al. 2011; Král et al. 2014; Paluch et al. 2015; Drössler et al. 2016; Parobeková et al. 2018). Structural differentiation occurs primarily at the finest spatial scales and diminishes rapidly (Král et al. 2010a; Alessandrini et al. 2011; Zenner et al. 2019). For example, using local deadwood and live tree distributions as discrimination criteria, Král et al. (2014) found mean patch size in a narrow interval between 570 and 800 m2. Drössler et al. (2016) indicated that homogeneous patches cover only a small proportion of the total forest area compared to heterogeneously structured stand fragments. In beech-dominated forests Glatthorn et al. (2018) showed that the point clouds of structural data are rather homogeneous and lacking clearly separated clusters interpretable as development stages of the forest development cycle and concluded that any discrimination between development stages relying on stand structural data requires splitting a continuous point cloud. Moreover, the forest cycle concept based on a temporal sequence of development phases (Korpeľ 1993; Leibundgut 1993) has been found to be inadequate compared to the real developmental trajectories of stand patches. The study by Král et al. (2018) suggests that over the long term, the proportion of acyclic transitions is significantly more frequent than cyclic transitions along theoretical pathways, and transitions among different multi-aged phases dominates. Thus, frequent small-scale disturbances followed by quick recovery results in development pathways that resemble a complex web with highly stochastic connections rather than a deterministic repeating cycle.

Despite the considerable efforts that have been undertaken to elucidate the structural attributes of the natural mixed-species forests typical of the mountainous regions of Central and Southeastern Europe, knowledge about inter-regional variation in the structural heterogeneity of these stands remains fragmentary. The substantial variation in necromass and biomass accumulation reported in the literature (Holeksa et al. 2009) raises the question how the differences in these cumulative characteristics are associated with the intra-stand pattern of structural heterogeneity. The Dinaric Mts. harbour one of the most volume-rich mountain primeval forests on the continent, with growing stocks commonly exceeding 1000 m3 per ha. In contrast, in the Western Carpathians, which are the northernmost region of occurrence of such forests composed of silver fir, Norway spruce and European beech, average stand volumes are only half of those reported in the Balkans. In previously published studies from the Western Carpathians (Paluch 2007; Paluch et al. 2015), a fine-scale clustering of under-canopy trees but close-to-random pattern of canopy trees was reported in a wide range of spatial scales ranging from an area corresponding to an average parcel of a canopy tree up to several hectares. In addition, it was found that distributions of local structural differentiation and basal areas of live and dead trees rapidly approach a symmetric unimodal distribution type as the spatial scale increases.

In this study, we hypothesized that the considerably higher level of biomass accumulation in Dinaric forests would be reflected in a different pattern of their structural heterogeneity. Specifically, we expected a higher proportion of stand patches with a high level of biomass accumulation and low structural differentiation, resulting in a strongly skewed distribution of basal areas of live trees. In addition, taking into account probable differences in disturbance regimes and mortality rates among canopy trees, we anticipated a lower ratio of dead to live tree basal area in the Dinaric region compared to that in the Western Carpathians.

Methods

Study area

The study was carried out in three close-to-primeval forests in Bosnia and Herzegovina, which are located in the central part of the Dinaric Mts. in Southeastern Europe. The investigated forests are situated in the northwestern (Lom reserve, 44° 27′ N, 16° 29′ E), middle (Janj reserve, 44° 08′ N, 17° 18′ E) and southeastern part (Perućica reserve, 43° 18′ N, 18° 43′ E) of the central Dinaric Mts. They have been officially strictly protected since the 1950s, but there is also evidence that strict nature protection in these areas goes back to the nineteenth century during the Austro-Hungarian Empire (Stupar and Milanović 2017). The core areas of Janj (57.2 ha) and Lom (55.8 ha) are surrounded by relatively large buffer zones (237.8 ha and 297.8 ha, respectively) in which only low-intensity salvage cutting has been performed. The strictly protected area in Perućica is much larger (1493 ha) and includes different forest associations.

One close-to-rectangular research area about 10 ha in size was established in each reserve in a mixed-species stand consisting of European beech, silver fir and Norway spruce. The selected stand fragments were as homogeneous as possible in terms of site conditions and terrain orography and showed no signs of direct human activity. All studied stands were situated in the montane belt between 1250 and 1480 m a.s.l. on western (Janj), northeastern (Lom) or southwestern (Perućica) slopes with inclinations between 5 and 25°. The mean annual temperature in the studied stands is 5–7 °C, and the annual precipitation ranges between 1400 and 1900 mm. Dolomite bedrock is mostly present in Janj, limestone predominates in Lom and Perućica (at the observed location), and brown soils prevail in all three studied stands. Tree species composition in the study area is shown in Table 1.

Table 1 Characteristics of live trees in the three stands studied
Full size table

The data from the Western Carpathians used later in the comparisons were collected in five forest reserves containing the best maintained remnants of beech–fir (–spruce) primeval forests in this region: Baniska (49º 27′ N 20º 37′ E), U źródeł Solinki (49º 6′ N 22º 33′ E), Łabowiec (49° 28′ N 20° 49′ E), Oszast (49° 25′ N 19° 11′ E) and Żarnówka (49º 35′ N 19º 33′ E). The research plots were situated in the upper zone of the lower montane belt between 900 and 1040 m a.s.l. The vegetation period is approximately 180–200 days, the mean annual temperature is 5–7 °C, and the annual precipitation is 900–1450 mm, reaching a maximum in the summer months. The prevalent soils are loamy dystric or eutric cambisols developed on flysch bedrock (sandstones and shales). The tree cover consists of beech, fir and spruce. More detailed data on the site and stand characteristics of these research objects can be found in previous papers (Paluch 2007; Paluch et al. 2015; Paluch and Bartkowicz 2020).

Field measurements

To enable comparability of results, the present study employed the same sampling scheme and computational approach used in the project conducted in the Western Carpathians (Paluch et al. 2015). In each research area, a grid of perpendicular axes separated by 20 m was established, and two concentric sample plots with radii of 7 and 15 m (and areas of 154 m2 and 708 m2, respectively) were placed on each of the 244–256 grid nodes. Within the smaller plot, the diameter at breast height (dbh) and species of all live trees of dbh ≥ 7 cm and dbh of all dead trees of an approximated dbh ≥ 20 cm were recorded together with their species when possible. Within the larger plot, the number and species of all live and dead trees of dbh ≥ 50 cm were registered. The dead trees were assigned to the sample plots based on the midpoint of their stem bases. Tree dbh was determined from two measurements (along and across the slope) with a calliper exact to the nearest 0.5 cm. In the case of missing bark, advanced decay or split logs, the dbh of dead trees had to be approximated. For each species, at least 100 height measurements in the entire range of dbh classes were carried out to 0.1 m accuracy with a Vertex IV hypsometer (Haglöf, Sweden). The dead tree definition included both standing and lying dead trees. These dead trees were divided into three decay categories: class A—recently dead trees (within ca. 5 years) that retained thin twigs and intact, adhering bark; class B—trunks partly without bark with slightly decomposing but still compact wood and branches of diameter ca. 10 cm present; and class C—snaps, fallen snaps and remnants of uprooted trees with soft, friable wood in advanced decay, with completely decomposed crowns and distinction between conifer species often impossible. Remnants that were almost entirely decomposed, with a blurred trunk outline, unidentifiable stem base position and recognizable only by ground microtopography, were discarded.

Data analysis

The measurements allowed determination of the number and basal area of live and dead trees on the plots (hereafter, BAL and BAD, respectively). As a measure of structural complexity, the empirical variance in tree heights was computed and scaled through comparison with a hypothetical variance of the uniform distribution (hereafter, the differentiation index, VARh, see Paluch et al. 2015 for more details):

$${\text{VAR}}_{\text{h}} = \tfrac{{S_{\text{h}}^{ 2} }}{{S_{\text{U}}^{2} }},$$

where S 2h —empirical variance in tree heights, S 2U  = (hmax − hmin)2/12—theoretical variance of the continuous uniform distribution for given minimal and maximal values, hmin and hmax. The uniform distribution describes an experiment where arbitrary outcomes lie between certain bounds and are equally probable. Tree heights (h) were approximated by the regression of the general form h = dbh2 (a0 + a1dbh+ a2dbh2)−1 + 1.3 with the coefficients fitted separately for each research area and species. The uniform distribution used in the calculation of the index was constrained by the values of the 1st and 99th percentile of the smoothed tree heights in the dataset. The aim of using smoothed tree heights (and not tree dbh directly) was to put more emphasis on the absolute dbh differentiation among smaller trees, which more strongly contributes to vertical canopy heterogeneity than dbh differentiation among large trees. The VARh index takes the non-negative values and attains the value of 1 for dbh structures corresponding to a uniform height distribution bounded by hmin and hmax and the values lower/higher than 1 for the under- or over-dispersed structures, respectively. Figure 1 displays examples of the differentiation indices VARh calculated for different dbh distributions.

Fig. 1
figure 1

Values of differentiation index (VARh) depending on tree frequencies in dbh classes (examples). Tree heights (h) used in the calculation of the index were approximated as h = dbh2(2.4472 + 0.9727dbh+ 0.01142dbh2)−1 + 1.3. The VARh index attains values close to 1 for dbh structures corresponding to a uniform height distribution (bounded by predefined minimal and maximal height values) and values lower/higher than 1 for under- or over-dispersed structures, respectively. In terms of the common definition of stand vertical structure, low index values indicate one-storey types, close to 1—multi-storey types, and significantly higher than 1—two-storey types

Full size image

To account for the fact that the distribution of trees of different sizes is potentially driven by partially distinct ecological processes or the same processes operating at a different (size-dependent) scale, the analysis was conducted separately for several dbh classes. For methodological coherence with the former study (Paluch et al. 2015), trees were divided into canopy trees (dbh ≥ 50 cm), mid-canopy trees (25 cm ≤ dbh < 50 cm) and under-canopy trees (7 cm ≤ dbh < 25 cm). In addition, the category of large canopy trees (dbh ≥ 75 cm) was also distinguished.

The spatial distribution patterns of trees were analysed using a two-step procedure. First, the dispersion index I was computed as the variance-to-mean ratio for the tree counts xi on the sample plots with radii of 7 and 15 m (i.e. for the spatial scale of 154 and 708 m2, respectively),

$$I = \left[ {\mathop \sum \limits_{i = 1}^{n} \left( {x_{i} - \hat{x}} \right)^{2} } \right]/\left[ {\hat{x}\left( {n - 1} \right)} \right].$$

The product I (n − 1), where n denotes the number of sample plots, was then used as a test statistic in a random labelling tests in which the recorded trees were randomly assigned to the sample plots in 1000 simulations. In the case of the larger plots (708 m2) which slightly overlapped for neighbouring grid nodes separated by 20 m, the empirical and simulated dispersion indices were calculated as means from two estimates obtained for two disjunct subsets with the nearest neighbour plots separated by \(20\sqrt 2\) m.

The second step was designed to take advantage of the known location of the sample plots. It adopted a paired-plot approach (Cressie 1991) in which the mean square for the plots separated by s grid units (= 20 m) is

$$V_{\text{s}} = \frac{1}{{n_{\text{s}} }}\mathop \sum \limits_{i = 1}^{{n_{\text{s}} }} \left( {A_{i} - B_{i} } \right)^{2} .$$

here Ai and Bi are the number of events in the ith pair of the plots s grid units apart, and the summation is over all distinct ns pairs of the plots separated by s units. Given the shape of the research areas, 200 m (i.e. s = 10) was the maximal inter-plot distance analysed. Under a random arrangement of plot counts, the variance estimates could be expected to be roughly constant (Cressie 1991). Low/high mean square values indicate under-/over-dispersion of attribute values and relative similarity/dissimilarity between the plots separated by s units (i.e. positive/negative spatial dependence). The empirical \(V_{\text{s}}\) values were compared to 95% simulation envelopes for 1000 randomisations of the plot counts. In this procedure the plot counts were randomly assigned to the sample plots, then \(V_{\text{s}}\) was calculated and the lower and upper simulation envelopes designated as the 25th and 975th based on the ascending order of \(V_{\text{s}}\) values. As a measure of the discrepancy between the empirical variance estimates \(\hat{V}_{\text{s}}\) and the expectations under the random arrangement of plot counts \(\bar{V}_{\text{s}}\), the following \(U\) measure for spatial dependence was defined:

$$U = \frac{1}{k}\sum \frac{{2\left( {\bar{V}_{\text{s}} - \hat{V}_{\text{s}} } \right)}}{{W_{\text{s}} }},$$

where the summation is over all s units for which \(\hat{V}_{s}\) is beyond the 95% simulation envelope, k is the number of s units in the sum, \(\bar{V}_{\text{s}}\) is the mean obtained from the simulations and \(W_{\text{s}}\) is the width of the estimated 95% interval under random plot assignment. Thus, U-values close to 0 denote a random pattern, and those greater than 1 or lower than − 1 denote a significant positive or negative spatial dependence between plot counts. As regards BAL, BAD and VARh, a paired-plots approach was used in which the plot counts were replaced by the basal area of live/dead trees or the differentiation index, respectively.

The empirical spatial patterns of the BAL, BAD and VARh were applied in simulations designed to model and compare expected variation of these variables depending on the spatial scale. In the two-step simulation procedure, first a random sample plot Xij with row and column indices i and j was drawn from the set of all sampling plots, and then N sample plots corresponding to the spatial scale being tested were randomly drawn (with replacement) from the surroundings of Xij defined by row and column indices m and n fulfilling |i − m| ≤ k and |j − n| ≤ k, respectively. The constant integer k was fitted to cover the spatial range of the analysis, and a toroidal shift method was applied to control for border effects. For full comparability, exactly the same simulation procedure was used for the data from the Western Carpathians and Dinaric Mts.

Results

Species composition, basal area and diameter distributions

In the three stands in the Dinaric Mts., the proportion of beech in tree number ranged between 32 and 72% and was considerably higher than its percentage in BAL (19–38%, Table 1). The dominance of beech over fir and spruce was particularly visible among mid-canopy trees (25 ≥ dbh < 50 cm) in the Lom reserve and among under-canopy trees (dbh < 25 cm) in the Janj and Lom reserves (Fig. 2). In both reserves the percentage of beech among canopy trees (dbh ≥ 50 cm) was insignificant. Overall, despite substantial differences in species composition, the stem densities and BAL values were similar (430–569 trees per ha and 56.8–65.9 m2 ha−1, respectively; Table 1). The dbh distributions for all the species pooled together had a negative exponential (Lom) or rotated sigmoidal form (Janj and Perućica), with the highest number of small trees and a slightly marked local maximum in dbh classes of about 60–75 cm (Fig. 2).

Fig. 2
figure 2

Diameter distributions of live trees in the three stands studied in the Dinaric Mountains

Full size image

The stem density of dead trees (with dbh ≥ 20 cm) ranged between 57 and 81 per ha, BAD values were between 18.0 and 24.9 m2 ha−1 (Table 2) and the dbh distributions had a rotated sigmoidal form (Fig. 3). The basal area of dead trees tended to increase with increasing proportion of coniferous species among live trees and was higher in the Janj and Lom than in Perućica reserve. The percentage of beech was considerably lower in BAD than in BAL (Tables 1, 2).

Table 2 Characteristics of dead trees in the three stands studied
Full size table
Fig. 3
figure 3

Diameter distributions of dead trees in the three stands studied in the Dinaric Mountains

Full size image

Distribution patterns of live and dead trees

In the single plots with a radius of 7 m (154 m2), the average number of live trees (dbh ≥ 7 cm) varied between 6.6 and 8.8. The 5th percentiles of the stem number were between 2 and 4, and the 95th percentiles between 12 and 15, respectively. At the spatial scale of single plots (154 m2), the distributions of all live trees deviated from randomness towards clumped patterns (Table 3). However, the effect of aggregation was attributable only to under-canopy trees (7 ≤ dbh < 25 cm), while the distribution patterns of larger trees were random (for diameter classes 25 ≤ dbh < 50 cm and 50 ≤ dbh < 75 cm) or regular (for trees of dbh ≥75 cm) (Table 3). The analysis carried out for live canopy trees (dbh ≥ 50 cm) in the plots with radii of 7 and 15 m revealed repulsion effects for both scales in the Lom and Perućica reserves and a shift from a regular towards a random pattern with increasing spatial scale in the Janj reserve (Table 3).

Table 3 Spatial pattern of live and dead trees on the 154 m2 circular sample plots (variance-to-mean ratios of the tree counts)
Full size table

An under-dispersion (i.e. positive spatial dependence, patchiness) of live tree densities (dbh ≥ 7 cm) was found in all the stands, and this tendency was detectable up to a distance of 20 m in the Lom and Perućica reserves and maintained up to 80 m in the Janj reserve (Table 4). This pattern resulted from under-dispersion of under-canopy trees observable in all the stands but also of medium-diameter trees (25 ≤ dbh < 50 cm) in the Janj and Lom reserves. For canopy trees (dbh ≥ 50 cm), the results were divergent: the densities were random in the Lom and Perućica reserves and under-dispersed at scales above 80 m in the Janj reserve. Nonetheless, the densities of larger trees with dbh ≥ 75 cm in all the stands were spatially independent (Table 4).

Table 4 Spatial pattern of live and dead trees densities in the distance range between 20 and 200 m: deviations between the empirical paired-plot variance estimates and expectations under a random arrangement of trees defined by the U value
Full size table

In the single plots (154 m2), the average number of dead trees (dbh ≥ 20 cm) was between 1.0 and 1.3. At this spatial scale, the distribution of dead trees from all decay classes pooled together was random in the Janj and Lom reserves and aggregated in the Perućica reserve (Table 3). Dead canopy trees (dbh ≥ 50 cm) also displayed the same types of distribution. For the larger scale, i.e. on plots with a 15 m radius (708 m2), the distribution of dead canopy trees (dbh ≥ 50 cm) was aggregated in all the stands, and this effect was attributable to dead trees representing different decay classes (Table 3). Deeper analysis indicated that dead canopy spruces and firs, but not beeches, exhibited the aggregated type of distribution (with variance-to-mean ratios between 1.15 and 1.67). At spatial scales between 20 and 200 m, the densities of dead canopy trees were random in the Lom and Perućica reserves and showed a slight tendency towards under-dispersion (i.e. positive spatial dependence) for distances less than 40 m in the Janj reserve (Table 4).

Spatial variation in basal area of live and dead trees

At the single plot scale, the frequency distributions of BAL took the form of a truncated bell-shaped distribution (Fig. 4) with a maximal frequency at 0.84 m2 per 154 m2. As the plot area increased, the BAL distributions rapidly approached a symmetric unimodal type and the spatial variation diminished (Fig. 4). At spatial scales between 20 and 200 m, the spatial variation of BAL (for trees of dbh ≥ 7 cm and ≥ 50 cm) was diverse: over-dispersed for plots separated by 20–40 m in the Janj reserve, random in the Lom reserve and under-dispersed for plots separated by 20 m in the Perućica reserve (Table 5).

Fig. 4
figure 4

Frequency distributions of the basal area of live trees (BAL) and dead trees (BAD) on the sample plots (154 m2) and obtained from random sampling for 1000 m2 and 10,000 m2 plot sizes (Dinaric Mountains). For comparison, the generalized BAL and BAD distributions from the Western Carpathians are given (Paluch 2007; Paluch et al. 2015)

Full size image
Table 5 Spatial pattern of the basal area of live trees (BAL) and dead trees (BAD) and index of vertical differentiation (VARh) in the distance range between 20 and 200 m
Full size table

At the stand level the BAD distribution substantially differed from the BAL distribution (Figs. 2, 3). At the single plot scale (154 m2), the highest proportion had the plots with low BAD values, below 0.05 m2ha-1. With increasing scale, the BAD distributions evolved towards a symmetric unimodal distribution (Fig. 4). In all the stands, the spatial variation of BAD (dbh ≥ 20 cm) was random, and a positive spatial dependence of dead canopy tree density (dbh ≥ 50 cm) was found only in the Janj reserve at scales between 20 and 40 m (Table 5).

Structural differentiation and its spatial variation

The frequency distributions of VARh on single sample plots (154 m2) had a similar, more or less asymmetric unimodal form (Fig. 5). Plots with index values in the range between 0.45 and 0.85 exhibited the highest frequency, indicating the dominance of patches of complex stand structures. The proportion of single plots with low tree height variation and close to a one-layered structure (i.e. with VARh < 0.20) was between 7 and 9%. However, such structures formed incidentally already for a patch size of 1000 m2, and for 10,000 m2 plots only complex structures close to unimodal or bimodal types were likely to occur (Fig. 5). In all the reserves the spatial pattern of VARh index values was random (Table 5). The variation coefficient of VARh decreased with increasing plot size at a similar rate (Fig. 6).

Fig. 5
figure 5

Frequency distributions of the differentiation index (VARh) on the sample plots (154 m2) in the Dinaric Mountains and obtained from random sampling for 1000 m2 and 10,000 m2 plot sizes. For comparison, the VARh distributions from the Western Carpathians are given (Paluch 2007; Paluch et al. 2015). Because of the inconsistent pattern of VARh distributions, empirical frequencies from the individual stands are shown (stands I–V)

Full size image
Fig. 6
figure 6

Variation in the basal area of live trees (BAL), dead trees (BAD) and the differentiation index (VARh) dependent on spatial scale in the Dinaric Mountains. For comparison, the relationships between variation coefficients and plot areas obtained in the Western Carpathians are given (mean value and 95% confidence interval for regression). Please note the logarithmic scale on the X axis

Full size image

Discussion

Methodological remarks

Most variables in natural environments display gradients, patches, trends or other complex patterns which can exist at many scales and correspond to the physical features of the environment or result from ecological processes. Perception of these spatial structures may be critically influenced by the sampling design and techniques used for assessing dispersion patterns (Bellehumeur et al 1997; Dungan et al. 2002). Nonetheless, there is no single “correct” or “optimal” scale for characterizing spatial heterogeneity, and comparison between landscapes using pattern indices must be based on the same spatial resolution and extent or better, involve multi-scale analysis (Wu 2004).

Although they are inferior to methods based on a complete census of individuals, the analytical methods used in this study (variance-to-mean ratio and paired-plot variance) have been recognized as powerful tools for classification of spatial patterns for incompletely mapped data (Goodall and West 1979; Perry and Mead 1979; Perry et al. 2002; Dale and Fortin 2014). The simulation experiment adjusted to the sampling design of this study confirmed the high statistical power of the approach based on the paired-plot variance to detect patchy structures (see Online Resource 1). Nonetheless, the sampling method applied (circular plots of a radius of 7 m in a 20 × 20 m grid) filtered out the spatial variation occurring at scales smaller than the size of the circular sample plots (154 m2). In addition, the spacing of neighbouring sample plots (20 m) also left a margin of uncontrolled dispersion patterns for individuals with dbh < 50 cm. As a consequence, there is the potential to fail to recognize spatial structures at a distance range corresponding to between-plot lags or more subtle spatial interactions (e.g. repulsion changing over to attraction with distance) detectable by analyses conducted in a continuous gradient of spatial scales. On the other hand, it seems that the hierarchical sampling of live and dead canopy trees (with dbh ≥ 50 cm) on regularly spaced plots with radii of 7 and 15 m in this case provided a solid basis for the assessment of dispersion patterns. Moreover, the range of spatial scales involved in our analyses (i.e. 154 and 708 m2 for radii of 7 and 15 m, and at larger scales adequate for the gradient of inter-plot distances of the 20 × 20 m sampling grid) covers the probable resolution of the patch mosaic pattern identified in former studies in similar forest associations (Král et al. 2010a, 2014; Alessandrini et al. 2011; Zenner et al. 2019).

Tree distribution patterns

The analysis showed that the distribution patterns of live canopy trees tend to shift from regular to random or aggregated types with increasing spatial scale. The repulsion effects consistently found for canopy trees at the scale of the single plots (154 m2) seem attributable to increasing competition for resources between individuals sharing the same space and are well documented in the literature (e.g. Stoll and Bergius 2005; Getzin et al. 2006). However, at larger scales (corresponding to inter-plot distances of 20 m or greater) the outcomes revealed different distribution patterns of live canopy trees: coarse-grained patchiness (positive spatial dependence) in Janj (for inter-plot distances above 80 m), random in Lom and fine-grained patchiness in Perućica (for inter-plot distances below 40 m). Moreover, taking into account basal areas instead of tree counts smoothed the spatial variation, which, in such cases, was classified as over-dispersed and random in Janj and Lom, and as under-dispersed (patchy) in Perućica, but only for inter-plot distances of 20 m. These inconsistent results suggest that the spatial distribution of canopy trees may fluctuate to some extent within the region. Besides site conditions (soils, terrain orography), the differences seem attributable to a wide range of agents which may cause tree death in mixed-species montane forests (wind, snow, insects, fungal pathogens, etc.) and inherent variation in the spatial extent and intensity of disturbances (Webb 1989; Holzwarth et al. 2013). Inter- and intra-stand variation in species composition may also play a crucial role, as the vulnerability of fir, spruce and beech to particular disturbance agents differs considerably (Nagel and Diaci 2006; Schütz et al. 2006). For example, larger canopy openings resulting from severe windstorms or a clumped pattern of spruce mortality associated with insect outbreaks are conducive to the formation of patch structures (Korpeľ 1993; Sproull et al. 2016; Parobeková et al. 2018; Orman and Dobrowolska 2017). In particular, climate warming and associated adaptive changes in species compositions of natural ecosystems acting in the direction of a decreasing proportion of conifer species, especially spruce, may shift disturbance regimes towards more severe ones (Seidl et al. 2017). In contrast, single standing dead individuals killed by insects or fungal pathogens create small gaps, favouring a single-tree replacement pattern and a close-to-random distribution of canopy trees. Nevertheless, mixed-species forests are subjected to a variety of interfering disturbing factors, and our study shows that classifying the spatial distribution patterns of canopy trees into one type is overly simplistic.

The attractive associations between under-canopy trees found in this study and also reported by other authors for similar forest ecosystems (Szwagrzyk and Czerwczak 1993) might be attributable to the spatially selective establishment of regeneration cohorts and/or their further recruitment to the understorey layer in specific habitats, e.g. in canopy gaps. In fact, the fine-scale aggregations of under-canopy trees in the Lom and Perućica reserves (≤ 20 m) corresponded with the distribution patterns of dead canopy trees in these stands, which, at a scale of 708 m2 (i.e. on the plots with a 15 m radius), showed a tendency to occur in small clumps. In the Janj reserve the positive spatial dependence found for under-canopy tree densities was more extensive (up to 80 m) and again corresponded with the coarse-grained patchy distribution of dead canopy trees (at scales below 40 m) but also live canopy trees (at scales above 80 m) in this stand. Nevertheless, an in-depth analysis (manuscript in preparation) reveals that the relationships between the local density and basal area of dead canopy trees in different decay stages and the density of advanced regeneration and under-canopy trees are very weak in the study stands (the Pearson’s correlation coefficients sporadically attain values above 0.25) and thus indicates that generation replacement is a highly stochastic process. Several authors have reported multiple releases followed by periods of suppression in the life history of single trees from the canopy layer (Splechtna et al. 2005; Zeibig et al. 2005; Nagel and Svoboda 2008; Firm et al. 2009; Trotsiuk et al. 2012).

At the single plot level, the structural differentiation expressed by the VARh index strongly varied, but the majority of plots represented multi-layered structures. Although we revealed instances of non-random distribution patterns of dead and live canopy trees and under-canopy trees, the structural differentiation was characterized by random spatial variation, i.e. neighbouring plots separated by 20 m were not more similar than subsets randomly drawn from the entire sample. This finding confirms the fine-scale structural heterogeneity of the study stands and suggests the dominance of the single-tree replacement pattern as reported in another studies (Parobeková et al. 2018; Keren et al. 2017). Some authors have developed methods of multivariate identification of structural units (stand patches) based on biomass accumulation, the proportion of live to dead trees, the presence of large trees and structural differentiation (Emborg et al. 2000; Podlaski 2008; Král et al. 2010b; Glatthorn et al. 2018). Nonetheless, because developmental trajectories of initially similar units may diverge and the boundaries between them may change, these units should not be regarded as integral to the dynamic processes of tree regeneration, growth and mortality (Král et al. 2018).

Inter-regional variation

Stand basal areas and growing stocks in the investigated Dinaric forests (56.8–65.9 m2 ha−1 and 903–1134 m3 ha−1, respectively) were almost twofold greater than those of the Western Carpathian forests (36.1–38.4 m2 ha−1 and 478–574 m3 ha−1, respectively) in the Baniska, U źródeł Solinki, Łabowiec, Oszast, and Żarnówka reserves investigated by using an analogous methodology (Paluch 2007; Paluch et al. 2015). With respect to data reported from natural mixed-species montane forests in Europe (Holeksa et al. 2009), the stands in the Janj, Lom and Perućica reserves undoubtedly attain one of the highest stand volumes, while the natural forests from the Western Carpathians attain one of the lowest. These differences were partly attributable to a higher proportion of conifers in the Dinaric forests, which in this region accounted for between 62 and 82% of the stand basal areas. In the Western Carpathians, the proportion of conifers ranged between 16 and 55% in the four studied stands and amounted to 67% in one stand (Oszast reserve). Moreover, the higher average heights reached by trees of comparable dbh (e.g., 3–7 m for dbh > 80 cm) clearly indicate better site conditions in the studied Dinaric forests compared to those in the Western Carpathians.

Despite substantial differences in biomass accumulation, the estimated basal areas of dead trees were similar in both regions and varied between 18.0 and 24.9 m2 ha−1 in the Dinaric Mts. and between 20.2 and 28.0 m2 ha−1 in the Western Carpathians. Considering the higher growing stocks and proportions of conifer species that are more resistant to decomposition (Korpeľ 1993; Lombardi et al. 2008; Müller-Using and Bartsch 2009; Přívětivý et al. 2018) in the Dinaric forests, these similar quantities may result from the quicker decay of woody debris under a wetter climate (Zhou et al. 2007; Rock et al. 2008; Přívětivý et al. 2016) or indicate enhanced tree mortality in recent decades in the Carpathians. Indeed, numerous studies have reported significant changes in species composition of unmanaged forests in this region due to a declining proportion of fir and spruce (Vrška et al. 2009; Saniga et al. 2011; Jaworski 2016).

The Dinaric and Carpathian forests showed considerable similarity in the frequency distributions of BAL. At the single plot scale (154 m2), the variation of BAL was smaller in the Dinaric Mts. than in the Western Carpathians, but for areas above 1000 m2 these differences quickly disappeared (Figs. 4, 6). In both cases the distributions had a truncated bell-shaped form with a low proportion of the smallest BAL values. This finding confirms earlier suggestions based on material from the Carpathians that in montane mixed-species forests, natural disturbance rarely reduces the stand basal area to a close-to-zero level even in small stand patches and that these biomass losses are rapidly balanced by the growth of neighbouring or under-canopy trees. Considering the greater overall biomass accumulation in Dinaric forests, a higher proportion of stand patches with extremely high BAL values is not surprising. However, at the same time, Dinaric reserves are characterized by a considerably lower proportion of stand patches with extremely low BAL values, which may suggest a generally less severe disturbance regime or better resilience potential in this region. The latter characteristic is directly associated with the small-scale structural heterogeneity of a stand, as multi-layered stand patches with a portion of the trees unaffected by the disturbance can rapidly balance biomass losses (Jaworski and Paluch 2002; Jaworski 2016). Indeed, the proportion of stand patches with lower VARh values (< 0.4) representing less diversified structures tended to be, at small spatial scales (< 1000 m2), higher in the Carpathians than in the Dinaric Mts. (Fig. 5).

A noteworthy result of our study is the considerable variation in VARh distributions observable at the inter-stand level, in particular in the Carpathian stands (Fig. 5). It seems that the most significant source of this variation is related to differences in the density of small-diameter trees. In the stands from the Dinaric Mts. and Western Carpathians considered in our study, the proportion of trees of dbh < 30 cm ranged between 57 and 84%. These dissimilarities may significantly affect the values of the VARh index because less steep diameter distributions representing a reversed J-shaped type generally result in a higher proportion of stand patches with VARh index values close to or higher than 1, which is characteristic of a uniform or bimodal type of vertical differentiation. It seems probable that the density of small-diameter trees varying at the regional level in similar site conditions echoes past large-scale disturbances of moderate intensity which triggered regeneration pulses and resulted in an increasing frequency of post-disturbance cohorts of young trees (Nagel et al. 2006; Zenner et al. 2015). In fact, disturbances of moderate intensity in larger areas of close-to-natural stands of similar species composition have been recorded in several studies (Splechtna et al. 2005; Firm et al. 2009; Diaci et al. 2011; Trotsiuk et al. 2012; Nagel et al. 2014; Jaloviar et al. 2017). Depending on the time perspective, however, the effect of such disturbances involving at least a short-term reduction in biomass accumulation may be indiscernible at present because of the increased growth reaction of the trees unaffected by the disturbance.

Conclusions

The study, which was conducted by using an identical methodology in several mixed-species montane forest relics of primeval character in the Dinaric Mts. and Western Carpathians, indicated substantial inter-regional differences in the growing stocks of live trees and proportion of live to dead trees. The lower proportion of stand patches of simple vertical structure and the lower proportion of live to dead tree basal areas in the Dinaric Mts. may suggest a less severe disturbance regime in this region compared to the Western Carpathians. Nevertheless, in both regions the basal area distributions have a similar bell-shaped form with a small proportion of extreme values, indicating that natural disturbance rarely reduces local basal area to a close-to-zero level and that these biomass losses are rapidly balanced by the growth of neighbouring or under-canopy trees. The dominant proportion of complex structures at the scale of very small stand patches indicates the prevalence of a single-tree replacement pattern. This seems remarkable, especially in the case of Dinaric forests distinguished from those in the Western Carpathians by relatively high growing stocks. Nonetheless, variation in the densities of under-canopy trees and frequency distributions of the differentiation index observed at the between-stand level highlight the stochastic nature of the generation replacement process and its possible connection with recruitment pulses caused by occasional intermediate disturbances.