Interactive comment on “ Validation of modelled forest biomass in Germany using BETHY / DLR ”

Abstract. We present a new approach to the validation of modelled forest Net Primary Productivity (NPP), using empirical data on the mean annual increment, or MAI, in above-ground forest stock. The soil-vegetation-atmosphere-transfer model BETHY/DLR is used, with a particular focus on a detailed parameterization of photosynthesis, to estimate the NPP of forest areas in Germany, driven by remote sensing data from VEGETATION, meteorological data from the European Centre for Medium-Range Weather Forecasts (ECMWF), and additional tree coverage information from the MODIS Vegetation Continuous Field (VCF). The output of BETHY/DLR, Gross Primary Productivity (GPP), is converted to NPP by subtracting the cumulative plant maintenance and growth respiration, and then validated against MAI data that was calculated from German forestry inventories. Validation is conducted for 2000 and 2001 by converting modelled NPP to stem volume at a regional level. Our analysis shows that the presented method fills an important gap in methods for validating modelled NPP against empirically derived data. In addition, we examine theoretical energy potentials calculated from the modelled and validated NPP, assuming sustainable forest management and using species-specific tree heating values. Such estimated forest biomass energy potentials play an important role in the sustainable energy debate.


Introduction
Models of carbon uptake by plants play an important role in answering questions concerning the mechanisms driving the carbon cycle and the roles of terrestrial carbon sinks and sources (Cox et al., 1999).Carbon uptake by plants, measured as Gross Primary Productivity (GPP), can be predicted by simple models that describe the physical, chemical, and plant physiological processes of plant development, as well as the interactions between plants and the atmosphere.Such "deterministic" models (sometimes also called "mechanistic" or "Monteith-type" models) calculate photosynthesis following the methods of Monsi and Saeki (1953) and Monteith (1965).
The idea behind these Monteith-type models is that the carbon uptake of sufficiently watered and fertilized plants is linearly correlated with the energy of the incident photosynthetically active radiation (PAR), or more precisely, the fraction of the PAR that is actually absorbed by the plants (fPAR).Following this approach, it is possible to calculate GPP for each vegetation type from the absorbed solar radiation (fPAR) and the light use efficiency (LUE) of the vegetation type.The LUE can be affected by environmental stress factors, particularly temperature, water limitation, and nitrogen availability.Species-specific fPAR values may be estimated by measurement of dry biomass accumulation, or may be derived from satellite data.
GPP, as estimated by such a model, can be converted to NPP by considering temperature-dependent maintenance respiration.Maintenance respiration can be estimated using allometric functions regarding leaf and root distribution following the approach of Ryan et al. (1995), or using the Leaf Area Index (LAI) of the vegetation following Running et al. (2000).In either case, NPP is defined as the remainder after plant maintenance respiration is subtracted from GPP.In a further step, Net Ecosystem Productivity (NEP) can be calculated by subtracting the heterotrophic respiration in an ecosystem from the ecosystem's NPP.
The Monteith-type model architecture has been used many times.For example, the C-Fix model, a Monteith-type parametric model by Veroustraete et al. (1994), was used by Verstraeten et al. (2006) to estimate net ecosystem fluxes for all of Europe.C-Fix is driven by vegetation type data of the Normalized Differenced Vegetation Index (NDVI) and meteorological data (temperature and daily incoming global radiation) obtained from about 800 weather stations administered by the World Meteorological Organization (WMO).Verstraeten et al. (2006) validated their results with eddy covariance flux tower measurements, obtaining an R 2 of 0.84 for pine forests and 0.59 for mixed deciduous forests.The Carnegie-Ames-Stanford Approach (CASA) model introduced by Potter et al. (1993) and validated by Potter et al. (2001Potter et al. ( , 2003) ) is another example of a Monteith-type model.The CASA model is driven by monthly NDVI data from the FASIR database of the Goddard Space Flight Center, monthly temperature and precipitation data from the International Institute for Applied Systems Analysis (IIASA), and monthly PAR data from the Goddard Institute for Space Studies.Validation of CASA was performed against atmospheric CO 2 concentration data from NOAA and the Geophysical Monitoring from Climate Change Flask Sampling Network, and obtained R 2 values between 0.09 and 0.67.
When the LUE approach is integrated into a coupled soilplant-atmosphere model, such as the Atmosphere-Land Exchange (ALEX) model, daily estimates of evapotranspiration and carbon assimilation fluxes can also be obtained (Anderson et al., 1997).Recently, more sophisticated models have begun to be developed that take this integrative approach even further.In computing the uptake of carbon by plants, these so-called "dynamic" models take into account the complex interactions between plants, soil, and the atmosphere, but also account for the carbon released by both plants and soil in a manner that respects the conservation of energy and momentum.At present, only a few dynamic models have been published.Examples are the Lund-Potsdam-Jena (LPJ) model developed by Prentice et al. (1992) and modified by Bondeau et al. (2007), the Equilibrium Terrestrial Biosphere Model (BIOME3) by Haxeltine and Prentice (1996), and the ORganizing Carbon and Hydrology in Dynamic EcosystEms (ORCHIDEE) model by Krinner et al. (2003).These global models are driven by meteorological input data, and phenology is calculated internally from those inputs using per-species physiological parameters.The spatial resolution for dynamic models can range from degrees, for global models such as Prentice et al. (1992) and Haxeltine and Prentice (1996), to kilometres, for regional models such as Wisskirchen (2005).Model outputs are typically GPP, NPP and NEP, total ecosystem respiration, and evapotranspiration.
This study used the Biosphere Energy Transfer Hydrology (BETHY/DLR) model, a dynamic model based on the Jena Scheme of Atmosphere Biosphere Coupling in Hamburg (JSBACH) by Knorr (1997), which was designed for global applications (see also Knorr and Heimann, 2001).It was modified by Wisskirchen (2005) for application to regional modelling.
Model validation is often conducted using data from devices called eddy covariance flux towers.The relationship between carbon and energy flux has been studied in international networks such as FLUXNET (Baldocchi et al., 2001) and AmeriFlux, as well as in projects such as EUROFLUX (Valentini, 2003) and CarboEurope.This research has shown that eddy covariance flux tower measurements can be used to quantify NEP at the spatial scale of the footprint of a tower (Baldocchi, 1997).As mentioned above, NEP may also be calculated by subtracting heterotrophic respiration from NPP.Therefore, robust methods have been developed to estimate heterotrophic respiration in order to convert NEP, as measured by eddy towers, into NPP (or, by considering plant maintenance respiration as well, GPP).
For example, the MODIS GPP product (MOD17, C4.5) was validated with eddy tower CO 2 flux estimates across diverse land cover types and climates (Heinsch et al., 2006).The main test areas were forest ecosystems in North America, but chaparral, oak savannah, northern grassland and Arctic tundra were also included in the investigation.It was found that MODIS overestimated GPP by about 20 % to 30 %, but this depended strongly on season and ecosystem type.Comparison of annual MODIS GPP (modelled with global meteorological data from NASA's Global Modeling and Assimilation Office) to tower-based GPP measurements yielded an R 2 of 0.72.
The primary objective of this study is to present a new approach to the validation of modelled NPP.We compare output from BETHY/DLR, run at 1 km 2 spatial resolution, with empirical measurements of the mean annual increment (MAI) in above-ground biomass (including bark) observed in forests in Germany.The MAI data are available at a regional scale called the NUTS-1 level; NUTS is an abbreviation for "Nomenclature des Unités Territoriales Statistiques," and is a system of hierarchically organised territorial units used for statistical purposes.The NUTS-1 MAI data were obtained from the National Forest Inventory (NFI) of Germany.
A secondary objective is to use our modelled and validated NPP to estimate theoretical energy potentials, given sustainable forestry practices, for the area of study.Sustainable energy potentials from forests are a key element in planning a sustainable energy economy for Germany (and, of course, the rest of the world), and so developing methods for estimating, and ultimately forecasting, these potentials is of great importance (BMVBS, 2010).

Model description
BETHY/DLR is a soil-vegetation-atmosphere-transfer (SVAT) model.SVAT models track the plant-mediated transformation of atmospheric carbon dioxide into energystoring hydrocarbons such as sugars, a process known as carbon fixation.BETHY/DLR models photosynthesis, and takes into account environmental conditions that affect it.Photosynthesis is parameterized following the combined approach of Farquhar et al. (1980) and Collatz et al. (1992), and treats separately the "light" and "dark" reactions of photosynthesis at the leaf level.A benefit of this design is that the photosynthetic rate can be mechanistically limited either by light availability or by the abundance of the carboxylation enzyme Rubisco, the key player in the Calvin cycle that fixes carbon.In addition, so-called C3 and C4 plants are distinguished in BETHY/DLR because significant differences exist between their carbon fixation physiologies; in particular, C4 plants such as sugar beet and corn can fix more atmospheric CO 2 at higher temperatures than can C3 plants such as barley and wheat.
In a second step, the rate of photosynthesis is extrapolated from leaf to canopy level taking into account both canopy structure and the interactions between soil, atmosphere, and vegetation.For closed canopies (trees, shrubs, and crops) the photosynthetic rate depends on the Leaf Area Index (LAI), a per-species metric of leaf upper surface area per unit ground area.To model self-shading, photosynthetic rate is reduced exponentially from the canopy top to the soil.In this approach, radiation absorption in the canopy is approximated in BETHY/DLR using the so-called "two-flux scheme" of Sellers (1985) with three canopy layers.
Besides photosynthesis, other energy transfers also need to be tracked.BETHY/DLR's energy balance model takes into account heat fluxes between the vegetation and the atmosphere above it, as well as the cooling effect of evapotranspiration from the soil and vegetation.Soil heat flux is also estimated, and the storage of heat in the canopy and in the air layer above the canopy is considered.
The coupling of these processes is of great importance, since temperature-dependent photosynthesis transforms light energy into chemical energy, and finally into carbohydrates, using water and CO 2 .Water is available for the plants from the soil, while evapotranspiration by plants and soil determine the water vapour deficit in the atmosphere, which is closely linked to the stomatal behaviour of leaves.Thus when considering the dynamic interaction of, for instance, the soil water balance and photosynthesis, the natural behaviour of vegetation can be reflected; this is the motivating idea of the SVAT approach.
The water cycle must also be modelled and included in this interaction scheme.In BETHY/DLR three water reservoirs are considered: soil water, snow, and "skin" or "intercepted" water on leaves and other parts of the vegetation.These reservoirs change in time and space depending on precipitation, temperature and evapotranspiration.A "bucket model" is used for calculating soil water dynamics, using the plant rooting depth as the depth of the soil core.In the bucket model, water availability for plants is governed by the soil water content above the permanent wilting point (at which water is so tightly bound by soil particles that it is unavailable to plants) and below the field capacity (at which the soil is full and further water added by precipitation or snowmelt runs off).The distribution of water within the soil is not modelled in BETHY/DLR, although this would be a worthwhile addition.Water limitation is considered by calculating the demand for evapotranspiration using the approach of Monteith (1965) with the criteria of Federer (1979), which state that evapotranspiration cannot be greater than the limit set by the soil water supply and the water uptake physiology of a plant's roots.
Autotrophic respiration is modelled in BETHY/DLR as the sum of maintenance respiration and growth respiration.Maintenance respiration is determined by plant-specific dark respiration rates, while growth respiration is proportional to the difference between GPP and maintenance respiration.For estimating NEP, heterotrophic soil microbe respiration is calculated as a function of temperature, scaled by the annual NPP and the effect of soil moisture (neglecting, for modelling purposes, other heterotrophic respiration).
The output of BETHY/DLR is a time series of NPP in daily steps, at the resolution and projection of the land cover classification.Here 1 km 2 resolution is used, in a latitudelongitude projection using the WGS84 (World Geodetic System 1984) datum.An overview of the input data used in this study and the internal model processes acting upon them is presented in Fig. 1.

Input data
The inputs of the BETHY/DLR model are three remote sensing datasets derived from SPOT-VEGETATION and MODIS, meteorological data provided by ECMWF, and two further datasets describing soil type and land elevation.

Remote sensing data
Per-pixel Leaf Area Index (LAI) time series were used to determine vegetation phenology.These time series were based on CYCLOPES (Carbon cYcle and Change in Land Observational Products from an Ensemble of Satellites) tenday composite datasets, which can be downloaded from the POSTEL (Pole d'Observation des Surfaces continentales par TELedetection) database.Criteria for the identification of gaps and outliers in the CYCLOPES datasets were defined and crosschecked against a hand-validated dataset.For each identified gap or outlier, harmonic analysis (HA), a type of superpositioning transformation, was applied to estimate a corrected value.Following this procedure, a global mean error of 9 % was found across a longer time series.Despite these gaps and outliers, the CYCLOPES dataset is the most consistent LAI dataset available (Garrigues et al., 2008); it was also chosen because it is available in the needed form of ten-day composites.
The CYCLOPES dataset also provides land cover and land use information for the year 2000, available as the product Global Land Cover 2000 (GLC2000).The Land Cover Classification System of the Food and Agriculture Organization of the United Nations (FAO), using 22 different land cover classes, was used to derive GLC2000 (Bartholome et al., 2002;DiGregorio and Jansen, 2001).The CYCLOPES LAI and GLC2000 datasets are both available as 10 • × 10 • maps in rectangular projection, with latitude and longitude using the WGS84 datum.Complete coverage of the study area (Germany) is available.
In order to make the GLC2000 land use/land cover classification useable for NPP modelling with BETHY/DLR, the GLC2000 vegetation classes were translated to BETHY/DLR's forest-related vegetation types.Each vegetation type in BETHY/DLR is described by biochemical parameters such as the maximum carboxylation rate and the maximum electron transport rate, encapsulating the maximum rates of, respectively, the light and dark reactions of photosynthesis.The internal parameterisation of BETHY/DLR allows a given GLC2000 vegetation class to be represented as a fraction of two BETHY/DLR vegetation types.In the context of this study, only four GLC2000 classes represented forest types found within Germany; the BETHY/DLR representation of these four classes is shown in Table 1.
To obtain information about fractional land cover, an additional dataset was used.The MODIS Vegetation Continuous Fields (VCF) dataset (DeFries et al., 2000;Hansen et al., 2002Hansen et al., , 2003) ) contains annual global data on percent tree cover at a spatial resolution of 500 m square, and is available for 2000 to 2005.For use in BETHY/DLR, the highresolution VCF data were aggregated to match the spatial resolution of the CYCLOPES data using a Geographic Information System (GIS).For GLC2000 class 6 ("Tree cover, mixed leaf type"), a weighting factor of 0.5 between deciduous and coniferous forest types was used (without scaling of those weights by the VCF fractional land cover).This approximation was assumed to be adequate, since only 15 % of German forest areas are described as mixed forest in the GLC2000, and since VCF contains no further information about the proportion of deciduous to coniferous trees.A similar approach was described by Jung (2008) using the AVHRR Continuous Field of Tree Cover dataset.

Meteorological data
In addition to remote sensing data, BETHY/DLR needs meteorological input data (Table 2).The ECMWF (European Centre for Medium-Range Weather Forecasts) provides this data at a spatial resolution of 0.25 From these data we calculated daily mean, minimum and maximum temperatures, daily mean cloud cover at three heights, and relative humidity.Daily temperature values were reconstructed at the 1 km 2 resolution of the land cover data, adjusting for the elevation difference between the ECMWF data and the elevation of each modelled pixel using the temperature gradient of the international standard atmosphere (−0.65 K per 100 m).
The daily average photosynthetically active radiation (PAR) was calculated from global radiation following the method of Burridge and Gadd (1974).This method calculates PAR from the incident sunlight for the given day and year, modified by atmospheric transmission, which depends on the degree of cloudiness.Daily average cloud cover was calculated using a weighted sum of each cloud layer.This approach leads to more exact results than the direct use of radiation forecast data, and is thus preferable to the direct use of ECMWF radiation data (Wisskirchen, 2005).Global radiation was calculated for each location at each one-hour time step.
Time series of the volumetric soil water content were needed to calculate the soil water budget of the model.Information about the soil type was taken from the International Soil Reference and Information Centre-World Inventory of Soil Emission Potentials (ISRIC-WISE) dataset, which has a resolution of 5 × 5 arcminutes.It is a harmonization of the global FAO-UNESCO Soil Map of the World (FAO-UNESCO 1974).

Eddy crosscheck
Before validating BETHY/DLR's modelled NPP for NUTS-1 regions across all of Germany, we performed a crosscheck of BETHY/DLR GPP results with eddy covariance flux tower measurements provided by FLUXNET.Two tower sites were selected, one in the Hainich forest and one in the Tharandt forest.The Hainich tower is to the west of Jena, Germany, in a mixed deciduous beech forest, while the Tharandt tower is south of Dresden in a coniferous forest.For both sites Level 4 data, providing information about GPP, are available for 2000 and 2001.GPP was calculated by subtracting the estimated ecosystem respiration from measured NEP) as described in Reichstein et al. (2005).These data were crosschecked against BETHY/DLR's modelled GPP, as annual sums at each station (Table 3).
Because BETHY/DLR underestimated annual GPP sums by 20 %-30 % for both stations over both years, a further analysis of monthly GPP sums was performed (Fig. 2).
BETHY/DLR's GPP estimates qualitatively follow the measured GPP for the coniferous forest at Tharandt over both years, but with a slight underestimation (Fig. 2a and b).This finding accords with the results of Wisskirchen (2005), who also observed good agreement between measured and modelled NEP at Tharandt.For the mixed deciduous forest of the Hainich station in 2000, BETHY/DLR's estimated GPP is good up to May and from September onwards.However, June, July and August exhibit a strong depression of the modelled GPP, responsible for the overall underestimation for this year (Fig. 2c).For Hainich in 2001 the modelled GPP is  overestimated in May, but underestimated in all other months in the growing season (Fig. 2d).
Overall, BETHY/DLR models monthly GPP well for coniferous forest, with an underestimation of yearly GPP of less than 30 %.For deciduous forest the underestimation of yearly GPP is also less than 30 %, but monthly GPP shows an RMSE of up to 25 (Table 3).In particular, BETHY/DLR seems to strongly underestimate the peak of productivity during the middle of summer, perhaps due to the more heterogeneous development patterns of the species in this vegetation group.
Comparisons of GPP calculated by other vegetation models, such as BIOME-BGC, ORCHIDEE and LPJ, with eddy covariance flux tower measurements revealed an overall RMSE of 30 % for forests in climate zones from boreal to Mediterranean (Jung, 2008).This magnitude of error corresponds well with our findings.

Validation strategy
To validate BETHY/DLR's modelled NPP, we used empirical data on the mean annual increment (MAI) of timbergrowing stocks.These data describe the change in the aboveground woody parts of trees with a diameter greater than 7 cm, given as solid tree volume (including bark), estimated during the second National Forest Inventory (NFI2) and the first National Forest Inventory (NFI1), divided by the time between the two inventories.We chose to use MAI instead of current annual increments (CAI) because for Germany no empirically-derived data with higher temporal resolution than the NFI is available.For Germany, these timber stock increment data are provided by the second National Forest Inventory, classified by forest, tree species, and age class.The aim of this ongoing large-scale survey is to establish statistically reliable central monitoring of the development of Germany's forests, to allow assessment of each forest's status and production potential.The NFI survey uses a permanent set of sampling sites, based on a nationwide 4 km × 4 km grid.Samples are taken at randomly chosen sites from this set, using a uniform procedure across all of Germany.This sampling procedure fulfils accuracy requirements at the national level, but more intensive sampling is conducted for greater accuracy in some smaller federal states, such as Bremen and Hamburg.From a statistical point of view the maximum error of the NFI survey regional level data is about 12 % for coniferous forest and below 8 % for deciduous forest, calculated from the data of BMELV (2004).
For each sampling site, the NFI surveys about 150 characteristics, such as tree species composition, timber volume, and growth, using an angle-count sampling method at each corner of the site.Furthermore, sampling circles with defined radii are drawn to survey tree species composition, tree heights, deadwood, ground vegetation and other characteristics.A more detailed description can be found in BMELV (2004).
The first NFI survey was conducted between 1986 and 1988 and was restricted to the ten states of the former West Germany.All data of the NFI1 are referenced to the year 1987.The second NFI survey was carried out in 2001 and 2002 across all federal states of Germany.The reference year for NFI2 is 2002.The NFI data are freely available at the NUTS-1 level (BMELV, 2004).
In order to validate the modelled NPP against these NFI surveys, the highly detailed NFI data had to be aggregated.In a first step, the tree species reported in the NFI statistics were grouped into BETHY/DLR's two temperate forest classes, coniferous and deciduous.Coniferous forest in Germany is mainly composed of spruce (∼57 %) and pine (∼33 %), while Germany's deciduous forest is dominated by beech (∼48 %) and oak (∼25 %).It was therefore assumed that all parameters (standing timber stock and mean annual increment of timber stock, in particular) of these two forest classes could be estimated as the sum of the metrics for the two principal species, plus an estimated value for all of the secondary species.
The NFI data were then used to calculate the MAI of the total above-ground biomass (MAI T ), which was then compared to the modelled NPP.To calculate MAI T for a NUTS region, the increments of the above-ground biomass of coniferous (MAI c ) and deciduous (MAI d ) trees were summed (Eq.1).
As described above, we calculated MAI c and MAI d as the sum of the increments of the total above-ground biomass of the two dominan species ( M) plus an estimate for the tertiary species ( M).Both MAI c and MAI d , represented simply as MAI, were thus calculated following Eq.( 2).
The index s represents the two dominant species of the forest class.
Since tree biomass depends upon age, M was calculated using per-species age-dependent biomass increments, β, for the ten age classes in the NFI data, as shown by Eq. ( 3).
The index a represents the age classes, and A represents the proportional area occupied by each age class (given in the NFI data).
Since values for β are not given by the NFI, the NFI's species-and age-dependent net increments of the outer bark volume V , expressed in m 3 a −1 area −1 , are used.To convert timber stock biomass into absolute dry biomass, speciesdependent conversion factors (ε) are needed (Table 4).Furthermore, species-and age-dependant conversion factors, P , are needed to estimate the total above-ground biomass.P represents the fraction of total above-ground biomass, including branches of less than 7 cm diameter and needles/leaves, that is outer bark (Table 4).Given these values, β was calculated as shown in Eq. (4).
Because the first NFI was conducted before 1989, values for V are only available for the ten states of the former West Germany.For the five states of the former East Germany, values for V were taken from yield tables (Erteld, 1963).
For this study values of ε were taken from Dieter and Englert (2001).Since only a single source of values for P was found, the accuracy of these values is unknown.
To calculate M it was assumed that a weighted average of the biomass increments of the two dominant tree species, M, would be representative for the tertiary species (Eq.5).
www  Statistical analysis indicated that the error of this approach is less than 6 %.
P and Ṽ represent the arithmetic averages of the weighted average means of P and V , respectively.Since both depend on tree age and species composition, they were thus calculated using Eq. ( 6) ( P being calculated analogously using P ).
Before the modelled results could be validated, they needed to be aggregated to the NUTS-1 level at which the NFI data are given.To accomplish this, the modelled NPP was first transferred to a GIS, taking into account the equi-rectangular WGS84 datum map projection.Then, in order to allow comparison of the datasets, the CAI of the above-ground biomass of the modelled NPP was calculated.Because the modelled NPP does not specify which parts of the plant contain the accumulated carbon, the below-ground carbon content had to be estimated and removed.Furthermore, the NPP (in units of carbon) was converted to above-ground biomass (in units of dry weight) by applying a conversion factor (Eq. 8).
R represents the ratio of the increment of below-to aboveground biomass, while λ is a conversion factor from NPP to total biomass.Species-specific values for R were taken from Pistorius and Zell (2005).Since the GLC2000 gives no information about tree species distribution, a mean value High NPP is shown in green, moderate NPP in yellow, and low NPP in red.Grey pixels represent areas which do not belong to the GLC2000 classes glc-2, glc-3, glc-4 or glc-6 (see also Table 1).Blue pixels represent pixels designated as forest in GLC2000, but that have less than 10 % forest cover according to VCF; their modelled NPP is therefore close to zero despite being considered forest.
of R for each of the two temperate forest classes modelled by BETHY/DLR (coniferous and deciduous) was calculated (deciduous: 0.19 ± 0.08, coniferous: 0.23 ± 0.04).To check the numbers upon which these calculations were based, the corresponding allocation factors for above-ground biomass were also calculated using the same dataset (deciduous: 0.81, coniferous: 0.84); these values agree closely with previously reported values (Zhou et al., 2006), supporting our estimated values for R. The value for λ was set to 2, which is seen as representative for both deciduous and coniferous trees (Houghton et al., 1997).
After these calculations and conversions, the CAI derived from the BETHY/DLR model output was compared to the MAI calculated from NFI statistics in order to validate BETHY/DLR's estimates of NPP.

Results and discussion
Figure 3 depicts annual modelled NPP for Germany in 2000 and 2001 at this study's spatial resolution of 1 km 2 .For both years the major forests of southern Germany (the Spessart, the Palatinate, and the Black Forests) are clearly identifiable from their high NPP values, whereas the northern forested areas, such as the Harz mountains, have substantially lower NPP.The mean annual NPP for 2000 was 139 (tC km −2 yr −1 ) with a maximum of 547 (tC km −2 yr −1 ); for 2001 the mean annual NPP was 137 (tC km −2 yr −1 ) with a maximum of 544 (tC km −2 yr −1 ).Total annual modelled NPP was thus 21.6 × 10 6 tC for 2000 and 21.3 × 10 6 tC for 2001.
Conversion of these NPP values to above-ground biomass as described above gives annual totals of 52.3 × 10 6 t for 2000 and 51.8 × 10 6 t for 2001.The value estimated from NFI's data is 82.7 × 10 6 t (for both of 2000 and 2001).Our modelled NPP thus shows an underestimation of 37 % for both years compared to empirical data.Furthermore, large areas with very low NPP can be identified, especially at the borders of larger forests such as the Black Forest of southwestern Germany.This is because these areas are considered to be forest according to GLC-2000, but MODIS VCF indicates very low forest cover fractions (down to one part per thousand).Such areas of conflicting land cover information are shown as blue pixels in Fig. 3.
Linear regressions of estimated above-ground biomass increment from modelled NPP (CAI) against the empirical data from the NFI (MAI) are presented in Fig. 4, with results for coniferous and deciduous trees separated for comparsion.
Figure 4 shows that BETHY/DLR underestimates the net increment of above-ground biomass for both deciduous and coniferous trees.The R 2 values of 0.74 and 0.76 for deciduous trees indicate a high degree of correlation, however.The correlation for coniferous trees is even stronger, with R 2 values of 0.95 and 0.93, but the underestimation is also higher here.In order to quantify the predictive accuracy of BETHY/DLR's NPP estimates, the root mean square error (RMSE) was calculated for all four panels; for deciduous trees the RMSE is 1.53 (2000) and 1.48 ( 2001), and for coniferous trees, 1.87 (2000) and 1. 93 (2001).
The MAI of above-ground biomass for deciduous trees in Germany for all NUTS-1 regions is 821.9 tons per km 2 (NFI), but the corresponding value estimated with BETHY/DLR is 530.9 tons per km 2 , 35 % less.For coniferous trees these values are 804.7 tons per km 2 (NFI) and 416.0 tons per km 2 (BETHY/DLR), a 48 % underestimate.
A reason for underestimation can be found in the land cover/land use classification used (GLC2000).Figure 5 presents a comparison of the forest areas derived from the NFI database, the forest areas of GLC2000, and the forest areas for the intersection of GLC2000 and MODIS VCF.
From Fig. 5b it is apparent that NFI and GLC2000 deciduous forest area estimates differ markedly; an underestimation of 66 % for Schleswig-Holstein (SH) and an overestimation of 106 % for Bavaria (BY) is observed, for example.For coniferous forest (Fig. 5a) the two area estimates are more comparable, with a mean difference of 20 %.These imbalances are reduced when looking at the total forest areas for deciduous and coniferous trees across all of Germany; total coniferous forest area estimates are 42 400 km 2 (NFI) and 47 100 km 2 (GLC2000), and for the deciduous forest, 60 800 km 2 (NFI) and 61 100 km 2 (GLC2000).It can also be seen in Fig. 5a and b that the GLC2000 underestimates forest areas for the northern states of Germany such as Schleswig-Holstein (SH) and Lower Saxony (NI), whereas it overestimates the forest areas for southern states such as Bavaria (BY) and Baden-Württemberg (BW).
In aggregate, then, GLC2000 represents forest area well, but its spatial distribution is not comparable with the NFI data.We hypothesize that the medium-scale forest structure found in most parts of Germany is not adequately described by the GLC2000, due to the difficulty of accurately classifying a heterogeneous land cover distribution even with a resolution of 1 km 2 (Mayaux et al., 2006).According to the Land Cover Classification System (DiGregorio and Jansen, 2001) used in deriving GLC2000, the GLC2000 class "Broadleaved Deciduous Closed to Open (100-40) % Trees" includes all forest areas with a forest fraction from 40 % to 100 % -a very wide range.In order to describe the forest cover fraction more precisely, the MODIS VCF product was combined with the GLC2000 to produce the area estimates used as inputs by BETHY/DLR.Figure 5a and b show the coniferous and deciduous forest areas that result from this combination of MODIS VCF and GLC2000.Clearly this approach led to underestimations of the forest area in Germany, both for coniferous (47 %) and deciduous (59 %) forest.This underestimation occurs because only areas reported as forested in the GLC2000 were carried forward to be combined with the MODIS VCF coverage data; areas designated as non-forest in GLC2000, but with a non-zero forest cover fraction in VCF, were treated as non-forested.As a result, those forest areas which were underestimated by GLC2000, such as Lower Saxony or Schleswig-Holstein, led to substantial underestimations of the increment of above-ground biomass.
Other land cover datasets with higher resolution, such as CORINE (100 m × 100 m) and MERIS GlobCover (300 m × 300 m), are available for Germany, and their land use structures show a better agreement with the NFI data.Since BETHY/DLR requires land cover and LAI inputs to be at the same spatial resolution, and since no higher-resolution LAI products are available yet for Germany, these finergrained land cover datasets unfortunately could not be used.Exploratory analysis shows, however, that the combination of GlobCover and VCF leads to an underestimation of forest area of 24 %, while the combination of CORINE and VCF yields an underestimation of only 7 %.This agrees with the findings of EEA (2006), which estimated the reliability of the CORINE classes 311 (coniferous forest) and 312 (deciduous   forest) at better than 85 %.We observe, therefore, that while area-wide land cover products at high resolution are needed and useful, high-resolution datasets for plant physiology parameters such as LAI must also keep pace if these products are to be of maximal utility.
Returning to Fig. 4, we note that when the number of observations is small, the slope of a regression line is very sensitive to outliers.In the case of deciduous forests in 2000 and 2001 two outliers can be identified that have a strong influence on the slope of the regression line: Bavaria and Baden Wuerttemberg, the two largest federal states of Germany.In Fig. 4 these states have the largest values (on both axes), because of their large areas.
To compensate for the large effect of both: the potential outliers and forest area underestimation, we normalized the CAI and MAI data for each NUTS-1 region by dividing these values by forest area, resulting in units of tons per km −2 per NUTS unit.For this the MAI data was devidid by the forest area reported in the NFI and the CAI data by the combined area of GLC2000 and VCF as used for the model run.The results are presented with linear regressions in Fig. 6.
Figure 6 shows that BETHY/DLR does not exhibit an underestimation for coniferous trees with these areanormalized metrics, indicating that the underestimation seen in Fig. 4 might indeed be explained by GLC2000's area underestimation, as discussed above.The R 2 values of 0.61 and 0.53 still indicate a reasonable degree of correlation.
The underestimated CAI for deciduous trees, observed in Fig. 4 and persisting in Fig. 6 after area-normalization, might be explained by BETHY/DLR's internal model parameters related to carbon uptake: maximum carboxylation rate and maximum electron transport rate.The values used for these parameters were taken from Knorr (1997) (see also Knorr and Heimann, 2001), where they were used for global carbon assessment; these values thus represent global mean values.However, forests in Germany are probably more productive than the "global mean" trees simulated by BETHY/DLR using these parameter values, because of their age.The last large reforestation programme in Germany followed World War II, to mitigate the deforestation experienced during the war.According to the NFI, the mean age of Germany's forests is about 67 yr (81 yr for deciduous and 54 yr for coniferous trees), an age class that is expected to exhibit a high rate of increase of timber biomass.Young and old trees differ in their carbon allocation and fixation strategies; in particular, carbon fixation and timber growth decreases with increasing tree age.In old trees, the maintenance respiration rate is nearly as high as the carbon uptake rate, and thus the large majority of GPP in older trees is dedicated to maintenance.The carbon uptake of young trees, on the other hand, is mainly used for growth.Studies show that the transition between these two metabolic regimes occurs at about 60 to 80 yr of age (Zhou et al., 2006).Therefore it is likely that the values used for the maximum carboxylation rate and the maximum electron transport rate are too low to accurately simulate the tree communities of Germany (see Zaehle et al., 2006 for further discussion of this issue).
Underestimations in the modelled NPP could also be the result of the neglect of nitrogen deposition in the model.Luysseart et al. (2010) showed an increase in modelled NPP of up to 30 % when nitrogen deposition is included in the model formulation.
Uncertainties for deciduous trees are higher (R 2 0.35 and 0.37) than for conifers, which might indicate higher structural variability among deciduous tree species.In particular, deciduous tree species exhibit greater variation in shape than do coniferous tree species in nature.
Finally, it should be taken into account that NFI statistics can only produce MAI values; these values were estimated from the difference between the first NFI survey, in 1987, and the second, in 2002.Until NFI conducts a third survey year, the effects of climatic variability cannot be captured by the NFI statistics.

Estimation of energy potentials
A further objective of this study is to derive energy potentials both from modelled NPP and from NFI data.The energy potential of forests is of considerable importance to the sustainable energy discussion and the development of sustainable energy policy.
To estimate theoretical energy potentials, species-specific lower heating values (H ) can be used to convert from absolute dry above-ground biomass.Heating values represent the maximum energy obtainable from the combustion of biogenic solid fuels, and are given in megajoules per kilogram.Since the GLC2000 gives no information about tree species, but does differentiate between deciduous and coniferous trees, mean heating values representative of each tree class were used.For this study, heating values for deciduous and coniferous trees were calculated for each NUTS-1 unit (Table 5), taking into account the relative abundance and age distribution of tree species in each region.Heating values for the main deciduous (oak and beech) and coniferous (pine and spruce) tree species can be found in many sources; see, for example, Kaltschmidt and Hartmann (2009) 3, areas which are not designated as forested by GLC2000 (see Table 1).
and Grammel (1989).We assume sustainable management of forest biomass; following this assumption, only the mean annual net increment of forest biomass is used to calculate theoretical energy potential.
Since B (above-ground biomass; see Eq. ( 8) is determined for absolute dry conditions, we calculated the energy potential (J ) as shown in Eq. ( 9).
Using this equation, theoretical energy potentials for 2000 and 2001 were estimated from the above-ground biomass previously calculated from modelled NPP (Fig. 7).Although the validation previously conducted had demonstrated a systematic underestimation of NPP, no correction was applied here to compensate for this, since such correction would have resulted in an incorrect spatial distribution of estimated energy potentials.
Comparison with Fig. 3 shows that those forest areas having the highest NPP values also have the highest theoretical energy potentials.This is also true, mutatis mutandis, for low NPP and energy potential, and is valid for both 2000 and 2001.Analysis revealed that the mean theoretical available energy potential for coniferous forest is 17.5 (TJ km −2 yr −1 ) for 2000 and 2001, while for deciduous forest these values are 25.0 (TJ km −2 yr −1 ) (2000) and 24.6 (TJ km −2 yr −1 ) (2001).Maximum values of 25.7 (TJ km −2 yr −1 ) (coniferous forest) and 25.4 (TJ km −2 y −1 ) (deciduous forest) were found.
The NFI data were also used to estimate empirical energy potentials, using Eq. ( 10) for all thirteen NUTS-1 units in Germany (Table 5).These estimates of energy potential are partitioned into the two main tree classes.
Since it was shown previously that underestimation existed both for the areas of forests (as derived from remote sensing) and for modelled NPP, it is unsurprising that the empirical energy potential estimated from NFI data is 37 % higher than the theoretical estimate from modelled NPP.

Conclusions
For this study we modelled the Net Primary Productivity (NPP) of German forests for 2000 and 2001 using the dynamic biomass model BETHY/DLR.We presented a new method for the validation of modelled NPP using empirical data related to MAI of above ground biomass.Modelled NPP was converted and aggregated to a net increment of above-ground biomass per NUTS-1 unit for comparison to these empirical data.With this method we showed a high degree of correlation between modelled and empirical NPP, although the modelled NPP underestimated the empirical NPP.A comparison with data from two eddy covariance flux towers revealed that BETHY/DLR represents annual productivity patterns well (particularly for coniferous forest) but with substantial underestimation.
In a second step, the sustainable theoretical energy potential of the above-ground biomass was estimated, using heating values to convert estimated above-ground biomass to energy units.For comparison, energy potentials were also calculated from empirical data, which revealed that modelled energy potentials are underestimated by 37 %, a consequence of the prior underestimation of modelled NPP.
Reasons for this pattern of underestimation were discussed; in particular, it was shown that GLC2000 does not represent the spatial distribution of forest areas well due to its limited resolution.We thus argue that 1 km 2 resolution is insufficient to describe the heterogeneous small-scale structure of mid-European forests.For future modelling, the use of higher-resolution land cover products might allow more accurate NPP estimation; this should be tested in future research.To facilitate the use of such products, however, there is a need for matching high-resolution datasets for vegetation metrics such as LAI.
Furthermore we hypothesize that the maximum carboxylation rate and maximum electron transport rate is agedependent and thus potentially responsible for the underestimation of modelled NPP, since BETHY/DLR does not take the tree age distribution into account.Further research in this area could thus lead to more exact results.
Modelled NEP, NPP and GPP are typically validated using eddy flux tower measurements, but such measurements are not available in many areas.Our presented validation method could therefore be helpful in the assessment of model outputs both at a broader spatial scale, and in less developed countries.Our method will, additionally, allow the development of a downscaling procedure for empirically derived NUTS-level data, allowing NUTS data to be partitioned into smaller spatial units.Our MAI-based validation method should, however, be tested in additional countries, and should be comprehensively compared to validation using eddy measurements, so that the benefits and drawbacks of each method are clearly understood.
This new MAI-based validation method will be useful in validating modelled NPP; as we demonstrated here, that also allows the further estimation of other metrics, such as bioenergy potentials.Such estimates of forest energy potentials play an important role in planning for a sustainable economy.More broadly, accurate and precise model results, crosschecked against empirical data, are needed for a better understanding of optimal forest management and the future possibilities of renewable energy.

Fig. 2 .
Fig. 2. Monthly GPP sums modelled with BETHY/DLR (solid lines) and measured with eddy flux towers (dashed lines) for the Tharandt and Hainich stations for 2000 and 2001.

Fig. 3 .
Fig. 3. Annual NPP of forest areas in Germany for 2000 and 2001.High NPP is shown in green, moderate NPP in yellow, and low NPP in red.Grey pixels represent areas which do not belong to the GLC2000 classes glc-2, glc-3, glc-4 or glc-6 (see also Table1).Blue pixels represent pixels designated as forest in GLC2000, but that have less than 10 % forest cover according to VCF; their modelled NPP is therefore close to zero despite being considered forest.

Fig. 4 .
Fig. 4. Estimated above-ground biomass increment from modelled NPP (CAI) versus empirical data from the NFI (MAI) for Germany's deciduous and coniferous trees for 2000 and 2001.Each cross represents one NUTS-1 region.Thick lines show linear regressions.Values are given in megatons per NUTS-1 unit per year.

Fig. 6 .
Fig. 6.Normalized estimated above-ground biomass increment from modelled NPP versus empirical data from the NFI for Germany's deciduous and coniferous trees for 2000 and 2001.Each cross represents one NUTS-1 region.Thick lines show linear regressions.Values are given in tons km −2 .

Fig. 7 .
Fig. 7. Sustainable theoretically available energy potential, in terajoules per 1 km 2 pixel, of forest areas in Germany for 2000 and 2001.Low energy potentials are shown in blue, intermediate potentials in beige, and high energy potentials in red.White represents, as in Fig.3, areas which are not designated as forested by GLC2000 (see Table1).

Table 1 .
Translation of GLC2000 vegetation classes to BETHY/DLR vegetation types with weighting factors.

Table 2 .
Summary of meteorological input data from ECMWF, including short names and code numbers.
resolution of up to four times per day.The data are produced by ECMWF's climate model analysis of meteorological station measurements of air temperature (at 2 m height), wind speed (at 10 m height), soil water content (in the four uppermost layers), and cloud cover.Daily precipitation values were obtained from the ECMWF INTERIM dataset.

Table 3 .
Annual GPP sums in (gC m −2 a −1 ) modelled with BETHY/DLR and measured with eddy flux towers at the Hainich and Tharandt stations for 2000 and 2001.The percent difference and RMSE are provided for each comparison.