Composition of live, dead and downed trees in Järvselja old-growth forest

The study area is in Järvselja Training and Experimental Forest Centre and is a part of the Järvselja nature conservation area. In 1924, an area of 12.8 ha (Raukas, 1967) in the southern part of a 19-hectare compartment was set aside to allow an experiment where all management actions were abandoned (Mathiesen, 1940). In 1936, the area was enlarged to cover the whole compartment and the site was included in the nature conservation registry, obtaining a state-controlled protection status. In 1959, the whole old-growth forest compartment was registered as a botanical-zoological conservation area (Krigul, 1971). Today the old-growth forest and its surrounding compartments are part of the Järvselja nature conservation area that covers altogether 188 ha. According to the site’s conservation plan, this area is divided into two special management zones (Ürgmets and Järvselja) and one, partially limited, management zone (Apna) (Laas et al., 2007).

The protected old-growth forest (known as compartment JS226) (Kollom, 2011; Laas et al., 2007) is the most important, scientifically interesting and most studied area in the conservation area. The old-growth forest compartment is divided into 15 sub-compartments, altogether covering 19.4 ha according to the forest inventory in 2016. There are three different Natura 2000 habitat types apparent in this compartment: 91D0 (bog woodlands), 9020 (Fennoscandian hemiboreal natural old broadleaved deciduous forests) and 9050 (Fennoscandian herb-rich forests with Picea abies (L.) H. Karst.). Several registered, protected habitats are registered in the compartment with records for the following species: Epipactis helleborine (L.) Crantz, Skeletocutis odora (Peck ex Sacc.) Ginns, Cypripedium calceolus L., Cinna latifolia (Trevir.) Griseb., Corallorhiza trifida Châtel., Corydalis intermedia (L.) Mérat, Amylocystis lapponica (Romell) Singer, Leucopaxillus salmonifolius M.M. Moser & Lamoure, and Tricholoma apium Jul. Schäff.

Compartment JS226 (58.28 N and 27.32 E) is divided according to three dominating vegetation types (Lõhmus, 1993): on the eastern side, fresh boreo-nemoral forests with rich groundcover species like Oxalis acetosella L., Hepatica nobilis Schreb., Rubus saxatilis L. and Maianthemum bifolium (L.) F.W. Schmidt. The western side is dominated by drained swamp forests, dominated by Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies) in the tree layer. The central part of the compartment is a fen, with a constantly high ground water table and very lush plant cover. Dominant tree species are black alder (Alnus glutinosa (L.) Gaertn.) and silver birch (Betula pendula Roth).

One of the earliest studies about the old-growth forest ground vegetation was carried out by Gustav Vilbaste in 1929. In that study three dominating forest site types were defined and described in the old-growth forest compartment: fresh boreo-nemoral forests, minerotrophic mobile water swamp forest and mixotrophic bog forests (Sepp & Rooma, 1993; Vilbaste, 1929).

There are numerous subsequent descriptive studies about the old-growth forest compartment. Mathiesen published an overview about the compartment in 1940 and Masing wrote about the old-growth forest conditions in 1960 (Masing & Rebane, 1995). Krigul carried out several studies in Järvselja old-growth forest from 1939 to 1971 (see Table 1). While earlier studies focused on volume description of living trees and coarse woody debris (CWD) (Krigul, 1940), later studies focused more extensively on stand health conditions and also presented additional stand-wide inventory data (Kollom, 2011).

Irdt and Rebane conducted one of the most important research studies in the period of 1954–1985 (Irdt & Rebane, 1985). They established a complex measurement plot transect in the old-growth forest and remeasured it regularly from 1954 until 1985. The measurements on this transect were carried out on 10 × 10 m plots. They assessed standing live tree diameters, heights, tree layer, geographical location, ground vegetation and soil microsites. Based on those descriptions, Irdt & Rebane (1985) partly described the stand dynamics. Heino Kasesalu published several studies about the old-growth forest, also describing the stand dynamics since 1984 up to 2001 (Kasesalu, 1993; Kasesalu, 2004).

In addition to the studies presented in Table 1, there are stand-wise forest inventory data available, covering a period of 94 years (from 1922 to 2016). Long-term dynamics based on the stand-wise forest management inventory data were analysed by Kollom (2011). In her study, the standing gross volume of all tree species in old-growth forest compartments increased up to 1962, and was then followed by a period of decrease up to 1983. Similar trends were reported by Pettai (1985) and Kasesalu (2004). In the period of 1983–2010, the average volume of the dominating layer showed low fluctuations, but a clear increase in the volume of the supressed spruce layer can be seen starting from 1993 (Kollom, 2011).

From spring until autumn 2013, an inventory was carried out with the aim of covering the study area with tree mapping data. To avoid potentially disruptive extensive measurements on more fragile sites and close to habitats under protection, two methodologies were followed: whole surface inventory (in sub-compartments 6, 9, 10, 11, 12, 13, 14, 15) and partial surface mapping based on permanent sample plots (in sub-compartments 1, 2, 3, 4, 5, 7, 8). To facilitate mapping, additional measurements on tree positions on the site were carried out to establish a network of geodetically accurate reference points (approx. 50 × 50 m interval, RTK GNSS method, with the location accuracy of 0.05–0.5 meters). This resulted in a 72-point grid and the grid points were marked on the ground with metallic rods (see Appendix A). Whole surface mapping covered 9.34 ha and partial mapping covered 10.03 ha. The study area included five different site types: Oxalis drained swamp site type 1.52 ha, Alder fen site type 2.31 ha, Myrtillus drained swamp site type 6.07 ha, Aegopodium site type 6.62 ha and Transitional bog forest site type 2.84 ha (forest site types described following Lõhmus (2004).

FieldMap technology (IFER-Monitoring and Mapping Solutions Ltd.) was used in conjunction with geodetically accurate reference points (see Appendix A), where tree distance was measured with a ForestPro impulse laser rangefinder. For tree positioning the Mapstar TruAngle angle encoder together with ForestPro impulse laser rangefinder was used. Tree heights to the live crown base and heights of the lowest dry branch on the tree stem were measured with a Vertex hypsometer to the nearest 0.1 m. Data were collected and stored on site with the FieldMap data collector software. Partial surface stem mapping was carried out on wetter and more fragile sites in parallel with full surface stem mapping elsewhere. Circular plots of various radii (15–20 m) were used in different sub-compartments. Permanent plots were located similarly with FieldMap measurements using the geodetically accurate reference points as shown in Appendix A.

Both full surface and partial surface tree mapping followed the design of the survey protocol of the Estonian Network of Forest Research Plots (ENFRP) (Kiviste et al., 2015). The measurement protocol was amended so that all tree heights were measured, where possible (in ENFRP a sub-sample of trees, e.g., every fifth tree, had tree height recorded). Live standing trees: all trees with DBH larger than 7 cm were measured. For every tree, measurements included tree position, tree species, tree height, live crown length, and diameter (in two perpendicular directions). The tree layer was divided according to tree classes into an upper storey (dominant layer and secondary cohort) and a lower storey (under-storey cohort and regeneration layer). For all trees, tree foliage vitality (vital, normal, poor) and dominance (dominant, co-dominant, suppressed) were determined. Standing and broken dead trees (snags: all standing dead or broken trees with DBH larger than 7 cm) were measured.

Other tree measurements included position, species, height, and diameter (in two perpendicular directions). If the height of the broken dead tree was under 1.3 m, then the diameter was measured 10 cm below the breakage point or from the root collar (marked accordingly). Standing and broken dead trees were also described according to five decay stages following studies by Fridman & Walheim (2000) and Köster et al. (2015): 1) Hard wood, recently dead; 2) Slightly decayed wood, outer layer of wood starts to soften; 3) Medium decayed, outer layers of stem wood soft, core still hard; 4) Well-decayed wood, wood soft through the tree stem; 5) Almost decomposed, wood soft, fragmented. Logs: All downed dead wood and tree stem parts with a top diameter larger than 6 cm were measured. Other measurements for every tree or stem included: top and trunk position, tree species, stem length or the length of a stem part, top and trunk diameter (in two perpendicular directions). The five decay stages of downed stems were also described according to the classification as described above for standing dead trees. The spatial arrangement of measured standing, dead and downed trees in Järvselja old-growth forest is shown in Figure 2.

Assessing the tree height and diameter (H-D) relation of old-growth forests using models that have been developed for managed forest stands has its own challenges (Nigul et al., 2021). Here, we employed several models and developed a model selection procedure, mostly to tackle problems like a very small number of measured trees per species and structural cohort, linearization and bias removal, or mean tree height dependent scaling for very small or very large trees (Figure 1). Our basic tool is the Näslund function with a controlled starting point at breast height (initiated 1.3, 0) (Näslund, 1936; Schmidt et al., 2011): (1) h=1.3+(da+b⋅d)3, h = 1.3 + {\left( {{d \over {a + b \cdot d}}} \right)^3}, where h refers to tree height, d is diameter at breast height (DBH), a and b are model parameters determining the location of the curve’s inflection point and the convergence behavior towards the height asymptote. Transforming Equation 1 by linearization yields: (2) dh−1.33=a+b⋅d. {d \over {\root 3 \of {h - 1.3} }} = a + b \cdot d. To avoid systematic bias in the parameter estimation, introduced by the transformation to the linear form, we applied adjusted weight functions as suggested by Artur Nilson (Nilson, 2002). We applied the following weight function: (3) W=(h−1.3)w, W = {\left( {h - 1.3} \right)^w}, where W = weight and h = tree height. w is a searchable parameter and can be chosen so that the regression residuals are minimized, and either are near or equal to 0. This criterion is achieved when w lies in the interval between 0.6 and 0.75.

Figure 1

We introduced an additional parameter c to obtain a scalable version of Equation 1: (4) h=1.3+c⋅(da+b⋅d)3, {\rm{h}} = 1.3 + {\rm{c}} \cdot {\left( {{{\rm{d}} \over {{\rm{a}} + {\rm{b}} \cdot {\rm{d}}}}} \right)^3}, where the variables and parameters are as given as in Equation 1. The new parameter c is either set to 1 in the case of a sample size n of measured trees is ≥ 9 (Figure 1 left side) or it is calculated from raw data by fixing the parameters a and b, using the default parameters from Appendix B, in the case that n < 9 (Figure 1 right side).

The selection procedure in the case of n ≥ 9 is as follows: from the weight function W, the scaling procedure creates a polynomial form cc. If cc = 0, the asymptotic height is set to 55 m and scaled parameters a and b are calculated including the weight function W. If cc ≠ 0, new parameters a and b are determined as given in Figure 1 and in the case of a < 0.08 new scaled parameters a and b are calculated where b depends on the mean height only. In all other cases, the calculated parameters a and b are then applied to Equation 4. The whole selection and scaling process is iterated until the sum of residuals of the regression is minimized and approaches 0. In short, the procedure either calculates the scaling parameter c or it scales the parameters a and b with respect to the weight function W.

Following the model selection procedure, H-D ratios were assessed, and model solutions calculated by using the measured data per stand element (Figure 1) – a combination of the tree species in a particular cohort in a particular site type and sub-compartment. There were 5 site types and 15 sub-compartments in the Järvselja old-growth forest compartment. Similar calculations were carried out for deciduous and coniferous trees, which resulted in 129 different height curves.

For comparison, we used the two best-performing height curve models from an earlier study by Nigul and colleagues (Nigul et al., 2021): a)

Chapman-Richards function with fitted parameters per species (Table 2): (5) h=1.3+α⋅(1−e−βd)γ, h = 1.3 + \alpha \cdot {\left( {1 - {e^{ - \beta d}}} \right)^\gamma }, where h= tree height; d= DBH; α, β and γ are constants per tree species;

To achieve comparable results, a scaling factor x, like c in Equation 4, was introduced: (6) h=1.3+x⋅α⋅(1−e−βd)γ, where x=∑(h−1.3)∑[α⋅(1−e−βd)γ], h = 1.3 + x \cdot \alpha \cdot {\left( {1 - {e^{ - \beta d}}} \right)^\gamma },\,{\rm{where}}\,x = {{\sum {\left( {h - 1.3} \right)} } \over {\sum {\left[ {\alpha \cdot {{\left( {1 - {e^{ - \beta d}}} \right)}^\gamma }} \right]} }}, (7) h=1.3+x⋅(da+b⋅d)3, where x=∑(h−1.3)∑[(da+b⋅d)3], h = 1.3 + x \cdot {\left( {{d \over {a + b \cdot d}}} \right)^3},\,{\rm{where}}\,x = {{\sum {\left( {h - 1.3} \right)} } \over {\sum {\left[ {{{\left( {{d \over {a + b \cdot d}}} \right)}^3}} \right]} }}, where h is tree height; d is DBH; α, β, γ, a, and b are constants for tree species (Table 4).

All three H-D curves were used to estimate single tree heights and to make comparisons with measured tree heights. Root mean square errors were calculated to assess the goodness of fit (Figure 3).

General information about measured live and dead standing trees is presented in Table 5, and the spatial allocation can be followed in Figure 2 (B). The measured tree height was used for the modelling of the tree stem volume (m³) and the tree stem lateral surface area (SLSA, m²). Missing tree heights were calculated using Equation 4. The calculation principle is described in more detail by Padari (2020). For species-specific tree functions a modification of Ozolinš taper curve formula (Ozolinš, 2002; Silava, 1988) was used for calculating diameter at any tree height: (8) γ(x)=1+(x2−0.01)⋅(p⋅(h−h0)+q⋅(d1.3−d0)), \gamma \left( x \right) = 1 + \left( {{x^2} - 0.01} \right) \cdot \left( {p \cdot \left( {h - {h_0}} \right) + q \cdot \left( {{d_{1.3}} - {d_0}} \right)} \right), (9) dl=d1.3⋅γ(lh)⋅(a0+a1⋅(lh)+a2⋅(lh)2+a3⋅(lh)3+a4⋅(lh)4+a5⋅(lh)5+a6⋅(lh)6)γ(1.3h)⋅(a0+a1⋅(1.3h)+a2⋅(1.3h)2+a3⋅(1.3h)3+a4⋅(1.3h)4+a5⋅(1.3h)5+a6⋅(1.3h)6), {d_l} = {d_{1.3}} \cdot {{\gamma \left( {{l \over h}} \right) \cdot \left( {{a_0} + {a_1} \cdot \left( {{l \over h}} \right) + {a_2} \cdot {{\left( {{l \over h}} \right)}^2} + {a_3} \cdot {{\left( {{l \over h}} \right)}^3} + {a_4} \cdot {{\left( {{l \over h}} \right)}^4} + {a_5} \cdot {{\left( {{l \over h}} \right)}^5} + {a_6} \cdot {{\left( {{l \over h}} \right)}^6}} \right)} \over {\gamma \left( {{{1.3} \over h}} \right) \cdot \left( {{a_0} + {a_1} \cdot \left( {{{1.3} \over h}} \right) + {a_2} \cdot {{\left( {{{1.3} \over h}} \right)}^2} + {a_3} \cdot {{\left( {{{1.3} \over h}} \right)}^3} + {a_4} \cdot {{\left( {{{1.3} \over h}} \right)}^4} + {a_5} \cdot {{\left( {{{1.3} \over h}} \right)}^5} + {a_6} \cdot {{\left( {{{1.3} \over h}} \right)}^6}} \right)}}, where γ(x) = perturbation coefficient, dl = diameter at the distance from stem base l, cm; p, h₀, q, and d₀ are parameters for the perturbation coefficient; d_1.3 = breast height diameter, cm; l = distance from stem base, m; h = tree height, m; a₀, a₁, …, a₆ are coefficients of the taper curve.

Figure 2

A volumetric transformation is needed for the tree stem volume calculation based on the taper curve. In the current study it was calculated using the following formula: (10) v=π40000∫0ldl2dl, v = {\pi \over {40000}}\int\limits_0^l {d_l^2dl,} where v = above bark stem volume up to the distance l from stem base (m³); l distance from stem base (for live tree equal to tree height) (m); dl =diameter at the distance l from stem base, cm (estimated by Equation 8 and 9).

SLSA above bark was calculated as follows: (11) S=π100∫0ldl dl, S = {\pi \over {100}}\int\limits_0^l {{d_l}\,dl,} where S = above bark SLSA up to the distance l from stem base 0, m²; l = distance from stem base (for live and standing dead trees) the tree height h, m; dl = diameter at the distance l from stem base l, cm (Equation 8 and 9).

The measured standing snags (see Table 5 for general statistics and Figure 2 (B) for spatial arrangement) were divided into two groups: snags that were at least 1.3 m or taller in height and snags lower in height than 1.3 m (stratification needed due to different calculation methods). For taller snags h= or ≥1.3 m, breast height diameter was measured. For shorter snags, the diameter was measured below the breakage point and needed to have DBH calculated. The breast height diameter was calculated by fitting the species-specific stem taper curve (Equation 8 and 9) to the measured diameter at snag height. For all measured snags the modelled tree height was further calculated using Equation 4 and applying the selection procedure as described before (see also Figure 1) to obtain the species-specific H-D curve. For species-specific tree taper functions in Estonian conditions, a modified Ozolinš taper curve (Equation 8 and 9) was used for calculating the diameters in the measured snag height above the bark. The calculation principle is described in more detail by Padari (2020). Snag volume (m³) was calculated based on the snag breast height diameter d_1.3, predicted height h and snag height l by Equation 10. Stem section lateral surface area above bark was calculated by Equation 11.

For CWD on the site, both co-ordinates and diameters of the bottom and top sides of tree stem sections were measured (see Table 5 for general statistics and Figure 2 (A) for spatial placement). For spatial assessment, a tree stem section was cut if the section was placed outside a particular plot or sub-compartment using GIS analysis. The diameters were interpolated for the cut point and only the sections inside a particular plot or sub-compartment were used for volume calculation. For stem section volume calculations, a truncated cone formula was used: (12) v=π⋅l⋅(d12+d1⋅d2+d22)12000, v = {{\pi \cdot l \cdot \left( {d_1^2 + {d_1} \cdot {d_2} + d_2^2} \right)} \over {12000}}, where v = stem section volume (m³), l = stem length (m), d₁ and d₂ = diameters for the stem, cm; accordingly, d₁ = diameter at the base of the stem section and d₂ = diameter at the top of the stem section.

The stem sections’ lateral surface area was calculated using the stem lateral surface area of a truncated cone as follows: (13) S=π⋅l2+(d1−d2)240000⋅d1+d2200, S = \pi \cdot \sqrt {{l^2} + {{{{\left( {{d_1} - {d_2}} \right)}^2}} \over {40000}}} \cdot {{{d_1} + {d_2}} \over {200}}, where S = stem section lateral surface area from section base to section top, m²; l = distance from stem section base to the section top, m; d₁ = diameter at stem section base, cm and d₂ = diameter at stem section top, cm.

For the estimation of wood density (wood dry mass per m³) and calculation of the respective C content of fresh (growing trees) timber and timber in different decomposition stages (dead and downed trees), data from different studies were employed. For the estimation of fresh (growing trees) wood densities in different tree species most of the densities used are available from Saarman & Veibri (2006); for grey alder the fresh wood density was available from (Uri et al., (2014) (see Table 4). Fresh wood C relative content was available for pine (Uri et al., 2019), birch (Uri et al., 2017) and black alder (Uri et al., 2014). In the case of fresh wood C relative content, we applied the pine value for spruce, the birch value for aspen and the black alder value for grey alder.

For the assessment of C content in dead and downed stems along the different decay classes (fresh wood and five decay classes), we used available measures on wood density and the relative C content for six tree species. Those are Scots pine, Norway spruce, silver and downy birch, black alder, common alder and common aspen measured by Köster et al. (2015). We then used these measurements (Köster et al., 2015) to derive the wood density dynamics as depending on the wood decay classes. The results of the regression analysis with parameter estimates are presented in Table 3.

Wood density along different decay classes was estimated as follows: (14) dDecay=d0+a1⋅d0⋅Decaya0, {d_{Decay}} = {d_0} + {{{a_1} \cdot {d_0} \cdot Decay} \over {{a_0}}}, where d_Decay = basic density in decay class, g/m³; d₀ – fresh wood basic density, g/m³ (Table 3); Decay – decay class (1 to 5); a₀, a₁ – model constants (Table 4).

For the estimation of the C content for all tree species represented in Järvselja old-growth forest the following substitutions were used: for rowan, the birch values were used, for linden and unknown species the aspen values were used, for hazel, bird cherry and other broadleaved species grey alder values were used and for fir the spruce values were used. The wood density estimations are calculated using Equation 14 and presented in Table 4.

Stem wood C relative content was estimated as follows: (15) C=V⋅d⋅C%100, C = {{V \cdot d \cdot C\% } \over {100}}, where C is carbon mass, kg; V is tree/snag/laying tree volume, m³; d is wood density (dry weight/raw volume), kg/m³; C% is carbon relative content, %.