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

(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, HR10450 Jastrebarsko, Croatia


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 longterm 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 earlyand latewood formations.The underlying idea behind this research was to test modelling performance by using treering 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.

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 narrowleaved 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.
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][27][28][29][30].Major improvements in BBGCMuSo include introducing a figURe 1. Geographical location of the research areas, Pokupsko Basin (west) and Spačva Basin (east), and the meteorological stations located in Jastrebarsko and Gradište.multilayer soil module with the possibility of using groundwater table information, management module, new plant pools, respiration acclimation, CO 2 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 CO 2 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 CO 2 concentration, nitrogen deposition and management.Finally, the normal run is done using current meteorology, CO 2 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 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 sitespecific 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].Furthermore, we used three optional input files: atmospheric CO 2 , nitrogen deposition, and management file.Atmospheric CO 2 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 figURe 2. Distribution of the number of model runs (simulations) and the number of tree cores according to age classes for forests of Pokupsko and Spačva basins.Note that each simulation is made at the forest compartment level.

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.
Tree cores were taken, on average, from 9.6 dominant and co-dominant trees per plot (min.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: 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.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 cm 2 •y -1 .To be able to assess the growth dynamics of the results from the modelling against the treering 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: , 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 figURe 3. Example of the sampling plot layout (yellow circles) within MODIS 1km pixels (red parallelograms).White lines mark borders of forest compartments (labels in italic) that are part of the management unit "Slavir", part of Spačva Basin.

Allocation parameter parameter value Reference
new fine root C : new leaf C 0.95 [52] new fruit C : new leaf C 0.14 [52] new woody stem C : new leaf C 1.42 [52] new coarse root C : new woody stem C 0.26 [54] TaBle

Carbon allocation parameters and values used in BBGCMuSo model
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 NPP w 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.

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 longterm 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(NPP w )) 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 NPP w (Pokupsko) and NPP w (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 NPP w 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 is >150 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 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).Multiobjective 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.We would like to thank two anonymous reviewers for their comments which greatly improved the manuscript.

figURe 4 .
figURe 4. Measured and simulated growth dynamics during the period 2000-2014 in two distinct locations of pedunculate oak forests (Pokupsko and Spačva Basin) based on standardized values z(BAI) from tree cores, and z(NPPw) from BBGCMuSo model.The trends and the corresponding equations for measured and modelled z-values are also shown.

figURe 5 .
figURe 5. Mean annual (October from the previous year -September of the current year) and growing season (April-September from the current year) air temperature (°C) and precipitation (mm) anomalies during the period 2000-2014 in two distinct locations of pedunculate oak forests (Pokupsko and Spačva Basin).

figURe 6 .
figURe 6. Assessing the predictive power of the model for two distinct locations of pedunculate oak forests at Pokupsko and Spačva Basin with standardized values (z-values) of BAI and NPPw.

2 =figURe 7 .
figURe 7. The correlation of residuals with seasonal (April-September) precipitation anomaly during the previous and the current year for Pokupsko (a) and Spačva Basin (b).

figURe 8 .
figURe 8.The correlation of residuals with seasonal (April-September) temperature anomaly during the previous and the current year for Pokupsko (a) and Spačva Basin (b).

TaBle 2 .
Parameter values adjusted to site-specific conditions.
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.

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

TaBle 5 .
Performance statistics (based on z-values) of Biome-BGCMuSo compared with the observed tree ring data.