Quantifying the carbon uptake by vegetation for Europe on a 1 km 2 resolution using a remote sensing driven vegetation model

Introduction Conclusions References


Introduction
Since the early 1950s the continuous rise of the atmospheric carbon dioxide concentration is known (Keeling, 1960).The linkage of higher concentrations of climatically active gases in the atmosphere with the threat of sustained global climate warming (Houghton et al., 1996) greatly boosted research activities to quantify the various components of the global carbon cycle.The terrestrial biosphere plays a prominent role in this system and thus can have major impacts on atmospheric CO 2 concentrations, with only minor changes in the terrestrial productivity (e.g.Fowler et al., 2009;Friedlingstein et al., 1994).
The quantification of the carbon exchange between biosphere and atmosphere on regional and global scales currently relies on modelling approaches.To this effect net primary productivity (NPP) is estimated describing the amount of CO 2 accumulated by vegetation taking into account the autotrophic respiration of the vegetation.The net carbon exchange between vegetation and atmosphere is calculated by taking also into account the heterotrophic respiration and NPP, called net ecosystem productivity (NEE).NEE can be determined by subtracting the heterotrophic respiration from NPP.
To compute NPP, currently two main modelling approaches are established.Models which use the concepts of Monteith (1965) and Monsi and Saeki (1953) assume that NPP can be calculated as the product of photosynthetic active radiation (fPAR) and the light-use efficiency.Examples are the models EPIC (Environment Policy Integrated Climate) by Williams et al. (1984), CASA (Carnegie-Ames-Stanford Approach) by Potter et al. (1993), and C-Fix by Veroustraete et al. (1994).Another important modelling approach follows the approach of Farquhar et al. (1980) who introduced a biophysical model of photosynthesis.This concept is often used for Soil-Vegetation-Atmosphere Transfer (SVAT) models and dynamic vegetation models as for instance the BIOME3 (equilibrium terrestrial biosphere) model by Haxeltine and Prentice (1996), the LPJ (Lund-Potsdam-Jena) model by Sitch et al. (2003), and the ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosystems) model by Krinner et al. (2005).
Most of the just mentioned models are designed for global applications and are run on coarse resolution (e.g. 1 • × 1 • ) due to the lack of spatial high-resolution input data or, for regional applications on e.g. 1 km 2 resolution.However, only few studies exist for Europe that simulate the carbon exchange of all vegetation types.Vetter et al. (2008) used seven models (five processed-based terrestrial ecosystem models: Biome-BGC, LPJ, ORCHIDEE, JULES and PIXGRO; and two data oriented models: MOD17+ and NETWORKANN) to investigate GPP and NEP time series in Europe .They put a special focus on the year 2003, for which they found consistent patterns in all seven model results.Veroustraete et al. (2002) used the C-Fix model to calculate NEP for the period 1997-1999 and compared their results with Euroflux data.They resulted in a total NEP for Europe of 2.7 PgC a −1 .Jung et al. (2007) performed a study to assess gross primary productivity (GPP) using three different process-oriented biosphere models.The study showed that uncertainties in meteorological input data have significant effects on the model response (typically 9 % relative variability).Total GPP values for Europe of 7.1 to 8.7 PgC a −1 (including parts of Russia) can be found in Jung et al. (2008) and of 3.9 to 5.8 PgC a −1 (excluding Russia, Belarus and Ukraine) in Beer et al. (2007).These findings indicate a high amount of variability with respect to the use of the model, the input data and the spatial resolution.
Uncertainties in model results and their sources were investigated by Lin et al. (2011).They focussed on boreal forests and found that biases in downward shortwave radiation data have most significant effects and result in a systematic error of the annual carbon uptake.However, the largest cumulative errors resulted from biospheric parameters controlling the light-use efficiency and respiration-temperature relationships.
A further possibility to quantify carbon fluxes and/or carbon stocks is to use remote sensing derived products, either directly or as drivers for vegetation models.Medium resolution remote sensing gives the opportunity to monitor the earth's surface, globally, every 3 days.Thus it has a potential of immense information gain on a global and regional level.The Moderate-resolution Imaging Spectroradiometer (MODIS) product MOD17 provides continuous time series of GPP and NPP on a 1 km 2 resolution with a near real-time 8 day interval (Zhao et al., 2005).Vegetation models driven by remote sensing as e.g.BETHY/DLR (Biosphere Energy Transfer Hydrology; Tum et al., 2011) bring further advantages, because they couple measured data with complex modelling approaches.With such approaches data assimilation is feasible in principle.
Model validation is often conducted using data from eddy covariance flux tower measurements.The relationship between carbon exchange and energy flux has been studied before in international networks such as FLUXNET (Baldocchi et al., 2001), 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 NEE at the spatial scale of the footprint of a tower (Baldocchi, 1997).
With this study we will contribute to the on-going discussion of the net carbon exchange balance in general and with special focus on Europe.We will investigate whether BETHY/DLR is suitable for monitoring carbon fluxes by comparing against independent data and present new GPP estimates for Europe for the period 2000-2007 on a 1 km 2 resolution.These were calculated with a modified version of our BETHY/DLR model, which was developed by Knorr (1997) (see also Knorr, 2000;Knorr and Heimann, 2001) and adapted for the use of remote sensing data by Wißkirchen (2005).Furthermore, we compare our results with GPP data taken from the FLUXNET database and state on the robustness and uncertainty of both: model results and measurements.

Input data
BETHY/DLR is designed to be driven by remote sensing data.The condition of vegetation is described by leaf area index (LAI) time series derived from optical sensors (e.g.SPOT-VEGETATION, ENVISAT-MERIS, Terra/Aqua-MODIS).Vegetation related parameters such as the Normalized Differenced Vegetation Index (NDVI), the Fraction of Absorbed Photosynthetic Active Radiation (FAPAR), or LAI are widely used in vegetation and ecosystem models (e.g.Wißkirchen, 2005;Veroustraete et al., 1994) and can easily be calculated using existing approaches (Gobron et al., 2008;Baret et al., 2007).For this study we used data from the CY-CLOPES database.This database contains global LAI time series for the period 1999-2007 on a 1 km 2 resolution, which are offered as 10 day composites in 10 • × 10 • tiles.The data are freely available.
Since time series data derived from optical satellites are often contaminated by gaps and outliers, time series analysis has to be applied to eliminate these data errors.We used a modified harmonic analysis (HA) for un-evenly spaced data sets based on the least-squares technique (Lomb, 1976;Bittner et al., 1994).Our HA decomposes time series by successively subtracting the highest peak in the power spectrum, then computing a new spectrum etc., resulting finally in a linear combination of trigonometric functions, i.e. sine and cosine oscillations.For each harmonic component, the amplitude, frequency, and phase are found by least-squares fitting approaches.Before applying HA, outliers are determined using criteria for unrealistic fluctuations of vegetation growth.In addition, large gaps (more than 5 missing composites) are filled using a land cover and regional-specific averaged LAI time series.This approach yields a spatio-temporally continuous data set as needed for BETHY/DLR.The general advantage of using remote sensing derived data sets in contrast to standardized growing functions is the possibility to better represent local phenological conditions.However, a combination of both (measured and modelled data) can easily be used for a data assimilation approach (Knorr et al., 2010).
Aside from the LAI time series, the CYCLOPES database also provides a land cover/land use classification (LCC) data set -the Global Landcover Classification 2000 (GLC2000), which was also derived from SPOT data.GLC2000 is representative for the year 2000 and contains 22 global land cover classes.The Land Cover Classification System (LCCS) of the FAO (Food and Agriculture Organization) was used (Bartholome and Belward, 2005;DiGregorio and Jansen, 2001) to derive these classes.With the current model setup it is in principle possible to apply any LAI time series and LCC product, irrespective of their source, but preferably with the same spatial resolution.Depending on the availability of input data sources the current model output is in 1 km 2 resolution but foreseen to be increased to 300 m × 300 m resolution.
BETHY/DLR is usually run with 10 day LAI data from SPOT/VEGETATION and GLC2000 land cover classification.However, any other data sources with adequate spatial and temporal resolution can be used.For example, the usability of MODIS LAI and MODIS land cover data sets was tested for modelling NPP in Kazakhstan (Eisfelder et al., 2013).The most relevant adaptations of the model include the frequency of LAI data update interval to 8 days (to meet the frequency of the MODIS LAI data), and the transformation of the land cover classes to BETHY/DLR internal vegetation classes.The usage of 16 day LAI input data for a test site in semi-arid central Kazakhstan indicated that this longer time step might be too coarse to capture variable meteorological conditions (Eisfelder et al., 2013).
Besides information about the condition of vegetation, the model is driven by meteorological time series, as precipitation, temperature, wind speed and radiation.Knorr (1997) used a weather generator on a daily basis to predict precipitation.Temperature was scaled linearly from monthly to daily averages.Wind speed was considered as constant (3 m s −1 ) and PAR was calculated following the approach of Weiss and Norman (1985) using the solar elevation, earth-sun distance and solar flux, computed from geographical latitude, Julian day and solar hour.
For our study we used data from the operational process chain of the European Centre for Medium-Range Weather Forecasts (ECMWF).The ECMWF re-analysis project (ERA-Interim) contains daily (4 times daily) data of the temperature at 2 m height, wind speed at 10 m height, soil water content (in the four uppermost soil layers) and cloud cover.Twice-daily data of precipitation is also available.All data sets have a spatial resolution of 0.25 • × 0.25 • .From these data sets, the daily mean, minimum, and maximum temperatures and the water vapour pressure are calculated.Temperatures are scaled linearly to 1 km 2 resolution, by using the difference between the ECMWF reference height and global SRTM (1 km 2 ) elevation data, and the temperature gradient of the US Standard Atmosphere, which is −0.65 K/100 m: T ECMWF represents the reference temperature of ECMWF, h ECMWF is the ECMWF reference elevation and h SRTM the elevation of the SRTM.
The ECMWF data set also includes estimates of daily PAR.However, this data set is not used.The analysis of Wißkirchen (2005) showed that calculating PAR following the approach of Burridge and Gadd (1974) yields very highly reliable results.Thus, daily fractions of high, medium, and low cloud cover from ECMWF are used to calculate atmospheric transmission.Using the geographical coordinates, the Julian day and year and the atmospheric transmission result in more accurate estimates of PAR at 1 km 2 spatial resolution (Wißkirchen, 2005).
Soil type information is also needed.We used the data set of Batjes (2006), which is based on the FAO-Unesco soil map of the world and the soil profile database of IS-RIC's (International Soil Reference and Information Centre) global WISE (World Inventory of Soil Emission potentials) database.It is freely available as grid file with the spatial resolution of 5 arc minutes × 5 arc minutes and contains information on 128 FAO soil types, including sand, silt, and clay content, layering and depth.

Processes
The approach to model GPP, NPP, and NEP with BETHY/DLR has previously been described in Knorr and Heimann (2001).However, for our study some processes were refined to make the model applicable for regional applications.Since most areas in Europe cannot be described With this method, carbon fluxes can also be computed for incomplete vegetation cover.Photosynthesis is calculated independently for each vegetation type fraction on an hourly basis, which is aggregated to daily values.For the case that more specific information on the land cover use (i.e.specific tree species distribution) is available, we introduced ten additional crops and tree species.For all species the following parameters are provided: maximum carboxylation rate, maximum electron transport rate, and maximum plant height (see Table 2).
For all vegetation types new maximum rooting depths were taken from Canadell et al. (1996).These values were used to limit the bucket size in its height for the soil water budget.The soil water content was initialized by the cumulative soil water content of the four uppermost layers of the ECMWF soil water on the first simulation day.After initialization, the soil water budget is calculated independently of further input data.Since a transition phase is needed, criteria were defined to reach equilibrium.At first, Wißkirchen (2005) introduced a 1 yr transition phase, assuming stable conditions after this time.However, further studies revealed that, depending on soil type and environmental conditions, transition phases of up to 4 yr are needed.The length of transition is calculated by comparing the soil water content of the first simulation day with the content 1 yr later.It is assumed that for stable conditions, the soil water content of these 2 days does not vary more than 5 %.Therefore, it is  Knorr, 1997).(1993) assumed that if less than 10 % of the simulated pixels show a variation of more than 5 % in their water content, stable conditions are reached.

Preparation of FLUXNET data
In order to compare our results with independent data, we used eddy covariance data provided by FLUXNET.
FLUXNET is an association of regional networks aiming to coordinate regional to global analyses of observations from micrometeorological tower sites.These flux tower sites use eddy covariance methods to determine fluxes as carbon dioxide, water vapour and energy flow between terrestrial ecosystems and the atmosphere by measuring meteorological data (e.g.sonic anemometer wind speed, temperature and densities of carbon dioxide and water vapor).Over 500 towers are operated worldwide from which around 150 are situated in Europe.The European towers were established to measure carbon fluxes over various vegetation types, including temperate conifers, deciduous and evergreen broad-leaved forests, croplands, grasslands, evergreen shrubs, and wetlands.The data are provided as half-hourly time series and available at http://www.europe-fluxdata.eu.It includes several pre-processed data sets starting from level-0 data to calculated values of GPP and ecosystem respiration, which are provided as level-4 products.GPP estimates are based on the measured net flux of carbon from the vegetation/soil to the atmosphere (called net ecosystem exchange, NEE) and on an estimation of the terrestrial ecosystem respiration (R eco ) according to (2) Reichstein et al. (2005) proposed a method to estimate ecosystem respiration by determining nighttime temperature sensitivities of ecosystem respiration and to extrapolate these estimates to daylight period after Lloyd and Taylor (1994): T 0 is the regression parameter with constant value of is the soil temperature and R ref the temperature dependent contribution of respiration.E 0 is the free parameter of activation energy to determine the temperature dependency (Reichstein et al., 2005).Thus, GPP is given by the difference of measured NEE and estimated daytime R eco .Various uncertainties are within this approach and have been previously discussed in literature.Papale et al. (2006) investigated uncertainties of algorithms and parameter estimations, especially during the preprocessing of the eddy flux data.They found a major concern in the heuristic low turbulence u * filtering, which introduces the largest uncertainties, while the quality of storage correction depends on the measured profile of CO 2 concentration.Furthermore, they found that spike removal does not affect directly the annual NEE but can affect the quality of gap filled data sets.Gap filling (Moffat et al., 2007) and partitioning (Reichstein et al., 2005;Desai et al., 2008) of NEE in the two components can also result in uncertainties.Reichstein et al. (2005) showed that temperature sensitivity of R eco derived from long-term data (annual) sets is different to the shortterm temperature sensitivity.Thus, if long-term temperature sensitivity is used for extrapolation to half-hourly daytime respiration of summer active vegetation, a systematic overestimation of R eco of more than 25 % might be realistic for annual timescales.On the other hand, for summer passive vegetation as e.g.found in the Mediterranean regions an underestimation of annual R eco is observed.Similar results were found by Richardson et al. (2006), who investigated the suitability of different respiration models for gap-filling techniques and for partitioning eddy flux measurements to respiration and GPP.They stated that "two of the most widely used models of ecosystem and soil respiration, the basic Q10 model and the "restricted" form of the Lloyd and Taylor model (Lloyd and Taylor, 1994) do a poor job of accounting for observed variation in ecosystem and soil respiration in comparison with other simple models."In addition they found little differences in the annual sum of respiration for some test sites among the models (∼ 75 gC m −2 a −1 ), but high ranges for other test sites (355 gC m −2 a −1 ), which is approximately 40 % of the mean annual respiration.Reichstein et al. ( 2007) reported typical uncertainties in eddy covariance flux measurements of less than 100 gC m −2 a −1 while the total systematic error, due to non-ideal observations and correction procedures, is below 200 gC m −2 a −1 .
For our analysis we requested GPP data from 74 towers, included in the GHG_Europe (Greenhouse Gas Management in European land use systems), CarboItaly and IMECC (Infrastructure for Measurements of the European Carbon Cycle) networks.The towers are located in 18 European countries and are maintained and operated for individual periods by local organisations.A list of all tower sites used in this study is presented in Table 3, and for better visual interpretation in Fig. 1.As can be seen, the complete simulation period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) was covered only by few tower sites (DK_SOR, DE_THA, and CH_DAV).According to the individual availability of FLUXNET data, we performed a model run with BETHY/DLR for each site.Thus a total of 273 complete, measured and modelled years were available.For the model run we translated the reported vegetation type to the corresponding BETHY/DLR type, assuming 90 % coverage for each site.Mixed forests were translated to a mixture of 50 % evergreen coniferous trees and 50 % broadleaved deciduous trees (Table 3).Since half-hourly time series of GPP are not directly comparable with BETHY/DLR's model output, the data was aggregated to monthly and annual sums.Since the model resolution of BETHY/DLR is on a 1 km 2 and the eddy covariance measurements are usually valid for some metres for small towers, or hundreds of metres for large towers, we assumed the grid cell which is situated closest to the geographical position of the eddy tower to be representative.In case that the tower is situated at a border of two grid cells, the western cell was used.

Results
The averaged modelled GPP results as calculated with BETHY/DLR for the period 2000-2007 are presented in Fig. 2. Furthermore, the yearly deviations from the 8 yr mean are given.It can be seen, that the GPP pattern is dependent on time and space and shows high variations for individual regions.This is especially notable for the area of Romania and the western part of the Ukraine where the GPP can vary from year to year up to ±250 gC m −2 a −1 from the long-term mean.To show the regional variation more clearly annual GPP fluxes per country are provided in Table 4. From this table, it can be seen that on average the CO 2 uptake by vegetation for Europe is about 2.5 PgC a −1 (±0.17 PgC a −1 ).The highest GPP for the area of Europe was calculated for 2004 (2.7 PgC a −1 ) and the lowest for 2001 (2.2 PgC a −1 ).In Table 4 we also calculated the coefficient of variation (CV) for each country.CV is the standard deviation divided by the average, which is the relative variability.The GPP's of Spain, Portugal and Turkey show the lowest relative variability (4.6, 6.5 and 7.2 %, respectively) indicating that photosynthesis of vegetation of these countries is relatively stable for the period 2000-2007.On the other hand, the relative variability of GPP of the Netherlands, Estonia, Latvia and Belgium is highest (46.4, 45.6, 43.8 and 43.1 %, respectively).It is obvious that our GPP maps with 1 km 2 spatial resolution allow deriving country-wide GPP even for small countries as e.g.Slovenia.According to FAO statistics, Slovenia has a land area of about 20 140 km 2 which can be represented by a rectangle of about 142 km length.The spatial resolution of other modelling exercises as e.g.those presented by Jung et al. (2007) is 0.25 • and 0.5 • .With this spatial resolution, the country-wide GPP of e.g.Slovenia is represented by only a few grid cells.
In order to compare our results with FLUXNET data, we first calculated annual GPP sums for all stations and modelled years.Statistical analysis revealed for 192 of 273 (70 %) measurements a higher GPP in the FLUXNET data.The range is 71-2766 gC m −2 a −1 (FLUXNET) and 182-2386 gC m −2 a −1 (BETHY/DLR), and is strongly dependent on the vegetation type.Annual mean values and standard deviations for the available years were also calculated (see Table 5).It can be found that on average the standard deviation of the FLUXNET data is 16.1 % of the mean GPP value.For the BETHY/DLR results this value is less than half as large (7.7 %), indicating less variability in the annual results.This finding corresponds with the assumption that data provided by FLUXNET is influenced by more environmental conditions and disturbances.
For a closer look to selected sites that represent main land cover types (i.e.deciduous-and evergreen-broadleaved forest, coniferous and mixed forest, grassland, and agriculture), monthly GPP values were calculated and presented as time series in Fig. 3.We chose results which are exemplary for the main vegetation types, as listed above, and preferably also cover the whole simulation period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007).Thus, results for the sites FI_Kaa, FI_Hyy, IT_Cpz and DK_Sor, which almost cover the 8 yr, are presented.In addition IT_Ro1, which contains missing data at the beginning and at the end, and DE_Geb with five years of data are shown.For all vegetation types a full coverage of the 2000-2007 period is exceptional, in particular for agriculture, no longer time series was available.Depending on the land cover type, it can be seen that the agreement of this comparison highly varies.
As can be seen from Fig. 3, our results for the grassland (a) and coniferous forest (b) sites are comparable with the FLUXNET data.Coefficients of determination (R 2 ) of 0.83 (grassland, GL) and 0.93 (coniferous forest, CF) combined with low RMSE of 14 (GL) and 39 gC m −2 month −1 (CF) were calculated.The average of these values for all grassland sites (23) are 0.71 (R 2 ) and 66 gC m −2 month −1 (RMSE), resulting from 95 site-years.The 19 available CF sites cover 89 site-years with a mean R 2 of 0.79 and a RMSE of 70 gC m −2 month −1 .Here FLUXNET estimates GPP consistently higher.These results show a good agreement of both the seasonal patterns, and the absolute GPP values.Good agreement in the seasonal pattern was also found for the mixed forest (MF) site of BF_Vie (Fig. 3c).Here a R 2 of 0.89 and a RMSE of 40 gC m −2 month −1 are obtained.Similar results are found for the other six MF sites (covering 39 siteyears), resulting in a mean RMSE of 67 gC m −2 month −1 and a R 2 of 0.82.Thus, it can be concluded that seasonal productivity of grassland and coniferous forest modelled by BETHY/DLR are close to the eddy covariance estimates.For comparison: the diagnostic biospheric carbon flux model VPRM (Vegetation Photosynthesis and Respiration Model) driven by satellite data from MODIS and a comparison between simulated and observed hourly NEE data from two eddy covariance flux measurement sites in Ontario and Québec yielded R 2 of 0.58 for mixed forest in Ontario and R 2 of 0.63 for Black Spruce in Québec (Lin et al., 2011).
The other three examples for agriculture, evergreen broadleaved forest (EF), and deciduous broadleaved forest (DF) (see Fig. 3d-f), show less satisfactory agreement.The agricultural sample (Fig. 3f) has a relatively high RMSE of 100 gC m −2 month −1 , but also a low R 2 of 0.65, which is mainly based on deviations in amplitude and length of the vegetative period.In the DE_Geb case BETHY/DLR predicts higher GPP values, compared to the FLUXNET measurements.This pattern is not consistent for all 13 agricultural sites, which we believe is caused by e.g.management practices (tillage, crop rotation, fertilization, etc.), which is not included in the BETHY/DLR scheme.In addition, since eddy towers on agriculture are usually only a few high, but BETHY/DLR's resolution is of 1 km 2 , differences might be due to the resolution.FLUXNET data was available for 46 site-years, resulting in a RMSE of 113 gC m −2 month −1 , combined with a R 2 of 0.58.The two broadleaved forest sites (Fig. 3d, e) also show low R 2 (0.6) with intermediate RMSE (EF: 54 gC m −2 month −1 , DF: 75 gC m −2 month −1 ).26 site-years were available for EF and 28 for DF.The agreement for all EF sites is very low (R 2 : 0.36, RMSE: 74 gC m −2 month −1 ).For most of the five sites, the modelled GPP values are consistently higher in the vegetative active phase and lower in the winter season compared to the estimates provided by FLUXNET.For other sites, however, the comparison shows a contrary pattern.Thus, no clear tendency for evergreen broadleaved forest vegetation can be found.
A closer look at the LAI time series which were used for the IT_Cpz site (Fig. 4) reveals distinct seasonal patterns, which were not expected for evergreen deciduous forest vegetation.It is thus questionable if CYCLOPES data are realistic for this type of vegetation.This finding was already discussed by Garrigues et al. (2008), who also stated that CY-CLOPES LAI is most realistic for grassland and agriculture.Significantly lower GPP estimates of BETHY/DLR during the winter season might thus be explained by uncertainties in the LAI data.Hence, for evergreen deciduous forests, eddy covariance measurements can be seen as more likely to predict realistic GPP values than our model does.
The DF site (IT_Ro1) again shows a clear trend of higher GPP reported by eddy covariance measurements, which is consistent for all eight stations.This results in a mean RMSE of 80 gC m −2 month −1 , combined with an R 2 of 0.78 for all DF sites.The better agreement for CF sites and less good agreement for deciduous sites stands in accordance to the findings of Tum et al. (2011), who used BETHY/DLR to predict solid wood increase for Germany's forests, and validated their results with empirical data of solid wood increase.

Discussion
A comparison of our modelled mean CO 2 uptake by vegetation for Europe (GPP ∼2.5 PgC a −1 ± 0.17 PgC a −1 ) with data from Jung et al. (2007) shows that our result is about 2.7 times lower than the mean multi-model GPP for Europe estimated with three terrestrial biosphere models (LPJ, www.geosci-model-dev.net/6/1623/2013/ORCHIDEE and Biome-BGC) using different drivers as e.g.different meteorological input data and different land cover maps.The mean multi-model GPP was calculated for the time from 1981 to 2000.Results of Jung et al. (2007) show a relative variability of about 11 % while the relative variability of our results is about 7 %.They also showed that differences of more than 2 PgC a −1 are observed if the same meteorological drivers, land cover and spatial resolution are used for different models (ORCHIDEE instead of Biome-BGC).Unfortunately, the comparison of the results from Jung et al. (2007) to ours is difficult because they extended the area of interest from Europe to the western part of Russia.Vetter et al. (2008) investigated the European 2003 carbon flux anomaly using a multi-model approach (as mentioned in the introduction).For determining the anomaly they calculated the long-term GPP average .We calculated the range of multi-model European GPP by including the subsets "north", "west" and "central" from 3.7 to 4.3 PgC a −1 .This range is still higher than our result for Europe.
Ciais et al. ( 2010) estimated the long-term carbon balance of European (EU-25) croplands and its component fluxes over the last two decades using process-based models (CASA, LPJ and ORCHIDEE-STICS) as well as inventories of yield statistics.Only the modified version ORCHIDEE-STICS was able to estimate GPP for European croplands (winter wheat and maize only, ignoring summer C3 crops) with 1360 gC m −2 a −1 .The cropland area for EU-25 is defined in Ciais et al. (2010) as 1.08 × 10 6 km 2 ignoring grassland.Thus a total European GPP of ∼1.5 PgC a −1 can be calculated.Luyssaert et al. (2010) presented a new synthesis, based on a suite of complementary approaches, of the primary production and carbon sink in forests of the European Union (EU-25) member states during 1990-2005.In this paper one can find mean GPP data derived from three models (Biome-BGC, LPJ and ORCHIDEE) representing the forest component fluxes.Together with the estimate of the European forest area (of 1.32 × 10 6 km 2 to 1.55 × 10 6 km 2 for EU-25) the total GPP is ∼ 1.2-∼ 1.4 PgC a −1 (LPJ), ∼ 1.3-∼ 1.5 PgC a −1 (Biome-BGC) and ∼ 1.6-∼ 1.8 PgC a −1 .Taking the result of Ciais et al. (2010) for cropland and of Luyssaert et al. (2010) forest, we find a total European GPP in the range from ∼ 2.7 to ∼ 3.3 PgC a −1 which is in good agreement with our result.
Both the eddy covariance approach to deriving GPP, and BETHY/DLR contain significant uncertainties.Both techniques involve models, their parameters and input data, while the main difference is that the eddy-covariance technique of deriving GPP involves local and BETHY/DLR largescale data, each with their own potential and limitations.Taken those principal uncertainties into account, we find that both the eddy covariance approach to deriving GPP, and BETHY/DLR arrive at comparable results at the sites presented in this study.Comparability of the two GPP data sets is in principal given, but strongly depends on the vegetation type.Concerning the uncertainty in the model results of BETHY/DLR and the eddy covariance flux measurements, further discussion is needed.
To estimate carbon fluxes with the eddy-covariance technique meteorological conditions have to be within a defined range.Unfavourable conditions, such as strong nonstationary and non-turbulence, cannot be used to calculate fluxes and thus need gap filling approaches (Chen et al., 2012).In addition, further environmental conditions (e.g.complex terrain and vegetation distribution) can negatively influence measurements.These influences have formerly been widely discussed by e.g.Chen et al. (2009), Göckede et al. (2004) and Sogachev et al. (2004).Foken and Wichura (1996) discussed a potential of error sources not only in the environmental condition of a tower site (i.e.internal boundary layers, surface layer height, gravity waves, etc.) but also in sensor configurations (e.g.flow distortion, sensor separation, and measuring height).For example, in Fig. 3a and e the FLUXNET data shows two unrealistic high peaks of vegetation activity in the winter seasons (FI_Kaa: month 24 and IT_Ro1: month 46).This pattern can be found for 28 of the 273 years (10 %), indicating errors caused either by the measurement or post-processing.
For some months eddy covariance measurements predict negative values (Fig. 3a, b and f), which is inconsistent with the definition of GPP as computed with BETHY/DLR.Here GPP is defined as the total carbon uptake by vegetation, which is necessarily positive.However, following Eq.( 2) GPP calculated from eddy flux measurements is the sum of measured NEE fluxes and extrapolated ecosystem respiration R eco .In consequence, if ecosystem respiration exceeds incoming CO 2 fluxes, negative GPP values are possible.These negative values are most likely caused by additional ecosystem components, such as anthropogenic or fauna caused inputs to the system, which are indirectly taken into account in the FLUXNET data.Similar findings were made by Mitchell et al. (2009), who found that aside from errors in NEE retrieval, fundamental problems in modelling approaches can lead to discrepancies in the data comparison.Thus, low correspondences of modelled and measured NEP (GPP) do not necessarily state on the validity of both data sources.Law et al. (2001) underpinned this finding and noted, "errors in the approaches to estimating NEE (eddy covariance approach) and NEP (biological approach) are large, but combining biological and eddy flux data is useful for model testing." A further error source which has to be accounted is the parameterization of BETHY/DLR's photosynthesis.We used parameters as described in Knorr and Heimann (2001) and Table 2. Since these are generalized, improvement could be archived when they are regionalized or scaled using assimilation techniques.Kato et al. (2013) showed for semi-arid woodland in Botswana that improvements in the terrestrial water and carbon simulations can be achieved if satellite and eddy covariance (FAPAR and latent heat flux, respectively) data are assimilated simultaneously.They found that the parameter describing photosynthesis, as e.g.maximum catalytic capacity of Rubisco (V 25 max ), and the standard ratio of CO 2 concentration inside and outside the leaf tissue for C3 plants as well as the parameter describing hydrology and the coupling of water and carbon fluxes, as e.g. the maximum plant-available soil moisture (W max ), showed high uncertainty reductions after assimilation resulting in lower RMSE between simulated data and observation-based values.It became obvious that V 25 max for the two PFTs (plant functional type) considered are reduced by assimilation to minimize the cost function.From a photosynthetic point of view a reduction of V 25 max reduces photosynthesis when the enzyme Rubisco is limiting; but due to the combined increase of W max GPP is increased during the vegetative phase.However, before assimilation is applied it has to be proven that general characteristics of carbon fluxes correspond with measured data.Friend et al. (2007) used the dynamic vegetation model ORCHIDEE and FLUXNET data (here latent and sensible heat as well as net radiation) to optimize 12 model parameters for a pine forest in France.They found that the optimized carboxylation rate V 25 max is about 1.5 times the a priori value of about 35.The a priori value is identical with the BETHY/DLR value for temperate broad leaved deciduous trees (PFT = 4).For most C3-PFTs we find that the electron transfer rate J m is highly correlated to V 25 max (J m ∼ 1.73 • V 25 max ; R 2 = 0.91).This finding is supported by Ishida et al. (1999), who expressed that V 25 max depends on the leaf/needle age and is for mature needles about 55.An increase of V 25 max and the combined increase of J m will increase GPP in general.Shirke (2001) also found that the total net photosynthetic rate of young leaves of P. julifora trees was just 36 % of that in mature leaves while old leaves assimilated 76 % of mature leaves.Kattge et al. (2009) investigated the relationship of photosynthetic capacity with leaf nitrogen by assimilating observations of carboxylation capacity V max and maximum photosynthesis rates A max into BETHY coupled with the climate-vegetation model ECHAM5/JSBACH.The optimized V 25 max of tropical trees is substantially lower than earlier estimates currently used in BETHY and the optimized V 25 max of temperate broadleaved trees is substantially higher than the values used in our study.This finding also supports our reduced GPP compared to other studies.
After discussing the uncertainties in measuring NEE using eddy covariance techniques, extrapolating RE and finally deriving GPP and discussing the uncertainties in modelling GPP using BETHY/DLR, it makes sense to embed the validation of our specific model into a wider context of the uncertainties of other model systems in order to elaborate differences or similarities in model behavior.For this reason we investigated several studies.Ichii et al. (2013) compared eight terrestrial biosphere model results with eddy-covariance observations from 24 CarboEastAsia sites from 1901 to 2010.They concluded that "in terms of pattern and magnitude, the carbon fluxes (i.e., gross primary productivity, ecosystem respiration, and net ecosystem exchange) at the temperate and boreal forest sites were simulated best, whereas the simulation results from the tropical forest, cropland, and disturbed sites were poor [. . .].These results indicate that the current model-based estimation of terrestrial carbon budget has large uncertainties." The multi-model analysis of Ichii et al. (2010) using eddy flux observations for validation and calibration came to the conclusion that "models using default model settings showed large deviations in model outputs from observation with large model-by-model variability."As an example, the slope of the regression line for GPP (observed data vs. modelled data) is 0.63 with a coefficient of determination of 0.91.Besides the large underestimation of model results with observations the shift of timing of the start of growing season for the model results are documented.
Huntzinger et al. ( 2012) reported about the North American Carbon Program (NACP) where 19 terrestrial biosphere models and observations at local to continental scales for the period 2000-2005 are evaluated and inter-compared.It was shown that GPP varied between 12.2 and 32.9 PgC a −1 and that "complications in modeling productivity leads to significant disagreement among the model estimates of GPP, with peak growing season differences of greater than 2 PgC month −1 in both Temperate and Boreal NA TransCom regions [. . .], and over 1000 gC m −2 yr −1 in regions of mixed and deciduous broadleaf forests and cultivated and managed lands."An interesting aspect is discussed, "this work suggests that precipitation and radiation data with higher temporal variability yield lower overall GPP [. . .] due to nonlinearities in the photosynthetic functions."Ito (2011) performed a historical meta-analysis of global terrestrial net primary productivity.He surveyed relevant literature from 1862 to 2011 and analyzed 251 estimates of total terrestrial NPP.The mean NPP, the standard deviation and the median were 56.2, 14.3 and 56.4 PgC a −1 , respectively.The relative variability of all estimates is about 25.4 %.When only the literature from 2000 on is investigated the global mean NPP is 59.5 PgC a −1 and the standard deviation is 8.9 PgC a −1 (resulting in a relative variability of about 15 %).The NPP estimates from 2000 on (in total 134) are mainly model-based and partly calibrated using observational data from flux towers.He also analyzed the global mean NPP and the standard deviation derived from different bio-geochemical models (BGC model) (as e.g.BETHY/DLR) and found for a global mean of 58.8 PgC a −1 with a standard deviation of 13.7 PgC a −1 (equivalent with a relative variability of about 23.3 %).From these findings one can conclude that the uncertainty of BGC models to estimate NPP on a global scale is more than 20 %.Several reasons for the uncertainty are discussed, e.g.parameterization of photosynthesis and stomata control, water budget or energy closure.For modelling NPP of forests, he recommended to introduce the stand age as an additional parameter.Also, the difference in input data especially land cover maps and meteorological data are possible error sources.In Wißkirchen (2005) inconsistencies between satellite derived land cover and LAI data sets had been shown to be significant error sources in modelling NEP on a regional scale compared to eddy covariance measurements.Ito and Sasai (2006) showed that modelled NPP is dependent on climate data sets and model parameterization.For two ecosystem models and three climate data sets they revealed global NPP in the range from 39.7 to 63.0 PgC a −1 .Woodward and Lomas (2004) have applied the Sheffield Dynamic Global Vegetation Model (SDGVM) on a global scale with 0.5 • × 0.5 • spatial resolution and daily time step.The calculation of GPP is linear related to LAI.A linear correlation between modelled-and observation-estimated GPP is found with a R 2 of 0.8.However, the slope of the regression line indicates an overall underestimation of about 25 % for the modelling results.It should be noted that Woodward and Lomas (2004) also mention "that estimated errors of NEP measurement by eddy covariance are in the range of ±30 and 180 gC m −2 yr −1 , which could bring many of the differences between simulation and observation within the ranges of measurement error."Baker et al. (2003) simulated and observed CO 2 fluxes at an eddy covariance tower in north-central Wisconsin, USA, using the revised version of the Simple Biosphere model (SiB2.5).This model is highly comparable with BETHY/DLR, using the same the parameterization of stomatal and canopy conductance.In addition, LAI is derived from satellite data (here NDVI) as well as the greenness fraction.SiB2.5 overestimated the magnitude of the net ecosystem exchange of CO 2 in both summer and winter.According to the authors, "mid-day maximum assimilation was well represented by the model, but late afternoon simulations showed excessive carbon uptake due to misrepresentation of within-canopy shading in the model.Interannual variability was not well simulated because only a single year of satellite imagery was used to parameterize the model."

Conclusions
In this study we introduced the BETHY/DLR vegetation model, which is an adapted version of the BETHY scheme by Knorr and Heimann (2001).BETHY/DLR is optimized to be driven with remote sensing derived products to calculate carbon fluxes.For this study we modelled annual gross primary productivity for the period 2000-2007 for Europe.On average we found annual GPP sums of 2.5 PgC a −1 (±0.17 PgC a −1 ).The GPP result for Europe of Ciais et al. (2010) for cropland and of Luyssaert et al. (2010) for forest are in the range from ∼2.7-∼ 3.3 PgC a −1 and in good agreement with our result.
To compare our model results with observation data we used monthly eddy covariance measurements taken from 74 FLUXNET stations distributed all over Europe.The criterion to include a site was data access and at least one consistent year of measurements.In total, data for eight vegetation types comprising 274 consistent years were available.Analysis showed good agreement between most of the main vegetation types.Especially for grassland, mixed forest and coniferous forest sites the model-data comparison shows comparable patterns (grassland: R 2 = 0.83, RMSE = 14 g m −2 month −1 ; mixed forest: R 2 = 0.89, RMSE = 40 g m −2 month −1 ; coniferous forest: R 2 = 0.93, RMSE = 39 g m −2 month −1 ).Intermediate agreement was found for agriculture (agriculture: R 2 = 0.65, RMSE = 101 g m −2 month −1 ) and higher differences for evergreen-broadleaved and deciduous-broadleaved forests (evergreen-broadleaved forests: R 2 = 0.6, RMSE = 54 g m −2 month −1 ; deciduousbroadleaved forests: R 2 = 0.6, RMSE = 75 g m −2 month −1 ).It can thus be concluded that the two approaches result in comparable patterns but with some improvements for evergreen-broadleaved and deciduous-broadleaved forests.Modelled maxima and minima of GPP for evergreenbroadleaved forests exceed the extremes of the observation while the maxima of modelled GPP do not fit the maxima of observed GPP and the length of the vegetative phase is not correctly fitted by model.However, differences in the approaches to calculate GPP have to be considered.These are mainly within the definition of estimated and/or modelled GPP, where the eddy covariance technique includes additional carbon fluxes, which are not considered in BETHY/DLR.Thus, if eddy covariance measurements are considered to be used to validate modelled GPP, it has to be taken into account that quantitative statements on the model accuracy cannot be made.However GPP estimates as provided by FLUXNET can be used to test models on their likelihood to predict e.g.seasonal vegetation patterns within reasonable degrees of uncertainty.Further studies will prove if our findings for Europe are transferable to other regions.In addition, data sets which can be related to biomass increase, such as empirical data on above ground biomass increase (i.e.cereal yields, stem wood increase), are seen as valuable for model validation and should be investigated for their applicability.
On the other hand, a literature survey confirmed our hypothesis that the default parameter of BETHY/DLR for modelling photosynthesis as e.g.V 25 max and J m are too low for temperate forest, due to not yet regarding either the dependence of these parameters on the nitrogen content of leaves/needles or on leaf/needle age.
This study shows that the validation of process-based modelling approaches is restricted by the data availability and comparability of measured data.Thus it is not only necessary to design comprehensive validation and calibration approaches, but also to gain knowledge of the uncertainty and reliability of the data used for comparison.Since complex process models already play an important role for understanding the dynamics of earth systems, and are particularly used to forecast future responses of vegetation to the noticeable climate change (e.g.temperature and water stress), more effort needs to be spent on collecting precise validation data to improve the significance of model results.Further studies should investigate if assimilation of e.g.eddy covariance and remote sensing derived data can further improve the model results.

1Figure 2 :Fig. 2 .
Figure 2: Annual difference map of Gross Primary Productivity for Europe on a 1 km² 2 resolution for the period 2000 -2007, and averaged GPP.Negative deviations are represented 3 in red; positive deviations in blue.The averaged GPP is represented in grey (low values) to 4 green (high values).White areas symbolize urban areas, water bodies and no data.5 6

Fig. 3 .
Fig. 3. Monthly sums of GPP for the period January 2000 to December 2007, for six selected FLUXNET stations.Black courses represent FLUXNET measurements, grey courses BETHY/DLR model results.BL stands for broadleaved.

Figure 4 :
Figure 4: Leaf Area Index time series as derived from CYCLOPES for the IT_Cpz FLUXNET site.Data is given as 10-day composites for the period 2000 -2007.

Fig. 4 .
Fig. 4. Leaf area index time series as derived from CYCLOPES for the IT_Cpz FLUXNET site.Data is given as 10 day composites for the period 2000-2007.

Table 2 .
Additional vegetation types of the BETHY/DLR model with vegetation parameters (after

Table 3 .
74 FLUXNET tower sites and corresponding vegetation types with associated BETHY/DLR translations and measuring times in between 2000 and 2007.

Table 4 .
Annual GPP sums for the period 2000-2007 per country as calculated by BETHY/DLR in MtC a −1 .CV stands for coefficient of variation [%].

Table 5 .
Average GPP values for the 2000-2007 period and corresponding standard deviation for all 74 tower sites in gC m −2 a −1 .