Agricultural Drought Detection with MODIS Based Vegetation Health Indices in Southeast Germany

Abstract

Droughts during the growing season are projected to increase in frequency and severity in Central Europe in the future. Thus, area-wide monitoring of agricultural drought in this region is becoming more and more important. In this context, it is essential to know where and when vegetation growth is primarily water-limited and whether remote sensing-based drought indices can detect agricultural drought in these areas. To answer these questions, we conducted a correlation analysis between the Normalized Difference Vegetation Index (NDVI) and Land Surface Temperature (LST) within the growing season from 2001 to 2020 in Bavaria (Germany) and investigated the relationship with land cover and altitude. In the second step, we applied the drought indices Temperature Condition Index (TCI), Vegetation Condition Index (VCI), and Vegetation Health Index (VHI) to primarily water-limited areas and evaluated them with soil moisture and agricultural yield anomalies. We found that, especially in the summer months (July and August), on agricultural land and grassland and below 800 m, NDVI and LST are negatively correlated and thus, water is the primary limiting factor for vegetation growth here. Within these areas and periods, the TCI and VHI correlate strongly with soil moisture and agricultural yield anomalies, suggesting that both indices have the potential to detect agricultural drought in Bavaria.

1. Introduction

Drought is a complex, globally occurring phenomenon that affects humans and nature alike [1]. It can be examined from different perspectives: meteorological (precipitation deficit), agricultural (soil moisture deficit), hydrological (runoff and water storage deficit), and socio-economic drought (consideration of water supply and demand) [2,3,4]. Drought events in Europe have far-reaching impacts and economic costs within different sectors of society [5,6]. Within the European Union and United Kingdom, the annual economic losses (1981–2010) from drought are estimated at EUR 9 billion year⁻¹ [7]. In Central Europe, and thus also in Bavaria, extreme drought events in recent years have repeatedly resulted in massive damage and losses in agriculture and forestry [8,9,10,11].

Since the beginning of the 21st century, Central Europe has experienced recurrent periods of exceptional drought, with the years 2003, 2015, and 2018, in particular, being extremely hot and/or dry [12,13,14,15]. While climate data from the past do not yet allow a consistent statement on drought trends in Central Europe [16,17,18,19], future projections mainly assume an increase in drought frequency and severity during the growing season [20,21,22,23]. In addition, there could be a significant increase in the risk of consecutive droughts and their affected areas in the future [15].

In this context, the detection and forecast of drought are becoming more and more important, a major challenge being the monitoring of the impacts of agricultural drought on vegetation. In addition to in-situ measurements, remote sensing based observations allow a broad scale detection and can identify and fill relevant knowledge gaps within this thematic area. A variety of approaches based on remote sensing exist to detect drought impacts on vegetation [24,25,26], some of which have already been applied in Central Europe, with different methodologies. The majority of existing studies made use of remotely sensed vegetation indices (VIs).

One possibility is the comparison of VIs with meteorological drought indices. It has been shown that the Normalized Difference Vegetation Index (NDVI) and Vegetation Optical Depth (VOD) could capture agricultural drought events in response to meteorological droughts within Europe [27]. In addition, it has been found that meteorological drought indices for shorter periods were associated with crop stress (defined via the Vegetation Condition Index (VCI) and Vegetation Health Index (VHI)), while longer accumulation periods correlated better with the vegetation status of forest areas. At the same time, the magnitude of regional differences in drought impacts within Europe was pointed out [28]. Equally, when the NDVI and Standardized Precipitation Index (SPI) data in spring and summer were compared to each other, high positive correlations were detected, especially in Eastern Europe and on the Iberian Peninsula [29].

Another approach within VI use is the analysis and comparison with soil moisture data. Comparing the NDVI with soil moisture indices such as the Palmer Drought Severity Index, the Self-Calibrating Palmer Drought Severity Index and the Normalized total depth Soil Moisture (NSM), clear differences in correlations have been observed across Europe [29]. Here, the NSM achieved the highest correlations (also in Eastern Europe and on the Iberian Peninsula), whereby this was explained primarily with the additional information included in NSM compared to the other indices about soil moisture in deeper soil layers. Using the NSM as a soil moisture variable, it was also found that in European areas with a warmer summer climate, the NDVI only responded to fluctuations in soil moisture, but not in temperature. In spring, this dependency was reversed. In colder regions, the NDVI was only dependent on temperature in both seasons. Furthermore, during the 2018 agricultural drought in the Netherlands, negative soil moisture anomalies occurred 2–3 weeks before the first VI reduction, using near-infrared reflectance of terrestrial vegetation and VOD as VI [30]. When comparing NDVI and Climate Change Initiative soil moisture data, an offset between soil moisture drought and vegetation drought in Europe was also found [31].

In Central Europe, analyses of VIs have been used to assess the impact of agricultural droughts on corn and winter wheat yields. Corn yields at the NUTS3 level in Germany correlated well with VCI and VHI in August, although clear spatial correlation differences could also be seen [28]. On the other hand, comparing Temperature Condition Index (TCI), VCI and VHI with winter wheat yields in northern and eastern Germany, VHI, in particular, was able to achieve higher correlations than meteorological indices [32]. Comparing the VIs NDVI and Enhanced Vegetation Index (EVI) with climatic datasets, it was shown that during the extreme drought in July 2018, 1.5 times more area within Europe was negatively impacted than during the extreme drought in August 2003. Differences within land-use types were also found [33].

Within remote sensing of agricultural drought, there are also studies which developed combined drought monitoring in Central Europe. Sepulcre-Canto et al. [34], for example, applied a combination of SPI, soil moisture anomaly, and the fraction of Absorbed Photosynthetically Active Radiation (fAPAR), while Trnka et al. [35] used a combination of soil moisture data (remotely sensed, modelled, and reported) and EVI in the Czech Republic and Slovakia. Both approaches were shown to be promising. Finally, approaches also exist within this topic to determine areas and time periods where water is the primary limiting factor for plant growth. In Europe, a correlation analysis between LST and NDVI has revealed that energy was a limiting factor mainly in northern Europe, at high altitudes, and in spring, while water limited plant growth mainly in southern Europe and in summer [36].

Despite the large numbers of existing studies on agricultural and vegetation-based drought in Central Europe using remote sensing, uncertainties and larger knowledge gaps still exist. For example, the spatial as well as temporal resolution of results indicating where and when water is the limiting factor for vegetation growth is still rather limited. There are also hardly any evaluations that deal with other explanatory environmental variables for vegetation stress caused by drought. In Bavaria in particular, there are still hardly any approaches for continuous drought monitoring via remote sensing-based VIs that can demonstrate an appropriate spatial and temporal resolution. Therefore, our study aims, for the first time, to provide results with a much better spatial and temporal resolution than currently exists as to where and when in Central Europe water is the limiting factor for vegetation growth. Second, the analysis of this limitation incorporates other environmental conditions to explain it, which has rarely been done before. Finally, an approach not yet practiced in Bavaria is presented to perform continuous drought monitoring via remote sensing-based VIs at an appropriate spatial and temporal resolution.

Addressing the aforementioned gaps in knowledge, the general objective of this study is to detect agricultural drought using remote sensing in Central Europe on water-limited areas for vegetation growth. Specifically, this means, with a spatial resolution of up to 250 m and a temporal resolution of up to 8 days:

• To determine, by means of a correlation analysis between Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI and LST, areas and time periods in which water is the limiting factor for vegetation growth. In addition, it will be examined whether and to what extent the factors of land cover and altitude influence these conditions;
• To carry out drought monitoring using the TCI, VCI, and VHI, as well as to evaluate their results by soil moisture and agricultural yields. The question to be addressed here is whether and to what extent these indices can be used to detect agricultural drought.

2. Materials and Methods

The general work-flow of this study is shown in Figure 1, including preprocessing of data, the classified correlation analysis between NDVI and LST, and the calculation and evaluation of the drought indices. All relevant analyses were carried out in R (version 4.0.3). More details on the following steps are provided in the respective subsections.

2.1. Study Area

The study area of Bavaria is located between 47°N and 50.5°N, and between 9°E and 14°E, in the southeastern part of Germany, Central Europe (Figure 2). The climate of the region is mainly influenced by strong relief, with higher elevations in the south (northern edge of the Alps) and east (Bavarian Forest and Fichtel Mountains; Figure S1 in Supplementary Materials). The mean annual temperature ranges from −3.3 to 11 °C, but in the majority of the territory, it is between 8 and 10 °C (Figure S2, based on data from the German Meteorological Service, 1991–2020). The mean annual precipitation sums range from 515 to 3184 mm, with wetter conditions in the southern part of Bavaria (Figure S3). The land cover in Bavaria is largely dominated by agriculture (54.58%) and forest (35.38%). The agricultural areas are mainly found in the northwest and southeast of Bavaria, while forest cover dominates towards the Alps and in the east of Bavaria. Open grasslands and larger water areas are mostly localized in the Alpine region and Alpine foothills. Munich and Nuremberg constitute the largest urbanized areas (Figure 2, based on Corine Land Cover (CLC) 2018 data). Bavaria is divided into 96 counties (based on the yield data for agricultural crops/Statistical Offices of the Federation and the states).

2.2. Data

MODIS NDVI, LST, and Land Cover (LC) products were used for both the correlation analysis between NDVI and LST and the calculation of drought indices. The chosen period for all MODIS products was from 2001 to 2020. The NDVI raster data (MOD13Q1) have a spatial resolution of 250 m and a temporal resolution of 16 days [37]. The LST raster data (MOD11A2) provide Daytime Land Surface Temperature [K] with a spatial resolution of 1 km and a temporal resolution of 8 days [38]. The LC product from MODIS (MCD12Q1) is issued annually with a spatial resolution of 500 m. We used the International Geosphere-Biosphere Programme (IGBP) classification (LC Type 1) [39].

Both the LST and LC data were resampled first to a spatial resolution of 250 m in order to be compatible with the NDVI data. In addition, an annual vegetation mask (Table S2) was created using the MODIS LC data and applied to the NDVI and LST data to preemptively avoid misinterpretation, especially with respect to NDVI. Finally, the 16-day NDVI grid data were linearly interpolated to an 8-day temporal resolution in order to perform both the correlation analysis between NDVI and LST and the calculation of drought indices in 8-day time steps.

In addition, we used the land cover product CLC 2018 and the European Digital Elevation Model (EU-DEM) v1.1 for a better and more differentiated evaluation of the correlation analysis between NDVI and LST. The CLC and EU-DEM datasets are produced within the framework of the EU Copernicus Service coordinated by the European Environment Agency (EEA) and have a spatial resolution of 100 [40] and 25 m [41], respectively. To combine the two datasets, the elevation model was aggregated to a spatial resolution of 100 m using the arithmetic mean.

To evaluate the calculated drought indices VCI, TCI, and VHI (see Section 2.4), soil moisture index data from 2001 to 2020, as well as yield data for agricultural crops in combination with land use data were used. The present work exploits the Soil Moisture Index (SMI) published by the Helmholtz Centre for Environmental Research (UFZ). SMI represents the percentile of soil moisture content in the topsoil up to a depth of 25 cm, relative to its frequency distribution in the period 1951–2019. The data have a monthly temporal resolution and a spatial resolution of 4 km [42,43]. In order to compare the data with the calculated raster data of the drought indices, the soil moisture index data were bilinearly resampled to a spatial resolution of 250 m.

The yield data for agricultural crops from 2001 to 2019 are published in the Regional Database Germany, which is operated by the Statistical Offices of the Federation and the states. The data are given as the annual mean yield per hectare in dt ha⁻¹ and are provided at county level [44]. The land use data for the agricultural holdings in Bavaria originate from the Bavarian State Office for Statistics and refer to the reference year 2016. The cultivated area of the respective field crop is also given for each county [45]. In the present paper, the area ratios of the individual cultivated crops within the counties are assumed to be relatively stable and thus refer to the observation period of the yield statistics (2001 to 2019). For the selection of yield data to be included in the analysis, data on annual harvest dates in Bavaria from 2001 to 2019 were used based on selected phenological observations of the respective crops (Phase ID 24—“Ernte”) from the German Meteorological Service (DWD) [46].

2.3. Correlation Analysis between NDVI and LST

The first step of the present study is a correlation analysis between NDVI and LST to determine when and where water is the limiting factor for vegetation growth. NDVI is one of the most widely used vegetation indices in remote sensing and is defined as follows [47,48]:

$$NDVI=\frac{{\rho }_{NIR}-{\rho }_R}{{\rho }_{NIR}+{\rho }_R},$$

(1)

where ${\rho }_{NIR}$ is the reflectance in the near-infrared band and ${\rho }_R$ is the reflectance in the red band. NDVI allows one to distinguish between healthy and stressed vegetation and to determine the growth status of vegetation, and thus a measure of general vegetation health [26,49,50]. In addition, NDVI can also be an effective indicator of the vegetation-moisture condition [51]. LST, in turn, mirrors soil moisture conditions, evapotranspiration, and vegetation water stress. Because evapotranspiration cools the earth’s surface, rising surface temperatures may be linked to decreasing water content or a soil moisture deficit and thus water stress in the vegetation canopy [52,53]. If both variables are well correlated, the limiting factor for vegetation growth can be determined depending on the season and location: The NDVI-LST correlation is negative if water is the limiting factor, but it is positive if energy is the limiting factor [36,54,55].

Here the correlation analysis between NDVI and LST was performed at pixel level, at a monthly time scale, following the method of Abdi et al. [56]. The Bravais-Pearson correlation coefficient was calculated for each pixel and month in the vegetation period from March to October using all available MODIS data from 2001 to 2020 (

Table 1. Monthly subdivision of Moderate Resolution Imaging Spectroradiometer (MODIS) day numbers for pixel-based monthly correlation analysis between NDVI and LST.

Month	MODIS Day Numbers
March	065–089
April	097–113
May	121–145
June	153–177
July	185–209
August	217–241
September	249–273
October	281–297

).

The calculated correlation coefficients were additionally analyzed by land cover, altitude, and their combination. For this, the monthly NDVI-LST correlation maps were resampled to 100 m spatial resolution. Subsequently, both the land cover data (CLC 2018) and the digital elevation model data (EU-DEM v1.1) were classified into different classes (

Table 2. Land cover (CLC 2018) and elevation (EU-DEM v1.1) classes and their IDs applied to the monthly NDVI-LST correlation maps.

ID	Land Cover Classes (CLC 2018)	ID	Altitude Classes (EU-DEM v1.1)
NAL	Non-irrigated arable land	AL1	<300 m
PAS	Pastures	AL2	300–500 m
BLF	Broad-leaved forest	AL3	500–800 m
CFF	Coniferous forest	AL4	800–1200 m
MXF	Mixed forest	AL5	>1200 m
NGL	Natural grasslands

). The selected land cover classes had to meet two criteria: (1) contain vegetation; and (2) cover at least 1% of the Bavarian land area. The altitude classification followed relevant papers with a reasonable number of classes [57,58,59,60]. The respective independent masks were then intersected with the monthly NDVI-LST correlation maps and finally characterized using boxplots.

When combining both classifications, the land cover masks were intersected with the altitude masks. All 30 (6 land cover classes × 5 altitude classes) combined masks were then intersected with the monthly NDVI-LST correlation maps, for which the median was determined for each class.

2.4. Calculation and Evaluation of the Drought Indices VCI, TCI, and VHI

The VHI consists of two components, i.e., VCI and TCI. The VCI is defined as follows [61,62,63]:

$$VCI=100\ast \frac{\left(NDVI-NDV{I}_{min}\right)}{(NDV{I}_{max}-NDV{I}_{min})}.$$

(2)

The NDVI is the smoothed 8-day NDVI, the NDVI_min, and the NDVI_max the corresponding multiyear absolute maximum and minimum. The VCI is accordingly scaled from 0 to 100 and takes on larger values the higher the NDVI is. The TCI is defined as follows [61,62,63]:

$$TCI=100\ast \frac{(T_{max}-T)}{\left(T_{max}-T_{min}\right)},$$

(3)

where T is the smoothed 8-day temperature, T_min and the T_max are the corresponding multiyear absolute maximum and minimum. The TCI is also scaled accordingly from 0 to 100 but takes on higher values the lower T is. In this work, T is defined via the LST. The combination of both described indices results in the VHI [61,62,64]:

$$VHI=\mathsf{\alpha }\ast VCI+\left(1-\mathsf{\alpha }\right)\ast TCI,$$

(4)

where α determines the share of VCI and TCI in the VHI. Since this share depending on location and time is unknown, the weight of both indices was assumed to be equal (α = 0.5) [64,65,66]. The value range of the VHI is also from 0 (severe vegetation stress) to 100 (very favorable conditions). The VHI is based on the assumption that there is a negative correlation between NDVI and LST: During a drought period, the NDVI tends to be low, while the LST tends to be high [54].

In this work, all three drought indices were calculated pixel-based for Bavaria in 8-day resolution within the growing season (March–October) from 2001 to 2020. Due to the VHI assumption of a negative correlation between NDVI and LST, it is important to exclude all areas with a non-negative NDVI-LST correlation in order to avoid misinterpretations. For this purpose, a threshold value of −0.1 was applied to all monthly NDVI-LST correlation maps, and all areas with a correlation coefficient >−0.1 were then removed from the respective drought index time series on a monthly basis.

For the evaluation of the drought indices time series, TCI, VCI, and VHI were compared with the SMI and the yield data. For the comparison with the monthly SMI, we calculated the monthly arithmetic means of the 8-days drought indices, based on Table 1. Subsequently, the Bravais-Pearson correlation coefficients were determined between the SMI and the drought index data-sets for the complete time series (monthly and total vegetation period). Agricultural yield data comprised winter wheat, winter rye, summer barley, oat, sugar beet, winter rapeseed, and silage corn, i.e., major crops for which harvest dates in Bavaria were available and which were on average harvested at the end of July at the earliest (Day of year > 200;

Table 3. List of selected crop yield data from the Regional Database Germany and the corresponding mean harvest dates (as the day of the year) from 2001 to 2019 in Bavaria in the evaluation study.

Selected Crops in The Yield Statistics (Regional Database Germany)	Mean Harvest Date (Day of The Year) in Bavaria from 2001 to 2019
Winter wheat	215
Winter rye	212
Summer barley	211
Oat	221
Sugar beet	283
Winter rapeseed	205
Silage corn	266

). The reason for this threshold value was the inclusion of at least one month (July) in the evaluation, in which water can be assumed as the limiting factor for vegetation growth over large areas in Bavaria.

To be compared with crop yield data, drought indices were averaged pixel-wise and annually for a crop-specific period of influence, i.e., all months from May to July for winter wheat, winter rye, summer barley, oat, and winter rapeseed or May to September for sugar beet and silage corn. Then, mean drought indices were determined for each county. A correlation analysis was then performed between the mean annual drought indices and the respective relative annual yield anomaly at county level in Bavaria. These correlations were accordingly also determined within the counties. In order to avoid misinterpretations, counties, where, for a specific cultivated crop, the area share of the total area of the seven selected crops was less than 3%, were excluded from the spatial correlation analysis.

In order to summarize the results across the individual crops, an additional analysis was carried out based on weighted relative crop yield anomalies. For this purpose, the share of the total area of the seven selected crops was determined for each county and each crop. Subsequently, the respective annual crop yield anomalies were weighted with the area shares for each county and each year and summed up. In any case, years in which data were not available for all crops grown in that county were removed to ensure comparability. Mean annual drought indices (May to July or May to September) were also weighted at the county level according to the corresponding area shares of the crop and added up to obtain a weighted drought index. Finally, an overall correlation analysis, as well as a spatial correlation analysis differentiated by county between the relatively weighted crop yield anomalies and the weighted drought indices, were carried out.

3. Results

3.1. Correlations between NDVI and LST

As evidenced by the monthly correlation analysis between NDVI and LST for the period 2001 to 2020, NDVI-LST relationships clearly depend on several factors. The first factor is the time of the year. Across Bavaria (Figure 3 and Table S3), positive NDVI-LST correlations occurred above all both at the beginning (March to May) and at the end (September and October) of the vegetation period. In contrast, in the middle of the vegetation period or summer season (June to August), negative correlations prevailed. Regardless of the season, areas in the Alps, the Bavarian Forest, and the Spessart region showed consistently positive correlations, while negative correlations prevailed for the main part of the vegetation period within the Main-Franconian Plates or the Isar-Inn Hills.

This seasonal influence on the sign of the correlation between NDVI and LST was also confirmed when being classified by land cover, altitude, and their combination. Similarly, the correlation coefficients were higher at the beginning and end of the growing season than in the summer months of June, July, and August. Negative NDVI-LST correlations occurred especially in these summer months, while in the marginal months of March, April, and October only positive correlations were observed in almost all classes (Figure 4 and Figure 5;

Table 4. Median pixel-based monthly correlation coefficients (Bravais-Pearson) between NDVI and LST in Bavaria from 2001 to 2020, by altitude and land cover classes. The percentage (in parentheses) of the respective altitude class indicates the area share of this class in the total area of Bavaria. The percentage of the respective land cover class indicates the area share of this land cover in the total area of the respective altitude class. The color-coding of the individual median values is based on Figure 3. Land cover classes according to CLC 2018 are non-irrigated arable land (NAL), pastures (PAS), broad-leaved forest (BLF), coniferous forest (CFF), mixed forest (MXF), and natural grasslands (NGL). Altitude levels are <300 m (AL1), 300–500 m (AL2), 500–800 m (AL3), 800–1200 m (AL4), and >1200 m (AL5; Table 2).

Altitude	CLC 2018	Number of Pixels	March	April	May	June	July	August	September	October
AL1	NAL	240,901 (48.69%)	0.35	0.00	−0.09	−0.23	−0.35	−0.33	−0.19	0.01
(7.01%)	PAS	77,745 (15.71%)	0.46	0.19	0.01	−0.12	−0.29	−0.36	−0.14	0.12
	NGL	1650 (0.33%)	0.46	0.09	−0.01	−0.08	−0.43	−0.41	−0.19	0.15
	BLF	38,559 (7.79%)	0.42	0.32	0.16	0.08	−0.06	−0.06	0.12	0.18
	CFF	24,399 (4.93%)	0.31	0.31	0.16	0.14	−0.03	−0.05	0.13	0.12
	MXF	22,339 (4.51%)	0.36	0.38	0.17	0.14	−0.03	−0.02	0.14	0.18
AL2	NAL	1,603,322 (43.47%)	0.44	0.07	−0.03	−0.12	−0.26	−0.29	−0.08	0.14
(52.28%)	PAS	609,244 (16.52%)	0.55	0.16	−0.01	−0.09	−0.26	−0.33	−0.01	0.18
	NGL	13,274 (0.36%)	0.44	0.19	0.09	−0.03	−0.23	−0.29	0.08	0.16
	BLF	253,590 (6.88%)	0.45	0.35	0.19	0.14	−0.01	0.04	0.20	0.18
	CFF	627,277 (17.01%)	0.40	0.22	0.15	0.10	−0.05	−0.07	0.10	0.12
	MXF	239,313 (6.49%)	0.44	0.31	0.17	0.13	−0.02	−0.01	0.16	0.17
AL3	NAL	511,891 (22.01%)	0.51	0.09	0.00	−0.07	−0.24	−0.21	−0.02	0.14
(32.96%)	PAS	646,484 (27.80%)	0.64	0.18	0.02	−0.06	−0.22	−0.16	0.10	0.20
	NGL	8681 (0.37%)	0.61	0.24	0.11	0.00	−0.18	−0.15	0.10	0.17
	BLF	92,655 (3.98%)	0.53	0.34	0.19	0.12	0.00	0.05	0.20	0.15
	CFF	595,578 (25.61%)	0.53	0.24	0.16	0.10	−0.02	0.03	0.11	0.13
	MXF	254,112 (10.93%)	0.53	0.30	0.18	0.09	−0.01	0.06	0.17	0.13
AL4	NAL	0 (0.00%)	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.
(5.51%)	PAS	95,015 (24.43%)	0.71	0.35	0.09	−0.02	−0.13	−0.02	0.21	0.23
	NGL	15,785 (4.06%)	0.60	0.54	0.27	0.03	−0.01	0.08	0.25	0.27
	BLF	37,420 (9.62%)	0.52	0.40	0.27	0.10	0.15	0.14	0.25	0.11
	CFF	129,960 (33.42%)	0.55	0.37	0.24	0.08	0.08	0.14	0.20	0.16
	MXF	94,833 (24.38%)	0.54	0.39	0.27	0.11	0.10	0.16	0.25	0.13
AL5	NAL	0 (0.00%)	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.	n.a.
(2.23%)	PAS	747 (0.47%)	0.34	0.69	0.43	0.13	0.13	0.16	0.29	0.40
	NGL	35,529 (22.58%)	0.34	0.66	0.46	0.17	0.10	0.16	0.36	0.42
	BLF	8214 (5.22%)	0.47	0.51	0.32	0.08	0.16	0.14	0.25	0.24
	CFF	58,431 (37.14%)	0.42	0.53	0.34	0.11	0.14	0.14	0.27	0.31
	MXF	14,883 (9.46%)	0.47	0.53	0.34	0.10	0.16	0.13	0.28	0.27

).

The second factor that significantly influences the NDVI-LST relationship is land cover (Figure 4) with a clear dichotomy of the results. Except for the marginal months of March and October, the medians of the correlation coefficients of the agricultural classes (non-irrigated arable land (NAL) and pastures (PAS)) were always significantly lower than those of the forest and grassland areas. While the agricultural areas showed predominantly negative correlation values, especially in summer (June, July, and August), none or only positive correlations occurred in the other classes.

When combining the factors land cover and altitude (Table 4), natural grasslands (NGL) at lower altitudes, in addition to the agricultural classes NAL and PAS, also showed significantly lower correlation coefficients than the forest areas except for March and October. Altitudes below 800 m which constitute more than 90% of the total area were characterized by this negative NDVI-LST correlation in the summer months, while forest areas showed no or a slightly positive correlation in summer. Above 800 m the effect of the dichotomy decreased significantly since forests predominated.

Thirdly, monthly NDVI-LST correlations vary with altitude in Bavaria (Figure 5). Except for the altitude level above 1200 m in March, the correlation coefficients increased continuously with altitude in all months. Negative correlations were on average only achieved at altitudes up to 800 m, and at altitudes below 300 m, the median of the correlation coefficients was clearly negative (from June to September). When combining altitude and land cover classes (Table 4), the results of the simple altitude classification were confirmed, whereby here even negative correlations could occur up to 1200 m in the median (July). Finally, it is important to note that hardly any strong relationships are observed within the NDVI-LST correlation analysis with respect to all described perspectives.

3.2. Drought Indices VCI, TCI, and VHI and Their Evaluation Results

Figure 6 is showing the prominent drought year 2003 as an example: The VHI is a clear combination of both indices, with weeks of low drought stress taking higher values (e.g., day 129) and weeks of high drought stress taking lower values (e.g., day 193–233). The TCI basically changes more dynamically than the VCI, both spatially and temporally. From the 20-year time series of the VCI, TCI, and VHI, all areas with a monthly mean NDVI-LST correlation coefficient >−0.1 had to be removed. While in the summer months of July and August a large part of the area was still retained (62.25 and 58.83%, respectively), in the marginal months of March (0.15%), April (8.00%), and October (5.95%) hardly any area was included in the analysis. The months of May (18.15%), June (34.31%), and September (24.04%) represented intermediate stages, in which mainly focal areas in the south-east and north-west of Bavaria were included in the analysis (Figure S6 and Table S3).

When evaluating for the above described sensitive areas the drought indices TCI, VCI, and VHI with the SMI (

Table 5. Bravais-Pearson correlation coefficients between Soil Moisture Index (SMI) and the MODIS-based drought index data sets for the complete time series from 2001 to 2020 (monthly and total). All correlation coefficients are statistically significant (p-value < 0.05). The corresponding scatterplots for the complete time series (March to October) are included in Figure S7.

SMI vs.	TCI	VCI	VHI
March	0.26	0.09	0.25
April	0.52	0.11	0.37
May	0.38	0.16	0.32
June	0.34	0.15	0.30
July	0.63	0.26	0.52
August	0.69	0.37	0.61
September	0.47	0.38	0.51
October	0.07	0.05	0.08
Total	0.54	0.27	0.48

), several observations could be made. Over the entire growing season, similar ranges of values were found for the correlation coefficients of SMI with TCI (0.54) or VHI (0.48), while the one between VCI and SMI was considerably lower (0.27). This tendency of lower VCI than TCI/VHI correspondence to SMI could be observed across all months, especially in the summer months of July and August. Generally, the MODIS-based drought indices corresponded much better to the SMI in summer than at the beginning and end of the growing season.

When comparing the MODIS-based drought indices with annual yield data for different crops, positive correlations of varying strength were found (

Table 6. Bravais-Pearson correlation coefficients (r) between crop-specific mean drought index values and the respective annual field crop or weighted yield anomalies (county level). All correlation coefficients are statistically significant (p-value < 0.05). The corresponding scatterplots between weighted yield anomaly and weighted annual mean drought index are included in Figure S8.

r	Winter Wheat	Winter Rye	Summer Barley	Oat	Sugar Beet	Winter Rapeseed	Silage Corn	Weighted Yield
TCI	0.12	0.34	0.10	0.28	0.38	0.10	0.62	0.47
VCI	0.21	0.25	0.19	0.25	0.50	0.09	0.51	0.45
VHI	0.16	0.33	0.14	0.28	0.46	0.10	0.61	0.49

). In contrast to SMI, the correlations of all the three drought indices with the yield anomalies were in a similar value range, e.g., across all field crops as indicated by the weighted yield; they ranged between 0.45 (VCI) and 0.49 (VHI). Additionally for the individual crops, no obvious differences were apparent among TCI, VCI, and VHI correspondence to yield anomalies. However, the correlation varied considerably for individual crops. While especially winter wheat, summer barley, and winter rapeseed showed low positive correlations, the correlations of sugar beet and silage corn were high.

The correlations between weighted crop yield anomalies and the corresponding TCI, VCI, and VHI at county level highlighted spatial patterns across Bavaria (Figure 7c and Figure S9e,f). Similar to Table 6, for each of the drought indices, almost all correlations were positive, with significant variations across counties. There were hardly any differences between the maps of the individual drought indices. However, all counties in the north of Bavaria tended to have higher correlation values than southern areas. Another focus area with positive correlations was the northwest of Bavaria, with statistically significant correlations >0.8. One county in the south, was also striking, showing a slight to medium negative linear correlation for all the three drought indices.

Within the correlation maps between the individual field crop yield anomalies and the drought indices TCI, VCI, and VHI, clear differences could be recognized depending on the crop type (Figure 7a,b and Figures S9–S12). While silage corn, oats, winter rye, and sugar beet showed a homogeneous picture of positive correlations, winter wheat, winter rapeseed, and summer barley showed a predominantly heterogeneous picture of positive and negative correlations. In addition, it became clear that the correlations in northern Bavaria generally had higher values than in the south. The strongly negative outlier in the southern county was also confirmed in the maps of individual crops. In contrast, there were hardly any notable differences between the individual drought indices. In conclusion, it must also be noted within this subsection that, as with the NDVI-LST correlation analysis, few strong relationships emerged in the evaluation of the drought indices.

4. Discussion

This study addresses three research questions that are not satisfactorily answered in the current literature. First, it provides results with much better spatial and temporal resolution than the current one to determine where and when water is the limiting factor for vegetation growth in Central Europe. Second, the analysis of this constraint incorporates other environmental conditions to explain it, which has rarely been the case before. Finally, an approach not yet practiced in Bavaria is presented to perform continuous drought monitoring using remote sensing-based VIs at appropriate spatial and temporal resolution.

Our study showed that the correlation between NDVI and LST depends on the season, land cover, and altitude, and that the TCI and VHI correlate well with both soil moisture and yield data in Bavaria. Considering that the NDVI-LST correlation is negative when water is the limiting factor for vegetation growth and positive when energy is the limiting factor, several conclusions can be drawn from the results presented. In Bavaria, this relationship is, firstly, dependent on the season: At the beginning and end of the vegetation period the correlation between NDVI and LST is predominantly positive, whereas in the summer months it is predominantly negative. Accordingly, the primary growth-limiting factor of energy is replaced by the factor of water, especially in July and August. This seasonal dependence of the NDVI-LST relationship has already been demonstrated globally, for example in China [67] and North America [54,55] and even both from temporal and spatial perspectives in Europe [36], although not in such a detailed spatial and temporal resolution.

Classifying the correlations according to land cover, it is noticeable that only agricultural land and grassland show negative NDVI-LST correlations in the summer months, while forest shows none or only positive correlations throughout the year. This indicates that forests have a greater water storage capacity and/or a higher resilience to drought stress than agricultural plants and grasses. Accordingly, the growth of forests in Bavaria in summer seems in this study primarily energy-limited, while water is the primary limiting factor for arable plants and grasses. However, it is important to emphasize that this is an averaging analysis over 20 years: Extreme drought years also affect forest areas in Central Europe and are associated with tree mortality [9,68].

The fact that the relationship between NDVI and LST varies depending on the land cover or vegetation type has also already been demonstrated in several regions of interest, e.g., on the Iberian Peninsula [69], in the Arctic [70], Mongolia [71], and North America [54]. The differences in forest areas and other vegetation types have only been investigated in Mongolia [71] and North America [54]. Similar correlation differences were found in Mongolia as well, although the analysis does not include seasonal differentiation. In North America, it was also shown that forest areas are less prone to negative NDVI-LST correlations than agricultural areas or grasslands. However, the different data, classification, and geographical conditions must be taken into account when making these comparisons.

The third factor that significantly influences the relationship between NDVI and LST is altitude. An increase in the correlation coefficient between NDVI and LST with altitude has similarly been observed in other parts of the world before, for example in North America [54], Mongolia [71], and Europe [36].

Additionally, in Bavaria the correlation coefficient increases with altitude in all months, with negative correlations only occurring below 800 m a.s.l. on average. Thus, especially in summer at low altitudes in Bavaria (<800 m), water is the primary limiting factor for vegetation growth, while above this altitude energy continues to be the primary factor.

Looking at the calculated drought index maps from TCI, VCI, and VHI (Figure 6), it can be seen that MODIS data can primarily capture the temporal and spatial dynamics of NDVI and LST. Comparing this with drought index calculations from other satellites (e.g., [72,73]), this is an advantage in drought monitoring, especially in the temporal dimension. However, some compromises have to be made within the spatial resolution. In general, it can thus be argued that MODIS data are particularly useful for near-real-time drought monitoring, while general classifications at selected time periods with high spatial resolution are more useful with other platforms.

When evaluating the drought indices TCI, VCI, and VHI, with the soil moisture index, two aspects are particularly relevant. Firstly, TCI and VHI show higher correlations with soil moisture than VCI, with TCI resulting in the highest correlation overall. It can be concluded from this that surface temperature generally reacts faster to changes in soil moisture than vegetation state. Further studies with similar methodology exist mainly for Asia and North America. Positive correlations between soil moisture and VHI have been observed also in [74,75,76]. In Mongolia, the correlation between TCI and soil moisture is higher than that for VCI, firstly, at different soil depths (but only with lagged measurements) [77] and secondly, in areas with predominant vegetation cover [78]. In China, soil moisture measurements and TCI also correlate higher than VCI in different months [79]. A similar picture emerges for rice fields in Vietnam, and only on forest land do both indices show a similar soil moisture correlation [80]. The results for Bavaria are thus largely supported in other regions of the world.

The second relevant aspect of the correlation analysis between soil moisture and the drought indices is the seasonal course. In summer, all indices correlate more strongly with soil moisture than at the beginning and end of the vegetation period. This suggests that soil moisture in summer has a stronger influence on both the surface temperature and the status of the vegetation. This can be explained by the fact that water is the primary limiting factor for plant growth in the warmer or hot summer months when the biomass to be sustained is high: If soil moisture shows negative or positive anomalies in summer, this has a stronger influence on the water balance of the vegetation and thus also on the temperature [81] than in the other seasons, in which energy is the primary limiting factor.

When comparing the agricultural yield data with the drought indices, three aspects are particularly central. Firstly, all three indices (TCI, VCI, and VHI) show similarly high positive correlations in general, i.e., as indicated in the weighted yields. Also, there is no general tendency for the superiority of any index for the individual crops. All three indices thus reflect the annual anomaly of agricultural yields in Bavaria well. Two aspects are noticeable here: VHI shows only slightly superior correlations here than its two components separately. On the other hand, VCI and TCI correlate with yield at the same level, although VCI is assumed to have a much more direct link to agricultural yield via NDVI than the TCI via LST. This suggests that drought, which is reflected in higher than normal vegetation thermal conditions, was the main cause of yield losses in the study area. Other correlation studies in Europe involving yield data have also shown positive correlations with the three drought indices almost without exception [61,64,82,83]. However, a direct comparison is usually difficult due to the different methodologies, data situations, selected crops, and geographical locations. When comparing with the most methodologically similar study by Bachmair et al. [28], the results are confirmed in terms of both correlation strength and spatial distribution, although in the other study the VCI shows slightly higher correlations than the VHI. The added value in the present paper compared to this study is mainly the analysis of several crops and the inclusion of the TCI.

The second aspect worth mentioning is the difference in correlation strength between the individual crops. Especially sugar beet and silage corn show higher correlations than the others, which can be explained primarily in the monitoring period: While most crops are harvested at the end of July/beginning of August, sugar beet and silage corn are harvested at the end of September/beginning of October (Table 3). Thus, for these two crops, more growing time and area (especially in August) is included in the correlation analysis, which is obviously reflected in the drought indices. In comparison with other studies, it is noticeable that winter wheat yields in north-eastern Germany [32] show higher correlations with the VIs used than in Bavaria. The decisive difference here, however, is in the methodological choice to include only relevant areas in north-eastern Germany in the monitoring or the correlation analysis, whereas in our study all vegetation areas with negative NDVI-LST correlation were included. Our results of the yields for wheat and barley in connection with the drought indices are also confirmed by a Europe-wide temporal correlation analysis between NDVI and wheat and barley yields [84]. Here, both correlations in Germany decrease significantly from June onwards (Day 153) and reach similar value ranges as in our analysis.

The third and last aspect to be noted is the spatial difference within the correlation analysis. The general impression is that for all indices and crops the correlation between yield and drought index is higher in northern/eastern Bavaria than in southern Bavaria. The reason for this can be traced in the different precipitation conditions (Figure S3): In drier northern Bavaria, the sensitivity of both the drought indices and agricultural yields to a change in the water balance is higher than that in wetter southern Bavaria. The higher sensitivity causes a higher correlation and thus a stronger relationship between the two variables.

Nevertheless, certain limitations of this study should be noted. The spatial resolution of the underlying data should be thoroughly considered. Due to the small-scale features of the landscape in Bavaria, mixed pixels often occur, especially in the remote sensing data, which can lead to inaccuracies both in the correlation analysis between NDVI and LST and in the index calculation. Within the evaluation data, this indication is equally true for both the soil moisture data (4 km) and the agricultural yield data (county level). In addition, it must be noted that the data in this paper show significant differences within temporal and spatial resolution. In order to match the data, they are adjusted on both levels (e.g., via resampling or interpolation). This causes additional inaccuracies and could be an explanation for partially low correlations within the results. In addition, the assumption of a 19-year identical distribution of crop shares within the counties introduces—to some extent—further uncertainties.

Another aspect worth discussing is the general use of VIs for estimating vegetation growth and the application of NDVI in particular. Thus, VIs only ever provide an indirect vegetation proxy, and several limitations must be considered here: On the one hand, detection on heterogeneous surfaces and at different canopy heights is problematic, and on the other hand, factors such as sensor calibration, sensor viewing conditions, solar illumination geometry and soil moisture, also influence the data quality of the indices [85,86,87]. NDVI in particular is one of the most widely used vegetation indices and thus offers a wide range of comparisons. It allows a large number of estimations of different vegetation properties and, in addition, it has a sensitivity to green vegetation even on areas with low vegetation cover. Nevertheless, the NDVI also brings uncertainties in the detection due to the influence of soil brightness, soil color, clouds and cloud shadows as well as leaf canopy shadows and a saturation problem in the presence of high vegetation diversity [86,88,89,90].

From a methodological point of view, two aspects, in particular, should always be kept in mind: Firstly, the application of the drought indices is in part severely limited in terms of area due to the threshold value set for the correlation coefficient between NDVI and LST (which refers to the entire period and partially masks individual drier years). The evaluation of the indices is thus only of limited significance, especially in the marginal months of the vegetation period and in counties with little area included. On the other hand, the area included in the analysis is always the same when evaluating the drought indices with the individual crop yields and no differentiation is made between the individual crops.

5. Conclusions

In this study, we aimed at identifying by means of a correlation analysis between NDVI and LST when and where in Bavaria vegetation growth is primarily water-limited and how this is related to the factors of land cover and altitude. In addition, we investigated whether and to what extent the remote sensing-based drought indices TCI, VCI, and VHI can capture agricultural drought and yield losses within these areas. The indices were calculated from 2001 to 2020 within the growing season and evaluated with both soil moisture and yield anomalies of agricultural crops.

We found that in Bavaria, especially in the summer months of July and August, on agricultural land and grassland and below 800 m, water is the primary limiting factor for plant growth. Within these areas and periods, the remote sensing-based drought indices TCI and VHI correlate strongly with soil moisture and agricultural yield anomalies. From both a soil- and a vegetation-based perspective, both indices have the potential to detect agricultural and vegetation-based drought, respectively.

However, there are also further research potentials within this thematic focus. With regard to the LST-NDVI relationship, it would be particularly interesting to include in the correlation analysis other relevant variables, such as evapotranspiration, precipitation, and global radiation and to specifically analyze individual drought years with regard to this relationship. Within the evaluation of the drought indices, further soil- and vegetation-related variables could be included in the analysis. In addition, a temporally differentiated (e.g., monthly) correlation analysis between the drought indices and the field crop yield anomalies would be desirable to determine possible seasonal focal points in the context.