SEEFOR 8 (2): 71-84
Article ID: 158
DOI: https://doi.org/10.15177/seefor.17-17

Original scientific paper

Biogeochemical Modelling vs. Tree-Ring Measurements - Comparison of Growth Dynamic Estimates at Two Distinct Oak Forests in Croatia

Maša Zorana Ostrogović Sever1, Elvis Paladinić1, Zoltán Barcza2,3,4, Dóra Hidy5, Anikó Kern6, Mislav Anić1, Hrvoje Marjanović7*

(1) Croatian Forest Research Institute, Division for Forest Management and Forestry Economics, Trnjanska cesta 35, HR-10000 Zagreb, Croatia;
(2) Eötvös Loránd University, Department of Meteorology, Pázmány P. st. 1/A, H-1117 Budapest, Hungary;
(3) Eötvös Loránd University, Faculty of Sciences, Excellence Center, H-2462 Martonvásár, Brunszvik u. 2., Hungary;
(4) Czech University of Life Sciences Prague, Faculty of Forestry and Wood Sciences, Kamýcká 129, 165 21 Prague 6, Czech Republic;
(5) Szent István University, MTA-SZIE Plant Ecology Research Group, Páter K. u.1., H-2103 Gödöllő, Hungary;
(6) Eötvös Loránd University, Department of Geophysics and Space Science, Pázmány P. st. 1/A, H-1117 Budapest, Hungary;
(7) Croatian Forest Research Institute, Division for Forest Management and Forestry Economics, Cvjetno naselje 41, HR-10450 Jastrebarsko, Croatia

* Correspondence: e-mail: hrvojem@sumins.hr

Citation: OSTROGOVIĆ SEVER MZ, PALADINIĆ E, BARCZA Z, HIDY D, KERN A, ANIĆ M, MARJANOVIĆ H 2017 Biogeochemical Modelling vs. Tree-Ring Measurements - Comparison of Growth Dynamic Estimates at Two Distinct Oak Forests in Croatia. South-east Eur for 8 (2): 71-84. DOI: https://doi.org/10.15177/seefor.17-17

Received: 28 Jun 2017; Revised: 29 Oct 2017; 29 Nov 2017; 12 Dec 2017; Accepted: 14 Dec 2017; Published online: 21 Dec 2017

Cited by:     Crossref  (1)     Google Scholar

Abstract

Background and Purpose: Biogeochemical process‑based models use a mathematical representation of physical processes with the aim of simulating and predicting past or future state of ecosystems (e.g. forests). Such models, usually executed as computer programs, rely on environmental variables as drivers, hence they can be used in studies of expected changes in environmental conditions. Process‑based models are continuously developed and improved with new scientific findings and newly available datasets. In the case of forests, long-term tree chronologies, either from monitoring or from tree-ring data, offer valuable means for testing modelling results. Information from different tree cores can cover a wide range of ecological and meteorological conditions and as such provide satisfactory temporal and spatial resolution to be used for model testing and improvement.
Materials and Methods:
In our research, we used tree-ring data as a ground truth to test the performance of Biome-BGCMuSo (BBGCMuSo) model in two distinct pedunculate oak forest areas, Kupa River Basin (called Pokupsko Basin) and Spačva River Basin, corresponding to a wetter and a drier site, respectively. Comparison of growth estimates from two different data sources was performed by estimating the dynamics of standardized basal area increment (BAI) from tree-ring data and standardized net primary productivity of stem wood (NPPw) from BBGCMuSo model. The estimated growth dynamics during 2000-2014 were discussed regarding the site-specific conditions and the observed meteorology.
Results: The results showed similar growth dynamic obtained from the model at both investigated locations, although growth estimates from tree-ring data revealed differences between wetter and drier environment. This indicates higher model sensitivity to meteorology (positive temperature anomalies and negative precipitation anomalies during vegetation period) than to site-specific conditions (groundwater, soil type). At both locations, Pokupsko and Spačva, BBGCMuSo showed poor predictive power in capturing the dynamics obtained from tree‑ring data.
Conclusions: BBGCMuSo model, similar to other process-based models, is primarily driven by meteorology, although site-specific conditions are an important factor affecting lowland oak forests’ growth dynamics. When possible, groundwater information should be included in the modelling of lowland oak forests in order to obtain better predictions. The observed discrepancies between measured and modelled data indicate that fixed carbon allocation, currently implemented in the model, fails in predicting growth dynamics of NPP. Dynamic carbon allocation routine should be implemented in the model to better capture tree stress response and growth dynamics.

Keywords: Pedunculate oak forest, basal area increment, net primary productivity, model testing

INTRODUCTION

In a changing environment, there is a growing need for estimating future forest productivity in order to forecast the impact of climate change on sustainability and adaptability of forests. Process-based modelling is a state-of-the-art technique used in predicting behaviour and future state of ecosystems with respect to environmental conditions [1, 2]. A variety of known ecophysiological and geochemical processes are implemented in these models, but continuous model development based on new knowledge is still needed [3, 4].

Process models, used for vegetation modelling, are complex and often have a high number of driving variables and parameters. This makes calibration (i.e. parameterisation) and validation of such complex vegetation models a challenging task [5]. Model calibration requires an extensive dataset of field measurements, as well as high computational skills and computing power. Most valuable source of field data, used in calibration and validation of vegetation process models, is high‑frequency (i.e. half-hour) eddy-covariance (EC) data [6]. A global network of EC flux measurement sites, such as FLUXNET [7], has great potential in facilitating means for better understanding of carbon dynamics in various biomes across regional and global scales [8]. However, for a particular or specific ecosystem, such as lowland forests of pedunculate oak, this dataset has limited use due to a relatively scarce spatial distribution of flux towers. In addition, even if flux measurements for the selected forest type do exist, single site measurements cover a relatively small area (few hundred meters to few kilometres). Taller towers are capable of covering even larger areas, such as the 82 m high tower in Hegyhátsal, Hungary, but in that case, fluxes reflect a multitude of different land covers [9].

Other sources of data that might be useful in assessing the results of modelling with process‑based models are long-term chronologies from monitoring or from tree-ring data. Databases of tree ring measurements usually cover a wide range of ecological and meteorological conditions. As such, these data contain information of satisfactory temporal and spatial frequency to be used for testing performance of complex process-based models. New knowledge gained through model comparison with various measurement datasets is a valuable source of information to be used for model improvements and further model developments [7, 10].

The dendrochronological approach provides a unique long-term understanding of the interplay between terrestrial ecosystems and external forcing agents [11, 12]. It is most suitable for trees of the temperate climate zone. Tree-ring width (TRW) data and its derived variables (e.g. tree basal area increment, BAI), reflect tree’s radial growth due to cambial activity. Tree’s stem growth at a given year often integrates the meteorological effects of the current year and several previous years, and it is further modified by site-specific conditions and management [13]. In this way, in their annual rings, trees preserve an archive of past growing conditions reflecting climate anomalies, competition, disturbance, soil characteristics or species-specific growth patterns [14, 15], as well as human-induced disturbances. Therefore, when using or interpreting tree-ring width data all these influencing factors should be kept in mind.

According to Hafner et al. [16], when analysing the response of lowland oak forests to climate conditions, it is important to consider the micro-environment (e.g. drier vs. wetter sites), but also to distinguish tree-ring data into early- and latewood formations. The underlying idea behind this research was to test modelling performance by using tree-ring data as ground truth, rather than to analyse the climatic influence on tree-ring formations. Therefore, we used a whole tree‑ring width as a proxy for realised annual growth to test the modelled growth dynamics at different locations. Simple visual interpretation of tree growth response to the observed meteorology was performed only with the purpose of providing additional insight into differences between investigated locations which could further be used in defining potential issues in the model logic.

In our research, TRW data, combined with dendrometric data (i.e. diameter at breast height), were used to assess the inter-annual variability of productivity in lowland oak forests. The aim of this study is to test modelling performance by comparing forest growth dynamics estimates from Biome-BGCMuSo model (BBGCMuSo) against the observed growth estimated from an extensive dataset of tree-rings. The observed differences will serve for defining modelling issues and indicating potential directions for further model improvements. There is evidence of growth decrease of pedunculate oak forests in Southeast Europe as a response to a change of water regime and climate [17]. Reliable model predictions are needed for the selection of appropriate adaptation measures for the preservation of those forests.

MATERIALS AND METHODS

Study Areas

The research was conducted in two distinct areas of managed pedunculate oak forest in Croatia, Kupa River Basin (also called Pokupsko Basin) located in western part of Croatia, approximately 35 km SW of Zagreb, and Spačva Basin located in eastern part of Croatia, approximately 35 km SE of Vinkovci (Figure 1).

The dominant tree species in both forest complexes is pedunculate oak (Quercus robur L.) with a significant share of common hornbeam (Carpinus betulus L.), and narrow-leaved ash (Fraxinus angustifolia Vahl.). Black alder (Alnus glutinosa (L.) Geartn.) is also present, but more abundantly in Pokupsko Basin where it holds a significant share in stock (Table 1). Oak forests in Croatia are managed as even-aged, with 140 years long rotations that end with two or three regeneration cuts during last 10 years of the rotation.

TABLE 1. Site description

Floodplain forests of Pokupsko Basin grow in the tectonic basin “Crna Mlaka” surrounded by hilly slopes of Samobor, Žumberak and Vukomerec hills on east, west and north side, and Kupa River on the south. The Basin lies between 15°32’ and 15°50’ longitude east, and 45°30’ to 45°42’ latitude north occupying mostly flat area, with altitude ranging from 107 to 115 m a.s.l. The climate in Pokupsko Basin is warm temperate with a mean annual air temperature of 10.6°C and precipitation of 962 mm·y-1 for the period 1981-2010 (data obtained from national Meteorological and Hydrological Service for nearest meteorological station in Jastrebarsko, 45°40’N, 15°39’E, 140 m a.s.l., approx. 2 km NW of the Pokupsko Basin forest). Soils are hydromorphic on clay parent material and according to the World Reference Base for Soil Resources [18], they are classified as luvic stagnosol (Table 1). Average groundwater table depth (based on the data from previous measurements until 1997 and those from 2008 onwards, made by the researchers of Croatian Forest Research Institute), is from 60 to 200 cm [19].

The forest complex of Spačva Basin lays at the most eastern part of Croatia, between Sava and Drava rivers, on the catchment area of Bosut River and its tributaries. Located between 18°45’ and 19°10’ longitude east and 44°51’ to 45°09’ latitude north, it occupies flat-curly basin of altitude ranging from 77 to 90 m a.s.l., which is intersected by numerous small rivers. According to Seletković [20], the climate in the eastern part of pedunculate oak distribution area in Croatia is warm temperate with maximum rainfall in June, without exceptionally dry months in summer, and with driest months occurring during cold period of the year. The mean annual air temperature is 11.5°C and precipitation is 686 mm·y-1 for the period 1981-2010 (data obtained from National Meteorological and Hydrological Service for nearest meteorological station Gradište, 45°10’N, 18°42’E, 89 m a.s.l., approx. located 4 km W of the Spačva forest). Majority of Spačva Basin forest soils are semi-terrestrial or hydromorphic soils on loamy-clay river sediments [21], and according to World Reference Base for Soil Resources [18], they are classified as chernozem. In the period from 1996 to 2012, an average observed groundwater table depth was ranging from 139 to 617 cm [22]. Differences between two forest complexes are summarized in Table 1.

Biome-BGCMuSo Model

Biome-BGCMuSo (BBGCMuSo) [23, 24] is a newer version of the original biogeochemical model Biome-BGC that simulates carbon, nitrogen, and water cycling in different terrestrial ecosystems [25]. In general, Biome-BGC is a process-based model widely used for estimating ecosystem productivity under current and changed environmental conditions [23, 26-30]. Major improvements in BBGCMuSo include introducing a multilayer soil module with the possibility of using groundwater table information, management module, new plant pools, respiration acclimation, CO2 regulation of stomatal conductance and transient run [23, 24].

There are two obligatory input datasets for running the model, namely meteorology and ecophysiological traits of the specific ecosystem, and several optional datasets, e.g. atmospheric CO2 concentration, nitrogen deposition, management data, groundwater table etc. Model simulation has three steps: spin-up, transient run and normal run. The purpose of spin-up is to bring the ecosystem to the steady state regarding soil carbon stocks using long-term local meteorological data. Transient run enables a smooth transition from spin-up phase to the normal phase as it slowly brings the ecosystem to steady state under current (changed) environmental conditions using varying data on CO2 concentration, nitrogen deposition and management. Finally, the normal run is done using current meteorology, CO2 concentration, nitrogen deposition, and management for the period of interest.

In this research, we simulated the productivity of selected stands in managed oak forests in two selected areas, Pokupsko Basin (4 forest management units with 947 forest compartments covering 11.1 kha in total) and Spačva Basin (13 forest management units with a total of 2918 forest compartments covering 47.8 kha in total). Only forest compartments categorised, according to the dominant tree species, in forest management plans as management class of pedunculate oak, older than 15 years in the year 2000, and for which regeneration harvests have not occurred between 2000 and 2014 were considered for simulation and tree coring (potentially 6.4 kha in 524 compartments in Pokupsko and 33.2 kha in 2083 compartments in Spačva Basin). For the selected 2607 forest compartments in both areas, we run the simulations. However, in the comparison of modelled and measured results only those forest compartments where we actually performed measurements were considered (36 in Pokupsko and 55 in Spačva Basin). More details on the selection of forest compartments are given in the next chapter.

Model simulation was performed on a forest stand level since each forest compartment corresponds to a single forest stand. To account for different management history of stands of different age (i.e. management compartments; for age class distribution see Figure 2), we set the spin-up and transient simulations to the period 1850-1999, while the normal run was done for the period of interest, i.e. 2000-2014. From the records in forest management plans for the year and volume of thinning / regeneration harvests, and the growing stocks at the forest compartment level we calculated the intensities. Specific intensity values were used for the simulation of each forest compartment in the normal run. For the spin-up and transient run we calculated the average intensity of thinning and regeneration harvests and applied those values to all forest compartments. However, the timing (i.e. the year) of the thinning was estimated from the existing records of stand age and year of thinning in the 10-year steps (e.g. thinning in 2009 in the records implies thinning in 1999, 1989, etc., until the final harvest of the previous stand).

Meteorological data used in the simulation was obtained from FORESEE database [31]. FORESEE is a gridded database with a spatial resolution of 1/6o x 1/6o  containing daily maximum/minimum temperature and precipitation fields for Central Europe. In addition, the FORESEE retrieval site (http://nimbus.elte.hu/FORESEE/map_query/index.html) offers the possibility for retrieval of core meteorological variables needed for running Biome-BGCMuSo, namely: the daily minimum, maximum and average daytime (from sunrise to sunset) temperature (°C), daily total precipitation (cm), daylight average vapour pressure deficit (Pa), shortwave radiant flux density (W·m-2) and day length from sunrise to sunset (seconds). By overlapping FORESEE database over a spatial distribution of the selected management compartments, a specific meteorological dataset was assigned to each forest compartment.Considering that this dataset covers the period from 1951, also taking into account that the current minimum prescribed rotation length for pedunculate oak is 140 years, we needed to approximate to the meteorology from 1851. For the purpose of our simulations, we assumed that meteorology for 1851-1950 was the same as it was during 1951-1970. Therefore, we used multiple times data from the period 1951-1970 without randomization. For ecophysiological traits, we used a parameter list for oak forests published in Hidy et al. [24], slightly adjusted to site-specific conditions (Table 2). The main difference between two investigated sites is a share of black alder. Black alder is a nitrogen-fixing species, and therefore a higher nitrogen fixation rate [32], relative to the share of black alder [33]), is used at Pokupsko site. Values of other adjusted parameters are set to default [24] (Table 2). Spatially explicit data on stand elevation was obtained from Croatian Forests Ltd. database containing all information on forest stands, as prescribed by the national regulation [34], while for stand latitude, we used the latitude of the corresponding FORESEE pixel [31, 35]. Site‑specific soil texture was calculated from previously collected soil data, resulting with one texture that was used in simulations at Pokupsko and the other at Spačva Basin [19, 36].

TABLE 2. Parameter values adjusted to site-specific conditions.

Furthermore, we used three optional input files: atmospheric CO2, nitrogen deposition, and management file. Atmospheric CO2 concentration data were obtained from Mauna Loa Observatory (available online at http://www.esrl.noaa.gov/gmd/obop/mlo/) and using relevant publications [37]. Nitrogen deposition data was based on Churkina et al. [38]. Management data (stand age, wood volume stocks by tree species, the volume of wood extracted with thinning or stand regeneration and year of the activities, soil type, etc.) was obtained from Croatian Forests Ltd. database. For transient run, management was reconstructed using available information on the stand age in each management compartment. The final cut year was set as the first year of stand development, and from that year onward thinning events were set at every 10 years using average thinning rate of 15%. For the normal run, the actual thinning rates were used for each forest compartment. Thinning rates were estimated from records of standing wood stock, estimated annual wood increment, year of thinning and volume of extracted wood available from the Croatian Forests Ltd. database. In order to use groundwater table information, the user should provide daily data. Unfortunately, daily data on groundwater table was not available for both investigated sites; therefore, this model feature was not used.

Tree-Ring Data

A field survey was conducted from spring 2015 to spring 2016. This research was part of the project EFFEctivity (http://www.sumins.hr/en/projekti/effectivity/) which had, in addition to the work presented here, also the goal of testing MODIS MOD17 annual Net Primary Productivity product [39]. This determined the design for the selection of the plot locations. In short, both forest areas, Pokupsko and Spačva, were overlaid with grid corresponding to MODIS 1 km resolution pixels. Only pixels with more than 90% forest cover and with homogenous age structure (forest compartments consisting of >70% of the pixel area had to have the stand age difference of less than 40 years) were selected and in each pixel four plots were installed. The location of the plot was at the centre of the MODIS pixels with 500 m resolution (each 1 km MODIS pixel can be subdivided into four 500 m pixels). The example of the plot layout is presented in Figure 3. In total, 109 temporary circular plots were placed within two investigated forest areas (41 plots in Pokupsko and 68 in Spačva). Sampling radius varied depending on the stand age and tree size of the sampled tree with larger trees being sampled using larger radius [40] (Table 3). On each plot diameter at breast height (dbh, 1.30 m above the ground) of all sampled trees, as well as tree location on the plot (i.e. distance from the plot centre and azimuth) were recorded.

TABLE 3. The radius of sampling with respect to stand age and tree size.

Tree cores were taken, on average, from 9.6 dominant and co-dominant trees per plot (min. 5, max. 13). In total, 1051 cores were collected, out of which 383 in Pokupsko Basin (247 Q. robur, 21 C. betulus, 34 F. angustifolia, 75 A. glutinosa, 6 other) and 668 in Spačva Basin (512 Q. robur, 44 C. betulus, 112 F. angustifolia).

One core per tree was taken at 1.30 m from the ground from the stem side facing the plot centre using increment borer (Haglof, Sweden) of 5.15 mm inner diameter. The collected cores were air-dried in the laboratory for several days and stored in the refrigerator at 4°C until further analysis. Following standard dendrochronological preparation methods, outlined in Stokes and Smiley [41]), the cores were glued to wooden holders and placed into a press for a day. Afterward, they were sanded with progressively finer grades of sandpaper (i.e. 120, 180, 240 and 320 grit). Finally, cores were scanned at high resolution (2400 DPI) and the scanned images were saved into the tree-core database on a local network drive for later TRW measurements.

Tree-ring widths were measured from scanned images to the nearest 0.001 mm using PC and specialized CooRecorder software v.7.8.1 (Cybis Elektronik & Data AB, Sweden). TRW measurements were corrected and underwent quality control through repeated cross-check routines for cross-dating and identification of measurement errors with COFECHA computer program [42, 43].

Data Analysis

The comparison of growth dynamics from 2000 to 2014 in two distinct locations, from two different data sources, was performed using basal area growth estimated from tree-rings and net primary productivity obtained from the BBGCMuSo model. To exclude age-related trend associated with TRW data we used basal area increment (BAI) as proposed in Biondi and Qeadan [44]. BAI was calculated using tree diameter at breast height, measured at the time of tree coring, and TRW data.

Net primary productivity (NPP) obtained from the BBGCMuSo model comprises net productivities of different tree parts. To be able to make a comparison of BAI with the model NPP data, we estimated the net primary productivity of stem wood (NPPw) using carbon allocation ratios (Table 4) in the following way:

NPP = NPPw + NPPl + NPPf + NPPfr + NPPcr

NPP = (1.42 + 1+ 0.14 + 0.95 + 0.26 · 1.42) · NPPl

NPP = 3.88 NPPl

NPPw = 1.42 NPPl

NPPw = 0.366 NPP

where NPP is net primary productivity, and subscripts w, l, f, fr and cr stand for wood, leaf, fruit, fine root and coarse root, respectively.

TABLE 4.Carbon allocation parameters and values used in BBGCMuSo model

Modelling results are area-based, i.e. NPPw is expressed in kg·C·m-2·y-1, while tree‑ring data are tree-based, i.e. BAI is expressed in cm2·y-1. To be able to assess the growth dynamics of the results from the modelling against the tree-ring data we calculated the average annual NPPw and the average annual BAI from all simulation runs (36 in Pokupsko basin and 55 in Spačva basin) and all tree cores for each of the forest areas, respectively. Then we standardized both NPPw and BAI. Standardization is introduced because a direct comparison of NPPw and BAI is not possible without introducing additional uncertainty. To make a comparison without the standardization we would need to calculate the net primary productivity of wood based on the tree core data. This would require the use of allometric functions for estimating wood volume, as well as the use of wood density and wood carbon content values. In addition, not all trees in the plots have been cored, and NPP of the uncored trees would have to be evaluated. All this would introduce additional errors. On the other hand, using standardized values, despite being somewhat more difficult to grasp, circumvents those problems and at the same time keeps the information on growth dynamics.

Standardized values (z-values) of BAI or NPPw, were calculated as:

$\fn_phv&space;\large&space;z=&space;\frac{\left&space;(&space;x_{i}-\bar{x}&space;\right&space;)}{\delta&space;_{x}}$

where  is the variable of interest (the average tree BAI or the average simulated NPPw) in the year  (i = 2000 to 2014) at the given area;  is the overall average of all tree BAI, or NPPw of all simulated forest compartments, at a given area during the entire 15-years long period of interest;  and   is the standard deviation of  during the period of interest (in our case years 2000 to 2014).

Before performing standardization, a Shapiro-Wilk W test for normality was conducted on a series of average annual tree BAI and average annual NPPw estimates for each forest area using procedure swilk in STATA 14 (StataCorp, College Station, TX, USA). Average annual BAI data series were normally distributed (W=0.9791, p=0.9630 for Pokupsko; W=0.9559, p=0.6224 for Spačva). Similarly, average NPPw data series were also normally distributed (W=0.9532, p=0.5756 for Pokupsko; W=0.9699, p=0.8563 for Spačva). Therefore, standardization of the data sets is allowed.

Meteorological data were analysed for the same period as for growth dynamics, from 2000 to 2014. For each location mean annual (October-September) and seasonal (April‑September) air temperature (°C) and precipitation (mm) anomalies were calculated as follows:

T = Ti - Tp

P = Pi - Pp

where T is air temperature anomaly, Ti is mean annual (Oct-Sep) or seasonal (Apr-Sep) air temperature in year i , Tp is mean annual/seasonal air temperature of the investigated period (2000-2014), P is precipitation anomaly, Pi is annual/ seasonal precipitation sum, and Pp is the mean of annual/ seasonal precipitation sums during the investigated period.

RESULTS AND DISCUSSION

Measured Growth Dynamics and Observed Meteorology

Tree-ring data revealed differences in growth dynamics between two investigated oak forests (Figure 4; green circles). Interestingly, at the wetter [19, 45, 46] oak site (Pokupsko), growth decreased in colder years (e.g. 2005-2006, Figure 5), which is in contrast to the common negative response of oak trees’ growth to temperature in spring and summer found in Čufar et al. [47], and to the positive response of oaks to rainy, humid and cloudy conditions during the current year’s summer [16]. Nevertheless, according to Renninger et al. [48], it is highly important to account for groundwater table information when interpreting the response of oak ecosystems to dry conditions. Due to low vertical water conductivity of gleysol soil, forests at Pokupsko site are partly flooded with stagnating water during winter and early spring. During a course of a vegetation season groundwater table is relatively high (Table 1), and therefore we can consider that at this particular site growth is rarely water‑limited, but can be rather sunlight-limited during colder cloudy years. For example, in Pokupsko basin in 2011 there was approx. 300 h more sun compared to 2010, or almost 17%. What is even more important, sun hours were in shortage during May and September of 2010 while the case of 2011 it was exactly the opposite (data for meteorological station Karlovac, http://klima.hr/klima_e.php?id=klima_elementi). Contrary to that, at the drier site (Spačva), growth decreased in warm and dry years (e.g. 2007), which is in line with Čufar et al. [47]. Forests at Spačva site can be considered to be water-limited during warm and dry years, especially after the prolonged drought from the ecological perspective, when groundwater table drops significantly (i.e. more than 1 m below the long-term average for a given month), although partial resupply of soil water reserves occurs laterally. In addition, forests in Spačva basin grow on fertile soil, with good water holding capacity, and are considered to be highly productive. Forests that grow on soil with high nutrient availability tend to have higher aboveground biomass and are found to be more susceptible to drought due to a predisposition to hydraulic failure [49]. At both sites model falsely indicated a drop in growth (z(NPPw)) for the dry 2011 (Figure 4). But in 2012, the growth decrease indicated by the model was evident also in tree core (z(BAI)). The observed reduction in growth was a consequence of two extremely warm and dry years in a row (i.e. 2011 and 2012) and is likely due to the carry-over effect of the drought [35, 50].

Evaluating the Predicting Power of the Model

Model results show some differences in growth dynamics between two investigated locations (Figure 4; red triangles). A strong reduction in simulated growth in years 2003 and 2012, observed at both sites, indicates high model sensitivity to dry conditions and high air temperature, i.e. negative precipitation anomalies and positive temperature anomalies during vegetation period (Figure 5). Similar modelling results at both locations indicate that the model is more sensitive to meteorology than to site‑specific conditions. Although two locations have somewhat different abundances of main tree species, soil characteristics and hydrology, as well as the measured growth dynamics (Figure 4), modelled growth shows similar dynamics (Pearson’s correlation coefficient between NPPw (Pokupsko) and NPPw (Spačva) is 0.563). Pappas et al. [3] obtained similar results when testing process-based LPJ‑GUESS model. Authors concluded that model has a very high sensitivity to photosynthetic parameters (i.e. light correlated parameters) and minor sensitivity to hydrological and soil texture parameters.

Quantitative comparison of the growth dynamics from the two different data sources (the measured tree rings and BBGCMuSo model) reveals a poor agreement (i.e. correlation) for both sites (Figure 6). Table 5 shows the results of statistical evaluation for the model-measurement agreement. The extremely dry year 2003 acted as an outlier, according to Tukey’s definition [51], for NPPw at wetter (Pokupsko) sites (Figure 6, left panel). It seems that a single extremely dry year, such as the year 2003, when strongly negative effects on vegetation productivity at European scale were recorded [52], has not significantly affected tree growth at the investigated sites (Figure 4). The ability of oak trees to overcome a single dry event could be explained by the presence of significant soil water reserves (e. g. records from groundwater monitoring in Spačva show that depth to groundwater in Spačva Basin can fluctuate from 0 to ~5 m, while the average water holding capacity of soils is150 mm·m-1 [46]) and/or large carbohydrate reservoirs (i.e. carbon storage pools in trees). The analysis of remote sensing data also indicates a different response of forests and other vegetation to drought [35], where the results suggest that drought in a given year might negatively affect growth in the consecutive years in case of forests, but not for other vegetation types. This is in line with the logic that due to stress carbohydrate reserves might be depleted because of decreased photosynthesis (due to stomatal closure) and/or increased respiration demand due to excess heat, which then has a legacy effect on the growth in the next year [53].

Significantly decreased growth in 2003 obtained from the model indicates that model routines, describing a biological response to the single drought event, have difficulties with predicting growth dynamics under such conditions. In stress conditions carbon allocation ratios (i.e. proportions of assimilates allocated to different plant pools/organs, as well as mobilization of reserves) change in order for the plant to successfully overcome stress [54]. The shortcoming of fixed carbon allocation ratios, currently implemented in the BBGCMuSo model, might become increasingly pronounced in the case of extreme events. Additionally, BBGCMuSo is a “source-driven” model, meaning that the current photosynthates are immediately allocated to tissues with fixed allocation ratios. According to new findings (e.g. [54]), this logic might not be completely applicable to forests, which can partly explain why the extreme event in a given year might have pronounced effect on plant growth in the forthcoming year(s). These issues are limiting the model to properly predict plant’s response to drought stress. The improvement of modelling results for both sites might be achieved if groundwater table information is used [24].

Residual analysis was performed to find possible sources of model-data discrepancy. According to Figures 7 and 8, the relationship between the studied meteorological variables (data from the previous year and the current year) and the model residuals was not significant for the two study sites. However, the trends, although not statistically significant, might be indicative. Positive/negative precipitation (Figure 7) anomaly in the current year seems to cause over/underestimation of NPPw at both sites, although data variability is high. This means that the role of water availability is more emphasized in the model than in reality. On the other hand, there is a difference between locations in the model performance for the current year with respect to the precipitation anomaly in the previous year (Figure 7). At the wetter (Pokupsko) location negative precipitation anomaly still seems to cause the underestimation of the NPPw by the model, but at the drier site (Spačva) the effect is reversed (negative precipitation anomaly seems to cause overestimation of NPPw by the model). Opposite to the effects of precipitation, positive/negative air temperature (Figure 8) anomaly in the current year seems to have different effects at each location. At the wetter (Pokupsko) location, positive temperature anomaly seems to cause underestimation of NPPw, while at the drier (Spačva) site it causes overestimation. Interestingly, positive air temperature anomaly in the previous year seems to cause underestimation of NPPw at both locations. The observed relations are not statistically significant, as we already emphasized, but are in accordance with the logic of buffered growth response against single drought events due to a large rooting depth of trees.

In a previous study [24] statistical evaluation of the observed Biome-BGCMuSo simulated carbon fluxes against measured eddy covariance data from Pokupsko Basin (Jastrebarsko site) showed much better agreement (see Table 5 in [24]) than in the current study. The explained variance of the observed gross primary production (GPP) and total ecosystem respiration (TER) reached 84% and 83%, respectively. This is in striking contrast to the negligible amount of explained variance (~2%) for the BAI dataset. It is important to note that biogeochemical models like Biome-BGCMuSo are typically calibrated/validated with data-rich eddy covariance based measurements. The application of BAI as validation data in the case of processed based models is rare in the literature. The use of BAI data or NPP estimated from BAI measurements, in Biome-BGCMuSo calibration might help to improve the allocation parameters and thus improve the predictive power of the model. This would be even more important when dynamic allocation will be implemented in the next version of Biome-BGCMuSo (Barcza, personal communication). Multi-objective calibration using both eddy covariance and BAI data (probably with additional variables such as leaf mass, leaf C:N ratio) might provide additional constraints to improve model performance in such cases. The calibration will be a challenging task where sophisticated calibration (e.g. Bayesian [55]) techniques will have to be used.

CONCLUSIONS

The comparison of modelling results with the observed tree-ring data revealed two important model issues related to its predictive power. The first one is the importance of including site-specific conditions (i.e. groundwater table information) with the purpose of enabling the model to be more case-sensitive. At both oak forest locations, Pokupsko and Spačva basins, BBGCMuSo showed poor predictive power in capturing the dynamics obtained from tree‑ring data. Using groundwater table information for modelling in lowland oak forests might improve model results.

The second issue is related to carbon allocation. Fixed carbon allocation ratios, currently used in BBGCMuSo model, do not enable the model to successfully predict plant response to stress conditions (e.g. drought). A dynamic carbon allocation routine might better capture tree stress response and growths dynamics. There is an urgent need to investigate and implement more sophisticated carbon allocation routines in the BBGCMuSo model.

Acknowledgments

The research has been supported by the Croatian Science Foundation (HRZZ UIP-11-2013-2492) and “Spačva” Project (OKFŠ, HŠ 2013-2016). The work of Z.B. and D.H. on the research was funded by the Széchenyi 2020 programme, the European Regional Development Fund and the Hungarian Government (GINOP-2.3.2-15-2016-00028) and Z.B. was also supported by the grant “EVA4.0”, No. CZ.02.1.01/0.0/0.0/16_019/0000803 financed by OP RDE.

We would like to thank two anonymous reviewers for their comments which greatly improved the manuscript.

REFERENCES

1. SÁNCHEZ-SALGUERO R, CAMARERO JJ, GUTIÉRRET E, GONZÁLEZ RF, GAZOL A, SANGÜESA-BARREDA G, ANDREU-HAYLES L, LINARES JC, et al. 2017 Assessing forest vulnerability to climate warming using a process-based model of tree growth: bad prospects for rear-edges. Glob Change Biol 23 (7): 2705-2719. DOI: https://doi.org/10.1111/gcb.13541
2. PETERS EB, WYTHERS KR, ZHANG S, BRADFORD JB, REICH PB 2013 Potential climate change impacts on temperate forest ecosystem processes. Can J For Res 43 (10): 939-950. DOI: https://doi.org/10.1139/cjfr-2013-0013
3. PAPPAS C, FATICHI S, LEUZINGER S, WOLF A, BURLANDO P 2013 Sensitivity analysis of a process-based ecosystem model: Pinpointing parameterization and structural issues. J Geophys Res Biogeosci 118 (2): 505-528. DOI: https://doi.org/10.1002/jgrg.20035
4. FATICHI S, LEUZINGER S, KÖRNER C 2014 Moving beyond photosynthesis: from carbon source to sink-driven vegetation modelling. New Phytol 201 (4): 1086-1095. DOI: https:// doi.org/10.1111/nph.12614
5. HARTIG F, DYKE J, HICKLER T, HIGGINS SI, O’HARA RB, SCHEITER S, HUTH A 2012 Connecting dynamic vegetation models to data – an inverse perspective. J Biogeogr 39 (12): 2240-2252. DOI: https:// doi.org/10.1111/j.1365-2699.2012.02745.x
6. BALDOCCHI DD 2003 Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: past, present and future. Glob Change Biol 9 (4): 479-492. DOI: https:// doi.org/10.1046/j.1365-2486.2003.00629.x
7. BALDOCCHI DD, FALGE E, GU L, OLSON R, HOLLINGER D, RUNNING SW, ANTHONI P, BERNHOFER C, et al. 2001 FLUXNET: A New Tool to Study the Temporal and Spatial Variability of Ecosystem-Scale Carbon Dioxide, Water Vapor, and Energy Flux Densities. B Am Meteorol Soc 82 (11): 2415-2434. DOI: https://doi.org/10.1175/1520-0477(2001)082<2415:FANTTS>2.3.CO;2
8. KUMAR J, HOFFMAN FM, HARGROVE WW, COLLIE RN 2016 Understanding the representativeness of FLUXNET for upscaling carbon flux from eddy covariance measurements. Earth Syst Sci Data Discuss. Manuscript under review. DOI: https://doi.org/10.5194/essd-2016-36
9. BARCZA Z, KERN A, HASZPRA L, KLJUN N 2009 Spatial representativeness of tall tower eddy covariance measurements using remote sensing and footprint analysis. Agr Forest Meteorol 149 (5): 795-807. DOI: https://doi.org/10.1016/j.agrformet.2008.10.021
10. EVANS MEK, FALK DA, ARIZPE A, SWETNAM TL, BABST F, HOLSINGER DKE 2017 Fusing tree-ring and forest inventory data to infer influences on tree growth. Ecosphere 8 (7): e01889. DOI: https://doi.org/10.1002/ecs2.1889
11. BÜNTGEN U, FRANK DC, WILSON R, CARRER M, URBINATI C, ESPER J 2008 Testing for tree-ring divergence in the European Alps. Glob Change Biol 14 (10): 2443-2453. DOI: https://doi.org/10.1111/j.1365-2486.2008.01640.x
12. ESPER J, FRANK DC, BÜNTGEN U, VERSTEGE A, LUTERBACHER J, XOPLAKI E 2007 Long-term drought severity variations in Morocco. Geophys Res Lett 34 (17): . DOI: https://doi.org/10.1029/2007GL030844
13. FRITTS HC 1976 Tree Rings and Climate. Academic Press, London, UK, 567 p
14. CARRER M, URBINATI M 2006 Long-term change in the sensitivity of tree-ring growth to climate forcing in Larix decidua. New Phytol 170 (4): 861-872. DOI: https://doi.org/10.1111/j.1469-8137.2006.01703.x
15. SCHARNWEBER T, MANTHEY M, CRIEGEE C, BAUWE A, SCHRÖDER C, WILMKING M 2011 Drought matters - declining precipitation influences growth of Fagus sylvatica and Quercus robur L. in north-eastern Germany. Forest Ecol Manag 262 (6): 947-961. DOI: https://doi.org/10.1016/j.foreco.2011.05.026
16. HAFNER P, GRIČAR J, SKUDNIK M, LEVANIČ T 2015 Variations in Environmental Signals in Tree-Ring Indices in Trees with Different Growth Potential. PLoS ONE 10 (11): e0143918. DOI: https://doi.org/10.1371/journal.pone.0143918
17. STOJANOVIĆ DB, LEVANIČ T, MATOVIĆ B, ORLOVIĆ S 2015 Growth decrease and mortality of oak floodplain forests as a response to change of water regime and climate. Eur J For Res 134 (3): 555-567. DOI: https://doi.org/10.1007/s10342-015-0871-5
18. WRB 2006 World reference base for soil resources 2006. A framework for international classification, correlation and communication. World Soil Resources Reports No. 103, FAO, Rome, Italy, 128 p
19. MAYER B 1996 Hydropedological relations in the region of lowland forests of the Pokupsko basin (in Croatian with English summary). Rad Šumar Inst 31 (1-2): 37-89
20. SELETKOVIĆ Z 1996 Climate in Pedunculate oak forests (in Croatian with English summary). In: Klepac D (ed) Hrast lužnjak u Hrvatskoj. HAZU - Centar za znanstveni rad Vinkovci i "Hrvatske šume" d.o.o., Vinkovci-Zagreb, Croatia, pp 71-82
21. KALINIĆ M 1975 Soils of forest communities in Spačva basin (in Croatian with English summary). In: Varićak T, Vidaković M, Švagelj D (eds) Sto godina znanstvenog i organiziranog pristupa šumarstvu jugoistočne Slavonije. HAZU Centar za znanstveni rad Vinkovci, Zagreb, Croatia, pp 413-432
22. PILAŠ I 2017 Climatic, hydro-pedologic and geomorphologic characteristics of Spačva basin. In: Final report of the project "Ecological-climate changes and regeneration challenges of Pedunculate oak forests in Spačva basin" (in Croatian). Croatian Forest Research Institute, Jastrebarsko, Croatia, pp 161-195
23. HIDY D, BARCZA Z, HASZPRA L, CHURKINA G, PINTÉR K, NAGY Z 2012 Development of the Biome-BGC model for simulation of managed herbaceous ecosystems. Ecol Model 226: 99-119. DOI: https://doi.org/10.1016/j.ecolmodel.2011.11.008
24. HIDY D, BARCZA Z, MARJANOVIĆ H, OSTROGOVIĆ SEVER MZ, DOBOR L, GELYBÓ G, FODOR N, PINTÉR K, et al. 2016 Terrestrial Ecosystem Process Model Biome-BGCMuSo v4.0: Summary of improvements and new modelling possibilities. Geosci Model Dev 9 (12): 4405-4437. DOI: https://doi.org/10.5194/gmd-9-4405-2016
25. THORNTON PE 2000 User’s Guide for Biome-BGC, Version 4.1.1. URL: ftp://daac.ornl.gov/data/model_archive/BIOME_BGC/biome_bgc_4.1.1/comp/bgc_users_guide_411.pdf
26. PIETSCH SA, HASENAUER H, KUCERA J, CERMÁK J 2003 Modeling effects of hydrological changes on the carbon and nitrogen balance of oak in floodplains. Tree Physiol 23 (11): 735-746
27. CIENCIALA E, TATARINOV FA 2006 Application of BIOME-BGC model to managed forests. 2. Comparison with long-term observations of stand production for major tree species. Forest Ecol Manag 237 (1-3): 252-266. DOI: https://doi.org/10.1016/j.foreco.2006.09.086
28. HLÁSNY T, BARCZA Z, FABRIKA M, BALÁZS B, CHURKINA G, PAJTÍK J, SEDMÁK R, TURCÁNI M 2011 Climate change impacts on growth and carbon balance of forests in Central Europe. Climate Res 47 (3): 219-236. DOI: https://doi.org/10.3354/cr01024
29. MERGANIČOVÁ K, MERGANIČ J 2014 The Effect of Dynamic Mortality Incorporated in BIOME-BGC on Modelling the Development of Natural Forests. J Environ Inform 24 (1): 24-31. DOI: https://doi.org/10.3808/jei.201400273
30. HLASNY T, BARCZA Z, BARKA I, MERGANIČOVÁ K, SEDMAK R, KERN A, PAJTIK J, BALAZS B, FABRIKA M, CHURKINA G 2014 Future carbon cycle in mountain spruce forests of Central Europe: Modelling framework and ecological inferences. Forest Ecol Manag 328: 55-68. DOI: https://doi.org/10.1016/j.foreco.2014.04.038
31. DOBOR L, BARCZA Z, HLÁSNY T, HAVASI Á, HORVÁTH F, ITTZÉS P, BARTHOLY J 2015 Bridging the gap between climate models and impact studies: The FORESEE Database. Geosci Data J 2 (1): 1-11. DOI: https://doi.org/10.1002/gdj3.22
32. BINKLEY D, CROMAC Kjr, BAKER DD 1994 Nitrogen Fixation by Red Alder: Biology, Rates, and Controls. In: Hibbs DE, DeBell DS, Tarrant RF (eds) The Biology and Management of Red Alder. Oregon State University Press, Corvallis. Or (xi) 256 p
33. MARJANOVIĆ H, OSTROGOVIĆ MZ, ALBERTI G, BALENOVIĆ I, PALADINIĆ E, INDIR K, PERESOTTI A, VULETIĆ D 2011 Carbon dynamics in younger stands of pedunculate oak during two vegetation periods (in Croatian with English summary). Sumar list 135 (11-12): 59-73
34. Rulebook on forest management 2015 Official Gazette of the Republic of Croatia 79/15.
35. KERN A, MARJANOVIĆ H, DOBOR L, ANIĆ M, HLÁSNY T, BARCZA Z 2017 Identification of Years with Extreme Vegetation State in Central Europe Based on Remote Sensing and Meteorological Data. South-east Eur for 8 (1): 1-20. DOI: https://doi.org/10.15177/seefor.17-05
36. MIKO S, HASAN O, KOMESAROVIĆ B, ILJANIĆ N, ŠPARICA MIKO M, ĐUMBIR AM, OSTROGOVIĆ SEVER MZ, PALADINIĆ E, MARJANOVIĆ H 2017 Promjena zaliha ugljika u tlu i izračun trendova ukupnog dušika i organskog ugljika u tlu te odnosa C:N (Soil carbon stock change and calculation of trends of total soil nitrogen, organic carbon, and C:N ratio - in Croatian), Study for the Croatian Agency for the Environment and Nature (Contract No. 10-14-1442/79). Zagreb: Croatian Geological Survey, Jastrebarsko: Croatian Forest Research Institute, Zagreb: Agricultural Land Agency, 81p
37. ETHERIDGE DM, STEELE LP, LANGENFELDS RL, FRANCEY RJ, BARNOLA JM, MORGAN VI 1996 Natural and anthropogenic changes in atmospheric CO2 over the last 1000 years from air in Antarctic ice and firn. J Geophys Res 101 (D2): 4115-4128. DOI: https://doi.org/10.1029/95JD03410
38. CHURKINA G, BROVKIN V, VON BLOH W, TRUSILOVA K, JUNG M, DENTENER F 2009 Synergy of rising nitrogen depositions and atmospheric CO2 on land carbon uptake moderately offsets global warming. Global Biogeochem Cy 23 (4): GB4027. DOI: https://doi.org/10.1029/2008GB003291
39. RUNNING SW, ZHAO M 2015 User’s Guide Daily GPP and Annual NPP (MOD17A2/A3) Products NASA Earth Observing System MODIS Land Algorithm, Ver. 3.0 for Coll. 6, 7 October 2015. URL: http://www.ntsg.umt.edu/files/modis/MOD17UsersGuide2015_v3.pdf (15 Nov 2017).
40. OSTROGOVIĆ MZ 2013 Carbon stocks and carbon balance of an even-aged Pedunculate Oak (Quercus robur) forest in Kupa river basin. Ph.D. thesis, University of Zagreb, Faculty of Forestry, 130 pp. (in Croatian with English summary). URL: https://dr.nsk.hr/en/islandora/object/sumfak%3A795 (Nov 2017).
41. STOKES MA, SMILEY TL 1968 An Introduction to tree-ring dating. Chicago: University of Chicago Press, USA, 73 p
42. HOLMES RL 1983 Computer-assisted quality control in tree-ring dating and measurements. Tree-Ring Bull 43: 69-78
43. GRISSINO-MAYER HD 2001 Evaluating crossdating accuracy: a manual and tutorial for the computer program COFECHA. Tree-Ring Res 57: 205-221
44. BIONDI F, QEADAN F 2008 A theory-driven approach to tree-ring standardization: defining the biological trend from expected basal area increment. Tree-Ring Res 64 (2): 81-96. DOI: https://doi.org/10.3959/2008-6.1
45. PERNAR N, BAKŠIĆ D 2005 The soil of floodplain forests. In VUKELIĆ J (ed) Floodplain forests in Croatia, Akademija šumarskih znanosti & Hrvatske šume & Grad Zagreb, Zagreb, Croatia, pp 71-85.
46. VRBEK B, PILAŠ I, PERNAR N 2011 Observed climate change in Croatia and its impact on the hydrology of lowlands. In: Bredemeier M, Cohen S, Godbold Dl, et al. (eds) Forest management and the water cycle: An ecosystem-based approach, Springer Netherlands, Dordrecht, Netherlands, pp 141-162. DOI: https://doi.org/10.1007/978-90-481-9834-4_8
47. ČUFAR K, GRABNER M, MORGÓS A, CASTILLO EM, MERELA M, DE LUIS M 2014 Common climatic signals affecting oak tree-ring growth in SE Central Europe. Trees 28 (5): 1267-1277. DOI: https://doi.org/10.1007/s00468-013-0972-z
48. RENNINGER HJ, CARLO N, CLARK KL, SCHÄFER KVR 2014 Physiological strategies of co-occurring oaks in a water- and nutrient-limited ecosystem. Tree Physiol 34 (2): 159-173. DOI: https://doi.org/10.1093/treephys/tpt122
49. GESSLER A, SCHAUB M, MCDOWELL NG 2017 The role of nutrients in drought-induced tree mortality and recovery. New Phytol 214 (2): 513-520. DOI: https://doi.org/10.1111/nph.14340
50. MOLEN MK, DOLMAN AJ, CIAIS P, EGLIN T, GOBRON N, LAW BE, MEIR P, PETERS W, et al. 2011 Drought and ecosystem carbon cycling. Agr Forest Meteorol 151 (7): 765-773. DOI: https://doi.org/10.1016/j.agrformet.2011.01.018
51. TUKEY JW 1977 Exploratory Data Analysis. Addison-Wesley. ISBN-13: 978-0201076165, 688 p
52. CIAIS PH, REICHSTEIN M, NICOLAS V, GRANIER A, OGÉE J, ALLARD V, AUBINET M, BUCHMANN N, et al. 2005 Europe-wide reduction in primary productivity caused by the heat and drought in 2003. Nature 437: 529-533. DOI: https://doi.org/10.1038/nature03972
53. ZHAO J, HARTMANN H, TRUMBORE S, ZIEGLER W, ZHANG Y 2013 High temperature causes negative whole-plant carbon balance under mild drought. New Phytologist 200 (2): 330-339. DOI: https://doi.org/10.1111/nph.12400
54. HAGEDORN F, JOBIN J, PETER M, LUSTER J, PRITSCH K, NICKEL UT, KERNER R, MOLINIER V, et al. 2016 Recovery of trees from drought depends on belowground sink control. Nat Plants 2: 16111. DOI: https://doi.org/10.1038/nplants.2016.111
55. VAN OIJEN M, ROUGIER J, SMITH R 2005 Bayesian calibration of process-based forest models: Bridging the gap between models and data. Tree Physiology 25 (7): 915-927. DOI: https://doi.org/10.1093/treephys/25.7.915
56. MORIASI DN, ARNOLD JG, VAN LIEW MW, BINGNER RL, HARMEL RD, VEITH TL 2007 Model Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations. T ASABE 50 (3): 885–900. DOI: https://doi.org/10.13031/2013.23153