Temporal Evolution of Carbon Stocks, Fluxes and Carbon Balance in Pedunculate Oak Chronosequence under Close-To-Nature Forest Management

Abstract

Under current environmental changes, forest management is challenged to foster contrasting benefits from forests, such as continuous wood supply while preserving biomass production, biodiversity conservation, and contribution to climate change mitigation through atmospheric carbon sequestration. Although being found as globally important, estimates of long-term forest C balance are still highly uncertain. In this context, the chronosequence experiments (space-for-time substitution) might fill this gap in even-aged forests, as they represent an approach that enables the assessment of forest net C balance in the long term. In this research, we explored the dynamics of C stocks and fluxes in different forest pools throughout the rotation period (140 years) of a Pedunculate oak (Quercus robur L.) forest in Croatia. For this purpose, we selected a chronosequence that was made up of seven forest stands with different age (5, 13, 38, 53, 68, 108, and 138 years). To address the issues of uncertainty in C balance estimates, we compared net ecosystem carbon balance (NECB) estimated while using two different approaches, which we name pool-change (from C stocks) approach and component-flux (from C fluxes) approach. Overall, the pool-change approach showed higher NECB estimate, with the greatest difference being observed in younger stands (<50 years). Component-flux approach showed significantly higher uncertainty. Throughout the rotation period, managed pedunculate oak stands become a C sink early in their development phase, between the age of 13 and 35 years according to pool-change and component-flux approach, respectively. During the 140 years, oak forest provided 187.2 Mg C ha⁻¹ (604 m³ ha⁻¹) through thinnings and 147.9 Mg C ha⁻¹ (477 m³ ha⁻¹) in the final cut, while preserving, on average, 88.9 Mg C ha⁻¹ in mineral soil down to 40 cm, 18.2 Mg C ha⁻¹ in dead wood, and 6.0 Mg C ha⁻¹ in the forest floor. Soil C stocks in our chronosequence did not show any age-related trend, indicating that current management practice has no negative effect on soil C stocks. Finally, under current close-to-nature forest management, Pedunculate oak forest showed to be sustainable in providing both economic and ecological ecosystem services.

1. Introduction

Forests provide many direct and indirect benefits for human well-being. They represent the most important terrestrial carbon (C) pool [1], which plays a crucial role in global C cycle and climate regulation [2], and they facilitate human demands for wood supply through biomass production, while providing a variety of other ecosystem services [3]. The future of this terrestrial sink, especially in temperate forests, is highly uncertain due to the strong forecasted anthropogenic and environmental pressures [4]. Therefore, investigating long-term forest C dynamics, especially in the temperate region and under various management practices, has been gaining great attention in the last decades [5,6,7,8,9,10,11,12,13,14].

Forest C dynamic is commonly investigated by assessing C stocks (i.e., C stored in various forest pools; in Mg C ha⁻¹) and/or C fluxes (i.e., C stock change within different pools; in Mg C ha⁻¹ y⁻¹). During time, both C stocks and fluxes change [7,15] due to stand age, environmental conditions [16,17,18], and/or disturbances [19,20], and they may result in either C accumulation or loss.

Net ecosystem productivity (NEP; Mg C ha⁻¹ y⁻¹) is considered as a common measure of yearly ecosystem carbon balance and it is calculated either from component-flux measurements or from Eddy covariance Net Ecosystem Exchange (NEE) flux [21]. However, in managed forest ecosystems, where harvests play a crucial role in regulating C dynamics [7,22,23], it is necessary to account for cumulative harvest C losses when calculating long-term net ecosystem carbon balance (NECB) [24,25].

Long-term NECB can be estimated while using two different approaches: (1) pool-change approach, where C balance is calculated as the cumulative difference between C stocks at different stand ages, or (2) component-flux approach, where the C balance is calculated as the difference between the cumulative sum of fluxes in different ecosystem components and cumulative harvest C losses. It should be highlighted that, for short-term (e.g., one year) estimates of NECB, special attention must be given to the changes in the carbon stocks of non-structural carbohydrates (NSC) of trees [26]. However, over a longer period (e.g., decades), the contribution from changes in NSC becomes less important, as annual NSC changes tend to cancel each other out, and the estimates on the NECB from the above two approaches should converge [26]. NECB calculated from pool-change approach reflects long-term dynamics that are mainly related to the past environmental conditions and forest management, while NECB that was calculated from the component-flux approach is related to the current environmental conditions (i.e., meteorology) and management. Under stable environmental conditions and unchanged management, both approaches should yield similar results. However, there are numerous evidences that growing conditions have been changing in the last few decades [4,27] and there are many environmental drivers that can affect forest C fluxes, both in positive (e.g., CO₂ fertilization, nitrogen deposition, extension of the growing season) and negative manner (drought, increased respiration due to higher temperatures, new pests, and diseases) [28]. Consequently, C balance estimates from stocks and fluxes can differ and observed differences can be variable across different stand ages.

The chronosequence method (space-for-time substitution) is widely used for studying age-related processes in forest ecosystems, e.g., forest growth [16], development of C stocks [6,7,11,29], and modelling ecosystem productivity [12]. This method applies a false time series of individual forest stands with similar structural, soil, and climate characteristics, with a basic assumption that the only difference among the stands is due to their age [30]. Even though there are many shortcomings that are related to this method [30,31], if appropriately applied, the chronosequence approach offers invaluable insights into the temporal dynamics of forest ecosystem processes.

In temperate region, oak forests are widely recognized as very productive ecosystems, which are valuable in both the economic and ecological sense, but at the same time are very vulnerable due to the great pressure in terms of timber exploitation [32] and serious oak dieback [33]. The first aim of our study is to describe the temporal evolution of carbon stocks by pools and carbon fluxes in the pedunculate oak forest, managed with the close-to-nature forest management, during a 140 years long rotation. The second aim is to investigate the uncertainty of estimates of long-term NECB by adopting two different estimate approaches, named here, pool-change approach, and component-flux approach. For that purpose, C stocks and fluxes are measured in pedunculate oak stands of different age, net ecosystem productivity is modelled as a function of age, and, taking into account cumulative carbon losses due to harvest removals, the estimates of NECB from two approaches are compared.

2. Materials and Methods

2.1. Study Area and Forest Management Practice

The study area is located in the continental part of Croatia, 35 km southwest from Zagreb, within the managed forests of river Kupa basin (Figure 1). The climate in the area is warm temperate with a mean annual temperature of 10.6 °C and precipitation of 962 mm y⁻¹ for the period 1981–2010 (Appendix A, Figure A1 and Figure A2). Relief is mostly flat with an altitude that ranges from 107 to 130 m a.s.l. Soil texture is loam and clay loam to clay. Soils are hydromorphic, hypogley, eugley, and pseudogley, on clay parent material. Groundwater table level during vegetation period typically varies from −60 to −200 cm [34]. The study site is an even-aged managed forest dominated by Pedunculate Oak (Quercus robur L.), with a significant share of Common hornbeam (Carpinus betulus L.), Black alder (Alnus glutinosa (L.) Geartn.), and Narrow-leaved ash (Fraxinus angustifolia Vahl.) in the tree layer. The understory layer is dominated by hazel (Corylus avellana L.) and common hawthorn (Crataegus monogyna Jacq.) [35].

In Croatia, Pedunculate oak forests cover 225 kha (9% of the total forest area) and they have a growing stock of 66.5 Mm³ (16% of the total growing stock) [36]. They are managed according to the 140 year-long rotations that may be extended if a stand is in a good condition and rotation extension would contribute at attaining uniform distribution of areas with respect to age-classes. Thinnings are conducted every ten years, ending with two or three regeneration cuts in a shelterwood-system during the last 10 years of the rotation. The first cut, termed Stand preparation cut, is aimed at enhancing acorn production of the standing oak trees. In the year of good acorn mast follows the Establishment cut, which aimed at promoting germination of acorns and establishment of young seedlings under the shelter of remaining mother-trees. Finally, 2–3 years after the establishment cut, the Final cut is made and remaining old trees are harvested, which resulted in a regenerated pedunculate oak stand. Furthermore, in the case of poor acorn mast, the acorns from designated seed areas are introduced into the stand before the final cut. On areas of the stand where seedling density remained too low after the final cut, planting of pedunculate oak seedlings of known origin (stand in-filling) must be made within three years after the final cut [37]. Furthermore, a given number of old trees (typically 1 per ha) must be left unharvested for biodiversity. In this way, a natural regeneration under continuous tree crown cover is ensured. After harvest, only coarse wood with a diameter equal to or greater than 7 cm is typically exploited, while all other woody debris and leaf/fruit litter are left at the site. Recently, however, due to demand for wood for energy, the thicker branches are also removed at some sites, while the thinnest one (<2 cm), as well as stumps and litter, are still left in the forest. This approach has been considered sustainable in terms of yield, biodiversity, and ecosystem stability [38], and in Croatia it dates back, more-or-less unchanged, to 1769, when the first forest management plan was made [39]. Regarding present forest management definitions [40], such an approach is classified as close-to-nature forest management.

2.2. Experimental Design and Field Measurements

The experiment was set up at the beginning of 2010. The chronosequence consisted of seven forest stands with different age (5, 13, 38, 53, 68, 108, and 138 years). In each stand, four permanent circular plots with a variable radius (Appendix A,

Table A1. Plot radius (m), dependent on the stand age and diameter at breast height of the sampled tree.

DBH (cm)	Age Class *
	I. (1–20 y)	II. (21–40 y)	III. (41–60 y)	IV. (61–80 y)	VI. (101–120)	VII. (121–140)
2–5	2	2	2	2	2	2
6–10	3.5	3.5	3.5	3.5	3.5	3.5
11–30	7	8	10	12	16	17
31–50	7	8	13	13	16	17
51–80	7	8	13	13	16	17
>80	56.4 (1 ha)	56.4	56.4	56.4	56.4	56.4

* There was no stand in the V. age class (80–100 y) within the Jastrebarsko forest that would meet the criteria required to be included in the chronosequence. (Note on the reading the table: e.g., for the III age class, all trees with DBH of 2–5 cm were measured within 2 m radius from the plot centre, and within 3.5 m radius all trees with DBH 6–10 cm, and within 10 m radius all trees with DBH of 11–30 cm were measured, etc.).

) were set up within a one-hectare area according to Terrestrial Carbon Observing (TCO) Protocols [41], with the exception of the youngest stand (five years), where five 3 × 3 m plots were set up. Permanent plots represent the base of chronosequence experiment, but they also facilitate other field experiments that are beyond of the scope of this study. In the 38 years old stand, a permanent eddy-covariance system was present [42]. Furthermore, to account for the overall spatial variability of the stand, again with the exception of the youngest stand, additional six to 10 temporary circular plots, with variable radii, were set in a 100 × 100 m grid around each set of permanent plots. The number of additional plots was limited by the size of the corresponding forest compartment (Figure 1).

Table 1. The main characteristics of forest stands in the chronosequence in 2011 (Mean ± Std. error).

Stand Age (y)	DBH (cm) b	H (m) b	N (Tree ha−1)	BA (m2 ha−1)	V (m3 ha−1)	Clay (%)	Silt (%)	Ts (°C)	SWC (m3 m−3)	FC (m3 m−3)
5 a	[0.91]	[0.53]	[120,000 ± 12,217]	[7.96 ± 0.83]	[1.4 ± 0.1]	0.22	0.48	11.05	0.392	0.52
13	4.7 ± 1.9	5.4 ± 0.8	8507 ± 1029	11.63 ± 1.42	36.7 ± 4.7	0.28	0.49	10.45	0.379	0.58
38	16.0 ± 5.4	17.6 ± 3.4	1684 ± 188	27.03 ± 1.49	267.0 ± 15.3	0.51	0.37	10.62	0.403	0.53
53	23.1 ± 7.6	21.2 ± 3.9	860 ± 215	26.24 ± 1.5	310.4 ± 13	0.33	0.44	10.31	0.394	0.54
68	28.2 ± 10.6	23.5 ± 5.2	406 ± 39	26.62 ± 1.53	376.5 ± 22.5	0.25	0.49	n.a.	n.a.	n.a.
108	37.5 ± 16.5	26.0 ± 6.3	390 ± 144	30.7 ± 2.79	492.3 ± 43.8	0.25	0.49	10.35	0.383	0.49
138	42.9 ± 23.2	24.3 ± 7.4	169 ± 22	28.57 ± 1.26	466.5 ± 20.1	0.26	0.50	10.63	0.437	0.58

DBH—mean diameter at breast height, H—mean tree height, N—number of trees per ha, BA—stand basal area, V—stand volume, Ts—yearly average soil temperature at 5 cm depth, SWC—yearly average soil water content at 30 cm depth, FC—SWC at field capacity, ^a Square brackets indicate that the root collar diameter and not DBH was measured; all saplings were counted. ^b Mean ± Standard deviation.

gives the main characteristics of each forest stand.

2.2.1. Aboveground Live Biomass Stocks

On each permanent and temporary plot, diameter at breast height (DBH) and tree height were measured for all live trees with DBH > 2 cm. In the youngest stand (five years old), all of the saplings per plot were counted, and on a subsample of 30 saplings per plot, a root collar diameter and a total sapling height were measured. The measurements were taken once a year, at the end of 2010 and 2011.

2.2.2. Dead Wood, Litterfall and Soil

Dead wood includes snags, stumps, and lying dead wood. Snags and stumps diameters and heights were measured, together with live trees. At the centre of each temporary plot, two 10 m crossed intersects were placed and all lying dead wood with a diameter greater than 1 cm was recorded according to the line intersect method [43]. For all dead wood, the decomposition class was estimated according to Hunter (1990) [44]. In the five years old stand, all lying dead wood was piled up after harvesting and deposited on the top of the stumps, therefore the line intersect method could not have been applied, as it assumes a random distribution of lying dead wood. Thus, we measured all stumps, counted all piles within a 50 × 50 m area, and based on the weight of one average pile, estimated that the average biomass of one pile was 20 kg. Dead wood measurements were taken once in 2011.

One litter trap was placed in the centre of each permanent plot and litterfall was collected several times a year both in 2010 and 2011. The samples from litter traps were transported to the laboratory for further processing on the day of collection.

Soil was sampled at four positions within each permanent plot, 5 m away from the plot centre in the directions N, E, S, and W. The forest organic soil layer (forest floor) was sampled with metal collar (d = 16 cm) and separated into two fractions: litter layer (OL) and fragmented and humus layer (OFH), according to the The International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects on Forests (ICP Forests) manual [45]. Mineral soil was sampled with split tube sampler (Eijkelkamp, Giesbeek, The Netherlands) and each soil core was cut into four samples according to a predefined depths (0–5, 5–10, 10–20, and 20–40 cm). The soil was sampled once in 2011. In total, 16 soil core and 16 forest floor samples were taken per stand, representing one age class. Samples of forest organic and mineral soil layers were further processed in the laboratory.

2.2.3. Meteorological Measurements, Soil Respiration and Dead Wood Decomposition Flux

Meteorological and soil respiration measurements were performed in six out of the seven stands (5, 13, 38, 53, 108, and 138 years). A small, custom made, meteorological station was set up inside 2 × 3 m fenced area within each stand. Since December 2010, stations recorded the temperature of air (Tair) and soil (Tsoil) at 5, 10, and 20 cm depths with thermometer sensors DS18B20+ (Maxim Integrated, San Jose, CA, USA) and volumetric soil water content (SWC) at 0–10 cm depth with ECH2O Soil Moisture Sensor 10HS (Decagon Devices, Pullman, WA, USA). At the eddy covariance site (38 years old stand), meteorological conditions were measured by the tower’s weather station [35,42,46]. The period between two successive measurements at all weather stations was 30 min.

Soil CO₂ efflux was measured while using an automatic chamber system previously described by Delle Vedove et al. (2007) [47]. In 2011, the soil respiration system (with two to three chambers) was continuously moved among the selected stands every 3–4 days. The measurements were always conducted at the same location within each stand. Inside the previously mentioned fenced area, metal collars, on top of which the flux chamber is fitted, were placed in autumn of 2010 and were left in the ground permanently in order to minimize soil disturbance during 2011 measurement season. When at the location, the system measured the soil CO₂ efflux every hour. At eddy covariance site soil respiration system consisted of three to five chambers that were permanently placed into the ground and measurements were performed cyclically, every four hours throughout the year. Data on soil temperature at 5 cm depth (Tsoil5) and soil water content (SWC) were used for soil respiration gap-filling (see chapter 2.5). SWC at field capacity (FC) was estimated to be equal to daily average SWC recorded in spring, before the bud-burst, three days after the SWC started dropping from saturation, similar to FC that was mentioned in the work of Reichstein et al. [48]. The fact that the soil is saturated was observed at the sites (footprints filled-in with water after walking) and corroborated with SWC records (SWC does not increase after the rain events).

In three stands (5, 38, and 108 years), a decomposition experiment was set up with the purpose of estimating decomposition rates for different tree species and dead wood size classes. A detailed description of the decomposition experiment is given in Ostrogović et al. (2015) [49].

2.3. Laboratory Analysis

All the laboratory analyses were performed at the Croatian Forest Research Institute in Jastrebarsko, Croatia. Carbon content in samples from different C pools was measured while using CNS 2000 elemental analyser. Litterfall was separated in leaves and fruits, dried at 105 °C, weighted, and leaves were analysed for C content. The soil was sieved and fine roots (d ≤ 2 mm) separated, washed, dried at 105 °C and weighted. Forest floor samples were pooled to one sample per plot for each fraction (i.e., OL and OFH separately), oven-dried, weighed, and analysed for C content. Mineral soil samples were sieved, oven-dried, weighed, analysed for C content, and soil bulk density was estimated for each soil sample and depth while using a known volume of the sample.

2.4. Carbon Stocks Calculations

Carbon stocks were estimated for live biomass (CLB), dead wood (CDW), forest floor (CFF), and mineral soil layer down to a depth of 40 cm (CM40) and were expressed in Mg C ha⁻¹. CLB was calculated, as described in details below, for aboveground wood, leaf and fruit biomass, belowground fine root biomass, and coarse root biomass while using data on tree diameter and height and applying local allometric equations [50,51,52], basic wood densities [53], biomass equations for twigs and branches [54], measured leaf and fruit biomass, age-dependent coarse root-to-shoot ratios [55,56], measured live fine root biomass, and C content in wood [6], and leaves and fruits [56]. CLB was equal to:

$$CLB=\left(W_{AGW}+W_{CR}+W_{FR}\right)·C{F}_W+W_L·C{F}_L+W_{FT}·C{F}_{FT}$$

(1)

where W_AGW is aboveground wood biomass, W_CR is coarse root biomass, W_FR is fine root biomass, W_L is leaf biomass, W_FT is fruit biomass, CF_W is C fraction in wood, CF_L is C fraction in leaf, and CF_FT is C fraction in fruit.

Aboveground wood biomass was calculated, as follows:

$$W_{AGW}=V_{D >3}·\rho +W_{D \le 3}$$

(2)

where V_{D > 3} is the volume of stems and branches (d > 3 cm on thinner end), ρ is species-specific basic wood density, W_{D ≤ 3} is biomass of twigs, and branches with d ≤ 3 cm estimated while using a biomass equation factor [54].

Coarse root biomass was calculated, as follows:

$$W_{CR}=W_{AG}·RS$$

(3)

where W_AG is total aboveground biomass (wood, leaf, and fruit), RS is literature age-dependent root-to-shoot ratio [55,56].

CDW was calculated for above and belowground dead wood biomass. Aboveground biomass of snags, stumps, and lying dead wood (i.e., coarse and fine woody debris) was estimated from the volume, while using wood density for different stages of decomposition [6]. The volume of snags and stumps was estimated in the same way as for living trees, and the volume of lying dead wood was estimated while using the transect method [43]. The belowground biomass of snags was calculated in the same way as coarse roots in live biomass. For calculating the belowground biomass of stumps, first the removed aboveground biomass was estimated while using quadratic polynomial function between diameter at 0.3 m (d_{0.3 m}) and aboveground biomass, i.e., $W_{stump}=0.0004\cdot d_{0.3 m}{}^2-0.0047\cdot d_{0.3 m}$ (p < 0.01, R² = 0.9829), calibrated on the subsample of 94 trees. Finally, CDW was calculated with the equation:

$$CDW=\left(V_{DW}\cdot {\rho }_{DEC}+V_{snag}·{\rho }_{DEC}·RS+W_{stump}·RS\right)·C{F}_W$$

(4)

where V_DW is the total volume of dead wood (snags, stumps, and lying dead wood), ${\rho }_{DEC}$ is wood density for different stages of decomposition, V_snag is the volume of snags, and W_stump is the aboveground biomass of stumps, RS is literature age-dependent root-to-shoot ratio [55,56], and CF_W is C fraction in wood.

CFF was estimated from the forest floor samples, i.e., OL and OFH fractions, and it was calculated, as follows:

$$CFF=W_{OL}·C{F}_L+W_{OFH}·C{F}_{OFH}$$

(5)

where W_OL is the dry weight of litter layer, W_OFH is the dry weight of fragmented to humified layer, while CF_L and CF_OFH are carbon fraction in leaves and in fragmented to humified layer, respectively.

CM40 was estimated from mineral soil cores. Each soil core consisted of four samples for different depths, and carbon stock in mineral soil down to 40 cm was calculated while using the following equation:

$$CM40={\displaystyle \sum }_{i=1}^4(C{F}_i·{\rho }_{Bi}·h_i)$$

(6)

where CF_i is carbon fraction in the ith soil layer, ${\rho }_{Bi}$ is the bulk density of the ith soil layer, and h_i is the height of the ith soil layer.

2.5. Carbon Fluxes Calculations and Modelling

Main C fluxes estimated in this research were net primary production (NPP), total soil respiration (SR), heterotrophic soil respiration (SRH) from dead coarse roots (SRH_CR), and from soil organic matter (SRH_SOM) and decomposition (respiration) of dead wood (R_DW).

Annual NPP was estimated as the sum of the annual production of wood (NPP_W), leaves (NPP_L), fruits (NPP_FT), and fine roots (NPP_FR). NPP_W was estimated as the difference between C stocks in wood biomass (aboveground and belowground) between two measurement points in time at each plot, i.e., winter 2010 and 2011. NPP_L and NPP_FT correspond to annual C stocks in leaves and fruits, respectively, assessed while using litter traps. NPP_FR was estimated starting from C stocks in fine roots and assuming that the live fine root fraction on total (live and dead) fine root was 50% [57,58] and that turnover rate for Pedunculate oak fine root was one year⁻¹ [59].

Total NPP was then modelled using gamma function as proposed in Tang et al. (2014) [60].

$$NPP=k_0·AG{E}^{k_1}·e^{k_2·AGE}$$

(7)

where k₀, k₁, and k₂ are fitting parameters, and AGE is a stand age.

Starting from measured fluxes, total annual soil respiration at each plot was calculated while using a modified soil respiration model [48,61]:

$$SR=\left(a_{AGE}+b_{AGE}·AGE\right)·e^{\left(a_{REW}+b_{REW}·RSWC\right)··\left(\frac{1}{T_{ref}-T_0}-\frac{1}{T_{soil5}-T_0}\right)}·\frac{RSWC}{RSW{C}_{1/2}+RSWC}$$

(8)

where AGE is a stand age, RSWC is a relative soil water content (ratio of actual SWC and SWC at field capacity), T_ref is a reference temperature (i.e., average temperature for research area in the study period, 11 °C), T₀ is the lower temperature limit for the soil respiration (−46 °C), T_soil5 is soil temperature at 5 cm depth, aAGE, bAGE, aREW, bREW, and RSWC_1/2 are model parameters (for more details see Appendix A,

Table A2. Parameter estimates of soil respiration model (Equation (8)) parameterized for stands older than 5 years (i.e., without complete crown cover).

Parameter	Value	Standard Error	t Value	p Value	Lower Conf. Limit	Upper Conf. Limit
R = 0.9069; R2 = 0.8225; n = 1313
aAGE	2.7652	0.0922	29.98	<0.0001	2.5843	2.9461
bAGE	−0.0047	0.0004	−11.16	<0.0001	−0.0055	−0.0039
aREW	157.0005	21.3813	7.34	<0.0001	115.0550	198.9459
bREW	376.0356	30.0037	12.53	<0.0001	317.1749	434.8963
RSWC1/2	0.0985	0.0227	4.34	<0.0001	0.0540	0.1430

and

Table A3. Parameter estimates of soil respiration model * parameterized for youngest 5 years old stand (i.e., without complete crown cover).

Parameter	Value	Standard Error	t Value	p Value	Lower Conf. Limit	Upper Conf. Limit
R = 0.9096; R2 = 0.8274; n = 192
Rref	4.8274	0.7813	6.18	<0.0001	3.2863	6.3685
E0	507.8437	25.1566	20.19	<0.0001	458.2199	557.4675
RSWC1/2	1.0342	0.2886	3.58	0.0004	0.4650	1.6035

*SR model used for the 5-year-old stand is $SR=R_{\mathit{ref}}\times e^{E0\times (\frac{1}{T_{\mathit{ref}}-T_0}-\frac{1}{T_{S5}-T_0})}\times \frac{\mathit{RSWC}}{{\mathit{RSWC}}_{1/2}+\mathit{RSWC}}$.

).

Annual SRH was estimated from total soil respiration measurements and from the decomposition of belowground dead wood biomass (i.e., coarse roots of snags and stumps). Following Bond-Lamberty and Thomson (2012) [62] and Subke et al. (2006) [63], we estimated the mean heterotrophic to total soil respiration ratio (SRH:SR) of 0.514 for temperate deciduous managed forests with mean annual precipitation range of 600–1200 mm, mean annual temperature range of 10–15 °C, and with total soil respiration in range of 3–13 t C ha⁻¹ y⁻¹. However, by using solely a constant ratio, we would have not been able to account for the increased share of SRH after management disturbance (i.e., final cut). Therefore, we divided SRH into heterotrophic respiration from coarse roots of stumps and snags (SRH_CR) and heterotrophic respiration from other soil organic matter (SRH_SOM). We assumed that the decomposition of coarse roots has a decreasing trend during the lifetime of the stand, as proposed in Harmon et al. (2011) [64], while the share of decomposition of soil organic matter in total soil respiration can be considered as constant during the time, as observed by Liu et al. (2015) [65]. SRH_CR was estimated as the decomposition of coarse roots biomass of stumps and snags, thus the equation proposed by Olson (1963) [66] was used:

$$SR{H}_{CR}=CD{W}_{BG}\left(1-e^{-kt}\right)$$

(9)

where CDW_BG is carbon in belowground dead wood biomass (i.e., coarse roots of stumps and snags), k is decomposition rate of 0.182, i.e., an average for fine woody debris (2 mm–7 cm) of naturally occurring tree species in pedunculate oak forest [49], and t is time of decomposition (i.e., 1 year).

SRH_SOM was estimated from remaining total annual soil respiration (i.e., the difference between SR and SRH_CR) using a share of heterotrophic to total soil respiration. The following equation was used:

$$SR{H}_{SOM}=(SR-SR{H}_{CR})·\left(SRH:SR\right)$$

(10)

where SR is total soil respiration and SRH:SR is heterotrophic to total soil respiration ratio of 0.514 (calculated as described above).

Annual R_DW was estimated based on aboveground dead wood biomass and decomposition rates (k) estimated from decomposition experiment [49] and previously published data [67] for fine and coarse woody debris, respectively. The following equation was used:

$$R_{DW}=W_{AGDW}\left(1-e^{-kt}\right)$$

(11)

where W_AGDW is aboveground dead wood biomass, k is decomposition rate, and t is time of decomposition (i.e., 1 year).

Total heterotrophic respiration (RH) was calculated as the sum of SRH_CR, SRH_SOM, and R_DW. The RH was modelled while using a logarithmic function, as follows:

$$RH=a·\ln \left(AGE\right)+b$$

(12)

where a and b are parameters, and AGE is a stand age. Finally, net ecosystem productivity (NEP) was estimated as the difference between NPP and RH.

2.6. Harvest Carbon Losses and Net Ecosystem Carbon Balance

The estimation of past harvests, particularly for old stands, is challenging not only because data on harvest might be missing (oldest stand is 138 years old, while management records are available only for the last 40 years), but also because the area of given forest compartment might have changed. Therefore, we used measured growing stock data to determine the site yield class of our selected stands in the chronosequence, and we then applied prescribed thinnings for specific site class from local species-specific yield tables [68]. Such an approach carries large uncertainty, but, on the other hand, it reflects the usual management practice. The volume of the final cut was estimated from the amount of C that was assumed to be removed from the oldest stands after the final cut by applying basic wood density for oak and C fraction of wood. Finally, harvest loss (V_HL) was obtained as the sum of thinnings and final cut. Cumulative harvest C loss was estimated, as follows:

$$C_{HL}={\displaystyle \sum }_{AG{E}_1}^{AG{E}_2}{V}_{HL}\left(AGE\right)·\rho ·C{F}_W$$

(13)

where V_HL(AGE) is the prescribed harvested wood volume at a given stand age, ρ is basic wood density for oak, and CF_W is the carbon fraction in wood. Typically, dying or standing dead trees are removed during thinnings or final cut. Therefore, C loss from tree mortality can be considered, at least to a great extent, accounted for under harvest C losses.

NECB throughout the rotation was estimated while using two different approaches: the pool-change approach and the component-flux approach. With the pool-change approach, NECB was calculated as the cumulative difference between C stocks in two different stand ages. C stocks in the first year were assumed to be equal to the stocks that remain in the forest after the final cut of the previous stand. For the five years old stand (the youngest stand in our chronosequence), we back-calculated stocks in aboveground live wood biomass, as it would be before the final cut, i.e., at age 0 (so-called “previous stand”), as follows:

$$CL{B}_{AGW}=\left(V_H·\rho ·(1+0.044\right)·CF)$$

(14)

where CLB_AGW is aboveground live wood biomass at age 0, V_H is harvested volume of five years old stand obtained from previous management plan, ρ is species-specific basic wood density, and 0.044 is an average share of biomass of small branches and twigs (W_{D ≤ 3}) in the total aboveground wood biomass (W_AGW) of all trees in 138 years old stand (the oldest stand in our chronosequence).

Applying the relative shares of C stocks, as obtained from 138 years old stand, to CLB_AGW at age 0 we estimated C stocks by all forest ecosystem pools. With the final cut, we assumed the following C stocks allocations: all CLB to 0; 4.4% of CLB_AGW to CDW_AG (representing the biomass of small branches and twigs that were left at the site after harvest); CLB_AGW of new stumps (calculated as CDW_AG in five years old stand while using an average $\rho$ for all occurring tree species instead of wood density for different stages of decomposition, i.e., 3.9 Mg C ha⁻¹) to CDW_AG; CLB_FR; and, CLB_CR to CDW_BG; CLB_L; and, CLB_FT to CFF_OL.

When considering that the final cut did not occur during our experiment, we had to estimate C stocks of a 140 years old stand after the final harvest. To do this, we used data from the 138 years and the 0 years old stands. In short, we estimated the two-year increment of the 138 years old stand, added it to the stock, and then performed a hypothetical final cut and wood removal. Distribution of the remaining carbon stocks by pools in such, 140 years old stand that had been harvested, was modelled based on the shares of carbon in pools of the 0 years old stand. The only difference was that in 140 years CLB_AGW of new stumps was estimated as C stocks of a bole with a diameter at the root collar (d_{0.3 m}) calculated from relation DBH = 0.6571 × d_{0.3 m} (p < 0.01, R² = 0.9953), calibrated on the subsample of 94 trees, with average stump height of 0.24 m (observed at the site) and species-specific basic wood densities.

With the component-flux approach, NECB was calculated as the difference between the cumulative sum of NEP and cumulative harvest C losses.

Finally, NECB from a pool-change approach, i.e., from C stocks (NECB_s), and NECB from component-flux approach, i.e., from C fluxes (NECB_f), were compared.

2.7. Statistical Analysis

Processing of meteorological data and raw data on soil CO₂ efflux measurements and modelling of soil respiration, net primary productivity (NPP), and heterotrophic respiration (RH) were performed with STATA 14 software (StataCorp, College Station, TX, USA) [69]. C stocks and fluxes were expressed, unless differently stated, as average values for each stand with standard errors. One-way ANOVA, followed by Tukey post-hoc test, was performed for C stocks and fluxes across the chronosequence.

3. Results

3.1. Carbon Stocks and Their Temporal Evolution

Table 2. Carbon stocks in the different pools expressed in Mg C ha⁻¹. Mean ± standard error. A number of replicates (No) is different for different pools and across chronosequence due to experimental design. Different letters indicate a significant difference among stand ages (p < 0.05).

C Stocks	Stand Age (Years)									Average (Share)
	0 *	5	13	38	53	68	108	138	140 *
No #	−	5	12	14	12	12	10	12	−	−
CLBAGW	0.0	0.4 ± 0 a	25.2 ± 3.2 a	93.8 ± 6.6 b	106.2 ± 5.1 b,c	129.9 ± 7.3 c	164.7 ± 14.6 d	158.5 ± 5.8 d	0.0	97.0 (40%)
CLBL	0.0	0.2 ± 0.1 a	1.5 ± 0.1 b	1.7 ± 0.1 b	1.9 ± 0.2 b	2.1 ± 0.1 b	2.1 ± 0.1 b	1.8 ± 0.3 b	0.0	1.6 (1%)
CLBFT	0.0	0	0	0.2 ± 0.1 a	1.0 ± 0.2 a,b	0.3 ± 0.1 a	0.5 ± 0.2 a	1.7 ± 0.4 b	0.0	0.5 (0%)
CLBCR&	0.0	0.5 ± 0.2	16.0 ± 1.8	28.2 ± 3.6	27.6 ± 3.2	30.5 ± 2.2	33.0 ± 3.0	31.8 ± 1.3	0.0	23.9 (10%)
CLBFR	0.0	6.0 ± 0.5 a	5.9 ± 0.3 a	2.3 ± 0.3 c	3.7 ± 0.5 b,c	3.3 ± 0.3 b,c	4.6 ± 0.3 a,b	2.8 ± 0.2 c	0.0	4.1 (2%)
CLB	0.0	7.2 ± 0.6 a	48.6 ± 4.8 b	126.2 ± 16.7 c	140.4 ± 15.8 c,d	166.1 ± 11.7 d	204.9 ± 18.1 e	196.6 ± 7.8 e	0.0	127.1 (53%)
CDWAG	13.3	3.7	2.1 ± 0.4 a	8.0 ± 1.2 b	3.9 ± 0.7 a,b	7.8 ± 2.4 a,b	3.0 ± 0.8 a,b	4.6 ± 1.3 a,b	19.5	4.7 (2%)
CDWBG	37.1	24.9	25.2 ± 3.6 a	12.5 ± 0.8 b	9.3 ± 1.1 b,c	4.1 ± 0.8 c	5.9 ± 0.6 b,c	11.1 ± 1.0 b,c	45.7	13.3 (6%)
CDW	50.4	28.6	27.3 ± 4.0 a	20.5 ± 1.8 a,b	13.2 ± 1.4 b	11.9 ± 2.3 b	8.9 ± 1.0 c	15.7 ± 2.0 b	65.2	18.0 (8%)
CFFOL	5.8	3.0 ± 0.4 a,b	3.1 ± 0.3 a,b	2.9 ± 0.3 a,b	3.5 ± 0.3 a,b	2.2 ± 0.1 a	4.2 ± 0.6 b	3.6 ± 0.3 a,b	7.1	3.2 (1%)
CFFOFH	4.7	1.9 ± 0.7 a	0.9 ± 0.2 a	1.8 ± 0.1a	2.8 ± 0.5 a,b	3.4 ± 0.8 a,b	2.7 ± 0.3 a,b	5.8 ± 1.6 b	5.8	2.8 (1%)
CFF	10.5	4.9 ± 0.8 a	4.0 ± 0.4 a	4.7 ± 0.3 a	6.3 ± 0.6 a,b	5.6 ± 0.8 a,b	6.9 ± 0.7 a,b	9.4 ± 1.6 b	12.9	6.0 (2%)
CM_0–5	18.1	14.9 ± 1.2	12.1 ± 0.4	19.4 ± 1.2	16.7 ± 0.5	16.8 ± 1.9	17.7 ± 1.0	22.3 ± 0.9	22.3	17.1 (7%)
CM_5–10	15.7	12.3 ± 1.1	14.3 ± 1.3	20.9 ± 0.9	16.2 ± 0.2	15.4 ± 1.8	12.4 ± 0.7	19.4 ± 1.0	19.4	15.8 (7%)
CM_10–20	27.2	16.6 ± 0.7	25.4 ± 1.2	31.4 ± 1.2	25.2 ± 1.4	27.6 ± 3.1	18.3 ± 1.3	33.6 ± 2.2	33.6	25.4 (10%)
CM_20–40	29.6	20.0 ± 1.4	35.7 ± 3.8	38.9 ± 1.6	28.0 ± 3.3	30.7 ± 1.8	23.6 ± 1.3	36.4 ± 1.2	36.4	30.5 (13%)
CM40	90.6	63.8 ± 2.2 a	87.5 ± 4.2 c	110.6 ± 2.5 d	86.1 ± 3.6 b,c	90.5 ± 4.4 c	72.0 ± 2.2 a,b	111.7 ± 2.8 d	111.7	88.9 (37%)
ΣC	151.5	104.4 ± 2.4 a	167.4 ± 6.7 a	262.0 ± 8.0 b,c	246.0 ± 7.1 b	274.1 ± 9.2 b,c	292.7 ± 6.0 c,d	333.4 ± 7.0 e	189.8	240.0 (100%)

Ecosystem pools: CLB—live biomass, CDW—dead wood, CFF—forest floor, CM—mineral soil at different depth in cm. Ecosystem sub-pools: AG—aboveground, BG—belowground, W—wood, L—leaf, FT—fruit, FR—fine root, CR—coarse root, OL—organic litter soil layer, OFH—organic fragmented and humus soil layer. Letters (a–e) next to numbers indicate significant difference between stands of different age (p < 0.05, One-way ANOVA). * Estimated (i.e., not measured) carbon stocks at the beginning of the first year of the stand development, and at the end of rotation period after the final cut (for calculation see chapter 2.6.); those values are not used in the calculations of the averages. # Number of replicates for CLB_AGW and CDW (variable due to different number of plots); No for CLB_L, CLB_FT, CFF and CM is 4, for CLB_FR is 16, and for CDW in 5 years old stand is 1. & One-way ANOVA was not performed for carbon stocks in coarse roots (CLB_CR) as this flux was not measured but estimated from aboveground carbon stocks.

reports the C stocks in the different pools. C stocks in live biomass increased with stand age and reached its maximum of 204.9 ± 18.1 Mg C ha⁻¹ at the age of 108 years. The C stocks in dead wood were higher in the younger stands (5, 13, and 38 years) than in the mature stand (108 years). C stocks in forest floor showed a linear increase with stand age (R² = 0.90, p < 0.004), while C stocks in the top 40 cm of soil mineral layer, although some significant differences were detected among age classes (Table 2), did not show any significant correlation with stand age (R² = 0.14, p > 0.40). In stands that were older than 50 years, C stocks in forest soil (forest floor and mineral layer) represented 34% of the total.

Most of the C was, on average, stored in live biomass (53%), even though the C distribution among the different pools changed with stand development (Figure 2). In fact, while in the youngest stand (five years), 61% of C was stored in mineral soil layer and the rest of C was mostly stored in dead wood, during stand development, the distribution of C stocks went in favour of live biomass, with the maximum share of 70% reached in the 108 years old stand.

3.2. Temporal Evolution of Carbon Fluxes and Net Ecosystem Productivity

Table 3. Forest carbon fluxes by main forest ecosystem pools and sub-pools, expressed in Mg C ha⁻¹ y⁻¹. Mean ± standard error.

C Fluxes	Stand Age (Years)							Average
	5	13	38	53	68	108	138
NPPAGW	0.12 ± 0.02 a	1.96 ± 0.25 b	3.94 ± 0.34 c	3.90 ± 0.56c	3.58 ± 0.29 c	2.56 ± 0.30 b,c	2.55 ± 0.16 b,c	2.66 ± 0.51
NPPL	0.16 ± 0.10 a	1.47 ± 0.13 b	1.68 ± 0.09 b	1.94 ± 0.20 b	2.11 ± 0.14 b	2.08 ± 0.10 b	1.80 ± 0.25 b	1.61 ± 0.26
NPPFT	0	0	0.16 ± 0.09 a	1.03 ± 0.24 a,b	0.29 ± 0.11 a	0.54 ± 0.18 a	1.70 ± 0.43 b	0.53 ± 0.24
NPPCR	0.10 ± 0.01	1.21 ± 0.12	1.23 ± 0.12	0.99 ± 0.14	0.82 ± 0.07	0.51 ± 0.06	0.52 ± 0.04	0.77 ± 0.16
NPPFR	3.02 ± 0.23 a	2.97 ± 0.15 a	1.13 ± 0.16 c	1.85 ± 0.26 b,c	1.63 ± 0.16 bc	2.29 ± 0.14 a,b	1.39 ± 0.11 c	2.04 ± 0.28
NPP	3.40 ± 0.23 a	7.61 ± 0.32 b	8.14 ± 0.40 b	9.71 ± 0.68 b	8.43 ± 0.36 b	7.98 ± 0.38 b	7.96 ± 0.47 b	7.60 ± 0.75
SR	−9.37	−8.32	−9.49	−8.18	n.a.(−8.00) #	−7.52	−7.19	−8.35 ± 0.38
SRHCR	−4.14 ± 0.23 a	−4.19 ± 0.64 a	−2.08 ± 0.18 b	−1.55 ± 0.20 b,c	−0.68 ± 0.14 c	−0.98 ± 0.11 b,c	−1.85 ± 0.20 b,c	−2.21 ± 0.54
SRHSOM	−2.69 ± 0.11 a,b	−2.12 ± 0.09 a	-3.81 ± 0.15 d	-3.41 ± 0.14 c,d	-3.76 ± 0.15 c,d	-3.36 ± 0.14 c,d	-2.74 ± 0.11 a,b	-3.13 ± 0.24
RDW	−0.21	−0.12 ± 0.03 a	−1.04 ± 0.15 d	−0.45 ± 0.08 a,b,c	−0.68 ± 0.16 c,d	−0.28 ± 0.05 a,b,c	−0.58 ± 0.11 b,c	−0.48 ± 0.12
RH	−7.04 ± 0.25 a	−6.43 ± 0.65 a,b	−6.93 ± 0.28 a	−5.41 ± 0.26 b	−5.12 ± 0.26 b	−4.62 ± 0.18 b,c	−5.17 ± 0.25 b	−5.82 ± 0.37
NEP	−3.64 ± 0.34 a	1.18 ± 0.72 b	1.21 ± 0.49 b	4.30 ± 0.73 c	3.31 ± 0.44 b,c	3.36 ± 0.42 b,c	2.79 ± 0.53 b,c	1.79 ± 1.0

NPP—net primary productivity, SR—total soil respiration, SRH—heterotrophic soil respiration, R_DW—respiration of dead wood in main forest ecosystem pools and sub-pools (AGW—aboveground wood, L—leaf, FT—fruit, CR—coarse root, FR—fine root, SOM—soil organic matter), RH—heterotrophic respiration, NEP—net ecosystem productivity. Different letters (a–d) next to the values indicate statistically significant (p < 0.05) difference between stands. # Value in brackets is linear interpolation of two neighbouring stands (i.e., 53 and 108 years old) and it has not been taken into account when calculating the average.

reports the main ecosystem C fluxes. Forest ecosystem in-coming carbon fluxes are considered to be positive (NPP), while out-going fluxes are considered as negative (SR and R_DW). Mean total NPP along the chronosequence was 7.60 ± 0.75 Mg C ha⁻¹ y⁻¹, with a maximum in the 53 years old stand. A significant difference was only observed for the youngest stand, which showed the lowest value. Leaf biomass production (NPP_L) was relatively uniform throughout the chronosequence, with the exception of the youngest stand, where the canopy was not completely closed. Fruit production showed significant differences among the stands, with the highest production in the oldest stand. Fine root production showed high variability, but no age-related trend was observed. SR showed a linear decrease with stand age (R² = 0.70, p < 0.05), although with low variability (i.e., 11%) among stands. On the contrary, R_DW showed significantly higher spatial variability (66%), but no age trend was observed. In the youngest stand, negative fluxes (i.e., SRH and R_DW) exceed positive NPP flux, which results with a net carbon loss from the ecosystem (i.e., negative NEP values). By the age of 13 years, NEP became positive and with small fluctuation remained positive throughout the chronosequence. The highest carbon sink was observed in 53 years old stand. Although sink seemingly declined with increasing stand age, the observed NEPs in stands of 53 years and older were not significantly (p < 0.05) different (Table 3).

Share of positive and negative carbon fluxes, as well as the share of below- and aboveground carbon fluxes differed between stands of different age (Figure 3). In the youngest stand, SRH_CR accounted for 40% of all fluxes and represented the greatest flux. Second flux in magnitude was fine root productivity that made almost 90% of NPP of the five years old stand. At the age of 13 years, SRH decreased and the share of positive and negative fluxes reached equilibrium at 60:40. By the age of 38 years, the share of above- and belowground NPP fluxes stabilized at 71:29. Share of R_DW was on average 3.6% of all fluxes.

The temporal evolution of the main C fluxes was modelled through time (Figure 4). NPP showed age-related trend and was modelled while using eq. 7 (R² = 0.98, RMSE = 1.2872, n = 28).

Table 4. Parameter values of the nonlinear model (Equation (7)) of net primary productivity (NPP).

Parameter	Value	SE	t-Value	p-Value	CI Lower (95%)	CI Upper (95%)
K0	2.1535	0.5246	4.10	<0.001	1.0730	3.2340
k1	0.4601	0.0870	5.29	<0.001	0.2809	0.6393
k2	−0.0074	0.0018	−4.12	<0.001	−0.0111	−0.0037

gives the model parameters. Total heterotrophic respiration (RH) also showed an age trend (Figure 5) and was modelled using Equation (10) (R² = 0.69, RMSE = 0.6259, n = 6). Parameters of the RH model are given in

Table 5. Parameter values of the logarithmic model (eq. 10) of RH.

Parameter	Value	SE	t-Value	p-Value	CI Lower (95%)	CI Upper (95%)
a	0.6555	0.2200	2.98	0.041	0.0447	1.2662
b	−8.2702	0.8248	−10.03	0.001	−10.5603	−5.9802

. Finally, NEP was estimated from NPP and RH and Figure 6 shows its temporal evolution.

3.3. Harvest Carbon Losses and Net Ecosystem Carbon Balance throughout the Rotation

The cumulative harvest carbon losses (CHL) amounted to 335.18 Mg C ha⁻¹, out of which 187.24 was through thinnings and 147.94 in the final cut. The distribution of CHL among age classes is given in Appendix A (

Table A4. Cumulative harvest carbon losses throughout the rotation.

Age Class	Harvest (m3 ha−1)	Cum. Harvest (m3 ha−1)	Harvest C Loss (Mg C ha−1)	Cum. Harvest C Loss (Mg C ha−1)	Source
20	0	0	0.00	0.00	Local species-specific yield tables (Špiranec 1975) Pedunculate Oak II site class
30	5	5	1.55	1.55
40	17	22	5.27	6.82
50	34	56	10.54	17.36
60	50	106	15.50	32.86
70	62	168	19.22	52.08
80	71	239	22.01	74.09
90	75	314	23.25	97.34
100	76	390	23.56	120.90
110	75	465	23.25	144.15
120	72	537	22.32	166.47
130	67	604	20.77	187.24
140	477	1081	147.94	335.18	Estimated *

* Based on the yield tables (thinnings) and the volume in the oldest, 138 years old, stand in the chronosequence (final harvest).

). NECB calculated from stocks (i.e., pool-change approach) and from fluxes (i.e., component-flux approach) is given in

Table 6. Net ecosystem carbon balance calculated from stocks (NECBs) and from fluxes (NECBf). Mean ± standard error.

Stand Age	ΣC (Mg C ha−1)	NECBs (Mg C ha−1)	Cum. Harvest C Loss (Mg C ha−1)	Cum. NEP (Mg C ha−1)	NECBf (Mg C ha−1)
0	151.5 ± 5.5	0		0.0	0.0
5	104.4 ± 2.4	−47.1 ± 6.0	0	−21.4	−21.4
13	167.4 ± 6.7	15.8 ± 8.7	0	−31.0	−31.0
38	262.0 ± 8.0	110.4 ± 9.7	1.55	11.8	10.3
53	246.0 ± 7.1	94.5 ± 9.0	17.36	59.8	42.4
68	274.1 ± 9.2	122.6 ± 10.7	32.86	113.3	80.4
108	292.7 ± 6.0	141.2 ± 8.2	120.9	253.5	132.6
138	333.4 ± 7.0	181.9 ± 8.9	187.24	340.0	152.7
140	189.8 ± 6.9	38.3 ± 8.9	335.18	345.0	9.8

, while Figure 7 shows its long-term dynamic. The two different approaches showed a clear difference in NECB estimates. Overall, NECB from the pool-change approach (NECBs) was higher than the component-flux approach (NECBf), with a significant difference observed in younger stands (< 50 years). Component-flux gave lower NECB estimates, but with higher uncertainty. NECBs reached greater C loss in the first five years, i.e., −47.1 Mg C ha⁻¹ compared to NECBf C loss of −31 Mg C ha⁻¹ reached at the age of 13 years. According to pool-change approach, the forest became C sink before the age of 13 years, while NECBf remains negative until 35 years. In older stands, the two approaches showed similar NECB values.

4. Discussion

4.1. Carbon Stocks in Different Ecosystem Pools and Their Temporal Evolution

The maximum value of C stocks in live biomass (CLB) observed in our study (204.9 Mg C ha⁻¹) is within the range of published data on deciduous tree species in Europe, e.g., 102 Mg C ha⁻¹ in an oak chronosequence in Austria [11], 172 Mg C ha⁻¹ in a beech chronosequence in Italy [12], 181 Mg C ha⁻¹ in an old-growth mixed forest in UK [13], and 231–233 Mg C ha⁻¹ in an even-aged beech forest in Germany [6]. Understandably, the observed differences in CLB can be linked to differences in stand basal area (BA), where greater BA implies greater CLB and vice versa. Furthermore, stand volume and mean tree height are both very important stand structure elements to be considered in the comparison of this kind. Higher values of CLB observed in Mund (2004) [6] can be attributed to 20% higher stand volume than in our research, while the lower values of CLB in Bruckman et al. (2011) [11] can be explained with lower mean tree height (15–20 %) than in our research. Soil nutrient availability and management differences are other key variables that likely contribute to the observed variability in CLB.

In temperate forest ecosystems, C stocks in dead wood (CDW) are extremely variable [15], and the main controlling factors are stand density, site quality, and forest management frequency and intensity. Average CDW estimated in our research (18 Mg C ha⁻¹) is significantly lower than CDW observed in a managed beech forests in Slovenia (28.4 Mg C ha⁻¹) [70]. Apart from management differences, the stand structure seems to be the likely reason for this. Namely, a significantly higher average number of stumps was observed in beech forest (app. 217 stumps ha⁻¹) compared to our oak chronosequence (average 127 stumps ha⁻¹) [61]. Research dealing with the estimation of CDW sometimes do not account for belowground carbon stocks, i.e., root biomass of stumps and snags [6]. When just aboveground CDW estimated in our research (4.7 Mg C ha⁻¹) are considered, our values are more in line with previous literature data, i.e., 5 Mg C ha⁻¹ in Mund (2004) [6] and 4.2 Mg C ha⁻¹ in Čavlović (2010) [71].

Average C stock in forest floor (CFF) estimated in our study (6 Mg C ha⁻¹) is significantly higher than data published for old black alder (Alnus glutionsa Gaertn.) stand near river Drava in Croatia (2.2 Mg C ha⁻¹) [72], but much lower than value published for 60-year old oak stand in the Netherlands (26.1 Mg C ha⁻¹) [73]. In alder stand, faster decomposition of dead organic matter in forest floor might be attributed to more favourable litter chemical composition (less tannin), higher nitrogen content and higher soil moisture. On the other hand, according to Schulp et al. (2008) [73], forests management activities in the Netherland are less intensive than in many other countries and this might influence CFF accumulation. At our study site, the average thickness of the forest floor was 1.5 cm [61], which is more than 4 times lower than 6.6 cm published in Schulp et al. (2008) [73].

Carbon stocks in mineral soil (CM) depend, among others, on the depth of the soil sampling, which varies among different studies. In this work, we have presented results according to four distinct layers (Table 2) to facilitate comparison with the results from similar works. For example, C stock in the top 15 cm of soil mineral layer (CM_0_15) observed in even-aged beech forests in Germany were 39 and 42 Mg C ha⁻¹ [6] as compared to the average of 45.6 Mg C ha⁻¹ in first 15 cm of soil at our study site. Values for the CM_0_20 (top 20 cm of soil) reported for the beech forest in Italy were disproportionally higher (105 Mg C ha⁻¹) [12] when compared to 58.3 Mg C ha^-1 in top 20 cm of soil in our research. A comparison of the average CM_0_20 in our research with the results from other oak forests showed, instead, good agreement. For example, for the sessile oak forest in Austria, Bruckman et al. (2011) [11] reported 38 Mg C ha⁻¹ for CM_0_20, while Hale (2015) [13] for mixed forests in the UK reported 43 and 41 Mg C ha⁻¹ (for old and young growth, respectively). Contrastingly, Schulp et al. (2008) [73] for the CM_0_20 in the 60-year old oak stand in the Netherlands reported 82 Mg C ha⁻¹, which is significantly higher value compared to the averages above, but closer to the value of 59.8 Mg C ha⁻¹ for the 68 years old stand in our chronosequence (Table 2).

The temporal evolution of C stocks distribution by pools was investigated in many studies [5,6,11,12,15,74,75]. The common result in all these studies is that younger stands, up to approx. 30 years of age, have the highest share of C stored in the soil. With time, in managed forests distribution of carbon stocks changes in favour of live biomass, while the C stocks in dead wood biomass show opposite behaviour with the highest values occurring in young stands and then gradually decline. In the case of the pedunculate oak forest (i.e., our chronosequence), most of the C was stored in soil up to approx. 20 years of age. After stand establishment, the share of CLB increases fast, reaching a peak at around 100 years, after when its share starts to decline. CDW has the highest values after the final cut, and then rapidly declines. Such behaviour is the result of management, where most of the dead or dying trees are harvested and extracted during the thinning. A slight increase in the share of CDW that was observed in the oldest stand is mostly the result of an increased number of larger oak stumps remaining after thinning, which can have a very long turnover time.

4.2. Carbon Fluxes and Net Ecosystem Productivity

Net primary productivity (NPP) is one of the main C fluxes in forest ecosystems, and it mostly depends on tree species [76], stand age [10,15,29], and climate [10,15]. The average NPP estimated in our study (7.6 ± 0.75 Mg C ha⁻¹ y⁻¹) is close to the global estimate for temperate deciduous forests published in Luyssaert et al. 2007 (7.38 ± 0.55 Mg C ha⁻¹ y⁻¹) [19] and Pregitzer and Euskirchen 2004 (7 Mg C ha⁻¹ y⁻¹) [15]. The parametrization of the NPP model, as proposed by Tang et al. (2014) [60], which was used in our study, yielded parameter estimates that were highly significant (p < 0.001). Therefore, the parameterized model may be considered reliable for estimating NPP of Pedunculate oak stands growing in similar conditions and managed with a similar management system.

Soil respiration (SR) is the second most important ecosystem C flux. It accounts for the respiration of live roots and microorganisms, and decomposition of forest floor, dead roots, and other soil organic matter. It is influenced by environmental conditions, with high temporal and spatial variability [77]. The temporal dynamic of SR is highly influenced by the dynamic of meteorological conditions, while spatial variability is largely due to soil properties, the activity of microorganism and of the soil fauna, and growth and turnover of fine root biomass [77,78]. In our study, the spatial variability of SR was around 11% across the chronosequence, which is relatively low value when compared to the broad range (10 to 100%) published in Campbell and Law (2005) [78] and Hanson et al. (1993) [79]. The average SR estimated in our research (8.35 ± 0.38 Mg C ha⁻¹ y⁻¹) is similar to soil respiration of oak stand in Belgium, 8.8 ± 0.5 Mg C ha⁻¹ y⁻¹ [80], and SR of oak chronosequence in China, 6.23–8.31 Mg C ha⁻¹ y⁻¹ [81].

Heterotrophic SR (SRH) is recognized as an important component of SR flux [62] with high likelihood for being affected by the climate change, therefore there is a growing interest for the improvement of SRH modelling [82]. Nevertheless, when estimating SRH under uniform environmental conditions, a simple empirical approach, as we used in this research, could provide satisfactory results. SRH estimated in our research (i.e., the sum of SRH_CR and SRH_SOM) is 5.34 ± 0.38 Mg C ha⁻¹ y⁻¹ and it is close to values that were published in Luan et al. (2011) [81], i.e., 4.31–6.78 Mg C ha⁻¹ y⁻¹. Moreover, the average share of SRH:SR estimated in the present study (65.5%) is within the range 63.2–68.6% published in Liu et al. (2015) [65].

The decomposition of dead wood depends on the debris’ size, soil temperature, and moisture, but also on substrate quality [49,77,83,84]. Substrate quality implies specific chemical composition of the dead wood sample (i.e., share of sugars, proteins, cellulose, lignin, and mineral compounds), which directly depends on tree species and sample size. Higher shares of sugars and proteins in dead wood characterise dead wood of higher substrate quality, which typically exhibits higher rates of decomposition. Decomposition flux also depends on the amount of dead wood C stocks, except for decomposition rates. We were interested in the decomposition flux of aboveground dead wood, as we already included belowground dead wood decomposition in the soil respiration estimates. Decomposition flux R_DW in our study ranged from 0.12 to 1.04 Mg C ha⁻¹ y⁻¹, with a mean of 0.48 Mg C ha⁻¹ y⁻¹ across the whole chronosequence. For comparison, the magnitude of R_DW corresponds to approx. 11% of SRH. The highest R_DW was observed in the 38 years old stand, although the highest C stock in the aboveground dead wood was observed in the oldest stand, where decomposition flux was half of what we estimated in the 38 years old stand. This can be explained with different decomposition rates. According to Ostrogović et al. (2015) [49], one of the main variables explaining the decomposition rates of woody debris is sample size: the smallest samples have the highest decomposition rate, and vice-versa. The higher decomposition flux in the 38 years old stand can be explained with a higher amount of fine woody debris and snags with smaller DBH, than in the other stands (data not shown).

Range of net ecosystem productivity (NEP) in our research spanned from −3.6 to 4.3 Mg C ha⁻¹ y⁻¹, which is close to the values of −2.5 to 5 Mg C ha⁻¹ y⁻¹ that were published for temperate forest ecosystems [7]. In the Pedunculate Oak chronosequence, NEP becomes positive quite early in the development phase, between the stand age of five and 13 years, which is a result of higher NPP fluxes than the total heterotrophic respiration (i.e., the sum of SRH and R_DW). Maximum NEP reached at 53 years, which corresponds to the age of highest NPP. At this age, stands still intensively grow in both diameter and height, but also they start with acorn production that accounts, on average, for 10% of total NPP. Later, with age, NEP slowly decreases as the ratio of positive and negative fluxes goes in favour of negative ones. Nevertheless, until the final cut (140 years), Pedunculate oak stands, managed with a close-to-nature management system, retain positive NEP. A similar trend of NEP was observed in Campbell et al. (2004) [7] and De Simon et al. (2012) [12].

4.3. Net Ecosystem Carbon Balance throughout the Rotation Estimated by Different Approaches

The long-term net ecosystem carbon balance (NECB) that was calculated while using the two different approaches showed relevant differences. The pool-change approach yielded in higher NECB estimates with a higher level of precision, although strongly driven with the first and the last data points (Figure 7). Maximum C loss observed for NECBs in the first five years relies on the accuracy of the C stocks estimated for age 0, while the C accumulation observed at the end of the rotation period, before the final cut, is based on C stocks in the oldest chronosequence stand. On the contrary, the component-flux approach showed lower NECB values, but with significantly higher uncertainty, which usually accompanies all modelling and predicting methods [28]. Nevertheless, if chronosequence method is appropriately applied (i.e., if the basic assumptions are met) and if reliable information on forest management activities is available, it might be an indicator that the observed differences between the two C balance estimates, i.e., NECBs and NECBf, could be due to changes in environmental conditions that have occurred over the last century.

Changes in environmental conditions are linked with an increase in both productivity as well as respiration [4,85,86,87]; therefore, their effect on forest C balance is specific to each forest ecosystem. While NECBs is an integrated measure of the long-term C accumulation under both recent (younger stands) and past (older stands) environmental conditions, NECBf revels C accumulation under current meteorology (in this study, 2011). Having that in mind, higher NECBs that were observed in younger stands might be linked with the recent environmental change (i.e., increased temperature, nitrogen deposition, and CO₂), favouring an increase in stand productivity rather than in heterotrophic respiration. On the contrary, meteorological conditions in 2011 might have had a stronger effect on respiration (i.e., extremely dry year) [88,89], which results results in lower NEP, and finally NECB. Although NECBf carries high uncertainty, it provides valuable information regarding the possible effects that future droughts could have on long-term forest C balances.

5. Conclusions

When carefully set up, chronosequence is a useful method for estimating the temporal evolution of C stocks and/or carbon fluxes. Nevertheless, when using this method for interpreting the net ecosystem carbon balance over the rotation period, one should be aware of the possible differences in estimates that originate from either the pool-change approach or the component-flux approach. When considering all of the shortcomings of the two approaches, in order to obtain a more reliable estimate of long-term C balance, if possible, a comparison of both approaches is recommended.

Under current close-to-nature forest management, Pedunculate oak forest showed to be sustainable in providing both economic and ecological ecosystem services. While assuming successful regeneration, this forest ecosystem becomes C sink very early in a development phase, between the age of 13 and 35, according to the pool-change and component-flux approach, respectively. During 140 years, this forest provides 604 m³ ha⁻¹ (i.e., 187.2 Mg C ha⁻¹) through thinnings and 477 m³ ha⁻¹ (i.e., 147.9 Mg C ha⁻¹) in the final cut, while preserving, on average, 88.9 Mg C ha⁻¹ in the mineral soil down to 40 cm, 18.2 Mg C ha⁻¹ in dead wood, and 6.0 Mg C ha⁻¹ in the forest floor. Finally, soil C stocks in our chronosequence did not show any age-related trend, which indicated that current management practice has no negative effect on soil carbon stocks.