Eddy Covariance vs. Biometric Based Estimates of Net Primary Productivity of Pedunculate Oak (Quercus robur L.) Forest in Croatia during Ten Years

Abstract

We analysed 10 years (2008–2017) of continuous eddy covariance (EC) CO₂ flux measurements of net ecosystem exchange (NEE) in a young pedunculate oak forest in Croatia. Measured NEE was gap-filled and partitioned into gross primary productivity (GPP) and ecosystem reparation (R_ECO) using the online tool by Max Planck Institute for Biogeochemistry in Jena, Germany. Annual NEE, GPP, and R_ECO were correlated with main environmental drivers. Net primary productivity was estimated from EC (NPP_EC), as a sum of −NEE and R_h obtained using a constant R_h:R_ECO ratio, and from independent periodic biometric measurements (NPP_BM). For comparing the NPP at the seasonal level, we propose a simple model that aimed at accounting for late-summer and autumn carbon storage in the non-structural carbohydrate pool. Over the study period, Jastrebarsko forest acted as a carbon sink, with an average (±std. dev.) annual NEE of −319 (±94) gC m⁻² year⁻¹, GPP of 1594 (±109) gC m⁻² year⁻¹, and R_ECO of 1275 (±94) gC m⁻² year⁻¹. Annual NEE showed high inter-annual variability and poor correlation with annual average global radiation, air temperature, and total precipitation, but significant (R² = 0.501, p = 0.02) correlation with the change in soil water content between May and September. Comparison of annual NPP_EC and NPP_BM showed a good overall agreement (R² = 0.463, p = 0.03), although in all years NPP_BM was lower than NPP_EC, with averages of 680 (±88) gC m⁻² year⁻¹ and 819 (±89) gC m⁻² year⁻¹, respectively. Lower values of NPP_BM indicate that fine roots and grasses contributions to NPP, which were not measured in the study period, could have an important contribution to the overall ecosystem NPP. At a seasonal level, two NPP estimates showed differences in their dynamic, but the application of the proposed model greatly improved the agreement in the second part of the growing season. Further research is needed on the respiration partitioning and mechanisms of carbon allocation.

1. Introduction

Currently global forests store approximately 30% of total anthropogenic CO₂ emissions [1], however, the potential of forests to act as a carbon sink in the future is uncertain due to possible saturation effect [2], or negative impact of changed environmental conditions on forest productivity [3]. The eddy covariance (EC) technique is a widely used, state-of-the-art method, and it has become a standard in the estimation and monitoring of high frequency (typically half-hourly) carbon and water fluxes within terrestrial ecosystems [4]. Long-term data series of net ecosystem carbon exchange (NEE) or net ecosystem productivity (NEP, NEP = −NEE), in combination with meteorological and biometric measurements, provide invaluable information on the response of forest ecosystem to environmental conditions and climate change [5,6,7,8].

EC data are essential for calibration of process-based models like Biome-BGC and its variants [9,10] and they are also used for calibration and testing of products derived from remote sensing with sensors, such as Moderate Resolution Imaging Spectroradiometer (MODIS) [11]. Unfortunately, the complexity and high maintenance costs of EC systems have limited the number of sites with measurements spanning a decade or longer. For example, there are 101 forest sites in the global network of micrometeorological tower sites (FLUXNET) 2015 database, out of which 44 datasets are ten or more years long, and only 13 of those are datasets for sites classified as Deciduous Broadleaf Forests [12]. The EC method is far from perfect. EC measurements often result with bad data (e.g., due to instrument malfunctions, ambient conditions not satisfying the EC method or precipitation) and even when data are good a series of corrections during the post-processing is needed [13] (see Section 2). Independent measurements (e.g., biometric, soil respiration, etc.) of the ecosystem fluxes within the footprint of the EC tower are highly needed for the understanding of the ecosystem being measured, as well as for the validation of the EC results.

The above and belowground plant growth, or net primary productivity (NPP) [14] of a forest stand, is typically assessed through biometric measurements. NPP can be estimated from sequential measurement of trees attributes (diameter at breast height—DBH, tree height, etc.) and the use of specific allometric equations and parameters (volume equations, height curves, basic wood densities, root-to-shoot ratios). However, some components of the forest NPP are rather difficult and/or costly to measure (e.g., fine roots production, herbivory losses, and changes in the non-structural carbohydrate reserves) and are not frequently assessed [15,16]. On the other hand, NPP cannot be directly measured with EC [14,17]. Additional data processing (e.g., flux partitioning) and measurements/assumptions (e.g., on the share of heterotrophic in the total ecosystem respiration) have to be made in order to estimate NPP from EC. Consequently, there are many studies comparing NEP (i.e., −NEE) from flux measurements with NEP from biometric measurements, where biometric NEP is calculated as a difference between biometric NPP and measured or estimated heterotrophic respiration [18,19,20,21]. On the other hand, some studies provide the comparison of EC and biometric NPP estimates as well [6,22,23].

Comparison of biometric and EC NEP estimates often show poor agreement on a year-to-year basis [18,20,22,24]. The inter-annual variability of forest productivity is a result of a direct response of trees to occurring meteorological conditions [25], but it is also affected by the postponed response due to the carry-over effect [26,27]. The postponed response in plants is often explained by the presence of carbon reserve pools in the form of non-structural carbohydrates (NSCs). This mechanism has evolved to help trees to cope with adverse effects, such as climate extremes and dampens the visible signal in trees growth [28,29,30,31]. Nevertheless, over a longer time period, a multi-year convergence of these two independent estimates has been observed [19,22,32,33].

In this paper, we provide a comprehensive overview of the measurements of CO₂ fluxes and biometric measurements that were conducted in a young, managed forest of pedunculate oak during the period from 2008 to 2017. We hypothesise that over the time period of 10 years there is a good correlation between EC and biometric annual NPP estimates. The aims of this study were to: (1) quantify variability of CO₂ fluxes from EC with respect to observed meteorology, (2) estimate NPP with biometric method and analyse growth dynamics for different tree species; and, (3) assess the agreement in NPP estimates from eddy covariance (NPP_EC) and biometric measurements (NPP_BM).

2. Materials and Methods

2.1. Research Area and Site Measurement Design

The research was carried out in pedunculate oak stands of Jastrebarsko forest, located approximately 30 km SW from Zagreb, Croatia, near the town of Jastrebarsko (Figure 1). Jastrebarsko forest is part of the Pokupski basin (i.e., river Kupa basin) forest complex, which is, with its area of approximately 13,600 ha, the second largest complex of pedunculate oak forests in Croatia. The river Spačva basin in the east, spanning over 40,000 ha, is the largest forest complex of pedunculate oak in Croatia and in Europe [34] (Figure 1). The area covered with forest in Croatia is 2.493 Mha or 44% of the total land area [35]. Pedunculate oak forests are among the most productive as well as ecologically and commercially most important stands, with the pedunculate oak alone having a standing stock of 48.4 Mm³ (11.6% of the total growing stock in Croatia) [35].

The terrain of the Pokupski basin is mainly flat, with the elevation ranging from 106 m above sea level at the central part of the basin up to 120 m and 130 m above sea level in the south-western and northern parts, respectively. The soil is deep (>5 m), acidic in the top 1 m layer (pH in H₂O 4.9–5.7), but it becomes neutral (pH in H₂O 6.9–7.5) at depths > 1.2 m. It is mainly composed of heavy clay (54% clay, 18% silt, 28% sand) with low vertical water conductivity, and according to the Food and Agriculture Organization World Reference Base (WRB) for soil resources [36], it is classified as luvic stagnosol. Parts of the forest are waterlogged or flooded with stagnating water during winter and early spring, while during summer the soil dries out. The basin is intersected with several natural creeks as well as a network of human-made drainage canals. The average groundwater table depth in the basin ranges from 0.6 m to 2 m, with the average of ~1.5 m at the study site, and seasonal oscillation in between 0 m and 4 m [37]. According to the Köppen classification, the climate of the area is temperate oceanic climate (Cfb), with a mean annual air temperature of 10.6 °C for the period 1981–2010 (data from the Croatian Meteorological and Hydrological Service for the nearest meteorological station in Jastrebarsko). Average annual precipitation during 1981–2010 is 962 mm year⁻¹, out of which around 500 mm falls during the vegetation season (April–September). The dominant tree species in the basin is pedunculate oak (Quercus robur L.), with a significant share of other species, namely black alder (Alnus glutinosa (L.) Geartn.), common hornbeam (Carpinus betulus L.), and narrow-leaved ash (Fraxinus angustifolia Vahl.). Frequently present in the understory are common hawthorn (Crataegus monogyna Jacq.) and hazel (Corylus avellana L.).

The study site is located in managed, young stands (during the time of the study 35–44 years old) dominated by pedunculate oak (Figure 1, left). Stands are the result of regeneration cuts of old (~140 years) pedunculate oak stands in the early 1970s.

The EC tower (coordinates: 45°37′10″ N and 15°41′16″ E) was erected in May 2007 as a part of the Carbon-Pro project [38] and from 21 September 2007 provides continuously meteorological and EC measurements (Figure 1, right). At the time of installation, the measurement height was 23 m (3–5 m above canopy). Because of tree growth, the tower was upgraded and the measurement height was elevated to 27 m on 1 April 2011.

A network of permanent circular plots was set up according to a 100 × 100 m grid (i.e., one plot per ha) around the EC tower during 2007 and winter months of 2008 (Figure A1). All trees were permanently marked with numbers and measured. Tree species, status (alive, dying, dead), and location on the plot were recorded. Stumps and coarse woody debris were also measured [39]. After a preliminary footprint analysis based on EC flux measurements taken from October to December 2007, 24 ‘dendrometer’ plots (marked with red colour on Figure A1) were selected for the installation of dendrometer bands and intensive monitoring. The radius of the 24 ‘dendrometer’ plots (henceforth referred to only as ‘plots’) is 8 m with an area of 201 m² per plot, which is equivalent to sampling intensity of 2% by area.

2.2. Eddy Covariance and Meteorological Measurements

2.2.1. EC and Meteorological Instrumentation

The EC measurement system is made up of a three-dimensional sonic anemometer 81,000 V (R. M. Young Company, Traverse City, MI, USA) and an open path infrared gas analyser—IRGA LI-7500 (LI-COR, Inc., Lincoln, NE, USA). Both instruments are operating at a sampling frequency of 20 Hz. Meteorological variables are measured at 30 s intervals and then averaged and recorded as half-hourly values by a CR1000 data-logger (Campbell Scientific, Inc., Logan, UT, USA). Air temperature and humidity are measured with a humidity and temperature probe (HMP45A, Vaisala Oyj, Vantaa, Finland), incoming and outgoing photosynthetic photon flux density—PPFD is measured with a quantum sensor (LI190SB, LI-COR, Inc., Lincoln, NE, USA), net radiation is measured with a net radiometer (NR Lite, Kipp & Zonen, Delft, The Netherlands), while incoming shortwave radiation is measured with the pyranometer (CMP 3, Kipp & Zonen, Delft, The Netherlands). Soil temperature is measured at three depths (5 cm, 15 cm, and 25 cm) with type T thermocouples. Average soil water content in the 0–30 cm soil layer is measured with two water content reflectometers (CS616, Campbell Scientific, Inc., Logan, UT, USA).

CO₂ storage term was measured with a profiler system for the CO₂ concentration measurement at six levels: 1, 2, 4, 8, 16, and 23 m (until 30 March 2011), and 27 m (from 1 April 2011). The profiler consisted of six tubes, pump, set of valves, 16 channel relay multiplexer (AM16/32B, Campbell Scientific, Inc., Logan, UT, USA), the CO₂ analyser (SBA4, PP System, Amesbury, MA, USA), and was controlled by the CR1000 data-logger. However, due to occasionally compromised data from the profiler (due to ants invading the filters of the tubes) and a final general failure of the profiler system in 2014, from 2014 until 2017 we decided to compute the storage term using only the CO₂ concentration measured at the top of the tower, as suggested by [40].

In the cases of gaps in meteorological measurements, we used half-hourly meteorological data (air temperature, global radiation, relative humidity, and precipitation) from a small, auxiliary meteorological station (WatchDog 2900ET, Spectrum Technologies, Aurora, IL, USA), which was installed within the same forest, approximately 2 km NE from EC tower.

2.2.2. EC Data Processing

Raw EC data were stored in one-hourly binary files (

Table A1. Percentage of half-hourly NEE files before processing, after quality control filtering procedures and after u* filtering procedure.

Year	Data Coverage before Processing [%]	Data Coverage after Quality Control [%]	Data Coverage after u* Filtering [%]
2008	95.5	46.6	39.4
2009	88.2	46.8	42.1
2010	91.9	47.0	45.4
2011	94.3	56.5	47.7
2012	83.3	49.8	42.7
2013	92.8	53.0	48.4
2014	95.7	51.4	47.3
2015	96.0	52.9	43.2
2016	83.8	45.1	41.6
2017	79.0	45.0	38.0

). Under normal circumstances, every one-hourly binary file has 72,000 records and files having less (or more) records than 72,000 have been discarded. The remaining one-hourly binary files had been further quality checked. For every one-hourly binary file mean value, the standard deviation of CO₂ and a difference (‘delta’) between two consecutive CO₂ concentration measurements was calculated. The one-hourly binary file was dismissed from further analysis if the standard deviation was greater than 30 ppm, or the mean value was less than 340 ppm or greater than 600 ppm, or the number of events when delta > 10 ppm was larger than 10. This procedure removed bad data caused by rain and system malfunctions alongside with non-physical shifts, dropouts, and discontinuities described in Vickers and Mahrt [41].

Fluxes were calculated using the free software EdiRe (http://www.geos.ed.ac.uk/homes/jbm/micromet/EdiRe) according to the project EUROFLUX methodology [42] and stored as half-hourly values. At first, to avoid the loss of flux due to the physical separation of wind speed and scalar concentration sensors, the time lag was removed. In open-path EC measurement systems gas analyser and three-dimensional (3D) anemometer are usually placed several decimetres (30 cm in our case) apart so the time lag is much smaller than in case of closed path systems. To avoid contamination of vertical wind component w by the other two wind components planar fit method [43] was performed after time lag removal. CO₂ flux was calculated as the average of the product of instantaneous fluctuations of vertical wind velocity (w’) and CO₂ concentration fluctuations (s’):

$$Fc=\overline{w^{\prime }{s}^{\prime }}$$

(1)

Calculated flux Fc was first corrected for the imperfect frequency response of the system. After the spectral corrections, the Webb-Pearman-Leuning (WPL) correction was performed to compensate for fluctuations in the density of carbon dioxide resulting from the fluctuations of air temperature and water vapour, which are not representatives of measured flux [42]. Only data that have passed the quality criteria described above were used in the calculation of fluxes.

To obtain NEE, the CO₂ storage flux (Sc) was added to the calculated flux Fc. Half-hourly NEE were checked for absolute limits. Values that were outside the range of −50 µmoL m⁻² s⁻¹ > NEE > 30 µmoL m⁻² s⁻¹ have been removed. High-frequency spikes, i.e., spikes affecting the single instantaneous measurement have been already removed in EdiRe. However, spikes were sometimes still occurring in the final half-hourly fluxes. They were detected and removed by the filtering procedure proposed by Papale et al. [44]. Half-hourly NEE values were categorized into night-time (solar radiation (Rg) < 20 Wm⁻²) and day-time NEE. The spike detection was based on the double-differenced time series, using the median of absolute deviation around the median (MAD):

$$MAD=\mathrm{median}\left(\left|d_i-M_d\right|\right)$$

(2)

where M_d denotes the median of the differences and value d_i was calculated for each half-hourly NEE as:

$$d_i=\left(NE{E}_i-NE{E}_{i-1}\right)-\left(NE{E}_{i+1}-NE{E}_i\right)$$

(3)

Half-hourly NEE was flagged as spike and removed from further processing if:

$$d_i<M_d-\left(\frac{z\times MAD}{0.6745}\right) \mathrm{OR}{d}_i>M_d+\left(\frac{z\times MAD}{0.6745}\right)$$

(4)

where z is a threshold value. As proposed in Papale et al. [44], the threshold was set to z = 4 during the night-time and to z = 5.5 during the day-time to allow for greater variability of daily fluxes.

During the stable conditions, when the stratification is stable and the wind speed is low, turbulence may not be fully developed. In such situations, the hypothesis underlying eddy covariance method cannot be met [42,45]. These occasions occurred mostly at night, when the ecosystem acts as a source of carbon. Eddy covariance method underestimates CO₂ flux under the stable conditions, which lead to the overestimation of NEE at the annual scale [42]. To check the reliability of night-time flux measurements at Jastrebarsko site, flux data was divided into 20 classes with respect to friction velocity u* (from 0 to 1 m s⁻¹ with increment of 0.05 m s⁻¹) and grouped according to the season (leaf-off, leaf-on; in order to account for differences in surface roughness), and air temperature class (2 °C wide temperature classes). For the each u* class average and median of NEE was calculated. Neither the average nor the median of night-time NEE showed significant trends with u* (data not shown). Therefore, we could not select a fixed u* threshold for the u* filtering, which is used for identifying low turbulence conditions. Instead, we used the value for the u* threshold calculated by the online gap-filling tool by the Department Biogeochemical Integration at the Max Planck Institute for Biogeochemistry in Jena, Germany [17,46].

The fluxes and meteorology data for Jastrebarsko oak forest (half-hourly eddy covariance and meteorological measurements, corrected daily precipitation, variable description, and processing codes) are provided in the S1 supplementary material online.

2.3. NEE Flux Partitioning and the Estimation of NPP_EC

Complex filtering procedures alongside with measurement system malfunctions can lead to a significant amount of gaps in time series. The typical data coverage within a year at an EC site is between 65% and 75% [47]. Average annual coverage with flux data of acceptable quality at Jastrebarsko site was even smaller (from 38.0% to 48.4%, see Table A1), among others due to frequent stable conditions and fog causing a condensation on windows of the LI-7500 IRGA. To estimate the annual carbon budgets missing or discarded flux values were replaced with estimates using standardized marginal distribution sampling (MDS) technique that was employed by FLUXNET, which is available as an online tool [17,46]. Gap-filled NEE fluxes were partitioned into gross primary production (GPP) and ecosystem respiration (R_ECO) while using the same online tool [46,48]. This tool uses the night-time based respiration model of Lloyd and Taylor [49] for the assessment of the R_ECO:

$$R_{ECO}\left(T\right)=R_{ECO}{}_{ref}\times e^{E_0(\frac{1}{T_{ref}-T_0}-\frac{1}{T-T_0})}$$

(5)

where T_ref marks reference temperature that is set to 15 °C, R_ECO,ref denotes respiration at the reference temperature, T₀ is set to −46.02 °C, T marks measured air or soil temperature, while the activation energy parameter, E₀, is allowed to vary [48]. Air temperature (T_a) varies more than soil temperature (T_s) and more variance in R_ECO can be explained by T_a [50]. Therefore, T_a was selected for partitioning of NEE into GPP and R_ECO. After R_ECO has been calculated, GPP is determined as:

$$GPP=-NEE+R_{ECO}$$

(6)

R_ECO consists of autotrophic (R_a) and heterotrophic (R_h) parts. At the EC site, R_h was estimated in 2008 and 2009 using soil respiration measurements [38] and the assumption that the share of R_h in soil respiration is 50% [19,38]. Using data for R_h published in [38] and R_ECO obtained from NEE flux partitioning for 2008 and 2009, the ratio R_h:R_ECO was calculated separately for 2008 and 2009. The ratio R_h:R_ECO from the data for 2008 and 2009 was then averaged and this ratio (39.19%) was subsequently used for the partitioning of R_ECO into R_a and R_h for all of the remaining years (2010 to 2017). Net primary production from eddy covariance (NPP_EC) was calculated as the sum of measured NEP (i.e., −NEE) and calculated heterotrophic respiration:

$$NP{P}_{EC}=NEP+R_h=-NEE+R_h=-NEE+0.3919\times R_{ECO}$$

(7)

2.4. Assessment of Flux Footprint

To estimate the position and spatial extent of the area that is contributing to the measured carbon flux, a Lagrangian stochastic particle dispersion model was used [51]. In the Lagrangian formulation, the diffusion of scalar released from the surface is assumed to be statistically equivalent to the dispersion of an ensemble of particles that impact the ground within the surface area source and thereafter transport the scalar [52,53]. The model uses stochastic backward trajectories of particles, which allows for the calculation of the footprint for a measurement point instead of an average over a sensor volume without a coordinate transformation that requires horizontal homogeneity of the flow [53]. Recently improved parameterisation of the model describes the flux footprint function’s upwind extent and the crosswind spread of the footprint [51]. Data used in simulation are wind velocity, wind direction, measurement height, surface roughness, friction velocity, Monin-Obukhov length, and zero-plane displacement. Model is available as an online tool [54]. Flux footprint analysis showed that the fetch for the 90% contribution to measured flux before the elevation of the tower in 2011 was ~350 m, and after the tower was elevated it was ~500 m (Figure 2).

In addition, the footprint analysis confirmed that the contribution to the measured fluxes arriving from areas outside of the forest (i.e., agricultural areas, highway, and the town of Jastrebarsko with distances from EC tower of 1.2 km W, 2.4 km NW, and 5 km NW, respectively) is likely negligible. All 24 circular plots with installed dendrometer bands on trees (see Section 2.5) lay within the footprint of EC tower.

2.5. Biometric Measurements and Estimation of NPP_BM

Estimation of forest stand NPP through the biometric method required the assessment of different components, namely the production of: stems and branches, leaves, flowers and fruits, coarse and fine roots, ground vegetation (grasses and non-woody plants), herbivory losses, pollen, volatile organic compounds (VOCs), and non-structural carbohydrates (NSC). Some of the mentioned NPP components are very difficult and/or costly to measure (VOCs, NSC, herbivory), while others (e.g., pollen) have a fairly insignificant contribution to the total NPP and could be disregarded. Within this study, we estimated the production of total woody biomass (NPP_WBt) by components: stem and branches (NPP_SB), twigs (NPP_T), and roots (NPP_R). NPP of leaves and fruits (NPP_LF) was estimated from litterfall measurements. The sum of NPP_WBt and NPP_LF is denoted as NPP_BM. Net primary production of a given component is expressed in units of gC m⁻², while net primary productivity is expressed in units of gC m⁻² year⁻¹. More details on measurements and calculations are given below.

2.5.1. NPP of Total Woody Biomass, Leaves and Fruits

From 2007 until 2017, at the end of growing season, the DBH of all trees (DBH > 2 cm) in each of the 24 plots (Figure 2 and Figure A1) were measured with a measuring tape at 1 mm precision. In addition, in April 2008, a total of 643 aluminium dendrometer bands were installed on all trees on the plots having DBH larger than 7.5 cm for the assessment of short-term (weekly to monthly) stem increments. Dendrometer bands have been self-made in the Croatian Forest Research Institute according to the method that was described by Keeland and Young [55]. During the vegetation season, cumulative stem circumference increments were measured on dendrometers with a precision of 0.01 mm, using small callipers with an electronic display. The frequency of measurements varied from weekly to monthly during the growing season, with lower frequency in most recent years due to the limitation in resources. The error of the measurement was estimated to be around 1% of the measured increment [56]. The effect of the initial ‘settling’ of the dendrometer band on a tree stem during 2008, which, if unaccounted for, could result with underestimation of increment during the first year, has been estimated and the increment was corrected, as described in [56].

Tree heights (h) of all trees were also measured each year (with the exception of 2011) using Vertex III hypsometers (Haglöf Sweden AD, Långsele, Sweden). This provided sample data for constructing species-specific, age-dependent height curves using adjusted Michailoff function [57] (

Table A2. Parameters for species-specific, age-dependent tree height curves ($h=\left({\alpha }_0+{\alpha }_1\times Age\right)\times e^{\frac{-\left({\beta }_0+{\beta }_1·Age\right)}{DBH}}+1.30$).

Species **	Parameters *
	${\mathit{\alpha }}_0$	${\mathit{\alpha }}_1$	${\mathit{\beta }}_0$	${\mathit{\beta }}_1$
Alnus glutinosa	10.31	0.37	0.59	0.1385
Carpinus betulus	6.60	0.59	−2.31	0.2863
Fraxinus angustifolia	6.25	0.58	−3.89	0.3275
Quercus robur	3.86	0.63	−1.61	0.2496

* Based on 35–44 years old stands. ** For all other tree species the parameters for Carpinus betulus were used.

). Volume (V _{≥ 3 cm}), over bark, of tree stem, and branches (d ≥ 3 cm), was assessed with the allometric equation of Schumacher and Hall [58] using local, species-specific parameters (

Table A3. Parameters of Schumacher and Hall function * ($V=b_0\times {DBH}^{b_1}\times h^{b_2}$) for the wood volume (including bark) of tree stem and branches **.

Species ***		Parameter		Source
	b0	b1	b2
Alnus glutinosa	4.23243∙10−5	2.002354	1.001300	[40]
Carpinus betulus	2.96400∙10−5	2.022705	1.102119	[41]
Fraxinus angustifolia	3.95282∙10−5	1.974875	1.001444	[42]
Quercus robur	4.96820∙10−5	2.048384	0.892124	[41]

* Unites: DBH (cm), h (m), V (m³). ** Wood from the parts of branches that are larger than 3 cm in diameter over bark. *** For the tree species that were present, but not listed here, the parameters of Carpinus betulus were used for the calculation of the wood volume.

) [59,60,61]. The volume of thin branches and twigs (d < 3 cm in diameter) was assumed to be 5% of V _{≥ 3 cm} [62]. The volume of below ground coarse roots (>2 mm) was calculated assuming a constant, species-invariant root-to-shoot (RS) ratio of 0.257. The value of RS ratio was obtained as the average of the RS ratios for stands with pedunculate oak published in [63]. Total tree volume was converted to biomass using species-specific basic wood densities (BWD), namely 450, 630, 570, and 580 kg m⁻³ for Alnus glutinosa, Carpinus betulus, Fraxinus angustifolia, and Quercus robur, respectively [64], and converted to carbon using a carbon fraction (CF) of 0.50 [65,66]. Carbon stock (C) at the time of measurement (t) was calculated for each component of tree i as:

$$\begin{array}{cc}\hfill \text{stem and branches}:& C_{S+B, i}\left(t\right) = V_{\ge 3 cm,i}\left(t\right) \times BWD \times CF\}\hfill \\ \hfill \text{twigs}:& C_{T, i}\left(t\right) = 0.05 \times C_{S+B, i}\left(t\right)\hfill \\ \hfill \text{roots}:& C_{R, i}\left(t\right) = RS \times \left(C_{S+B, i}\left(t\right)+C_{T, i}\left(t\right)\right)\hfill \\ \hfill \text{total tree woody biomass}:& C_{WBt, i}\left(t\right) = \left(C_{S+B, i}\left(t\right)+C_{T, i}\left(t\right)+C_{R, i}\left(t\right)\right).\hfill \end{array}\}$$

(8)

The NPP was calculated from changes in C_i at a given plot. For example, the NPP_WBt is calculated as:

$$\begin{array}{c}NP{P}_{WBt}=\overline{NP{P}_{WBt}\left(plot\right)}=\\ =\frac{1}{N_{plots}}{\displaystyle \sum }_{plot}\frac{1}{A\times \left(t_2-t_1\right)}\times \left[{\displaystyle \sum }_{i=1}^{N_{alive}}\left(C_{WBt, i}\left(t_2\right)-C_{WBt, i}\left(t_1\right)\right)+{\displaystyle \sum }_{j=1}^{N_{dead or removed}}\left(C_{WBt, j}\left(t_j\right)-C_{WBt, j}\left(t_1\right)\right)\right]\end{array}$$

(9)

where N_plots is the number of plots, A is the plot area, t₁ and t₂ are times of the beginning and the end of the period of interest, and $C_{WBt,j}(t_j)$ is estimated carbon stock of tree j at the time $t_j$ of its death or removal ($t_1\le t_j\le t_2$).

From 2008 until 2015, 16 litterfall baskets (45 cm in diameter) that were placed at centres of randomly selected plots, were used for the assessment of leaves and litter production. In the early spring of 2016, additional litterfall baskets were installed resulting with a total of 39 installed baskets. Litterfall was collected from baskets several times a year. On several occasions, some of the litterfall baskets were damaged by wild boars (probably attracted by acorns that fell in a basket). In such cases, for that particular season, data from the damaged basket were discarded.

The collected litterfall was sorted into three fractions: leaves, small twigs (d < 1 cm), and fruits. Samples were dried and weighed with 0.01 g precision. Value of biomass (dry weight) was divided by the area of the sampler and multiplied by the carbon fraction of 0.50 to obtain the production from leaves (${NPP}_L$), fruits (${NPP}_F$), and their sum (${NPP}_{LF}$).

2.5.2. Modelling Seasonal Dynamics of NPP_BM

The stem increment of deciduous trees growing in the continental parts of Europe starts in the spring before the new leaves develop. Part of that stem increment is due to tree re-hydration and it should not be considered as growth in terms of increase in carbon stock. Furthermore, although all of the carbon in leaf biomass is usually accounted as part of the current season’s NPP, a part of that carbon actually comes from tree reserves of the previous season(s) and not from carbon that is fixed in the current season [28]. On the other hand, stem increment almost ceases in summer, particularly after the mid of August [6]. However, from the CO₂ flux measurements we can observe that the forest continues to accumulate carbon; apparently more than the amount required for the growth of stems and branches (see Section 3). The fixed carbon is likely stored in tree reserves [6]. These facts should be addressed when comparing seasonal dynamics of NPP estimates from eddy covariance and biometric measurements.

Since the measurement and assessment of changes in trees’ reserves is very challenging, and it was beyond the scope of this paper, in order to account for the seasonal NPP_LF we assumed a simplistic model where the amount of carbon mobilised from trees’ reserves and used in early spring for the formation of new tissues (primarily that of leaves), is approximately equal to the NPP_LF. Consequently, the NPP_LF of the current season should not be accounted for in the spring, but only later in the season, when the reserves are being replenished. The underlying idea is to account, in a simple and transparent way, for the effects of NSC storage changes. The actual mechanism of carbon reserve management in trees is significantly more complex and is currently still being researched [6,16,23,25,67,68,69].

In our simple model, we assumed that the production of every tree component under the average growing conditions exhibits the pattern, which can be described by the logistic curve of the form:

$$NP{P}_i\left(DOY\right)=\frac{NP{P}_{i\_annual}}{1+e^{-k\left(DOY-DO{Y}_{max}\left(i\right)\right)}}$$

(10)

where NPP_i(DOY) is the cumulative net primary production of the i-th component (e.g., stem and branches, leaf, fruit) from the first until a given Day-Of-Year (DOY), NPP_i__annual is the annual NPP of the component i, DOY_max(i) is the DOY when the rate-of-change of NPP(i, DOY) is at its maximum (i.e., the inflexion point of the annual NPP curve), and k is a parameter that is related to the maximum rate-of-change of NPP (r_max) where r_max is the quotient of the maximum daily net primary production (dNPP_max) and the annual net primary production (NPP_annual):

$$r_{max}=\frac{dNP{P}_{max}}{NP{P}_{annual}}.$$

(11)

If r_max is known (or estimated) the parameter k can be calculated using the following equation:

$$k=4\times atanh\left(r_{max}\right)=2\times ln\left(\frac{1+r_{max}}{1-r_{max}}\right)$$

(12)

Inversely,

$$r_{max}=tanh\left(\frac{k}{4}\right)=\frac{e^{k/4}-e^{-k/4}}{e^{k/4}+e^{-k/4}}$$

(13)

Using the nl routine for nonlinear least-squares estimation in STATA 14 statistical software (StataCorp LP, College Station, TX, USA), we fitted the NPP_WBt with the Equation (10) for each year separately and obtained estimates for the parameters NPP_WBt,annual, DOY_max__WBt and k (

Table A4. Parameters of the NPP_WBt fit with Equation (10).

Parameter	Year	Estimate	Std. Err.	t	p > |t|	CI95lower	CI95uper
DOYmax_WBt	2008	151.6	0.5	317.43	0.000	150.7	152.6
	2009	154.1	0.7	227.61	0.000	152.7	155.5
	2010	158.3	1.1	138.99	0.000	155.9	160.7
	2011	144.7	1.6	91.91	0.000	141.3	148.0
	2012	143.2	1.3	108.31	0.000	140.3	146.1
	2013	152.0	2.0	75.94	0.000	147.6	156.4
	2014	156.3	1.2	126.73	0.000	153.5	159.0
	2015	155.3	1.7	93.04	0.000	151.5	159.1
	2016	154.4	2.7	57.64	0.000	146.9	161.8
	2017	163.5	0.3	590.49	0.001	160.0	167.0
k	2008	0.0476	0.0011	42.77	0.000	0.0453	0.0499
	2009	0.0426	0.0012	36.81	0.000	0.0402	0.0450
	2010	0.0397	0.0015	26.67	0.000	0.0365	0.0428
	2011	0.0431	0.0027	15.83	0.000	0.0373	0.0490
	2012	0.0356	0.0015	23.36	0.000	0.0323	0.0389
	2013	0.0287	0.0017	16.5	0.000	0.0249	0.0325
	2014	0.0322	0.0010	31.4	0.000	0.0299	0.0345
	2015	0.0404	0.0026	15.63	0.000	0.0346	0.0463
	2016	0.0293	0.0019	15.24	0.000	0.0240	0.0346
	2017	0.0307	0.0001	248.18	0.003	0.0291	0.0322
NPPWBt annual	2008	504.9	2.3	216.26	0.000	500.1	509.8
	2009	550.9	3.7	150.55	0.000	543.3	558.6
	2010	571.5	7.5	76.41	0.000	555.7	587.3
	2011	513.1	9.7	52.89	0.000	492.3	533.9
	2012	457.0	5.3	85.69	0.000	445.3	468.6
	2013	540.4	9.1	59.64	0.000	520.5	560.4
	2014	486.5	4.8	102.17	0.000	475.9	497.1
	2015	409.9	6.3	65.16	0.000	395.7	424.1
	2016	526.9	11.6	45.46	0.000	494.8	559.1
	2017	302.0	0.3	900.4	0.001	297.7	306.2

). Seasonal dynamics of NPP_LF was estimated assuming that the shape parameter k for the NPP_LF is the same as for NPP_WBt. Parameter NPP_LF,annual equals the measured NPP_LF. Looking into the characteristics of the logistic curve, we propose that the DOY_max__LF occurs on DOY, at which the time derivative of NPP_WBt has the second inflexion point (Figure A2A), resulting with the general behaviour of NPP_BM analogous to that shown in Figure A2B. In that case, DOY_{max_LF} must be the value in which the second derivative of the logistic function has its minimum and the third derivative is zero.

An analytical solution for DOY_{max_LF} in this case exists and it is:

$$DO{Y}_{max\_LF}=DO{Y}_{max\_WBt}-\frac{\ln \left(2-\sqrt{3}\right)}{k}.$$

(14)

Finally,

$$NP{P}_{BM}\left(t\right)=NP{P}_{WBt}\left(t\right)+NP{P}_{LF}\left(t\right)$$

(15)

3. Results

3.1. Meteorological Conditions during the Period 2008–2017

Overall mean air temperature at the study site during the period 2008–2017 was 11.24 °C. The year with the lowest mean annual air temperature (10.22 °C) was 2010, which was the only year with a value below the 30-year (1981–2010) average (10.62 °C) (data from the Croatian Meteorological and Hydrological Service—CMHS for the meteorological station Jastrebarsko, 6 km NW of EC tower). The year 2014 had the highest average annual air temperature (12.11 °C;

Table 1. Average tree age (years), mean annual precipitation (P, mm), mean annual air temperature (T_a, °C), mean annual soil temperature (T_s, °C), mean annual global radiation (R_g, W m⁻²), mean annual soil water content of the 0–30 cm soil layer (SWC, vol_H2O vol⁻¹), difference between mean soil water content in May and in September (SWC_May–SWC_Sep, vol_H2O vol⁻¹), start of growing season (SGS, day of year (DOY)), end of growing season (EGS, DOY), annual carbon fluxes (net ecosystem exchange (NEE), gross primary productivity (GPP), ecosystem reparation (R_ECO), from eddy covariance measurements, gC m⁻² year⁻¹), and autotrophic and heterotrophic respiration (R_a and R_h, respectively, gC m⁻² year⁻¹). NPP_EC and NPP_BM are two independent estimates of the annual net primary productivity (NPP, gC m⁻² year⁻¹). NPP_EC is calculated from eddy covariance measurements and from heterotrophic respiration estimates (R_h, see text). NPP_BM is calculated from biometric measurements and it is the sum of NPPs of individual tree components: SB (stems and branches), T (twigs), R (coarse roots), L (leaves), F (seeds), SE is the standard error and SD is the standard deviation of the annual average (Aver.) during the period of the study. Note: Values of SD are in brackets in order to be distinguishable from SE.

Year	Age	P	Ta	Ts	Rg	SWC	ΔSWCMay–Sep	SGS	EGS	NEE ± SE	GPP	RECO	Ra	Rh	NPPEC	NPPSB ± SE	NPPT	NPPR	NPPL ± SE	NPPF ± SE	NPPBM
2008	35	900	11.34	11.12	153	0.448	0.233	77	295	−352 ± 13	1469	1117	679	438	790	387 ± 15	19	104	164 ± 20	1 ± 0	675
2009	36	940	11.25	11.23	155	0.420	0.197	98	293	−496 ± 15	1622	1126	685	441	937	413 ± 18	21	111	207 ± 10	4 ± 1	756
2010	37	1255	10.22	10.26	143	0.518	0.072	86	287	−286 ± 14	1615	1329	808	521	807	430 ± 18	21	116	207 ± 15	9 ± 4	783
2011	38	576	11.05	10.63	160	0.403	0.195	88	295	−353 ± 14	1642	1289	784	505	858	395 ± 17	20	107	199 ± 11	18 ± 5	739
2012	39	987	11.49	11.06	162	0.419	0.071	90	296	−261 ± 14	1642	1382	840	541	802	351 ± 16	18	95	165 ± 8	1 ± 0	630
2013	40	1238	10.92	10.52	148	0.493	0.124	111	291	−356 ± 13	1630	1275	775	500	855	403 ± 26	20	109	204 ± 12	20 ± 6	756
2014	41	1755	12.11	11.82	140	0.586	−0.011	90	272	−232 ± 14	1522	1290	785	506	738	359 ± 17	18	97	174 ± 8	6 ± 2	654
2015	42	1260	11.47	11.10	152	0.502	0.212	86	278	−373 ± 13	1775	1402	852	549	923	310 ± 13	16	84	162 ± 13	60 ± 25	632
2016	43	744	11.11	11.01	146	0.512	0.273	101	282	−339 ± 13	1644	1305	794	512	850	393 ± 19	20	106	166 ± 5	4 ± 1	689
2017	44	927	11.39	10.76	164	0.487	0.049	100	292	−147 ± 13	1384	1236	752	484	632	227 ± 16	11	61	183 ± 5	1 ± 0	483
Aver.		1058	11.24	10.95	152	0.479	0.141	93	288	−319	1594	1275	775	500	819	367	18	99	183	12	680
SD		(332)	(0.48)	(0.43)	(8)	(0.053)	(0.089)	(10)	(8)	(94)	(109)	(94)	(57)	(37)	(89)	(60)	(3)	(16)	(19)	(18)	(88)

).

The coldest months were December, January, and February, with the average monthly temperatures during the study period of 1.68 °C, 0.56 °C, and 2.46 °C, respectively. The absolute minimum air temperature was recorded on 9 February 2012 (−20.68 °C). July and August were the hottest months with the average daily air temperatures typically reaching values higher than 25 °C on several days and average monthly temperatures during the study period of 21.42 °C and 20.76 °C, respectively. The absolute maximum air temperature was recorded on 8 August 2013 (38.92 °C).

Average annual precipitation during the study period was 1058 mm (1981–2010 average, measured at the Jastrebarsko town meteorological station, was 962 mm). The majority of the precipitation falls in summer and early autumn. With an annual precipitation of 576 mm, the year 2011 was the driest year in the study period, while the year 2014 had the highest annual precipitation (1755 mm). Interestingly, the precipitation during the 2014 vegetation season exceeded the 30 years average for the annual precipitation, while the summer rainfall alone was 727 mm (

Table A5. Seasonal meteorological data; average air and soil temperature, average global radiation and sums of precipitation by season.

	Season				Annual	Season				Annual
Year	DJF	MAM	JJA	SON	Values	DJF	MAM	JJA	SON	Values
	Ta (°C)					Ts (°C)
2008	2.55	11.30	20.08	10.90	11.29	3.65	9.69	18.09	12.25	11.12
2009	0.87	12.45	20.02	11.49	11.28	3.70	10.53	17.65	12.56	11.24
2010	0.75	10.91	20.10	10.13	10.27	3.80	9.59	17.69	11.42	10.31
2011	0.24	11.65	20.47	10.14	11.04	2.20	9.82	18.12	11.23	10.63
2012	0.63	12.16	21.81	11.89	11.49	3.03	10.31	18.58	12.73	11.06
2013	0.73	10.36	20.62	11.26	10.92	2.01	9.63	17.89	12.07	10.52
2014	3.85	12.42	19.30	12.20	12.10	4.86	11.19	17.93	12.90	11.82
2015	2.74	11.67	21.27	10.75	11.54	3.43	10.29	18.78	11.73	11.11
2016	2.91	11.36	20.19	10.82	11.11	3.78	9.99	18.51	11.96	11.23
2017	−0.40	12.57	21.62	10.65	11.57	2.06	10.48	18.49	11.98	10.83
	Rg (Wm−2)					P (mm)
2008	62	189	256	107	153	106	279	245	223	899
2009	48	193	262	113	155	228	164	345	204	940
2010	41	174	257	93	143	292	306	305	382	1255
2011	49	209	265	112	160	100	94	193	166	576
2012	63	205	281	96	162	170	184	137	526	987
2013	43	178	270	98	148	310	207	197	355	1024
2014	45	181	246	86	140	247	325	527	603	1755
2015	51	194	268	96	152	215	230	288	325	991
2016	48	168	261	106	146	102	95	338	208	744
2017	57	209	286	98	164	130	85	99	501	927

). This influenced soil water content in the top 30 cm (SWC_{0–30 cm}), which was, for the best part of 2014, higher than the SWC at field capacity (Figure 3), implying that the soil water was in excess during the whole year.

3.2. NEE

Mean diurnal variation of NEE by months is shown in Figure 4. On a daily basis, a net carbon uptake (i.e., negative values of NEE) during the vegetation season usually started soon after sunrise, while the net release of CO₂ started about 30 min before the nightfall.

The exchange of CO₂ showed a clear seasonal pattern (Figure 4 and Figure 5). During dormant season NEE was moving around zero. Some smaller uptake could be occasionally noticed at midday also during the winter months, probably due to mosses, grasses and other herbaceous vegetation. Carbon uptake started in early spring with the development of grasses and the budding of leaves on trees. Noticeable CO₂ uptake started usually at the end of March. During May, carbon uptake increased strongly and reached a peak in June. However, there were also days when the forest lost carbon to the atmosphere during the vegetation season, which indicates that respiration was larger than the CO₂ uptake due to unfavourable meteorological conditions. With the start of leaf senescence in early autumn carbon uptake rapidly declined. During winter, when vegetation was in a state of dormancy, the absolute values of NEE reached their minimal values and ecosystem was a net source of carbon, which indicates active respiration even during the cold winter days. Growing season lasted mainly until the middle of October. The average (±SD) length of the growing season (GSL) was 195 (±13) days, where GSL was defined as the number of days between the first and the last day of the period when integrated three-day NEE was negative [40] (Table 1).

Ten years of EC measurements showed the high inter-annual variability of net ecosystem CO₂ exchange (Table 1, Figure 5). Annual sums of NEE (±95% confidence interval) range from −147 ± 13 gC m⁻² year⁻¹ in 2017 to −496 ± 15 gC m⁻² year⁻¹ in the year 2009, with an average of −319 gC m⁻² year⁻¹. Linear regression of annual NEE with years reveals a small positive, even though it was not significant (p = 0.106), trend in NEE (16.8 ± 18.3 gC m⁻² year⁻¹).

Net carbon balance of the vegetation is controlled, to a great extent, by the environmental variables. However, annual means of the environmental variables are rather poor predictors of annual NEE at Jastrebarsko forest. The matrix of cross correlation coefficients between main annual fluxes and the mean annual values of the environmental variables is given in

Table A6. Correlation coefficients for the main annual fluxes (NEE, GPP, R_ECO) and the mean annual values of the main driver variables (R_g—global radiation, T_a—air temperature, T_s—soil temperature, P—precipitation, SWC—soil water content, ΔSWC_May-Sep—difference of the mean monthly SWC in May and SWC in September, SGS and EGS—day of year for the start and the end of the growing season, respectively). Significant (p < 0.05) coefficients are marked with asterisk.

	NEE	GPP	RECO	Rg	Ta	Ts	P	SWC	ΔSWCMay–Sep	SGS	EGS
NEE	1
GPP	−0.578	1
RECO	0.327	0.582	1
Rg	0.107	−0.184	−0.107	1
Ta	0.209	−0.213	−0.038	0.105	1
Ts	−0.030	−0.108	−0.155	−0.176	0.902 *	1
P	0.255	−0.004	0.250	−0.644 *	0.328	0.398	1
SWC	0.421	−0.118	0.281	−0.778 *	0.191	0.272	0.776 *	1
ΔSWCMay–Sep	−0.707 *	0.403	−0.240	0.089	−0.219	−0.096	−0.645 *	−0.422	1
SGS	−0.021	0.221	0.238	−0.075	0.054	0.019	0.160	0.025	−0.228	1
EGS	−0.165	−0.189	−0.382	0.672 *	−0.458	−0.550	−0.749 *	−0.834 *	0.247	0.131	1

. Figure 6 provides an overview of the relationships (or lack thereof) between the NEE and the annual means of the main meteorological variables as well as soil water content. Correlations between NEE and average annual global radiation ($\overline{R_g}$), air temperature ($\overline{T_a}$), and precipitation ($\overline{P}$) are not statistically significant (Figure 6A–C). On the other hand, NEE showed strong negative correlation (r = −0.71, p = 0.02) with the difference (i.e., the decline) in soil water content from May to September (ΔSWC_May–Sep, Figure 6D). Here, we introduced a ΔSWC_May–Sep as a proxy variable that describes the dynamic of soil moisture. In the case of excess, or of a shortage of water in the soil during the growing season, the ΔSWC_May–Sep will be smaller. On the other hand, in the case of optimal meteorological and environmental conditions, there is neither a shortage of water, nor an excess of it. Hence, the depletion of the soil water can be expected and ΔSWC_May–Sep will be larger.

The multiple regression analysis of the relationship between the annual NEE and environmental variables requires addressing the issue of collinearity among predictors. One way of removing the correlation among predictor variables is by using Principal Components Analysis (PCA). We conducted PCA with a subset of environmental variables that correlate, at least slightly (|r| > 0.1, see Table A6), with the annual NEE. Results of the PCA (

Table A7. Principal Components Analysis (a and b) and the resulting correlation between NEE and Principal Components (c).

(a) Eigenvalues of the Correlation Matrix
Principal Component (PC)	Eigenvalue	Difference	Proportion of Variance	Cumulative Proportion
PC1	3.0422	1.3715	0.507	0.507
PC2	1.6706	0.7073	0.278	0.786
PC3	0.9633	0.7660	0.161	0.946
PC4	0.1973	0.0945	0.033	0.979
PC5	0.1028	0.0791	0.017	0.996
PC6	0.0237	-	0.004	1.000
(b) Relation between the Physical Driver Variables and the Principal Components
	Principal Components (Eigenvectors)
Variable	PC1	PC2	PC3	PC4	PC5	PC6	Unexplained
P	0.5313	−0.1376	−0.1502	−0.4509	0.6785	0.1121	0
Ta	0.3025	0.6464	0.0065	0.2386	0.1597	−0.6389	0
Ts	0.3513	0.5572	0.2789	−0.2349	−0.3487	0.5581	0
Rg	−0.3844	0.4151	−0.4819	0.2968	0.4023	0.4446	0
SWC	0.4925	−0.2747	0.1317	0.7710	0.0573	0.2587	0
ΔSWCMay–Sep	−0.3352	0.0694	0.8062	0.0468	0.4769	0.0565	0
(c) Correlation between the NEE and the Principal Components (PCs)
	NEE	PC1	PC2	PC3	PC4	PC5	PC6
NEE	1
PC1	0.3391	1
PC2	−0.0284	0	1
PC3	−0.6233 #	0	0	1
PC4	0.5962 #	0	0	0	1
PC5	−0.1664	0	0	0	0	1
PC6	−0.0355	0	0	0	0	0	1

^# Statistically significant at p < 0.1.

a,b) and the obtained correlation between NEE and Principal Components (PC, Table A7c), reveals that the first two components (PC1 and PC2) explain 78.6% of variance among the predictor variables, but they do not correlate with NEE (Table A7c). The PC3 and PC4 explain only 16.1% and 3.3% of the variance, but they have correlation coefficients with NEE of r = −0.6233 and r = 0.5962, respectively, although statistically significant only at p < 0.1 level. The PCs with low variances are usually considered to be unimportant and they are sometimes discarded from further analysis, but in our case, similarly to [70], PC3 and PC4 are very important. Comparison of the model of NEE with R_g and ΔSWC_May–Sep as predictors (adj. R² = 0.395, root-mean-square error (RMSE) = 73.3, p = 0.07) against the model of NEE with PC3 and PC4 as predictors (adj. R² = 0.671, RMSE = 54.0, p = 0.01) shows much better performance of the later one (Figure A3 and Figure A4).

3.3. GPP, R_ECO and NPP_EC

Carbon fluxes GPP, R_ECO and NPP_EC followed the seasonal pattern of NEE. Cumulative fluxes (NEE, GPP, R_ECO, NPP) during the period of the study are shown in the Appendix A (Figure A5). During the dormant season, the amount of carbon fixed by photosynthesis (GPP) had daily values that are close to zero (Figure 7). During the early spring, with the development of leaves, GPP started to rise rapidly and peaked in June. After a peak in early summer, GPP started to decline slowly and dropped again to small values after leaf fall. R_ECO followed the yearly pattern of air temperature (Figure 8). The highest daily values of R_ECO were achieved in the summer, while during the dormant season the daily values of R_ECO were significantly lower. The yearly pattern of NPP_EC followed the yearly pattern of GPP. The highest daily values of NPP were also achieved in June and the lowest during the dormant season (Figure 9). Negative NPP indicates that autotrophic respiration exceeded GPP. Such a situation occurred mainly during the cold part of the year. Annual sums of carbon fluxes are shown in Table 1. Both the annual GPP and annual R_ECO do not exhibit a statistically significant trend with time. Similarly to annual NEE, the annual values of GPP and R_ECO also do not correlate significantly with the average annual R_g, T_a and P. Furthermore, they do not correlate with ΔSWC_May–Sep either.

The highest annual GPP (1775 gC m⁻² year⁻¹) was recorded in the year 2015. The highest recorded annual R_ECO (1402 gC m⁻² year⁻¹) occurred in 2015 as a result of drought and high soil and air temperatures during the vegetation period. Despite a high radiated energy, the lowest recorded annual value of GPP (1384 gC m⁻² year⁻¹) occurred in 2017, while the lowest annual value of R_ECO (1117 gC m⁻² year⁻¹) occurred in the year 2008. NPP_EC was highest in the year 2009 (937 gC m⁻² year⁻¹), while the year 2017 had the lowest annual NPP (632 gC m⁻² year⁻¹). Temporal analysis showed a slightly negative (−10.7 gC m⁻² year⁻²) but not statistically significant (p = 0.303) trend in NPP_EC.

3.4. Stand Characteristics and NPP_BM

3.4.1. Stand Development

Temporal evolution of the forest stand structure during the study period 2008–2017, for all species combined, is shown in

Table 2. Temporal evolution of d100, h100 stand density, basal area and wood volume of the forest in the area around the eddy covariance tower (mean ± std. error).

DBH * (cm)	Year (at the End of the Season)
	2007		2008		2009		2010		2011		2012		2013		2014		2015		2016		2017
Stand age (years)
		35		36		37		38		39		40		41		42		43		44		45
d100—quadratic mean diameter of 100 thickest trees per hectare (cm)
		25.8		26.5		26.4		26.7		27.4		28.5		28.5		29.8		29.4		30.0		30.3
h100—top height ** (m)
		20.8		21.3		21.6		21.9		22.3		22.7		23.1		23.5		23.9		24.3		24.6
N—stand density (trees ha−1)
<10	946	±109	913	±108	888	±108	839	±104	808	±95	784	±93	782	±93	766	±94	729	±89	708	±86	686	±84
≥10	1085	±53	1118	±55	1114	±59	1065	±62	1054	±62	1070	±64	1070	±64	1046	±64	1023	±65	977	±64	969	±66
All	2031	±136	2031	±136	2002	±136	1903	±136	1862	±128	1854	±129	1852	±129	1812	±128	1752	±125	1686	±120	1655	±119
G—basal area (m2 ha−1)
<10	2.8	±0.4	2.6	±0.4	2.6	±0.4	2.4	±0.4	2.3	±0.3	2.2	±0.3	2.2	±0.3	2.2	±0.3	2.0	±0.3	1.9	±0.2	1.8	±0.2
≥10	22.5	±1.0	23.8	±1.0	24.3	±0.9	24.5	±1.0	25.4	±1.0	26.3	±1.1	27.3	±1.1	27.7	±1.1	27.9	±1.1	28.0	±1.2	28.3	±1.2
All	25.3	±1.0	26.5	±1.0	26.9	±1.0	26.9	±1.1	27.7	±1.1	28.5	±1.1	29.5	±1.2	29.9	±1.2	29.9	±1.2	29.9	±1.3	30.1	±1.3
V—wood *** volume (m3 ha−1)
<10	15	±2	14	±2	14	±2	13	±2	13	±2	12	±2	12	±2	12	±2	10	±2	10	±1	9	±1
≥10	214	±10	231	±10	240	±9	247	±11	261	±11	274	±11	290	±12	300	±12	306	±13	313	±14	321	±14
All	230	±10	245	±10	253	±10	260	±11	274	±11	286	±12	302	±12	311	±13	317	±13	323	±14	330	±15

* DBH—Stem diameter at 1.30 m above the ground. ** Top height is defined as the average height of 100 thickest trees per hectare [71]. *** Tree trunk and branches down to 3 cm diameter at thinner end, including bark.

. The results disaggregated according to tree species are given in the Appendix A (

Table A8. Temporal evolution of stand density, basal area and wood volume, according to tree species and size class, of the forest in the area around the eddy covariance tower (mean (std. error)).

Species	DBH	Year (at the End of the Season)
	(cm)	2007		2008			2009		2010		2011		2012		2013		2014		2015		2016		2017
N—stand density (trees ha−1)
Quercus robur	<10	119	±25	104	±24		100	±23	85	±22	81	±21	76	±20	74	±20	67	±18	61	±17	59	±17	52	±15
	≥10	502	±61	516	±63		505	±62	488	±61	486	±61	491	±61	491	±61	480	±57	464	±53	453	±51	447	±48
	All	620	±64	620	±64		605	±64	574	±60	568	±59	566	±59	564	±59	548	±56	525	±51	512	±50	499	±47
Fraxinus angustifolia	<10	91	±34	89	±33		81	±30	73	±28	68	±26	64	±24	64	±24	60	±22	56	±20	56	±20	54	±19
	≥10	90	±25	92	±26		100	±29	82	±25	84	±24	86	±25	86	±25	86	±25	84	±25	80	±24	80	±24
	All	181	±54	181	±54		181	±54	155	±46	153	±44	151	±43	151	±43	146	±41	140	±39	136	±38	134	±37
Carpinus betulus	<10	489	±97	474	±95		468	±94	457	±92	455	±92	447	±91	447	±91	444	±91	431	±88	417	±87	406	±86
	≥10	154	±80	169	±80		175	±81	175	±81	175	±81	183	±84	183	±84	185	±82	193	±82	197	±82	205	±81
	All	643	±145	643	±145		643	±145	633	±140	630	±139	630	±139	630	±139	628	±138	624	±136	614	±132	612	±130
Alnus glutinosa	<10	90	±46	90	±46		86	±45	74	±38	56	±27	52	±24	52	±24	50	±25	39	±20	35	±17	33	±15
	≥10	325	±67	325	±67		315	±68	305	±67	293	±66	293	±67	293	±67	278	±65	267	±63	232	±57	222	±56
	All	415	±99	415	±99		401	±99	378	±92	349	±84	345	±84	345	±84	328	±81	307	±75	268	±67	255	±65
Other broadleaf	<10	158	±52	156	±52		154	±53	149	±51	147	±51	145	±50	145	±50	145	±50	141	±49	141	±49	141	±49
	≥10	15	±6	17	±6		19	±8	15	±6	15	±6	17	±6	17	±6	17	±6	15	±7	15	±7	15	±7
	All	172	±53	172	±53		172	±53	164	±51	162	±51	162	±51	162	±51	162	±51	156	±51	156	±51	156	±51
All species	<10	946	±109	913	±108		888	±108	839	±104	808	±95	784	±93	782	±93	766	±94	729	±89	708	±86	686	±84
	≥10	1085	±53	1118	±55		1114	±59	1065	±62	1054	±62	1070	±64	1070	±64	1046	±64	1023	±65	977	±64	969	±66
	All	2031	±136	2031	±136		2002	±136	1903	±136	1862	±128	1854	±129	1852	±129	1812	±128	1752	±125	1686	±120	1655	±119
G—basal area (m2 ha−1)
Quercus robur	<10	0.5		±0.1	0.4	±0.1	0.4	±0.1	0.4	±0.1	0.4	±0.1	0.3	±0.1	0.3	±0.1	0.3	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.0
	≥10	11.7		±1.3	12.4	±1.4	12.7	±1.4	13.2	±1.4	13.9	±1.5	14.5	±1.5	15.3	±1.6	15.7	±1.6	15.9	±1.6	16.4	±1.7	16.6	±1.7
	All	12.2		±1.3	12.9	±1.4	13.1	±1.4	13.6	±1.4	14.3	±1.5	14.9	±1.5	15.6	±1.6	16.0	±1.6	16.2	±1.6	16.6	±1.7	16.8	±1.7
Fraxinus angustifolia	<10	0.3		±0.1	0.3	±0.1	0.3	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1
	≥10	2.2		±0.6	2.4	±0.6	2.6	±0.7	2.2	±0.6	2.3	±0.6	2.4	±0.6	2.5	±0.7	2.6	±0.7	2.7	±0.7	2.7	±0.7	2.7	±0.7
	All	2.5		±0.7	2.7	±0.7	2.9	±0.7	2.5	±0.6	2.6	±0.6	2.6	±0.7	2.8	±0.7	2.8	±0.7	2.9	±0.7	2.8	±0.7	2.9	±0.7
Carpinus betulus	<10	1.3		±0.3	1.3	±0.3	1.3	±0.3	1.3	±0.3	1.3	±0.3	1.3	±0.3	1.3	±0.3	1.3	±0.3	1.3	±0.3	1.2	±0.2	1.2	±0.2
	≥10	2.3		±1.4	2.5	±1.4	2.6	±1.5	2.7	±1.5	2.8	±1.6	3.0	±1.6	3.0	±1.7	3.1	±1.7	3.2	±1.7	3.3	±1.7	3.4	±1.7
	All	3.6		±1.5	3.8	±1.6	3.9	±1.6	4.0	±1.7	4.2	±1.7	4.3	±1.8	4.4	±1.8	4.4	±1.8	4.5	±1.8	4.6	±1.8	4.6	±1.8
Alnus glutinosa	<10	0.4		±0.2	0.4	±0.2	0.4	±0.2	0.3	±0.2	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.1	±0.1	0.1	±0.1	0.1	±0.0
	≥10	6.0		±1.2	6.2	±1.3	6.0	±1.3	6.1	±1.3	6.1	±1.3	6.1	±1.3	6.1	±1.4	6.0	±1.4	5.8	±1.3	5.4	±1.3	5.3	±1.3
	All	6.4		±1.3	6.6	±1.4	6.4	±1.4	6.4	±1.4	6.3	±1.4	6.3	±1.4	6.3	±1.4	6.1	±1.4	5.9	±1.4	5.5	±1.3	5.4	±1.3
Other broadleaf	<10	0.2		±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1
	≥10	0.2		±0.1	0.3	±0.1	0.3	±0.1	0.2	±0.1	0.3	±0.1	0.3	±0.1	0.3	±0.1	0.3	±0.1	0.2	±0.1	0.2	±0.1	0.2	±0.1
	All	0.5		±0.2	0.5	±0.2	0.5	±0.2	0.4	±0.1	0.4	±0.1	0.5	±0.1	0.5	±0.2	0.5	±0.2	0.4	±0.1	0.4	±0.1	0.4	±0.1
All species	<10	2.8		±0.4	2.6	±0.4	2.6	±0.4	2.4	±0.4	2.3	±0.3	2.2	±0.3	2.2	±0.3	2.2	±0.3	2.0	±0.3	1.9	±0.2	1.8	±0.2
	≥10	22.5		±1.0	23.8	±1.0	24.3	±0.9	24.5	±1.0	25.4	±1.0	26.3	±1.1	27.3	±1.1	27.7	±1.1	27.9	±1.1	28.0	±1.2	28.3	±1.2
	All	25.3		±1.0	26.5	±1.0	26.9	±1.0	26.9	±1.1	27.7	±1.1	28.5	±1.1	29.5	±1.2	29.9	±1.2	29.9	±1.2	29.9	±1.3	30.1	±1.3
V—wood ** volume (m3 ha−1)
Quercus robur	<10	3		±1	3	±1	3	±1	2	±1	2	±1	2	±1	2	±1	2	±0	2	±0	2	±0	1	±0
	≥10	115		±13	125	±14	130	±14	139	±15	149	±16	159	±17	171	±18	179	±18	185	±19	194	±20	200	±20
	All	118		±13	128	±14	133	±14	142	±15	152	±16	161	±17	173	±18	181	±18	187	±19	196	±20	201	±20
Fraxinus angustifolia	<10	2		±1	2	±1	1	±1	1	±1	1	±1	1	±0	1	±0	1	±0	1	±0	1	±0	1	±0
	≥10	19		±5	21	±6	24	±6	20	±5	22	±6	23	±6	25	±6	26	±7	27	±7	27	±7	28	±7
	All	21		±5	23	±6	25	±7	22	±5	23	±6	24	±6	26	±7	27	±7	28	±7	28	±7	29	±8
Carpinus betulus	<10	7		±1	7	±1	7	±1	7	±1	7	±1	7	±1	7	±1	7	±2	7	±1	7	±1	6	±1
	≥10	20		±13	22	±13	24	±14	26	±15	27	±16	29	±17	30	±17	31	±18	33	±19	35	±19	36	±19
	All	27		±13	29	±14	31	±15	33	±16	34	±17	36	±17	37	±18	39	±19	40	±19	42	±20	42	±20
Alnus glutinosa	<10	3		±2	3	±2	3	±2	2	±1	2	±1	1	±1	1	±1	1	±1	1	±1	1	±0	1	±0
	≥10	58		±12	60	±12	59	±12	60	±13	61	±13	61	±13	62	±14	61	±14	60	±14	56	±14	56	±14
	All	60		±13	63	±13	61	±13	62	±14	62	±14	62	±14	63	±14	62	±14	61	±14	56	±14	56	±14
Other broadleaf	<10	1		±0	0	±0	0	±0	0	±0	0	±0	0	±0	0	±0	0	±0	0	±0	0	±0	0	±0
	≥10	2		±1	2	±1	3	±1	2	±1	2	±1	2	±1	2	±1	2	±1	1	±1	1	±1	1	±1
	All	3		±1	3	±1	3	±1	2	±1	2	±1	2	±1	2	±1	2	±1	1	±1	1	±1	1	±1
All species	<10	15		±2	14	±2	14	±2	13	±2	13	±2	12	±2	12	±2	12	±2	10	±2	10	±1	9	±1
	≥10	214		±10	231	±10	240	±9	247	±11	261	±11	274	±11	290	±12	300	±12	306	±13	313	±14	321	±14
	All	230		±10	245	±10	253	±10	260	±11	274	±11	286	±12	302	±12	311	±13	317	±13	323	±14	330	±15

** Tree trunk and branches down to 3 cm diameter at the thinner end.

). Changes of trees’ DBH distribution with time, according to tree species (Figure A6), and species-specific, age-dependent tree height curve parameters (Table A2) are given in the Appendix A.

3.4.2. NPP_WBt by Tree Species

Most important part of the stand production, at least from the perspective of timber production, comes from the formation of new xylem in tree stems and branches. While it is relatively straightforward to assess the contribution of a given tree species to the part of the NPP that is attributed to aboveground woody growth, the situation is much more complex for non-woody plant components. Therefore, we present only the results of the NPP_WBt (total woody above- and below-ground) partitioned according to main tree species (Figure 10), while the annual production of leaf litter (NPP_L) and trees’ fruits (NPP_F) is given in Table 1.

The NPP_WBt ranged from 299 gC m⁻² year⁻¹ (2017) to 567 gC m⁻² year⁻¹ (2010), with the mean of 484 gC m⁻² year⁻¹ during the period 2008–2017. Disaggregation with respect to tree species clearly shows that Quercus robur is the key species (Figure 10A) with an average contribution to NPP_WBt of 66.8%. The average contributions of other important tree species are much smaller; namely, Alnus glutinosa 10.5%, Carpinus betulus 11.8%, Fraxinus angustifolia 9.9%. Interestingly, the trend of total productivity of woody biomass with time is significantly negative (−17.8 gC m⁻² year⁻² or −3.7% year⁻¹, p = 0.030). It is negative for all main species except for Quercus robur, for which the trend is also slightly negative but not significant (−1.8% year⁻¹, p = 0.342) (Figure 10B). The most pronounced negative trend was observed in Fraxinus angustifolia (−11.6% year⁻¹, p < 0.001), followed by Alnus glutinosa (−8.4% year⁻¹, p < 0.001) and Carpinus betulus (−3.9% year⁻¹, p = 0.037). The year 2017, with unfavourable meteorological conditions, is the last year in the used dataset and it has a strong effect on this trend. Without data from 2017, a statistically significant (95% confidence), negative trend in NPP_WBt is observed only for Fraxinus angustifolia (−11.1% year⁻¹, p < 0.001) and Alnus glutinosa (−7.5% year⁻¹, p = 0.006), while other tree species, as well as total NPP_WBt, do not display statistically significant trends.

Seasonal dynamics of NPP_WBt was reconstructed by fitting Equation (10) for each year separately and obtained estimates for the parameters NPP_WBt,annual, DOY_max__WBt, and k are given in Appendix (Table A4). Recalling Equations (11) and (13), the measured annual NPP_WBt and the existing estimates of k, maximum daily net primary production (dNPP_max) was calculated for each year. The dNPP_max ranged from 6.1 gC m⁻² day⁻¹ in 2008 to 2.3 gC m⁻² day⁻¹ in 2017 and it has exhibited a steady downward trend (−0.37 gC m⁻² day⁻¹ year⁻¹, p < 0.001) during the period of the research. DOY_max__WBt ranged from 143 in 2012 to 165 in 2017 with the average of 153 and it has not exhibited a significant trend with time.

3.4.3. NPP_BM

The NPP_BM ranged from 483 gC m⁻² year⁻¹ (2017) to 783 gC m⁻² year⁻¹ (2010), with the mean of 680 gC m⁻² year⁻¹ (Table 1). Although neither the annual NPP_L, nor the NPP_F trends were significant (p = 0.309 and p = 0.501, respectively), the overall trend of NPP_BM was significantly negative (−18.7 gC m⁻² year⁻² or −2.7% year⁻¹, p = 0.046) due to a major contribution of NPP_WBt. The data for the last year in the dataset, namely the year 2017 with unfavourable meteorological conditions, strongly influenced the trend. Without data from 2017, NPP_BM does not show a statistically significant trend (p = 0.239). On average, NPP_LF contributed 29.1% of NPP_BM, out of which NPP_F had a relatively small contribution of 1.8%, except in 2015, a typical oak mast year, when the share of NPP_F in NPP_BM was outstandingly high 9.5%.

Seasonal dynamics of NPP_LF were estimated assuming that the shape parameter k for the NPP_LF is the same as for NPP_WBt (Table A4). The DOY_{max_LF}, as calculated using Equation (14) and parameters from Table A4 (see Section 2.5.2), ranged from 175 in 2011 to 206 in 2017, while the average was 190.

3.5. Comparison of NPP Estimates from Eddy Covariance and Biometric Measurements

In order to assess the agreement between the estimates of NPP from eddy covariance and from biometric measurements, the annual NPP_BM was regressed on NPP_EC (Figure 11). Although the correlation is statistically significant (p = 0.03), the measured NPP_BM in all years was lower than NPP_EC (slope of the regression: 0.67). Paired samples t-test revealed that NPP_BM is significantly (p < 0.001) smaller than NPP_EC. A comparison of the seasonal dynamics of NPP_EC with NPP_BM and its components is shown in Figure 12. In each year, NPP_BM started earlier than NPP_EC in the spring and ended earlier in the autumn. NPP_EC shows negative values at the beginning of the year until the DOY 94 (on average), which is result of the calculation method (Equation (7)). During vegetation season, NPP_EC exceeds NPP_BM approximately after DOY 203 and continues to increase until it reaches its maximum at, depending on the year, between DOY 295 and 311, On the other hand, NPP_BM ends earlier in the autumn, achieving its 95% in average around the DOY 254.

4. Discussion

4.1. Variability of CO₂ Fluxes from EC

The observed seasonal variability of carbon fluxes at Jastrebarsko forest shows a pattern typical for the temperate broadleaved forest in the northern hemisphere [5,40]. Annual sums of carbon fluxes are within the range of values for temperate deciduous broadleaved forests published in a recent review by Baldocchi et al. [7] (

Table 3. Carbon fluxes in temperate deciduous broadleaved forests (average ± standard deviation).

Country, Site	Genus	Period	NEE	GPP	RECO	Reference
			gC m−2 year−1
Denmark, Sorø	Fagus	1996–2009	−156 ± 103	1727 ± 136	1570 ± 97	[72]
France, Hesse	Fagus	1995–2005	−386 ± 171	1397 ± 192	1011 ± 138	[6]
Germany, Hainich	Fagus, Fraxinus	2003–2012	−483 ± 70	1498 ± 83	1015 ± 51	[73]
UK, Straits Inclosure	Quercus	1999–2010	−486 ± 115	2034 ± 228	1548 ± 192	[40]
US, Duke Forest	Quercus, Carya	2001–2008	−402 ± 96	1982 ± 300	1580 ± 237	[74]
US, Harvard Forest	Acer, Quercus	1992–2004	−242 ± 100	1400 ± 164	1153 ± 105	[26]
Croatia, Jastrebarsko	Quercus	2008–2017	−319 ± 94	1594 ± 109	1275 ± 94	this study

). If we make a comparison at a site-to-site level, we notice differences in fluxes. The site Denmark, Sorø [72] has smaller NEE but larger GPP and R_ECO, while the site France, Hesse [6] has similar NEE but smaller GPP and R_ECO. Both sites are beech forests with distinction that Sorø was ~80 years old, while France, Hesse is ~40 years old, similar to Jastrebarsko. A higher NEE at Hainich site [73], with respect to Jastrebarsko site, might be related to the difference in forest structure. Hainich forest is characterised as an old-growth forest (up to 250 years) with a highly diverse horizontal and vertical structure [73]. Also, while average GPPs are similar, higher NEE at Hainich might be due to lower R_ECO as a result of significantly lower average air temperature at Hainich site (8 °C) as compared to Jastrebarsko (11.2 °C). There are two sites (the UK, Straits Inclosure [40] and US, Duke Forest [74]) where all fluxes are noticeably higher than the ones that were observed at Jastrebarsko site. The UK, Straits Inclosure site [40] is an 80-years old oak plantation in a climate that typically does not exhibit summer drought, which might contribute to the higher annual carbon accumulation than our natural forest in a climate with hot and, comparatively, dry summers. Nevertheless, relatively large fluxes observed at UK, Straits Inclosure site are still interesting and more thorough comparison would be needed. On the other hand, fluxes observed at the US, Duke Forest might not be directly comparable due to some differences in methodology, namely gap-filling and partitioning [74,75]. Finally, all fluxes at US, Harvard Forests [26] are somewhat smaller but not far from fluxes at Jastrebarsko forest. The aim of this simple comparison with other sites was not to explain the reported differences in fluxes, but to emphasize the multitude of possible issues and the very limited number of comparable sites with long time-series, which might help in explaining the observed variability.

Carbon fluxes are strongly driven by environmental variables [76,77]. Their inter-annual variability is found to be influenced by meteorological conditions, e.g., light [21,22,25,40,78,79], temperature [78,79,80], length of growing season [77,79], soil water status [6], this study, as well as with stand characteristics, e.g., successional change in forest composition [26]. However, our results indicate that, at the Jastrebarsko site, the average annual values are poor predictors of annual carbon fluxes (Figure 6A–C).

In the present study, we found a significant linear relationship only between NEE and ΔSWC_May–Sep (Figure 6D). The low NEE values that were observed in 2012 and 2017 were clearly the result of prolonged dry conditions from previous 2011 and 2016, represented by low ΔSWC_May–Sep values (Figure 6D). On the contrary, the low NEE that was observed in 2014 was probably due to cloudy weather with low mean annual R_g (Figure 6A). In 2014, the low ΔSWC_May–Sep did not imply dry, but extremely wet conditions. In fact, SWC did not change much from May to September due to excessive rain, which caused SWC to be high throughout the year (Figure 3 bottom, and Figure 6C). The highest NEE values that were recorded in 2009 and 2015 were due to favourable meteorological conditions with sufficient water and optimal light during the growing season. The relationship between NEE, ΔSWC_May–Sep, and light is given in a 3D scatter plot (Figure A3). The poor correlation that was observed between NEE and mean annual air temperature and precipitation (Figure 6B,C) is perhaps most noticeable if we look at the years 2008 and 2017. The year 2017 was, on average, 0.14 °C warmer with only 12 mm less precipitation compared to 2008. However, the difference in annual NEE was 205 gC m⁻² year⁻¹. Although 2008 and 2017 did not differ much with respect to mean annual air temperature and precipitation, the difference in the pattern of precipitation is striking. There is almost a complete lack of precipitation during the spring and summer of 2017 followed by a sudden cooling and excess precipitation in September 2017 (Figure 3). Similarly to findings of Zscheischler et al. [81], our results also indicate that the intrayear weather patterns and short-term favourable weather conditions are more important than the annual means in control of interannual variability of carbon fluxes.

4.2. Dynamic of NPP Estimated with Biometric Method

In his definition of ecological succession, Odum [82], among others, stated that the succession “culminates in a stabilized ecosystem in which maximum biomass (or high information content) and symbiotic function between organisms are maintained per unit of available energy flow”. Characteristic of NEP dynamic throughout the succession, according to Odum [82,83], is a rapid increase during the early stage of succession that reaches maximum during the mid-succession, followed by decline, and eventually reaching zero NEP at the late stage of succession. However, a recent publication by Curtis and Gough [83] question Odum’s hypothesis of steep decline and zero NEP at late succession in the case of temperate broadleaf forests. It appears that the wood NPP in some forests, transitioning from early to middle succession, exhibits resilience that is attributed to changes in species composition [84].

In temperate forests of Western and Central European lowlands, the pedunculate oak is a typical species of the final stage in forest succession [85]. Forest management in Croatia is considered to be close to nature, striving at securing environmental and economic functions of forests [34]. In other words, current, close-to-nature management of pedunculate oak forests (tending, thinning, and regeneration cuts) is aiming at continuously maintaining pedunculate oak as the main species. The absence of regeneration cuts, e.g., due to conservation of old pedunculate oak stands, may lead to ecosystem regression and can result in forest-like stand of Crataegus monogyna [86].

Jastrebarsko forest at the EC site is in the middle successional stage. Woody biomass NPP (NPP_WBt) showed a negative trend during the investigated period and it could be an indicator of the declining NPP, which would be in line with the theoretical framework of Odum [82]. However, the observed trend is also affected by the unfavorable meteorological conditions in 2017 (the last year of our study), which had a strong negative effect on NPP and consequently impacted the trend. Looking into other possible reasons, we should consider reports pointing out that decreasing soil nutrient availability and increasing stomatal limitation, leading to reduced photosynthetic rates, are the two main causes of NPP decline in forest stands [87]. On the other hand, a recent study in an oak forest [88] pointed out that tree mortality was the primary cause of the decline in aboveground biomass accumulation with stand age, and not a reduction in individual tree growth.

Excluding data from 2017 and performing further analysis with respect to tree species (Section 3.4.2) revealed that only Fraxinus angustifolia and Alnus glutinosa have a statistically significant decrease in NPP_WBt. In the last several years, the spread of Hymenoscyphus fraxineus fungi (more widely known by its older name Chalara fraxinea) played an important part in the decline of Fraxinus angustifolia [89]. On the other hand, the decline in Alnus glutinosa is more likely age-related, due to its shorter maximum lifespan of 100–160 years [90]. For example, the typical period of rotation for stands dominated by Alnus glutinosa in Croatia is 70 years, as compared to 140 years long rotation for Quercus robur stands. Hence, the increased mortality of Alnus glutinosa and Fraxinus angustifolia probably contributed to the observed trend of NPP_BM. This can be corroborated if we look at the evolution of species-specific densities at Jastrebarsko forest (Table A8). We can note that the Alnus glutinosa had the lowest survival rate of 61.4%, when compared to the average survival rate of 81.5%.

While a decreased share of Fraxinus angustifolia probably did not affect nutrient availability, the reduction in the share of Alnus glutinosa, which is a nitrogen-fixing species, could have had a negative effect on nutrient availability and NPP_BM trend. A model analysis showed that the effects of soil denitrification (during flood periods) and N-fixation by Alnus glutinosa have a significant role in Jastrebarsko forest [10].

4.3. Comparison of Biometric and Eddy Covariance NPP Estimate

In studies that are similar to this one, it has been common to compare NEP estimates from EC and BM methods [8,18,19,20,21,22,32,91,92,93], because such estimates are methodologically different and have independent measurement errors [19]. The biometric method of the NEP assessment requires the assessment of R_h. The R_h is usually assessed from soil respiration measurements, using methods such as root exclusion by trenching, tree girdling, or isotope labelling [94,95], or assuming, e.g., based on the work of Hanson et al. [94], that R_h is usually 50% of the soil respiration flux [19,38]. Alternatively, R_h can be calculated from the R_ECO if the share of R_h (or R_a) in R_ECO is known [6]. In our case, R_h is estimated from R_ECO and a fixed R_h:R_ECO ratio, where R_ECO is obtained with EC method by the partitioning of NEE fluxes. This rendered that the comparison of NEPs would no longer be valid. Therefore, in order to preserve the independence of two measurements, we had to compare estimates of NPPs. A similar approach of comparing NPPs was already used in other studies [6,22,23].

Comparison of annual NPP_EC and NPP_BM showed good correlation (R = 0.680, adj. R² = 0.396, p = 0.03). However, the observed correlation is obtained only when including the last year of the experiment (year 2017) when unfavorable meteorological conditions negatively affected both estimates. Our results are in line with previous studies stating that only long-term monitoring of C fluxes can reveal a correlation between biometric and EC estimates [18,22]. EC estimate was always higher than the biometric one (Figure 12), with an average difference of 140 gC m⁻² year⁻¹, i.e., 17%. One of the possible reasons for the observed bias is a fact that NPP_EC includes the production of whole ecosystem (trees, understory and ground vegetation) with all of its above- and belowground carbon pools, while NPP_BM (in our study) includes only the production of trees, excluding fine roots and grasses. Consequently, we can assume that the higher NPP_EC, as compared to NPP_BM, is due to the estimate from which productions of fine roots and grasses are missing. Productions of fine roots and grasses can be significant contributors to total biometric NPP estimate, ranging from 10% [6] to 41% [22] for fine roots or 45% for fine roots and grasses together [96]. Underestimation of biometric NPP has been identified as key potential source of bias in similar studies [93].

One of the main assumption, which we made for the purpose of NPP_EC calculation, is that the ratio of R_h:R_ECO is fixed in time because R_h estimates were available only for the years 2008 and 2009 [38]. The fixed R_h:R_ECO ratio assumption is questionable, and probably not true [97,98], but we think that it is not fully unjustified. Namely, the forest in the footprint was thinned in 2006 and 2007 with an average intensity of 8% by volume (see caption of Figure A1). Thinning residues constituted a one-time addition to the litter pool probably causing slight increase in R_h in several years following thinning. In addition, Bond-Lamberty et al. [98] recently provided evidence that global R_h is rising due to shifts in soil organic carbon (SOC) forms and enhanced SOC mineralization driven by rising global temperatures. At the same time, according to Mori et al. [99], continuous accumulation of living biomass during the study period should have resulted with an increase in R_a. When considering both studies [98,99], R_ECO should have experience an increase with time. At our site, R_ECO exhibited a small, positive, but not statistically significant trend (p = 0.140). All of these facts do not point unanimously toward conclusion that there should be a significant trend in R_h:R_ECO during the study period. Therefore, in the absence of sufficient evidence, we considered that the use of the constant R_h:R_ECO ratio is the only justifiable approach. Nevertheless, further research at Jastrebarsko site is needed on the partitioning of the ecosystem respiration into autotrophic and heterotrophic in order to improve the estimates of NPP from flux measurements, which could also be beneficial for the NPP estimates from remote sensing [100,101,102].

Before comparing the seasonal dynamic of NPP from two independent estimates (Figure 12), it is important to recall the processes which drive each estimate. EC NPP estimate is primarily driven by canopy photosynthesis [4], therefore it reflects current accumulation of atmospheric C that can further be partitioned into structural growth and labile C storage [23]. On the other hand, biometric NPP estimate represents biomass growth that uses carbohydrates from current assimilation as well as previously stored NSC [103]. NSC is used in spring for the growth of leaves, roots and even wood in deciduous broadleaf trees [104]. e.g., Barbaroux et al. [105] reported that mean quantities of NSC, mobilized between October of the previous year and June of the current year, reached 2.8% and 1.2% of total carbon biomass in sampled, approx. 40-year-old, trees of Quercus petrea and Fagus sylvatica, respectively.

Figure 12 reveals clear difference in seasonal NPP dynamic from two independent estimates. In the spring, NPP_BM starts before NPP_EC on the account of C reserves stored in previous years [78,106]. Later in the vegetation season, stem growth slows down and even ceases, but the forest ecosystem continued to absorb carbon, most likely in the NSC pool [6,19,94].

NSC represents very important storage/reserve pool of carbon with distinct seasonal pattern [23]. To account for seasonal changes in this carbon pool we used a simple modelling approach (see Section 2.5.2) aimed at addressing the replenishing of NSC pool in late summer and autumn. NPP estimated in this way greatly improved the agreement of NPP_BM with NPP_EC (Figure 12). We also note that the NPP_LF, with its average contribution of 29.1% to the total NPP_BM, is close to findings of Gough et al. [22] which reported that >25% of net annual photosynthetic C assimilation occurred after growth had stopped in the autumn. Nevertheless, we are aware of the limitations of our approach due to the complexity of the NSC and wood formation dynamics [69,78,107]. Further development of used model is highly needed as the understanding and modelling of seasonal carbon allocation to the reserve pool is still being investigated [16,67,68].

5. Conclusions

Eddy covariance measurements of CO₂ flux over Jastrebarsko pedunculate oak forest during the 2008–2017 period show that the forest was carbon sink in every year. NEE exhibited high interannual variability, but it was poorly correlated with the average annual R_g, T_a, and total annual P. On the other hand, NEE did correlate significantly with the May-to-September decrease in soil water content. Independent estimation of NPP with eddy covariance and biometric method showed good overall agreement (R² = 0.463, p = 0.03), although in all years NPP_BM was lower than NPP_EC. This could be because the NPP of fine roots and grasses were not measured as part of NPP_BM, or due to the use of constant R_h:R_ECO ratio in calculation the of NPP_EC. The use of simple theoretical model for the replenishment of carbon reserves in the late season greatly improved the seasonal agreement of NPP_BM and NPP_EC estimates. Further research is needed on the contribution to NPP of fine roots and understory vegetation, partitioning of the ecosystem respiration into autotrophic and heterotrophic, as well as on the non-structural carbohydrates dynamics.