Eddy Covariance vs. Biometric Based Estimates of

Full text

Turn on search term navigation

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). 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=^w′ ^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₂ 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=median(|_di−_Md|)

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

_di=(NE_Ei−NE_Ei−1)−(NE_Ei+1−NE_Ei)

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

_di<_Md−(z×MAD0.6745) OR _di>_Md+(z×MAD0.6745)

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:

_RECO(T)=_RECO _ref×^{e_E0(1_Tref−_T0−1T−_T0)}

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+_RECO

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_PEC=NEP+_Rh=−NEE+_Rh=−NEE+0.3919×_RECO

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). 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⁻³ 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:_{CS+B, i}(t) = _{V≥3 cm,i}(t) × BWD × CF}twigs:_{CT, i}(t) = 0.05 × _{CS+B, i}(t)roots:_{CR, i}(t) = RS ×(_{CS+B, i}(t)+_{CT, i}(t))total tree woody biomass:_{CWBt, i}(t) = (_{CS+B, i}(t)+_{CT, i}(t)+_{CR, i}(t)).}

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

NP_PWBt=NP_PWBt(plot)¯==1_Nplots∑plot1A×(_t2−_t1)×[∑i=1_Nalive(_{CWBt, i}(_t2)−_{CWBt, i}(_t1))+∑j=1_{Ndead or removed}(_{CWBt, j}(_tj)−_{CWBt, j}(_t1))]

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_{CWBt, j}(_tj)is estimated carbon stock of tree j at the time_tjof its death or removal (_t1≤_tj≤_t2).

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 (_NPPL), fruits (_NPPF), and their sum (_NPPLF).

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_Pi(DOY)=NP_{Pi_annual}1+^{e−k(DOY−DO_Ymax(i))}

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

_rmax=dNP_PmaxNP_Pannual .

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

k=4×atanh(_rmax)=2×ln(1+_rmax1−_rmax)

Inversely,

_rmax=tanh(k4)=^ek/4−^e−k/4^ek/4+^e−k/4

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:

DO_{Ymax_LF}=DO_{Ymax_WBt}−ln(2−3)k.

Finally,

NP_PBM(t) = NP_PWBt(t)+NP_PLF(t)

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

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). 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. 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 (_Rg¯), air temperature (_Ta¯), 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.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. 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.

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

Supplementary Materials

The following are available online at https://www.mdpi.com/1999-4907/9/12/764/s1, Supplementary Material S1: The EC and ancillary data (Jastrebarsko oak forest half-hourly eddy covariance and meteorological measurements, corrected daily precipitation, variable description, and processing codes).

Author Contributions

M.A, M.Z.O.S. and H.M. wrote the paper with equal contributions, prepared all data and made all tables and figures. H.M. conceptualized the experiment. A.P. and G.A. designed the EC system. A.P., G.A. and H.M. set up the EC system. H.M. maintained and is the PI of the EC site at Jastrebarsko forest. H.M., M.Z.O.S. and E.P. set up sampling plots, constructed and installed dendrometer bands. H.M., M.Z.O.S., E.P., I.B., M.A., G.T. with the help of technicians from Croatian Forest Research Institute performed field samplings and tree measurements. M.A. and H.M. processed the EC and meteorological data and sample plot measurements. D.V. and H.M. secured funding for the research. All Authors contributed to the interpretation of the results and to the discussion.

Funding

This research was funded by projects: Carbon balance drafting and new resources management tools according to Kyoto Protocol—CARBON-PRO, funded by the EU programme INTERREG-IIIB CADSES (ID 131669); Carbon balance and carbon cycling in young Pedunculate oak (Quercus robur L.) stands, funded by the Croatian Forest Ltd. (OKFŠ 2009–2010); Sustainability of carbon storage in managed forest of pedunculate oak, funded by the Croatian Forest Ltd. (OKFŠ 2011–2013); Estimating and forecasting forest ecosystems productivity by integrating field measurements, remote sensing and modelling—EFFEctivity, funded by the Croatian Science Foundation (HRZZ-UIP-11-2013-2492).

Acknowledgments

We would like to express our gratitude to the technicians of the Croatian Forest Research Institute who helped in setting-up, maintenance and measurements of sampling plots as well as to Diego Chiabà and Michel Zuliani for the help in setting the station. We are grateful to Croatian Forests Ltd. and Ivo Marjanović for their help during building of the tower infrastructure. We are grateful to Anikó Kern and Zoltán Barcza for their valuable suggestions and comments. This work M.A. was partly supported by a STSM Grant from COST Action FP1304. We are grateful to Academic Editor and four anonymous Reviewers for their valuable comments that helped improve the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Figure A1.Aerial image of study site within ‘Jastrebarsko forest’ where circular plots have been set up; Plots with dendrometer bands are marked with a red colour; flag shows the location of EC tower; forest compartment numbers are shown in white font. The management history: Since the establishment, stand tending were performed in the 1980s, 1990s with the last tending in 2002 in the south-western forest compartment (40A). The first pre-commercial thinning was made in 1999 in the north-western compartment (36A) with the thinning intensity of 11% by volume, followed by another one of only 1.5% in early 2008 (before the start of the vegetation season). The first thinning in the south-eastern forest compartment (41A) was made in 2006 with the thinning intensity of 7% by volume, followed by a thinning made in 2007 in the north-eastern compartment (37A) with the thinning intensity of 9% by volume. During the vegetation seasons 2008 until 2017 (i.e., the period of this study) there was no further thinning or harvesting activities in the above mentioned forest stands.

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.3D scatter plot of NEE versus ΔSWCMay–Sep and mean annual solar radiation Rg.

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

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

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

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

[ Table omitted. See PDF. ]

Table

Table A2.Parameters for species-specific, age-dependent tree height curves (h=(_α0+_α1×Age)×^{e−(_β0+_β1·Age)DBH}+1.30).

[ Table omitted. See PDF. ]

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

Table

Table A3.Parameters of Schumacher and Hall function * (V=_b0×^DBH_b1×^h_b2) for the wood volume (including bark) of tree stem and branches **.

[ Table omitted. See PDF. ]

* Unites: DBH (cm), h (m), V (m3). ** 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.

Table

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

[ Table omitted. See PDF. ]

Table

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

[ Table omitted. See PDF. ]

Table

Table A6.Correlation coefficients for the main annual fluxes (NEE, GPP, RECO) and the mean annual values of the main driver variables (Rg—global radiation, Ta—air temperature, Ts—soil temperature, P—precipitation, SWC—soil water content, ΔSWCMay-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 omitted. See PDF. ]

Table

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

[ Table omitted. See PDF. ]

# Statistically significant at p <0.1.

Table

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

[ Table omitted. See PDF. ]

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

References

1. Canadell, J.G.; Le Quere, C.; Raupach, M.R.; Field, C.B.; Buitenhuis, E.T.; Ciais, P.; Conway, T.J.; Gillett, N.P.; Houghton, R.A.; Marland, G. Contributions to accelerating atmospheric CO2 growth from economic activity, carbon intensity, and efficiency of natural sinks. Proc. Natl. Acad. Sci. USA 2007, 104, 18866–18870.

2. Nabuurs, G.J.; Lindner, M.; Verkerk, P.J.; Gunia, P.J.; Deda, P.; Michalak, R.; Grassi, G. First signs of carbon sink saturation in European forest biomass. Nat. Clim. Chang. 2013, 3, 792–796.

3. Hanewinkel, M.; Cullman, D.A.; Schelhaas, M.J.; Nabuurs, G.J.; Zimmermann, N.E. Climate change may cause severe loss in the economic value of European forest land. Nat. Clim. Chang. 2013, 3, 203–207.

4. Baldocchi, D. Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: Past, present and future. Glob. Chang. Biol. 2003, 9, 1–14.

5. Carrara, A.; Janssens, I.A.; Curiel Yuste, J.; Ceulemans, R. Seasonal changes in photosynthesis, respiration and nee of a mixed temperate forest. Agric. For. Meteorol. 2004, 126, 15–31.

6. Granier, A.; Breda, N.; Longdoz, B.; Gross, P.; Ngao, J. Ten years of fluxes and stand growth in a young beech forest at Hesse, North-Eastern France. Ann. For. Sci. 2008, 65, 704–717.

7. Baldocchi, D.; Chu, H.; Reichstein, M. Inter-annual variability of net and gross ecosystem carbon fluxes: A review. Agric. For. Meteorol. 2018, 249, 520–533.

8. Teets, A.; Fraver, S.; Hollinger, D.Y.; Weiskittel, A.R.; Seymour, R.S.; Richardson, A.D. Linking annual tree growth with eddy-flux measures of net ecosystem productivity across twenty years of observation in a mixed conifer forest. Agric. For. Meteorol. 2018, 249, 479–487.

9. Running, S.W.; Hunt, E.R.J. Generalization of a forest ecosystem process model for other biomes, BIOME-BGC, and an application for global-scale models. In Scaling Physiological Processes: Leaf to Globe; Ehleringer, J.R., Field, C., Eds.; Academic Press: San Diego, CA, USA, 1993; pp. 141–158. ISBN 012233440X.

10. Hidy, D.; Barcza, Z.; Marjanović, H.; Ostrogović Sever, M.Z.; Dobor, L.; Gelybo, G.; Fodor, N.; Pinter, K.; Churkina, G.; Running, S.; et al. Terrestrial Ecosystem Process Model Biome-BGCMuSo v4.0: Summary of improvements and new modelling possibilities. Geosci. Model. Dev. 2016, 9, 4405–4437.

11. Running, S.W.; Nemani, R.R.; Heinsch, F.A.; Zhao, M.; Reeves, M.; Hashimoto, H. A Continuous Satellite-Derived Measure of Global Terrestrial Primary Production. BioScience 2004, 54, 547–560.

12. FLUXNET2015 Dataset. Available online: http://fluxnet.fluxdata.org/data/ (accessed on 7 September 2018).

13. Burba, G. Eddy Covariance Method for Scientific, Industrial, Agricultural and Regulatory Applications. A Field Book on Measuring Ecosystem Gas. Exchange and Areal Emission Rates; LI-COR Biosciences: Lincoln, NE, USA, 2013; p. 332. ISBN 978-0-615-76827-4.

14. Schulze, E.-D. The Carbon and Nitrogen Cycle of Forest Ecosystems. In Carbon and Nitrogen Cycling in European Forest Ecosystems, Ecological Studies; Schulze, E.-D., Ed.; Springer: Berlin/Heidelberg, Germany, 2000; Volume 142, pp. 3–13. ISBN 978-3-642-57219-7.

15. Brunner, I.; Godbold, D.L. Tree roots in a changing world. J. For. Res. 2007, 12, 78–82.

16. Dietze, M.C.; Sala, A.; Carbone, M.S.; Czimczik, C.I.; Mantooth, J.A.; Richardson, A.D.; Vargas, R. Nonstructural carbon in woody plants. Annu. Rev. Plant Biol. 2014, 65, 667–687.

17. Reichstein, M.; Falge, E.; Baldocchi, D.; Papale, D.; Aubinet, M.; Berbigier, P.; Bernhofer, P.; Buchmann, N.; Gilmanov, T.; Granier, A.; et al. On the separation of net ecosystem exchange into assimilation and ecosystem respiration: Review and improved algorithm. Glob. Chang. Biol. 2005, 11, 1424–1439.

18. Ohtsuka, T.; Saigusa, N.; Koizumi, H. On linking multiyear biometric measurements of tree growth with eddy covariance-based net ecosystem production. Glob. Chang. Biol. 2009, 15, 1015–1024.

19. Curtis, P.S.; Hanson, P.J.; Bolstad, P.; Barford, C.; Randolph, J.C.; Schmid, H.P.; Wilson, K.B. Biometric and eddy-covariance based estimates of annual carbon storage in five Eastern North American deciduous forests. Agric. For. Meteorol. 2002, 113, 3–19.

20. Hanson, P.J.; Edwards, N.T.; Tschaplinski, T.J.; Wullschleger, S.D.; Joslin, J.D. Estimating the Net Primary and Net Ecosystem Production of a Southeastern Upland Quercus Forest from an 8-Year Biometric Record. In North American Temperate Deciduous Forest Responses to Changing Precipitation Regimes; Ecological Studies (Analysis and Synthesis), Volume 166; Hanson, P.J., Wullschleger, S.D., Eds.; Springer: New York, NY, USA, 2003.

21. Rocha, A.V.; Goulden, M.L.; Dunn, A.L.; Wofsy, S.C. On linking interannual tree ring variability with observations of whole-forest CO2 flux. Glob. Chang. Biol. 2006, 12, 1378–1389.

22. Gough, C.M.; Vogel, C.S.; Schmid, H.P.; Su, H.B.; Curtis, P.S. Multi-year convergence of biometric and meteorological estimates of forest carbon storage. Agric. For. Meteorol. 2008, 148, 158–170.

23. Gough, C.M.; Flower, C.E.; Vogel, C.S.; Dragoni, D.; Curtis, P.S. Whole-ecosystem labile carbon production in a north temperate deciduous forest. Agric. For. Meteorol. 2009, 149, 1531–1540.

24. Saleska, S.R.; Miller, S.D.; Matross, D.M.; Goulden, M.L.; Wofsy, S.C.; da Rocha, H.R.; de Camargo, P.B.; Crill, P.; Daube, B.C.; de Freitas, H.C.; et al. Carbon in amazon forests: Unexpected seasonal fluxes and disturbance-induced losses. Science 2003, 302, 1554–1557.

25. Gough, C.M.; Vogel, C.S.; Schmid, H.P.; Curtis, P.S. Controls on annual forest carbon storage: Lessons from the past and predictions for the future. Bioscience 2008, 58, 609–622.

26. Urbanski, S.; Barford, C.; Wofsy, S.; Kucharik, C.; Pyle, E.; Budney, J.; McKain, K.; Fitzjarrald, D.; Czikowsky, M.; Munger, W. Factors controlling CO2 exchange on timescales from hourly to decadal at Harvard Forest. J. Geophys. Res. 2007, 112, G02020.

27. Van der Molen, M.K.; Dolman, A.J.; Ciais, P.; Eglin, T.; Gobron, N.; Law, B.E.; Meir, P.; Peters, W.; Phillips, O.L.; Reichstein, M.; et al. Drought and ecosystem carbon cycling. Agric. For. Meteorol. 2011, 2011 151, 765–773.

28. Carbone, M.S.; Czimczik, C.I.; Keenan, T.F.; Murakami, P.F.; Pederson, N.; Schaberg, P.G.; Xu, X.; Richardson, A.D. Age, allocation and availability of nonstructural carbon in mature red maple trees. New Phytol. 2013, 200, 1145–1155.

29. Kern, A.; Marjanović, H.; Dobor, L.; Anić, M.; Hlásny, T.; Barcza, Z. Identification of Years with Extreme Vegetation State in Central Europe Based on Remote Sensing and Meteorological Data. South.-East Eur. For. 2017, 8, 1–20.

30. Ostrogović Sever, M.Z.; Paladinić, E.; Barcza, Z.; Hidy, D.; Kern, A.; Anić, M.; Marjanović, H. Biogeochemical Modelling vs. Tree-Ring Measurements—Comparison of Growth Dynamic Estimates at Two Distinct Oak Forests in Croatia. South.-East Eur. For. 2017, 8, 71–84.

31. Körner, C. Carbon limitation in trees. J. Ecol. 2003, 91, 4–17.

32. Babst, F.; Bouriaud, O.; Papale, D.; Gielen, B.; Janssens, I.A.; Nikinmaa, E.; Ibrom, A.; Wu, J.; Bernhofer, C.; Kostner, B.; et al. Above-ground woody carbon sequestration measured from tree rings is coherent with net ecosystem productivity at five eddy-covariance sites. New Phytol. 2014, 201, 1289–1303.

33. Barford, C.C.; Wofsy, S.C.; Goulden, M.L.; Munger, J.W.; Pyle, E.H.; Urbanski, S.P.; Hutyra, L.; Saleska, S.R.; Fitzjarrald, D.; Moore, K. Factors controlling long- and short-term sequestration of atmospheric CO2 in a mid-latitude forest. Science 2001, 294, 1688–1691.

34. Klepac, D.; Fabijanić, G. Management of Pedunculate Oak Forests. In Pedunculate Oak in Croatia, 1st ed.; Klepac, D., Dundović, J., Gračan, J., Eds.; Croatian Academy of Arts and Sciences, Centre for Scientific Work, Vinkovci and Croatian Forests Ltd.: Zagreb, Croatia, 1996; pp. 452–457. ISBN 953-154-079-9.

35. Ministry of Agriculture of the Republic of Croatia, Forest Management Area Plan of the Republic of Croatia; Croatian Forests Ltd.: Zagreb, Croatia, 2017; p. 859. Available online: http://www.mps.hr/hr/sume/sumarstvo/sumskogospodarska-osnova-2016-2025 (accessed on 12 September 2018).

36. IUSS Working Group WRB. World Reference Base for Soil Resources 2014 International Soil Classification System for Naming Soils and Creating Legends for Soil Maps, Update 2015; World Soil Resources Reports No. 106; FAO: Rome, Italy, 2015; p. 192. ISBN 978-92-5-108369-7.

37. Mayer, B. Hydropedological relations in the region of lowland forests of the Pokupsko basin. Rad. Šum. Inst. Jastrebarsko 1996, 31, 37–89. (In Croatian)

38. Marjanović, H.; Ostrogović, M.Z.; Alberti, G.; Balenović, I.; Paladinić, E.; Indir, K.; Peressotti, A.; Vuletić, D. Carbon dynamics in younger stands of Pedunculate oak during two vegetation periods. Šum. List (Spec. Issue) 2011, 135, 59–73, (In Croatian with English summary).

39. Ostrogović Sever, M.Z. Carbon Stocks and Carbon Balance of an Even-Aged Pedunculate Oak (Quercus robur L.) Forest in Kupa River Basin. Ph.D. Thesis, Faculty of Forestry, University of Zagreb, Zagreb, Croatia, 2013.

40. Wilkinson, M.; Eaton, E.L.; Broadmeadow, M.S.J.; Morison, J.I.L. Inter-annual variation of carbon uptake by a plantation oak woodland in south-eastern England. Biogeosciences 2012, 9, 5373–5389.

41. Vickers, D.; Mahrt, L. Quality Control and Flux Sampling Problems for Tower and Aircraft Data. J. Atmos. Ocean. Technol. 1997, 14, 512–526.

42. Aubinet, M.; Grelle, A.; Ibrom, A.; Rannik, U.; Moncrieff, J.; Foken, T.; Kowalski, A.S.; Martin, P.H.; Berbigier, P.; Bernhofer, C.; et al. Estimates of the annual net carbon and water exchange of forests: The euroflux methodology. Adv. Ecol. Res. 2000, 30, 113–175.

43. Wilczak, J.M.; Oncley, S.P. Sonic anemometer tilt correction algorithms. Bound.-Layer Meteorol. 2001, 99, 127–150.

44. Papale, D.; Reichstein, M.; Aubinet, M.; Canfora, E.; Bernhofer, C.; Kutsch, W.; Longdoz, B.; Rambal, S.; Valentini, R.; Vesala, T.; et al. Towards a standardized processing of Net Ecosystem Exchange measured with eddy covariance technique: Algorithms and uncertainty estimation. Biogeosciences 2006, 3, 571–583.

45. Baldocchi, D.; Valentini, R.; Running, S.; Oechel, W.; Dahlman, R. Strategies for measuring and modelling carbon dioxide and water vapour fluxes over terrestrial ecosystems. Glob. Chang. Biol. 1996, 2, 159–168.

46. Eddy Covariance Gap-Filling & Flux-Partitioning Tool. Available online: https://www.bgc-jena.mpg.de/bgi/index.php/Services/REddyProcWeb (accessed on 1 June 2018).

47. Falge, E.; Baldocchi, D.; Olson, R.; Anthoni, P.; Aubinet, M.; Bernhofer, C.; Burba, G.; Ceulemans, R.; Clement, R.; Dolman, H.; et al. Gap filling strategies for defensible annual sums of net ecosystem exchange. Agric. For. Meteorol. 2001, 107, 43–69.

48. Wutzler, T.; Lucas-Moffat, A.; Migliavacca, M.; Knauer, J.; Sickel, K.; Sigut, L.; Menzer, O.; Reichstein, M. Basic and extensible post-processing of eddy covariance flux data with REddyProc. Biogeosciences 2018, 15, 5015–5030.

49. Lloyd, J.; Taylor, J.A. On the temperature-dependence of soil respiration. Funct. Ecol. 1994, 8, 315–323.

50. Van Dijk, A.J.I.M.; Dolman, A.J. Estimates of CO2 uptake and release among European forests based on eddy covariance data. Glob. Chang. Biol. 2004, 10, 1445–1459.

51. Kljun, N.; Calanca, P.; Rotach, M.W.; Schmid, H.P. A simple two-dimensional parameterisation for Flux Footprint Prediction (FFP). Geosci. Model. Dev. 2015, 8, 3695–3713.

52. Flesch, T.K.; Wilson, J.D.; Yee, E. Backward-time lagrangian stochastic dispersion models and their application to estimate gaseous emissions. J. Appl. Meteorol. Climatol. 1995, 34, 1320–1332.

53. Kljun, N.; Rotach, M.W.; Schmid, H.P. A three-dimensional backward Lagrangian footprint model for a wide range of boundary-layer stratifications. Bound.-Lay Meteorol. 2002, 103, 205–226.

54. Flux Footprint Prediction (FFP) Online Data Processing. Available online: http://geography.swansea.ac.uk/nkljun/ffp/www/ (accessed on 13 April 2018).

55. Keeland, B.D.; Young, P.J. Installation of Traditional Dendrometer Bands; National Wetlands Research Center: Lafayette, LA, USA, 2015. Available online: https://www.nwrc.usgs.gov/topics/Dendrometer/ (accessed on 15 August 2018).

56. Marjanović, H. Modelling tree development and elements of stand structure in young stands of Pedunculate oak (Quercus Robur L.). Ph.D. Thesis, Faculty of Forestry, University of Zagreb, Zagreb, Croatia, 10 May 2009; p. 213, (In Croatian with English Summary).

57. Michailoff, I. Zahlenmäßiges Verfahren für die Ausführung der Bestandeshöhenkurven. [Numerical algorithm for the implementation of stand height curves]. Forstw. Cbl. u Thar. Jahrb. 1943, 6, 273–279. (In German)

58. Schumacher, F.X.; Hall, F.D.S. Logarithmic expression of timber-tree volume. J. Agric. Res. 1933, 47, 719–734.

59. Cestar, D.; Kovačić, Đ. Wood volume tables for Black Alder and Black Locust. Rad. Šum. Inst. Jastrebarsko 1982, 49, 1–149, (In Croatian with English Summary).

60. Špiranec, M. Wood volume tables. Rad. Šum. Inst. Jastrebarsko 1975, 22, 1–262, (In Croatian with German Summary).

61. Cestar, D.; Kovačić, Đ. Wood volume tables for Narrow-leaved Ash (Fraxinus parvifolia Auct.). Rad. Šum. Inst. Jastrebarsko 1984, 60, 1–178, (In Croatian with English Summary).

62. Balboa-Murias, M.A.; Rojo, A.; Álvarez, J.G.; Merino, A. Carbon and nutrient stocks in mature Quercus robur L. stands in NW Spain. Ann. For. Sci. 2006, 63, 557–565.

63. Cairns, M.A.; Brown, S.; Helmer, E.H.; Baumgardner, G.A. Root biomass allocation in the world’s upland forests. Oecologia 1997, 111, 1–11.

64. Šumarska Enciklopedija, 1st ed.; Leksikografski zavod FNRJ: Zagreb, Croatia, 1959.

65. IPCC. Good Practice Guidance for Land Use, Land-Use Change and Forestry; Penman, J., Gytarsky, M., Hiraishi, T., Kruger, D., Pipatti, R., Buendia, L., Miwa, K., Ngara, T., Tanabe, K., Wagner, F., et al., Eds.; IPCC/IGES: Hayama, Japan, 2003.

66. IPCC. 2006 IPCC Guidelines for National Greenhouse Inventories; Eggleston, H.S., Buendia, L., Miwa, K., Ngara, T., Tanabe, K., Eds.; IGES: Kobe, Japan, 2006.

67. Martínez-Vilalta, J. Carbon storage in trees: Pathogens have their say. Tree Physiol. 2014, 34, 215–217.

68. Richardson, A.D.; Keenan, T.F.; Migliavacca, M.; Ryu, Y.; Sonnentag, O.; Toomey, M. Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. Agric. For. Meteorol. 2013, 169, 156–173.

69. Locosselli, G.M.; Buckeridge, M.S. Dendrobiochemistry, a missing link to further understand carbon allocation during growth and decline of trees. Trees-Struct. Funct. 2017, 31, 1745–1758.

70. Jolliffe, I.T. A note on the use of principal components in regression. J. R. Stat. Soc. Ser. C Appl. Stat. 1982, 31, 300–303.

71. Pretzsch, H. Forest Dynamics, Growth and Yield: From Measurement to Model; Springer: Berlin, Germany, 2009.

72. Pilegaard, K.; Ibrom, A.; Courtney, M.S.; Hummelshøj, P.; Jensen, N.O. Increasing net CO2 uptake by a Danish beech forest during the period from 1996 to 2009. Agric. For. Meteorol. 2011, 151, 934–946.

73. Herbst, M.; Mund, M.; Tamrakar, R.; Knohl, A. Differences in carbon uptake and water use between a managed and an unmanaged beech forest in central Germany. For. Ecol. Manag. 2015, 355, 101–108.

74. Novick, K.A.; Oishi, A.C.; Ward, E.J.; Siqueira, M.B.S.; Juang, J.-Y.; Stoy, P.C. On the difference in the net ecosystem exchange of CO2 between deciduous and evergreen forests in the southeastern United States. Glob. Chang. Biol. 2015, 21, 827–842.

75. Sulman, B.N.; Roman, D.T.; Scanlon, T.M.; Wang, L.; Novick, K.A. Comparing methods for partitioning a decade of carbon dioxide and water vapor fluxes in a temperate forest. Agirc. For. Meteorol. 2016, 226–227, 229–245.

76. Malhi, Y.; Baldocchi, D.; Jarvis, P.G. The carbon balance of tropical, temperate and boreal forests. Plant Cell Environ. 1999, 22, 715–740.

77. Law, B.E.; Falge, E.; Gu, L.; Baldocchi, D.D.; Bakwin, P.; Berbigier, P.; Davis, K.; Dolman, A.J.; Falk, M.; Fuentes, J.D.; et al. Environmental controls over carbon dioxide and water vapor exchange of terrestrial vegetation. Agric. For. Meteorol. 2002, 113, 97–120.

78. Delpierre, N.; Berveiller, D.; Granda, E.; Dufrene, E. Wood phenology, not carbon input, controls the interannual variability of wood growth in a temperate oak forest. New Phytol. 2016, 210, 459–470.

79. Froelich, N.; Croft, H.; Chen, J.M.; Gonsamo, A.; Staebler, R.M. Trends of carbon fluxes and climate over a mixed temperate-boreal transition forest in southern Ontario, Canada. Agric. For. Meteorol. 2015, 211, 72–84.

80. Saigusa, N.; Yamamoto, S.; Murayama, S.; Kondo, H. Inter-annual variability of carbon budget components in an AsiaFlux forest site estimated by long-term flux measurements. Agric. For. Meteorol. 2005, 134, 4–16.

81. Zscheischler, J.; Fatichi, S.; Wolf, S.; Blanken, P.D.; Bohrer, G.; Clark, K.; Desai, A.R.; Hollinger, D.; Keenan, T.; Novick, K.A.; et al. Short-term favorable weather conditions are an important control of interannual variability in carbon and water fluxes. J. Geophys. Res.-Biogeosci. 2016, 121, 2186–2198.

82. Odum, E.P. Strategy of ecosystem development. Science 1969, 164, 262–270.

83. Curtis, P.S.; Gough, C.M. Forest aging, disturbance and the carbon cycle. New Phytol. 2018, 219, 1188–1193.

84. Gough, C.M.; Vogel, C.S.; Hardiman, B.; Curtis, P.S. Wood net primary production resilience in an unmanaged forest transitioning from early to middle succession. For. Ecol. Manag. 2010, 260, 36–41.

85. Peterken, G.F. Woodland Conservation and Management, 2nd ed.; Chapman & Hall: London, UK, 1993.

86. Ortmann-Ajkai, A.; Csicsek, G.; Lukacs, M.; Horvath, F. Regeneration patterns in a pedunculate oak (Quercus robur L.) strict forest reserve in southern Hungary. Šum. List 2017, 141, 39–46.

87. Gower, S.T.; McMurtrie, R.E.; Murty, D. Aboveground net primary production decline with stand age: Potential causes. Trends Ecol. Evol. 1996, 11, 378–382.

88. Xu, C.Y.; Turnbull, M.H.; Tissue, D.T.; Lewis, J.D.; Carson, R.; Schuster, W.S.F.; Whitehead, D.; Walcroft, A.S.; Li, J.B.; Griffin, K.L. Age-related decline of stand biomass accumulation is primarily due to mortality and not to reduction in NPP associated with individual tree physiology, tree growth or stand structure in a quercus-dominated forest. J. Ecol. 2012, 100, 428–440.

89. Barić, L.; Županić, M.; Pernek, M.; Diminić, D. First records of Chalara fraxinea in Croatia—A new agent of ash dieback (Fraxinus spp.). Šum. list 2012, 136, 461–469.

90. Claessens, H.; Oosterbaan, A.; Savill, P.; Rondeux, J. A review of the characteristics of black alder (Alnus glutinosa (L.) Gaertn.) and their implications for silvicultural practices. Forestry 2010, 83, 163–175.

91. Black, K.; Bolger, T.; Davis, P.; Nieuwenhuis, M.; Reidy, B.; Saiz, G.; Tobin, B.; Osborne, B. Inventory and eddy covariance-based estimates of annual carbon sequestration in a sitka spruce (Picea sitchensis (Bong.) Carr.) forest ecosystem. Eur. J. For. Res. 2007, 126, 167–178.

92. Peichl, M.; Brodeur, J.J.; Khomik, M.; Arain, M.A. Biometric and eddy-covariance based estimates of carbon fluxes in an age-sequence of temperate pine forests. Agric. For. Meteorol. 2010, 150, 952–965.

93. Campioli, M.; Malhi, Y.; Vicca, S.; Luyssaert, S.; Papale, D.; Penuelas, J.; Reichstein, M.; Migliavacca, M.; Arain, M.A.; Janssens, I.A. Evaluating the convergence between eddy-covariance and biometric methods for assessing carbon budgets of forests. Nat. Commun. 2016, 7.

94. Hanson, P.J.; Edwards, N.T.; Garten, C.T.; Andrews, J.A. Separating root and soil microbial contributions to soil respiration: A review of methods and observations. Biogeochemistry 2000, 48, 115–146.

95. Kutsch, W.; Bahn, M.; Heinemeyer, A. Soil Carbon Dynamics: An Integrated Methodology; Cambridge University Press: Cambridge, UK, 2009; p. 286. ISBN 978-0-521-86561-6.

96. Ohtsuka, T.; Mo, W.H.; Satomura, T.; Inatomi, M.; Koizumi, H. Biometric based carbon flux measurements and net ecosystem production (NEP) in a temperate deciduous broad-leaved forest beneath a flux tower. Ecosystems 2007, 10, 324–334.

97. Subke, J.A.; Inglima, I.; Cotrufo, M.F. Trends and methodological impacts in soil CO2 efflux partitioning: A meta-analytical review. Glob. Chang. Biol. 2006, 12, 921–943.

98. Bond-Lamberty, B.; Bailey, V.L.; Chen, M.; Gough, C.M.; Vargas, R. Globally rising soil heterotrophic respiration over recent decades. Nature 2018, 560, 80–83.

99. Mori, S.; Yamaji, K.; Ishida, A.; Prokushkin, S.G.; Masyagina, O.V.; Hagihara, A.; Hoque, A.; Suwa, R.; Osawa, A.; Nishizono, T. Mixed-power scaling of whole-plant respiration from seedlings to giant trees. Proc. Natl. Acad. Sci. USA 2010, 107, 1447–1451.

100. Turner, D.P.; Ritts, W.D.; Cohen, W.B.; Maeirsperger, T.K.; Gower, S.T.; Kirschbaum, A.A.; Running, S.W.; Zhao, M.; Wofsy, S.C.; Dunn, A.L.; et al. Site-level evaluation of satellite-based global terrestrial gross primary production and net primary production monitoring. Glob. Chang. Biol. 2005, 11, 666–684.

101. Neumann, M.; Zhao, M.; Kindermann, G.; Hasenauer, H. Comparing MODIS net primary production estimates with terrestrial national forest inventory data in Austria. Remote Sens. 2015, 7, 3878.

102. Kern, A.; Marjanović, H.; Barcza, Z. Evaluation of the Quality of NDVI3g Dataset against Collection 6 MODIS NDVI in Central Europe between 2000 and 2013. Remote Sens. 2016, 8, 955.

103. Gaudinski, J.B.; Torn, M.S.; Riley, W.J.; Swanston, C.; Trumbore, S.E.; Joslin, J.D.; Majdi, H.; Dawson, T.E.; Hanson, P.J. Use of stored carbon reserves in growth of temperate tree roots and leaf buds: Analyses using radiocarbon measurements and modeling. Glob. Chang. Biol. 2009, 15, 992–1014.

104. Pallardy, S.G. Physiology of Woody Plants, 3rd ed.; Elsevier: Amsterdam, The Netherlands; London, UK; Boston, MA, USA, 2008.

105. Barbaroux, C.; Breda, N.; Dufrene, E. Distribution of above-ground and below-ground carbohydrate reserves in adult trees of two contrasting broad-leaved species (Quercus petraea and Fagus sylvatica). New Phytol. 2003, 157, 605–615.

106. Palacio, S.; Camarero, J.J.; Maestro, M.; Alla, A.Q.; Lahoz, E.; Montserrat-Martí, G. Are storage and tree growth related? Seasonal nutrient and carbohydrate dynamics in evergreen and deciduous Mediterranean oaks. Trees 2018, 32, 777–790.

107. Begum, S.; Kudo, K.; Rahman, M.H.; Nakaba, S.; Yamagishi, Y.; Nabeshima, E.; Nugroho, W.D.; Oribe, Y.; Kitin, P.; Jin, H.O.; et al. Climate change and the regulation of wood formation in trees by temperature. Trees-Struct. Funct. 2018, 32, 3–15.

AuthorAffiliation

¹Department of Forest Management and Forestry Economics, Croatian Forest Research Institute, Cvjetno Naselje 41, 10450 Jastrebarsko, Croatia

²Department of Agricultural, Food, Environmental and Animal Sciences, University of Udine, Via Delle Scienze 206, 33100 Udine, Italy

³Institute of Biometeorology, National Research Council of Italy, Via Caproni 8, 50145 Firenze, Italy

⁴Department of Geophysics, Faculty of Science, University of Zagreb, Horvatovac 95, 10000 Zagreb, Croatia

^*Author to whom correspondence should be addressed.

Word count: 14369

Show less

© 2018. This work is licensed under https://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Abstract

Translate

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.

Details

Title

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

Author

Anić, Mislav; Maša Zorana Ostrogović Sever; Alberti, Giorgio; Balenović, Ivan; Paladinić, Elvis; Peressotti, Alessandro; Tijan, Goran; Večenaj, Željko; Vuletić, Dijana; Marjanović, Hrvoje

Publication year

2018

Publication date

Dec 2018

Publisher

MDPI AG

e-ISSN

19994907

Source type

Scholarly Journal

Language of publication

English

DOI

https://doi.org/10.3390/f9120764

ProQuest document ID

2174412781

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

Jump to:

Full text

Abstract

Details

Suggested sources