Dynamic height growth models for highly productive pedunculate oak (Quercus robur L.) stands: explicit mapping of site index classification in Serbia

Abstract

Key message

We applied the generalized algebraic difference approach (GADA) to develop dynamic models of height growth for pedunculate oak (Quercus robur L.) in Serbia. According to the dominant heights, the studied region comprises some of Europe’s most productive sites for pedunculate oak. Therein, we have generated a map showing the current site index class of stands. Such a map could be used to enhance forest management and evaluate climate change impacts.

Context

Although sustainable forest management requires reliable prediction of forest site productivity, such indicators are currently unavailable for pedunculate oak sites in Serbia. The site index (SI) curves represent the most commonly used indirect scale for site productivity classification. The dynamic equations derived by the Generalized Algebraic Difference Approach (GADA) are the state-of-the-art approach in growth modeling, but they have not been widely applied for studying the height dynamics of pedunculate oak.

Aims

The main objectives of this study were to develop the first dynamic site index curves for pedunculate oak in Serbia and subsequently to provide stand-level maps with predicted site indices.

Methods

We have tested five flexible polymorphic equations with variable asymptotes derived by the GADA approach. Models were calibrated using artificially established growth trajectories obtained from 3636 detailed temporary sample plots. The selection of the most suitable model was accomplished according to (1) quantitative measures of goodness of fit, (2) the analysis of residual scattering, and (3) the biological plausibility of obtained height growth curves.

Results

After correcting the error terms with a continuous first-order autoregressive structure and conducting a three-stage performance analysis, the GADA dynamic site index model derived from the Hossfeld base equation shows the best overall properties. Insight into the oscillations of relative error suggested that 100 years is the most suitable age for site index referencing. Comparison with existing height growth models revealed greater flexibility and a considerably better representation of the height growth dynamic of pedunculate oak in the studied region. Additionally, we have produced a spatially explicit map showing the expected SI₁₀₀ for 1907 stands with pedunculate oak within 22 management units.

Conclusion

Dynamic SI-curves based on GADA will serve forest practitioners to update management plans and serve as a reference point for benchmarking the impact of climate change and for developing adaptation strategies. The utilized approach allowed unbiased estimation of SI₁₀₀ across all age classes so that the results could be mapped at a broader scale. This study provides the second known application of the dynamic model for pedunculate oak in Europe but the first that includes some of the most productive sites in the species distribution range.

1 Introduction

Serbia is a hotspot of climate change in Europe (Kržič et al. 2011; Đurđević, V.; Vuković, A.; Vujadinović-Mandić 2018; Erić et al. 2021). In the last two decades, almost the whole country has experienced an increase of mean maximum temperature during the summer season above 2.5 °C, while given projection up to the year 2100 under the RCP8.5 scenario assumes increasing of 4.7 °C in the June–August period (Đurđević, V.; Vuković, A.; Vujadinović-Mandić 2018). Forest management thus has to focus on developing strategies to mitigate the environmental and socio-economic consequences of climate change on forests. To secure the sustainability of wood production accompanied by all other ecosystem services (Löf et al. 2016; Brockerhoff et al. 2017), adaptation strategies should arise from the interpretation of foreseen environmentally driven shifts in (1) species occurrence and (2) growing dynamics of trees and stands (Pretzsch et al. 2012; Mette et al. 2013; Dyderski et al. 2018; Castaño-Santamaría et al. 2019). The present distribution of pedunculate oak (Quercus robur L.) forests in Serbia is tied to belts around large rivers, where underground water abundantly compensates for the lack of rainfalls necessary for the trees to grow. Decreasing river water levels due to drought events thus directly affects the growth of pedunculate oak. Oak decline and dieback in this area have been registered for a long time now (Dubravac and Dekanić 2009; Medarević et al. 2009; Bauer et al. 2013), and have recently been directly related to the decrease in water availability (Stojanović et al. 2013, 2015; Kostić et al. 2021b, a).

Although covering a relatively small total area of 32,400 ha (Banković et al. 2009), stands with pedunculate oak represent one of the most valuable Serbian forest areas. Around 90% of these forests (29.000 ha) are growing alongside the riverbanks of the Sava River in the Ravni Srem region (Rađević et al. 2020), situated in the north-western part of the country, in Vojvodina province, where forest occupies just above 7% of total area (Banković et al. 2009). Due to the high economic and ecological value of the floodplain forest of pedunculate oak, they have been studied from many different aspects (Bobinac 2008a, b; Ivanišević and Knežević 2008; Medarević et al. 2009; Bauer et al. 2013; Stojanović et al. 2013, 2015; Rađević et al. 2020; Kostić et al. 2021b, a). Forest management has a long tradition in this area starting in the year 1755 with the first governance and protection guidelines (Bobinac 2008a). Driven by the historical circumstances, the several different phases of utilization and regeneration of pedunculate oak stands in the past were documented (Medarević et al. 2009; Rađević et al. 2020). The present renewal procedure is carried out by regenerative cutting on larger areas, in combination with additional acorn sowing and intensive forestry tending and protection measures in the youngest stages of development.

However, there is a severe lack of reliable productivity indicators for pedunculate oak stands that are essential for prescribing suitable management targets and appropriate silviculture measures. The most commonly used indirect scale for site productivity classification is obtained by sectioning the height-age oscillation range with the desired number of the site index curves (Skovsgaard and Vanclay 2013; Brandl et al. 2018). The development of site index curves until now has been done only for beech stands in the eastern part of Serbia (Stajić et al. 2016) and subsequently expanded to include beech stands in the central region (Stajić et al. 2021), but not for other domestic tree species including the pedunculate oak. From neighboring countries, we found Špiranec’s (1975) height growth model for pedunculate oak in Croatia using the guide curve method. The utilized dataset was mainly collected in the Spačva forest area along the Sava River, with almost identical growing conditions to those in Ravni Srem, so these two regions form a unique cross-state complex of floodplain pedunculate oak forests. Contrary to other important tree species in Europe, modeling the height growth for pedunculate oak based on the generalized algebraic difference approach (GADA) has not been widely conducted. Currently, a dynamic height growth model is only available for low-productivity sites (Barrio-Anta and Diéguez-Aranda 2005). Therefore, our objective was to complement the height growth curves of pedunculate oak using state-of-the-art GADA models by covering a wider range of growth conditions including growth dynamics in medium to highly productive conditions.

Reliable assessment of the expected height growth trajectories and resulting site indices depends on the suitability of a modeling methodology and the integrity of underlying data. Ideally, the growth modeling is based on long-time series that are evenly distributed across a complete spectrum of age classes and site conditions (Brandl et al. 2018). However, varying properties of the utilized datasets and their impact on modeling outcomes are a matter of continuous discussions, that often lead to rather contrasting conclusions. Socha et al. (2019) reported the four potential data origins for site index calibration in their work, although it seems that data from stem analysis of unsuppressed trees grown at preselected sites were most frequently used as a data source for height growth modeling (Barrio-Anta and Diéguez-Aranda 2005; Bravo-Oviedo et al. 2008; Martín-Benito et al. 2008; Nunes et al. 2011; Seki and Sakici 2017; Matisons et al. 2018; Castaño-Santamaría et al. 2019; Socha et al. 2020). Due to gathering complexity, the site index calibration is slightly less often based on the height–age data pairs from successive surveys of the permanent sample plots (Cieszewski et al. 2007; Nord-Larsen et al. 2009; Álvarez-González et al. 2010; Sharma et al. 2019; Manso et al. 2021). Besides, some authors have also tested the fusion of different data sources (Sharma et al. 2011; Socha and Tymińska-Czabańska 2019) or relied on their combinations in the fitting and validation phase (Martín-Benito et al. 2008; Nunes et al. 2011).

Although most desirable data for developing site index curves come from stem analysis or re-measured permanent sample plots (Cieszewski et al. 2000), this data can still exhibit severe deficiencies. Despite high material and/or time costs, these datasets potentially underrepresent the growth trends at lower productive sites (Nord-Larsen et al. 2009). Namely, locations for setting up permanent experimental plots or cutting the trees for growth analysis are often selected in favor of easily approachable, well-managed stands on better sites (Sharma et al. 2011). In other circumstances, the sampling design was biased toward low-productivity sites causing unexpectedly high asymptotes (Martín-Benito et al. 2008; Socha and Tymińska-Czabańska 2019). The unbalanced sampling along a site quality gradient is sometimes also caused by truncated growth series at both ends (Adame et al. 2006; Nord-Larsen et al. 2009). Bias in the models can be introduced by data originating from the trees that were not in dominant status over their whole lifespan (Cherubini et al. 1998; Bravo-Oviedo et al. 2007).

The development and spread of the Generalized Algebraic Difference Approach- GADA (Cieszewski and Bailey 2000) resolved the substantial issues from former techniques of creating site curves index. GADA solves the basic growth equation after expansion with the two parameters related to theoretical variable X. Unlike the preceding base-age-specific models, reference points are treated locally as site-specific parameters, whereas all other parameters are estimated globally. Mathematically sound derivations of GADA models are finished by self-referencing, which resulted in base-age invariant and parsimonious growth models (Cieszewski 2004). While the complexity and robustness of the dynamic model guarantee fitting performances, accompanying path-invariance properties ensure applicability for both forward and backward predictions at any time intervals. Such desirable properties allow generalized equations to simultaneously reach polymorphism with different asymptotes, securing the ability to assess specific growth trends across the complete range of site conditions (Cieszewski and Bailey 2000; Cieszewski and Strub 2008). Therefore, in the last decades, GADA derivations have been increasingly used for the development of site index models, as well as for describing diameter growth (Sharma et al. 2017) stand basal area growth (Barrio-Anta et al. 2006; Castedo-Dorado et al. 2007) and tree biomass of the individual trees (Tang et al. 2017).

Producing spatially explicit maps represents an effective way to disseminate scientific results to practitioners in forestry and an efficient tool to support the decision-making process (Brus et al. 2012; Brandl et al. 2014; Lombardi et al. 2016; Parresol et al. 2017; Baumbach et al. 2019; Schumacher et al. 2020). Besides its importance for ongoing forest management, the expected site indices maps are suitable for assessing the effects of climate change on heights and related growth rates in the future. In that sense, the current state represents a valuable reference point for comparison with the values predicted by climate-sensitive growth models like those presented by Albert and Schmidt (2010), Sharma et al. (2015), and Socha et al. (2021). Moreover, the outputs of the site index model can be directly utilized for training the machine learning algorithms, which based on the selected environmental predictors, are capable of mapping the current height occurrences, as well as expectations under various climate change scenarios (Weiskittel and Kershaw 2010; Castaño-Santamaría et al. 2019). The risk analysis based on foreseeing the future growth dynamics in modified environmental conditions is an important part of management adaptation strategies (Hanewinkel and Peyron 2014; Sperlich et al. 2020; Martinez del Castillo et al. 2022). Also, with such kind of productivity maps we aim to provide ground for a better understanding of the present variations of pedunculate oak’s radial increment (Stojanović et al. 2013, 2015; Kostić et al. 2021b, a).

Having in mind the previously mentioned, our specific objectives in this study were the following:

• (i) to calibrate the first state-of-the-art dynamic site index equation (SI) for pedunculate oak based on the empirical data collected across the wide range of environmental conditions, from low to medium and high productivity including some of the most favorable sites where pedunculate oak reaches its maximal heights, and.

• (ii) to create a spatially continuous map with the expected site index for 1907 pedunculated oak stands within 22 management units in Ravni Srem.

The results will help forest managers to adjust ongoing management with the growth dynamics of species. At the same time, an explicit mapping of the SI classification will serve as a reference point for the evaluation of climate change's impact on height growth and one of the required tools to develop adaptation strategies.

2 Material and methods

2.1 Study area

The distribution of pedunculate oak (Quercus robur L.) in Serbia concentrates in the Sava River basin located in the Ravni Srem region (West Serbia), occupying the 42.000 ha of forests, which are managed by the public forest enterprise (PE) “Vojvodinasume” (Fig. 1). Ravni Srem is physically divided into the “lower Srem” (grey box at Fig. 1) and the “upper Srem” (black box at Fig. 1). In addition to pedunculate oak, herein is reported the presence of Turkey oak (Quercus cerris L.), ash (Fraxinus angustifolia Vahl), willows (Salix alba L. and Salix cinerea L.), poplar (Populus alba L.), black alder (Alnus glutinosa L.), elms (Ulmus minor Mill. and Ulmus effusa Willd.) and hornbeam (Carpinus betulus L.).

Fig. 1

The studied forest complex in the Ravni Srem region in the north-western part of Serbia, with spatial positions of 25 management units (differently colored and labeled with unique state-wide ID) and the closest meteorological stations (blue symbols)

The investigated terrain is flat with elevations mostly below 100 m a.s.l. The inner part of the forest complex is protected from river flooding by the embankment built alongside the riverside around 1930 (Stojanović et al. 2013). According to Ivanišević and Knežević (2008), the most represented soil types are humofluvisol, humogley, and eutric cambisol, occupying 38%, 33%, and 17% of the total forest area. The remaining 11% of stands are growing on pseudogley-gley, eugley, pseudogley, alluvial-diluvial soil, and chernozem.

The climatic records in the monthly resolution were available from 8 surrounding meteorological stations shown in Fig. 1. The stations enclose the area between 45.133º–44.467º latitude and 19.233º–22.283º longitude, covering the altitude range from 79 to 113 m a.s.l. Based on the climate time series (1949–2021), the observed climatic conditions are homogeneous across the studied region. The mean annual temperature ± standard deviation was 11.4 ± 0.5 °C with 0.6 ± 0.5 °C in January and 22.0 ± 0.4 °C in July. The mean annual precipitation was 659 ± 33 mm m⁻² with the highest amounts during May–June (73 ± 7.5 mm m⁻²). Besides rainfall, the underground water represents an additional water source for the plants in the alluvial plane due to the vicinity of rivers and irrigation canals. Distribution and growth of pedunculate oak are vulnerable to fluctuation in the water regime (Čater and Batič 2006; Cedro 2007; Pretzsch et al. 2013; Stojanović et al. 2015). Soil properties can thus significantly modify the reaction to water scarcity (Kostić et al. 2021a).

2.2 Data collection

PE “Vojvodinasume” provided the database for the height growth modeling based on the enterprise’s regular stand inventories. The shade intolerance predetermined the even-age structure of pedunculate oak stands, while exceptional interest in pedunculate oak forests affected current management plans to contain reliable stand age records. Generally, the stand’s age information comes from historical records that were checked and updated by analyzing the trees felled down during silviculture interventions. From the network of circular temporary sample plots (TSP) that underlie stand inventories, we extracted and used 3711 height-age pairs in total. To simulate the age spatial variability, the mean stand age at each TSP was modified by adding values up to ± 2 years randomly drawn with uniform probabilities.

The top height for each sample plot was calculated as the mean squared height of trees with the largest diameter, which is proportional to the 100 thickest trees per hectare (García 1998; Ochal et al. 2017; Manso et al. 2021). The height-age relationship was then statistically and visually analyzed for erroneous data. All outlier values for a given age exceeding the first and third quantile for more than ± 1.5 times the interquartile range were excluded from further analysis. After excluding 1.5% of the total sample, the remaining 3656 data pairs were used for establishing the absolute age–height variation spectrum for pedunculate oak in the Ravni Srem region (Fig. 2). To replicate the expected height growth trajectories, the artificial chronosequences were assembled following the procedure proposed by Socha and Tymińska-Czabańska (2019). Firstly, the percentile threshold values in 5% steps were calculated for each year along the time axis. Subsequently, the single points within the same percentile ranges were matched across ages except for two marginal classes. This resulted in 19 artificial age–height chronosequences (Fig. 2). The remaining data from the percentile classes were used as a validation dataset. Both parameterization and validation data sets are structured in 20-year-wide age classes (Table 1). According to frequency distribution (N), it is evident that both samples show no bias in distribution, nor unwonted patterns in minimal (min), mean ($\overline{X }$), maximal (max) values, or standard deviations (SD).

Fig. 2

Top heights over stand age calculated for temporary sample points (TSP) and 19 artificial heights growth trajectories (points + lines). Different percentile ranges are colored according to the given scale

Table 1 The number of observations (N), the mean age ($\overline{{\varvec{X}} }$), standard deviation (SD), minimal (min), and maximal (max) age in age classes of datasets used for calibration and verification of height growth models

By overlapping the results of detailed pedological and phytocoenological research conducted along an existing ecological gradient of Ravni Stem, a comprehensive list of specific forest types was defined (Jović et al. 1994) and afterward updated (Ivanišević and Knežević 2008). Even though established typological classification omits detailed productivity research, it is noticed that some of the described forest types demand specific silviculture and regeneration approaches (Bobinac 2008b). From the available database, we have succeeded in identifying 17 of 19 described forest types where pedunculate oak appeared as the main or admixed species. In Table 2, we provide the list of herein studied forest types alongside the area they cover, the number of delineated stands and associated the TSP.

Table 2 The investigated forest types, the area they cover, the number of delineated stands, and associated the temporary sample points (TSP)

2.3 GADA fitting, calibration, and validation

In order to adequately represent the natural properties of different height growth trajectories, we have applied the functions derived based on the generalized algebraic difference approach- GADA (Cieszewski 2001). Unlike the basic algebraic difference approach-ADA (Bailey and Clutter 1974), the given generalization method allows for two selected parameters to vary between sites, securing the possibility that the expanded model simultaneously achieves the desired properties of polymorphism and multiple asymptotes. The approach transforms basic equations into adjustable dynamic forms, where the theoretical site variable ${X}_{0}$ is used to relate the local (site-specific) parameters to a vector of global (fixed) parameters ($\widehat{\beta }$) and just one varying parameter. The system of height-age curves for the different site qualities is created by letting the expected height (${H}_{0}$) at base age (${T}_{0}$) vary in the desired frame of site index values. Therefore, the calculation of the height ($H$) at a given age ($T$) from any GADA model $g()$ can be schematically shown as $H=g\left({H}_{0},{T}_{0},T \right| \widehat{\beta })$. To select the GADA model that best describes height growth on a range of site quality, we have tested five dynamic forms of the equation with two site-specific parameters. The models used in this study (M1- M5) are shown in Table 3.

Table 3 Base models and tested dynamic equations derived from the generalized algebraic difference approach (GADA) (M1–M5)

The local and global parameters were estimated simultaneously as a part of nonlinear fixed-effects models by using the generalized non-linear least squares (GNLS) of the NLME package (Pinheiro et al. 2021). In the time series where the past values predetermine the succeeding one, the autocorrelation in residuals is expected. To handle interdependence in longitudinal data during the fitting phase, we have done autoregressive modeling of the residuals. The autocorrelation was modeled by the CAR(x), which is a continuous-time autoregressive error structure and recommended to account for irregular and unbalanced data (Zimmerman and Núñez-Antón 2001; Seki and Sakici 2022; Zimmerman and Nunez-Anton 2001; Seki and Sakici 2022). The homoscedasticity or any other pattern of dependence within model error was examined by plotting standardized residuals versus (1) predicted values of height, (2) residuals at different time lags, (3) given age, and (4) observed site index.

To avoid possible subjective selection of the most suitable model for representing the height growth of pedunculate oak, a three-step procedure was used to compare the performance of the GADA candidates. In the first phase, the goodness-of-fit was evaluated using these four statistics: estimation bias ($\overline{e }$), root mean square error (${\text{RMSE}}$), adjusted coefficient of determination (${{R}_{{\text{adj}}}}^{2}$), and the Akaike information criterion (${\text{AIC}}$). The randomness of residual scattering is assessed with Durbin-Watson statistics (d-w). Since the quality of fit does not simultaneously reflect the quality of prediction, a validation study using independent data is irreplaceable (Nunes et al. 2011; Manso et al. 2021). In the second step, we assessed the prediction performance of the models by calculating the estimation bias (Bias), the mean absolute difference (${\text{MAD}}$), and model estimation efficiency (${\text{MEF}}$). The last part of the adopted selection procedure implicates checking the biological plausibility of obtained height growth curves. For that purpose, we analyzed the practical acceptability of estimated upper asymptotes (Asy.), together with the current annual height increment (CAI) at the time of culmination (Age).

The direct application of the GADA models is sensitive to the choice of the age (${T}_{0}$) to which the site index will be referenced. For the proper selection of reference age that will reduce the uncertainties of predictions, it is necessary to follow objective criteria rather than rely on individual judgment. Therefore, we have followed the approach for base age selection proposed by Diéguez-Aranda et al. (2005). The applied method involves determining the reference ages that minimize the relative error.

Furthermore, using the selected height growth model, we have created the operational map that shows expected heights in reference age for all studied stands with pedunculate oak in Ravni Srem. To process the data and derive the desired outcomes, we have written a devoted script in the R environment (Team R Core 2020).

3 Results

3.1 Calibration and validation of GADA models

During the calibration of all five tested GADA models, the applied optimization algorithm each time gave proper convergence. Table 4 summarizes the parameter estimates obtained after expanding the error terms with a continuous first-order autoregressive structure (CAR1). Based on given t-values, all presented parameters were found to be highly statistically different from zero with a p value < 0.001 (***) or p < 0.01 (**).

Table 4 Parameter estimates of the fitted models (M1–M5) and their standard errors, t values and significance (p)

Table 5 shows three groups of indicators used for the evaluation of the model’s performance and representativeness: (i) fitting statistics, (ii) validation statistics, and (iii) growth elements. Fitting statistics suggest that studied models achieved a good fit to the artificially established growth trajectories, each explaining at least 96% of variations in the empirical dataset. The Durbin–Watson statistics (d-w) confirmed that autocorrelation problems in residuals were successfully solved, generating similar but rather high AIC values, especially for M2-M5 models, while M1 had a notably lower AIC coefficient. Models’ predictions versus calibration data showed the absence of any bias $(\overline{{\text{e}} })$. Although RMSE varies in the relatively narrow range, its values suggest that the M1-M3 equations provide slightly lower prediction errors in contrast to models M4 and M5.

Table 5 Indicators of the model’s goodness-of-fit

The analyzed statistical parameters in the validation phase (ii) have quite similar (Bias) or identical values (MAD, MEF) for each of the tested models. These results suggest that the GADA functions gave precise and accurate predictions for the independent dataset of top heights. For the growth elements, indicators were used to reflect the biological realism of predictions for the sites with the most favorable growth conditions (SI₁₀₀ = 38 m). To evaluate the asymptotic behavior of studied models (Asy.), the expected height is estimated for an age of 200 years. Statistics in the growth elements complement fitting and validation properties because mathematical models have to guarantee practical usefulness through a realistic depiction of inflection point occurrences. The selection of the most appropriate model has to represent a compromise between the statistical parameters and growth realism (Barrio-Anta and Diéguez-Aranda 2005), especially when statistics fail to filter out the best candidate as shown by the analyses of growth parameters (Socha et al. 2020).

3.2 Model selection

According to the presented results, the M5 and M4 estimate the expected heights slightly above empirical records from the oldest stands, so the application of these models cannot be considered adequate. The other two GADA models, M2 and M3, were excluded from further consideration since they forecasted the first-year culmination of the current annual increment (CAI), which is biologically unrealistic. Eventually, the only reasonable prediction of CAI culmination age (11) and values (1.13 m), accompanied by a relatively high yet realistic asymptote (40.4 m), is delivered from the M1 model. Thus, this dynamic GADA model with two site-specific parameters was chosen as the most appropriate for developing the site index curves for pedunculate oak in the Ravni Srem region (Fig. 3). The selected equation M1 with substituted parameters has the following form:

Fig. 3

Development of top height (a) and current annual increment (CAI) (b) over stand age with model M1 of the generalized algebraic difference approach (GADA) for 9 site index classes (SI₁₀₀) from red to green. Empirical height measurements from temporary sample points (TSP) for model fitting are displayed in grey points

$$H={H}_{0}\cdot \frac{{T}^{1.690559 }\left({T}_{0}^{1.690559 } {X}_{0}+{e}^{9.105472}\right)}{{T}_{0}^{1.690559 }\left({T}^{1.690559 }{X}_{0}+{e}^{9.105472}\right)}$$

$${X}_{0}={H}_{0}-22.21846+\sqrt{{({H}_{0}-22.21846)}^{2}+\frac{2{H}_{0}{\cdot e}^{9.105472 }}{{T}_{0}^{1.690559}}}$$

Predicted site index curves based on model M1 over observed TSP height values (grey points) provide ecologically meaningful properties such as polymorphism, multiple asymptotes, and sigmoid growth patterns (Fig. 3). Also, the selected model reliably depicted the anticipated shifts in growth rate as the site's quality declines. The corresponding CAI curves are shown in the lower part of the plot. Presented growth trajectories were created for anticipated top heights in the reference age of 100 years by dividing the entire heights spectrum from 22 m to 38 into 2 m site indexes. Compared to the sites with the most favorable growth conditions, the inflection point at low-quality areas occurs 11 years later, at the age of 22 years, with a CAI of 0.36 m. The maximal predicted height is approximately 14 m lower compared to the 40.4 m that is expected at SI₁₀₀ = 38 m.

Graphical examination of error scattering across the range of predicted top heights (Fig. 4a) shows slightly increased variation for values above 20 m, which can be attributed to higher replication of the sample in this range. No temporal autocorrelation was found between standardized residuals and own values on a first lag (Fig. 4b). In addition, the symmetrical dispersion of residuals is registered across different ages (Fig. 4c) and predicted site index (Fig. 4d).

Fig. 4

Distribution of residuals for predicted heights (a), first lag residuals (b), temporary sample points (TSP) age classes (c), and site index classes (SI₁₀₀) (d)

To determine the optimal reference age, the change of relative error (RE%) for the whole range of ${T}_{0}$ was tested (Fig. 5). Among several points with evidently low RE% amounts, the vertical dotted line intersected the value obtained for the base age of 100 years. The occurrence of a relatively small prediction error for a reference age of 100 years has confirmed the eligibility of the nominally used base age (Fig. 6). Such selection completely coincides with the conventionally accepted reference age for species with longer rotation periods.

Fig. 5

Relative errors (RE%) in top height prediction for different base ages and the number of observations. The vertical dotted line represents reference age 100 for site index classification

Fig. 6

Comparison of modeled top heights based on model M1 of the generalized algebraic difference approach (GADA) for 9 site index classes (SI₁₀₀) with Špiranec’s (1975) height growth models for pedunculate oak in Croatia. The solid black line represents relative site index classes (I, II, III) of Špiranec’s (1975), where dotted lines designate lower and upper thresholds of classes

3.3 Site index mapping

Based on the average value of age and heights from TSP, we have created maps with expected SI₁₀₀ for each of the 1907 studied stands. Two main parts of this map are displayed in Fig. 7, where the analyzed stands are painted according to the color scale given as legend. The top map shows the expected values of site indices in the eastern (lower) Srem area, which are reasonably lower compared to the SI₁₀₀ ascertained for the stands in the western part. The determined site indices do not exhibit any significant correlation with stand age. However, certain discrepancies for site indices within different parts of the same stands are discovered. This fact may indicate that the current stand delineation does not thoroughly reflect the existing differences in the productivity potential of sites.

Fig. 7

Mapping of SI₁₀₀ for stands with pedunculate oak in eastern (a) and western (b) part of Ravni Srem based on model M1 using the generalized algebraic difference approach (GADA)

4 Discussion

Herein we have compared the performance of five dynamic GADA formulations that meet the crucial requirement of polymorphism and multiple asymptotes. Calibration was done with suggested procedures for correction of empirical autocorrelation of the residuals (Barrio-Anta and Diéguez-Aranda 2005; Bravo-Oviedo et al. 2008). For the selection of the most appropriate height growth model, some authors rely more on goodness-of-fit criteria and validation statistics (Adame et al. 2006; Diéguez-Aranda et al. 2006; Sharma et al. 2011), while others underline the importance of the model prediction realism, solely or in a compromise with statistical criteria (Barrio-Anta and Diéguez-Aranda 2005; Nunes et al. 2011; Socha et al. 2020). The good values of used statistical parameters indicate the overall suitability of tested GADA derivations of Hossfeld IV, Chapman-Richards, Lundqvist-Korf, and Weibull growth functions. Among all tested forms, only model M1 derived by Cieszewski (2004) simultaneously showed good statistical properties, sensible asymptotes, and adequate characteristics of the inflection point at the most productive sites. Therefore, M1 was selected for developing the site index curves for pedunculate oak in the Ravni Srem region. The visual inspection of residuals confirmed the suitability of the chosen model, showing no significant trends along any of the four observed predictors (Fig. 4). The identical GADA derivation has also been found to be the most suitable for describing observed height-age patterns for other economically important tree species (Nord-Larsen and Johannsen 2007; Bravo-Oviedo et al. 2008; Nord-Larsen et al. 2009).

The occurrences of systematically biased results are less likely with the well-designed forest inventories conducted across broader spatial scales. The main advantage of such databases is that they typically enclose a wide range of age classes and site conditions (Stimm et al. 2021). However, the suitability of inventory datasets for growth modeling of natural stands highly depends on the reliability of age records (Sharma et al. 2011). When reliable age records are available at plot or stand level, inventory surveys evolve into the comprehensive data source for height growth modeling (Huuskonen and Miina 2007; Brandl et al. 2018). Such circumstances significantly reduced the concern that some sites are better represented than others. The familiar fact that light-demanding species are naturally preconditioned to even-aged structures, together with an enormous commercial interest in pedunculate oak wood, has affected the management plans in Ravni Srem to contain accurate records of the stand's age. Yet, tracking the height changes over time is a far more challenging task when the forest management plans are established based on a stand-wise network of temporary sample plots (TSP). To overcome the lack of long-term observation, Socha et al. (2019) proposed a pragmatic approach that largely relies on the assumption of theoretically expected height growth patterns along an ecological gradient. The given cost-effective method assumes the calibration of dynamic models based on artificially established site-specific height growth trajectories. We interpret that here-in presented results support the findings of Socha et al. (2019), where under favorable circumstances, data from TSPs can be used as another valuable source of data for creating site index curves.

The applied procedure of reference age selection showed that a wider range of ages could be safely used for dominant height prediction. Among several available criteria for the final choice, particularly practical seems the opinion of Diéguez-Aranda et al. (2005) that the reference age should be much earlier than the final harvesting age to support management decisions about silvicultural treatment. Having in mind that the common rotation period for pedunculate oak in Serbia is 160 years (Rađević et al. 2020), we believe that an age of 100 years could be conveniently used for site index referencing. The obtained results follow the reference age used for oak in Poland (Socha et al. 2020) and various other tree species with prolonger rotation periods (Seki and Sakici 2017; Matisons et al. 2018; Sharma et al. 2019). Contrary to that, the reference age for oaks in Spain and Denmark is 60 and 50 years, but the oldest stands rarely exceed 25 m in height (Barrio-Anta and Diéguez-Aranda 2005; Nord-Larsen et al. 2009).

Since we just developed the first site index curves for pedunculate oak in Serbia, a comparison was only possible with models from other European countries. The GADA equations for oak were calibrated based on the data collected across the entire territory of Denmark (Nord-Larsen et al. 2009) and Poland (Socha et al. 2020). However, direct comparison with these results will be omitted since such considerations could be highly biased because the authors did not specify which oak species they modeled. In Spain, the regional-wide GADA models depicted the height growth of pedunculate oak in low-quality stands situated on the southwest margins of the species distribution (Barrio-Anta and Diéguez-Aranda 2005). The expected heights in the north-west part of Spain at the reference age of 100 years on the best sites reach 28 m which is 10 m lower than in Serbia. For comparison with pedunculate oak growing in less favorable site conditions than floodplains (Bauhus et al. 2017), we have used the figures from the current yield tables of Baden-Württemberg, Germany (1993). Anticipated heights over there range between 20 and 34 m in 100 years, meaning relatively higher values are recorded in Serbia at both margins of the oscillation spectrum. Accordingly, it seems that sites with underground water can provide better-growing conditions than those where rainfalls represent the principal water source. Even though it is reported that pedunculate oaks are less drought-sensitive than other native tree species in Central Europe (Mette et al. 2013; Pretzsch et al. 2013), higher water balance apparently has positive effects on height growth, what is also found by Stimm et al. (2021).

Based on unclear criteria, Špiranec (1975) divided the observed height spectrum of pedunculate oak in Croatia into three mean growth patterns, which seem incompatible with presented site index curves and underlying empirical data (Fig. 6). These historical curves are static and include binomial relative coefficients, so performance comparisons are limited only to visual exploration. Given overlapping of models demonstrated a quite inferior flexibility performance of Špiranec’s growth model. Besides, the old curves at the best sites predict maximal reachable heights far above all empirical measurements in the Ravni Srem region of Serbia. From the perspective of ecological conditions, we are dealing with one forest complex in two countries, and therefore given high predictions do not reflect reality. Due to the comparable growing conditions and the comprehensiveness of the data used for model calibration, the developed height growth models—here in particular M1—can potentially be used in a wider application area beyond our study region. We see great potential, especially for the Spačva forest systems in Croatia, but also for other pedunculate oak forests growing in the Sava River basin.

Mapping the present site index of pedunculate oak stands provides a reliable impression of the possible productivity at the studied sites. Although in coarse resolution, this kind of map is a practical way to communicate the results and represents a necessary reference point for the assessment of expected changes in dominant heights growth under the influence of climate changes (Albert and Schmidt 2010; Weiskittel and Kershaw 2010; Castaño-Santamaría et al. 2019) and oscillation in the CO2 sequestration (Kim et al. 2017). Such investigations are essential to adapt forest ecosystems to the effects of climate change and environmental changes, while shown results are the first step in perceiving the related risk (Brockerhoff et al. 2017). Contrary to several other European studies, the significant age-related trend of site indices variation was not detected in this research (Albert and Schmidt 2010; Sharma et al. 2011; Socha and Tymińska-Czabańska 2019; Stimm et al. 2022). That could be somehow expected since TSP data already infiltrated the potential age-related trends into the height growth model. It is important to highlight that such independent evaluation of stand indices across age classes is another significant advantage of the used modeling approach and represents a highly desirable property for mapping height growth potential across broader areas. On the other side, the appropriateness on a wider scale is probably a trade-off for reducing the performances on the local level, which means that the actual growth of a particular stand may deviate from the model-predicted trajectory. The maps of expected heights for the entire investigated region prove the existence of a southwestern gradient in site indices variability, where predicted SI₁₀₀ is considerably lower in the eastern part of the studied area called lower Srem. Altogether, observed trends and the first site productivity indicators should significantly improve current silviculture practices and facilitate decision-making in forest management. The presented results could complement the existing typological classification of forest types and, in perspective, offer a new stand delineation based on site productivity classification.

Derived site indices are the transitional solution which should be perceived as a reliable tool for site classification. In further research, the appropriateness of presented site index models needs to be more closely verified and compared with models fitted to longitudinal data. To comprehensively observe the growth dynamics of the individual trees and stands from the first years to maturity, it is necessary to conduct stem analysis and set up permanent sample plots. Furthermore, the appearance of regional growth differences should be considered as a result of the variation in local growing conditions, so their inclusion in models will improve the understanding and prediction of pedunculate oak growth in this area. Using the models with built-in environmental predictors, the effects of site conditions change on the growth dynamics could be evaluated and eventually compared with here-in-produced maps. Depending on the type and complexity of the employed model, the derived site index curves represent a vital component of classical yield tables or modern growth simulators.

5 Conclusion

In this research, we presented the first height growth models and the very first implementation of the GADA methodology for pedunculate oak in Serbia. Moreover, the established site index seems to be the second found application of the dynamic model for pedunculate oak in Europe, and the first one that depicts the maximal heights this tree species can reach. Site indices are calibrated based on an artificially knitted growth series established using the management inventory data. Underlying data from TSPs is more desirable than other commonly used sources since it is the most cost-effective data source that simultaneously represents the broad range of growth conditions and age classes across wider regions. From five tested GADA models, the best compromise between statistical and biological considerations is demonstrated by a polymorphic equation with multiple asymptotes from Cieszewski (2004). The minimal relative error (RE%) for the developed base-age invariant site index model is registered for a reference age of 100 years. Compared to Špiranec’s (1975) height growth models, the presented site index curves demonstrated greater flexibility and sound predictions over the scope of empirical data. Finally, we have provided spatially continuous maps of the expected site index in a reference age for studied stands with the recorded presence of pedunculate oak in the Ravni Srem Region. Due to the utilized modeling approach, estimations of SI₁₀₀ are independent of age classes, representing an advantageous property for mapping at a broader scale.

Obtained results will impact current management practice as the first productivity figures for pedunculate oak in Serbia and serve as a reference point for a better understanding of future climate-induced shifts in height growth. In further research, the calibration of the site index should be checked with longitudinal data and expanded with environmental variables to enable climate-sensitive predictions.