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

We analysed 10 years (2008–2017) of continuous eddy covariance (EC) CO2 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 (RECO) using the online tool by Max Planck Institute for Biogeochemistry in Jena, Germany. Annual NEE, GPP, and RECO were correlated with main environmental drivers. Net primary productivity was estimated from EC (NPPEC), as a sum of −NEE and Rh obtained using a constant Rh:RECO ratio, and from independent periodic biometric measurements (NPPBM). 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−2 year−1, GPP of 1594 (±109) gC m−2 year−1, and RECO of 1275 (±94) gC m−2 year−1. Annual NEE showed high inter-annual variability and poor correlation with annual average global radiation, air temperature, and total precipitation, but significant (R2 = 0.501, p = 0.02) correlation with the change in soil water content between May and September. Comparison of annual NPPEC and NPPBM showed a good overall agreement (R2 = 0.463, p = 0.03), although in all years NPPBM was lower than NPPEC, with averages of 680 (±88) gC m−2 year−1 and 819 (±89) gC m−2 year−1, respectively. Lower values of NPPBM 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.


Introduction
Currently global forests store approximately 30% of total anthropogenic CO 2 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 2 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 2 fluxes from EC with respect to observed meteorology, (2)

estimate
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 3 (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 H2O 4.9-5.7),but it becomes neutral (pH in H2O 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- 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 2 O 4.9-5.7),but it becomes neutral (pH in H 2 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 −1 , 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 2 per plot, which is equivalent to sampling intensity of 2% by area.

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 2 storage term was measured with a profiler system for the CO 2 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 2 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 2 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.

EC Data Processing
Raw EC data were stored in one-hourly binary files (Table A1).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 2 and a difference ('delta') between two consecutive CO 2 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 2 flux was calculated as the average of the product of instantaneous fluctuations of vertical wind velocity (w') and CO 2 concentration fluctuations (s'): 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 2 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 −2 s −1 > NEE > 30 µmoL m −2 s −1 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 −2 ) 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): where M d denotes the median of the differences and value d i was calculated for each half-hourly NEE as: Half-hourly NEE was flagged as spike and removed from further processing if: 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 2 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 −1 with increment of 0.05 m s −1 ) 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.

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 : where T ref marks reference temperature that is set to 15 • C, R ECO,ref denotes respiration at the reference temperature, T 0 is set to −46.02 • C, T marks measured air or soil temperature, while the activation energy parameter, E 0 , 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: 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: Forests 2018, 9, 764 7 of 35

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).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.

Biometric Measurements and Estimation of NPPBM
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 (NPPWBt) by components: stem and branches (NPPSB), twigs (NPPT), and roots (NPPR).NPP of leaves and fruits (NPPLF) was estimated from litterfall measurements.The sum of NPPWBt and NPPLF is denoted as NPPBM.Net primary production of a given component is expressed in units of gC m −2 , while net primary productivity is expressed in units of gC m −2 year −1 .More details on measurements and calculations are given below.

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 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.

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 −2 , while net primary productivity is expressed in units of gC m −2 year −1 .More details on measurements and calculations are given below.From 2007 until 2017, at the end of growing season, the DBH of all trees (DBH > 2 cm) in each of the 24 plots (Figures 2 and 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).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) [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 −3 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: stem and branches : The NPP was calculated from changes in C i at a given plot.For example, the NPP WBt is calculated as: where N plots is the number of plots, A is the plot area, t 1 and t 2 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 ≤ t j ≤ 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 ).

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 2 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: 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 ): If r max is known (or estimated) the parameter k can be calculated using the following equation: Inversely, 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).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: Finally, 3. Results

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).
The coldest months were December, January, and February, with the average monthly temperatures during the study period of 1.68 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).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.

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 2 started about 30 min before the nightfall.
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 −2 ), mean annual soil water content of the 0-30 cm soil layer (SWC, vol H2O vol −1 ), difference between mean soil water content in May and in September (SWC May -SWC Sep , vol H2O vol −1 ), 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 −2 year −1 ), and autotrophic and heterotrophic respiration (R a and R h , respectively, gC m −2 year −1 ).NPP EC and NPP BM are two independent estimates of the annual net primary productivity (NPP, gC m −2 year −1 ).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.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).This influenced soil water content in the top 30 cm (SWC0-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.The exchange of CO2 showed a clear seasonal pattern (Figures 4 and 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 CO2 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 CO2 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).The exchange of CO 2 showed a clear seasonal pattern (Figures 4 and 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 2 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 2 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 CO2 exchange (Table 1, Figure 5).Annual sums of NEE (±95% confidence interval) range from −147 ± 13 gC m −2 year −1 in 2017 to −496 ± 15 gC m −2 year −1 in the year 2009, with an average of −319 gC m −2 year −1 .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 −2 year −1 ).
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. 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 (R g ), air temperature ( ), and precipitation (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 (ΔSWCMay-Sep, Figure 6D).Here, we introduced a ΔSWCMay-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 ΔSWCMay-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 ΔSWCMay-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 A7a,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.6233and r = 0.5962, respectively, although statistically significant only Ten years of EC measurements showed the high inter-annual variability of net ecosystem CO 2 exchange (Table 1, Figure 5).Annual sums of NEE (±95% confidence interval) range from −147 ± 13 gC m −2 year −1 in 2017 to −496 ± 15 gC m −2 year −1 in the year 2009, with an average of −319 gC m −2 year −1 .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 −2 year −1 ).
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. 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 (R g ), air temperature (T a ), and precipitation (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 A7a,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.6233and 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 2 = 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 2 = 0.671, RMSE = 54.0,p = 0.01) shows much better performance of the later one (Figures A3 and A4).

GPP, RECO and NPPEC
Carbon fluxes GPP, RECO and NPPEC followed the seasonal pattern of NEE.Cumulative fluxes (NEE, GPP, RECO, 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.RECO followed the yearly pattern of air temperature (Figure 8).The highest daily values of RECO were achieved in the summer, while during the dormant season the daily values of RECO were significantly lower.The yearly pattern of NPPEC 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, 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.
Sep either.
The highest annual GPP (1775 gC m −2 year −1 ) was recorded in the year 2015.The highest recorded annual RECO (1402 gC m −2 year −1 ) 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 −2 year −1 ) occurred in 2017, while the lowest annual value of RECO (1117 gC m −2 year −1 ) occurred in the year 2008.NPPEC was highest in the year 2009 (937 gC m −2 year −1 ), while the year 2017 had the lowest annual NPP (632 gC m −2 year −1 ).Temporal analysis showed a slightly negative (−10.7 gC m −2 year −2 ) but not statistically significant (p = 0.303) trend in NPPEC.2. The results disaggregated according to tree species are given in the Appendix A (Table A8).Changes of trees' DBH distribution with time, according to tree species    2. The results disaggregated according to tree species are given in the Appendix A (Table A8).Changes of trees' DBH distribution with time, according to tree species  The highest annual GPP (1775 gC m −2 year −1 ) was recorded in the year 2015.The highest recorded annual R ECO (1402 gC m −2 year −1 ) 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 −2 year −1 ) occurred in 2017, while the lowest annual value of R ECO (1117 gC m −2 year −1 ) occurred in the year 2008.NPP EC was highest in the year 2009 (937 gC m −2 year −1 ), while the year 2017 had the lowest annual NPP (632 gC m −2 year −1 ).Temporal analysis showed a slightly negative (−10.7 gC m −2 year −2 ) but not statistically significant (p = 0.303) trend in NPP EC .

Stand Development
Temporal evolution of the forest stand structure during the study period 2008-2017, for all species combined, is shown in Table 2.The results disaggregated according to tree species are given in the Appendix A (Table A8).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. * 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.

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 −2 year −1 (2017) to 567 gC m −2 year −1 (2010), with the mean of 484 gC m −2 year −1 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 −2 year −2 or −3.7% year −1 , 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 −1 , p = 0.342) (Figure 10B).The most pronounced negative trend was observed in Fraxinus angustifolia (−11.6% year −1 , p < 0.001), followed by Alnus glutinosa (−8.4% year −1 , p < 0.001) and Carpinus betulus (−3.9% year −1 , 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 −1 , p < 0.001) and Alnus glutinosa (−7.5% year −1 , p = 0.006), while other tree species, as well as total NPP WBt , do not display statistically significant trends.
G-basal area (m ha ) <10 2.8 ±0. 4  * 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.

NPPWBt 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 NPPWBt (total woody above-and belowground) partitioned according to main tree species (Figure 10), while the annual production of leaf (NPPL) and trees' fruits (NPPF) is given in Table 1.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 A (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 −2 day −1 in 2008 to 2.3 gC m −2 day −1 in 2017 and it has exhibited a steady downward trend (−0.37 gC m −2 day −1 year −1 , 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.

NPP BM
The NPP BM ranged from 483 gC m −2 year −1 (2017) to 783 gC m −2 year −1 (2010), with the mean of 680 gC m −2 year −1 (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 −2 year −2 or −2.7% year −1 , 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.

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 < Forests 2018, 9, 764 18 of 35 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.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).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 RECO, while the

Variability of CO 2 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).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 −2 year −1 .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.

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].

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 2 = 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 −2 year −1 , 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].

Conclusions
Eddy covariance measurements of CO 2 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 2 = 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.* Unites: DBH (cm), h (m), V (m 3 ).** 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.

Figure 1 .
Figure 1.Map (left) shows a distribution of pedunculate oak stands in Croatia (marked in green) with the enlarged areas of the Pokupski basin and the location of the study site with the installed eddy covariance tower (right).

Figure 1 .
Figure 1.Map (left) shows a distribution of pedunculate oak stands in Croatia (marked in green) with the enlarged areas of the Pokupski basin and the location of the study site with the installed eddy covariance tower (right).

Figure 2 .
Figure 2. Flux footprint; (a) before the elevation of eddy covariance (EC) tower (23 m); (b) after the elevation of EC tower (27 m); little cyan circles denote the location of 24 circular plots.

Figure 2 .
Figure 2. Flux footprint; (a) before the elevation of eddy covariance (EC) tower (23 m); (b) after the elevation of EC tower (27 m); little cyan circles denote the location of 24 circular plots.

Forests 1 .
NPP of Total Woody Biomass, Leaves and Fruits

Figure 3 .
Figure 3. Top panel: Average daily air (Ta) and soil (Ts) temperatures; bottom panel: daily precipitation (P, left y-axis) and soil water content (SWC, right y-axis) in the top 30 cm soil layer.Horizontal lines indicate the values of SWC at wilting point (dotted line), at field capacity (short-dashed line) and at saturation point (long-dashed line) for the location of the EC site [37].

Forests
2018, 9, x FOR PEER REVIEW 13 of 403.2.NEEMean diurnal variation of NEE by months is shown in Figure4.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 CO2 started about 30 min before the nightfall.

Figure 4 .
Figure 4. Mean diurnal variation of NEE by months during the years 2008 to 2017; blue dashed lines represent standard deviation among years; vertical dashed lines mark the time of sunrise and sunset on the 15 day in a month.

Figure 4 .
Figure 4. Mean diurnal variation of NEE by months during the years 2008 to 2017; blue dashed lines represent standard deviation among years; vertical dashed lines mark the time of sunrise and sunset on the 15 day in a month.

Forests 2018, 9 ,
x FOR PEER REVIEW 15 of 40 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 Rg and ΔSWCMay-Sep as predictors (adj.R 2 = 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 2 = 0.671, RMSE = 54.0,p = 0.01) shows much better performance of the later one (Figures A3 and A4).

Figure 6 .
Figure 6.Annual NEE with respect to main meteorological variables: (A) mean annual global radiation; (B) mean annual air temperature; (C) annual precipitation; and, (D) ΔSWCMay-Sep, the difference in the mean monthly soil water content between the months of May and September (n.s.not significant, CI-confidence interval).

Figure 6 .
Figure 6.Annual NEE with respect to main meteorological variables: (A) mean annual global radiation; (B) mean annual air temperature; (C) annual precipitation; and, (D) ∆SWC May-Sep , the difference in the mean monthly soil water content between the months of May and September (n.s.-not significant, CI-confidence interval).

Figure 9 .
Figure 9. Daily values of NPPEC in Jastrebarsko forest during 2008-2017.3.4.Stand Characteristics and NPPBM 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 Table2.The results disaggregated according to tree species are given in the Appendix A (TableA8).Changes of trees' DBH distribution with time, according to tree species

Figure 9 .
Figure 9. Daily values of NPPEC in Jastrebarsko forest during 2008-2017.3.4.Stand Characteristics and NPPBM 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 Table2.The results disaggregated according to tree species are given in the Appendix A (TableA8).Changes of trees' DBH distribution with time, according to tree species

Figure 10 .
Figure 10.NPP WBt of total woody biomass (above-and below-ground) according to tree species (A); anomaly (relative deviation from the mean) of NPP WBt for the main tree species during the period 2008-2017 (B).

Forests 2018, 9 ,
x FOR PEER REVIEW 20 of 40 at, depending on the year, between DOY 295 and 311, On the other hand, NPPBM ends earlier in the autumn, achieving its 95% in average around the DOY 254.

Figure 11 .
Figure 11.Comparison of NPP estimates from eddy covariance (NPP EC ) and biometric measurements (NPP BM ) at Jastrebarsko forest.

Figure 12 .
Figure 12.Comparison of cumulative NPP determined with EC method (red line) and with biometric measurements (blue line) in Jastrebarsko pedunculate oak forest for years 2008 until 2017 (NPP EC -eddy covariance; NPP BM -biometric method; NPP WBt -total woody biomass including coarse roots; NPP T -twigs; NPP SB -stem and branches, points indicate Day-Of-Year (DOY) of dendrometer measurements).

Figure A2 .
Figure A2.Proposed theoretical framework for calculation of NPPBM.Logistic curve model of NPPWBt (Y) has an upper asymptote at NPPWBt annual and the inflexion point and the maximum of Y' at DOYmax_WBt.DOYmax_LF is assumed to occur at the second inflexion point of Y' i.e., the point where third derivative Y''' = 0 (panel A); annual dynamics of NPPWBt and NPPLF modelled using the same parameter k (see text), and their sum NPPBM (panel B).

Figure A3 .
Figure A3.3D scatter plot of NEE versus ΔSWCMay-Sep and mean annual solar radiation Rg.

Figure A2 .
Figure A2.Proposed theoretical framework for calculation of NPP BM .Logistic curve model of NPP WBt (Y) has an upper asymptote at NPP WBt annual and the inflexion point and the maximum of Y' at DOY max_WBt .DOY max_LF is assumed to occur at the second inflexion point of Y' i.e., the point where third derivative Y"' = 0 (panel A); annual dynamics of NPP WBt and NPP LF modelled using the same parameter k (see text), and their sum NPP BM (panel B).

Forests 2018, 9 , 40 Figure A2 .
Figure A2.Proposed theoretical framework for calculation of NPPBM.Logistic curve model of NPPWBt (Y) has an upper asymptote at NPPWBt annual and the inflexion point and the maximum of Y' at DOYmax_WBt.DOYmax_LF is assumed to occur at the second inflexion point of Y' i.e., the point where third derivative Y''' = 0 (panel A); annual dynamics of NPPWBt and NPPLF modelled using the same parameter k (see text), and their sum NPPBM (panel B).

Figure A3 .
Figure A3.3D scatter plot of NEE versus ΔSWCMay-Sep and mean annual solar radiation Rg.

Figure A3 .
Figure A3.3D scatter plot of NEE versus ∆SWC May-Sep and mean annual solar radiation R g .Forests 2018, 9, x FOR PEER REVIEW 28 of 40

Figure A4 .
Figure A4.3D scatter plot of NEE versus the principal component PC3 and PC4.

Figure A4 .
Figure A4.3D scatter plot of NEE versus the principal component PC3 and PC4.

Figure A5 .
Figure A5.Cumulative fluxes (NEE, GPP, RECO and NPP; panels A-D, respectively) estimated with EC during the period of the study.

Figure A5 .
Figure A5.Cumulative fluxes (NEE, GPP, R ECO and NPP; panels A-D, respectively) estimated with EC during the period of the study.
• 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).

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).

Table A1 .
Percentage of half-hourly NEE files before processing, after quality control filtering procedures and after u* filtering procedure.
Figure A6.Evolution, from the year 2008 (age 36) until 2017 (age 45), of tree's DBH distribution, by tree species, for stands around eddy covariance tower.

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

Table A3 .
Parameters of Schumacher and Hall function * (V = b 0 × DBH b 1 × h b 2 ) for the wood volume (including bark) of tree stem and branches **.

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

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.

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

Correlation between the NEE and the Principal Components (PCs)
# Statistically significant at p < 0.1.

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)).