Creating a Regional MODIS Satellite-Driven Net Primary Production Dataset for European Forests

Net primary production (NPP) is an important ecological metric for studying forest ecosystems and their carbon sequestration, for assessing the potential supply of food or timber and quantifying the impacts of climate change on ecosystems. The global MODIS NPP dataset using the MOD17 algorithm provides valuable information for monitoring NPP at 1-km resolution. Since coarse-resolution global climate data are used, the global dataset may contain uncertainties for Europe. We used a 1-km daily gridded European climate data set with the MOD17 algorithm to create the regional NPP dataset MODIS EURO. For evaluation of this new dataset, we compare MODIS EURO with terrestrial driven NPP from analyzing and harmonizing forest inventory data (NFI) from 196,434 plots in 12 European countries as well as the global MODIS NPP dataset for the years 2000 to 2012. Comparing these three NPP datasets, we found that the global MODIS NPP dataset differs from NFI NPP by 26%, while MODIS EURO only differs by 7%. MODIS EURO also agrees with NFI NPP across scales (from continental, regional to country) and gradients (elevation, location, tree age, dominant species, etc.). The agreement is particularly good for elevation, dominant Remote Sens. 2016, 8, 554; doi:10.3390/rs8070554 www.mdpi.com/journal/remotesensing Remote Sens. 2016, 8, 554 2 of 18 species or tree height. This suggests that using improved climate data allows the MOD17 algorithm to provide realistic NPP estimates for Europe. Local discrepancies between MODIS EURO and NFI NPP can be related to differences in stand density due to forest management and the national carbon estimation methods. With this study, we provide a consistent, temporally continuous and spatially explicit productivity dataset for the years 2000 to 2012 on a 1-km resolution, which can be used to assess climate change impacts on ecosystems or the potential biomass supply of the European forests for an increasing bio-based economy. MODIS EURO data are made freely available at ftp://palantir.boku.ac.at/Public/MODIS_EURO.


Introduction
Net primary production (NPP), the difference between Gross Primary Production (GPP) and plant autotrophic respiration, is the net carbon or biomass fixed by vegetation through photosynthesis.NPP represents the allocation rate of photosynthetic products into plant biomass and can be used to measure the quantity of goods provided to society by ecosystems [1][2][3].NPP of forest ecosystems is essential to estimate the potential supply of biomass for bioenergy, fiber and timber supply.NPP is also a key variable to assess environmental change impacts on ecosystems [4] since any variation in the growing conditions influences the carbon cycle due to changes in carbon uptake and/or respiration.As interest grows in utilizing forests for a "bio-based economy" [5,6], more accurate and realistic forest productivity estimates become increasingly important.In addition, competing forest ecosystem services, such as biodiversity or and nature conservation, need to be considered to ensure sustainable use of our forests and to avoid unsustainable over-exploitation of renewable resources.
Within the EU-28 160.9 million ha or 37.9% of the total land area are covered with forests [7].These forests provide resources for the timber industry, the energy sector (24.3% of the energy in the EU-28 is generated from renewable sources of which 64.2% consists of forest biomass and waste [8]), but also for non-timber ecosystem services such as clean air, water, biodiversity or protection against natural hazards.Accurate and consistent forest information is a precondition for assessing the production and harvesting potential of forest resources in Europe.
There are conceptually different data sources and methods to assess forest productivity like: (i) The MODIS algorithm MOD17 uses remotely sensed satellite-data and climate data to predict spatially and temporally continuous NPP and GPP (Gross Primary Production or carbon assimilation) based on an ecophysiological modelling approach [2].In addition to satellite reflectance data and climate data, it requires the biophysical properties of land cover types, which are stored in the Biome Property Look-Up Tables (BPLUT) [9].(ii) National forest inventory data can be used to assess the timber volume stocks as well as volume increment and removal, if repeated observations are available [10].This terrestrial bottom-up approach collects forest information by measuring sample plots arranged on a systematic grid design across larger areas.In combination with biomass expansion factors or biomass functions, volume or tree information can be converted into biomass or carbon estimates to account for differences in wood densities, the carbon fraction and different allocation into compartments [11,12].(iii) Flux towers record the gas-exchange in plant-atmosphere interactions [13], which can be used to derive GPP from Net Ecosystem exchange (NEE).NEE is estimated using eddy covariance data, climate measurements and other ancillary data [14].
Net Primary Production (NPP) from (i) top-down satellite-driven MOD17 algorithm and (ii) bottom-up NPP estimates using terrestrial forest inventory data were compared in a pilot study for Austria on national scale [15].Top-down and bottom-up refer to the level of scaling of the primary recorded information (for MOD17 1-km remote sensing products and for Terrestrial NPP single tree observations).Our definition for top-down differs from traditional carbon cycle modelling [16].This study wants to extend and test this concept for Europe on a continental scale.
For this purpose, we obtain two wall-to-wall spatially-explicit and consistent MODIS NPP datasets by acquiring the global dataset using global climate driver and by creating a regional dataset MODIS EURO using 1-km European climate data.We evaluate these two datasets by comparing with the NPP derived from forest inventory data from 12 European countries.We assess the reliability and potential discrepancies of the MODIS satellite-driven top-down versus the terrestrial bottom-up NPP estimates from continental to national scale and across different gradients like location, elevation or stand density.This will provide a better understanding of the reliability of remote sensing based NPP estimates, which could be used also for regions, where no terrestrial measurements are available.

Materials and Methods
We used two conceptually different methods to estimate NPP, (i) the MODIS NPP algorithm MOD17 and (ii) terrestrial forest inventory data and tree carbon estimation methods.Both have their respective strengths and weaknesses.MODIS NPP has the advantage of providing spatially continuous estimates with a consistent methodology, which is important for any large-scale studies.It incorporates biogeochemical principles in mechanistic modelling environment and the vegetation feedback to climate conditions through changes in Leaf Area Index and absorbed radiation [17].It does not distinguish between different vegetation apart from general Land Cover types, has a coarse spatial resolution and might not be able to represent specific local conditions due to its calibration to global conditions.In contrast, terrestrial forest inventory NPP assesses the actual carbon allocation by trees and captures local small-scale effects (e.g., site conditions, tree age or forest management) as well as regional differences in estimating tree carbon [12,18].It covers only the increment of trees assessed by the inventory system and might not capture local specifics of litter fall and fine root turnover very well, since broad model assumptions have to be used.

MODIS NPP
Since the year 2000, the MOD17 product provides spatially and temporally continuous NPP estimates across the globe [17].The algorithm behind uses the reflectance data from the sensor MODIS (MODerate resolution Imaging Spectroradiometer) of the TERRA and AQUA satellites operated by National Aeronautics and Space Administration of the United States (NASA).MOD17 provides GPP and NPP estimates at a 1-km resolution [2,17] and incorporates basic biogeochemical principles adopted from Biome-BGC [19].It integrates a light use efficiency logic using remotely sensed vegetation information to estimate GPP (Equation ( 1)) with a maintenance and growth respiration module to derive NPP (Equation ( 2)).
GPP " LUE max ˆfTmin ˆfvpd ˆ0.45 ˆSW rad ˆFPAR NPP " GPP ´RM ´RG LUEmax is the maximum light use efficiency, which get adjusted by f Tmin and f vpd to address water stress due to low temperature (Tmin) and vapor pressure deficit (VPD).SWrad is short wave solar radiation load, of which 45% is photosynthetically active.FPAR is the fraction of absorbed photosynthetic active radiation.R M is the maintenance respiration and is estimated using LAI (Leaf Area Index), climate data and biome-specific parameters.R G is the growth respiration and is estimated to be approx.25% of NPP.The complete algorithm is documented in [18] and more details are found in the cited literature therein.
The MOD17 algorithm requires climate data, FPAR and LAI (leaf area index) data as well as land cover data, which is derived from MODIS reflectance data [20].We obtained the global MODIS NPP product (MOD17A3 Version 055) provided by the Numerical Terradynamic Simulation Group (NTSG) at University of Montana available at ftp://ftp.ntsg.umt.edu/pub/MODIS/NTSG_Products/.This data set (hereafter called MODIS GLOB) covers the period of 2000 to 2012, which is the time period covered by our terrestrial data (see next chapter), and provides the annual NPP in gC¨m ´2¨year ´1.
The source of FPAR and LAI input is MODIS15 LAI/FPAR Collection 5, which was temporally gap filled to close data gaps due to unfavorable atmospheric conditions such as cloudiness or heavy aerosol presence [9].For Land cover, we used the land cover product MOD12Q1 Version 4 Type 2 [21] representing the conditions in year 2001.
Climate data are important input into the MODIS NPP algorithm and climate data have a strong impact on the MODIS NPP results [15,22].MODIS GLOB uses the global climate data set NCEP2 [23] described in the following Section 2.2.In Europe, we have high quality daily climate data, the E-OBS data set [24], which was recently downscaled to a 1-km resolution [25].
We next ran the MOD17 algorithm with the downscaled European climate data [25] and obtained an additional MODIS NPP estimate for the period 2000-2012 (hereafter called MODIS EURO), which differ from MODIS GLOB provided by NTSG only in the used daily climate input data.We used the same FPAR, LAI and Land cover input, as used for the global NPP product, MODIS GLOB.MODIS EURO covers our study region, the EU-28 including Norway, Switzerland and the Balkan states (see Figure 1) and is made available under ftp://palantir.boku.ac.at/Public/MODIS_EURO.The source of FPAR and LAI input is MODIS15 LAI/FPAR Collection 5, which was temporally gap filled to close data gaps due to unfavorable atmospheric conditions such as cloudiness or heavy aerosol presence [9].For Land cover, we used the land cover product MOD12Q1 Version 4 Type 2 [21] representing the conditions in year 2001.
Climate data are important input into the MODIS NPP algorithm and climate data have a strong impact on the MODIS NPP results [15,22].MODIS GLOB uses the global climate data set NCEP2 [23] described in the following Section 2.2.In Europe, we have high quality daily climate data, the E-OBS data set [24], which was recently downscaled to a 1-km resolution [25].
We next ran the MOD17 algorithm with the downscaled European climate data [25] and obtained an additional MODIS NPP estimate for the period 2000-2012 (hereafter called MODIS EURO), which differ from MODIS GLOB provided by NTSG only in the used daily climate input data.We used the same FPAR, LAI and Land cover input, as used for the global NPP product, MODIS GLOB.MODIS EURO covers our study region, the EU-28 including Norway, Switzerland and the Balkan states (see Figure 1) and is made available under ftp://palantir.boku.ac.at/Public/MODIS_EURO.

Climate Data
As outlined, the two MODIS NPP estimates, MODIS GLOB and MODIS EURO, differ only in the daily climate data input: MODIS GLOB employs the global NCEP2 climate data set [23] and MODIS EURO uses European downscaled climate data [25].We provide here a brief overview of the two climate data sets.
The NCEP2 data set (NCEP-DOE Reanalysis 2) is a reanalyzed global daily climate data set with a spatial resolution of 1.875° × 1.875°.This corresponds to approx.220 km at the equator at latitude 0° (approx.136 × 220 km at latitude 50°).To compensate the coarse spatial resolution, for MODIS GLOB the climate data for the 1 km MODIS pixels was deduced with an bilateral interpolation method based on the neighboring NCEP2 pixels [9].

Climate Data
As outlined, the two MODIS NPP estimates, MODIS GLOB and MODIS EURO, differ only in the daily climate data input: MODIS GLOB employs the global NCEP2 climate data set [23] and MODIS EURO uses European downscaled climate data [25].We provide here a brief overview of the two climate data sets.
The NCEP2 data set (NCEP-DOE Reanalysis 2) is a reanalyzed global daily climate data set with a spatial resolution of 1.875 ˝ˆ1.875 ˝.This corresponds to approx.220 km at the equator at latitude 0 (approx. 136 ˆ220 km at latitude 50 ˝).To compensate the coarse spatial resolution, for MODIS GLOB the climate data for the 1 km MODIS pixels was deduced with an bilateral interpolation method based on the neighboring NCEP2 pixels [9].
The downscaled climate data used for MODIS EURO provide daily climate data on a 0.0083 ˝ˆ0.0083 ˝resolution (approx. 1 ˆ1 km at the equator and approx.0.6 ˆ1 km at 50 ˝latitude) [25].This data set was developed out of the E-OBS gridded climate data set (0.25 ˝resolution, using data from 7852 climate stations) [24] in conjunction with the WorldClim data set [26].

Terrestrial NFI NPP
Terrestrial forest data such as national forest inventory (NFI) data assess accumulated carbon on a systematic grid using a permanent plot design.From repeated observations of diameter at breast height (DBH) and/or tree height (H) in combination with biomass functions or biomass expansion factors the carbon accumulation of trees is estimated.Since this method is based on single tree measurements and local biomass studies, NPP derived from forest inventory data incorporates local effects such as weather patterns, climate anomalies, stand age, differences in biomass allocation, site and soil effects and different forest densities due to forest management [15,27].
We obtained 196,434 forest inventory plots covering 12 European countries.In Europe, each country has its own National Forest Inventory (NFI) system, which all have different measurement periods, sampling designs and methodologies [10] (Table S1 in the Supplementary Material).Thus, we first had to develop a harmonized and consistent terrestrial dataset for estimating Terrestrial NPP.We calculated NPP using the forest inventory data according to Equation (3).
CARB INC is the carbon increment of trees (gC¨m ´2¨year ´1).FR TO is the carbon used for fine root turnover [28,29].Fine root turnover FR TO is assumed to be equal to the carbon flow into litter C LF [27,30].Both processes are controlled by the same factors and the assumption of similarity between the above-and belowground turnover of short-living plant organs is supported by recently collected European data on fine root turnover [29] and litter fall [31].C LF is the flow of carbon into litter (gC¨m ´2¨year ´1) estimated using a climate-sensitive and species-dependent model [31] and is calculated as: CF is the carbon fraction of dry biomass which is set equal to 0.5 [11].T is the mean annual temperature from the year 2000 to 2012 ( ˝C).P is the mean annual precipitation 2000 to 2012 [mm].For temperature and precipitation we use the European climate data [25] to capture important small-scale regional effects such as elevation or topography in a more realistic way.Equation ( 4) is applied for all plots where broadleaf species contribute most to total basal area and Equation ( 5) is used for coniferous-dominated plots (see Table S2 of the Supplementary Material).
We used data from nine National Forest Inventories (Austria, Czech Republic, Germany, France, Finland, Norway, Poland, Romania, Spain), and three Regional Forest Inventories (Belgium, Estonia, Italy).We grouped our 12 countries in four geographic regions, North Europe, Central-West Europe, Central-East Europe and South Europe [7], to address the large environmental, elevational and climatic gradients in Europe.Countries within a region should have similar climatic and edaphic conditions as well as similar tree allometries and allocation patterns [32].The original locations of the inventory plots were falsified to the nearest pixel of the MODIS grid to guarantee the locations of the plots remain unknown.Temporal consistency with the MODIS data (available since year 2000) was ensured by using only inventory data, which provide CARB INC (Equation ( 3)) for the time period 2000 to 2012. Figure 1 shows our study region with the four geographic regions completely covered by MODIS EURO, and the 12 countries, where we have NFI NPP.
Although all our terrestrial forest inventory data assess properties of trees, there are different sampling methods and increment calculation by country in place, which may strongly affect the resulting estimates [33,34].Four different methods to estimate tree carbon increment CARB INC are used in our data: (1) repeated observations of fixed area plots (used in Norway, Poland, Belgium); (2) repeated angle count sampling (for Austria, Germany, Finland); (3) increment cores (France, Romania, Italy); as well as increment predictions from (4) tree growth models (Czech Republic, Estonia, Italy).Tree growth model predictions were used if no increment observations, either from repeated observations or from increment cores, were available.
In the Supplementary Material, we provide all details for our 12 inventory data sets, the local sampling system, the available data and the used increment method (Table S1 in the Supplementary Material).
The tree carbon results for determining carbon increment CARB INC (Equation ( 3)) were estimated using the carbon calculation method applied by the local forest inventory organization and compiled in [32].Local biomass functions and biomass expansion factors were used to derive tree biomass and carbon fractions to convert biomass into carbon.In the Supplementary Material, we provide a detailed description on processing the NFI data, the tree carbon estimates and stand variables to describe the represented forests (e.g., mean age, basal area or stand density index).
Using this methodology, we processed the forest inventory data from the 12 countries (Table S1) and derived harmonized carbon stocks for all inventory plots.The forest inventory data set consists of 196.434 plots, harmonized across 12 European countries.We applied the carbon increment method for each country and calculated NPP by inventory plot (hereafter called NFI NPP) using Equations ( 3)-(5).

Analysis of NPP Results
We thus have three NPP sources: two using the MOD17 algorithm with different daily climate data: (i) MODIS GLOB produced by the Numerical Terradynamic Simulation Group (NTSG) at University of Montana and (ii) MODIS EURO by running the original MOD17 algorithm and the latest BPLUTs parametrized by [9] with downscaled daily climate data from Europe [25] as well as (iii) Terrestrial NFI NPP using forest inventory data from the 12 countries (Table S1) and local carbon estimation methods [32].
We compared the three NPP datasets across Europe, by our 4 regions (Figure 1) and the 12 countries to analyze our results across different spatial scaling.We extracted for each forest inventory plot at the corresponding MODIS cell the average NPP from MODIS GLOB and MODIS EURO for 2000 to 2012.We next computed for all plots the difference between the two MODIS NPP estimates and the Terrestrial NFI NPP (∆NPP GLOB = MODIS GLOB minus NFI NPP and ∆NPP EURO = MODIS EURO minus NFI NPP).
We used each NFI plot separately and did not compute average values for MODIS pixels.This avoided smoothing effects due to different spacing between inventory grid points and the plot clusters used in some countries (Table S1).
To analyze the effect of gradients on the NPP results, we collected potentially meaningful meta-information such as plot location (Longitude and Latitude in WGS1984), Elevation (EU-DEM 30 m resolution), MODIS Land Cover type or forest characteristics (dominant tree species, mean age, stand density, tree height, etc.) and analyzed patterns of ∆NPP GLOB and ∆NPP EURO across these gradients.
Terrestrial and remote sensing NPP estimates exhibited discrepancies in previous research [15,18] and as explanation the authors suggested changes in stand density, which are commonly caused by forest management and disturbances [15,18].Since major parts of the forests in Europe are managed [7] and affected by natural disturbances such as wind damage or forest fire [35], they should have experienced changes in stand density as compared to unmanaged forests.Stand density directly affects terrestrial NPP estimates by its impact on the development of DBH and H of the remaining trees after forest management operations until canopy closure is reached.On the other hand MODIS NPP is based on the "big leaf" concept and assumes a full coverage of forest area.We thus use Stand density index (SDI) [36] in the analysis of our NPP estimates.

Results
NPP estimated using the MOD17 algorithm has the advantage of providing spatial-and temporal-continuous NPP estimates across Europe on a 1-km resolution and Figure 2   Our NFI dataset covers the full elevational and latitudinal range of forest conditions in Europe including different site conditions, tree species, development stages or management practices.For most countries we have more than 5000 inventory plots (exception: Belgium with 512 plots) and in most cases a plot spacing of at least 4 by 4 km (Table S1).This dataset also provides information on forest properties such as tree age, carbon stocks or stand density and Table 2 indicates that these characteristics vary across Europe.Terrestrial NFI NPP is driven by forest information collected by field crews.Thus it provides NPP and the carbon accumulation by forest stands during a certain time period.Table 1 gives a summary of the forest inventory results by country, by region and the whole dataset, with the terrestrial NFI NPP at the right side.Our NFI dataset covers the full elevational and latitudinal range of forest conditions in Europe including different site conditions, tree species, development stages or management practices.For most countries we have more than 5000 inventory plots (exception: Belgium with 512 plots) and in most cases a plot spacing of at least 4 by 4 km (Table S1).This dataset also provides information on forest properties such as tree age, carbon stocks or stand density and Table 2 indicates that these characteristics vary across Europe.

NPP Estimates across Different Scales
Comparing all our three NPP estimates on a European scale allowed us to explore the general behaviour and evaluate the agreement of the two remote sensing driven NPP products, MODIS GLOB and MODIS EURO, with the terrestrial driven NFI NPP estimates (Figure 3).
Re-running the MOD17 algorithm with local climate data reduced the remotely sensed MODIS NPP in terms of median, mean and variation as compared to the global climate driver (Figure 3).NFI NPP is close to MODIS EURO regarding median and mean, but show larger variation.In addition, Figure 3 confirms that our data is clearly right-skewed (NFI NPP in particular).
Zooming in and examining the different NPP estimates by ecoregion and country allowed us to analyze our results on a higher spatial resolution and to assess local effects such as different regional growing conditions, the impact of local biomass allometries or tree species composition [32] as well as the potential effect of different forest management practices in Europe [7].
We provide in Table 2 the median NPP for the three NPP sources (MODIS GLOB, MODIS EURO and NFI NPP) and the differences between MODIS and NFI NPP (∆NPP GLOB and ∆NPP EURO ), both in absolute values in gC¨m ´2¨year ´1 and normalized in relation to NFI NPP (Rel.∆NPPi in %).Results are given in Table 2 for Europe, by country and for the four eco-regions [7].
At the European level, the MODIS GLOB gives an NPP of 680 gC¨m ´2¨year ´1, the MODIS EURO resulted in 577 gC¨m ´2¨year ´1, and the NPP from the NFI data exhibit a value of 539 gC¨m ´2¨year ´1.The differences in NPP (∆NPP GLOB ) using the global dataset MODIS GLOB are larger than ∆NPP EURO using the regional dataset MODIS EURO (+26% vs. +7%).The same pattern is evident across all four regions and most countries.Only for Poland and Germany ∆NPP GLOB is smaller than ∆NPP EURO .∆NPP GLOB is positive for most countries (negative in only 2 countries), while the discrepancy of MODIS EURO is more randomly distributed in Europe and the 4 regions (∆NPP EURO positive in 5 countries and negative for 7 countries).In addition, Table 2 shows that Rel.∆NPP EURO is smaller than 10% for all countries except five (France, Germany, Czech Republic, Poland and Spain).

NPP Estimates across Different Scales
Comparing all our three NPP estimates on a European scale allowed us to explore the general behaviour and evaluate the agreement of the two remote sensing driven NPP products, MODIS GLOB and MODIS EURO, with the terrestrial driven NFI NPP estimates (Figure 3).Re-running the MOD17 algorithm with local climate data reduced the remotely sensed MODIS NPP in terms of median, mean and variation as compared to the global climate driver (Figure 3).NFI NPP is close to MODIS EURO regarding median and mean, but show larger variation.In addition, Figure 3 confirms that our data is clearly right-skewed (NFI NPP in particular).
Zooming in and examining the different NPP estimates by ecoregion and country allowed us to analyze our results on a higher spatial resolution and to assess local effects such as different regional growing conditions, the impact of local biomass allometries or tree species composition [32] as well as the potential effect of different forest management practices in Europe [7].
We provide in Table 2 the median NPP for the three NPP sources (MODIS GLOB, MODIS EURO and NFI NPP) and the differences between MODIS and NFI NPP (∆NPPGLOB and ∆NPPEURO), both in absolute values in gC•m -2 •year -1 and normalized in relation to NFI NPP (Rel.∆NPPi in %).Results are given in Table 2 for Europe, by country and for the four eco-regions [7].
At the European level, the MODIS GLOB gives an NPP of 680 gC•m −2 •year −1 , the MODIS EURO resulted in 577 gC•m −2 •year −1 , and the NPP from the NFI data exhibit a value of 539 gC•m −2 •year −1 .The differences in NPP (∆NPPGLOB) using the global dataset MODIS GLOB are larger than ∆NPPEURO using the regional dataset MODIS EURO (+26% vs. +7%).The same pattern is evident across all four regions and most countries.Only for Poland and Germany ∆NPPGLOB is smaller than ∆NPPEURO.∆NPPGLOB is positive for most countries (negative in only 2 countries), while the discrepancy of MODIS EURO is more randomly distributed in Europe and the 4 regions (∆NPPEURO positive in 5 countries and negative for 7 countries).In addition, Table 2 shows that Rel.∆NPPEURO is smaller than 10% for all countries except five (France, Germany, Czech Republic, Poland and Spain).This suggests that the discrepancy between MODIS EURO and NFI NPP is smaller than for MODIS GLOB and NFI NPP and we wanted to confirm this along the NPP gradient by showing the country medians in Figure 4.  We used in Figure 4 the aggregated NPP of all inventory plots of one country, since the spatial This suggests that the discrepancy between MODIS EURO and NFI NPP is smaller than for MODIS GLOB and NFI NPP and we wanted to confirm this along the NPP gradient by showing the country medians in Figure 4.
Figure 4 provides the results by country of the NPP estimates resulting from the NFI data versus MODIS EURO with an R 2 0.68, a residual standard error (RSE) of 52.0 gC¨m ´2¨year ´1 or 9.7% of median of the NFI NPP.Aside from Germany and Poland MODIS EURO and NFI NPP are similar across the NPP gradient for the analyzed countries.The results for MODIS GLOB in the right corner exhibit consistent overestimation of NFI NPP, smaller agreement (R 2 = 0.59) and larger error (RSE 80.6 gC¨m ´2¨year ´1 equal 15.0% of median NFI NPP).
We used in Figure 4 the aggregated NPP of all inventory plots of one country, since the spatial coverage and thus the error structure of the two NPP sources are very different (one MODIS pixel covering 1 km 2 or 100 ha and the size of an NFI plot ranging from approx.0.01 to 0.2 ha; Table S1).A direct plot-to-pixel comparison is provided in Figure S1 in the Supplementary Material.

NPP across Elevational, Latitudinal and Longitudinal Gradients
From Figures 3 and 4 as well as Table 2, we can see that the top-down MODIS EURO NPP estimates are consistent with the bottom-up terrestrial driven forest inventory NPP estimates at the European, regional and country level.Next, we investigated whether any patterns across gradients between MODIS EURO and NFI NPP may exist.For this purpose, we showed here ∆NPP EURO for selected gradients, Elevation, Latitude and Longitude.We chose these gradients, since they have a strong effect on environmental and climatic conditions such as growing season length or weather patterns, but also on tree allometries and species composition, and are irrespective of country borders.
We aggregated our results into classes to increase the readability and show Figure 5 the results for whole Europe (results on the different regions are available in Figures S2-S5 in the Supplementary Material).Images for additional gradients like tree age, tree height, MODIS land cover and dominant tree species are provided in Figures S6-S9 in the Supplementary Material.patterns, but also on tree allometries and species composition, and are irrespective of country borders.
We aggregated our results into classes to increase the readability and show Figure 5 the results for whole Europe (results on the different regions are available in Figures S2-S5 in the Supplementary Material).Images for additional gradients like tree age, tree height, MODIS land cover and dominant tree species are provided in Figures S6-S9 in the Supplementary Material.Grouping by elevation in Figure 5a does not indicate striking differences and shows, that the agreement between MODIS EURO and NFI NPP is consistent across the elevational gradients.At certain latitude and longitude classes however local discrepancies exist, which may correspond to the findings in Table 2 and Figure 4.

Stand Density Effects
We analyzed ∆NPP EURO (differences in NPP between MODIS EURO versus NFI NPP) by SDI (Stand Density Index [36] calculated with Equation (S10) in the Supplementary Material) for all of Europe (Figure 6).Grouping by elevation in Figure 5a does not indicate striking differences and shows, that the agreement between MODIS EURO and NFI NPP is consistent across the elevational gradients.At certain latitude and longitude classes however local discrepancies exist, which may correspond to the findings in Table 2 and Figure 4.

Stand Density Effects
We analyzed ∆NPPEURO (differences in NPP between MODIS EURO versus NFI NPP) by SDI (Stand Density Index [36] calculated with Equation (S10) in the Supplementary Material) for all of Europe (Figure 6).∆NPP shows in Figure 6 a significant trend by stand density index SDI (using linear regression; R 0.31; ∆NPP = 103.1 − 0.247 × SDI; p < 0.001), which confirms that differences in stand density have an effect in our data from the 12 European countries.MODIS EURO NPP estimates are higher than NFI NPP at low SDI classes, while at intermediate SDI classes no discrepancies are evident (Figure 6).At high SDI classes MODIS EURO are lower than NFI NPP.
We analyzed the effect of SDI for each country, since SDI could be an explanation for the discrepancies visible in Table 2, Figures 4 and 5. Local effects of forest management intensity, disturbances or differences in the local inventory data design and methodology (Table S1) could lead to differences in SDI.We performed similar graphical analysis as shown in Figure 6 for each country and present here as examples two "extreme" countries: (i) France-positive ∆NPP +10%, with MODIS EURO overestimating NFI NPP; and (ii) Germany-negative ∆NPP −16%, where MODIS EURO underestimates NFI NPP.
For France, MODIS EURO and NFI NPP results agree at high stand density and show ∆NPP shows in Figure 6 a significant trend by stand density index SDI (using linear regression; R 0.31; ∆NPP = 103.1 ´0.247 ˆSDI; p < 0.001), which confirms that differences in stand density have an effect in our data from the 12 European countries.MODIS EURO NPP estimates are higher than NFI NPP at low SDI classes, while at intermediate SDI classes no discrepancies are evident (Figure 6).At high SDI classes MODIS EURO are lower than NFI NPP.
We analyzed the effect of SDI for each country, since SDI could be an explanation for the discrepancies visible in Table 2, Figures 4 and 5. Local effects of forest management intensity, disturbances or differences in the local inventory data design and methodology (Table S1) could lead to differences in SDI.We performed similar graphical analysis as shown in Figure 6 for each country and present here as examples two "extreme" countries: (i) France-positive ∆NPP +10%, with MODIS EURO overestimating NFI NPP; and (ii) Germany-negative ∆NPP ´16%, where MODIS EURO underestimates NFI NPP.
For France, MODIS EURO and NFI NPP results agree at high stand density and show discrepancies at low stand density (Figure 7a).Apparently, MODIS EURO does well in capturing the NPP of stands with high densities, but does not agree with NFI NPP from very open stands.The same patterns are also visible for other countries, where MODIS EURO overestimates NFI NPP such as Spain or Czech Republic (not shown).

Discussion
Top-down satellite driven MODIS NPP (Net Primary Production) estimates using local European daily climate data (MODIS EURO) exhibit smaller differences from the bottom-up terrestrial forest inventory NFI NPP estimates (Table 1) than the original MODIS GLOB estimates using global climate data (Figure 3; Table 2).This confirms that the output from the climate sensitive MOD17 algorithm can be substantially improved by using enhanced daily climate data [22] and supports the findings of the pilot study in Austria [15] by extending the focus to a continental scope.The local European daily climate data [25] used for MODIS EURO reduced across scales from continental (Figure 3) to national scale (Figure 4) substantially the differences between NPP using the MOD17 algorithm and terrestrial forest inventory data (Table 2).Both NPP estimates are also consistent across various gradients (elevation, latitude and longitude in Figure 5 and tree age, tree For Germany on the other hand, MODIS EURO and NFI NPP are similar at low stand density classes, but show increasing deviations with increasing stand densities (Figure 7b).We see the same result for other countries as well, where MODIS EURO underestimates NFI NPP such as Poland (not shown).This may be seen as an indication that besides stand density an additional driver might cause discrepancies between MODIS EURO versus terrestrial NFI NPP.

Discussion
Top-down satellite driven MODIS NPP (Net Primary Production) estimates using local European daily climate data (MODIS EURO) exhibit smaller differences from the bottom-up terrestrial forest inventory NFI NPP estimates (Table 1) than the original MODIS GLOB estimates using global climate data (Figure 3; Table 2).This confirms that the output from the climate sensitive MOD17 algorithm can be substantially improved by using enhanced daily climate data [22] and supports the findings of the pilot study in Austria [15] by extending the focus to a continental scope.The local European daily climate data [25] used for MODIS EURO reduced across scales from continental (Figure 3) to national scale (Figure 4) substantially the differences between NPP using the MOD17 algorithm and terrestrial forest inventory data (Table 2).Both NPP estimates are also consistent across various gradients (elevation, latitude and longitude in Figure 5 and tree age, tree height, MODIS Land cover type and dominant species in Figures S6-S9).
In this study we evaluated MODIS EURO in comparison to the global MODIS NPP dataset [17] using our terrestrial NFI NPP.The specific methodologies and differences of our forest inventory data sets (Table S1) and missing information on fine roots and litter fall do not permit a proper validation of NPP.Since the forest inventory data was collected with a different purpose [10], it contains a different error structure due to the small sample plot size and large grid spacing (one or very few plots within a MODIS pixel) as compared to the continuous 1-km MODIS grid.
The large variations and local discrepancies apparent in this study (Figure 3; Table 2) are also reflected in a study on evaluating NPP and GPP (Gross Primary Production) from the MOD17 algorithm for North and South America [37].While the authors reported no general bias in the MODIS NPP product, they found over-as well underestimation especially for certain locations and forest biomes of more than 30%.This study shows that in Europe discrepancies between MODIS EURO and terrestrial NFI NPP exceeds 10% in three out of twelve countries (Table 2).
This study improves the knowledge on explaining discrepancies between remote sensing and terrestrial NPP estimates by highlighting the effect of stand density index (SDI).Forests with stand density of 200 or lower are expected to have gaps, canopy cover below 100% and low competition between trees.Under such conditions the NFI NPP is substantial lower than MODIS NPP (Figures 6 and 7).This can be explained that at low stand density a substantial share of NPP is undetected by the forest inventory system (gaps filled with young trees or shrubs below diameter threshold), while MODIS NPP is able to capture these gaps via leaf area index provided by the satellite [15].Figure S10 in Supplementary Material confirms that the stand density related trend of ∆NPP in Figures 6 and 7 is mainly caused by NFI NPP, which shows a stronger increase with SDI than MODIS NPP.
Since we tested this effect with MODIS GLOB as well, we can conclude that any MODIS productivity estimates irrespective from the used climate input cannot detect such important effects adequately.The relatively large pixel size of 1-km apparently does not allow MODIS NPP to capture small scale patterns such as clear-cuts, thinning operations or disturbance events, while a forest inventory can detect them better.This confirms the findings of the pilot study in Austria [15] and indicates that differences in stand density needs consideration also on the much larger European scale.
MODIS EURO agrees very well with NFI NPP at average stand densities (Figure 6).This could be explained with the calibration of the BPLUT tables used in the MOD17 algorithm [9] using large-scale global terrestrial NPP data [27].The calibration data most likely represents average forest conditions and may not capture very open or very dense forests adequately.The NFI NPP on the one hand represents the conditions of the (small) area covered by an inventory plot, while MODIS NPP provides a smoothed average NPP of a 1-km pixel.A consistent stand density map at 1-km resolution would be needed to test this hypothesis.
But NFI NPP estimates capture not only differences in stand density and forest management, they are also strongly influenced by local tree allometries and local carbon estimation methods [38].For Germany, stand density cannot explain the observed discrepancies satisfactory in Figure 7b.
In fact the results are quite different compared to whole Europe (Figure 6), France (Figure 7a) or our pilot study [15].Germany is planning to modify the currently used tree biomass estimation methodology [39] which is used in this study, for future carbon assessments.Following reanalysis of existing data [40] and collection and analysis of new sample data [41], improved biomass functions were developed for Germany [42].This new updated methodology results in approx.5% lower aboveground biomass estimates.Thus future German NFI NPP estimates will be lower as well, which will most likely reduce the gap between MODIS and NFI NPP observed for this country in this study.This suggests that, interpretation of discrepancies between NPP estimates needs consideration of the tree carbon estimation methods, since they directly affect increment estimates.
However, there might be other potential drivers leading to inconsistencies both in MODIS EURO and NFI NPP, that could be analyzed in future studies.
Concerning NFI NPP, few countries do not consider adequately the contribution of small trees to the NPP of a forest, either by not considering the ingrowth of small trees [33] or a particular large diameter threshold in some countries (Table S1).This could explain, why in Spain and France MODIS EURO is higher than NFI NPP, as we were not able to include ingrowth here and thus the French and the Spanish NFI NPP estimates might not represent the NPP of their forests sufficiently.
The accuracy of the litter fall and fine root estimates for NFI NPP (Equations ( 3)-( 5) need further research as well.The litter fall models used in this study were derived in a meta-analysis using Eurasian litter fall data [31].They have substantial variation in the used input data and might contain potential inaccuracy, when applied in certain regions.In addition, the estimates for litter fall and fine roots are driven by the same climate data than MODIS EURO.Although the specific climate input differs (periodic average climate used in Equations ( 4) and ( 5) for NFI NPP versus daily maximum, minimum temperature and precipitation used in MOD17), it cannot be ruled out yet that the climate source explains the better match of MODIS EURO and NFI NPP.Thus, the performance of the currently used approach and alternative options for instance by using Foliage mass and Leaf longevity [43] needs to be tested using European litter fall data.
Potential errors in the MODIS EURO product could involve wrong classification of forest biomes by MODIS Land cover [44], limitations of the global parameters of the MOD17 algorithm capturing European forest conditions (see discrepancies in NPP for evergreen broadleaf forests in Figure S3), mismatches in LAI and FPAR by region or forest fragmentation [45].

Conclusions
In this study we created a regional Net Primary Production (NPP) dataset by running the MOD17 algorithm with local European climate data on 1-km resolution for the years 2000 to 2012 (MODIS EURO).We additionally obtained the global MODIS NPP product (MODIS GLOB) and evaluated the two MODIS NPP datasets with bottom-up forest inventory driven NPP (NFI NPP).We thus compared two conceptually different methods for assessing forest productivity across Europe, and test whether local climate data enhances the ability of the MOD17 algorithm to capture European forest conditions.
Running the MOD17 algorithm with local daily climate data substantially improves the quality of MODIS satellite-driven NPP across Europe as compared to the global NPP product (MODIS GLOB).Top-down satellite-driven MODIS EURO and bottom-up NFI NPP agree by regions and by countries, across gradients by longitude, latitude and elevation, if potential discrepancies by stand density due to forest management or the used carbon estimation methods are addressed.
This newly created MODIS EURO dataset is a consistent, continuous, spatial and temporal explicit forest productivity measure of the European forest area providing realistic estimates, which compare well with forest inventory information.This is important since reliable wall-to-wall forest productivity estimates are increasingly important for the growing bio-economy or for increasing our knowledge on other forest ecosystem services such as carbon sequestration.
As long as the MODIS program (based on Satellite "Terra" launched in 1999 and "Aqua" in 2002) is operational and local climate data is available, we can obtain reliable large-scale forest productivity measures for European forests.Since the lifetime of the satellites carrying the MODIS sensor is unknown, we strongly suggest the implementation and testing of this concept in the upcoming European satellite technologies such as the Copernicus Programme to ensure consistent and realistic productivity estimates also in the future.

Supplementary Materials:
The following are available online at www.mdpi.com/2072-4292/8/7/554/s1,Table S1: Summary of the properties of the different forest inventory datasets, Table S2: Tree species groups used in this study, description and selected tree species, Figure S1

Figure 1 .
Figure 1.Our study region separated into four regions, countries with forest inventory data for estimating terrestrial National Forest Inventory (NFI) Net Primary Production (NPP) are marked with dots.

Figure 1 .
Figure 1.Our study region separated into four regions, countries with forest inventory data for estimating terrestrial National Forest Inventory (NFI) Net Primary Production (NPP) are marked with dots.
illustrates this by showing MODIS EURO for the years 2000 to 2012.Note that MODIS EURO also covers not-forest land cover types such as crops, shrub-or grassland.Remote Sens. 2016, 8, 554 8 of 18

Figure 3 .
Figure 3.Comparison of MODIS GLOB and MODIS EURO with NFI NPP: The box represent the Median and the 25th and 75th percentile, the diamond give the arithmetic mean, the whiskers extend to 1.5 of the interquartile range, values outside this range are indicated by circles, on the bottom the number of values represented by the boxplots are given.The number of observations is different since climate data is missing for certain pixels to compute MODIS NPP.To enhance the interpretability of the image, NFI NPP results larger 2100 gC•m −2 •year −1 (445 observations) are not shown, but are included in the boxplot.

Figure 3 .
Figure 3.Comparison of MODIS GLOB and MODIS EURO with NFI NPP: The box represent the Median and the 25th and 75th percentile, the diamond give the arithmetic mean, the whiskers extend to 1.5 of the interquartile range, values outside this range are indicated by circles, on the bottom the number of values represented by the boxplots are given.The number of observations is different since climate data is missing for certain pixels to compute MODIS NPP.To enhance the interpretability of the image, NFI NPP results larger 2100 gC¨m ´2¨year ´1 (445 observations) are not shown, but are included in the boxplot.Remote Sens. 2016, 8, x 10 of 18

Figure 4 .
Figure 4. Comparison of MODIS EURO using European climate data and NFI NPP (MODIS GLOB vs. NFI NPP in the subplot in the bottom-right corner), we present median by country, solid line is 1:1 line, dashed line represents the linear trend of the 12 countries, Coefficient of determination R 2 , Residual standard error (RSE) and the trend function are given.

Figure 4
Figure 4 provides the results by country of the NPP estimates resulting from the NFI data versus MODIS EURO with an R 2 0.68, a residual standard error (RSE) of 52.0 gC•m -2 •year −1 or 9.7% of median of the NFI NPP.Aside from Germany and Poland MODIS EURO and NFI NPP are similar across the NPP gradient for the analyzed countries.The results for MODIS GLOB in the right corner exhibit consistent overestimation of NFI NPP, smaller agreement (R 2 = 0.59) and larger error (RSE 80.6 gC•m -2 •year -1 equal 15.0% of median NFI NPP).We used in Figure4the aggregated NPP of all inventory plots of one country, since the spatial

Figure 4 .
Figure 4. Comparison of MODIS EURO using European climate data and NFI NPP (MODIS GLOB vs. NFI NPP in the subplot in the bottom-right corner), we present median by country, solid line is 1:1 line, dashed line represents the linear trend of the 12 countries, Coefficient of determination R 2 , Residual standard error (RSE) and the trend function are given.

Figure 5 .
Figure 5. NPP Difference (∆NPP) MODIS EURO minus NFI NPP by Elevation classes (a), by Latitude (b) and by Longitude (c), properties of illustration analogous to Figure 3, on the top the number of values represented by the boxplots are given.

Figure 5 .
Figure 5. NPP Difference (∆NPP) MODIS EURO minus NFI NPP by Elevation classes (a), by Latitude (b) and by Longitude (c), properties of illustration analogous to Figure 3, on the top the number of values represented by the boxplots are given.

:
Direct pixel-to-plot comparison of MODIS EURO and NFI NPP, Figure S2: For North Europe ∆NPP grouped by Elevation, Latitude and Longitude, Figure S3: For Central-West Europe ∆NPP grouped by Elevation, Latitude and Longitude, Figure S4: For Central-East Europe ∆NPP grouped by Elevation, Latitude and Longitude, Figure S5: For South Europe ∆NPP grouped by Elevation, Latitude and Longitude, Figure S6: Difference ∆NPP grouped by age classes, Figure S7: Difference ∆NPP grouped by tree height classes, Figure S8: Difference ∆NPP grouped by MODIS Land cover types, Figure S9: Difference ∆NPP grouped by dominant species, Figure S10: MODIS EURO and NFI NPP by Stand density Index (SDI) classes.

Table 2 .
NPP and ∆NPP (always using median) for the whole dataset ("All Countries"), for each country separately and for each region (MODIS NPP using global climate data-MODIS GLOB; MODIS NPP using local European climate data-MODIS EURO and NPP using forest inventory data-NFI NPP); ∆NPP and Rel.∆NPP both for MODIS GLOB and MODIS EURO.Positive differences indicate that MODIS NPP overestimates NFI NPP and vice versa.NPP and ∆NPP (gC•m −2 •year −1 ) MODIS MODIS ∆NPP Rel.∆NPP [%] GLOB EURO NFI NPP GLOB EURO GLOB EURO

Table 1 .
[36]ary of the forest inventory results: Number of plots with data, Time period covered by NFI NPP, Mean elevation (range Minimum-Maximum) in meter above sea level (EU-DEM 30 m resolution).For the following plot statistics we provide mean and standard deviation: Mean quadratic DBH (cm), Mean Tree height (m), Basal area at 1.3 m height (m 2 ¨ha ´1), Stem number (ha ´1), Tree carbon per hectare (gC¨m ´2), Median age class, SDI Stand Density Index[36](for details on this variables see Supplementary Material), NPP is the NFI Net primary production (gC¨m ´2¨year ´1) according to Equation (3), For Czech Republic we only have country means.Empty cells (-) indicate that this variable is not available from the NFI data set.At the end of each section, statistics of the region are given and at the bottom of the table summary statistics for whole Europe.

Table 2 .
NPP and ∆NPP (always using median) for the whole dataset ("All Countries"), for each country separately and for each region (MODIS NPP using global climate data-MODIS GLOB; MODIS NPP using local European climate data-MODIS EURO and NPP using forest inventory data-NFI NPP); ∆NPP and Rel.∆NPP both for MODIS GLOB and MODIS EURO.Positive differences indicate that MODIS NPP overestimates NFI NPP and vice versa.