Spatial and Seasonal Patterns in Vegetation Growth-Limiting Factors over Europe

: Water and energy are recognized as the most inﬂuential climatic vegetation growth-limiting factors. These factors are usually measured from ground meteorological stations. However, since both vary in space, time, and scale, they can be assessed by satellite-derived biophysical indicators. Energy, represented by land surface temperature (LST), is assumed to resemble air temperature; and water availability, related to precipitation, is represented by the normalized di ﬀ erence vegetation index (NDVI). It is hypothesized that positive correlations between LST and NDVI indicate energy-limited conditions, while negative correlations indicate water-limited conditions. The current project aimed to quantify the spatial and seasonal (spring and summer) distributions of LST–NDVI relations over Europe, using long-term (2000–2017) MODIS images. Overlaying the LST–NDVI relations on the European biome map revealed that relations between LST and NDVI were highly diverse among the various biomes and throughout the entire study period (March–August). During the spring season (March–May), 80% of the European domain, across all biomes, showed the dominance of signiﬁcant positive relations. However, during the summer season (June–August), most of the biomes—except the northern ones—turned to negative correlation. This study demonstrates that the drought / vegetation / stress spectral indices, based on the prevalent hypothesis of an inverse LST–NDVI correlation, are spatially and temporally dependent. These negative correlations are not valid in regions where energy is the limiting factor (e.g., in the drier regions in the southern and eastern extents of the domain) or during speciﬁc periods of the year (e.g., the spring season). Consequently, it is essential to re-examine this assumption and restrict applications of such an approach only to areas and periods in which negative correlations are observed. Predicted climate change will lead to an increase in temperature in the coming decades (i.e., increased LST), as well as a complex pattern of precipitation changes (i.e., changes of NDVI). Thus shifts in plant species locations are expected to cause a redistribution of biomes.

biomes except for the evergreen forests, where the correlation switched to positive. Based on the same relations, Nemani and Running [29] proposed a simple land-cover classification scheme of various land-cover types across the continental USA. Tateishi and Ebata [30] used the LST-NDVI relation to differentiate between phenology patterns on a global scale. They showed a negative LST-NDVI correlation in most of the equatorial latitudes and a positive one in the high altitudes. Julien et al. [31] used LST and NDVI values to estimate changes in vegetation in the European continent between 1982 and 1999. They were able to distinguish between trends from the arid areas in southern Europe and those from the moister northern regions. Over Mongolia, Karnieli et al. [32] demonstrated that the slope of the LST-NDVI correlation changed across biomes and latitudes. Significant negative correlations were found in the desert and desert steppe regions, while significant positive correlations were shown in the high mountains and taiga.
Longer temporal scales-months and seasons-were also examined. From the intra-annual perspective, LST-NDVI relations were found to be seasonally dependent. Over the USA, Sun and Kafatos [33] revealed positive correlations during winter and negative correlations during summer. This general pattern can be modified in shorter periods (months) and for smaller spatial land-cover units [15]. Moreover, significant changes in the monthly correlations were identified in conjunction with land covers. Karnieli et al. [34] investigated a long LST and NDVI time series over North America. They revealed positive correlations from April-May, over most of the domain except for the southern latitudes. Specifically, the evergreen needleleaf forest biome deviated from these positive relations. In contrast, the LST-NDVI slopes gradually shifted to negative values in mid-summer, while cropland, grassland, and shrubland exhibited significant negative correlations in June and July. Similar seasonal dependencies were observed around Berlin, Germany [35]. Six different land-cover/land-use categories were defined and analyzed for all four seasons. In autumn, winter, and spring seasons all LST-NDVI correlations were positive but statistically non-significant. During summer, all correlations were negative, but again generally not significant, except in green urban areas and forests where they were found to be significant.
Based on the above literature review, one can infer that energy and water control vegetation growth, but they are mostly time-and region-dependent. Earth observation datasets of biophysical satellite-derived products such as LST and NDVI enable the modeling of vegetation productivity with the ability to cover large areas worldwide. It is hypothesized that positive LST-NDVI correlations indicate energy-limited conditions, while negative correlations appear for water-limited conditions. Consequently, the direction and magnitude of the correlation may change seasonally. The overarching goal of the work presented here was to quantify the spatial and seasonal (spring and summer) distributions of LST-NDVI relations as main indicators for vegetation growth-limiting factors (water and energy constraints) in the European domain. Specifically, based on long-term (2000-2017) LST and NDVI Earth observation data, the study objectives were twofold: (1) to explore long-term, pan-European correlations between LST and NDVI during the spring and summer seasons, based on 18 years of Earth observation data; and (2) to analyze LST-NDVI patterns for each European biome at the seasonal scale.

Study Area
The study area was the European continent, as identified by its standard geographic definition. It is bordered by the Arctic Ocean in the north, the Atlantic Ocean in the west, and the Mediterranean Sea in the south. In the east, Europe is bordered by the watershed that divides the Ural Mountains and the Ural River, the Caspian and Black Seas. The territory of Turkey and some regions south of the Caucasus Mountains were included in the analysis, but Cyprus was excluded. Additionally, Iceland, Severny-Yuzhny, and other smaller northern islands were excluded since they are almost entirely covered by glaciers with very little vegetation cover variability.
The climate of Europe varies with respect to latitude, longitude, and altitude. Solar radiation is a function of latitude. Consequently, the European climate changes from the hot, dry Mediterranean summer climate in the south, through the colder and wetter middle region with no dry season, to polar climate in the north. Under the influence of the Gulf Stream, Western Europe has mild winters and relatively higher annual temperatures. Parts of Eastern Europe are characterized by dry and even arid climate, affected by continental air masses that originate in Asia. The climatic gradient also exists at the main mountainous regions such as the Alps, the Carpathians, and the Caucasus mountains. All these factors determine European biomes.
The study domain included seven terrestrial biomes; here, we use the biome definition provided by the World Wide Fund for Nature (WWF) [36]. The biomes were categorized based on their climatic and geologic features, as well as their evolutionary history, gridded to a 0.05-degree spatial resolution, and their locations within the study area are shown in Figure 1. The biomes are as follows: (1) Tundra, the most northern Arctic belt, is characterized by low temperatures and short vegetation growing season in the summer. The plants mainly consist of dwarf shrubs and lichens. (2) Boreal forest/taiga is a subarctic biome where the forests are mainly evergreen coniferous species with deciduous trees dominating in the oceanic climate. (3) The temperate forest is the largest biome in Europe, spreading from west to east across the domain. It is characterized by cool winters and warm moist summers, mainly with deciduous trees. (4) Temperate conifer forest is a cold and windy biome that is mainly confined to high elevated ridges such as the Alps, Carpathians, and Pontic Mountains, as well as in Scotland. (5) Temperate steppe covers large parts of Eastern Europe, north of the Black Sea and the Caspian Sea, as well as south of the Caucasus Mountains. This biome is characterized by a warm, dry climate, and since it receives a small amount of rain, shrubs are common and trees are sparse. (6) Mediterranean forest, woodland, and scrub, located in southern Europe, is generally characterized by hot dry summers and mild to cool, rainy winters. (7) Desert and xeric shrubland is located in the east of the domain, near the Caspian Sea, north and south of the Caucasus Mountains. This biome is characterized by less than 250 mm of annual rainfall and high evaporation rates, with dispersed woody plants. The biome map was downloaded from the WWF Terrestrial Ecoregions of the World website (https://www.worldwildlife.org/publications/terrestrial-ecoregions-of-the-world).

MODIS Data Sources and Analyses
The daytime LST and NDVI were products of MODIS Collection 6 datasets. Data were extracted from the NASA/USGS Land Resources Distributed Active Archive Center (LP DAAC) (https://lpdaac. usgs.gov/) as MOD13C2 and MYD13C2 data products for NDVI (from Terra and Aqua satellites, respectively) and MOD11C3 and MYD11C3 for LST. The data were downloaded at a spatial resolution of 0.05 degrees, for the 18 years between 2000 and 2017. The monthly products included images between March and August. This seasonal window allows examination of the main European vegetation growing season while avoiding the massive winter snow cover. The daytime LST values were converted from kelvin to degrees Celsius. Both NDVI and LST were averaged for three periods: the entire period (March to August) and two sub-periods, representing springtime (March-May) and summer (June-August). The March-August period was selected since it represents a large amplitude of vegetation phenology, in most of the European domain [37]. Each season included 18 years of LST and NDVI images, cropped to fit the study area. Pixels that included no data, or were covered with clouds, were removed. Moreover, water bodies and urban areas were masked out of the images, according to land-cover maps produced by the European Environmental Agency's COoRdinate INformation on the Environment (CORINE) in 2012 (https://land.copernicus.eu/pan-european/corine-land-cover). characterized by hot dry summers and mild to cool, rainy winters. (7) Desert and xeric shrubland is located in the east of the domain, near the Caspian Sea, north and south of the Caucasus Mountains. This biome is characterized by less than 250 mm of annual rainfall and high evaporation rates, with dispersed woody plants. The biome map was downloaded from the WWF Terrestrial Ecoregions of the World website (https://www.worldwildlife.org/publications/terrestrial-ecoregions-of-theworld).

LST-NDVI Relation Analyses at Pan-European and Biome Scales
The pre-processed LST and NDVI products were analyzed to quantify the relationships between these factors for the entire European domain, and at the biome scale. The relationships between LST and NDVI were analyzed by plotting the averaged values of these factors, for each defined timeframe using the 18 images against each other, for each of the timeframes specified in Section 2.2. Scatterplots were computed for the entire European domain, one for each timeframe. Similarly, biome-specific scatterplots were generated.
The correlation between the two variables throughout the examined period was quantified using the Pearson correlation (r lst-ndvi ) for each individual pixel. For each pixel, values of LST and NDVI were examined for each three months for 18 years, giving a total of 54 observations per pixel. For the entire period of six months, the total number of observations was 108. The resulting r map illustrates the strength and sign of the LST-NDVI relations. This analysis was conducted for each of the three periods specified in Section 2.2.
For each period, the significance of the LST-NDVI correlation for each pixel was also computed using the t-statistic. Significant (α = 0.05) positive correlation was defined when r > 0 and t > 1.96 and significant negative correlation when r < 0 and t < −1.96, according to the t-distribution. Correlations where −1.96 < t < 1.96 were defined as non-significant. Each pixel was then assigned to its biome category. General quantification of the percentage of negatively, positively, or non-significantly correlated pixels within each biome for each period was generated to estimate the overall LST-NDVI trend distribution across the biomes (positive, negative, or non-significant) for each timeframe.

LST-NDVI Correlation Dynamics
In order to explore long-term pan-European correlations between LST and NDVI, their values were firstly mapped over the spatial domain of interest. Long-term (2000-2017) averages of higher daytime LST and NDVI for the entire seasonal study period (March-August) and the two sub-seasons  Figure 2A,E,I generally indicate higher mean daytime LST values along the Mediterranean region in southern Europe that gradually decreased toward the higher latitudes. From June-August, the mean temperature range was 0.5 to 55 • C (with a mean of 26 • C), compared to a range of −12 to 37.5 • C (with a mean of 12.4 • C) from March-May. NDVI ( Figure 2B,F,J) was characterized by overall higher values from June-August (range = −0.18 to 0.91, mean = 0.66) than from March-May (range = −0.19 to 0.81, mean = 0.42), with a substantial increase in northern and eastern Europe. Negative NDVI values indicate snow-covered pixels. Pairs of LST and NDVI maps were used to create density scatterplots ( Figure 2C,G,K). In these plots, red and yellow colors represent higher pixel concentrations. For the entire period ( Figure 2C), two distinct trends can be identified: one with a negative slope and a second with a positive slope. The strengths of these opposite correlations changed with the sub-period considered. A dominant positive correlation was observed during the spring months ( Figure 2G). In contrast, during the summer months, most of the domain was characterized by a negative correlation ( Figure 2K). The spatial distribution of r lst-ndvi also varied between seasons ( Figure 2D,H,L). For the entire period, most of the European region was characterized by a significant positive LST-NDVI relationship, with negatively-correlated areas only in the Mediterranean region and the Caucasus Mountains. From March-May, the positive relationship was even more pronounced across the continent. The warmer June-August sub-season reflected a different spatial pattern: the southern half of the study region showed mostly negative r lst-ndvi values. However, mountainous areas (i.e., the Alps and the Carpathian Mountains) and the northern latitudes showed non-significant r lst-ndvi correlation, with patches of positive values.

Biome-Specific LST-NDVI Relations
As noted above, seven biomes were involved in the current study: tundra; boreal forest/taiga; temperate forest; temperate conifer forest; temperate grassland, savanna, and shrubland; Mediterranean forest, woodland, and scrub; and desert and xeric shrubland (Figure 1). The biome-specific relations between LST and NDVI were highly diverse among the various biomes and throughout the entire study period (Figure 3). The northern biomes with the lowest temperature values (Tundra and Taiga) showed a constant positive relationship between LST and NDVI. The temperate biomes (broadleaf and mixed forests, conifer forests, and grasslands) were characterized by shifting trends between seasons (positive relationships from March-May that turned negative from June-August). The biomes with higher temperatures (Mediterranean forest and desert) typically had a negative LST-NDVI relationship throughout the entire studied timeframe. Table 1 shows r lst-ndvi relations as percent of pixels that displayed each trend type (positive, negative, and non-significant). Overall, during the entire period (March-August), most of the biomes-65% of the study area-displayed a significant positive relationship. The Mediterranean and desert biomes, with a negative correlation, covered only 19% of the area ( Figure 2D). During the spring season (March-May), 80% of the European domain-across all biomes-showed the dominance of significant positive relations. Exceptions were the Mediterranean forest and desert biomes, where the positive correlation covered about half of the area with relatively high coverage of non-significant pixels ( Figure 2H). However, the summer season (June-August) featured a shift when 56% of the European area showed negative correlation values, 19% had positive values, and 25% of the pixels were non-significant ( Figure 2L). Tundra and boreal forest/taiga were the only biomes that demonstrated a majority of significantly positive LST-NDVI relationships (57% and 60% of their total areas, respectively). The rest of the biomes were characterized by a large areal coverage of negative r lst-ndvi values, except the temperate conifer forest biome which typically showed a larger portion of non-significant r lst-ndvi values (45% of its area).

Discussion
This project explored the spatial and seasonal distributions of LST-NDVI relations as indicators of VGLFs, while distinguishing between water and energy limitation factors at the pan-European scale. Further, the project intended to link the LST-NDVI relations and the European biomes. To accomplish this goal, long-term (2000-2017) MODIS-derived LST and NDVI variables were obtained. Since the study covered a rather large continental scale, positive and negative correlations between LST and NDVI were observed in different regions of the domain and at different sub-seasonal periods. From a detailed analysis of these relationships, it can be inferred that when water is the limiting factor for vegetation growth during the summer season in the study area's low latitudes, the LST-NDVI correlations are typically negative. However, when energy is the limiting factor for vegetation growth, as is the case at the northern latitudes or high elevations areas (especially at the beginning of the growing season), LST and NDVI correlations are positive. The correlation signs-either negative, positive, or non-significant-varied between different seasons and biomes. For the spring months, March-May, most of the areas of all biomes showed a positive correlation. During the summer months, June-August, most of the areas of the biomes changed to negative correlations, but the northern biomes of tundra and taiga mostly continued to display positive correlations. Large areas of non-significant correlations were also observed. These findings led to the following discussions.

Energy Balance Perspective
The spatial distribution of LST-NDVI relationships can also be examined from the perspective of the energy balance [38,39]. In rural areas, vegetation generally has a lower albedo than non-vegetated terrain in the same region, due to its spectral characteristics (vegetation cover) and its 3D structure. This 3D structure due to tree canopies leads to multiple scattering, increasing the shortwave radiation. Thus, longwave radiation is trapped, resulting in increased net longwave radiation [40]. Therefore, because of their lower albedo, the vegetated areas store more energy in the form of net radiation than the surrounding non-vegetated areas [38]. Part of this energy is stored in the surface, increasing its temperature (LST). This leads to a temperature gradient between the surface and lower atmosphere, and part of this energy is used to evaporate water, increasing turbulent latent heat flux. Evaporative cooling depends on the available energy [41], the landscape heterogeneity [42], as well as the type of vegetation and the soil moisture conditions. While evaporative cooling over grassland might exceed that over forest at times of ample soil moisture, the reverse is likely to occur under conditions of low soil moisture [43][44][45].
Generally, in southern latitudes (especially in arid/semi-arid areas) and the mid-latitudes during the summer period, solar radiation is high, leading to higher air temperature. During the daytime, high LST (and therefore high saturation water vapor pressure) is observed, whereas atmospheric water vapor pressure is usually low. Therefore, a vegetated surface in these cases is more likely to use a great part of the incoming energy to evaporate water. The higher the water availability (e.g., precipitation, irrigation), the higher the latent heat flux. Since the dominant type of vegetation is broadleaf, stomatal resistance is low, further contributing to high latent heat flux. Therefore, in these cases, the incoming energy in vegetated areas is mainly used to evaporate water, despite the lower albedo of vegetated areas. Conversely, the incoming energy in non-vegetated areas is stored in the respective surface, leading to the increase of temperature. For this reason, the LST in non-vegetated areas is usually higher than the LST of the surrounding vegetated areas. Pixels corresponding to vegetated areas (high NDVI) are expected to have a lower LST than the pixels that correspond to non-vegetated areas (low NDVI), which leads to a revised LST-NDVI correlation ( Figure 2D,H,L).
Conversely, in northern latitudes, both the solar radiation and the air temperature are generally low. In general, where low LST (and therefore lower stomatal resistance) is observed, the relative humidity is consequently higher, and the air is close to saturation. Furthermore, in northern latitude areas, where boreal forest/taiga and tundra dominate, the respective stomatal resistance is expectedly high. Therefore, a relatively low latent heat flux is expected for vegetated areas in these regions. The net incoming energy in the form of radiation is stored in the vegetation canopies, increasing the respective LST. A small part of this energy is transported to the atmosphere as sensible heat flux due to the low-temperature gradient between the surface and air temperature in these latitudes. The mechanism is the same for the surrounding non-vegetated areas (no latent heat flux and low sensible heat flux). Thus, the incoming energy is stored in the surfaces, with higher storage in vegetated areas because the net incoming energy is higher due to the lower albedo. The higher heat storage in vegetated areas results in a higher LST than in non-vegetated ones. Consequently, pixels corresponding to vegetated areas (high NDVI) are expected to have a higher LST than the pixels corresponding to non-vegetated areas (low NDVI), which leads to a positive LST-NDVI correlation ( Figure 2D,H,L). Finally, in the mid-latitudes and transition zones, as well as regions that are characterized by mixed-land-cover types, a combination of the above cases is expected, where both seasonality and water availability effects are most prominent.

European Biomes and Vegetation Growth-Limiting Factors
Two concerns were raised in this work regarding climate change and land-use changes: how will biome extents vary and how will the LST-NDVI relationships be affected. Parmesan [46] reviewed the ecological effects of climate change on phenology and species range shifts, and Watson et al. [47] estimated the vulnerability of ecoregions under climate change. In turn, Eigenbrod et al. [48] showed that the vulnerability of ecosystems to climate changes that induce biome shift depends on habitat intactness and on the scale on which intactness is estimated. In a study of European vegetation, Thuiller et al. [49] revealed that species distribution shifts triggered by climate change threaten plant biodiversity. In particular, species loss and turnover were found to be mainly determined by just two climatic variables: temperature and precipitation. In our approach, these are represented by LST and NDVI, respectively. For the coming decades, climate projections generally indicate an increase of temperature (thus, of LST) and a complex pattern of precipitation changes (thus, of NDVI). Therefore, plant species are expected to shift within and between biomes.
Since the relative strength of the different VGLFs might change, a major issue regards estimating future vegetation stress. Summer droughts in particular are expected to increase in the coming decades in Europe [50], possibly modifying the role of the different VGLFs and strengthening water limitation (and thus negative LST-NDVI relations). In turn, this is expected to lead to more severe summer fires and a drastic increase in the summer burned area [51]. Such estimates, based on empirical models, are rooted in the effects of the different VGLFs and can thus be affected by changes in the role of the different limiting factors, as expressed, for example, by the LST-NDVI correlations studied here. During the past years, for example, unexpected summer fires plagued northern European regions. A deeper understanding and modeling of the response of VGLFs and of the vegetation itself to the changing climate conditions is thus needed.

Relation to LST-NDVI Models
Previous studies, mostly limited in area and observation period, revealed a strong negative correlation between LST and NDVI. Based on this prevailing assumption, several vegetation/drought/stress indices were developed and widely applied. The hypothesis behind this assumption was the synergetic effect of water availability and energy effect on vegetation: less water and higher temperature promote vegetation stress and vice versa [26]. Hence, this naive model used only the ratio between LST and NDVI, acquired from the NOAA/AVHRR imagery, for mapping drought over Papua New Guinea and Mongolia [52,53]. A widely-used index based on the negative LST-NDVI relation is the Vegetation Health Index (VHI), developed by Kogan [54] as a proxy characterization of vegetation health or a combined estimation of moisture and thermal conditions. The index originally used long-term data from the NOAA-AVHRR spacecraft and has recently been used with the Visible Infrared Imaging Radiometer Suite (VIIRS) spectral bands [55] as well. VHI has been applied on a regional scale in many areas around the world [56][57][58], but has also been proposed as a global drought index [55,[59][60][61]. Global coverage of the VHI is updated weekly on an official NOAA webpage (https://www.star.nesdis.noaa.gov/smcd/emb/vci/VH/vhbrowse.php).
Another spectral index that relies on an empirical parameterization of negative relations between LST and NDVI is the Temperature Vegetation Dryness Index (TVDI) [62,63]. This index is based on the location of a pixel in the LST-NDVI triangle or trapezoid-shaped scatterplot that associates higher LST with less evaporation, while higher NDVI values are identified with higher vegetation cover. This index has been applied to a vast number of sites all over the world, including Europe, as listed in Petropoulos et al. [64]. Wang and Dickinson [5] have reservations about using TVDI in northern latitudes and cold areas where temperatures, as a major control of latent heat flux, are mostly positively related to evaporation. Garcia et al. [65] examined the performance of the TVDI in a spatially heterogeneous region in southern Spain using MODIS imagery. They revealed that the index successfully assessed the water pressure deficit only under water-limited conditions and failed under energy-limited conditions.
Since droughts occur mostly in southern latitudes, these indices were developed with the implicit assumption that NDVI and LST are always and everywhere negatively correlated. Due to their simplicity, such as LST-NDVI-based models, they are widely accepted and commonly used. Furthermore, the availability of long-term remote sensing images at different spatial and temporal resolutions enabled the easy application of these models. In areas where water is ultimately the limiting factor for vegetation growth throughout the year, this assumption is apparently correct. However, these indices have been applied globally, even in regions and periods where the correlation is positive or non-significant.

Conclusions
In the current work, LST is a proxy for near-surface air temperature and NDVI is related to precipitation. The analysis of the correlations between these two biophysical variables provides information on how they vary in space and time and how they are related to each other. When energy is the limiting factor for vegetation growth, as in the case of northern latitudes and high elevations, significant positive correlations exist. Under water availability constraints, on the other hand, significant negative correlations are observed. During spring months, most of the domain is characterized by significant positive correlations that become significantly negative toward the summer, except for the northernmost regions. This behavior can be explained by taking into account the energy partitioning between latent and sensible heat flux, mainly regulated by surface albedo and temperature, as well as by air temperature.
The spatial distribution of the LST and NDVI correlation was found to be in line with previous studies that addressed large-scale climatic constraints but primarily relied on ground-based measurements and models. Long-term and global coverage of Earth-observed imagery from space provides a reliable dataset for biome mapping and ecological models. Predicted climate change will lead to an increase of temperature in the coming decades (i.e., increased LST), as well as a complex pattern of precipitation changes (i.e., changes of NDVI), thus shifts in plant species locations are expected causing redistribution of biomes.
The current project demonstrates that the drought/vegetation/stress indices, based on the prevalent hypothesis of a negative correlation between LST and NDVI, are also spatially and temporally dependent. Consequently, they are not valid in regions where energy is the limiting factor or during specific periods of the year. Furthermore, their results are doubtful when the correlations are statistically non-significant. Consequently, it is essential to re-examine this assumption and restrict applications of such an approach only to areas and periods in which negative correlations are observed.