Next Article in Journal
Change Alignment-Based Image Transformation for Unsupervised Heterogeneous Change Detection
Next Article in Special Issue
A Record-Setting 2021 Heat Wave in Western Canada Had a Significant Temporary Impact on Greenness of the World’s Largest Protected Temperate Rainforest
Previous Article in Journal
Msplit Estimation Approach to Modeling Vertical Terrain Displacement from TLS Data Disturbed by Outliers
Previous Article in Special Issue
Delineating Fire-Hazardous Areas and Fire-Induced Patterns Based on Visible Infrared Imaging Radiometer Suite (VIIRS) Active Fires in Northeast China
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Critical Climate Periods Explain a Large Fraction of the Observed Variability in Vegetation State

1
Space Research Group, Department of Geophysics and Space Sciences, Institute of Geography and Earth Sciences, ELTE Eötvös Loránd University, H-1117 Budapest, Hungary
2
Faculty of Forestry and Wood Sciences, Czech University of Life Sciences Prague, 165 21 Prague 6, Czech Republic
3
Institute of Geography and Earth Sciences, Department of Meteorology, ELTE Eötvös Loránd University, H-1117 Budapest, Hungary
4
Excellence Center, Faculty of Science, ELTE Eötvös Loránd University, H-2462 Martonvásár, Hungary
5
Doctoral School of Environmental Sciences, ELTE Eötvös Loránd University, H-1117 Budapest, Hungary
6
Centre for Agricultural Research, Agricultural Institute, H-2462 Martonvásár, Hungary
7
Lechner Knowledge Center, H-1149 Budapest, Hungary
8
Croatian Forest Research Institute, Department of Forest Management and Forestry Economics, HR-10450 Jastrebarsko, Croatia
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(21), 5621; https://doi.org/10.3390/rs14215621
Submission received: 30 September 2022 / Revised: 1 November 2022 / Accepted: 5 November 2022 / Published: 7 November 2022

Abstract

:
Previous studies have suggested that a major part of the observed variability in vegetation state might be associated with variability in climatic drivers during relatively short periods within the year. Identification of such critical climate periods, when a particular climate variable most likely has a pronounced influence on the vegetation state of a particular ecosystem, becomes increasingly important in the light of climate change. In this study, we present a method to identify critical climate periods for eight different semi-natural ecosystem categories in Hungary, in Central Europe. The analysis was based on the moving-window correlation between MODIS NDVI/LAI and six climate variables with different time lags during the period 2000–2020. Distinct differences between the important climate variables, critical period lengths, and direction (positive or negative correlations) have been found for different ecosystem categories. Multiple linear models for NDVI and LAI were constructed to quantify the multivariate influence of the environmental conditions on the vegetation state during the late summer. For grasslands, the best models for NDVI explained 65–87% variance, while for broad-leaved forests, the highest explained variance for LAI was up to 50%. The proposed method can be easily implemented in other geographical locations and can provide essential insight into the functioning of different ecosystem types.

Graphical Abstract

1. Introduction

The leaf development stage is a major determinant of gross photosynthesis, carbon balance and biomass accumulation of plants. The amount and condition of green leaves are influenced by genetic factors, environmental conditions including biotic and abiotic factors, and other disturbances, such as management [1,2,3,4]. Observations revealed that vegetation greenness and productivity have large interannual and spatial variability, which is associated with climatic fluctuations [4,5,6,7], such as large-scale drought events or heat waves [8,9]. Understanding and quantifying the causes of the observed variability is a major challenge for the scientific community due to the interaction between the multiple drivers [10].
It has been demonstrated that the variability of annual productivity and vegetation state can be associated with environmental conditions during shorter periods within the year, instead of the annually aggregated climate conditions [11,12,13,14]. This indicates that not only the environmental effect itself but also the timing of the climatic fluctuation can be critical factors for the annual productivity and vegetation state [15,16,17]. It was also demonstrated that the vegetation state (i.e., intra-annual variations) can be affected by the environmental conditions in the antecedent periods with different time lags spanning from months to years [6,17,18,19,20,21,22,23]. These short periods during the year, when the environmental conditions have a profound role in the forthcoming vegetation state and overall annual productivity, are referred to as critical climate periods [11,12,24,25,26]. The existence of critical climate periods within the year and the associated time lags are less studied for different climate zones and ecosystem types.
Climate change will have many facets, and the most obvious one will be the increase in the frequency of extreme events and their devastating consequences on various ecosystems [5]. However, even smaller changes, such as local shifts in the precipitation pattern, changes in temperature, radiation, vapor pressure deficit and soil water content, etc. will affect vegetation [6,15]. Plants will probably attempt to adapt to the shifts in the meteorological conditions, but the ability of different plant species and ecosystems to adapt varies, and is still greatly unknown [7,10]. Therefore, the identification of the periods during the year when a particular climate variable is most likely to have pronounced influence on the vegetation state of a particular ecosystem has high relevance.
Satellite-based remote sensing is a powerful tool to study the vegetation state and overall productivity of ecosystems at large spatial scales with high accuracy [8,10,20,27,28,29,30,31]. The identification of ecosystem-specific critical climate periods under diverse climatic conditions has been hindered so far because detailed ecosystem maps have seldom been available for large areas, such as landscapes or countries. Therefore, remote-sensing-based studies have usually focused on generic land cover types (e.g., deciduous broadleaf forests, grasslands; [32]) instead of particular ecosystems (e.g., habitat types or dominant species). In recent years, new generation ecosystem maps have been published in several countries of Central Europe, such as the Ecosystems in Slovakia dataset [33], the dataset of Czechia [34], the Croatia dataset [35], or the NÖSZTÉP database of Hungary [36]. These datasets offer new perspectives in the vegetation-related, remote-sensing studies focusing on ecosystem-specific productivity, carbon balance, stability, and vulnerability in a way that is unprecedented, and can provide insight to a better understanding of plant functioning.
In this study, we present a comprehensive analysis of the Hungarian forests and grasslands using the remote-sensing-based Normalized Difference Vegetation Index (NDVI) and Leaf Area Index (LAI) as indicators of plant functioning, provided by the measurements of the Terra/MODIS sensor [37] with 500 m spatial resolution. NDVI and LAI were used to identify critical climate periods on an intra-annual time scale during 2000–2020. Studying the response of the ecosystems to extreme environmental conditions is particularly important in light of the ongoing climate change, and Hungary, as one of the driest countries in Central Europe, offers an ideal study area.
The present study focuses on different ecosystem categories within deciduous broadleaf forests and grasslands, representing semi-natural ecosystem categories in Hungary (located in the Pannonian biogeographical region). Semi-natural in this context means that those ecosystems might include plantations and are typically affected by some kind of management, but they represent otherwise perennial ecosystems not affected by land-use change and are not part of the annual cropping systems. Deciduous broadleaf forests have a ~23% share in the country [36] and are subject to disturbances due to pests and extreme weather, with some of them located at the receding edge of the distribution, which indicates vulnerability to climate change [38]. Grasslands cover about 11% of the total area of Hungary [36,39] and are characterized by broad biodiversity due to the special pedo-climatic conditions, forming the oldest managed grasslands in Europe [40]. Grassland ecosystems are drought-prone during summer in the majority of the country [41], which highlights the issue of their vulnerability and long-term stability with climate change.
In our study, we address the following questions: (1) What is the long-term mean NDVI and LAI and the intra-annual and interannual variability for the investigated forest and grassland ecosystems? (2) What is the relationship between the foliar development stage and the antecedent environmental conditions with a relatively short time lag? (3) Which are the critical climate periods with short time lags during the year that are predominantly responsible for the observed interannual variability of NDVI and LAI?

2. Materials and Methods

2.1. Study Area

The present study focuses on Hungary, located in Central Europe in the Pannonian biogeographical region (BGR) [40], which almost completely corresponds to the area of Hungary. The climate of the area is mostly continental (with warm or hot summers) according to Köppen’s classification scheme, with high variability in the meteorological conditions, and a tendency toward summer droughts [41]. Based on the meteorological dataset used (see Section 2.3), the mean annual temperature during 2000–2020 was ~11.7 °C, while the mean annual precipitation ranged from 750 mm in the south-west down to less than 500 mm in the central part of the study area, which is one of the driest regions of Central Europe. The mean elevation of Hungary based on SRTM data [42] is 150 m above sea level.

2.2. Remote-Sensing-Based NDVI and LAI Datasets

Vegetation indices (VI) are designed to provide information on the greenness of the vegetation, which is closely related to foliar development, and the overall state of vegetation [29]. The Normalized Difference Vegetation Index (NDVI) and Leaf Area Index (LAI) at 500 m spatial and 8-day temporal resolution, derived from the data of the MODIS sensor on-board satellite Terra [37], were used for the 2000–2020 study period. Collection 6 MODIS products were downloaded from NASA LP DAAC for the tile h19v04 [43].
NDVI was derived from the surface reflectances of Band-1 (ρRED; visible red) and Band-2 (ρNIR; near-infrared) of the MODIS sensor, included in the MOD09A1 product [44], as (ρNIRρRED)/(ρNIR + ρRED). The quality filtered and smoothed NDVI time-series at an 8-day temporal grid was created based on the work of Kern et al. [30,31]. To avoid any misleading and non-realistic sudden decrease in the NDVI due to the presence of unrecognized and incorrectly quality flagged atmospheric effects (such as high aerosol content, presence of thin cirrus, sub-pixel clouds, etc.), the Best Index Slope Extraction (BISE) method [45] was applied at the pixel-level. The quality-checked and gap-filled NDVI dataset was then resampled into daily resolution using linear interpolation, taking into account the actual Julian dates of the measurements, indicated in the temporal composite products. Finally, the Savitzky–Golay filter [46] was also applied with a 30-day window and a second-degree polynomial smoothing to gain smoothed daily datasets. From the daily NDVI time-series, a dataset with a regular 8-day resolution was created, corresponding to the 8-day equidistant resolution of the MODIS LAI dataset. We refer to these 8-day long intervals of the annual cycle as periods. The temporal grid of the resulting 46 data per year served as a basic temporal resolution in our research.
Unlike NDVI, which is calculated only from the reflectances, LAI is obtained with more sophisticated calculation methods, which include modeling elements [47]. In this study, we used the MODIS LAI [m2 m−2] dataset, contained within the MOD15A2H data product [48]. The MODIS LAI values represent a one-sided green leaf area per unit ground area in broadleaf canopies and one-half the total needle surface area per unit ground area in coniferous canopies. Both are the results of the applied algorithm, which uses daily L2G-lite surface reflectance, and also biome type [47].
For the pre-processing of the LAI dataset, the QC flag and the additional “Extra QC flag” information within the HDF files were used. Only data with QC flags 0 and 32 along with Extra QC flags 0, 1, 8, 9, 128, 129, 136, and 137, and the snow/ice covered versions (values shifted with 4) were accepted. Allowance of the snow-/ice-covered pixels (similarly to the method for NDVI) was important to gain reliable spring phenology in mountainous areas, where otherwise no valid data point would have been present. Data with other QC flags were discarded and temporally filled by linear interpolation at the pixel-level. After QC filtering, the BISE method was applied to eliminate the remaining sudden decreases in the LAI courses.

2.3. Meteorological and Environmental Datasets

To study the effects of the weather on the state of the vegetation, we used meteorological data for the period 2000–2020. Daily minimum (Tmin) and maximum (Tmax; °C) temperature, precipitation (Prec; mm day−1), daylight average vapor pressure deficit (VPD; Pa), and daylight average shortwave radiative flux (i.e., global radiation, Rad; W m−2) were obtained from the FORESEE v3.2 dataset [49,50]. Daily mean soil water content (SWC2; m3 m−3) was calculated from the hourly volumetric soil water content data for the second soil layer (0.07–0.28 m depth) of the ERA5-Land dataset [51,52].
The daily data stored at the original FORESEE grid (1/6° × 1/6°) or ERA5-Land grid (0.1° × 0.1°) were resampled to the finer spatial grid of the MODIS products and averaged to the regular 8-day temporal grid. In the cases of temperature and precipitation, the resampling was performed based on the methodology of Kern et al. [53], where the elevation of the pixels was also taken into account.
In our study, the climate variables were used to detect the climatology of plant-weather interactions (see 2.6 and 2.7). Although co-linearity might exist between the climate variables (primarily between Tmin and Tmax; Tmax and VPD; Rad and SWC2; as it was investigated in our previous study [13]), we assume that all variables affect plant processes differently during different stages of plant growth. This means that it is not unreasonable to investigate them separately.

2.4. Ecosystem Category and Land Cover Datasets

The classification into ecosystem categories was based on the species/habitat type distribution information available in the Ecosystem Map of Hungary of the NÖSZTÉP project [34,54]. This map (hereafter referred to as the NÖSZTÉP dataset) has 56 different Level 3 categories at 20 m spatial resolution, covering the entire area of Hungary. Resampling of the dataset to the grid of the QKM (250 m) and HKM (500 m) MODIS products resulted in values of the actual shares (in percent) of each ecosystem category for every pixel [31].
We selected reliable pixels with pure land cover types from homogeneous areas. The determination of the pure pixels was based on at least a 90% share threshold of a given NÖSZTÉP category within a pixel. To ensure that the selected pixels were situated in a largely homogenous environment, we applied an additional condition. That is, each of the eight neighboring QKM pixels around each HKM pixel should have contained at least 60% of the same category. This criterion was fundamental to reducing the noise introduced by the potential presence of other ecosystem categories in the neighboring pixels, which would affect the recorded reflectance due to the known geolocation inaccuracy of the MODIS measurements [55], or due to the artefacts of the gridding procedure [56] and possible re-projection inaccuracies.
To reduce the effects of possible land cover change during 2000–2020 relative to the base year of the NÖSZTÉP dataset (2015–2017), the CLC2000 and CLC2012 [57] land cover datasets resampled to the HKM MODIS grid were used, where only pixels that met the criteria of the same share threshold were retained for further analyses. In the case of broad-leaved forests, an additional NDVI threshold (0.8) was also used to discard unstable pixels with a growing season mean NDVI of less than 0.8 in any of the years during the study period. Finally, only the ecosystem categories that had at least 10 pixels were used (Table 1 and Figure 1). We also created a group called All pedunculate oak (n = 15) by merging multiple oak categories in which pedunculate oak was also present, due to the importance of pedunculate oaks and the lack of a pure category. The selected ecosystems, with the exception of All pedunculate oak, followed a distinct precipitation-temperature gradient from wetter & cooler to dryer & warmer (Figure S1).
All of the selected ecosystem categories are of significant socio-economic importance, not only for Hungary but also at the Central European scale, as some of them are fairly unique (e.g., Open sandy grassland). Therefore, despite the low pixel numbers in some cases, investigation thereof is highly relevant.

2.5. Statistical Analysis of the Remote Sensing-Based Products

Different statistical measures were calculated from NDVI and LAI for each ecosystem category. Area-averaged (AAy,p) values were derived for each year y and period p based on the selected pixels (n) of the given ecosystem category (Equation (1)).
A A p , y = i = 1 n I n d e x i , p , y n ,           I n d e x { N D V I ,   L A I }  
Multiannual mean NDVI and LAI curves (MAMp,i) were calculated separately for each pixel i and each period p (Equation (2)) and also at the level of ecosystem categories (MAMp, Equation (3)) as the average of all pixel-level multiannual means for each 8-day period and for the 21 years:
M A M i , p = y = 1 21 I n d e x i , p , y 21 ,         I n d e x { N D V I ,   L A I }
M A M p = i = 1 n M A M i , p n = y = 1 21 A A p , y 21 ,
The interannual variability (IAVp) for every period p was defined by Equation (4) as the standard deviation ( σ A A p , y ) of the yearly country averaged mean values (AAp,y):
I A V p = σ A A p , y = y = 1 21 ( A A p , y M A M p ) 2 21 1 ,
The standard error of the MAMp values ( S E M A M p ) were also determined (Equation (5)) at the period level based on the standard deviations ( σ M A M p ) of the yearly country-averaged mean value (AAp,y) to indicate the interannual variability in each period for each ecosystem category:
S E M A M p = σ A A p , y 21 ,
Spatial variability (SVp) for every period p was defined (Equation (6)) as the standard deviation ( σ M A M i , p ) of the pixel-level MAMi,p values:
S V p = σ M A M i , p = i = 1 n ( M A M i , p M A M p ) 2 n 1 ,
We tested the difference in the MAM values between the different ecosystem categories at each period during the year through a one-way analysis of variance (ANOVA).

2.6. Critical Climate Period Analysis

Plant physiology is strongly influenced by environmental conditions, and also by the timing and intensity of human intervention (i.e., management). Figure 2 presents a conceptual view of the Central European ecosystems in terms of the typical effects of environmental conditions and management during the phenophases of broad-leaved forests and grasslands. The environmental conditions presumably affect the particular parts of the phenological cycle differently, which must be considered during the interpretation.
To quantify the effects of the environmental conditions on vegetation state (NDVI and LAI), statistical analysis was performed using moving-window correlation analysis at 8-day temporal resolution. The same analysis was performed for NDVI and LAI separately. For simplicity, hereafter we refer only to NDVI. For each period in a year, a 32-day (four periods) moving mean area-averaged NDVI along with climate data anomaly values were calculated consecutively. The 32-day-long area-averaging proved to be a robust method for the identification of relationships at large spatial scales [13]. After that, linear correlation coefficient (Pearson’s R) values were calculated between the obtained 32-day mean area-averaged anomaly values of NDVI and the 32-day mean climate variables, similarly to previous studies [6,14,18,20].
Since the vegetation state is affected by preceding environmental conditions, we used different lags between them to identify the most influential climate periods. For each 8-day period, five distinct R-values were calculated using five different time lags between the 32-day mean of NDVI and the 32-day mean of a given climate variable, where the time lag is the difference between the starting day of the climate variable and the starting day of the vegetation state interval (Figure 3). The first R-value corresponds to a time lag of three periods (24 days) with one period (8 days) overlap between the vegetation state and the climate variable. The additional four R-values were calculated by incrementally increasing the time lag by two periods to capture any possible lag effect, with a maximum time lag of almost three months (Figure 3). Since the study focused on the short-term effects, larger time lags (i.e., carry-over effects spanning many months and one or more years) were not included in the analysis. A significant R-value (at p ≤ 0.01 or p ≤ 0.05) indicates that the vegetation state co-varies with the environmental conditions of the corresponding previous time interval. This represents a statistical relationship, but does not necessarily indicate causality (it might be a statistical by-product caused by the co-linearity of the climate data variability).
The most influential 8-day climate period for a given ecosystem category was identified as the period with the highest number of instances when the 32-day mean climate variable (overlapping that 8-day period and for any of the five predetermined lags) and the corresponding 32-day mean vegetation state correlated significantly (at p ≤ 0.01 or p ≤ 0.05, see Figure 3). The importance of a given climate period increases with the number of significant correlation instances, i.e., with the significant correlation frequency (SCF) (Figure 3). Periods with high SCF likely have a determinant role in the upcoming ecosystem dynamics, thus they could be considered influential or even critical climate periods, as described in the literature [11,12,24,25]. In our case, the maximum number of instances when a period could be influential is 20 (except during the start and the end of the year), due to the five different lags and four periods (32 days) used in the correlation calculation.
The average time lag between the NDVI and a given climate variable for each ecosystem category was calculated in the following way. For every 8-day period, the time lag associated with the largest significant correlation coefficient was identified. In case none of the correlations was significant, no time lag was identified. The mean time lag was calculated as the mean of the identified time lags during the year.
To quantify the relative effects of climate variables on the state of the vegetation (NDVI and LAI), the relative effect (REp) values were calculated (Equation (7)) for every ecosystem category and for every period based on the first R-curve (time lag of 24 days):
R E p = S p f M A M p 100 ,
where Sp is the slope of the linear regression for a given period p in the years 2000–2020 between the 32-day long area-averaged NDVI values of the following 32 days (starting at the period p) and the area-averaged climate variable in the preceding 32 days (starting at the period p-3, where 3 is the lag expressed in the number of 8-day periods), MAMp is the multiannual mean of the 32-day area-averaged NDVIp starting at the period p, and f is a factor depending on the climate variable used to express the relative effect in percent within a sensible physical range. The value of f was set as 1 °C, 1 °C, 5 mm, 0.01 m3 m−3, 10 W m−2 and 50 Pa for Tmin, Tmax, Prec, SWC2, Rad and VPD, respectively. Therefore, the calculated relative effect is the normalized slope (normalized by the mean NDVIp value), expressing the sensitivity of the given ecosystem category to the investigated climate variable. The calculated value of the relative effect (Equation (7)) is the percent change of NDVI for a predefined unit change in a given environmental condition (when other effects are not considered).

2.7. Model Construction

The relative effect of a climate variable on NDVI (Equation (7)) does not take into account the interactions and possible co-linearity between different climate variables. This may mask the true effects of that climate variable [13,58]. In order to further analyze the effects of the climate fluctuations on NDVI anomalies in a selected period of interest (PI), multiple linear regression models were constructed for the mean NDVI and LAI, using climate variables during 2000–2020 as predictors. For the exact timing of the PI, the late summer was selected based on the magnitude of the NDVI variability (see Section 3.1). The models were built based on area-averaged values for each ecosystem category separately. The general form of the considered model to estimate mean NDVI (NDVIPI) and in the PI was:
N D V I P I = A P ( α A P · T m a x A P + β A P · T m i n A P + γ A P · P r e c A P + δ A P · S W C 2 A P + ε A P · R a d A P + ζ A P · V P D A P ) + η · Y e a r + c o n s t a n t ,
where TminAP, TmaxAP, PrecAP, SWC2AP, RadAP, and VPDAP, are the minimum and maximum temperature, precipitation, soil water content, radiation, and VPD for a given averaging period (AP) before the PI, respectively; Year is the calendar year of the observation accounting for a possible trend in the NDVIPI; constant is the constant term of the fit; and αAP, βAP, γAP, δAP, εAP, ζAP, and η are model coefficients. The time lags (i.e., climate variables of the previous periods) were introduced to account for possible short-term lag effects.
To reduce the large number of independent variables in the model (Equation (8)), the Boruta method [59] was used to select the climate variables in which the five former periods considered had relevant roles in affecting the vegetation state during the predefined period. The Boruta method quantifies the importance of the driving climate variables using random forest classification. In our case, Boruta was used with the six mean climate variables averaged over 32 days for the same five periods with different time lags, similarly to the earlier described R-curves (Section 2.6).
The constructed models were calibrated and then validated by the application of the leave out one year (LOOY) cross-validation technique [60]. In the LOOY validation, the model was calibrated with a dataset from which one calendar year of the data was omitted year by year. Model predictions were then tested against the observations in the year that was left out, and the procedure was repeated for all years of the validation dataset.
Processing of all the datasets for the present study and execution of the calculations were performed using the Interactive Data Language (IDL) version 8.6 (Harris Geospatial Solutions, USA, Broomfield, CO), STATA 14.2 (StataCorp., USA, College Station, TX, USA) and R [61].

3. Results

3.1. Multiannual Mean (MAM) and Interannual Variability (IAV) Curves

Based on the MAM curves of the investigated ecosystem categories (Figure 4) the category Beech has the highest MAM. The category Black locust dominated plantation shows a much slower green-up during the spring, due to its long-lasting flowering. In the case of herbaceous vegetation, the category Open sandy grassland shows very similar MAM to the Grassland on saline soil during the growing season, but not during the spring green-up. The ranges showing the standard error (SE) of the MAM indicate much higher interannual variability for all grasslands than for forests. The ranges of the ±1.96 SE (corresponding to a 95% confidence interval) are also shown for all periods. The ANOVA-based results of the assessments for the differences between the MAM curves (not shown) were in agreement with the periods when the error ranges of the individual MAM curves overlapped. For forests, those overlapping periods primarily corresponded to the spring green-up and autumn senescence of the vegetation.
In the case of forests, the IAV shows two peaks corresponding to the overall timing of the green-up and senescence (Figure 4) in accordance with our earlier study [30], where the mean start and end of the green-up were at DOY 100 and DOY 123, respectively. Those variabilities were much higher (even up to 3–4 times) than the variability at any other time during the growing season. On the contrary, for the grasslands, the IAV was higher during the late summer–early autumn periods (the highest occurred typically between DOYs 210–270) than in the transitional periods, which was especially striking in comparison to forests. The highest variabilities were associated with the category Grassland on saline soil (Figure 4).
Based on the results, we focused on a 32-day-long period of interest during the late summer (DOY 225–256, 13 August–13 September), when the effects of forest green-up and senescence were not present and the investigated grassland categories showed the highest interannual variability in their NDVI-based vegetation state. During the selected period of interest (PI), the observed interannual NDVI variability of the investigated forest categories, albeit comparatively small, was likely linked to climate effects. This PI served for model constructions (see Section 3.4 and Section 3.5).

3.2. Relationships between the NDVI/LAI and the Climate Variables

Based on the annual curves of the Pearson’s R-values, it is clear that Tmin and Tmax have different effects on the vegetation state of different ecosystem categories described by NDVI (Figure 5) and LAI (Figure S3), demonstrating the importance of their joint application.
Tmin was important for all investigated ecosystem categories during the green-up (as higher daily Tmin seemed to promote the green-up), while for forests it was also important during the senescence (as higher daily Tmin values contributed to a prolonged and longer senescence). In the case of the grasslands during the summer and early autumn, significant negative R-values were found for Tmax, while no significant relationships were found for the forests. The unambiguous role of the SWC2 for grasslands was obvious, as indicated by the stable and high positive R-values during the whole growing season starting from May. Weaker but also positive SWC2-related correlation was present during the summer for the category Black locust dominated plantation, while for forests, no significant positive relationships were found in the same period. The significant negative R-values of the forests in May imply a determinant role of SWC2 in leaf development, but this could also be the effect of co-linearity between SWC2 and other climate variables. That is, due to the likely high SWC2 and its comparatively slow temporal change in spring, SWC2 probably better represented the overall environmental conditions. The R-values representing the correlation between NDVI and Prec and NDVI and SWC2 showed a similar pattern, although for Prec, a stronger variability was present, probably due to the sporadic nature of precipitation (timing and amount). R-curves for VPD were similar to those for Tmax. In the case of Rad, the R-values were mostly negative, but generally not significant. Only the forest ecosystem categories were associated with significant positive R-values during the spring green-up, and only the Beech during the autumn senescence, while for grasslands, the negative R-value was significant during parts of summer. The maximum and minimum values of the significant R-values of the first R-curve (time lag of 24 days) during the year are presented in Table 2 for NDVI (Table S1 for LAI) for each climate variable and ecosystem category, in accordance with Figure 5. For LAI, the R-curves were similar to those made with NDVI, with some differences (Figure S3), mostly affecting the non-significant values. Hereafter, we focus primarily on NDVI, which is the most common vegetation-related characteristic in the literature.
The annual courses of the R-values referring to different time lags separately for each ecosystem category clearly show that the application of longer time lags can provide a similar or even greater correlation (Figures S4–S9). For example, considering Tmin/Tmax and VPD for the category of Black locust dominated plantation, the NDVI anomalies in late July–early August were not affected by the temperatures of the directly antecedent period, but were more affected by the earlier temperatures, with a time lag of three months (i.e., temperatures in May). The number of cases when the first R-curve did not have the highest (absolute) R-value (implying a longer lagged effect) was small. These figures also show the situation when the state of the vegetation was influenced not only by the immediately preceding conditions but also by the earlier ones. For grasslands, the strongest influential period of SWC2 was always the one just before the investigated period, but the SWC2 of the earlier periods also contributed to the NDVI anomaly of the investigated period (Figure S6). On the other hand, in the case of the forests, the category Beech showed the strongest negative correlation during the spring (April–May), with effects that could be long-lasting, even in June (Figure S7). This may be important, considering that climate change will likely affect precipitation amounts and patterns, and consequently also spring SWC2.

3.3. Critical Climate Periods

A graphical overview of the distribution of significant correlation frequencies at different periods during the year between NDVI and a given climate variable enables easy identification of critical climate periods (Figure 6). In the case of NDVI and Tmin, a high frequency of positive correlations can be observed before the spring green-up for all of the studied ecosystem classes, and in the case of forests, also before the autumn senescence (Figure 6a). A similar pattern can be observed also for NDVI and Tmax (Figure 6b), but somewhat less pronounced before the spring leaf unfolding. Additionally, for forests, the distribution was somewhat prolonged (positively skewed), likely reflecting the importance of higher temperatures during the entire green-up period. A high frequency of significant (at p ≤ 0.01) negative correlation of NDVI and Tmax was visible for grasslands and the Black locust dominated plantation, with peaks in May and August. In grasslands, the critical periods of the Prec and SWC2 (Figure 6c,d) were somewhat similar and very long. The very strong and long-lasting positive correlation between NDVI and Prec (and SWC2) reflects the importance of precipitation for grasslands and their productivity during the summer. For the forests, the negative effects of increased precipitation were present in April. The results for Black locust dominated plantation present an interesting mixture of the results for grasslands and forests for all climate variables. The difference in behavior of the black locust is interesting, especially considering the fact that this species is not native to Europe and was introduced to Hungary in the early 18th century [62]. In the cases of NDVI and Rad, periods with negative correlations were mostly detected, albeit with weaker significance (with p ≤ 0.05 as opposed to p ≤ 0.01; Figure 6e) for forests. Not surprisingly, the most frequent and highly significant (p ≤ 0.01) negative correlation was found in summer for Grassland on saline soil. The exception to this was the green-up period of the forests, when NDVI and Rad showed a positive correlation. For NDVI and VPD (Figure 6f), the critical periods were similar to those obtained with Tmax, though there were no relevant critical periods during the winter (January–February). The critical climate periods for LAI are presented in Figure S10.
The average time lags between the observed change in the NDVI of each ecosystem category and the climate variables used are presented in Table 3 (and for LAI in Table S2), with the highly significant (p ≤ 0.01) positive and negative correlations shown separately. The share of time with significant correlation (STSC) between the NDVI and the climate variable for the period from 25 Jan to 10 Dec (forty 8-day periods; the peak of the winter was discarded to minimize the snow effects) is also presented separately based on the sign of the correlation. The highest STSC was present in the case of the SWC2 and Prec for grasslands, where ~70% of the time the correlation between NDVI and Prec or SWC2 was significantly positive. Interestingly, this was not the case for forest categories, where NDVI and those climate variables appeared to have a more ambiguous relationship and were significantly correlated during a relatively short period (STSC ≤ 25%). The next highest STSC values were present for the herbaceous ecosystem categories and the Black locust dominated plantation, with correlations with both signs of NDVI and Tmax and of NDVI and VPD.

3.4. Relative Effects of the Climate Variables on Each Dominant Ecosystem

While the significant R-values identify periods during the year when the vegetation state was sensitive to the environmental conditions (correlation is significant), they do not provide complete information on the magnitude of NDVI/LAI change caused by the change in the climate variable (assuming all other climate variables are constant). Based on the relative effects (see Section 2.6) of the climate variables on NDVI, there are noticeable differences in the size of the responses of the investigated ecosystem categories (Figure 7). Although in the case of Tmax and VPD, the R-values of the herbaceous ecosystem categories had similar magnitude during late summer (period 30, referring to DOY 233–264 in the state of vegetation), Grassland on saline soil stood out, as it showed notably higher relative change (8.5% vs. −5.8% and −3.6% vs. −2.7%, respectively) for the same unit change in climate variables in comparison to other herbaceous types.
Similarly to the R-curves (Section 3.5), the relative effect curves for LAI are similar to the ones based on NDVI (Figure S11).

3.5. Important Climatic Variables during the Selected Period of the Year

Climatic control of the vegetation state was further analyzed for the selected period (PI) during the late summer (see Figure 4). The Boruta method enabled the identification of the periods when a given climate variable was important for explaining the observed late summer NDVI (Figure 8) and LAI variability (Figure S12). It should be emphasized that in this case, by using Boruta, the effects of all six of the climate variables used were considered (unlike in the case of the relative effect calculation, when only one variable was considered, with all others assumed constant).
Considering the identified important periods for a given climate variable, the results are in very good agreement with the previously presented R-curves (Figures S4–S9). In the case of herbaceous vegetation, generally, the period directly before the investigated late summer period has the highest importance. SWC2 had a determinant effect even in the earlier periods, emphasizing its prominent role in the overall state of the herbaceous vegetation. For forests, earlier periods (i.e., those with longer time lags) might have a stronger effect in determining the state of the vegetation, in line with the observed higher R-values.
Beyond the timing, the results of the Boruta method also provide information about the relative importance of the different climate variables relative to each other for the vegetation state during PI (Figure 8). For example, for the selected late summer period, the important role of SWC2 was clear for almost all the studied ecosystem categories. In the case of the herbaceous vegetation, an additional role of the VPD and/or Tmax was also visible. Contrastingly, for forests, the considered climate variables were less important during the late summer, or, for the category Turkey oak, not determinant at all. For the categories Beech, Turkey oak, and All pedunculate oak, the variable Year also had an important role (with “importance” values of 6, 11.57, and 4.73, respectively), implying the existence of a trend in the NDVI dataset.

3.6. Modeling NDVI and LAI in the Selected Late Summer Period

Based on the LOOY cross-validation, the constructed NDVI models explained 65–87% and 32–40% of the observed NDVI variability for grasslands and forests, respectively (Table 4 and Figure S13). In the case of the modeled LAI, the explained variances of the forests with significant R were higher (33–50%) than for the modeled NDVI. Generally, the category Grassland on saline soil had the highest explained variance (81–87%) in the variability of the NDVI and LAI. For forests, the variability of the NDVI was not captured at a significant p level. For the category All pedunculate oak, the variability was not explained by any of the corresponding models.

4. Discussion

4.1. Methodological Aspects

The present study was based on area-averaged NDVI, LAI, and climate data during 2000–2020 for eight different ecosystem categories, with the aim of gaining insight into the overall behavior of ecosystems in one specific biogeographical region (namely in the relatively small Pannonian BGR). The area-averaging approach used in the study is justified by the fact that the investigated ecosystems were close to natural and grew in areas with similar environmental conditions (Figure S1). The climatic gradients in the study area were relatively small, without any pronounced geographical or high elevation gradients, which additionally supported the approach used in our analysis.
The area-averaging approach proved to be a robust method in the identification of relationships at large spatial scales in our earlier studies [13,30,63]. Area-averaged vegetation indices were used in some studies [18,28], even at a continental scale [64]. The present research aimed to identify periods when a given climate variable most notably affected an ecosystem in general (i.e., we were seeking for robust responses). We did not focus on the possible individual, pixel-level responses where the local micro-conditions could confound the apparent response of the ecosystem. Such confounding effects at fine spatial scales can arise from the issues associated with the driving environmental data (such as the uncertainty of the climate data used, especially in the cases of minimum temperature, convective precipitation, and reanalysis-based SWC); local topography (elevation, aspect, slope, etc.) and soil properties [58,65,66]; management [3,67,68]; species composition [69]; longer-lagged effects [22,23]; recovery of previous disturbances; and other disturbances, such as insects or diseases [31]. The confounding effects of human intervention (management) are also relevant, especially for grasslands. A large fraction of the Hungarian grasslands is affected by regular mowing and/or grazing, which cause a change in the foliar mass. This might result in false attribution of the observed changes to environmental conditions. However, since the management throughout the study period did not change much, the general pattern of the average management effects for the whole area might be similar from year to year. This again speaks in favor of the area-averaging approach.
The presented area-averaged MAM curves (Figure 4) can be considered as “fingerprints” characterizing the investigated ecosystem categories during the 2000–2020 time period and can provide a reference for future studies. MAM curves for a studied ecosystem group were shown in some other studies based on one or more vegetation-related characteristics [28,66], but for more species or ecosystem groups, they are rarely compared. In the case of grasslands, the drop in the value of the vegetation index shown in the MAM curves around DOY 200 might be due to the combined effects of human management (e.g., mowing) and natural processes (e.g., drought).

4.2. Moving-Window Correlation Analysis and Critical Climate Periods

To quantify the relationship between the climate variables and the foliar state or annual productivity, Pearson’s, Spearman’s, partial correlation coefficients [6,14,18,20,26], partial least squares regression [11,12], and forward stepwise regression [24] were used. The studies differ in terms of temporal resolution of the detected critical periods, ranging from daily to seasonal scale, where the lengths of the applied time lags were highly variable [6,18].
The applied moving-window correlation analysis of the present study provided an easy and intuitive tool to quantify the relationship between the vegetation state and the climate variables. The R-curves used were calculated for six climate variables, with a maximum time lag of almost three months, which enabled us to obtain distributions of frequency of significant correlations between NDVI/LAI and selected climate variables (Figure 6 and Figure S10). Such distribution provided a clear overview of the direction (positive/negative), importance (magnitude of the frequency), timing, and duration of critical climate periods for the vegetation state of a given ecosystem category.
Analysis of climatic control based on vegetation indices derived from remote sensing data is common in the literature. Most studies rely only on mean temperature, precipitation, and occasionally radiation [6,11,12,18,28] as the most important climatic drivers of plant growth and senescence. Here, we also used SWC and VPD, since many plant processes are essentially governed by the available water in the soil [25,70,71], and also depend on atmospheric drought [70,72]. SWC is an important environmental variable in ecosystems prone to summer drought, such as those in the Pannonian BGR. We also considered the possibility of the different biophysical effects of Tmin and Tmax, governing different processes during plant growth [13,73,74]. The need for their joint application is emphasized by the fact that seasonal differences exist in the vegetation response to different daytime and night-time warming [75,76].
Based on the presented NDVI-based R-curves (with a time lag of 24 days, Figure 5), the strongest positive correlations were associated with SWC2 and Prec during the late summer for the category of Closed grassland on hard mountainous ground and Grassland on saline soil (see also Table 2). The strongest negative correlation was found also for grasslands between NDVI and Tmax, and also between NDVI and VPD during the summer, implying a strong dependence of NDVI on the antecedent environmental conditions. In contrast, forests showed the strongest relationships during the green-up and the autumn senescence, indicated by a moderately positive correlation with Tmin, Tmax, Rad, and VPD, and by a negative correlation in the case of Prec and SWC2. The latter result is interesting, as it points to the possible role of the increase in soil thermal capacity (due to increase in SWC caused by Prec) along with a likely decrease in the available soil oxygen (due to soil pores being more filled with water) in slowing down forest greening in spring. However, during the summer, some of the forest categories were also affected primarily positively by Prec and SWC2 (Sessile oak with hornbeam and Black locust dominated plantation), although not always at the significance level of p ≤ 0.01. We assume that the weaker R-values of the grasslands in July (Figure 5) might also be attributed to other effects (management and drought consequences), as discussed above (Section 4.1), also affecting the calculated critical climate periods.
The presented critical climate periods were determined based on the short-term response of the foliar state (in contrast to the overall yearly photosynthetic production) using 32-day-long moving-windows. According to the present study, the dominant climate variables from the short-term antecedent period influencing the foliar state of the vegetation state during the summer in the Pannonian BGR were primarily the SWC and the Prec. These results are in accordance with other studies using NDVI, both at the regional scale, such as Zhang et al. [77] for China, or at the global scale [12]. In the case of the grasslands, Tmax and VPD had also a significant role, as demonstrated in previous remote-sensing-based studies, as well for other regions [72,78].
The results reveal a difference in the importance of late-summer meteorology between grasslands and forest ecosystems, which is indicative of the inherent difference in the approach to dealing with drought stress (“strong resilience” approach–grasslands vs. “strong resistance” approach–forests) [79]. However, according to the projected climate changes in the region, an increase in drought frequency/duration and heat stress should be expected, which will lead to an increased frequency of carbon starvation (depletions of the non-structural carbohydrates) and hydraulic failure in trees, thus putting the strong resistance approach in forests to the test.
These findings and the identified critical climate periods are in agreement with previous studies focusing on the productivity or foliar development of temperate vegetation [11,12], emphasizing the existence of climatic key periods in relation to production and also to plant functioning in the diversity of the ecosystem types [80]. Our study does not address the issue of the effects of the extreme years and trends in environmental variables due to climate change. However, the critical periods and the relative effects of the environmental variables on NDVI and LAI identified in this study can provide an insight into the intrinsic traits of the investigated ecosystem categories and streamline efforts in future studies.

4.3. Modeling

The identified critical climate periods and the associated relationships between NDVI (LAI) and environmental drivers inspire the use of simple, multivariate linear models. Process-based models, such as Biome-BGC or its variants [81,82], are excruciatingly complex, frequently calibrated for whole biomes or only for certain ecosystems, and as such are not adequate for studies such as the present one. Therefore, (multiple) linear models, with all of their limitations, are still frequently in use [13,14]. In this study, we focused on the part of the vegetation season (13 August – 13 September) that is a sensitive time interval for Hungarian grassland/herbaceous ecosystems [83], but outside of the periods of leaf-unfolding or senescence, in order to investigate the effects of climate variables on ecosystem categories during periods of fully developed foliage and possible drought.
Models containing a large number of predictors in various periods are neither practical nor explanatory. The selection of the independent variables based exclusively on R-curves might not be justifiable, due to the possible co-linearity between the climate variables. The application of the Boruta method in selecting the important variables enabled the construction of linear regression models with only a few variables that were statistically the most relevant [59]. The timings of the climate variables indicated by the Boruta method corresponded to the periods of the highest R-value of the five R-curves (Figure 5). Building more sophisticated models was out of the scope of the present paper, as we only addressed the causes behind the interannual variability of the vegetation state during the late summer using climate variables.
Based on our results for grasslands, the NDVI models performed the best (with an explained variance of 65–87% of the observed NDVI variability), while for forests, the models constructed to simulate LAI variability were the best (33–50% of explained variance, see Table 4). According to the work of Kern et al. [31], modeling the NDVI of the oak forests in the same late summer period gave similar results (with a maximum R2 of 0.51), but in those models, the NDVI of the preceding period was also used as a proxy variable. Verbesselt et al. [84] also reported similar results (R2 = 0.48) during their model evaluation to forecast insect-induced tree mortality, based on MODIS and long-term daily climate data averages.
Our results corroborated the important effects that environmental conditions have on the state of the grasslands directly before the investigated period. Both for forests and grasslands, the longer-term lagged effects are well-known [3,22,23]; they lead to a delayed reaction and a lag before the change in the foliar state after the change in the environmental conditions. However, in the case of forests, the length of the time lag was found to be shorter when investigated by remote-sensing-based datasets at the ecosystem-level compared to the in situ tree-ring studies, due to the complex physiological processes related to the post-drought upregulation of photosynthesis and repair of canopy damage [85,86]. Despite the well-known existence of the time lag, the investigation of short-term effects (by neglecting the longer-term legacy effects) is common in the literature [6,18,20]. In our study, we used time lags of a maximum of three months; therefore, our models did not account for possible effects with longer time lags, such as the depletion of the non-structural carbohydrate storage [87] or the occurrence of hydraulic failure in trees [88].

4.4. Uncertainty Issues: The Importance of the Land Cover Dataset

Differences between the ecosystem categories regarding their MAM (Figure 4) and their relationships with the antecedent environmental conditions (Figure 5, Figure 6, Figure 7 and Figure 8) confirmed the importance of distinguishing between different ecosystem categories.
Our work was based on pixels with a high share (>90%) of the given ecosystem categories within the pixels, but also within the side-neighboring pixels (>60%), similarly to the criterion applied by e.g., de Beurs and Townsend [89] and Kern et al. [31], which guaranteed that our results were not significantly affected by other land cover types. In addition, we used only those pixels that were free from any significant land cover change according to the NÖSZTÉP and CLC datasets. Furthermore, in the case of forests, only pixels with stable yearly growing season phenological curves were used, similarly to the work of Kern et al. [31]. While the maximum-allowed share of 10% of the other ecosystem categories in a pixel might contribute to uncertainties in our results, its significance was small due to the low mean share of the non-target ecosystem categories within a pixel (<5%, Table 1). The NÖSZTÉP land cover dataset [36] with its fine spatial resolution enabled the unique and accurate categorization of the ecosystem categories in every MODIS pixel.
The differences in the results for different ecosystem categories show not only the diversity in behaviors and responses of different ecosystems, but also emphasize the importance of knowledge of the exact type of vegetation and mixture of species present within an area.

4.5. Uncertainty Issues: Trends

The issue of trends in remotely sensed vegetation indices and productivity during longer timescales is well known [4,27,90,91,92], although its sign and magnitude depend not only on the time of year, ecosystem type, geographical location, pixel size, and investigated period, but also on the applied methodology for detecting long term changes and breakpoints [93,94]. Therefore, trend analysis is a sensitive and challenging scientific topic, where several methods have been elaborated during the last decades [94,95,96,97]. Despite the existing knowledge and ongoing efforts, the trend in the remotely sensed vegetation indices is still a matter of debate and the subject of investigation [10,71]. Compared to the interannual variability, the magnitude of the NDVI trend is generally small (globally 0.0013 yr−1 based on the period 1982–2011 [2], or 0.001 yr−1 based on the period 2001–2015 [98]).
It is worth noting that during the selection of the predictors during the construction of the NDVI model, the input dataset of the Boruta method also contained the calendar Year as a variable, serving as the proxy variable for a trend. The Year variable was identified as an important driver only in the case of forests, but the trends were relatively small. Consequently, even if a trend exists, such a small trend in the dataset does not modify our critical climate period results significantly, but probably leads to a minor amplification or dampening of the signals.

5. Conclusions

In this work, the ecosystem-specific responses of the Pannonian biogeographical region’s vegetation were studied with respect to the climate variables during 2000–2020, based on MODIS, NDVI, and LAI. The presented novel method, using the moving-window correlation analysis and time lags between the NDVI/LAI and climate variables of interest of up to three months, proved to be an easy and intuitive tool supporting the identification, timing and duration of the critical climate periods. Forests showed a pronounced positive correlation with spring and autumn temperatures, but a negative correlation with May precipitation, while for grassland ecosystems, the strongest positive correlations were associated with soil water content and precipitation during the late summer, and negative correlations were associated with radiation.
By using multivariate statistical modeling in combination with the Boruta method for the selection of the most important climate variables, it was possible to estimate the response of NDVI/LAI for different ecosystem categories to the changes in climate variables during late summer. Results reveal a difference in the importance of late summer meteorology between grasslands and forest ecosystems, which is indicative of the inherent difference in the approach to dealing with drought stress.
The present study provides insight into the behavior and key periods relevant to the functioning of ecosystems in the Pannonian biogeographical region. The presented methodology can easily be applied in other geographical regions under markedly different climates. The results could be relevant for areas where increases in the frequency and severity of summer drought and heat are expected due to climate change.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs14215621/s1. Figure S1: Thermopluviogram of the investigated ecosystem categories in Hungary; Figure S2: Geographical location of all grasslands and forests of the studied categories from the NÖSZTÉP ecosystem map dataset in Hungary, based on the native dataset with 20 m × 20 m resolution; Figure S3: Pearson’s R values between the area-averaged LAI and the area-averaged climate variables during 2000–2020; Figure S4: Pearson’s R values between the NDVI and the climate variables (temperatures) for grasslands; Figure S5: Pearson’s R values between the NDVI and the climate variables (temperatures) for woody vegetation (forests); Figure S6: Pearson’s R values between the NDVI and the climate variables (precipitation and SWC) for grasslands; Figure S7: Pearson’s R values between the NDVI and the climate variables (precipitation and SWC) for woody vegetation (forests); Figure S8: Pearson’s R values between the NDVI and the climate variables (radiation and VPD) for grasslands; Figure S9: Pearson’s R values between the NDVI and the climate variables (radiation and VPD) woody vegetation (forests); Figure S10: Critical climate periods within the year, when Tmin, Tmax, Prec, SWC2, Rad, and VPD significantly influence the state of the vegetation (LAI); Figure S11: Relative effects of Tmin, Tmax, Prec, SWC2, Rad, and VPD on the investigated ecosystem groups based on LAI; Figure S12: Relevant climate variables affecting the vegetation state (expressed by LAI) during the PI based on the results of the Boruta algorithm; Table S1: Maximum and minimum significant (p ≤ 0.01) R values representing correlation between LAI and the climate variables used during the year, based on the first R-curve; Table S2: The average length of the time lags (lag, expressed in number of days) between LAI and the climate variables for the time lags with the greatest significant R.

Author Contributions

Conceptualization, A.K., Z.B. and H.M.; methodology, A.K. and H.M.; software, A.K. and R.H.; validation, A.K.; formal analysis, A.K.; investigation, A.K.; resources, A.K., Z.B. and H.M.; data curation, A.K. and E.B.; writing—original draft preparation, A.K., Z.B. and H.M.; writing—review and editing, A.K., Z.B. and H.M.; visualization, A.K.; supervision, A.K. and H.M.; funding acquisition, A.K., Z.B. and H.M. All authors have read and agreed to the published version of the manuscript.

Funding

The research has been supported by the Hungarian Scientific Research Fund (OTKA FK-128709); by the Croatian Science Foundation project MODFLUX (HRZZ IP-2019-04-6325); by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences (grant no. BO/00254/20/10); by the National Multidisciplinary Laboratory for Climate Change, RRF-2.3.1-21-2022-00014 project; by the Hungarian National Research, Development and Innovation Office under Grant number TKP2021-NVA.29; and by the grant “Advanced research supporting the forestry and wood-processing sector’s adaptation to global change and the 4th industrial revolution”, No. CZ.02.1.01/0.0/0.0/16_019/0000803 financed by OP RDE”. The research was prepared with the professional support of the Doctoral Student Scholarship Program of the Co-Operative Doctoral Program of The Ministry of Innovation and Technology, financed from the National Research, Development and Innovation Fund.

Data Availability Statement

Not applicable.

Acknowledgments

We thank NASA for producing and distributing the MODIS products. CORINE datasets were retrieved through the Copernicus land portal on a principle of full, open, and free access, as established by the Copernicus data and information policy Regulation (EU) No 1159/2013 of 12 July 2013. We are grateful to the Lechner Knowledge Center in Hungary for producing and freely disseminating the National Ecosystem Base Map land cover dataset of Hungary, created within the NÖSZTÉP Project. Special thanks to Róbert Lehoczki and Márta Belényesi. We are grateful for the anonymous reviewers for their constructive comments, which helped us to improve the quality of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Anderegg, W.R.L.; Hicke, J.A.; Fisher, R.A.; Allen, C.D.; Aukema, J.; Bentz, B.; Hood, S.; Lichstein, J.W.; Macalady, A.K.; McDowell, N.; et al. Tree mortality from drought, insects, and their interactions in a changing climate. New Phytol. 2015, 208, 674–683. [Google Scholar] [CrossRef] [PubMed]
  2. Huang, K.; Xia, J.; Wang, Y.; Ahlström, A.; Chen, J.; Cook, R.B.; Cui, E.; Fang, Y.; Fisher, J.B.; Huntzinger, D.N.; et al. Enhanced peak growth of global vegetation and its key mechanisms. Nat. Ecol. 2018, 2, 1897–1905. [Google Scholar] [CrossRef] [PubMed]
  3. Catorci, A.; Lulli, R.; Malatesta, L.; Tavoloni, M.; Tardella, F.M. How the interplay between management and interannual climatic variability influences the NDVI variation in a sub-Mediterranean pastoral system: Insight into sustainable grassland use under climate change. Agric. Ecosyst. Environ. 2021, 314, 107372. [Google Scholar] [CrossRef]
  4. Piao, S.; Wang, X.H.; Park, T.; Chen, C.; Lian, X.; He, Y.; Bjerke, J.W.; Chen, A.; Ciais, P.; Tommervik, H.; et al. Characteristics, drivers and feedback of global greening. Nat. Rev. Earth Environ. 2020, 1, 14–27. [Google Scholar] [CrossRef] [Green Version]
  5. IPCC. Climate Change 2022: Impacts, Adaptation, and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change; Pörtner, H.-O., Roberts, D.C., Tignor, M., Poloczanska, E.S., Mintenbeck, K., Alegría, A., Craig, M., Langsdorf, S., Löschke, S., Möller, V., et al., Eds.; Cambridge University Press: Cambridge, UK, 2022; p. 3056. Available online: https://www.ipcc.ch/report/sixth-assessment-report-working-group-ii/ (accessed on 20 August 2022).
  6. Linscheid, N.; Estupinan-Suarez, L.M.; Brenning, A.; Carvalhais, N.; Cremer, F.; Gans, F.; Rammig, A.; Reichstein, M.; Sierra, C.A.; Mahecha, M.D. Towards a global understanding of vegetation–climate dynamics at multiple timescales. Biogeosciences 2020, 17, 945–962. [Google Scholar] [CrossRef]
  7. Seyednasrollah, B.; Young, A.M.; Li, X.; Milliman, T.; Ault, T.; Frolking, S.; Friedl, M.; Richardson, A.D. Sensitivity of deciduous forest phenology to environmental drivers: Implications for climate change impacts across North America. Geophys. Res. Lett. 2020, 47, e2019GL086788. [Google Scholar] [CrossRef]
  8. Zhao, M.; Running, S.W. Drought-Induced Reduction in Global Terrestrial Net Primary Production from 2000 through 2009. Science 2010, 329, 940–943. [Google Scholar] [CrossRef] [Green Version]
  9. Zscheischler, J.; Mahecha, M.D.; von Buttlar, J.; Harmeling, S.; Jung, M.; Rammig, A.; Randerson, J.T.; Schölkopf, B.; Seneviratne, S.I.; Tomelleri, E. A few extreme events dominate global interannual variability in gross primary production. Environ. Res. Lett. 2014, 9, 035001. [Google Scholar] [CrossRef] [Green Version]
  10. Piao, S.; Liu, Q.; Chen, A.; Janssens, I.A.; Fu, Y.; Dai, J.; Liu, L.; Lian, X.; Shen, M.; Zhu, X. Plant phenology and global climate change: Current progresses and challenges. Glob. Change Biol. 2019, 25, 1922–1940. [Google Scholar] [CrossRef]
  11. Guo, L.; Cheng, J.; Luedeling, E.; Koerner, S.E.; He, J.S.; Xu, J.; Gang, C.; Li, W.; Luo, R.; Peng, C. Critical climate periods for grassland productivity on China’s Loess Plateau. Agric. For. Meteorol. 2017, 233, 101–109. [Google Scholar] [CrossRef]
  12. Chen, C.; He, B.; Guo, L.; Zhang, Y.; Xie, X.; Chen, Z. Identifying critical climate periods for vegetation growth in the northern hemisphere. J. Geophys. Res. Biogeosci. 2018, 123, 2541–2552. [Google Scholar] [CrossRef]
  13. Kern, A.; Barcza, Z.; Marjanovic, Z.; Árendás, T.; Fodor, N.; Bónis, P.; Bognár, P.; Lichtenberger, J. Statistical modelling of crop yield in Central Europe using climate data and remote sensing vegetation indices. Agric. For. Meteorol. 2018, 260–261, 300–320. [Google Scholar] [CrossRef]
  14. Vogel, J.; Rivoire, P.; Deidda, C.; Rahimi, L.; Sauter, C.A.; Tschumi, E.; van der Wiel, K.; Zhang, T.; Zscheischler, J. Identifying meteorological drivers of extreme impacts: An application to simulated crop yields. Earth Syst. Dyn. 2021, 12, 151–172. [Google Scholar] [CrossRef]
  15. Hovenden, N.W. Seasonal not annual rainfall determines grassland biomass response to carbon dioxide. Nature 2014, 511, 583–586. [Google Scholar] [CrossRef] [PubMed]
  16. Hoover, D.L.; Rogers, B.M. Not all droughts are created equal: The impacts of interannual drought pattern and magnitude on grassland carbon cycling. Glob. Change Biol. 2016, 22, 1809–1820. [Google Scholar] [CrossRef]
  17. Huang, M.; Wang, X.; Keenan, T.F.; Piao, S. Drought timing influences the legacy of tree growth recovery. Glob. Change Biol. 2018, 24, 3546–3559. [Google Scholar] [CrossRef] [Green Version]
  18. Wang, J.; Rich, P.M.; Price, K.P. Temporal responses of NDVI to precipitation and temperature in the central Great Plains, USA. Int. J. Remote Sens. 2003, 24, 2345–2364. [Google Scholar] [CrossRef]
  19. Delgado-Balbuena, J.; Arredondo, J.T.; Loescher, H.W.; Pineda-Martínez, L.F.; Carbajal, J.N.; Vargas, R. Seasonal precipitation legacy effects determine the carbon balance of a semiarid grassland. J. Geophys. Res. Biogeosciences 2019, 124, 987–1000. [Google Scholar] [CrossRef]
  20. Pei, Z.; Fang, S.; Yang, W.; Wang, L.; Wu, M.; Zhang, Q.; Han, W.; Khoi, D.N. The Relationship between NDVI and Climate Factors at Different Monthly Time Scales: A Case Study of Grasslands in Inner Mongolia, China (1982–2015). Sustainability 2019, 11, 7243. [Google Scholar] [CrossRef] [Green Version]
  21. Anderegg, W.R.L.; Trugman, A.T.; Badgley, G.; Konings, A.G.; Shaw, J. Divergent forest sensitivity to repeated extreme droughts. Nat. Clim. Change 2020, 10, 1091–1095. [Google Scholar] [CrossRef]
  22. Kannenberg, S.A.; Schwalm, C.R.; Anderegg, W.R.L. Ghosts of the past: How drought legacy effects shape forest functioning and carbon cycling. Ecol. Lett. 2020, 23, 891–901. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Bastos, A.; Orth, R.; Reichstein, M.; Ciais, P.; Viovy, N.; Zaehle, S.; Anthoni, P.; Arneth, A.; Gentine, P.; Joetzjer, E.; et al. Vulnerability of European ecosystems to two compound dry and hot summers in 2018 and 2019. Earth Syst. Dyn. 2021, 12, 1015–1035. [Google Scholar] [CrossRef]
  24. Craine, J.M.; Nippert, J.B.; Elmore, A.J.; Skibbe, A.M.; Hutchinson, S.L.; Brunsell, N.A. Timing of climate variability and grassland productivity. Proc. Natl. Acad. Sci. USA 2012, 109, 3401–3405. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Eckes-Shephard, A.H.; Tiavlovsky, E.; Chen, Y.; Fonti, P.; Friend, A.D. Direct response of tree growth to soil water and its implications for terrestrial carbon cycle modelling. Glob. Change Biol. 2021, 27, 121–135. [Google Scholar] [CrossRef]
  26. Teasdale, J.R.; Cavigelli, M.A. Meteorological fluctuations define long-term crop yield patterns in conventional and organic production systems. Sci. Rep. 2017, 7, 688. [Google Scholar] [CrossRef] [Green Version]
  27. Zhou, L.; Tucker, C.J.; Kaufmann, R.K.; Slayback, D.A.; Shabanov, N.V.; Myneni, R.B. Variations in northern vegetation activity inferred from satellite data of vegetation index during 1981 to 1999. J. Geophys. Res. Atmos. 2001, 106, 20269–20283. [Google Scholar] [CrossRef]
  28. Butterfield, Z.; Buermann, W.; Keppel-Aleks, G. Satellite observations reveal seasonal redistribution of northern ecosystem productivity in response to interannual climate variability. Remote Sens. Environ. 2020, 242, 111755. [Google Scholar] [CrossRef]
  29. Zeng, L.; Wardlow, B.D.; Xiang, D.; Hu, S.; Li, D. A review of vegetation phenological metrics extraction using time-series, multispectral satellite data. Remote Sens. Environ. 2020, 237, 11511. [Google Scholar] [CrossRef]
  30. Kern, A.; Marjanović, H.; Barcza, Z. Spring vegetation green-up dynamics in Central Europe based on 20-year long MODIS NDVI data. Agric. For. Meteorol. 2020, 287, 107969. [Google Scholar] [CrossRef]
  31. Kern, A.; Marjanović, H.; Csóka, G.; Móricz, N.; Pernek, M.; Hirka, A.; Matošević, D.; Paulin, M.; Kovač, G. Detecting the Oak lace bug infestation in oak forests using MODIS and meteorological data. Agric. For. Meteorol. 2021, 306, 108436. [Google Scholar] [CrossRef]
  32. Bonan, G.B.; Levis, S.; Kergoat, L.; Oleson, K.W. Landscapes as patches of plant functional types: An integrating concept for climate and ecosystem models. Glob. Biogeochem. Cycles 2002, 16, 5-1–5-23. [Google Scholar] [CrossRef] [Green Version]
  33. Černecký, J.; Gajdoš, P.; Špulerová, J.; Halada, Ľ.; Mederly, P.; Ulrych, L.; Ďuricová, V.; Švajda, J.; Černecká, Ľ.; Andráš, P.; et al. Ecosystems in Slovakia. J. Maps 2019, 16, 28–35. [Google Scholar] [CrossRef] [Green Version]
  34. Vačkář, D.; Grammatikopoulou, I.; Daněk, J.; Lorencová, E. Methodological aspects of ecosystem service valuation at the national level. One Ecosyst. 2018, 3, e25508. [Google Scholar] [CrossRef]
  35. Bardi, A.; Papini, P.; Quaglino, E.; Biondi, E.; Topić, J.; Milović, M.; Pandža, M.; Kaligarič, M.; Oriolo, G.; Roland, V.; et al. Karta prirodnih i poluprirodnih ne-šumskih kopnenih i slatkovodnih staništa Republike Hrvatske. AGRISTUDIO Srl TEMI Srl TIMESIS Srl Hrvat. Agencija Za Okoliš I Prir. Zagreb. 2016, p. 56. Available online: http://www.haop.hr/sites/default/files/uploads/dokumenti/03_prirodne/projekti/NIP-projekt_zavrsno_izvjesce.pdf (accessed on 3 July 2022).
  36. Tanács, E.; Belényesi, M.; Lehoczki, R.; Pataki, R.; Petrik, O.; Standovár, T.; Pásztor, L.; Laborczi, A.; Szatmári, G.; Molnár, Z.; et al. A national, high-resolution ecosystem basemap: Methodology, validation, and possible uses. Term. Közl. 2019, 25, 34–58. [Google Scholar] [CrossRef]
  37. Justice, C.O.; Vermote, E.; Townshend, J.R.G.; Defries, R.; Roy, D.P.; Hall, D.K.; Salomonson, V.V.; Privette, J.L.; Riggs, G.; Strahler, A.; et al. The Moderate Resolution Imaging Spectroradiometer (MODIS): Land remote sensing for global change research. Trans. Geosci. Remote Sens. 1998, 36, 1228–1249. [Google Scholar] [CrossRef] [Green Version]
  38. Mátyás, C.; Berki, I.; Bidló, A.; Csóka, G.; Czimber, K.; Führer, E.; Gálos, B.; Gribovszki, Z.; Illés, G.; Hirka, A.; et al. Sustainability of Forest Cover under Climate Change on the Temperate-Continental Xeric Limits. Forests 2018, 9, 489. [Google Scholar] [CrossRef] [Green Version]
  39. Barcza, Z.; Bondeau, A.; Churkina, G.; Ciais, P.; Czóbel, S.; Gelybó, G.; Grosz, B.; Haszpra, L.; Hidy, D.; Horváth, L.; et al. Modeling of biosphere-atmosphere exchange of greenhouse gases—Model based biospheric greenhouse gas balance of Hungary. In Atmospheric Greenhouse Gases: The Hungarian Perspective; Haszpra, L., Ed.; Springer: Berlin/Heidelberg, Germany, 2010; pp. 295–330. ISBN 978-90-481-9949-5. [Google Scholar] [CrossRef]
  40. EC, European Commission, Directorate-General for Environment; Sundseth, K. Natura 2000 in the Pannonian Region; Publications Office: Luxemburg, Belgium, 2010; ISBN 978-92-79-22586-8. [Google Scholar] [CrossRef]
  41. Lakatos, M.; Bihari, Z.; Izsák, B.; Szentes, O. Globális és hazai éghajlati trendek, szélsőségek változása: 2020-as helyzetkép. (Global trends and climate change in Hungary in 2020). Sci. Et Secur. 2021, 2, 164–171. [Google Scholar] [CrossRef]
  42. Jarvis, A.; Reuter, H.I.; Nelson, A.; Guevara, E. Hole-Filled Seamless SRTM Data V4. International Centre for Tropical Agriculture (CIAT). 2008. Available online: http://srtm.csi.cgiar.org (accessed on 20 May 2021).
  43. LP DAAC (Land Processes Distributed Active Archive Center). MOD09A1, Collection 6. NASA EOSDIS Land Processes DAAC, USGS Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota. Available online: https://lpdaac.usgs.gov (accessed on 14 February 2021).
  44. Vermote, E. MOD09A1 MODIS/Terra Surface Reflectance 8-Day L3 Global 500m SIN Grid V006 [Data Set]; NASA EOSDIS Land Processes DAAC: Reston, VA, USA, 2015. [Google Scholar] [CrossRef]
  45. Viovy, N.; Arino, O.; Belward, A.S. The Best Index Slope Extraction: A method for reducing noise in NDVI time-series. Int. J. Remote Sens. 1992, 13, 1585–1590. [Google Scholar] [CrossRef]
  46. Chen, J.; Jönsson, P.; Tamura, M.; Gu, Z.; Matsushita, B.; Eklundh, L. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. Remote Sens. Environ. 2004, 91, 332–344. [Google Scholar] [CrossRef]
  47. Knyazikhin, Y.; Glassy, J.; Privette, J.L.; Tian, Y.; Lotsch, A.; Zhang, Y.; Wang, Y.; Morisette, J.T.; Votava, P.; Myneni, R.B.; et al. MODIS Leaf Area Index (LAI) and Fraction of Photosynthetically Active Radiation Absorbed by Vegetation (FPAR) Product (MOD15) Algorithm Theoretical Basis Document. 1999. Available online: https://lpdaac.usgs.gov/documents/90/MOD15_ATBD.pdf (accessed on 15 July 2021).
  48. Myneni, R.; Knyazikhin, Y.; Park, T. MOD15A2H MODIS/Terra Leaf Area Index/FPAR 8-Day L4 Global 500m SIN Grid V006 [Data set]; NASA EOSDIS Land Processes DAAC: Reston, VA, USA, 2015. [Google Scholar] [CrossRef]
  49. Dobor, L.; Barcza, Z.; Hlásny, T.; Havasi, Á.; Horváth, F.; Ittzés, P.; Bartholy, J. Bridging the gap between climate models and impact studies: The FORESEE Database. Geosci. Data J. 2015, 2, 1–11. [Google Scholar] [CrossRef]
  50. FORESEE Database. Available online: https://nimbus.elte.hu/FORESEE/ (accessed on 20 May 2021).
  51. CCCS. Copernicus Climate Change Service: ERA5-Land Hourly Data from 2001 to Present [Data Set]; ECMWF: Reading, UK, 2019. [Google Scholar] [CrossRef]
  52. Muñoz-Sabater, J.; Dutra, E.; Agustí-Panareda, A.; Albergel, C.; Arduini, G.; Balsamo, G.; Boussetta, S.; Choulga, M.; Harrigan, S.; Hersbach, H.; et al. ERA5-Land: A state-of-the-art global reanalysis dataset for land applications. Earth Syst. Sci. Data 2021, 13, 4349–4383. [Google Scholar] [CrossRef]
  53. Kern, A.; Marjanović, H.; Barcza, Z. Evaluation of the quality of NDVI3g dataset against Collection 6 MODIS NDVI in Central-Europe between 2000 and 2013. Remote Sens. 2016, 8, 955. [Google Scholar] [CrossRef] [Green Version]
  54. NÖSZTÉP Ecosystem Map of Hungary. Available online: https://alapterkep.termeszetem.hu/ (accessed on 14 November 2020).
  55. Wolfe, R.E.; Nishihama, M.; Fleig, A.J.; Kuyper, J.A.; Roy, D.P.; Storey, J.C.; Patt, F.S. Achieving sub-pixel geolocation accuracy in support of MODIS land science. Remote Sens. Environ. 2002, 83, 31–49. [Google Scholar] [CrossRef]
  56. Tan, B.; Woodcock, C.E.; Hu, J.; Zhang, P.; Ozdogan, M.; Huang, D.; Yang, W.; Knyazikhin, Y.; Myneni, R.B. The impact of gridding artifacts on the local spatial properties of MODIS data: Implications for validation, compositing, and band-to-band registration across resolutions. Remote Sens. Environ. 2006, 105, 98–114. [Google Scholar] [CrossRef]
  57. EEA. 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): Copenhagen, Denmark, 2016; Available online: http://land.copernicus.eu/pan-european/corine-land-cover/clc-2012 (accessed on 24 February 2020).
  58. Zhou, H.; Shao, J.; Liu, H.; Du, Z.; Zhou, L.; Liu, R.; Bernhofer, C.; Grünwald, T.; Dušek, J.; Montagnani, L.; et al. Relative importance of climatic variables, soil properties and plant traits to spatial variability in net CO2 exchange across global forests and grasslands. Agric. For. Meteorol. 2021, 307, 108506. [Google Scholar] [CrossRef]
  59. Kursa, M.B.; Rudnicki, W.R. Feature selection with the Boruta package. J. Stat. Softw. 2010, 36, 1–13. [Google Scholar] [CrossRef] [Green Version]
  60. Santos, J.A.; Malheiro, A.C.; Karremann, M.K.; Pinto, J.G. Statistical modelling of grapevine yield in the Port Wine region under present and future climate conditions. Int. J. Biometeorol. 2011, 55, 119–131. [Google Scholar] [CrossRef]
  61. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2018; Available online: https://www.R-project.org/ (accessed on 9 January 2020).
  62. Keresztesi, B. (Ed.) The Black Locust; Akadémiai Kiadó: Budapest, Hungary, 1988; 197p, ISBN 963-05-4696-5. [Google Scholar]
  63. Kern, A.; Marjanović, H.; Dobor, L.; Anic, M.; Hlásny, T.; Barcza, Z. Identification of Years with Extreme Vegetation State in Central Europe Based on Remote Sensing and Meteorological Data. South East Eur. For. SEEFOR 2017, 8, 1–20. [Google Scholar] [CrossRef] [Green Version]
  64. Zeng, H.; Jia, G.; Epstein, H. Recent changes in phenology over the northern high latitudes detected from multi-satellite data. Environ. Res. Lett. 2011, 6, 045508. [Google Scholar] [CrossRef] [Green Version]
  65. Buttler, A.; Mariotte, P.; Meisser, M.; Guillaume, T.; Signarbieux, C.; Vitra, A.; Preux, S.; Mercier, G.; Quezada, J.; Bragazza, L.; et al. Drought-induced decline of productivity in the dominant grassland species Lolium perenne L. depends on soil type and prevailing climatic conditions. Soil Biol. Biochem. 2019, 132, 47–57. [Google Scholar] [CrossRef]
  66. Kowalski, K.; Okujeni, A.; Brell, M.; Hostert, M. Quantifying drought effects in Central European grasslands through regression-based unmixing of intra-annual Sentinel-2 time series. Remote Sens. Environ. 2022, 268, 112781. [Google Scholar] [CrossRef]
  67. Reinermann, S.; Asam, S.; Kuenzer, C. Remote sensing of grassland production and management—A review. Rem. Sens. 2020, 12, 1949. [Google Scholar] [CrossRef]
  68. Wood, D.J.; Powell, S.L.; Stoy, P.C.; Thurman, L.L.; Beever, E.A. Is the grass always greener? Land surface phenology reveals differences in peak and season-long vegetation productivity responses to climate and management. Ecol. Evol. 2021, 11, 11168–11199. [Google Scholar] [CrossRef] [PubMed]
  69. Hofer, D.; Suter, M.; Haughey, E.; Finn, J.A.; Hoekstra, N.J.; Buchmann, N.; Lüscher, A. Yield of temperate forage grassland species is either largely resistant or resilient to experimental summer drought. J. Appl. Ecol. 2016, 53, 1023–1034. [Google Scholar] [CrossRef] [Green Version]
  70. Zhang, Q.; Ficklin, D.; Manzoni, S.; Wang, L.; Way, D.; Phillips, R.; Novick, K.A. Response of ecosystem intrinsic water use efficiency and gross primary productivity to rising vapor pressure deficit. Environ. Res. Lett. 2019, 14, 074023. [Google Scholar] [CrossRef]
  71. Feng, X.; Fu, B.; Zhang, Y.; Pan, N.; Zeng, Z.; Tian, H.; Lyu, Y.; Chen, Y.; Ciais, P.; Wang, Y.; et al. Recent leveling off of vegetation greenness and primary production reveals the increasing soil water limitations on the greening Earth. Sci. Bull. 2021, 66, 1462–1471. [Google Scholar] [CrossRef]
  72. Yuan, W.; Zheng, Y.; Piao, S.; Ciais, P.; Lombardozzi, D.; Wang, Y.; Ryu, Y.; Chen, G.; Dong, W.; Hu, Z.; et al. Increased atmospheric vapor pressure deficit reduces global vegetation growth. Sci. Adv. 2019, 5, eaax1396. [Google Scholar] [CrossRef] [Green Version]
  73. Peng, S.; Piao, S.; Ciais, P.; Myneni, R.; Chen, A.; Chevallier, F.; Dolman, A.J.; Janssens, I.A.; Peñuelas, J.; Zhang, G.; et al. Asymmetric effects of daytime and night-time warming on northern hemisphere vegetation. Nature 2013, 501, 88. [Google Scholar] [CrossRef]
  74. Xia, J.; Chen, J.; Piao, S.; Ciais, P.; Luo, Y.; Wan, S. Terrestrial carbon cycle affected by non-uniform climate warming. Nat. Geosci. 2014, 7, 173–180. [Google Scholar] [CrossRef]
  75. Tan, J.; Piao, S.; Chen, A.; Zeng, Z.; Ciais, P.; Janssens, I.A.; Mao, J.; Myneni, R.; Peng, S.; Peñuelas, J.; et al. Seasonally different response of photosynthetic activity to daytime and night-time warming in the northern hemisphere. Glob. Change Biol. 2014, 21, 377. [Google Scholar] [CrossRef]
  76. Wu, X.; Liu, H.; Li, X.; Liang, E.; Beck, P.S.A.; Huang, Y.M. Seasonal divergence in the interannual responses of northern hemisphere vegetation activity to variations in diurnal climate. Sci. Rep. 2016, 6, 19000. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  77. Zhang, Q.; Kong, D.; Singh, V.P.; Shi, P. Response of vegetation to different time-scales drought across China: Spatiotemporal patterns, causes and implications. Glob. Planet Change 2017, 152, 1–11. [Google Scholar] [CrossRef] [Green Version]
  78. Ding, J.; Yang, T.; Zhao, Y.; Liu, D.; Wang, X.; Yao, Y.; Peng, S.; Wang, T.; Piao, S. Increasingly important role of atmospheric aridity on Tibetan alpine grasslands. Geophys. Res. Lett. 2018, 45, 2852–2859. [Google Scholar] [CrossRef]
  79. De Keersmaecker, W.; Lhermitte, S.; Tits, L.; Honnay, O.; Somers, B.; Coppin, P. Global vegetation resistance and resilience. Glob. Ecol. Biogeogr. 2015, 24, 539–548. [Google Scholar] [CrossRef]
  80. Gherardi, L.A.; Sala, O.E. Enhanced interannual precipitation variability increases plant functional diversity that in turn ameliorates negative impact on productivity. Ecol. Lett. 2015, 18, 1293–1300. [Google Scholar] [CrossRef]
  81. White, M.A.; Thornton, P.E.; Running, S.W.; Nemani, R.R. Parameterization and sensitivity analysis of the Biome-BGC terrestrial ecosystem model: Net primary production controls. Earth Interact. 2000, 4, 1–85. [Google Scholar] [CrossRef]
  82. Hidy, D.; Barcza, Z.; Marjanovic, H.; Ostrogovic Sever, M.Z.; Dobor, L.; Gelybó, G.; Fodor, N.; Pintér, K.; Churkina, G.; Running, S.; et al. Terrestrial ecosystem process model Biome-BGCMuSo v4.0: Summary of improvements and new modeling possibilities. Geosci. Model Dev. 2016, 9, 4405–4437. [Google Scholar] [CrossRef] [Green Version]
  83. Nagy, Z.; Pintér, K.; Czóbel, S.; Balogh, J.; Horváth, L.; Fóti, S.; Barcza, Z.; Weidinger, T.; Csintalan, Z.; Dinh, N.Q.; et al. The carbon budget of a semi-arid grassland in a wet and a dry year in Hungary. Agric. Ecosyst. Environ. 2007, 121, 21–29. [Google Scholar] [CrossRef]
  84. Verbesselt, J.; Robinson, A.; Stone, C.; Culvenor, D. Forecasting tree mortality using change metrics derived from MODIS satellite data. For. Ecol. Manag. 2009, 258, 1166–1173. [Google Scholar] [CrossRef]
  85. Kannenberg, S.A.; Novick, K.A.; Alexander, M.R.; Maxwell, J.T.; Moore, D.J.; Phillips, R.P.; Anderegg, W.R. Linking drought legacy effects across scales: From leaves to tree rings to ecosystems. Glob. Change Biol. 2019, 25, 2978–2992. [Google Scholar] [CrossRef]
  86. Wong, C.Y.S.; Young, D.J.N.; Latimer, A.M.; Buckley, T.N.; Magney, T.S. Importance of the legacy effect for assessing spatiotemporal correspondence between interannual tree-ring width and remote sensing products in the Sierra Nevada. Remote Sens. Environ. 2021, 265, 112635. [Google Scholar] [CrossRef]
  87. Hartmann, H. Carbon starvation during drought-induced tree mortality—Are we chasing a myth? J. Plant Hydraul. 2015, 2, e005. [Google Scholar] [CrossRef]
  88. Ruehr, N.K.; Grote, R.; Mayr, S.; Arneth, A. Beyond the extreme: Recovery of carbon and water relations in woody plants following heat and drought stress. Tree Physiol. 2019, 39, 1285–1299. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  89. De Beurs, K.M.; Townsend, P.A. Estimating the effect of gypsy moth defoliation using MODIS. Remote Sens. Environ. 2008, 112, 3983–3990. [Google Scholar] [CrossRef]
  90. Nemani, R.R.; Keeling, C.D.; Hashimoto, H.; Jolly, W.M.; Piper, S.C.; Tucker, C.J.; Myneni, R.B.; Running, S.W. Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science 2003, 300, 1560–1563. [Google Scholar] [CrossRef]
  91. Keenan, T.; Gray, J.; Friedl, M.A.; Toomey, M.; Bohrer, G.; Hollinger, D.; Munger, J.; O’Keefe, J.; Peter Schmid, H.; SueWing, I.; et al. Net carbon uptake has increased through warming-induced changes in temperate forest phenology. Nat. Clim. Chang. 2014, 4, 598–604. [Google Scholar] [CrossRef]
  92. Zhu, Z.; Piao, S.; Myneni, R.B.; Huang, M.; Zeng, Z.; Canadell, J.G.; Ciais, P.; Sitch, S.; Friedlingstein, P.; Arneth, A.; et al. Greening of the Earth and its drivers. Nat. Clim. Chang. 2016, 6, 791–795. [Google Scholar] [CrossRef]
  93. De Jong, R.; Verbesselt, J.; Schaepman, M.E.; de Bruin, S. Trend changes in global greening and browning: Contribution of short-term trends to longer-term change. Glob. Chang. Biol. 2012, 18, 642–655. [Google Scholar] [CrossRef]
  94. Forkel, M.; Carvalhais, N.; Verbesselt, J.; Mahecha, M.D.; Neigh, C.S.R.; Reichstein, M. Trend Change Detection in NDVI Time Series: Effects of Inter-Annual Variability and Methodology. Remote Sens. 2013, 5, 2113–2144. [Google Scholar] [CrossRef] [Green Version]
  95. Verbesselt, J.; Hyndman, R.; Newnham, G.; Culvenor, D. Detecting Trend and Seasonal Changes in Satellite Image Time Series. Remote Sens. Environ. 2010, 114, 106–115. [Google Scholar] [CrossRef]
  96. De Jong, R.; de Bruin, S.; de Wit, A.; Schaepman, M.E.; Dent, D.E. Analysis of monotonic greening and browning trends from global NDVI time-series. Remote Sens. Environ. 2011, 115, 692–702. [Google Scholar] [CrossRef] [Green Version]
  97. Ives, A.R.; Zhu, L.; Wang, F.; Zhu, J.; Morrow, C.J.; Radeloff, V.C. Statistical inference for trends in spatiotemporal data. Remote Sens. Environ. 2021, 266, 112678. [Google Scholar] [CrossRef]
  98. Zhang, Y.; Song, C.B.; Band, L.E.; Sun, G.; Li, J. Reanalysis of global terrestrial vegetation trends from MODIS products: Browning or greening? Remote Sens. Environ. 2017, 191, 145–155. [Google Scholar] [CrossRef]
Figure 1. Map of the selected HKM MODIS pixels within Hungary with the dominant NÖSZTÉP ecosystem categories. The map of the used ecosystem categories with all pixels and with the native 20 m × 20 m resolution is presented in Figure S2.
Figure 1. Map of the selected HKM MODIS pixels within Hungary with the dominant NÖSZTÉP ecosystem categories. The map of the used ecosystem categories with all pixels and with the native 20 m × 20 m resolution is presented in Figure S2.
Remotesensing 14 05621 g001
Figure 2. Conceptual figure of illustrative phenological profiles of broad-leaved forests (Sessile oak with hornbeam) and grasslands (Grassland on saline soil), with the indication of the most important phenophases. Phenological profiles are presented for the multiannual means (MAM) and for two selected years.
Figure 2. Conceptual figure of illustrative phenological profiles of broad-leaved forests (Sessile oak with hornbeam) and grasslands (Grassland on saline soil), with the indication of the most important phenophases. Phenological profiles are presented for the multiannual means (MAM) and for two selected years.
Remotesensing 14 05621 g002
Figure 3. Methodology of the critical climate period analysis based on moving-window correlation. An example of the five R-curves used based on five different time lags is presented on the left. The different intervals of the climate variables relative to the interval of the vegetation state and the corresponding time lags between them are shown in the top right corner. One period is an 8-day long interval.
Figure 3. Methodology of the critical climate period analysis based on moving-window correlation. An example of the five R-curves used based on five different time lags is presented on the left. The different intervals of the climate variables relative to the interval of the vegetation state and the corresponding time lags between them are shown in the top right corner. One period is an 8-day long interval.
Remotesensing 14 05621 g003
Figure 4. Multiannual mean (MAM) and interannual variability (IAV) curves of each ecosystem category based on NDVI and LAI, derived from the selected pixels. The range of the ±1.96 SEMAM based on the yearly country averaged curves are also shown for the MAM curves. Pink rectangles/columns indicate the period of interest (PI) in late summer (DOY 225–256, 13 August–13 September).
Figure 4. Multiannual mean (MAM) and interannual variability (IAV) curves of each ecosystem category based on NDVI and LAI, derived from the selected pixels. The range of the ±1.96 SEMAM based on the yearly country averaged curves are also shown for the MAM curves. Pink rectangles/columns indicate the period of interest (PI) in late summer (DOY 225–256, 13 August–13 September).
Remotesensing 14 05621 g004
Figure 5. Pearson’s R-values between the area-averaged NDVI (after the indicated period) and the area-averaged climate variables (before the indicated period) separately (af). Both the NDVI and the climate variables refer to 32-day-long intervals, with an 8-day overlap. Pink rectangles indicate the period of interest (PI) in late summer.
Figure 5. Pearson’s R-values between the area-averaged NDVI (after the indicated period) and the area-averaged climate variables (before the indicated period) separately (af). Both the NDVI and the climate variables refer to 32-day-long intervals, with an 8-day overlap. Pink rectangles indicate the period of interest (PI) in late summer.
Remotesensing 14 05621 g005
Figure 6. Critical climate periods within the year, when Tmin (a), Tmax (b), Prec (c), SWC2 (d), Rad (e), and VPD (f) significantly influenced the state of the vegetation (NDVI). The significant correlation frequency (i.e., the numbers of the instances, maximum 20; five time lags & four 8-day periods; see methods) associated with significant R is shown for all 8-day periods during the year. Periods corresponding to positive R-values are indicated with green, while those with negative R-values are indicated with red color. Significant values at p ≤ 0.01 are shown in darker shades, while those with p ≤ 0.05 are shown in lighter shades.
Figure 6. Critical climate periods within the year, when Tmin (a), Tmax (b), Prec (c), SWC2 (d), Rad (e), and VPD (f) significantly influenced the state of the vegetation (NDVI). The significant correlation frequency (i.e., the numbers of the instances, maximum 20; five time lags & four 8-day periods; see methods) associated with significant R is shown for all 8-day periods during the year. Periods corresponding to positive R-values are indicated with green, while those with negative R-values are indicated with red color. Significant values at p ≤ 0.01 are shown in darker shades, while those with p ≤ 0.05 are shown in lighter shades.
Remotesensing 14 05621 g006
Figure 7. Relative effects of Tmin (a), Tmax (b), Prec (c), SWC2 (d), Rad (e) and VPD (f) for the investigated ecosystem groups based on NDVI and the first R-curves with a time lag of 24 days. It expresses what percentage change the NDVI showed for one reference unit change of the climate variable, where the reference units were 1 °C, 1 °C, 5 mm, 0.01 m3 m−3, 10 W m−2, and 50 Pa in the case of Tmin, Tmax, Prec, SWC2, Rad, and VPD, respectively. Pink rectangles indicate the period of interest (PI) in late summer.
Figure 7. Relative effects of Tmin (a), Tmax (b), Prec (c), SWC2 (d), Rad (e) and VPD (f) for the investigated ecosystem groups based on NDVI and the first R-curves with a time lag of 24 days. It expresses what percentage change the NDVI showed for one reference unit change of the climate variable, where the reference units were 1 °C, 1 °C, 5 mm, 0.01 m3 m−3, 10 W m−2, and 50 Pa in the case of Tmin, Tmax, Prec, SWC2, Rad, and VPD, respectively. Pink rectangles indicate the period of interest (PI) in late summer.
Remotesensing 14 05621 g007
Figure 8. Relevant climate variables affecting the vegetation state (described by NDVI) during the 13 August–13 September period based on the results of the Boruta method. The colors express the “importance” value as the result of the Boruta feature selection. The * and the following colored rectangles indicate that in the case of the given ecosystem category, the Year variable was also found to be relevant, indicating the existence of a trend.
Figure 8. Relevant climate variables affecting the vegetation state (described by NDVI) during the 13 August–13 September period based on the results of the Boruta method. The colors express the “importance” value as the result of the Boruta feature selection. The * and the following colored rectangles indicate that in the case of the given ecosystem category, the Year variable was also found to be relevant, indicating the existence of a trend.
Remotesensing 14 05621 g008
Table 1. Spatial, areal, and altitudinal descriptive statistics of the selected HKM MODIS pixels regarding the different ecosystem categories of the NÖSZTÉP dataset. The “total native share” gives information about the total presence of the given ecosystem category within Hungary based on the native resolution of the LC dataset (at 20 × 20 m).
Table 1. Spatial, areal, and altitudinal descriptive statistics of the selected HKM MODIS pixels regarding the different ecosystem categories of the NÖSZTÉP dataset. The “total native share” gives information about the total presence of the given ecosystem category within Hungary based on the native resolution of the LC dataset (at 20 × 20 m).
Ecosystem CategoriesNumber of the PixelsMean of the
Pixel-Level
Share of the
Given Category within the Pixels
Share in Hungary Based on the Used Pixels at 500 m (and all Possible at 20 m)Altitudinal
Distribution
(Median,
5 and 95
Percentiles)
(1) Open sandy grassland1199.4%0.01% (0.68%)117 m [106–144 m]
(2) Grassland on saline soil108699.4%1.01% (2.27%)86 m [83–93 m]
(3) Closed grassland on hard mountainous ground 11798.9%0.11% (4.91%)219 m [81–286 m]
(4) Beech19298.5%0.18% (1.49%)586 m [262–879 m]
(5) Sessile oak with hornbeams21098.6%0.19% (1.74%)400 m [242–530 m]
(6) Turkey oak11597.6%0.11% (2.83%)293 m [216–435 m]
(7) Black locust dominated plantation4297.7%0.04% (4.87%)139 m [112–210 m]
(8) All pedunculate oak (merged category)1595.3%0.01% (1.67%)115 m [86–154 m]
Table 2. Maximum and minimum significant (p ≤ 0.01) R-values representing the correlation between NDVI and the climate variables used during the year, based on the first R-curve (see also Figure 5). The months associated with the maximum and minimum significant R-values are also indicated, referring to the 32-day long period of the investigated vegetation state. The R-value with the greatest absolute value of each climate variable is marked in bold. In the case of no significant correlation (p > 0.01), the R-values are not indicated.
Table 2. Maximum and minimum significant (p ≤ 0.01) R-values representing the correlation between NDVI and the climate variables used during the year, based on the first R-curve (see also Figure 5). The months associated with the maximum and minimum significant R-values are also indicated, referring to the 32-day long period of the investigated vegetation state. The R-value with the greatest absolute value of each climate variable is marked in bold. In the case of no significant correlation (p > 0.01), the R-values are not indicated.
Ecosystem CategoriesTminTmaxPrecipitationSWC2RadiationVPD
RMonthRMonthRMonthRMonthRMonthRMonth
Positive correlation with NDVI
Open sandy grassland0.68III--0.76IX0.83VII–VIII----
Grassland on saline soil0.72III0.71II–III0.83VIII-IX0.88VII--0.61II–III
Closed grassland on hard m. ground0.76II–III0.72II–III0.85IX0.91IX–X--0.55III–IV
Beech0.80X0.81X----0.67IV–V0.75IV–V
Sessile oak with hornbeam0.78X0.81IV–V0.58X–XI--0.62IV–V0.84IV–V
Turkey oak0.71X0.81IV–V----0.65IV–V0.81IV–V
Black locust dominated plantation0.78II–III0.73II–III0.68IX0.66IX–X--0.61III
All pedunculate oak0.78III–IV0.72III------0.67IV–V
Negative correlation with NDVI
Open sandy grassland--−0.75VIII–IX----−0.66VI–VII−0.83VIII–IX
Grassland on saline soil--−0.80VIII–IX----−0.69VIII–IX−0.83VIII–IX
Closed grassland on hard m. ground--−0.71V–VI----−0.59VIII–IX−0.73V–VI
Beech----−0.77V−0.72V–VI----
Sessile oak with hornbeam----−0.75V−0.60V–VI−0.68XII--
Turkey oak----−0.76V−0.60V–VI----
Black locust dominated plantation--------−0.61II–III−0.63IX
All pedunculate oak------------
Table 3. The average length of the time lags (lag, expressed in a number of days) between NDVI and the used climate variables for the time lags with the most significant R and the share of time between 25 Jan and 10 Dec (forty 8-day periods; the peak of the winter was discarded to minimize the snow effects) when the correlation was highly significant (p ≤ 0.01; STSC–share of time with significant correlation), according to the sign of the correlation.
Table 3. The average length of the time lags (lag, expressed in a number of days) between NDVI and the used climate variables for the time lags with the most significant R and the share of time between 25 Jan and 10 Dec (forty 8-day periods; the peak of the winter was discarded to minimize the snow effects) when the correlation was highly significant (p ≤ 0.01; STSC–share of time with significant correlation), according to the sign of the correlation.
Ecosystem CategoriesCorrelationTminTmaxPrecipitationSWC2RadiationVPD
Lag
(Days)
STSCLag
(Days)
STSCLag
(Days)
STSCLag
(Days)
STSCLag
(Days)
STSCLag
(Days)
STSC
Open sandy grasslandnegative-0%42.532%-0%-0%34.334%38.449%
positive28.428%-0%26.955%24.068%-0%-0%
Grassland on saline soilnegative-0%36.040%-0%-0%41.330%27.453%
positive30.413%28.010%36.070%24.673%-0%32.05%
Closed grassland on hard m. groundnegative-0%37.330%-0%-0%40.738%12.635%
positive26.020%30.413%37.770%26.970%-0%24.03%
Beechnegative-0%-0%44.825%42.318%80.05%-0%
positive34.745%35.040%-0%-0%32.010%35.738%
Sessile oak with hornbeamnegative-0%-0%29.315%24.05%40.018%-0%
positive29.828%28.933%40.08%50.78%32.05%24.023%
Turkey oaknegative-0%-0%46.020%24.03%88.03%-0%
positive28.825%25.823%-0%40.03%52.010%42.433%
Black locust dominated plantationnegative82.78%80.015%-0%-0%53.718%76.825%
positive30.153%43.623%33.118%38.425%-0%32.05%
All pedunculate oaknegative-0%-0%725%5613%64.010%-0%
positive4040%25.823%-0%-0%-0%24.013%
Table 4. The performance metrics of the cross-validation modeling NDVI and LAI in the selected late summer period (13 August–13 September, period 29–32).
Table 4. The performance metrics of the cross-validation modeling NDVI and LAI in the selected late summer period (13 August–13 September, period 29–32).
NDVILAI
Ecosystem CategoriesR2pRMSEbiasR2pRMSEbias
Open sandy grassland0.65 *0.0000.046−0.0010.58 *0.0000.1880.000
Grassland on saline soil0.87 *0.0000.036−0.0010.81 *0.0000.154−0.002
Closed grassland on hard m. g.0.76 *0.0000.0370.0000.77*0.0000.1510.000
Beech0.36 *0.0040.0080.000#
Sessile oak with hornbeam0.100.1610.015−0.0010.43 *0.0010.1770.006
Turkey oak0.40 *0.0020.0130.0000.50 *0.0000.176−0.001
Black locust dominated p.0.32 *0.0070.032−0.0000.33 *0.0060.4700.006
All pedunculate oak0.020.5280.0150.001#
* Significant correlation coefficient (at p ≤ 0.01) of the multiple linear regression model. # In the case of LAI for category Beech and All pedunculate oak, the Boruta method did not identify any variable as important.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Kern, A.; Barcza, Z.; Hollós, R.; Birinyi, E.; Marjanović, H. Critical Climate Periods Explain a Large Fraction of the Observed Variability in Vegetation State. Remote Sens. 2022, 14, 5621. https://doi.org/10.3390/rs14215621

AMA Style

Kern A, Barcza Z, Hollós R, Birinyi E, Marjanović H. Critical Climate Periods Explain a Large Fraction of the Observed Variability in Vegetation State. Remote Sensing. 2022; 14(21):5621. https://doi.org/10.3390/rs14215621

Chicago/Turabian Style

Kern, Anikó, Zoltán Barcza, Roland Hollós, Edina Birinyi, and Hrvoje Marjanović. 2022. "Critical Climate Periods Explain a Large Fraction of the Observed Variability in Vegetation State" Remote Sensing 14, no. 21: 5621. https://doi.org/10.3390/rs14215621

APA Style

Kern, A., Barcza, Z., Hollós, R., Birinyi, E., & Marjanović, H. (2022). Critical Climate Periods Explain a Large Fraction of the Observed Variability in Vegetation State. Remote Sensing, 14(21), 5621. https://doi.org/10.3390/rs14215621

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop