Drought Extent and Severity on Arable Lands in Romania Derived from Normalized Difference Drought Index (2001–2020)

Abstract

The aim of this study was to evaluate the frequency and severity of drought over the arable lands of Romania using the Normalized Difference Drought Index (NDDI). This index was obtained from the Moderate Resolution Imaging Spectro-Radiometer (MODIS) sensor of the Terra satellite. The interval between March and September was investigated to study the drought occurrence from the early stage of crop growth to its harvest time. The study covered a long period (2001–2020), hence it is able to provide a sound climatological image of crop vegetation conditions. Corine Land Cover 2018 (CLC) was used to extract the arable land surfaces. According to this index, the driest year was 2003 with 25.6% of arable land affected by drought. On the contrary, the wettest year was 2016, with only 10.8% of arable land affected by drought. Regarding the multiannual average of the period 2001–2020, it can be seen that drought is not a phenomenon that occurs consistently each year, therefore only 11.7% of arable land was affected constantly by severe and extreme drought. The correlation between NDDI and precipitation amount was also investigated. Although the correlations at weekly or monthly levels are more complicated, the annual regional mean NDDI is overall negatively correlated with annual rainfall. Thus, from a climatic perspective, we consider that NDDI is a reliable and valuable tool for the assessment of droughts over the arable lands in Romania.

1. Introduction

Drought is generally considered the most complex meteorological phenomenon [1], given that many factors contribute to its onset, such as precipitation amount, soil characteristics, terrestrial water accessible to plants, soil/air temperature and humidity or wind speed. Other factors that define the characteristics of the active surface, or the physiological peculiarities of the plants, in addition to the anthropogenic influence on the environment, are highly significant in its occurrence [1,2]. Despite the policies and efforts to reduce extreme weather effects, droughts will remain unavoidable.

As present, because regional development is considered to be one of the main factors contributing to economic growth, the better prepared a region to cope with adverse weather conditions, such as drought, the more the region can contribute to the development of the whole country [3,4]. Therefore, strategies and actions to mitigate drought impacts in less developed countries are essential, because these regions are among the most vulnerable, with the lowest financial and technical capacity to adapt or to mitigate the effect of this extreme phenomenon [5].

Scientific research has led to a wide variety of results and applications in the field of monitoring and control of drought effects [2]. However, several scientific problems and challenges remain, such as finding solutions to mitigate drought effects and improve living standards, which can be solved through a better understanding of meteorological extremes, with the aim of preventing possible future impacts [1,2].

The relevance of this topic derives, in particular, from the fact that in the context of climate change, the impact of this climatic phenomenon is expected to become more pronounced because periods of drought are expected to become longer in numerous regions of the world [5,6].

Classically, Wilhite and Glantz [7] defined four types of drought: meteorological (lack of precipitation), hydrological (lack of water supply), agricultural (crop water deficit) and socio-economic (combined effect of drought on human activities). Regardless of its type, the drought phenomenon has been investigated using a variety of indices which have been developed over time, such as the Palmer Drought Sensitivity Index (PDSI) [8], the Standardized Precipitation Index (SPI) [9,10], or the Standardized Precipitation Evaporation Index (SPEI) [11]. In addition, some indices based on the relationship between land surface temperature (LST) and land cover have been developed, such as the Temperature Condition Index (TCI) [12], Vegetation Temperature Condition Index (VTCI) [13] or Temperature Vegetation Dryness Index (TVDI) [14,15].

Taken into account the complicated definition of drought [7,16], vegetation indices appear to provide more appropriate estimations of drought occurrence and intensity and the associated agriculture impacts. In this regard, remote sensing products such as the Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), Leaf Area Index (LAI) and Fraction of Photosynthetically Active Radiation (FPAR) have also been developed [17,18,19,20,21,22]. Furthermore, other indices have been investigated, such as the Vegetation Condition Index (VCI) [12], Crop Water Stress Index (CWSI) [23], Vegetation Health Index (VHI) [12], Global Vegetation Moisture Index (GVMI) [24], Soil Water Index (SWI) [25] or Remote Sensing Drought Risk Index (RSDRI) [26,27]. Among these, one of the most flexible and useful indices to be applied for drought investigation on arable lands is the Normalized Difference Drought Index (NDDI) [20].

In Europe, many studies have focused on drought and its characteristics using various methods, indices and satellite products. For example, a study focusing on the north-east of the Iberian Peninsula developed by Vicente-Serrano [28] used vegetation indices derived from AVHRR (Advanced Very High Resolution Radiometer) images. In Spain and the Mediterranean region, Vicente-Serrano et al. [29] used AVHRR images and the NDVI index to study the drought impact over agricultural lands located in the Middle Ebro valley, one of the most arid regions in Europe. In addition, Gouveia et al. [30] analyzed the drought impacts on vegetation over the entire Mediterranean basin, using NDVI and SPEI indices, with the purpose of determining the stage at which vegetation is more impacted by drought. Drought was analyzed by Sepulcre-Canto et al. [31], who combined the SPI, the anomalies of soil moisture and the anomalies of the FPAR. Furthermore, Dalezios et al. [32] applied a number of drought indices based on NOAA (National Oceanic and Atmospheric Administration)-AVHRR and the Reconnaissance Drought Index (RDI) in Thessaly, central Greece, which is a drought-prone agricultural region characterized by vulnerable agriculture. The driest years, such as 2000, 2003 or 2008, were analyzed in many studies aiming to monitor the drought severity in Europe. For this, Sea-Viewing Wide Field-of-view Sensor (SeaWiFS) and Medium Resolution Imaging Spectrometer (MERIS) instruments [33] or MODIS vegetation indices (VIs), NDVI and enhanced vegetation index (EVI) [34] have been used. It should be noted that the assessment of drought at the continental scale identified the Carpathians region as subject to an increase in drought for 1950–2012 [35].

With the increase in the availability of remote sensing products and indices, their use has also become more frequent in Romania. These have been used to analyze and detect floods [36], land cover changes [37,38], landslide-prone hilly areas in Moldova and various other areas [39] and Saharan dust intrusion [40,41], and it is clear that they can also help to precisely monitor drought.

A comprehensive analysis of the correlation between NDVI and SPEI, aiming to evaluate the response of vegetation’s photosynthetic activity to drought conditions from 1998 to 2014 over Romania and the Republic of Moldova, was undertaken by Páscoa et al. [11]. Changes in the forest ecosystems in south-western Romania, due to global climate change and anthropogenic impacts during the past three decades, and correlated with the evolution of aridization, were assessed by Prăvălie et al. [42] using the NDVI and the UNEP (United Nations Environmental Programme) aridity index. In addition, remote sensing data have been used in several studies to monitor the summer surface urban heat island of the city of Bucharest [43], Cluj-Napoca [44], and Galați [45]. Due to climate change, some studies have discussed the topic of desertification [46], which was studied in the south-east of Romania using LANDSAT TM images and LST.

Recently, Angearu et al. [47] analyzed drought severity using the Drought Severity Index (DSI) in Romania and its validation based on meteorological data, soil moisture content and agricultural production. In addition, drought assessment based on a multi-temporal analysis and trends of the DSI obtained from Terra MODIS satellite images was undertaken.

In this context, the main objective of our study was to enhance the knowledge on drought frequency and severity in Romania by providing a long-term (2001–2020) and comprehensive view of its impact over arable lands as reflected by the NDDI. Our results show that NDDI is a highly reliable and valuable tool for the assessment of droughts over the arable lands in Romania.

2. Materials and Methods

2.1. Study Area and Its Geographical Features

With a territory of about 238,500 km², Romania is the largest country in south-eastern Europe (Figure 1). The distribution of major landforms, with 31% mountainous area, 33% hills and sub-mountainous areas and 36% plains and meadows, provides its territory with significant climatic diversity [48].

First, regarding latitude, the average annual temperature decreases by 3 °C between the south and the north of Romania (from 11 to 8 °C), and regarding altitude it decreases by about 14 °C from the lowlands to the highest mountain peaks (from 11 to −3 °C). The amounts of precipitation are also strongly influenced by the topography. Compared to the average values in the plain areas from the west (about 600 mm) and east (about 400 mm) of the country, in the high mountain areas average precipitation rises on the slopes exposed to the advection of humid air masses at more than 1400 mm [49].

Second, the western and central regions of Romania are particularly impacted by the cyclones formed in the Atlantic and Mediterranean Seas (following a Pannonian track). These cyclones produce significantly more precipitation compared to the cyclones of Mediterranean origin, following a trans-Balkanian track, which produce more significant precipitation in the southern and eastern regions of the country [41,49,50,51,52,53].

It should be noted that the multiannual amount of precipitation at the country scale remained generally stable during the last interval [50], except for some areas in the north and north-west of Romania, with a positive trend, and from the east, south and south-east of the country, with a negative trend [54]. At the same time, there was an increase in evapotranspiration in Romania which led to increased aridity [46,50], which led some authors [46] to discuss a possible ongoing process of desertification, similarly to that in other regions of southern Europe, such as Spain [55], Italy [56] and Greece [57].

The arable lands in Romania cover an area of about 11 million hectares, of which more than 50% are located in the plain regions of the south, west and extreme north-east of the country, where the arable lands represent the main type of land cover (Figure 1,

Table 1. Territorial extension of arable lands according to CLC 2018 in the main regions of Romania (thousands of hectares).

Territorial Extension of Arable Lands in Romania (Source: CLC 2018)
Carpathians	Transylvanian Plateau	Subcarpathians	Pannonian Plain	Crisana and Banat Hills	Moldavian Plateau	Romanian Plain	Getic Plateau	Dobrogea Plateau and Danube Delta
270	750	313	2600	363	1390	3.070	680	1680

).

Moreover, due to the large extent of arable lands, it should be noted that Romania represents one of the main producers of maize and wheat in Europe [50], but given the small ratio of irrigated land, representing less than 10% of all arable lands [58], the annual rank in cereal production at European level is highly variable depending on year-to-year weather conditions.

2.2. Data Used

2.2.1. Moderate Resolution Imaging Spectroradiometer (MODIS) Data

The MODIS Vegetation Indices (MOD13Q1) Version 6, based on a 16-day composite at 250 m spatial resolution [59,60,61], were used to calculate the NDDI. The MOD13Q1 product includes 12 layers, of which only three layers were used: band 1 (16 days Red reflectance—620 to 670 nm), band 2 (16 days Near Infrared reflectance—841 to 876 nm) and band 7 (16 days Middle Infrared reflectance—2105 to 2155 nm) (

Table 2. Composite periods of MOD13 used in the study [61].

Composite	Day of Year (No.)	Starting Day during Non-Leap Years	Starting Day during Leap Years *	Season
1	65	6-Mar	05-Mar	Spring
2	81	22-Mar	21-Mar
3	97	7-Apr	06-Apr
4	113	23-Apr	22-Apr
5	129	9-May	08-May
6	145	25-May	24-May
7	161	10-Jun	09-Jun
8	177	26-Jun	25-Jun	Summer
9	193	12-Jul	11-Jul
10	209	28-Jul	27-Jul
11	225	13-Aug	12-Aug
12	241	29-Aug	28-Aug

* Due to numerous corrections made to the acquired data, both in leap years and in non-leap years, 16-day Moderate Resolution Imaging Spectroradiometer (MODIS) composites have the same data range [61].

and

Table 3. Layers of MOD13Q1 product [59]; the bands used to calculate the Normalized Difference Drought Index (NDDI) (1, 2 and 7) are given in bold.

Layer Name	Description	Units	Data Type	Fill Value	No Data Value	Valid Range	Scale Factor
250 m 16 days NDVI	16 day NDVI	NDVI	16-bit signed integer	−3000	N/A	−2000 to 10,000	0.0001
250 m 16 days EVI	16 day EVI	EVI	16-bit signed integer	−3000	N/A	−2000 to 10,000	0.0001
250 m 16 days VI Quality	VI quality indicators	Bit Field	16-bit unsigned integer	65535	N/A	0 to 65534	N/A
250 m 16 days red reflectance	Surface Reflectance Band 1	N/A	16-bit signed integer	−1000	N/A	0 to 10,000	0.0001
250 m 16 days Near Infrared reflectance	Surface Reflectance Band 2	N/A	16-bit signed integer	−1000	N/A	0 to 10,000	0.0001
250 m 16 days blue reflectance	Surface Reflectance Band 3	N/A	16-bit signed integer	−1000	N/A	0 to 10,000	0.0001
250 m 16 days Middle Infrared reflectance	Surface Reflectance Band 7	N/A	16-bit signed integer	−1000	N/A	0 to 10,000	0.0001
250 m 16 days view zenith angle	View zenith angle of VI Pixel	Degree	16-bit signed integer	−10000	N/A	0 to 18,000	0.01
250 m 16 days sun zenith angle	Sun zenith angle of VI pixel	Degree	16-bit signed integer	−10000	N/A	0 to 18,000	0.01
250 m 16 days relative azimuth angle	Relative azimuth angle of VI pixel	Degree	16-bit signed integer	−4000	N/A	−18,000 to 18,000	0.01
250 m 16 days composite day of the year	Day of year VI pixel	Julian day	16-bit signed integer	−1	N/A	1 to 366	N/A
250 m 16 days pixel reliability	Quality reliability of VI pixel	Rank	8-bit signed integer	−1	N/A	0 to 3	N/A

). These MODIS products were extracted from the Land Processes Distributed Active Archive Center (LP DAAC), using Application for Extracting and Exploring Analysis Ready Samples (AppEEARS), ver. 2.42.1 [62]. Thus, 721 Surface Reflectance Bands for 16-day composite images (1, 2 and 7 bands) were used to calculate the NDDI.

The 16-day composite was chosen because it minimizes the errors induced by clouds. However, these errors generated by the clouds were still present and were eliminated as further explained in Section 2.3.1.

2.2.2. Corine Land Cover (CLC) Data for Arable Lands

The arable lands from Romania were extracted from CLC 2018. CLC 2018 is a dataset produced within the framework of the Copernicus Land Monitoring Service that refers to the land cover status of 2018 in Europe. The CLC service has a long heritage (formerly known as the “CORINE Land Cover Program”), and is coordinated by the European Environment Agency (EEA). It provides consistent and thematically detailed information on land cover and land cover changes in Europe. This project has a regular update time period of approximately six years [63]. In this classification, arable lands represent a distinct land cover category.

2.2.3. Precipitation Data

The relationship between atmospheric precipitation amount and NDDI was also analyzed, in order to understand how the crop vegetation responds to this major atmospheric driver of drought occurrence. The precipitation time series used in this study were taken from the ENSEMBLE project gridded data set E-OBS [64]. These data were obtained using a kriging interpolation procedure from the European Climate Assessment and Dataset (ECA&D) time series at meteorological stations [65]. The E-OBS version used in this study was version 21 (release date: May 2019) and covers Europe with a spatial resolution of 0.1° from 1950 to 2019 at a daily time step. The precipitation amount for the year 2020 was added from the ERA-5 land reanalysis hourly dataset of the Copernicus Climate Change Service C3S Climate Data Store (CDS) [66]. Reanalysis was conducted by combining the model data with observations across the world into a complete and consistent dataset. The spatial resolution of the ERA-5 land reanalysis is 0.1°, thus no resampling was needed. The precipitation amount was calculated for each 16-day composite (12 composites between March and September) and for each year (20 years), similar to that of the NDDI.

2.3. Methodology

2.3.1. Gap-Filling of MODIS Images

The quality of remote sensing products is highly influenced by weather conditions. Among these factors, cloud cover can frequently induce gaps in the time series of the satellite optical imagery. To fill these gaps retrieved from MODIS over Romania, the Data Interpolating Empirical Orthogonal Functions (DINEOF) procedure was applied [67,68,69]. Recently, this method has been used in various remote sensing products such as for the reconstruction of total suspended matter [70] sea surface salinity [71], sea surface temperature derived from MODIS [68,72], MODIS-Aqua chlorophyll products [73] or LST over Bucharest [74]. More appropriate to our study, Filliponi et al. [75] applied the DINEOF algorithm to the reconstruction of the MODIS Fraction of Green Vegetation around the world.

The DINEOF procedure was run to gap-fill all MODIS NDDI composites and reconstruct the missing pixels. A full completeness (100% availability) composite (August 2012) was set as a profile mask, thus only pixels within the arable lands were filled. The DINEOF gap-filling method was applied using the rtsa R package version 0.3 for Raster time series Analysis (https://github.com/ffilipponi/rtsa/blob/master/DESCRIPTION (accessed on 15 February 2020) [75]. To evaluate the DINEOF method, artificial gaps were created for a full completeness composite and compared with the original. The difference between pixels of the original composite and DINEOF gap-filled pixels are shown in

Table 4. The statistical parameters of the original and gap-filled data pixels (//–not applicable).

Summary Statistical	Original Data	Gap-Filled Data	Original vs. Gap-Filled Data
Minimum	0.02	−0.08	//
1st Quartile	0.34	0.31	//
Median	0.46	0.44	//
Average	0.48	0.46	//
3rd Quartile	0.59	0.57	//
Maximum	1.75	1.72	//
Root-Mean-Square Error	//	//	0.03
Mean-Absolute-Error	//	//	0.09
R-squared	//	//	0.91

. For more details regarding this procedure see Figures S1–S3 in the Supplementary Materials.

2.3.2. NDDI Calculation and Drought Assessment

Generally, no index can fully describe the complexity of drought at both temporal and spatial levels. Hence, it is recommended to combine several parameters, indicators or indices (including remote sensing data) in a single product for drought classification [20].

First, band 1 (B1), band 2 (B2) and band 7 (B7) from 16-day composite images at 250 m resolution of the MOD13Q1 product were extracted. Using these bands, the NDVI was calculated using B1 and B2 bands [20], by applying Equation (1), and NDWI was calculated from B2 and B7 [20,76], by applying Equation (2).

The NDVI is a classical index which measures the development and the density of vegetation and has values from −1.0 to 1.0. Negative values indicate clouds or water, whereas positive values indicate soil without vegetation (values near to zero), and dense green vegetation (values equal or higher than 0.6) [77,78]. The NDVI is widely used to evaluate the main parameters of vegetation, induced mainly by climate conditions, human activities and other anthropic or natural causes. NOAA/AVHRR, SPOT (French: Satellite Pour l’Observation de la Terre), MODIS or LANDSAT [9,10] imagery can be used to achieve these products. It should be noted that the NDVI is included in this product and can be downloaded already calculated. The NDVI was calculated using the following formula [20]:

$${\mathrm{NDVI}}_{\mathrm{Modis}}=\frac{B2-B1}{B2+B1},$$

(1)

where B1 and B2 refer to MODIS band 1 and band 2, respectively.

Additionally, the NDWI is an index which measures the water content of leaves and is used for detecting and monitoring vegetation humidity. The NDWI is influenced by plant dehydration, and it is considered to be a better indicator for drought monitoring than the NDVI [77]. The NDWI also has values from −1.0 to 1.0. The common range for green vegetation is −0.1 to 0.4. This index increases with vegetation water content or from dry soil to free water [78,79]. Both the NDVI and NDWI have been used in different studies to observe their relationship with LAI for the study of the characteristics of vegetation that covers different regions, including arable lands [80,81,82] or FPAR [19,83]. The NDWI was calculated using the following formula [20]:

$${\mathrm{NDWI}}_{\mathrm{Modis}}=\frac{B2-B7}{B2+B7},$$

(2)

where B2 and B7 refer to MODIS band 2 and band 7, respectively.

It has been found that the NDWI is more sensitive than the NDVI to drought conditions, providing information about the amount of water that enters the plant [20]. In addition, it has been found that the average in cases of drought is below 0.5/0.3 in the case of the NDVI/NDWI, whereas in periods without drought, the NDWI/NDVI has values above 0.4/0.6 [6].

Although meteorological drought has been and can be well studied further using the NDVI and NDWI, a new index combines the information provided by the NDVI and NDWI. The NDDI was recently developed [20], and has been used to monitor the drought parameters in different regions of the world [6,84,85,86,87], representing a sensitive drought assessment tool for agriculture [88]. The first research to analyze the potential of this drought monitoring index was conceived with good results for the Flint Hills region of eastern Kansas and north-eastern Oklahoma [89].

Using the NDVI and NDWI results, the NDDI was calculated according to the equation below [20]:

$$\mathrm{NDDI}=\frac{\mathrm{NDVI}-\mathrm{NDWI}}{\mathrm{NDVI}+\mathrm{NDWI}},$$

(3)

The resulting values of the NDDI range generally from 0 (no drought) to >1.0 (extreme drought).

The NDDI has a stronger response to summer drought conditions than a simple difference between the NDVI and NDWI, and is more sensitive indicator of drought in grasslands and arable lands than the NDVI alone. Because NDWI values decrease more than NDVI values during summers with severe drought, suggesting that the NDWI is more sensitive than the NDVI to drought conditions, the calculation of the NDDI is considered to be a more complex calculation compared to the NDWI and NDVI [20].

Finally, it was found that the NDDI combines well the information provided by the NDVI and NDWI, that it has a wider range of values than a simple NDVI–NDWI differentiation, and that it can be used, based on MODIS images at a good resolution, for the analysis of drought at local scales [90]. Thus, it has been used increasingly often in different regions to study the extent and severity of drought, especially during the vegetation period [6,22,86,91]. Because the NDDI is more sensitive and more accurate, drought-affected territories will be identified more often compared to using the NDVI or NDWI, with differences of up to 5% [6,90]. In Romania, this index has also been used based on MODIS images but to a small extent [83,84,92].

Using the final NDDI products for the 12 composites for each year between 2001 and 2020, a total of 240 16-day composites were derived. In these composite images, the absolute and relative frequency of NDDI values indicating drought (>0.5) were computed for each pixel. To assess the drought severity, the drought frequency for each pixel was considered, taking into account the NDDI classes [84,89] higher than 0.5, which indicate moderate drought (0.5–0.6), severe drought (0.6–1.0) and extreme drought (>1.0). Furthermore, these results were aggregated at annual and multiannual levels.

All of the cartographic presentations of the drought spatial distribution in this paper were constructed using ArcGIS software, version 10.3, produced by ESRI (Environmental Systems Research Institute).

2.3.3. Spearman’s Correlation Analysis between NDDI and Precipitation Amount

To check the relationship between the NDDI and the precipitation amount—one of the major drivers of drought conditions—Spearman’s correlation was applied. Using an iteration for each pixel, Spearman’s correlation was first applied between each 16-day composite precipitation amount and the corresponding NDDI values for the same period, and then for the entire analyzed period (2001–2020). In this analysis a significance level of p-value <0.10 was used.

However, the relationship between atmospheric precipitation and the NDDI is far more complex. Therefore, in addition to a direct Spearman correlation between the two parameters, we determined the inertial effect of atmospheric precipitation on the state of crop vegetation. For this, a lagged correlation was applied between the composite n of atmospheric precipitation and the composite n + 1 of the NDDI (1lagged composite correlation). Similarly, we applied the correlation for 2composites lagged (n composite of atmospheric precipitation with n + 2 composite of NDDI) to fully cover the inertial response of the NDDI to precipitation input.

3. Results and Discussions

3.1. Drought Extent and Severity According to NDDI

The results of our analysis indicate firstly that, as a multiannual mean, 17.2% of arable land was affected by drought during the analyzed period (2001–2020), with a larger extent during the very dry years (25.6% in 2003, 24.1% in 2012, 23.0% in 2002 and 21.9% in 2020). By comparison, during very humid years, the arable lands were clearly less affected by drought (10.8% in 2016, 11.0% in 2014, 11.1% in 2018 and 11.2% in 2010). Thus, as a general feature, drought has constantly affected the territory of Romania, although only some regions have been severely impacted (

Table 5. Relative frequency (%) of NDDI drought classes between 2001 and 2020 (ToD—Type of drought, Nd—No drought, Md—Moderate drought, Sd—Severe drought, Ed—Extreme drought and Ad—All types of drought) and the precipitation amount (mm) accumulated between 5 of March and 13 of September.

Normal Difference Drought Index (NDDI)
ToD	Nd	Md	Sd	Ed	Ad	Precipitation Amount (mm)
Range	<0.5	0.5–0.6	0.6–1	>1	>0.5
2001	80.2	6.4	10.5	2.9	19.8	370.9
2002	77.0	7.2	12.3	3.5	23.0	310.0
2003	74.4	7.4	13.6	4.7	25.6	207.2
2004	81.3	5.5	9.8	3.5	18.7	319.8
2005	81.5	5.4	10.1	3.0	18.5	458.7
2006	83.0	4.9	9.8	2.3	17.0	375.3
2007	78.8	7.9	10.9	2.3	21.2	295.1
2008	87.0	4.5	6.9	1.5	13.0	289.9
2009	85.5	5.6	7.7	1.2	14.5	276.1
2010	88.8	4.0	6.1	1.1	11.2	381.0
2011	84.6	5.1	8.4	1.9	15.4	270.3
2012	75.9	7.2	13.4	3.5	24.1	267.3
2013	85.4	4.9	8.2	1.5	14.6	339.3
2014	89.0	3.8	6.0	1.3	11.0	405.5
2015	82.4	5.8	9.6	2.2	17.6	262.3
2016	89.2	3.7	5.9	1.2	10.8	329.4
2017	84.2	5.2	8.5	2.1	15.8	316.6
2018	88.9	4.2	5.9	1.0	11.1	340.9
2019	80.5	5.5	10.6	3.4	19.5	286.3
2020	78.1	7.2	11.8	2.9	21.9	360.1
2001–2020	82.8	5.6	9.3	2.4	17.2	323.1

and Figure 2). Taking into account the fact that more than 70% of the arable lands in Romania are found in plain regions, we observe that these are the regions constantly facing the risk of drought occurrence and the associated impacts on agriculture production.

Additionally, the trend analysis applied to all drought types (moderate, severe, extreme) between 2001 and 2020 (Table 4) indicates a decreasing trend (R² = 0.14, p-value < 0.10), in contrast to other trends identified using other methods for drought assessment [46]. More details regarding the trend analysis are given in Figures S4 and S5 in the Supplementary Materials. Moreover, it can be observed that the first part of the analyzed period (2001–2010) recorded more years in terms of drought extent and severity than the second interval (2011–2020). Therefore, we can assume that the discussed increase in the frequency of drought events derived mainly from standardized precipitation index [93,94,95] is not firmly supported from the perspective of crop vegetation conditions, assessed using the NDDI.

In terms of spatial distribution of drought for the 12 annual composites (Figure 3), and from year to year (Figure 4), we can see that the NDDI drought classes are present in all of the arable regions of Romania. As a main feature, a general country scale asymmetry, with the south-eastern regions more prone to drought than the western side of the country, can be clearly seen. This asymmetry is explained partially by the asymmetry in the field of atmospheric precipitation between these regions, but probably also by a plethora of other factors. For instance, it is known that the NDDI can also be influenced by the prevailing crop type (more wheat in the south-east of Romania), by some types of soils that easily lose the water reserve in the upper layers (such as those in the south and south-east of Romania developed on loess and sands), and by the level of underground water. For instance, analyzing the results for each composite between March and September (Figure 3), it can be seen that drought is less pronounced from April to June, an interval corresponding both to the peak in annual precipitation amount in Romania and to the maximum crop vegetation development. After the end of June, the drought increases in spatial frequency as a combined effect of the decrease in precipitation amount and the depletion of the arable lands by vegetation, due to the harvest of some important crops, such as wheat.

In general, the south-east of Romania was affected by drought even during the wettest years, such as in 2010, 2014, 2016 or 2018 (Figure 4). Normally, if we discuss the manifestation of the desertification process in Romania [46], one should expect this phenomenon in those regions that are affected by drought even during the wettest years. However, our analysis does not support the so-called theory of desertification in Romania, at least from the perspective of crop vegetation conditions.

To underline the variability of drought occurrence at the country scale, we identified, for each pixel on the map of arable lands in Romania, the year recording the highest value of the NDDI, indicating either severe or extreme drought conditions. The results were simplified by grouping them into four classes that indicate the 5 year interval in which these values were recorded (Figure 5).

Although the southern and western parts of Romania recorded the highest values of the NDDI, mainly between 2001 and 2005 (due to 2003 and 2002 drought events), and secondly during 2006–2010 (as an effect of the 2007 drought episode), in the central and north-western parts of Romania, over large parts of the arable lands, the interval 2011–2015 had the highest drought impact (particularly due to the 2012 drought event). For the north-eastern part of Romania, even if 2001–2005 recorded most of the maximum values of the NDDI, no interval appears to be dominant for the most severe drought conditions as a whole: 2016–2020 prevailed over the southern part of the region together with the Bărăgan region, 2006–2010 prevailed over the central part of the region, and 2011–2020 was more present over the northern part, in particular (Figure 5).

Consequently, the most important region subject to drought in Romania (Figure 6), as derived from the NDDI values, extends over the southern part of the Moldavian Plateau and the eastern part of the Romanian Plain (the so-called Bărăgan Plain, one of the most important agricultural regions in Romania). However, within this region, drought was not pronounced along the valley of the rivers (Siret, Buzău, Ialomița). In addition, we can distinguish a compact strip of arable land oriented from north to south, located east of the Romanian Plain and west of the Dobrogea Plateau, along the Danube, between the cities of Galați and Călărași, which is not severely affected by drought. This region is represented by lands between the Danube branches, recording a high degree of soil humidity and benefiting from well-developed irrigation systems [47]. A second important region impacted by drought is located in the north-eastern part of Romania, within the Moldavian Plain. In addition to these two important regions, drought is also common in the Dobrogea region, in the southern part of the Romanian plain and in the western part of the country, in a region that is relatively distant from the extremity of the Pannonian Plain.

It is important to underline that all of these regions that are affected by drought have at least three elements in common: the low amount of precipitation with a multiannual mean less than 300 mm between March and September [49], a very high level of groundwater vulnerability [96,97] and the prevalence of wheat crops in the arable lands.

For a more comprehensive view on drought, in addition to the frequency of drought classes for the entire period from 2001 to 2020, we selected the upper/lower third (that is, 7 years) of the driest/wettest years from 2001 to 2020, according to the results shown in Table 2. We thus present more synthetically the distribution of the frequency of drought classes of the NDDI, selecting the most relevant years for drought/humid conditions (Figure 7 and Figure 8).

The driest 7 years between 2001 and 2020, as shown by the NDDI values (Table 5), were by 2003, 2012, 2002, 2020, 2007, 2001 and 2019. Clearly 2003, with 7.4% of the arable land affected by moderate drought and 19.0% affected by severe and extreme drought, represents a record year from this point of view. The second driest year was 2012 with 7.2% of arable land affected by moderate drought, and 16.9% affected by severe and extreme drought.

The spatial extent of drought as shown by the NDDI values for these years presents a similar pattern as for multiannual mean, but extreme drought conditions are extended over the Dobrogea Plateau, the north-eastern part of the Moldavian Plateau and the southern central part of the Romanian Plain (Figure 7). For the western part of Romania during these years, the differences from the multiannual mean are not high, indicating a less pronounced impact of drought in this region.

In contrast, the wettest 7 years between 2001 and 2020, as shown by the NDDI values, were 2016, 2014, 2018, 2010, 2008, 2009 and 2013. The wettest years were 2016 and 2014 with only 10.8% and 11%, respectively, of arable land affected by drought (Figure 8, Table 5). Regarding the territorial extension of the drought during these intervals, the most affected areas remained the eastern and southern parts of Romania. Certainly, during these years the drought manifested with a lower degree of severity. However, even in these general humid conditions the drought remained a threat in some regions pf eastern Romania, particularly in the northeast and southeast of the Moldavian Plateau, in the eastern part of the Romanian Plain and the Dobrogea Plateau. Because these are the regions in Romania that are constantly affected by drought, drought in these regions should not be considered an extreme meteorological event, but a common climate feature. That is, in these regions agriculture is not possible without intensive irrigation.

3.2. The Relationship between Atmospheric Precipitation and NDDI

The relationship between the atmospheric precipitation amount and the NDDI is governed by a logical inverse correlation (−0.37 Pearson coefficient, statistically significant for p < 0.10). In addition, the difference (70.4 mm) between the means of the precipitation amount for the five most dry/humid years, defined using the NDDI (Table 5), is statistically significant at p < 0.10.

Generally, we can observe that the 1-composite lagged correlation is prevalent for the entire interval with a spatial frequency of pixels with significant correlation reaching the maximum (ca. 30% of all the arable lands) from May to July (Figure 9). Moreover, for March and April (with a minimum between 23 April and 8 May) the direct correlation is very weak, as also observed by Potopová et al. [98] in the Republic of Moldova for the correlation between the NDVI and SPEI. The thermal increase in this period most likely causes a rapid increase in vegetal activity and so, although the precipitation amounts are not high, the plants develop rapidly (Figure 9a). In fact, the correlation is weak for this period, and also for the 1- and 2- composite lagged correlations, supporting the same explanation of the rapid increase in vegetal mass triggered by the temperature increase. Moreover, this is the period when the maize crop, which covers 25% of the arable lands in Romania [99] and is a highly demanding crop in terms of water supply, reaches its maximum vegetal activity. After this period, starting with the second 10-day period of May, the 1-composite lagged correlation increases in significance (more than 30% of the arable lands) with its peak at the beginning of June. This can represent the fact that the precipitation amount reaches its annual maximum during this period. The 2-composite lagged correlation was the most important towards the end of the study period, with its peak between 29 August and 13 of September. In brief, we can observe that the response of crop vegetation to atmospheric precipitation amount is fasterduring the maximum development phase of crops and early summer and slower at the end of crop development.

Overall, at the annual level the 1-composite lagged correlation reached the highest score of frequency of statistical significance on the arable lands, and was close to the direct correlation (Figure 9b). We note that for the annual level the correlation was aggregated for the p-value <0.05.

Therefore, we showed the 1-composite lagged correlation for each composite separately (Figure 10) and for annual conditions (Figure 11). For each composite separately, one can first observe that the significant correlations are mainly negative, reflecting the logical and expected relationship between low/high precipitation amount and drought/excessive humid conditions.

However, some pixels, particularly those located in the south-eastern part of Romania, present a positive significant correlation. This is mainly a result of the vegetation development despite the weak amount of atmospheric precipitation, especially when the plants benefit from a consistent reserve of soil humidity. In addition, the correlation is thoroughly influenced by the crop phenology and the year-to-year crop rotation.

The annual conditions (Figure 11) resume the overall negative correlation between precipitation amount and NDDI values. Interestingly, this negative correlation is stronger over the arable lands in regions with higher precipitation, such as the regions in the northern half of the country. This is induced mainly by the prevalence in these regions of crops with a summer maximum in their development that corresponds with the annual peak in atmospheric precipitation, such as maize.

By comparison, the lack of statistical significance is specific for some spots in the south-eastern part of Romania, but also along the Danube and sparsely in the remainder of the country, and correspond mainly with areas with extensive irrigation.

4. Conclusions

Spatial and temporal characteristics of the NDDI were described in an attempt to analyze and monitor the frequency and severity of drought on the arable land in Romania for a relatively long period (2001–2020). The drought assessment was based on the calculation of the relative frequency of NDDI values indicating drought conditions (NDDI >0.5). The analysis focused on the March–September period to cover the crop development from the greening period of vegetation in early spring, to the stage of maximum vegetation and harvest in late summer and early autumn. In addition, the NDDI values were correlated with atmospheric precipitation, with the aim to determine the role of this important driver on drought frequency and severity.

The main conclusion of our study is that the NDDI represents a valuable tool to assess drought from both temporal and spatial perspectives. However, drought assessment using the NDDI should be interpreted with caution because the NDDI is a result of both meteorological and non-meteorological conditions, which are constantly changing over time. This aspect is underlined, in particular, by the weak significant correlation that we obtained between the NDDI and precipitation amount at the country scale. Moreover, when used on longer time period, as in our study, the NDDI offers a comprehensive view of the degree of aridity for the analyzed region.

Regarding the spatial distribution of drought using NDDI drought classes, it was observed that the most affected regions are the eastern and southern parts of the Romanian Plain, the entire Dobrogea region, and the north-eastern part of the Moldavian Plateau. In addition, most of the plain in the west, and the central and western parts of the Romanian Plain, are affected by moderate drought. Overall, we highlight that these are the arable regions in Romania for which drought should not be considered an extreme event, but a main climate feature, that can be handled by permanent irrigation.

Our results also do not indicate a clear trend regarding the multiannual frequency of drought, as expressed by NDDI values, during the analyzed period between 2001 and 2020. Therefore, based on our results, the discussion regarding a possible desertification manifesting in some parts of Romania is doubtful, at least from the point of view of vegetation conditions.