- Download PDF-
SEEFOR 8 (1): 1-20
Article ID: 134
Original scientific paper
Identification of Years with Extreme Vegetation State in Central Europe Based on Remote Sensing and Meteorological Data
Anikó Kern1*, Hrvoje Marjanović2, Laura Dobor3, Mislav Anić2, Tomáš Hlásny3, Zoltán Barcza4,5
(1) Eötvös Loránd University, Department of Geophysics and Space Science, Pázmány P. st. 1/A, H-1117 Budapest, Hungary; (2) Croatian Forest Research Institute, Cvjetno naselje 41, HR-10450 Jastrebarsko, Croatia; (3) Czech University of Life Sciences, Faculty of Forestry and Wood Sciences, Department of Forest Protection and Entomology, Kamýcká 129, 165 21 Prague 6, Czech Republic; (4) Eötvös Loránd University, Department of Meteorology, Pázmány P. st. 1/A, H-1117 Budapest, Hungary; (5) Eötvös Loránd University, Faculty of Sciences, Excellence Center, H-2462 Martonvásár, Brunszvik u. 2.
* Correspondence: e-mail: email@example.com
Citation: KERN A, MARJANOVIĆ H, DOBOR L, ANIĆ M, HLÁSNY T, BARCZA Z 2017 Identification of Years with Extreme Vegetation State in Central Europe Based on Remote Sensing and Meteorological Data. South-east Eur for 8 (1): 1-20. DOI: https://doi.org/10.15177/seefor.17-05
Received: 15 Feb 2017; Revised: 10 Apr 2017; Accepted: 12 Apr 2017; Published online: 25 Apr 2017
Background and Purpose: Determination of an extreme year from the aspect of the vegetation activity using only meteorological data might be ambiguous and not adequate. Furthermore, in some ecosystems, e.g. forests, the response is not instantly visible, but the effects of the meteorological anomaly can be seen in the following year. The aim of the present paper is to select and characterize typical and anomalous years using satellite-based remote sensing data and meteorological observations during the recent years of 2000-2014 for Central Europe, based on the response of the vegetation.
Materials and Methods: In the present study vegetation characteristics were described using remotely sensed official products of the MODerate resolution Imaging Spectroradiometer (MODIS), namely NDVI, EVI, FPAR, LAI, GPP, and NPP, with 8-day temporal and 500 meter spatial resolution for the period of 2000-2014. The corresponding mean temperature and precipitation data (on the same grid) were derived from the Open Database for Climate Change Related Impact Studies in Central Europe (FORESEE) daily meteorological dataset. Land cover specific anomalies of the meteorological and vegetation characteristics were created and averaged on a country-scale, where the distinction between the main land cover types was based on the synergetic use of MODIS land cover and Coordination of Information on the Environment (CORINE) Land Cover 2012 datasets.
Results: It has been demonstrated that the anomaly detection based solely on basic meteorological variables is ambiguous since the strength of the anomaly depends on the selected integration time period. In contrast, the effect-based approach exploiting the available, state-of-the-art remote sensing based vegetation indices is a promising tool for the characterization of the anomalous behaviour of the different land cover types. The selection of extreme years was performed in an explicit way using percentile analysis on pixel level.
Conclusions: Plant status in terms of both positive and negative anomalies shows strong land cover dependency in Central Europe. This is most likely due to the differences in heat and drought resistance of the vegetation, and species composition. The selection of country-specific extreme years can serve as a basis for forthcoming research.
Keywords: remote sensing, anomalous vegetation conditions, phenology, MODIS
Variability is one fundamental feature of the Earth’s climate system [1, 2]. Climate fluctuations result in considerable interannual variability of the meteorological parameters like temperature, precipitation, cloud cover, radiation and others. Climate change research typically calculates meteorological parameters in annual basis to track global changes and patterns . Percentiles as the measures of annual deviations from the long term mean are also used to quantify expected return-time of extreme weather events .
Meteorological parameters exhibit strong intra-annual variability as well with extreme events like heatwaves, drought spells, flash floods, extreme cold periods etc. In meteorology, quantification of the extreme events has long tradition [1, 4, 5]. Climate indices are calculated typically from daily meteorological records to estimate e.g. strength of heatwaves, number of days with extreme precipitation, length of dry periods etc. [1, 6, 7].
Plant processes and productivity have strong economic impacts due to its direct relationship with crop yield, wood production, animal welfare, fodder quantity and quality, ecosystem services, etc. . Thus, quantifying plant production and understanding the processes behind variability of plant status on annual scale is of high importance [9, 10].
Studies focusing on ecosystem processes have problems with the existing climate indices and annual means due to the complex interactions between meteorological parameters and the vegetation state . Intra-annual variability of the meteorological parameters might exert strong impact on plant growth and productivity . Actual vegetation state is the result of beneficial and negative effects from the past. In other words, vegetation status is the integrator of all effects from soil water status, temperature, variability of radiation, etc. Therefore extreme indices and annual means are not directly applicable to study plant growth and explain its interannual variability.
One logical way to study the impact of climate variability on plant processes is the effect-based approach, where the plant response is characterized first, then the meteorological cause of the effect is sought (see e.g. [13, 14]). Due to economic and political reasons these studies should practically focus on the country scale (see e.g. ). However, quantification of plant status on large spatial scales is not straightforward due to the small spatial representativeness and amount of observation based, in situ data (e.g. biomass data, phenology observations, leaf area index measurements, eddy covariance data, crop census data, etc.).
Remote sensing (RS) provides a very convenient solution to this problem. Data obtained from sensors on board Earth Observing satellites like the Moderate resolution Imaging Spectroradiometer (MODIS) provides unique information  which can be linked to plant status via vegetation indices like Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI) , Leaf Area Index (LAI), Fraction of Photosynthetically Active Radiation (FPAR)  or modelled plant production like Gross Primary Production (GPP), and the Net Primary Production (NPP) . They are produced at convenient spatial scale (with finest 250 m – 1 km), and today already comprise an exceptional 17-years-long data record.
The spatial resolution of the MODIS Collection 6 products for NDVI, EVI, LAI, FPAR, GPP, and NPP is 500 m or finer (for NDVI and LAI with 250 m resolution) which is sufficiently small to isolate large number of pixels that contain a single land cover (LC) type. This enables that the analysis of plant behaviour can be studied for each land cover type separately, which is needed due to the substantially different dependence of plant functional types on the environmental conditions . However, the classification of pixels into LC categories is not trivial due to land use change and land cover classification errors. For best results, other sources of land cover type are needed, like national/regional forest management maps, agriculture land use monitoring systems for payment of subsidies (e.g. in EU countries), or CORINE Land Cover database . Here we propose a novel method that would enable selection and identification of the anomalous years based on RS data and RS-derived products (NFVI, EVI, LAI, FPAR, GPP, and NPP) according to different land cover types and spatial coverage. This method overcomes the issues that are related with the study of the meteorological anomalies alone. Based on the area for which the method is implemented, it could help to assess differences at spatial and temporal scale for a given LC within the area of interest. It can be applied at a national level to study the behaviour of each LC type as well as to assess the magnitude of the anomaly that has occurred in terms of “repeat time”.
In the present paper country averages are studied to get robust results and to support future research. At smaller spatial scales climate is not necessarily the dominant driver of plant development and vigour. For example forest production can be adversely affected due to pathogens or insects . Cropland management clearly affects plant status at smaller scales; for example winter wheat is typically harvested in the region by June-July which affects overall greenness. Management practices like fertilization, harvest and irrigation, and also soil type clearly modulate the effect of meteorological conditions in grasslands and croplands. Country averages by land cover types are expected to be driven by climate fluctuations and provide robust and clear signal about the overall productivity in the given country (see e.g. [15, 22] for such approaches).
The aim of the paper is to present a method to select years that can be characterized as anomalous based on observed plant status and greenness. Using multiple vegetation indices we also test the similarity/dissimilarity of the different vegetation metrics in terms of their usability to detect anomalous plant status. To get robust and easily interpretable results we used country-means. Bosnia and Herzegovina, Croatia and Hungary (and in Supplement Czech Republic, Slovakia and Slovenia) were selected to represent Mediterranean, continental and alpine climates. The reason behind such choice is the relatively small geographical extent and the fact that extreme weather events typically affect large areas and sometimes the whole country . Countries with transitional economies are particularly vulnerable to climate extremes . It means that understanding the cause-effect relationships might support the prevention of adverse effects on plant state and this can support economic growth and human welfare.
MATERIALS AND METHODS
The target area of this study is the broader region of the Carpathian Basin, determined by the coverage of the applied meteorological dataset (see below). Bosnia and Herzegovina, Croatia and Hungary were selected to represent Mediterranean and continental climate as well, with high biodiversity and variability in the meteorological conditions. Results for Czech Republic, Slovakia and Slovenia are presented in the Supplemental Material. The reason behind chosen countries is the spatial domain of the applied meteorological database, covering fully only the selected countries.
To investigate the effects of weather variability on the vegetation activity and greenness, we used the FORESEE database . This freely available meteorological database contains observed and projected daily maximum/minimum temperature and precipitation fields for Central Europe on a regular grid with a spatial resolution of 1/6° × 1/6°, covering the years between 1951 and 2100. For the 1951-2014 period, FORESEE provides observation-based interpolated meteorological fields for the wide region of the Carpathian Basin.
In order to match the temporal and spatial resolution of the applied remote sensing datasets and of the FORESEE, 8-day mean temperature values and precipitation sums were calculated on the finer grid of the MODIS products based on the methodology of Kern et al. . Country averaged mean and anomaly values were calculated for the above mentioned six countries within the study area.
Mean temperature and precipitation fields (respectively) calculated for the entire domain area of the FORESEE database are presented in Figure 1 and 2 for the study period of 2000-2014. Yearly anomaly fields, shown in Figure 3-6, are illustrations of years which might be considered as extreme from the meteorological perspective. Figure 3 and 4 show the anomaly maps of temperature and precipitation (respectively) for 2011, while Figure 5 and 6 shows the anomaly maps of temperature and precipitation (respectively) for 2014 (note that the temperature and precipitation anomaly maps have different legends).
Vegetation Related MODIS Products
In the present study we used the latest version (Collection 6) MODIS NDVI, EVI, FPAR, LAI, GPP, and NPP with 500 × 500 m spatial resolution, as part of the MOD13A1, MOD15A2H, MOD17A2H and MOD17A3H official products [26, 27] derived from the measurements of MODIS sensor on board satellite Terra. The longest possible datasets (covering the period of 2000-2014) were chosen to match the availability of the MODIS products and the temporal coverage of the applied meteorological dataset.
Quality filtering of each dataset was performed using the quality flag information included in the datasets, based on the method described by Kern et al. . Besides this, to filter out unrealistic sudden increases and decreases in the state of the vegetation the so called Best Index Slope Extraction (BISE) method  was applied afterwards on pixel level for NDVI and EVI data . While FPAR, LAI and GPP has the same 8-day temporal resolution, NDVI and EVI are 16-day composite products, therefore temporal resampling based on the Julian date information included in the MOD13 datasets was also necessary to create NDVI and EVI dataset with the same 8-day temporal resolution .
It has to be noted that although most of the MODIS products are state-of-the-art, standardized, well documented data sets, the quality of the model based products is expected to be lower and that affects their applicability . Most notably, MODIS annual NPP products suffers from a major error (with large unexpected positive bias for the first four years (2000-2003) as it is presented in Figure 7 and should not be used in the FORESEE domain. Therefore, the anomaly values from the yearly NPP products were calculated relative to the period of 2004-2014.
FIGURE 7. Country-averaged annual mean MODIS NPP values for the investigated countries in Central Europe (note that due to technical issues the first 4 years (2000-2003) are not applicable for further processing).
Methods and Metrics for Defining Anomalous Years
From the quality (and in the case of MOD13A1 data BISE-filtered) 8-day temporal resolution datasets country averaged multiannual means and yearly anomaly values were calculated for the whole vegetated area of the selected countries within the Carpathian Basin. Land cover specific country averaged values were also derived for broadleaf and coniferous forests, croplands and grasslands. The yearly anomalies were calculated without the first 40 and last 40 days of the year to avoid any misleading result originating from the effect of snow cover. Using the derived yearly anomaly values, relative anomaly values were calculated by dividing the yearly anomaly values by the maximum of the absolute anomaly values during the investigated 2000-2014 period.
To analyse the spatial distribution of the yearly anomaly fields, we calculated the 0.5, 2, 9, 25, 50, 75, 91, 98 and 99.5 percentiles values of all yearly mean anomaly values for all land cover specific pixels separately within a given country. The selected percentiles define thresholds for the classification of anomalies into the following 10 categories:
The derived percentile values are appropriate to des-cribe the total distribution of the anomaly values of a given land cover type (within a given country) during the study period of 2000-2014. Based on the proposed categorization we calculated for every year separately the number of the pixels (as the percentage for a given LC type) within all the created percentile ranges. Note that there might be situations when e.g. strong positive anomaly was not occurring in the entire time period for a given LC type at all. In this study the maximum anomaly is used independent of its magnitude.
The advantages of the proposed categorization is the following: (1) the method we applied is the same for all LC types, but it yields LC specific results; (2) it enables the identification of areas within the selected geographical area (e.g. country), where mean anomaly is more positive or more negative with respect to the country average, indicating that some areas are more/less productive, more/less prone or more/less sensitive to meteorological anomalies (where a special attention should be paid to croplands due to the possible yearly crop type change); (3) it is appropriate to be done for any other vegetation related characteristics.
In this paper we propose the following classification of the years, which is based on the percentage of the area showing positive and negative anomalies of a given vegetation characteristics (NDVI, EVI, FPAR, LAI, GPP, and NPP). Using the yearly anomaly fields of given vegetation characteristic, for a specific land cover type, the chosen year can be graded based on the thresholds presented in Table 1. It should be noted that the grade for a given year at a given location depends on the selected reference period (e.g. MODIS era, in our case 2000-2014) and the area representing a domain (e.g. region, country or a continent).
Synergy of the MODIS and CORINE Land Cover Databases
In order to study the response of the various vegetation types to meteorological anomalies we distinguished the main land cover types based on the synergistic use of two land cover datasets. The so called MCD12 land cover products based on MODIS observations  and the CORINE 2012 database  were used to identify pixels with high probability of being broadleaf forests, coniferous forests, croplands and grasslands.
From the five types of land cover classification contained in the MCD12Q1 product  we used the Type-1 (International Geosphere Biosphere Programme - IGBP) to identify croplands and grasslands and the Type-3 (MODIS-derived LAI/fPAR scheme) classifications for broadleaf and coniferous forests. The categories of the land cover classification for the used MODIS pixels are given in the translation matrix between the different LC schemes (Table 2). The reason behind the usage of Type-3 for forests was that the widely used Type-1 IGBP land cover classification suffers from well-documented errors [32, 33] with 75% overall accuracy  including misclassification of the mixed forest types as well . We found that the misclassification is especially evident for the lowland forests along rivers, for example like the Gemenc forest by the Danube in Hungary, or the river Spačva basin in eastern Croatia, extending also to Vojvodina in Serbia. According to MODIS Type-1 land cover classification both forests have been categorized as mixed forest, while they are broadleaf forests in reality. In fact Spačva forest with its area of 40 kHa size is the largest complex of pedunculate oak forest in Europe . In contrast to this, the Type-3 classification doesn’t have mixed forests category, resulting in all forested pixels categorized as either broadleaf or coniferous category. We have to note that in Type‑3 classification croplands and grasslands are categorized together as “Grasses/Cereal crops”, therefore Type-3 classification cannot be used to discriminate croplands and grasslands.
From the yearly MCD12Q1 land cover datasets with 500 m × 500 m spatial resolution covering the period of 2001-2013 we selected the pixels (hereafter called stable LC pixels) which had no land cover change during the entire period. The usage of the temporal stability selection criteria excluded pixels which underwent any kind of land cover change during the investigated period, either actual, or resulting from the error in classification. Implementation of such a strict rule resulted with dropping significant number of the pixels. The percentage of stable LC pixels was only ~56% of the whole study area based on Type-1 classification. Specifically, 32%, 34%, 53%, 69%, and 28% of the pixels remained for the categories of broadleaf forests, coniferous forests, mixed forests, croplands, and grasslands, respectively. Based on Type-3 classification the percentages of stable LC pixels in the whole study were 48% and 42% for broadleaf and coniferous forests, respectively. The percentages were calculated relatively to the mean number of the pixels during the 13 years of the total dataset.
In order to increase the reliability of the applied LC categorization we used the CORINE land cover dataset as a reference LC dataset for the year 2012. The accuracy of CORINE is 87.0 ± 0.7% , which is significantly better than the 75% accuracy of MODIS . Using GIS software [QGIS 2.16], we intersected the vector-based CORINE layer with the grid of the MODIS pixels at 500 m × 500 m spatial resolution and obtained the share of every of the 44 CORINE land cover types present within each of the MODIS grid-cells (i.e. pixel). Due to differences in land cover categorization between MODIS and CORINE, we used translation matrix (Table 2) to unify the classification. Based on the area share information of CORINE we selected so called CLC2012 pure pixels for deciduous forests, coniferous forests, croplands and grasslands, which contained at least 90% area share from the given LC type (Table 2, column 1). In the case of croplands and grasslands the 90% threshold was applied for the sum of the area shares of the different CLC categories (separately for croplands and grasslands) listed in Table 2 (column 3).
Using the derived set of MODIS stable pixels and CORINE pure pixels we selected “reliable” pixels which were classified by both land cover database to the same vegetation type of broadleaf forests, coniferous forests, croplands and grasslands. Figure 8 shows the location of the remaining reliable pixels of the main land cover types based on MODIS IGBP (MCD12 Type-1), LAI/FPAR scheme (MCD12 Type-3) classifications and CLC2012 database. This land cover map shows the location of the pixels which were finally used in the present study as reliable pixels, being constant during the study period and having probably right classification. Based on the information of the CORINE database regarding the percentage values of the presence of a given land cover type the mean percentage for the selected broadleaf forests, coniferous forests, croplands and grasslands pixel was 99.0%, 99.0%, 99.3% and 98.2%, respectively. The problem of using the MODIS Type-1 classification for forests is illustrated in Figure 9. The presented map shows the location of the pixels which are classified as broadleaf forests and coniferous based on the MODIS Type-3 and CORINE classifications, but not based on MODIS Type-1 classification. The number of the forested pixels (in the study area) which are not represented well in Type-1 classification is considerable: 59.5% and 66.3% for broadleaf and coniferous forests, respectively, relatively to the number of the derived reliable broadleaf and coniferous pixels, which were selected finally to our further investigations.
FIGURE 9. Location of the broadleaf and coniferous forests pixels according to both the MODIS Type-3 and CORINE classifications which are classified by MODIS Type-1 classification in some other LC category.
After the synergetic use of the two MODIS datasets with CORINE dataset the number of the selected pixels was 6.3% for broadleaf, 3.1% for coniferous forest, 17.1% for cropland and 0.1% for grassland pixels relative to the all vegetated pixels in the whole domain. It means, with the final pixel selection we were taking into account 26.6% of the whole domain in a total value. The country specific pixel numbers are presented in Table 3.
TABLE 3. Number of the pixels at 500 m × 500 m spatial resolution which complied with the selection criteria separately to the main LC types and countries (selection criteria: MODIS LC type for a given pixel did not change during 2001-2013 and its LC classification with 90% share corresponds to that of CORINE 2012 - see Table 2).
RESULTS AND DISCUSSION
Problems with the Selection of Extreme Years Based Solely on Meteorological Data
In order to study the dependence of the strength of meteorological anomalies on the selected integration time period, first we present country averaged mean temperature and precipitation anomalies in Figure 10 (left and right side images, respectively) during 2000-2014 for Bosnia and Herzegovina, Croatia and Hungary (and in Figure S1 in the Supplementary Material for Czech Republic, Slovakia and Slovenia).
While the columns show whole year anomalies in Figure 10, the curves indicate periods with various length (starting at different date and all ending with 7th of October within the year), representing gradually only the predominant part of the growing season. 7th of October was used as the end of the integration period as we can hypothesize that meteorology in the dormant season affects plant state to a lesser extent. The beginning was varied with ~2 weeks periods and ranging from 1st January to 26th June.
FIGURE 10. Temperature (left side images) and precipitation (right side images) anomalies of time periods with different length within the years during 2000-2014 based on the FORESEE database for Bosnia and Herzegovina, Croatia and Hungary (curves represent the selected time periods ending in 7th October, while columns indicate anomalies for the entire year).
The figure shows that even if there is an extreme temperature or precipitation anomaly for a given year, it does not necessarily mean that the vegetation was exposed to similar meteorological anomalies during the growing season. A good example is 2003, when the country-averaged annual temperature anomaly was negative in the presented countries (especially for Hungary), but a strong, mostly positive anomaly was detectable for the growing season. On the contrary, year 2004 or 2005 (with similar negative annual temperature anomalies) showed only negative temperature anomalies during the shorter time periods until the end of the growing season. In terms of annual temperature anomaly year 2007 was very similar to year 2000 for all three selected countries, but with much higher intra-annual variability, on average having cooler than usual period in summer and early autumn. On the contrary, in 2012 (with annual temperature anomaly similar to the one in 2007) the summer-early autumn period was characterized with the highest positive anomaly during the year, when the mean temperature of the countries was continuously higher than the average. Finally, year 2014 is worth mentioning as well. Though it had the largest annual anomaly, its temporal evolution was not consistent. Precipitation showed similar features, however the intra-annual variability was much lower. It is important to mention the two consecutive years of 2010 and 2011, which was characterized by very strong positive (2010) and then negative precipitation anomaly (2011) all over the Carpathian-Basin.
These results illustrate the problems in defining temperature and precipitation anomalies as they strongly depend on the selected integration period. In other words, we cannot unambiguously select extreme years based on the yearly (or monthly/seasonal) mean meteorological conditions that can be related with e.g. vegetation state, crop yield, forest productivity, outbreak of insects or other phenomena. To find the extreme years in plant greenness and describe the response of the vegetation greenness, it is more straightforward to study anomalous behaviour of the plant state (or other phenomena that is related with ongoing climate anomalies).
Selection of Anomalous Years Based on the Overall Impact on Vegetation
Yearly relative anomalies of the vegetation related characteristics (such as NDVI, EVI, FPAR, LAI, GPP, and NPP) during 2000-2014 are presented in Figure 11 for croplands and grasslands. Figure 12 shows the same for broadleaf and coniferous forests for Bosnia and Herzegovina, Croatia and Hungary (Figure S2 and S3 in the Supplementary Material illustrate results for the Czech Republic, Slovakia and Slovenia). The yearly, land cover specific anomalies of the meteorological variables (temperature and precipitation) are also presented in Figure 11 and 12.
FIGURE 11. Yearly relative anomalies of the vegetation related characteristics (such as NDVI, EVI, FPAR, LAI, GPP and NPP) and of the meteorological variables (temperature and precipitation) during 2000-2014 for croplands and grasslands of Bosnia and Herzegovina, Croatia and Hungary.
FIGURE 12. Yearly relative anomalies of the vegetation related characteristics (such as NDVI, EVI, FPAR, LAI, GPP and NPP) and of the meteorological variables (temperature and precipitation) during 2000-2014 for broadleaf and coniferous forests of Bosnia and Herzegovina, Croatia and Hungary.
Based on the results, the relative anomalies for the different vegetation characteristics could be described as quite consistent for croplands and grasslands, showing the same direction of relative anomaly in most of the cases. It is also notable that GPP and NPP show opposite character for some years compared to other indexes. (Anomalies of NPP for the year 2000-2003 are not shown; see section 2.4. in the Materials and Methods for the explanation).
For a given land cover type in some area (country) within the reference period 2000-2014, the calculated relative anomalies facilitate straightforward comparison between years, as well as between different land cover types within the same year. For example, we can select years when croplands and grasslands were both anomalous (like 2000, 2003 and 2012 with negative anomalies, or 2010 and 2014 with positive anomalies), but in some years (such as 2001 and 2004) croplands and grasslands showed different response to the environmental conditions. These differences might be related to land cover type specific features, species distribution and human disturbances (harvest, fertilization), and maybe other factors like soil water holding capacity, etc. Additionally, croplands in some countries cover C4 plants with higher drought resistance, which might affect the results. Also, croplands and grasslands might exhibit different critical time periods when the ongoing weather affects their productivity . Clearly, it is not only the yearly anomalies of the meteorological parameters that affect the plant status, but the within-year temporal pattern of the meteorological elements. Note that the number of pixels (see Table 3) was low for grasslands in all countries except Bosnia and Herzegovina and Croatia, which can affect the reliability of the results for grasslands in those countries.
While in the case of croplands and grasslands the vegetation related characteristics are mostly consistent, forests reveal different behaviour (Figure 12). As it is obvious from the figure, in the case of forests the different characteristics are not showing similarity for the anomalies, as it was found for croplands and grasslands. It is especially true for GPP and NPP, which are in fact models and not indices that are derived directly from the reflectance data . It means that GPP and NPP results must be treated with caution. For example, in 2012 the reason behind the large disagreement between the mean anomaly of NDVI, EVI, LAI and FPAR, and the mean behaviour of GPP and NPP might be related with the extreme meteorological conditions during the growing season, causing significant drought in the Carpathian-Basin and a consequently simulated lower productivity by the GPP-model. In addition, forests with typically deeper root zone do not react in the same way to changing meteorological conditions as the shallow-rooted croplands and grasslands . Forest productivity depends on the available soil water content in the topsoil to a lesser extent, because deeper groundwater supplies can be accessible to trees. Trees store a relatively large amount of non-structural carbohydrate – which is the product of the previous year – to mitigate negative effect of shortage in nutrients and photosynthetic carbon uptake. Therefore they are likely to be less affected by the shorter term weather anomalies, but showing stronger exposure to the longer term changes. Based on the observed results, we could not select the same years as anomalous years for herbaceous vegetation, which were highlighted in the case of croplands and grassland over the studied countries. The role of the carbohydrate reserves seems to be important which might impose lagged effect . This can be recognised in 2013 as the effect of the previous, extreme year (evident in strong negative anomaly for croplands, see Figure 11) which has contributed to a negative anomaly in broadleaf forests (see Figure 12).
Co-variation between the meteorological and vegetation anomalies does exist and it can give us some basic information at the annual scale. However, the main aim of the present study was not to execute in-depth analysis of cause-effect of environmental variables on the detectable plant state anomalies, but to select anomalous years. In order to quantify their co-variance Pearson’s r values (correlation coefficients) were calculated and compared. The correlation coefficients between the yearly anomalies of the vegetation characteristics and the yearly meteorological anomalies were calculated for the different plant types separately. Table 4 present the correlation coefficients calculated between the NDVI anomalies and the temperature and precipitation anomalies, where the statistically significant (P<0.05) values are indicated. In the case of temperature, the calculated r values show the strongest correlation for forests (especially with broadleaf forests), while in the case of precipitation the largest correlation is present for croplands and grasslands. Both the direction and the strength of the correlation depend strongly on the land cover type, indicating complex relationship between the vegetation and meteorological conditions. It also implies that large scale and land cover specific studies should not be made together, but separately and at finer temporal and spatial scale.
TABLE 4. Correlation coefficients (r) between the yearly NDVI and meteorological (temperature and precipitation) anomalies for the main land cover type in the case of the different countries (note, for Czech Republic and Slovenia the number of the selected grassland pixels were less than 5 (Table 2)).
Selection of Anomalous Years Based on the Magnitude and Spatial Extent of NDVI Anomalies
The investigation of the country averaged mean yearly anomaly does not provide information about its spatial distribution. Therefore, pixel-level investigations were performed using percentiles analysis for the area of a given land cover type affected with severe anomaly (see Section 2.4 in the Materials and Methods). Here we present the results based only on the NDVI, while the results for all other metrics are provided in the Supplementary material (S6, S8, S10, S12).
Different land cover types react with varying intensity to meteorological anomalies as it can be inferred by careful inspection of Figure 13. It shows that, for example, broadleaf forests (top left) maintain rather stable mean annual NDVI and the corresponding anomalies are in the range of +/- 0.05 of the mean in 99% of the cases. Coniferous forests (Figure 13, top right), on the other hand, show higher variability of the mean NDVI compared to broadleaf forests, but the magnitude is still lower than for grasslands (Figure 13, bottom right), and considerably lower than in the case of croplands (Figure 13, bottom left).
FIGURE 13. Categorization of mean annual NDVI anomalies according to the land cover for the selected countries and for the whole FORESEE domain (numbers at the top are percentiles (p; left of 0.5p - most extreme negative anomaly; 0.5p - 2p - extremely negative anomaly; 2p - 9p - very negative anomaly; 9p - 25p - negative anomaly; 25p - 75p - normal range; 50p - median; 75p - 91p positive anomaly; 91p - 98p very positive anomaly; 98p-99.5p extremely positive anomaly normal; right of 99.5p - most extreme positive anomaly); BIH - Bosnia and Herzegovina, CZE - Czech Republic; HRV- Croatia; HUN - Hungary; SVK - Slovakia; SVN - Slovenia).
Differences between countries exist, but they are rather small, with increasing discrepancies toward the extremes (2p and less or 98p and more). However, this could be the result of relatively small number of pixels with extreme values and should be interpreted with caution. This is particularly the case for grasslands, where the number of pixels was large enough only for Bosnia and Herzegovina and Croatia (see Table 3). We would also like to point out the differences among countries in the case of the croplands, and in particular to Slovenia. It can be seen that the high values for the anomalies in NDVI are occurring in Slovenia less frequently than e.g. in Slovakia or in Czech Republic. The reason behind this could be climate (relatively large amount of precipitation in Slovenia - see Figure 2) but also it could be due to the differences in crop type or cultivation method. Somewhat similar behaviour can be observed also for the coniferous forests in Hungary, where the observed magnitude of their positive anomaly is apparently lower than for coniferous forests in the other countries. However, here the caution is needed, because it could be an artefact due to the modest number of available pixels and this would require further investigation. In any case, these results corroborate the logic of country level analysis, where the considered area is large enough to reduce the effects of random noise (caused e.g. by species or management differences), but small enough to retain the information on the exiting differences among the countries. The relatively uniform distance between the values of the mean anomalies for the shown percentiles (Figure 13) indicates that the choice of the percentiles, representing the border between the classes of anomalies (see Section 2.4 in the Materials and Methods for the list) has been made properly.
In an effort to identify specific years that can be characterized as extremes according to a given RS characteristic, we performed a land cover and country specific analysis. The results for Bosnia and Herzegovina, Croatia and Hungary based on NDVI are given in Figures 14-16 and Table 5 while the results for other countries and other RS characteristics are available in the Supplementary materials (Figures S7, S9, S11, S13 and Tables S2-S5). Distribution of shares of NDVI anomalies by categories according to the land cover and years (Figures 14-16) clearly shows the dominant effects of meteorological conditions. However, the responses of various land cover types differ significantly. What immediately strikes the attention is the fact that croplands and grasslands have almost unanimous agreement on which year was the worst (2003) and which one was the best (2014). However, year 2003 cannot be classified as extreme for croplands and grasslands in all countries, although they are relatively close geographically (e.g. compare Hungary and Croatia). When we focus on forests, the situation is rather clear regarding the best year (2014 was the best for all forest type and in all countries except for coniferous forest of Hungary and Slovenia). Here the case of Slovenia is somewhat easier to interpret. Namely, in January-February 2014 large part of Slovenia and part of Croatia was experiencing severe case of freezing rain that caused extensive damage to forests of the affected region [37, 38].
Regarding the negative extreme in the case of forests, the situation is more complex. Unlike for croplands and grasslands, year 2003 was not the worst in all cases except for broadleaf forests in Hungary and coniferous forests in Hungary and Croatia. Interestingly, year 2004 and 2005 are indicated as worst or second worst years for forests in different countries. This could imply that the negative effects of year 2003 (which was in fact an extremely warm and dry year all over Europe  on forest ecosystems were not immediately visible, but the consequences of the unfavourable weather was “remembered” by the ecosystem and manifested itself in the following years. This “dampening” effect of the negative anomalies in the case of forests is very important. It shows how difficult is to make straightforward conclusions based only on the current observation without taking into account the events of the past.
Table 6 shows the summary of country-specific extreme years for the different land cover types and different vegetation indices (based on the approach given in Table 1). Note that NPP is not used here due to problems with the temporal coverage of the dataset (see Section 2.3). The table shows that for herbaceous vegetation and for coniferous forests in Central Europe year 2003 is undoubtedly the most important year that affected vegetation status adversely. Besides year 2003, year 2000 is also notable. Meteorological conditions during those years need deeper evaluation. In contrast, for broadleaf forests year 2005 and 2012 might be notable, but as it was mentioned earlier, there is no clear evidence for unanimous selection of extremely unfavourable year in terms of overall forest development status. Considering above-average plant performance year 2014 is clearly noticeable for herbaceous vegetation and coniferous forests (to a lesser extent). For broadleaf forests the situation is similar to the bad year case, namely there is no unambiguous single good year. 2009, 2011 and 2014 might be studied in some countries to gain deeper understanding of the cause of the positive anomalies. The presented results corroborate that the different vegetation characteristics provide relatively consistent results in terms of anomalous plant behaviour.
TABLE 6. Summary of country-specific extreme years for the different land cover types and different vegetation indices (based on the approach given in Table 1). Note that NPP is not used here due to problems with the temporal coverage of the dataset (see Section 2.3).
In this study we do not attempt to discriminate the different RS-based indices, which means that the indices are not ranked or qualified in any way. NDVI, EVI, LAI and FAPAR are all related to plant greenness and leaf development status, which is in turn related to photosynthetic capacity and plant productivity. GPP and NPP are based on a series of model assumptions and they provide measures that are directly related to productivity. Therefore, we assume that they all serve as proxies of plant processes, and as such we can expect from them to respond to the climate fluctuations. The selection of the relatively large number of indices was done due to their different complexity in their algorithm and assumptions which were used during their calculations.
For example, NDVI and EVI are indices obtained directly from reflectances provided by the polar orbiting MODIS instrument, which can already be related to the potential plant productivity. On the other hand, FPAR and LAI (which are also obtained from remotely sensed data) are also related to plant productivity and plant development status, but unlike the former two, FPAR and LAI are calculated taking into account the biome type derived from the land cover information . Thus they embedded additional information but possibly additional errors as well (either random due to e.g. misclassification of the pixel, or bias due to the possible bias in biome parameters used in FPAR and LAI calculation). Finally, MOD17 products, namely GPP and NPP, are likewise not independent from FPAR  but they also carry additional information (and probably additional errors) due to the meteorology information used for their calculation.
Considering the above, it is clear that all of the indices or metrics potentially have advantages and setbacks and it is difficult to select the best one. Therefore, in our work we decided to investigate all and treat them as an ensemble, similarly to the ensemble technique used in numerical weather prediction.
The selection of land cover specific anomalous years highlighted that plant response to unusual climate conditions strongly depend on the land cover type. Intra-species differences also exist but in this study this was not addressed. The advantage of the country-mean studies is the robustness which was demonstrated earlier in other studies [15, 22]. Using country means the spatial differences are most likely diminished to some extent, and in fact regions with considerably higher anomalous behaviour might exist within the countries that might require further studies.
We need to mention that in the present approach anomalies were defined for the entire growing season. It is clear that the growing season (from start of season to cessation) might be split into different time periods that needs further investigation. These anomalies associated with smaller temporal scale might provide additional information on the nature of anomalies and their meteorological driver.
The presented results might provide invaluable information for researchers associated with plant production (ecologists, agronomists, foresters, etc.). As ecosystem services are closely tied to plant productivity , stakeholders might also find the presented information useful. With the availability of additional years, the time series should be extended and thus the trend in the strength of the anomalies might be estimated.
The present study highlighted the difficulties related with the selection of appropriate time periods and climate indices to define extreme weather from the point of view of ecosystems. The effect-based approach relying on RS data resolves this difficulty as the observed anomalies act as integrators of the past environmental conditions thus they are good indicators of unusual conditions. Explicit selection of anomalous years indicated that the unusual plant state is not independent of the plant functional type. It is also clear that large differences exist between herbaceous and woody vegetation, where the latter is also associated with legacy effect from the previous year . This legacy effect is one major issue which needs further research. In any case, legacy effect (which is well documented) alone questions the pure meteorological approach for the selection of anomalous years, as climate anomalies in the last year might exert stronger impact on the current plant growth than the present year. We propose to further refine the effect-based approach separately for the different land cover types to better understand the main drivers of plant growth in Central Europe and also worldwide.
We thank NASA, for producing and distributing the MOD13 NDVI data. Earth Observing System Data and Information System (EOSDIS) 2009, Earth Observing System ClearingHouse (ECHO) / Reverb Version 10.91.5 [online application], Greenbelt, MD: EOSDIS, Goddard Space Flight Center (GSFC) National Aeronautics and Space Administration (NASA), URL: https://wist.echo.nasa.gov/api/. CORINE Land Cover (CLC) 2012 Database is the property of European Commission - Directorate-General for Internal Market, Industry, Entrepreneurship and SMEs (DG-GROW) and European Environment Agency is data custodian. The CLC database has full, open and free access in line with the Copernicus delegated regulation (EU) No 1159/2013 of 12 July 2013, supplementing Regulation (EU) No 911/2010 of the European Parliament and licensing conditions for GMES users and defining criteria for restricting access to GMES dedicated data and GMES service information. The research has been supported by the Hungarian Scientific Research Fund (OTKA PD-111920 and K-104816), the Croatian Science Foundation (HRZZ UIP-11-2013-2492) and the Széchenyi 2020 programme, the European Regional Development Fund and the Hungarian Government (GINOP-2.3.2-15-2016-00028).
Supplementary material 1. (here)
- IPCC 2012 Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation. In: Field CB, Barros V., Stocker TF, Qin Q, Dokken DJ, Ebi KL, Mastrandrea MD, Mach KJ, Plattner G-K, Allen SK, Tignor M, Midgley PM (eds) A Special Report of Working Groups I and II of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK, and New York, NY, USA, 582 p. URL: https://www.ipcc.ch/pdf/special-reports/srex/SREX_Full_Report.pdf (14 Feb 2016)
- IPCC 2013 Climate Change 2013: The Physical Science Basis. In: Stocker TF, Qin D, Plattner G-K, Tignor M, Allen SK, Boschung J, Nauels A, Xia Y, Bex V, Midgley PM (eds) Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 p. URL: http://www.climatechange2013.org/ (14 Feb 2016)
- IPCC 2014 Climate Change 2014: Synthesis Report. In: Pachauri RK, Meyer LA (eds) Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. IPCC, Geneva, Switzerland, 151 p. URL: https://www.ipcc.ch/report/ar5/syr/ (14 Feb 2016)
- Easterling DR, Evans JL, Groisman PYa, Karl TR, Kunkel KE, Ambenje P 2000 Observed Variability and Trends in Extreme Climate Events: A Brief Review. B Am Meteorol Soc 81 (3): 417-425. DOI: https://doi.org/10.1175/1520-0477(2000)081<0417:OVATIE>2.3.CO;2
- Heim R 2015 An overview of weather and climate extremes – Products and trends. Weather and Climate Extremes 10: 1-9 DOI: https://doi.org/10.1016/j.wace.2015.11.001
- Klein Tank AMG, Können GP 2003 Trends in indices of daily temperature and precipitation extremes in Europe, 1946–99. J Climate 16 (22): 3665-3680. DOI: https://doi.org/10.1175/1520-0442(2003)016<3665:TIIODT>2.0.CO;2
- Bartholy J, Pongrácz R 2007 Global Planet Change 57 (1-2): 83-95. DOI: https://doi.org/10.1016/j.gloplacha.2006.11.002
- Millennium Ecosystem Assessment 2005 Ecosystems and Human Well-being: Synthesis. Island Press, Washington, DC, USA, 137 p. URL: http://www.millenniumassessment.org/documents/document.356.aspx.pdf (25 Sep 2016)
- Hanewinkel M, Cullmann DA, Schelhaas M-J, Nabuurs G-J, Zimmermann NE 2012 Climate change may cause severe loss in the economic value of European forest land. Nat Clim Change 3 (3): 203–207. DOI: https://doi.org/10.1038/nclimate1687
- Hlásny T, Trombik J, Dobor L, Barcza Z, Barka I 2016 Future climate of the Carpathians: climate change hot-spots and implications for ecosystems. Reg Environ Change 16 (5): 1495-1506. DOI: https://doi.org/10.1007/s10113-015-0890-2
- Hlásny T, Trombik J, Bosel aM, Merganic J, Marusák R, Seben V, Stepánek P, Kubista J, Trnka M 2017 Climatic drivers of forest productivity in Central Europe. Agr Forest Meteorol 234: 258-273. DOI: https://doi.org/10.1016/j.agrformet.2016.12.024
- Hovenden MJ, Newton PCD, Wills KE 2014 Seasonal not annual rainfall determines grassland biomass response to carbon dioxide. Nature 511: 583-586. DOI: https://doi.org/10.1038/nature13281
- Reichstein M, Bahn M, Ciais P, Frank D, Mahecha MD, Seneviratne SI, Zscheischler J, Beer C, Buchmann N, et al. 2013 Climate extremes and the carbon cycle. Nature 500: 287-295. DOI: https://doi.org/10.1038/nature12350
- Yi C, Pendall E, Ciais P 2015 Focus on extreme events and the carbon cycle. Environ Res Lett 10 (7): 70201. DOI: https://doi.org/10.1088/1748-9326/10/7/070201
- Bognár P, Kern A, Pásztor Sz, Lichtenberger J, Koronczay D, Ferencz Cs 2017 Yield estimation and forecasting for winter wheat in Hungary using direct broadcast MODIS data. Int J Remote Sens 38 (11): 3394-3414. DOI: https://doi.org/10.1080/01431161.2017.1295482
- Justice CO, Vermote E, Townshend JRG, Defries R, Roy DP, Hall DK, Salomonson VV, Privette JL, et al. 1998 The Moderate Resolution Imaging Spectroradiometer (MODIS): land remote sensing for global change research. Trans Geosci Remote Sens 36 (4): 1228-1249. DOI: https://doi.org/10.1109/36.701075
- Didan K, Munoz AB, Solano R, Huete A 2015 MODIS Vegetation Index User’s Guide (MOD13 Series) Version 3.00, June 2015 (Collection 6). URL: https://vip.arizona.edu/documents/MODIS/MODIS_VI_UsersGuide_01_2012.pdf (25 May 2015)
- Knyazikhin Y, Glassy J, Privette JL, Tian Y, Lotsch A, Zhang Y, Wang Y, Morisette JT, PVotava, et al. 1999 MODIS Leaf Area Index (LAI) and Fraction of Photosynthetically Active Radiation Absorbed by Vegetation (FPAR) Product (MOD15) Algorithm Theoretical Basis Document. URL: http://eospso.gsfc.nasa.gov/atbd/modistables.html (15 Sept 2016)
- Running SW, Zhao M 2015 User’s Guide - Daily GPP and Annual NPP (MOD17A2/A3) Products NASA Earth Observing System MODIS Land Algorithm, Ver. 3.0 for Coll. 6, URL: http://www.ntsg.umt.edu/sites/ntsg.umt.edu/files/modis/MOD17UsersGuide2015_v3.pdf (14 Jan 2015)
- Teuling AJ, Seneviratne SI, Stöckli R, Reichstein M, Moors E, Ciais P, Luyssaert S, van den Hurk B, et al. 2010 Contrasting response of European forest and grassland energy exchange to heatwaves. Nat Geosci 3: 722-727. DOI: https://doi.org/10.1038/ngeo950
- Büttner G, Maucha G 2006 The Thematic Accuracy of CORINE Land Cover 2000. Assessment Using LUCAS (Land Use/Cover Area Frame Statistical Survey), European Environment Agency, Copenhagen, EEA Technical Report 7, 85 pp. URL: http://land.copernicus.eu/user-corner/technical-library/technical_report_7_2006.pdf (1 Feb 2016)
- Kern A, Bognár P, Pásztor Sz, Barcza Z, Timár G, Lichtenberger J, Ferencz Cs 2015 Monitoring the state of vegetation in Hungary using 15 years long MODIS Data. Geophysical Research Abstracts Vol. 17, 13613. URL: http://meetingorganizer.copernicus.org/EGU2015/EGU2015-13613.pdf (14 Nov 2016)
- Zscheischler J, Mahecha MD, von Buttlar J, Harmeling S, Jung M, Rammig A, Randerson JT, Schölkopf B, et al. 2014 A few extreme events dominate global interannual variability in gross primary production. Environ Res Lett 9 (3): 035001. DOI: https://doi.org/10.1088/1748-9326/9/3/035001
- Dobor L, Barcza Z, Hlásny T, Havasi Á, Horváth F Ittzés P, Bartholy J 2014 Bridging the gap between climate models and impact studies: The FORESEE Database. Geosci Data J 2 (1): 1-12. DOI: https://doi.org/10.1002/gdj3.22
- Kern A, Marjanović H, Barcza Z 2016 Evaluation of the quality of NDVI3g dataset against Collection 6 MODIS NDVI in Central-Europe between 2000 and 2013. Remote Sens 8 (11): 955. DOI: https://doi.org/10.3390/rs8110955
- LP DAAC (Land Processes Distributed Active Archive Center). MOD13A2 and MYD13A2, Collection 6. NASA EOSDIS Land Processes DAAC, USGS Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota. URL: https://lpdaac.usgs.gov (20 Oct 2016)
- EOSDIS, Earth Observing System ClearingHouse (ECHO)/Reverb Version 10.91.5 [online application]. Greenbelt, MD: EOSDIS, Goddard Space Flight Center (GSFC) National Aeronautics and Space Administration (NASA). URL: http://reverb.echo.nasa.gov/reverb (20 Oct 2016)
- Viovy N, Arino O, Belward AS 1992The Best Index Slope Extraction: A method for reducing noise in NDVI time-series. Int J Remote Sens 13 (8): 1585-1590. DOI: https://doi.org/10.1080/01431169208904212
- Strahler A 1999 MODIS Land Cover Product Algorithm Theoretical Basis Document (ATBD) Version 5.0. MODIS Land Cover and Land-Cover Change. URL: https://modis.gsfc.nasa.gov/data/atbd/atbd_mod12.pdf (13 Jun 2016)
- EEA 2016 Co‑ORdinated INformation on the Environment (CORINE) Land Cover 2012, Version 18.4. European Commission - Directorate-General for Internal Market, Industry, Entrepreneurship and SMEs (DG-GROW, data owner), European Environment Agency (EEA, data custodian). URL: http://land.copernicus.eu/pan-european/corine-land-cover/clc-2012 (24 Apr 2016)
- Friedl MA, Sulla-Menashe D, Tan B, Schneider A, Ramankutty N, Sibley A, Huang X 2010 MODIS Collection 5 global land cover: Algorithm refinements and characterization of new datasets. Remote Sens Environ 114 (1): 168-182. DOI: https://doi.org/10.1016/j.rse.2009.08.016
- Pérez-Hoyos A, García-Haro FJ, San-Miguel-Ayanz J 2012 Conventional and fuzzy comparisons of large scale land cover products: Application to CORINE, GLC2000, MODIS and GlobCover in Europe. ISPRS J Photogramm 74: 185-201. DOI: https://doi.org/10.1016/j.isprsjprs.2012.09.006
- Hao G, Gen-Suo JIA 2014 Assessing MODIS Land Cover Products over China with Probability of Interannual Change. Atmos Ocean Sci Lett 7: 564-570. DOI: https://doi.org/10.1080/16742834.2014.11447225
- Herold M, Mayaux P, Woodcock CE, Baccini A, Schmullius C 2008 Some challenges in global land cover mapping: an assessment of agreement and accuracy in existing 1km datasets. Remote Sens Environ 112 (5): 2538-2556. DOI: https://doi.org/10.1016/j.rse.2007.11.013
- Klepac D 2000 The largest pedunculate oak forest in Croatia, Spačva (in Croatian). Croatian Academy for Science and Arts, The Center for Scientific Work in Vinkovci: Vinkovci, Zagreb, Croatia, 2000, 116 p
- Semenov MA, Shewry PR 2011 Modelling predicts that heat stress, not drought, will increase vulnerability of wheat in Europe. Scientific Reports 1: 66. DOI: https://doi.org/10.1038/srep00066
- Veselič Ž, Grecs Z, Kolšek M, Oražem D, Matijašić D, Beguš J 2015 Icebreak in Slovenian forests and salvation consequences (In Slovenian with English summary). UJMA 29: 188-194. URL: http://www.sos112.si/slo/tdocs/ujma/2015/188_194.pdf (22 Dec 2016)
- Vuletić D, Kauzlarić Ž, Balenović I, Krajter Ostojić S 2015 Assessment of Forest Damage in Croatia Caused by Natural Hazards in 2014. South-East Eur For 5 (1): 65-79. https://doi.org/10.15177/seefor.14-07
- CIAIS P, REICHSTEIN M, VIOVY N, GRANIER A, OGÉE J, ALLARD V, AUBINET M, BUCHMANN N, et al. 2005 Europe-wide reduction in primary productivity caused by the heat and drought in 2003. Nature 437: 529-5 DOI: https://doi.org/10.1038/nature03972
- Frank D, Reichstein M, Bahn M, Thonicke K, Frank D, Mahecha MD, Smith P, van der Velde M, et al. 2015 Effects of climate extremes on the terrestrial carbon cycle: concepts, processes and potential future impacts. Glob Change Biol 21 (8): 2861-2880. DOI: https://doi.org/10.1111/gcb.12916
© 2017 by the Croatian Forest Research Institute. This is an Open Access paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0).