Overview of methods and tools for evaluating future woody biomass availability in European countries

The lifespan of many tree species is over several hundred years; therefore, knowledge regarding tree and forest stand dynamics, and long-term impacts due to environmental change is largely incomplete. Retrospective tree ring analyses can only partly close this knowledge gap, as this method indeed offers insights into tree growth, but not into stand dynamics. Forest inventories primarily assess managed forests, where the influences of climate change and thinning might be coalesced, and difficult to differentiate. Models are frequently used as a means to circumvent data collection and subsequent analyses1. However, modelling is no substitute for underlying field data, and the full potential of any modelling approach is only fulfilled when feedback between modelling studies and empirical analyses are achieved. A unique source of information, however, was provided by long-term experimental plots established in approximately 1872, the founding year of the International Union of Forest Research Organisations2. These plots, which were surveyed 10–20 times until the present day, provide the longest existing time series data on forest stand dynamics available, approximately 140 years. Originally, the study stands were established to examine stand growth principles3, but not for growth trend analysis2. In addition to the per se uniqueness of the sites, the locations are in Central European regions where the longest time series on driving variables (precipitation and temperature) dates back to 1781. The original data acquisition objective was to support sustainable forestry at a local scale; however, we subsequently used these unique records to quantify and characterize changes in Central European forest growth. We chose Norway spruce (Picea abies (L.) Karst.) and European beech (Fagus sylvatica L.) as the study species. These taxa dominate Central Europe’s forests occupying 30%, that is, a total area of 14 × 106 ha of all forest areas. The plots selected for this study represent pure, even-aged stands, which were established by planting or seeding. Site conditions varied broadly, and soils ranged from dry silty sands to moist deep silts. Since the first site observations and records in 1870, the plots were maintained under continuous scientific control, and surveyed on a single tree basis. Therefore, investigators excluded plots and reports impacted by disturbances, including storms or bark beetle infestations. We included only unmanaged, or at most moderately thinned, but always fully stocked plots. This selection resulted in a unique survey data set from 36 spruce and 22 beech plots. Based on these data we show that both species currently exhibit significantly faster tree growth, stand volume growth and standing stock accumulation than still in 1960 and the decades before. Self-thinning lines remain constant, while growth rates increase indicating the stock of resources have not changed, while growth velocity and turnover have altered. This means stands still follow similar general allometric rules, but proceed more rapidly through usual trajectories. As we can demonstrate, this results in stands currently having lower tree numbers per unit area than past stands at the same age. Our data also reveal that the growth acceleration is stronger on fertile sites, which is supported by scenario runs with an ecophysiological growth model. Results Changes in stand dynamics and environmental conditions First, we pooled our data and compared them to standard yield tables4 5 (Fig. 1). Yield tables, common forestry tools that tabulate stand growth age dependently, were developed primarily from 1795 to 1965. Yield table data were derived from long-term plot field survey data; and thus they served to represent past growth conditions in this comparison. We found stand growth rates and standing stocks after 1960 (empty symbols in Fig. 1) exceeded the yield table ranges by 50–100%, which called the validity of yield table range data into question, and suggested essential changes in stand dynamics. For the same period our plot surveys spanned, we compiled available data on environmental variables reported to drive forest growth dynamics (Fig. 2). For Central Europe, forest environmental and growing conditions exhibited significant changes since the first experimental forest plots were established in 1870 (Fig. 2). During the addressed period, the atmospheric CO2 concentration rose from 295 p.p.m. in 1900 to approximately 390 p.p.m. in 2010 (refs 6, 7; see Matyssek and Sandermann8 for the possible effects of atmospheric composition on trees). This means an increase of more than 30% within nearly one century. Wet N-deposition increased by 0.5–1.0 kg ha−1 per decade7 9. Throughout Central Europe, average total N-deposition increased from approximately 2.5 kg ha−1 per year in 1900 to more than 9 kg ha−1 per year in the first decade of the twenty-first century6. Global average temperature has increased by roughly 0.7 °C (ref. 7) within the twentieth century. During the same period, the average temperature in Europe has risen by 0.95 °C (ref. 7). In Germany, the mean annual air temperature increased by 1.0 °C (ref. 10) during the twentieth century; and the sum annual precipitation increased by 9% during the same period. However, the annual distribution varied. During the winter months, precipitation increased by 19% during the last century, and rainfall in summer decreased by 3% on average. If only the second half of the twentieth century is measured, summer precipitation shows a 16% reduction. In addition, the rise in atmospheric CO2 concentration, the higher N-deposition and the increase in air temperature were two to three times higher in the second half of the twentieth century compared with the first half. However, the strong temperature increase during the last 50 years reported for all of Europe7, but was not reported for Germany. The annual mean temperature increase for 1950–2000 was equal to the century average at 1.0 °C; only winter temperatures showed a higher increase in the second half of the twentieth century10. Higher temperatures will also extend the growing season11 12. Menzel and Fabian11 reported the average annual growing season has been extended by 10.8 days, since the early 1960s. Chmielewski and Rötzer12 also found the vegetation period was lengthened between 0.6 and 6.3 days per decade in different European natural regions during 1969–1998. Based on temperature data from the four climate stations used in this study (Supplementary Table 1), the length of the growing season, defined as the number of days annually with temperatures above 10 °C, was calculated for the last 110 years. Averaged over the climate stations, the growing season was extended by 22 days. The main increase, however, was detected over the last 50 years (Fig. 2). This suggests notable wood volume growth rate increases at the stand level over the last 100 years (Fig. 1), coinciding with an increase in resource supply (CO2, N), together with an extended growing season accompanied by changes in other climatic variables (Fig. 2). These observations justified statistical analyses of growth trends, and model-based examination of the underlying mechanisms. Growth trends of key stand variables First, we employed linear mixed models (LMMs) to determine whether the most important stand characteristics were dependent on only stand age, or also on calendar year. Standing wood volume (V), mean diameter (dq), dominant height (ho=mean height of the 100 tallest trees per hectare) and mean tree volume ( ) currently grow significantly faster than in the past (Figs 3, 4 and Supplementary Tables 2,3). Under the environmental conditions of the year 2000, any given mean diameter was attained following stand establishment up to more than one decade earlier than its counterpart in 1960 or before (Fig. 3). Stand volume growth, expressed as periodic annual increment of volume (PAIV) changed from 1960 to 2000 by, respectively, 10% and 30% for Norway spruce and European beech (Table 1). Most stands continued to accumulate volume, and have not reached a final constant yield plateau (Figs 1 and 3). For a 60-year-old Norway spruce stand, we expected a maximum standing volume V of 760 m3 ha−1 in 1960, and 810 m3 ha−1 in 2000. In a 130-year-old European beech stand, the expected maximum volume in 1960 was 630 m3 ha−1 compared with 700 m3 ha−1 in 2000 (Fig. 3). A consequence of accelerated stand development was a more rapid decrease in stand tree number N per unit area (Fig. 3), and change in tree mortality rate, MORT (Fig. 4). A comparison between 1960 and 2000 showed a 17% decrease in Norway spruce tree number and a 21% decrease in European beech. We did not detect a significant change in Norway spruce mortality rate, however, European beech exhibited a −17% change (Table 1). The calendar year effect on size and stand growth, and volume accumulation were significantly positive, and the calendar year effects on tree numbers were significantly negative. Significance levels obtained with LMM were at least P 5,000) to assure only fully stocked stands were included. All stands are monospecific and even-aged; they originated from planting or seeding. We used such long-term observational plots to determine whether growth rate changed over time, and how changes affected tree growth, standing volume, stand density, mortality and other stand attributes. The selected plots represent growth conditions in the plains and highlands of Middle and Central Germany (Supplementary Fig. 1). The localities occur from 330 to 843 m above sea level. Although European beech plots dominate the Atlantic plains and highlands climate, Norway spruce plots are located in submontane and montane highlands, and the pre-alpine mountain zone (Supplementary Tables 6–8). The long-term mean temperature and annual precipitation exhibits a broad range across both species (5.7–8.5 °C and 605–1,369 mm per year), as well as for each species separately (for example, for Norway spruce, 6.9–8.3 °C and 812–1,256 mm per year). The plot distribution over 13 eco-regions and 11 geological zones is reflected by the broad spectrum of soil types. The poorest soils are podsols derived from sandstone and cretaceous material in the Palatinate and Upper Palatinate regions; the most fertile soils are parabrown soils from diluvial loess-loam in the pre-alpine highlands. The majority of stands are supported on soils of mediocre fertility on Triassic, Jurassic and Cretaceous formations between Frankfurt and Munich (Supplementary Fig. 1, see Supplementary Tables 6 and 8 for further details). The data set comprises plots in present day mature stands surveyed up to 18 times since 1870, but also in young stands established in the last decade and only surveyed twice. Hence, the plots cover both historic and recent growth behaviour under respective environmental conditions. The broad variation of stand age (21–188 years), dominant height (10.7–44.4 m), tree number per hectare (133–11,238 trees ha−1) and quadratic mean tree diameter (5.4–54.4 cm) show the plots represent a rather wide range of stand developmental stages. Total yield (50−2,459 m3 ha−1), standing volume (50–1,637 m3 ha−1) and periodic annual volume increment (7.1–41.5 m3 ha−1 per year), as well as SI (19.8–43.1 m) emphasize the wide spectrum of site conditions and productivity levels (cf. Supplementary Table 9). Survey and evaluation of long-term observational plots We based our analyses on the International Union of Forest Research Organisations2 standard variables, which quantified above ground mean tree and stand stem volume, rather than single tree volume or biomass. Therefore, additional assumptions for scaling from volume to mass were avoided; however, it required the following variable definitions: (i) all stand variables relate to a unit area of 1 ha (104 m2); (ii) tree diameter, dq (cm), refers to the quadratic mean diameter at breast height (1.30 m) for all trees per plot; (iii) dominant height, ho (m) is the mean height of the 100 tallest trees per hectare; (iv) mean tree volume, (m3) is the arithmetic mean stem volume; (v) annual tree volume growth, (m3 yr−1) is the mean annual volume growth of the mean trees with volume ; (vi) PAIV (m3 ha−1 per year) refers to the entire stand’s mean annual stem volume growth during a period between two surveys; (vii) standing stand volume, V (m3 ha−1) is the accumulated stem volume per hectare; (viii) the successive surveys of remaining, dead and harvested trees generates the current tree number per unit area, N (ha−1), and enables the calculation of the annual tree mortality rate, MORT (% per year); and (ix) the plots’ SI is expressed as measured or expected stand height at age 100 years. We used the prevalent yield tables by Assmann and Franz4, and Schober5 for Norway spruce and European beech, respectively. For complementary details, see Supplementary Methods. Dependency of stand variables on age and calendar year We examined whether stand development on observational plots reflects any long-term growth trends by modelling stand characteristics dependent on stand age and calendar year. Certainly, stand characteristics from successive surveys (for example, PAIV, standing volume (V) and tree mortality rate (MORT) depend on age. An additional calendar year effect on stand characteristics would indicate a growth trend; if stands at a defined age perform differently in different calendar periods or decades, this will indicate a change in growth and site conditions. For investigating this, the following basic LMM structure was the most appropriate: Variable Y represents the stand characteristic of interest (for example, PAIV, V and so on), untransformed or logarithmized, depending on whether the logarithmic transformation rendered a better model fit. Similarly, A represents stand age, untransformed or its logarithm. Choosing appropriate combinations of Y and A logarithmic and untransformed values allowed us to sufficiently cover nonlinear age-dependent relationships with a linear regression model. The second explanatory variable, calendar year, corresponding to a given observation, is indicated by the variable year. The indices i, j and t represent the location an observational plot is included, the plot itself and the point of time a plot survey has occurred. Fixed effects parameters are β 0–β 3, whereas b i and b ij are location and plot random effects (b i ~N(0,τ 1 2 ), b ij ~N(0,τ 2 2 )). Including these random effects, we avoid biased results due to the plot-specific and possibly also location-specific autocorrelation among the observations. Finally, ε ijt denotes i.i.d. errors (ε ijt ~N(0,σ 2)). The calendar year effect and its interaction with age (represented by β 2 and β 3 parameters) were only maintained in the model when they were statistically significant. Otherwise, the model was reduced accordingly and fitted again. If the interaction was significant, but not the isolated year effect, both were maintained in the model44. High stand ages in our data primarily occurred with recent calendar years only, therefore, we excluded specific observations beyond a certain age to develop a balance of age-calendar year combination data set. Our models were fitted for 60 years and younger stand ages in Norway spruce, and 130 years and younger in European beech. For most stand characteristics as response variables this resulted in a sample size of n=157 and n=225 for Norway spruce and European beech, respectively. For the growth variables PAIV and the sample size reduced to n=141 (spruce) and n=217 (beech) as there is no increment information available for the plots’ last surveys. The mortality rate, MORT, could be meaningfully analysed for the completely unthinned plots only which results in n=90 (spruce) and n=119 (beech). All models were fitted by maximizing the restricted maximum likelihood criterion (cf. Zuur et al.44). Allometric relationships of stand growth and size variables The relationships between mean tree growth and mean tree size ( versus ), and tree number per unit area and mean size (N versus ) are cornerstones of allometric theory45 46 47 48 49. In the double logarithmic scale, both relationships follow a straight line ln(y)=a+b × ln(x) (equivalent to y=e a × x b ) with rather general and species-overarching values for the slope b. However, it is widely accepted that line levels, represented by intercept a, depend on environmental conditions and species21 50 51. We used a LMM, which is very similar to the basic model shown above to test the extent both allometric relationships are influenced by calendar-year-dependent trends: where y and x represent and , or N and , respectively. The variable and index names are defined the same as in equations 1. We employed exactly the same data used in fitting the age trend models, including only stands younger than 61 (Norway spruce) and 131 (European beech) years. Site dependency of the relationship’s temporal trend Results of the previous analyses suggested an upward shift with time (significant parameter β 2 in equation 2) in the allometric relationship between and as the common mechanism underlying the observed growth trends. With the same data, we tested any change in allometry dependent on site conditions by formulating the following linear mixed regression model: where SI is the respective plot’s SI, expressed as an expected stand height at an age of 100 years (see above). The other variable meanings and names are defined exactly the same as above. If parameter β 3 differs significantly from zero, this indicates the allometric shift depends on site quality. All statistical analyses were performed with R 3.0.2 (ref. 52). Growth trend and changed arrival age at threshold values To quantify how stand characteristics changed, we choose stand age 75 years, calculated how the stands perform at that age in 2000, and divided that by stand performance at the same age in 1960. For this purpose, both values were derived for the respective stand variables from the fitted model equations using the fixed effects parameter estimates (equations 1, 2 and Supplementary Tables 2–4) while setting the random effects to zero. For PAIV, for example, this procedure yields PAIV age 75,2000 and PAIV age 75,1960, and the ratio RPAIV age 75,2000/1960=PAIV age 75,2000/PAIV age 75,1960, which reflects the growth trend since 1960. The age when a mean tree variable or a stand characteristic arrives at a defined threshold value is a practical and relevant measure. Assume the fixed effects parameters of equations 1 have been estimated, we have the following equation for estimating , which is the general expected value for a given stand characteristic or its logarithm: This can be rearranged as which allows us to estimate A, the mean age (or its logarithm) when a certain threshold value is reached under the environmental conditions of a given year. Process-based modelling The physiological growth model BALANCE53 54 we used for scenario analyses simulates the three-dimensional development of individual trees in a stand, and estimates the consequences of environmental influences. Tree development is calculated as a response to individual environmental conditions, and as environmental conditions change with individual tree development, the influences of competition, stand structure, species mixture and management options can be assessed with the model (Supplementary Fig. 2). Initial tree biomass is calculated from the dimensional variables tree height, height to crown base, diameter at breast height, tree position and crown radii. Biomass increase is the result of the interaction between physiological processes, which are dependent on the physical and chemical microenvironment. These are in turn influenced by stand spatial structure. Asymmetric crown shapes are included, and generate a spatially explicit representation of the environment. The calculation levels vary from stand level to individual trees, from tree components (that is, foliage, branches, stems, and fine and coarse roots) to crown and root layers, which are spatially subdivided into segments. Consequently, an increase in biomass is simulated based on the carbon and nitrogen uptake from each segment, depending on its energy supply and resource availability. By using weather data, microclimate and water balance are simulated for each layer and segment, respectively. Air temperature and radiation within the stand is calculated for every crown layer of every tree on the basis of leaf area distribution for the respective tree and its competitors. The spatial distribution of light and water availability is estimated on a daily basis. Water balance simulation examines soil conditions in different soil layers, where vertical and horizontal water flows between rooted and non-rooted fractions are considered. Based on the Penman-Monteith55 approach, potential evapotranspiration is estimated, from which the actual evapotranspiration of a tree is calculated using maximum water uptake derived from water content within soil volume pervaded by fine roots. Water can be exchanged between rooted and un-rooted soil layers. Total soil water content is reduced by drainage, which is equivalent to percolation from the deepest soil layers. Foliage biomass and leaf area as well as light availability and photosynthetically active radiation (PAR) absorption change with the onset of bud burst. A tree’s foliage emergence date determines its assimilation and respiration rate, but also alters the environmental conditions in the immediate surrounding area. Bud burst of a tree species is estimated using an air temperature sum model, whereas foliage senescence is simulated depending on the respiration sum for each segment of a tree. Based on the aggregated driving variables, all physiological processes, that is, assimilation, respiration, nutrient uptake, growth, senescence and allocation can be calculated for each individual tree. Nutrient uptake is the result of demand, supply and absorption capacity, whereby demand is based on the difference between the actual nitrogen concentration, and a given optimal concentration. Supply is defined by soil characteristics of the rooted volume, uptake capacity by root surface, and its specific absorption rate. Physiological processes are calculated in 10-day time steps using aggregated results of daily environmental conditions. Gross primary production is estimated depending on leaf surface, absorbed PAR, temperature, CO2 concentration, water and nitrogen supply (Supplementary Fig. 3). Total respiration is the sum of maintenance losses and growth respiration. Maintenance respiration is calculated for each segment as a function of biomass, specific respiration rate and temperature. Growth respiration is estimated as a constant fraction of maximum photosynthesis. The fixed carbon not required for respiration is distributed into plant compartments, including foliage, branches, stems and roots. The available carbon for allocation is apportioned into different compartments according to growth and respiration demands. Carbon allocation is defined by the relationships between the compartments according to the functional carbon balance theory56, and the pipe model theory57. Consequently, all tissues within a segment, that is, foliage and branches, or fine and coarse roots, as well as the amount of stem wood, are mechanistically linked to each other. Dimensional tree growth is estimated annually, based on biomass accumulation during that year. Volume expansion depends on the necessary amount of twigs and transport branches, and the amount of coarse roots within root segments. Therefore, crown development is preferred in the direction of best assimilation conditions during the previous year. If net assimilation rates are negative, the crown segment is regarded as dead. If no segments contain living biomass, the tree is assumed dead and removed from calculations. BALANCE has been extensively validated for basic micro-meteorological and physiological processes, for water balance, annual tree development and entire stand development54 58. Detailed descriptions of BALANCE, and single modules can be obtained from Grote and Pretzsch53, or Rötzer et al.54 Scenarios calculated with BALANCE Climate data from four German climate stations with daily time series for more than 100 years formed the foundation of the growth simulations using BALANCE. Supplementary Table 1 shows the geographical coordinates, mean air temperature and precipitation values for the chosen simulation periods. These four stations are representative of most climate regions in Central Europe. We chose a sandy loam soil type with an available field capacity of 186 mm to a maximum rooting depth of 1 m. Field capacity decreases from 35 to 25 vol% with increasing soil depth, whereas the wilting point was set constant at 8 vol%. Plant available N was assumed low at the beginning of the simulation. This was justified because historically, forests in Germany were displaced by agriculture to nutrient poor sites. Growth development of a 30-year-old Norway spruce and 35-year-old European beech stand (Supplementary Table 11) was simulated for time spans 1901–1930 and 1981–2001. The first simulation scenario (reference) was the 1901–1930 period using the daily climate record obtained from each station, and continuously increasing atmospheric CO2 concentrations from 295 p.p.m. to 307.p.p.m., and N-depositions from 6 kg N ha−1 per year to 7 kg N ha−1 per year. The second scenario reproduced the recent climate conditions from 1981 to 2010 while keeping N-depositions and CO2 concentration on the previous level. In scenario 3, an increase of atmospheric CO2 concentrations from 338 p.p.m. from 1981 to 389 p.p.m. in 2010, and increased N-deposition from 15 kg N ha−1 per year in 1981 to 20 kg N ha−1 per year in 2010 was included in the model. Scenario 3 thus integrates all recent environmental conditions. Atmospheric CO2 concentration and N-deposition data were derived from Churkina et al.6 In this way, the single and overall influences of climate, CO2 and N-deposition were analysed. Author contributions H.P. initiated the study, interpreted the data and wrote the paper. P.B. performed statistical analyses, interpreted the data and wrote the paper. G.S. compiled the data. E.U. interpreted the data and revised the manuscript. T.R. performed and interpreted simulation runs and wrote the paper. Additional information How to cite this article: Pretzsch, H. et al. Forest stand growth dynamics in Central Europe have accelerated since 1870. Nat. Commun. 5:4967 doi: 10.1038/ncomms5967 (2014). Supplementary Material Supplementary Information Supplementary Figures 1-3, Supplementary Tables 1-11, Supplementary Notes 1, Supplementary Methods and Supplementary References. Show more