Decadal assessment of agricultural drought in the context of land use land cover change using MODIS multivariate spectral index time-series data

ABSTRACT

Using a multivariate drought index that incorporates important environmental variables and is suitable for a specific geographical region is essential to fully understanding the pattern and impacts of drought severity. This study applied feature scaling algorithms to MODIS time-series imagery to develop an integrated Multivariate Drought Index (iMDI). The iMDI incorporates the vegetation condition index (VCI), the temperature condition index (TCI), and the evaporative stress index (ESI). The 54,474 km² Vietnamese Central Highlands region, which has been significantly affected by drought severity for several decades, was selected as a test site to assess the feasibility of the iMDI. Spearman correlation between the iMDI and other commonly used spectral drought indices (i.e. the Drought Severity Index (DSI–12) and the annual Vegetation Health Index (VHI–12)) and ground-based drought indices (i.e. the Standardized Precipitation Index (SPI–12) and the Reconnaissance Drought Index (RDI–12)) was employed to evaluate performance of the proposed drought index. Pixel-based linear regression together with clustering models of the iMDI time-series was applied to characterize the spatiotemporal pattern of drought from 2001 to 2020. In addition, a persistent area of LULC types (i.e. forests, croplands, and shrubland) during the 2001–2020 period was used to understand drought variation in relation to LULC. Results suggested that the iMDI outperformed the other spectral drought indices (r > 0.6; p < 0.005). The analysis revealed an increase in drought risk in some provinces of the Central Highlands including Gia Lai, Kon Tum, and Dak Lak. It was also found that changes in LULC patterns could minimize (reforestation) or exacerbate (deforestation) the impacts of drought. Our study suggests that applying a multivariate drought index enables a better understanding of drought patterns at the local scale. This provides valuable information for the development of appropriate land and environmental management practices that can affect and mitigate climate change effects.

KEYWORDS:

• Evapotranspiration
• integrated multivariate drought index
• drought severity index
• vegetation health index
• ground standardized precipitation index
• ground reconnaissance drought index

1. Introduction

Drought is a natural disaster, which can cause serious impacts on human health and well-being, economy, and the environment (Sherval, Askew, and McGuirk 2014; NDMC 2022; WMO 2022). The phenomenon occurs when rainfall is lower than average over an extended period (WMO 2016). This process results in significant water shortages, including the supply of surface water and groundwater. Drought severity has manifold impacts on the environment and socio-economic development in more than 130 countries worldwide (Son et al. 2021; UNDRR 2021). The United Nations estimated that about 55 million people worldwide are affected by severe drought and predicted worsening droughts could displace up to 700 million people by 2030 (UNDRR 2021). They also reported that the potential impact of drought on people’s livelihoods might be comparable with the damage caused by the COVID-19 pandemic (UNDRR 2021). Of those negative effects of drought, agricultural production, especially crop and livestock sectors, are the most affected groups. For instance, loss of agricultural productivity originated from drought occupies almost 86% of the total impacts of climate-related disasters (Wu, Ma, and Yan 2020; FAO, F. and A.O. of the U.N. 2021). Hence, understanding drought severity is valuable for the management of drought-related practice.

Drought can be classified into different types, subject to the causes and effects of drought, geographical regions, or field of study (WMO and GWP 2016; West, Quinn, and Horswell 2019; Wu, Ma, and Yan 2020). Thornthwaite (1947), for example, categorized droughts into four types: permanent (water demand is always higher than supplied precipitation), seasonal (rainy and dry), contingent (irregular and changeable rainfall), and invisible (occurs at any time). Wilhite and Glantz (1985), on the other hand, classified drought into four groups: meteorological drought, hydrological drought, agricultural drought, and socio-economic drought. In addition, several studies suggested other types of drought such as ecological drought (Crausbay et al. 2017; Bradford et al. 2020) and flash drought (Christian et al. 2019; Hu et al. 2021). In practice, all drought types are connected to each other and it is hard to identify the onset and end of each drought type (West, Quinn, and Horswell 2019; Son et al. 2021; Ha et al. 2022). For instance, a decline in precipitation over time is the primary factor affecting meteorological drought. If the precipitation deficiencies persist, the stream flow, lake, and groundwater levels will be reduced, which will lead to hydrological drought. Soil moisture deficits are caused by low precipitation and water supply and are associated with agricultural drought. Finally, a socioeconomic drought occurs when water demand for socioeconomic activity exceeds water supply because of a weather-related deficit, including the impacts of meteorological, agricultural, or hydrological droughts (Elhag and Zhang 2018).

Conventionally, ground meteorological data were collected to monitor and predict drought events. However, these ground-based stations are often sparsely distributed making regional drought monitoring challenging (Jiao, Wang, and McCabe 2021; Ha et al. 2022). In the recent years, remote sensing has been heavily utilized in predicting drought severity improving the drawbacks of traditional approaches in relation to cost, time, and spatial distribution (West, Quinn, and Horswell 2019; Jiao, Wang, and McCabe 2021). Remote sensing that is capable of providing multiple spatial and temporal resolution of data acquiring from various sensors at different spectral regions (wavelengths), has significantly contributed to drought studies over the previous several decades (West, Quinn, and Horswell 2019; Jiao, Wang, and McCabe 2021). Remotely sensed products of drought (mostly using drought indices) have commonly utilized information of land surface temperature (LST), rainfall, vegetation greenness, evapotranspiration (ET), and soil moisture (WMO and GWP 2016; Jiao, Wang, and McCabe 2021) [see comprehensive reviews by Jiao, Wang, and McCabe (2021), Mu et al. (2013), West, Quinn, and Horswell (2019), Wu, Ma, and Yan (2020) and Yu et al. (2019)]. According to these studies, each remotely sensed drought index is generated to investigate a specific type of drought based on its source, which may be utilized as input elements. For other types of drought, meteorological drought often utilizes climatic data (e.g. precipitation, air temperature) (Spinoni et al. 2019; Yu et al. 2019), whereas hydrological drought employs streamflow or groundwater data (Myronidis et al. 2018; Jiao, Wang, and McCabe 2021). Furthermore, vegetation condition or evapotranspiration are commonly used to characterize drought patterns in agricultural drought studies (West, Quinn, and Horswell 2019; Sun et al. 2020). Although drought indices have been produced for several decades, the majority of them are region-specific, and so these indices may not be relevant to covering various geographical settings and environmental conditions such as climate zones, land cover types, and biome types (Hao and AghaKouchak 2013; Son et al. 2021).

Several integrated drought indices have been developed to better monitor and characterize drought conditions (Hao, Zhang, and Yao 2015; Cao et al. 2019; Shah and Mishra 2020). The Palmer Drought Severity Index (PDSI) may be the first multivariate indicator that combined precipitation, temperature, and soil moisture data to detect agricultural drought (defined as soil moisture deficit) patterns (Guttman 1998; Mu et al. 2013), but its parameters heavily depend on regional geographical features (Guttman 1998; Paulo and Pereira 2006; Yu et al. 2019). Kogan (1995) proposed Vegetation Health Index (VHI) that combines Vegetation Condition Index (VCI) and Temperature Condition Index (TCI) to examine the dryness condition associated with the water-energy balance of the land surface. However, utilizing VHI to examine drought conditions may not be appropriate due to the difference in spatial distribution, contribution, and intensity of moisture and temperature to vegetation phenology, which are major factors influencing the weight of VCI and TCI (Karnieli et al. 2006; Bento et al. 2020; Sun et al. 2020). Additionally, VHI does not directly respond to the evapotranspiration process, which is one of the major factors influencing drought conditions on cultivated land in semi-arid areas (Bento et al. 2020; Zeng et al. 2022). Mu et al. (2013) developed the global-scale Drought Severity Index (DSI), a well-known index, using standardized combinations of the MOD16 ET/PET (Evapotranspiration/potential ET) and the MOD13 NDVI (the Normalized Difference Vegetation Index) products of MODIS (the Moderate Resolution Imaging Spectroradiometer) to monitor drought conditions. The ratio of actual ET and PET behaved similarly to precipitation-based indices (Anderson et al. 2011) using evaporation from canopy interception,¹ wet and moist soil surfaces, and transpiration from canopy stomata (Mu et al. 2007; Mu, Zhao, and Running 2011; Mu et al. 2013). However, DSI does not consider the use of land surface temperature (LST) which may lead to missed information of the land surface balance of water, energy, and carbon dioxide. To fill this gap, Tran et al., (2019) added LST to Enhanced Drought Severity Index (EDSI) and applied it to monitor drought at regional scale of the Mekong River Delta. In this study, we integrated additional potentially related factors such as vegetation, temperature, and evapotranspiration to better quantify and understand drought patterns.

The use of time series remotely sensed data may provide an additional temporal dimension to further decipher the complex drought pattern (Frazier et al. 2018; Jönsson et al. 2018; Ghil et al. 2002). However, the differences in ranges of features (i.e. the difference between the lowest and highest values) within an input time-series dataset are relatively high and, in most cases, data used to develop drought indices are measured in different units (e.g. millimeter versus degree Celsius) (Heo and FitzHugh 2000; Lokman et al. 2019; Zakaria 2019). As a result, feature scaling (e.g. data normalization) – a method used to normalize the range of independent variables or features of data – is widely employed as a pre-processing stage to standardize or improve the quality of the input dataset (Heo and FitzHugh 2000; Lokman et al. 2019). Among several feature scaling techniques that have been employed such as linear regression, logistic regression, or neural networks, normalization and standardization are commonly used methods for drought monitoring and assessment (McKee, Doesken, and Kleist 1993; Kogan 1995; Rhee, Im, and Carbone 2010; Mu et al. 2013; Hao, Zhang, and Yao 2015; Tran et al. 2019). The primary difference between normalization and standardization is that the standardized data follows a Gaussian distribution (Aniruddha 2020). Also, in multivariate linear regression, all features need to ensure having similar ranges of values for standardized performance (Dzierżak 2019; Wan 2019). However, for DSI (Mu et al. 2013) and EDSI (Tran et al. 2019), there is a difference in range between indicators of NDVI and ET/PET ([−1; 1] and [0;1], respectively), while maximum LST value is not defined. Therefore, applying feature scaling algorithms is necessary to have a similar range before the parameter integration.

Our study aims to develop an integrated multivariate drought index (iMDI) from MODIS data for monitoring agricultural drought. The iMDI employed a combination of vegetation condition, LST, and evapotranspiration based on feature scaling algorithms (i.e. normalization and standardization). The normalized approach helps to homogenize the range of variables between 0 and 1 before combining these parameters (Shumway and Stoffer 2017; Aniruddha 2020). The standardized method assists in rescaling the integrated values such that they are centered around the mean with a unit standard deviation, as in a normal distribution (Shumway and Stoffer 2017; Aniruddha 2020). The approach allows filling data gaps in the proposed drought indices of DSI (Mu et al. 2013) and EDSI (Tran et al. 2019) with respect to variable min – max values. The combination of three previously published drought indices (i.e. VCI, TCI, and ESI) enables us to understand vegetation and moisture conditions in relation to water stress conditions. As such, it will be an effective index for agricultural drought monitoring. The Central Highlands, one of the major agricultural regions of Vietnam, was selected as a case study to demonstrate how the proposed drought index can be applied to quantify the pattern of drought severity over time (i.e. 2001–2020) (Nguyen et al. 2016; Dinh et al. 2017).

2. Materials and methods

2.1. Study site

The Central Highlands in Vietnam (11°12”00”–15°27”15” N and 107°12–108°59’37” E) consists of five administrative provinces including Kon Tum, Gia Lai, Dak Lak, Dak Nong, and Lam Dong with an area of 54,474 km² (Figure 1). It is considered as a “green lung” of Vietnam with a total forested area of nearly 26,000 km², which is 17.5% of the country’s forested area (Trieu, Pham, and Dao 2021). The study region is comprised of several plateaus with elevation ranging from 500 m to 1,500 m. These plateaus are covered by basaltic soils and bordered by high mountain ranges and massifs to the east (known as the South Annamite Range). The study area has a typical tropical monsoon climate with a rainy (May-October) and a dry (November-April) season (Figure 2). The study region receives 90% of rainfall in the rainy season with a minimum of 100 mm/month. Severe drought events commonly occur in the dry season. Mean annual temperature is about 20 ºC with low seasonality for tropical region. Nevertheless, due to topography, the climate are characterized into three types of sub-regions, which are North Highlands, Central Highlands, and South Highlands (Figure 3) (Toan and Dac 1975).

Figure 1. (A) Location of the study area, with (b) elevation, provincial boundaries, and meteorological stations. The source of the digital elevation model is from ASTER Global Digital Elevation Model.

Figure 2. Monthly mean temperature (T, the primary y-axis) and precipitation (P, the secondary y-axis) acquired from five meteorological stations in the study area from 50 years of climate data (1972 –2021) from the National Centre for Hydro-Meteorological Forecasting. Consecutive months receiving 100+ mm of precipitation (above the standard lines) is defined as the rainy season.

Figure 3. The flowchart of integrated multivariate drought index computation based on MODIS data.

Coffee plantation occupies 85% of the cultivated area of the Central Highlands with the remaining crops being rubber (5.1%), pepper (1.4%), and others (8.5%). Coffee plantations in this area account for 90% of the production for the entire Vietnam (GSO 2021). The flowering season of coffee trees begins in March requiring sufficient water for producing seeds, but rainfall could only provide 25% of the demand and additional water supply from irrigation is sometimes necessary (Amarasinghe et al. 2015). However, deforestation created 20% decline of natural forests between 1990 and 2019 (from 60% to less than 40%) (Chuyen 2020; MARD, M. of A. and R.D. 2021), leading to an increase in surface runoff, reduction in the infiltration of rainwater, and decrease in the amount of stored groundwater to provide irrigation water for crops in the dry season (Chuyen 2020). Consequently, drought has severe impacts on water supply for coffee, rice and other cash crops (Dinh et al. 2017; Nguyen et al. 2017; Vu, Vu, and Trinh 2019). However, the limitation on the number of meteorological stations, where each province only has one station, has created a challenge for broadly monitoring drought in the Highlands. Therefore, to overcome these issues in the region, there is an emerging need for large-scale, long-term, and cost-effective monitoring and mapping techniques that employ remotely sensed data that are available at no or low cost.

2.2. Data preparation

The monthly and annual time-series MODIS LST (MOD11A1), NDVI (MOD09 GA) and ET/PET (MOD16A2) data acquired between 2001 and 2020, were obtained via Google Earth Engine (GEE). GEE, a cloud-based platform, is commonly used to perform deep learning approaches in image analysis with the “JavaScript API (application programming interface) and an online Integrated Development Environment (IDE) code editor” (Gorelick et al. 2017). The GEE is a flexible platform which enables users to perform online analysis of satellite images without downloading the images (https://earthengine.google.com/). This platform allows the users to utilize all GEE data and functions to perform high-performance computing analysis to detect changes and monitor trends on the Earth’s surface using big data (e.g. long-term time series satellite images) (Erickson 2014; Gorelick et al. 2017). In this study, we retrieved the monthly mean, maximum, and minimum values of LST, NDVI, ET, and PET in the 2001–2020 period (Figure 3). Afterward, all the parameters were exported at a spatial resolution of 1000 m for further processing (see the workflow in Figure 3). In addition to the satellite data, the concurrent ground precipitation and temperature data, provided by the National Hydrometeorological Center, was collected from five meteorological stations in the study region (Figure 1) to quantify ground-based annual drought (12-month drought) for consideration of the iMDI performance. Notably, the 500-m annual MODIS Land Cover Type 1 Product (MCD12Q1) for 2000–2019 was utilized to monitor land use and land cover change and access its impacts on drought patterns as anthropogenic factors in the study area (Sulla-Menashe and Friedl 2018) (please see description of land classes in Table 1).

Table 1. Description of the land cover classes derived from MODIS product (Sulla-Menashe and Friedl 2018).

Name	MCD12Q1 Land Cover	Description
Evergreen Broadleaf forests	2	Tree cover of evergreen broadleaf and palmate trees is greater than 60% and their canopy are higher than 2 m.
Deciduous Broadleaf forests	4	At least 60% of the tree cover is neither deciduous nor evergreen, with a canopy height above 2 m.
Perennial Croplands	14	40–60% of the area is either cultivated, found in natural trees, or vegetation.
Annual Croplands	12	Cultivated cropland are ≥ 60%.
Shrubs	9	Canopy height is more than 2 m, and tree cover ranges from 10 to 30%.
Residence	13	Impervious surface (i.e. building materials, asphalt, and vehicles) ≥ 30%.
Water Bodies	17	Permanent water bodies (≥ 60%).

2.3. Integrated multivariate drought index (iMDI)

In this study, iMDI is proposed as an extension of VHI through an addition of ESI for agricultural drought characterization (see Figure 1 for the workflow of iMDI), and details follow:

1) The VCI (eq. 1(1) $VC{I}_i=\frac{V_{i,mean-}{V}_{min}}{V_{max}-V_{min}}$(1) ) was obtained from pixel-based normalization of NDVI (Kogan 1995; Elhag and Zhang 2018):

(1) $VC{I}_i=\frac{V_{i,mean-}{V}_{min}}{V_{max}-V_{min}}$(1)

where V_mean, V_max and V_min are the mean, maximum and minimum NDVI of a period (2001–2020) respectively, and I is the year. This vegetation greenness metric accounts for local differences in ecosystem productivity because the drought implications for vegetation cannot be easily detected from NDVI directly (Kogan 1995; Tran et al. 2017). VCI assumes that climate variance influences vegetation vigor. Extreme drought may decrease vegetation growth and result in the lowest NDVI in a multi-year observation in certain climate region and season. Conversely, the highest NDVI represents an optimal climatic condition for plant growth. Therefore, VCI is superior to NDVI in drought stress detection (Tran et al. 2017; Elhag and Zhang 2018). The range of VCI is from 0 to 1, and when there is a drought in a particular year, vegetation vigor is low (low V_i,mean); hence, the numerator of the eq.1(1) $VC{I}_i=\frac{V_{i,mean-}{V}_{min}}{V_{max}-V_{min}}$(1) is low (low VCI_i).

(2) The TCI (eq. 2(2) $TC{I}_i=\frac{T_{max-}{T}_{i,mean}}{T_{max}-T_{min}}$(2) ) is used to assess vegetation stress induced by high temperatures and excessive precipitation (Kogan and Sullivan 1993; Kogan 1995):

(2) $TC{I}_i=\frac{T_{max-}{T}_{i,mean}}{T_{max}-T_{min}}$(2)

where T_mean, T_max and T_min are the mean, maximum and minimum LST of a period (2001–2020) respectively, and i is the year. In relation to the effects on vegetation growth, high LST in the growing season of vegetation presents unfavorable conditions (i.e. drought), whereas low LST mostly shows favorable conditions (Singh, Roy, and Kogan 2003). This demonstrates the opposite effect compared to the VCI, which is generated from NDVI. The range of TCI is limited in [0;1] and when there is a drought in a particular year, temperature in a particular year is high and hence TCI_i is low.

(3) Drought indices involved in evapotranspiration estimation are mainly based on the ratio of ET/PET, also known as Evaporative Stress Index (the ESI, eq. 3(3) $ES{I}_i=\frac{E{T}_i}{PE{T}_i}$(3) ):

(3) $ES{I}_i=\frac{E{T}_i}{PE{T}_i}$(3)

The ET and PET, directly derived from MODIS ET products, has been demonstrated to be effective in several drought studies (Anderson et al. 2011; Mu et al. 2013; Dorjsuren, Liou, and Cheng 2016) and is rigorously validated with daily ET estimations based on flux tower eddy covariance measurements on a worldwide scale (Jung et al. 2010). MODIS ET also includes NASA Soil Moisture Active Passive (SMAP) within its algorithm framework which is effective for agricultural drought monitoring (Running, Mu, and Zhao 2019; Brust et al. 2021). As such, the ESI (ET/PET) can be used as an approximate proxy for soil moisture which can provide a direct interpretation for a water stress condition. Using other methods such as precipitation data or hydrologic modeling is not efficient as it often requires extensive information (e.g. management strategies and water table distribution). In contrast to precipitation indices, we will likely be able to acquire information on where stress is reduced by active water management or other non-precipitation water inputs (Anderson et al. 2007a, 2007b, 2011). In addition, precipitation and soil moisture data that are available on a global scale are too coarse or uncertain for regional analysis, and thus, were not used for iMDI computation in our study (Mu et al. 2013; Kerr et al. 2016; Peng et al. 2017). The range of ESI is limited in [0;1] and when ESI is low, there is likely to be an association with low soil moisture or dryness, with consequential negative impacts on vegetation.

Lastly, (4) we used mathematical functions to combine these parameters as an integrated index to come up with the proposed index (iMDI) (eq. 4(4) $iMD{I}_i=VC{I}_i+TC{I}_i+ES{I}_i$(4) ) since their values range from 0 to 1 without any unit. Afterward, we standardized iMDI as an anomaly to reduce the dissimilarity with missing data of local mean values (eq. 5(5) $Z_{iMDI}=\frac{iMD{I}_i-iMD{I}_{mean}}{{\delta }_{iMDI}}$(5) ).

(4) $iMD{I}_i=VC{I}_i+TC{I}_i+ES{I}_i$(4)

(5) $Z_{iMDI}=\frac{iMD{I}_i-iMD{I}_{mean}}{{\delta }_{iMDI}}$(5)

where ${\delta }_{iMDI}$and iMDI_mean are the iMDI’s standard deviation and average, respectively. The $Z_{iMDI}$score is dimensionless and varies in $\left(-\mathrm{\infty },+\mathrm{\infty }\right)$ for dryness to wetness conditions, respectively. When there is a drought in a particular year all three indices are low, hence the iMDI_i is low. Conditions of drought are shown by values less than zero, while wet conditions are presented by values above zero. Similar to drought indices that were developed in previous studies (e.g. DSI, EDSI, and SPI), the iMDI anomaly is available to detect abnormally wet periods, whereas our study mainly focused on drought detection. To provide a consistent interpretation of the drought assessment, we classified the drought conditions in the study area into five categories that are similar to the classification in the SPI (McKee, Doesken, and Kleist 1993) and SDI (Mu et al. 2013). Specific threshold values for each category were adjusted to meet the local conditions (Table 2).

Table 2. The drought categories based on iMDI dynamic range and relative categories (McKee, Doesken, and Kleist 1993; Mu et al. 2013; Tran et al. 2019).

Category	Description	SPI	SDI	iMDI
D0	Non-drought	More than 0	More than −0.3	More than 0
D1	Near Drought	0.0 to −0.99	−0.3 to −0.89	0.0 to −0.80
D2	Moderate drought	−1 to −1.49	−0.90 to −1.19	−0.80 to −1.10
D3	Severe drought	−1.5 to −1.99	−1.20 to–- 1.49	−1.10 to −1.80
D4	Extreme drought	−2 or less	−1.50 or less	−1.40 or less

The combination of VCI, TCI, and ESI, obtained from LST, NDVI, and ET products offer valuable information on drought effects in connection to local conditions. Both VCI and TCI are efficient for early drought warning because these indices incorporate spectral profiles of vegetation, which can reveal vegetation response to the influence of environmental conditions (Seiler, Kogan, and Sullivan 1998; Ontel et al. 2021; Zeng et al. 2022). The ESI enables the moisture content to be detected (Running, Mu, and Zhao 2019), and this index is also a valuable supplement that can be incredibly useful in areas where rainfall data is sparse or inconsistent (Anderson et al. 2011). Therefore, an integration of three proposed indices provides effective information about vegetation and water content as an efficient indicator of agricultural drought monitoring. Notably, there are two critical aspects to emphasize in this phase: (i) standardization improves minor fluctuations or seasonal variation, and (ii) standardization does not eliminate any existing trend in the time series. Furthermore, it does not change the statistical significance of trends in any sense (de Carvalho Júnior et al. 2015; Orvos et al. 2015). As only five meteorological stations have been used to observe weather conditions in the Central Highlands across an area of 54,474 km², remotely sensed data should be applied for drought monitoring.

2.4. Performance assessment

To assess the performance of iMDI, Spearman correlation analyses were conducted to test the correlation between iMDI and the annual field-derived Standardized Precipitation Index (SPI–12) (McKee, Doesken, and Kleist 1993) and Reconnaissance Drought Index (RDI–12) (Tsakiris and Vangelis 2005; Tsakiris, Pangalou, and Vangelis 2007). These indices are calculated from data measured from weather stations and have been widely used to evaluate the performance of spectral drought indices (Cao et al. (2019), Shah and Mishra (2020), and Wu, Ma, and Yan (2020). The annual mean air temperature and rainfall data from five weather stations in the Central Highlands from 2001 to 2020 were used to calculate the ground-based drought indices (i.e. SPI and RDI). Correlation between ground-based indices and some commonly used multivariate drought indices that are VHI (Kogan 1995) and DSI (Mu et al. 2013) were carried out to evaluate the performance of iMDI compared to other spectral drought indices. Results from the Spearman correlation analysis between remotely sensed drought indices and ground-based drought information from 2001 to 2020 provided a basis to justify the effectiveness of the proposed drought index.

2.4.1. Ground-based drought indices

The SPI (eq.6(6) $SPI=S\frac{t-\left(c_{2}t+c_1\right)t+c_0}{\left[\left(d_3+d_2\right)t+d_1\right]t+1}$(6) ) is based on the standardized approach to rainfall probability over a certain time period (McKee, Doesken, and Kleist 1993). It describes the change in precipitation using the joint ($\mathrm{\Gamma }$) probability distribution and normalizes the rainfall of skewed probability distribution (McKee, Doesken, and Kleist 1993; Tigkas, Vangelis, and Tsakiris 2013; Um, Kim, and Park 2018). The formula for SPI calculation is as follows:

(6) $SPI=S\frac{t-\left(c_{2}t+c_1\right)t+c_0}{\left[\left(d_3+d_2\right)t+d_1\right]t+1}$(6)

In equation 6(6) $SPI=S\frac{t-\left(c_{2}t+c_1\right)t+c_0}{\left[\left(d_3+d_2\right)t+d_1\right]t+1}$(6) , c₀ = 2.515517; c₁ = 0.802853; c₂ = 0.010328; d₁ = 1.432788; d₂ = 0.189269; and d₃ = 0.001308 are the calculation parameters, $t=\sqrt{ln\frac{1}{G{(x)}^2}}$, where $G(x)$ is the rainfall distribution probability associated to the $\mathrm{\Gamma }$ function; $x$ presents the precipitation sample; and S indicates the positive and negative probability density coefficient. When $G(x)$is more than 0.5 and S is equal to 1, or when $G(x)$ is equal and less than 0.5 and S is equal to-1, the probability density integral formula of $\mathrm{\Gamma }$ distribution function is used to determine $G(x)$(eq.7(7) $G(x)=\frac{1}{{\beta }^{\gamma }\mathrm{\Gamma }(\mathrm{x})}\underset{0}{\overset{x}{\int }}{x}^{\gamma -1}{e}^{\frac{-x}{\beta }}dx,\hspace{0.17em}x>0$(7) ):

(7) $G(x)=\frac{1}{{\beta }^{\gamma }\mathrm{\Gamma }(\mathrm{x})}\underset{0}{\overset{x}{\int }}{x}^{\gamma -1}{e}^{\frac{-x}{\beta }}dx,\hspace{0.17em}x>0$(7)

where, γ and β are the shape and scale parameters, respectively, of the $\mathrm{\Gamma }$ distribution function.

The RDI (eq.8(8) $RD{I}_i=\frac{\sum _{j=1}^k{R}_{\mathrm{i}\mathrm{j}}}{\sum _{j=1}^{k}PE{T}_{ij}},i=1(1)N\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{j}=1(1)\mathrm{k}$(8) ) is conceptualized as a water deficit, representing the water balance deficit between input (precipitation) and output (reference evapotranspiration) (Tsakiris and Vangelis 2005).

(8) $RD{I}_i=\frac{\sum _{j=1}^k{R}_{\mathrm{i}\mathrm{j}}}{\sum _{j=1}^{k}PE{T}_{ij}},i=1(1)N\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{j}=1(1)\mathrm{k}$(8)

(9) $PET=\hspace{0.17em}0.0023\ast \mathrm{R}\mathrm{A}\ast \left(\mathrm{T}-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}\hspace{0.17em}+\hspace{0.17em}17.8\right)\ast \mathrm{T}\mathrm{D}0.5$(9)

where R_ij and PET_ij are the rainfall and potential evapotranspiration (eq. 8(8) $RD{I}_i=\frac{\sum _{j=1}^k{R}_{\mathrm{i}\mathrm{j}}}{\sum _{j=1}^{k}PE{T}_{ij}},i=1(1)N\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{j}=1(1)\mathrm{k}$(8) ) of the j-th month of the i-th year, respectively; k and N is the total number of months and years of the available data, respectively. The PET in the equation 8(8) $RD{I}_i=\frac{\sum _{j=1}^k{R}_{\mathrm{i}\mathrm{j}}}{\sum _{j=1}^{k}PE{T}_{ij}},i=1(1)N\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{j}=1(1)\mathrm{k}$(8) was obtained using Hargreaves equation (eq. 9(9) $PET=\hspace{0.17em}0.0023\ast \mathrm{R}\mathrm{A}\ast \left(\mathrm{T}-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}\hspace{0.17em}+\hspace{0.17em}17.8\right)\ast \mathrm{T}\mathrm{D}0.5$(9) ) (Tigkas, Vangelis, and Tsakiris 2013, 2015) that requires monthly average geo-datasets of mean temperature (T – mean, degree Celsius), daily temperature range (TD, degree Celsius), and extra-terrestrial radiation (RA, radiation on top of atmosphere expressed in mm/month as equivalent of evaporation). In this study, the annual field drought indices (i.e. SPI–12 and RDI–12) and PET (eq.9(9) $PET=\hspace{0.17em}0.0023\ast \mathrm{R}\mathrm{A}\ast \left(\mathrm{T}-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}\hspace{0.17em}+\hspace{0.17em}17.8\right)\ast \mathrm{T}\mathrm{D}0.5$(9) ) from 2001 to 2020 were measured using DrinC software (Tigkas, Vangelis, and Tsakiris 2013, 2015).

2.4.2. Remotely sensed drought indices

The VHI (eq. 10(10) $VH{I}_i=\alpha VC{I}_i+\left(1-\alpha \right)TC{I}_i$(10) ) (Kogan 1995) and DSI (eqs. 11(11) $Z_i=\frac{R{T}_i-R{T}_{mean}}{{\delta }_{RT}}+\frac{V_i-V_{mean}}{{\delta }_V}$(11) and 12(12) $d_{euc}(x,y)=\sqrt{\sum _{i=1}^{20}{(x_i-y_i)}^2}$(12) ) (Mu et al. 2013) were calculated to provide information for a comparison between iMDI and these spectral drought indices. The equations for these indices are presented as follows:

(10) $VH{I}_i=\alpha VC{I}_i+\left(1-\alpha \right)TC{I}_i$(10)

where α is the weighted coefficient between VCI and TCI. The value of α varies according to moisture and temperature conditions. In the case of uncertain moisture conditions, α is set to 0.5, indicating that VCI and TCI are equally weighted for annual VHI computation (Kogan 1995; Tran et al. 2017). For the annual DSI,

(11) $Z_i=\frac{R{T}_i-R{T}_{mean}}{{\delta }_{RT}}+\frac{V_i-V_{mean}}{{\delta }_V}$(11)

$DSI=\frac{Z_i-Z_{mean}}{{\delta }_Z}$

where Z_i is sum of the RT and NDVI and i is the year; RT_i and V_i are ET/PET and NDVI, respectively; ${\delta }_Z$, ${\delta }_{RT}$, ${\delta }_V$, Z_mean, RT_mean, and V_mean are the standard deviations and means of Z, RT, and NDVI in the 2001–2020 period, respectively.

2.5. Time-series clustering analysis

A time-series clustering analysis consists of two parts: measuring similarity and performing clustering algorithm (Aghabozorgi, Shirkhorshidi, and Wah 2015). Clustering is commonly used in data mining, and details of cluster analyses can be found in Aghabozorgi, Shirkhorshidi, and Wah (2015), Guijo-Rubio et al. (2020), and Zhang et al. (2014). In this study, two commonly utilized clustering algorithms, K-Means and Hierarchical clustering, were applied to group drought datasets of the Central Highlands (Guyet and Nicolas 2016). K-Means clustering is defined as an algorithm that “uses the Euclidean distance to separate n observations into k clusters” (Guijo-Rubio et al. 2020)(eq. 11(11) $Z_i=\frac{R{T}_i-R{T}_{mean}}{{\delta }_{RT}}+\frac{V_i-V_{mean}}{{\delta }_V}$(11) ).

(12) $d_{euc}(x,y)=\sqrt{\sum _{i=1}^{20}{(x_i-y_i)}^2}$(12)

The distance measurements compare time point i of time-series x with the same time point i of time-series y, and thus, both time series must be of equal length (n = 20). The measurement process of K – Means clustering is followed by the steps described below (Grabusts and Borisov 2009): (i) choose an iMDI value for K; (ii) initialize the K cluster centers at random (if necessary); (iii) determine the class memberships of the N objects by assigning them to the nearest cluster center or centroid; (iv) re-estimate the K cluster centers’ (means) by assuming the memberships discovered above are valid; (v) exits if the N objects in the previous iteration is at 95% or if some nominated threshold of members remaining stable is exceeded. Otherwise, return to step (iii) to repeat the process.

In addition to the K-Mean clustering, the L-moments based heterogeneity measure was applied to analyze whether the regions might be considered homogeneous. According to Hosking and Wallis (1997), L – moments summarize the properties or geometries of theoretical probability distributions and observed samples. First, these data are examined in order to identify some problems with the homogeneity of the single time series. Then, 500 homogeneous regions with population parameters equal to the regional average sample L – moment ratios were generated using a four-parameter Kappa distribution-based Monte Carlo simulation approach. Finally, the target region’s characteristics were compared with those of the homogenous simulation region. For the sample and simulated regions, the heterogeneity measures H statistic and V statistic can be defined as follows eqs.12(12) $d_{euc}(x,y)=\sqrt{\sum _{i=1}^{20}{(x_i-y_i)}^2}$(12) and 13(13) $H=\frac{V-{\mu }_v}{{\sigma }_v}$(13) (Yoo et al. 2012):

(13) $H=\frac{V-{\mu }_v}{{\sigma }_v}$(13)

(14) $V_i=\sqrt{\frac{\sum _{i=1}^n{n}_i{({\tau }^i-{\tau }^R)}^2}{\sum _{i=1}^n{n}_i}}$(14)

Where ${\mu }_v$ and ${\hspace{0.17em}}^{\sigma }\mathrm{V}$ are the average and the standard deviation of simulated V_i values, respectively; n_i and ${\tau }^i$ are the record length and the sample L-moment at site i, respectively; and ${\tau }^R$ presents the regional averaged sample L-moment. The L-coefficient of variation, the L-skewness, and the L-kurtosis are used as the L-moment in eq. (12)(12) $d_{euc}(x,y)=\sqrt{\sum _{i=1}^{20}{(x_i-y_i)}^2}$(12) to calculate each homogeneity [H(1), H(2), and H(3)], respectively. A region is regarded as acceptably homogeneous if H ≤ 2, possibly heterogeneous if 2 < H ≤ 3, and certainly heterogeneous if H > 3 (Yoo et al. 2012).

2.6. Trend analysis

A pixel-based time-series trend analysis of iMDI was performed to estimate the drought spatial variation in the study area from 2001 to 2020. The slope of the linear regression line was calculated using eq. 13(13) $H=\frac{V-{\mu }_v}{{\sigma }_v}$(13) .

(15) $a=\frac{n\sum _{i=1}^{n}iMD{I}_{i}Y{}_i-\sum _{i=1}^{n}iMD{I}_i\sum _{i=1}^n{Y}_i}{n\sum _{i=1}^{n}iMD{I^2}_i-{\left[\sum _{i=1}^{n}iMD{I}_i\right]}^2}$(15)

where, a is slope of the regression line; iMDI_i and Y_i represent the drought index as a response variable and year as an explanatory variable, respectively, during the period of 20 years. Afterward, only pixels with values that are statistically significant at the 95% level are investigated using a 2-tailed t-test. Negative values of the slope indicate an environmentally stressed or increasing trend in drought intensity since the drought indices decrease toward negative infinity. Conversely, a slope value greater than 0 suggests a decline in drought level over time. To present the true drought conditions in the study site, we multiplied the resulting drought slope values by minus 1. This enables providing basis information to identify whether each of these measurements has a significant change either positively or negatively (increasing or decreasing trend of drought) over time.

3. Results

3.1. The performance of the iMDI

The results from Spearman correlation analysis between iMDI, DSI, VHI, and in-situ drought measurements (i.e. SPI and RDI) at five weather stations are shown in Table 3. Although all spectral indices are significantly correlated with SPI and RDI at the five weather stations (p < 0.005), r–values obtained from iMDI are significantly higher than that of VHI and DSI. For instance, correlation coefficients in iMDI range between 0.16–0.68, while r-values in VHI and DSI are only between 0.10–0.55 and 0.05–0.67, respectively (Table 3). This information in addition to the highest r-value indicates the performance of iMDI and therefore, demonstrates the suitability of using this index for the examination of drought condition in the study area.

Table 3. The correlation coefficients between in situ indices and remotely sensed indices upon meteorological stations (p <0.005) during the 2001–2020 period.

	Kon Tum		Pleiku		Buon Me Thuot		Dak Nong		Da Lat
	SPI	RDI	SPI	RDI	SPI	RDI	SPI	RDI	SPI	RDI
iMDI	0.40	0.31	0.57	0.53	0.63	0.54	0.61	0.68	0.20	0.16
VHI	0.23	0.08	0.33	0.55	0.17	0.25	0.52	0.52	0.17	0.10
DSI	0.40	0.38	0.05	0.07	0.24	0.35	0.60	0.67	0.21	0.05

3.2. Trend analysis of each cluster

The time-series clustering analysis showed five clusters of drought patterns (i.e. A, B, C, D, and E) in the study area (Figure 4). The values of heterogeneity measures (H1, H2 and H3) among five clusters were less than 2 (Table 4) and the probability density boxplot of each cluster showed a normal distribution with the mean value of iMDI centered around zero (Figure 4b). These information demonstrates that grouping drought into five clusters is appropriate to present the drought pattern in the study area (Zhang et al. 2014).

Figure 4. (A) the time-series clusters of drought patterns (A, B, C, D, and E) and (b) the distribution of iMDI values for each cluster in the study region. Boxes, whiskers, bar line, and crosses are range, maximum and minimum, median, and mean of iMDI values, respectively.

Table 4. Descriptive statistics for average time-series of the iMDI (Figure 4b).

Cluster ID	A	B	C	D	E
Counts	14,563	13,145	12,977	6,775	6,332
Slope	−1.40	7.93	10.44	−8.16	3.77
Median	−0.034	−0.023	0.06	−0.014	−0.002
Max	1.55	1.84	1.49	1.27	1.63
Variance	0.33	0.28	0.38	0.24	0.27
Skewness	−0.15	0.08	−0.46	0.03	−0.03
Kurtosis	0.98	0.93	0.72	−0.28	0.49
p-value	0.1624	0.0000	0.0000	0.0000	0.0002
H1	−0.044	−0.039	0.060	−0.028	−0.003
H2	−0.210	0.198	−0.842	0.093	−0.053
H3	1.767	1.803	1.059	−0.54	0.939
Trend	Not Significant	Increasing	Increasing	Decreasing	Increasing

Figure 4 shows that clusters A and D were mainly concentrated in Gia Lai and Kon Tum provinces. Cluster B was found in Gia Lai and Lam Dong provinces. Whereas clusters C and E were located in the middle and east of the Central Highlands, respectively. Comparing the trend of changes among five clusters (Table 4 and Figure 5), clusters B, C, and E presented an increase in drought value (slope >0), whereas cluster D showed a decrease trend (slope <0). This indicates that the risk of drought declined in B, C, and E regions, but increased in D region. In cluster A, an insignificant trend was investigated since the p-value was more than 0.1 despite the slope being less than zero.

Figure 5. The temporal variation of each cluster (A, B, C, D, and E) of monthly iMDI time-series during the 2001–2020 period.

3.3. Spatiotemporal variations of drought pattern during the 2001–2020 period

The annual patterns showed both decreasing and increasing trends in drought risk, with a dominance of decreasing trend across the study area. However, the change varies across space and season, showing that an increase in drought risk was mainly found in Gia Lai and Kon Tum provinces (Figures 2 and 6a). In the early dry season (November–December), an increase in drought risk was observed in the central part of the study area (e.g. a contiguous region between Kon Tum, Dak Lak, and Gia Lai). This trend expanded toward the east part of the study area in the following months and reached the maximum extent in April (the end of the dry season) (Figures 6b and 7). The increasing trend continued in the east region of the Central Highlands, and this reduced significantly in August (Figure 6b). In contrast, the lowest extent of decreasing trend was recorded in the middle of rainy season (August) and the highest value was discovered in December (Figures 6b and 7).

Figure 6. (A) Annual trend and (b) Monthly trend of drought (iMDI) patterns in the study area.

Figure 7. Monthly variation of affected-drought area trend in area extent over the period 2001–2020.

The annual affected area by drought categories between 2001 and 2020 is shown in Figure 8. It is noted that D0 was not included in the analysis because it is a non-drought category. The temporal change of the iMDI demonstrates that the drought intensity reached a peak in 2004 (Severe Drought) and 2005 (Extreme Drought), which is illustrated by the largest affected area in moderate and above drought categories (D2-D4). The areas influenced by severe and extreme drought (D3-D4) were also high in 2002, 2008, 2011, and 2016. A decrease in drought risk pattern was seen in the periods of 2006–2010 and 2012–2015, and the drought-affected area was the lowest in 2017.

Figure 8. Areal extent by drought categories between 2001 and 2020. Note: D1 = near drought; D2 = moderate drought; D3 = severe drought; D4 = extreme drought.

In general, drought events occur every year, but mostly concentrated in the dry season. The near drought pattern (D1) started from August/September surrounded the South Annamite Range and gradually expanded to other regions in accordance with increasing drought levels (i.e. D2, D3, and D4) in the following months (Figure 9). The drought intensity tended to be higher and more severe toward the end of the dry season (April). During this time, the drought risk is found to be more severe in the areas between Gia Lai and Dak Lak, south of Kon Tum, Dak Nong, and Lam Dong.

Figure 9. Seasonal change of iMDI from 2001 to 2020. The displayed years respond to the peak of affected-drought area in Figure 8 but excluding 2020 as final year of period (2002, 2004, 2005, 2011, 2016, and 2020). Please note that D0 represents a non-drought category.

3.4. The pattern of drought variation through land use and land cover

The LULC map shows the dominant cover of forests, shrubs, and croplands, whilst residential areas and water bodies occupy only a small proportion (Figure 10). From 2001 to 2020, forested land cover (i.e. evergreen, and deciduous broadleaf forests) declined by 684,444 km² (i.e. 657,478 km² of evergreen and 26,966 km² of deciduous), whereas other LULC types that are residence, shrub, and water bodies increased 3,383 km², 495,020 km², and 9,406 km²respectively (Figure 10 and Table 5). In terms of agricultural land, the area of annual croplands increased by 83,451 km², while perennial croplands decreased 538 km² in 2000–2010 period but increased by 93,165 km² in 2011–2019.

Figure 10. The LULC of Central Highlands in 2000, 2010, and 2019.

Table 5. Areas of LULC in the study area between 2001 and 2019, performed from LULC data in Figure 10 (km²).

Land Use	2001	2010	2019
Evergreen Broadleaf Forests	2,350,703	1,910,464	1,693,225
Deciduous Broadleaf Forests	121,823	86,504	94,857
Annual Croplands	466,213	549,664	525,500
Perennial Croplands	452,034	451,496	569,311
Residence	17,805	19,751	21,188
Shrubs	2,047,975	2,435,829	2,542,995
Water Bodies	4,922	7,864	14,328

Although LULC change occurred in most of LULC types, a large area remained stable over time. To exclude the drought trends that were impacted by LULC change, we focused on the persistent LULC areas where there was not having a conversion from one LULC type to another. An intersection of LULC layers of 20 years was employed to consider the persistent LULC category at a resolution of 1 km². Five vegetated LULC classes (i.e. evergreen broadleaf forests, deciduous broadleaf forests, annual croplands, perennial croplands, and shrubs) were considered to assess the variations of annual and seasonal drought index by LULC. The change in drought patterns and its relation to LULC between 2001 and 2020 is presented in Figure 11. Most lines tend to have a negative trend from 2001 to 2006 and then a positive trend afterward. However, an abnormal pattern was found in deciduous broadleaf forests, followed by annual croplands and perennial croplands. The iMDI < − 1.4 (extreme drought, Table 3) was explored in connection to these land uses in 2004 and 2005, while the iMDI >1 was found in 2017. The drought patterns in evergreen broadleaf forest and shrub ranged from −1 and above.

Figure 11. The variations of annual drought index associated with persistent LULC.

The increasing trend in drought-affected areas in deciduous broadleaf forest types was not explored for all months (Figure 12b), while this trend varied across seasons in evergreen broadleaf forests and shrublands classes (Figure 12). For the croplands (i.e. perennial and annual), an increase in drought-affected areas was observed in rainy season, mainly concentrated in July and August, while the decreasing pattern was found in the dry season, especially in December (Figure 12d).

Figure 12. The trends in monthly drought area based on persistent LULC classes between 2001 and 2020.

In addition, a decreasing trend was discovered in an area of 8,131 km² of the evergreen broadleaf forests, followed by shrubs and perennial croplands with a decrease of 7,062 km² and 1,113 km², respectively (Table 6). Whereas the increasing trend in drought-affected areas was mainly seen in annual croplands (171 km²) and shrubs (1,510 km²). Additionally, the difference between the increased and decreased areas was found to be greater in forests and perennial croplands than that of annual croplands and shrubs.

Table 6. Increasing or decreasing trend in annually affected-drought area (km²) for unchanged land use classification in the study area during 2001–2020 period.

Land Use/Land Cover	Increase	Decrease
Evergreen Broadleaf Forests	91	8,131
Deciduous Broadleaf Forests	0	0
Annual Croplands	171	609
Perennial Croplands	84	1,113
Shrubs	1,510	7,062

5. Discussion

Various drought indices have been developed to detect drought intensity and severity worldwide (West, Quinn, and Horswell 2019; Ha et al. 2022). Despite their valuable contributions, the suitability and performance of these indices for regional drought analysis are different and significantly dependent on the climate process, local LULC characteristics, and associated anthropogenic activities (Hao and AghaKouchak 2013; Hao, Zhang, and Yao 2015). Although the integration of multiple factors such as rainfall, soil moisture and temperature, and vegetation greenness is seen as an effective approach to develop a regional drought index, a comprehensive drought index that considers all these factors is neglected in previous studies (Rhee, Im, and Carbone 2010; Yao et al. 2010; Anderson et al. 2011; Hao, Zhang, and Yao 2015; Cao et al. 2019; Sun et al. 2020). Compared to drought indices that were developed and applied in previous studies, the iMDI developed in our study is seen to be more effective as it integrates more factors such as vegetation, temperature, and evapotranspiration to quantify and understand drought patterns at local scale. An integration of three published drought indices (i.e. VCI, TCI, and ESI), obtained from these environmental factors, enables the vegetation, soil moisture, and water deficit conditions to be taken into account. The information derived from LST and NDVI reflects the vegetation conditions in relation to temperature, while ESI presents moisture status and actual water stress conditions (Kogan 1995; Mu et al. 2013; Running, Mu, and Zhao 2019; Tran et al. 2019). The use of normalization algorithms in this study created a standardize measurement (from 0 to 1 scale) for each input factor before combining these to create drought index allows the iMDI to be appropriately used for multivariate linear regression (Wan 2019). Notably, using iMDI and time-series analysis methods enable better understanding of the changes in drought intensity in relation to LULC pattern. For example, a significant reduction in natural drought severity was found to be associated with a rapid expansion of perennial croplands (Table 6).

In the Central Highlands, a few studies applied remotely sensed data to monitor and assess drought intensity (Dinh et al. 2017; Le et al. 2021; Thuong et al. 2022). Dinh et al. (2017) applied the Normalized Difference Drought Index (NDDI) to examine the drought variation in the past three decades in the study area. The NDDI was calculated using NDVI and NDWI that permits more detailed information to be observed in respect to the drought pattern (i.e. higher spatial resolution) than ground-based drought indices (e.g. SPI). However, the aforementioned study only focused on drought analysis in March, while drought intensity should be continuously investigated by season or year. Bushra et al. (2019) suggested that using NDWI is unable to differentiate the stress to vegetation canopies occurred due to drought. Alternatively, Le et al. (2021) and Thuong et al. (2022) used remote sensing- based VHI to examine the drought variation in the study area in a 20-year time period. The VHI has significantly improved the major NDDI drawbacks in relation to applications of NDVI and LST based on the normalized algorithm. However, some limitations of VHI have been described within the introduction section (Bento et al. 2020; Sun et al. 2020; Zeng et al. 2022). More importantly, previous studies in the region have not included LULC, while changes in LULC significantly influence drought conditions.

In this study, we examined both annual and seasonal patterns of the drought conditions and integrated the trend of drought and pattern of LULC in drought-affected areas. The trend and pattern of drought within the LULC types can support the development of appropriate land use management plans to minimize the effects of drought. Our study discovered that evergreen broadleaf forests help to reduce affected-drought areas because trees could store water in rainy seasons and provide water balance in dry seasons. High density forest cover or trees with large canopy enable a significant decrease in heat flux, resulting in the reduction of soil evapotranspiration (Rhee, Im, and Carbone 2010). However, a decrease in forest area (i.e. evergreen broadleaf forests and deciduous broadleaf forests), mostly replaced by croplands, was exposed from 2000 to 2019. Changing from woody vegetation cover to cropland could negatively affect the level of water table in relation to reduction of infiltration processes due to increase in surface flow resultant of deforestation. Hence, converting land use types in the Central Highlands should be of concern for forest protection. Within the deciduous broadleaf forests, the seasonal change in drought pattern is clearly presented because these forests shed their leaves in dry season and grow back in rainy season. This change facilitates explaining the insignificant pattern of cluster A (Figure 4) because this region is mostly covered by shrubs and deciduous broadleaf forest types.

Besides, understanding the long-term trends in drought patterns under the impacts of human activities also provides valuable information to develop appropriate land and environmental management practices which can adapt and mitigate the climate change effects. It has been shown that the drought intensity becomes significantly high in El Niño years (e.g. 2002, 2004), while in the La Nina events (e.g. 2008, 2011), the iMDI value represents non-drought conditions (Figure 9). However, developing irrigation systems to produce crops in the dry season can reduce drought severity, as a result of increasing soil moisture and expanding leaf area (Rey, Holman, and Knox 2017). The drought severity was not significant in the dry season of 2015–2016 despite a historical El Niño event because of a substantial enhancement to the irrigation system for crop cultivation that enabled a greater supply of water in the dry season in the study area (Vu, Vu, and Trinh 2019). In fact, the croplands often consume a large amount of water, especially in the dry season, in order to achieve high productivity. The amount of water supply for perennial croplands (e.g. coffee) is much higher than that of annual croplands (e.g. rice). For example, coffee normally consumes 390–520 litters/root with a cycle of 25–27 days/time in the granulation stage, while paddy rice needs only 175–200 litter/root for the whole growth stage (Hieu and Hao 2000; Tan et al. 2018). Also, coffee trees require adequate soil moisture conditions during the flowering and fruiting stages, so that a significant amount of water provided by irrigation often added to the soil in dry season (Tan et al. 2018). As such, this helps minimize the drought severity in the areas where perennial croplands are cultivated.

Considering that major findings obtained from this study provide spatially explicit information (e.g. the trend and pattern of drought severity are presented on maps), spatial-based management and mitigation strategies can be developed to better adapt to local variations in drought and environmental conditions. The applications mentioned above demonstrate that our study is directly linked to the Goal 13 (climate action) and Goal 15 (life on land) among 17 of UN Sustainable Development Goals adopted in 2015 (General Assembly of the UN 2015; Honeck et al. 2018). Supporting these SDGs is considered as an efficient solution to achieve the long-term sustainable development of the region.

6. Conclusions

This research proposed an integrated multivariate drought index (iMDI) using monthly normalized ESI (ET/PET), LST, and NDVI from the MODIS products during 2001–2020 period. The iMDI synthesized the moisture deficits, soil thermal stress, and vegetation growth status in drought processes and makes it favorable for comprehensive drought monitoring. The results from this study revealed that the influence of drought severity mainly focused on the boundary regions of Gia Lai and Kon Tum provinces. Additionally, changes in drought conditions upon five persistent land use types (i.e. evergreen broadleaf forests, deciduous broadleaf forests, annual croplands, perennial croplands, and shrubs) were explored for 20 years in the study area. An increasing trend in affected-drought areas was discovered in cropland and shrubland areas.

This study demonstrated that the proposed approach that uses multi-source remote sensing data enables the achievement of more accurate drought information. This approach does not only cover the meteorological drought information but also reveals the impacts of drought on agriculture. The natural hazard associated with drought is a complex phenomenon and it is a slow process, not an event, that always begins with a precipitation deficit and progresses to drought conditions if the shortage persists for an extended period of time. Various effects of drought have been discovered, including a reduction in soil moisture, an increase in LST, and a negative impact on vegetation growth. As such, it is important that information related to evapotranspiration, soil, and vegetation is incorporated in the monitoring and assessment of drought condition to better understand this process, we must consider the spheres of evapotranspiration, soil, and vegetation.

Given that necessary MODIS data that are used to develop iMDI are globally available, the approach proposed here can be applied to any other areas in the world. In order to achieve better results for local scale application, we suggest that future studies employ higher spatial resolution data to develop drought indices using a similar approach we developed here. For instance, the continental monthly product for AET available in Australia (Guerschman et al. 2022; McVicar et al. 2022), in addition to Landsat LST and NDVI that can be utilized to develop iMDI at 30 m resolution for this region. In addition, a higher spatial resolution of LULC (e.g. 30 m, 10 m) should be used to analyze the impacts of anthropogenic activities on drought conditions. It is also important that future research considers other physical and socio-economic factors in addition to LULC such as topography, irrigation systems, soil type, household income, population density, and agricultural practices to fully explain the drought variations.

Acronyms

DSI	=	Drought Severity Index
EDSI	=	Enhanced Drought Severity Index
ESI	=	Evaporative Stress Index
ET	=	Actual evapotranspiration
iMDI	=	integrated Multivariate Drought Index
LST	=	Land Surface Temperature
LULC	=	Land Use Land Cover
MODIS	=	MODerate-resolution Imaging Spectroradiometer
NDVI	=	Normalized Different Vegetation Index
PET	=	Potential EvapoTranspiration
RDI	=	Reconnaissance Drought Index
SPI	=	Standardized Precipitation Index
TCI	=	Temperature Condition Index
UNDRR	=	United Nations Office for Disaster Risk Reduction
VCI	=	Vegetation Condition Index
VHI	=	Vegetation Health Index

Author Contributions

T.V.T., D.X.T. conceived the idea of the study; D.B.N., T.V.T. obtained and processed data from GEE for measurement; T.V.T., D.X.T. analyzed the data. T.V.T., D.B, C.H., and S.W.M. wrote manuscript. All authors contributed to the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study (i.e. normalized difference vegetation index, land surface temperature, evapotranspiration, land use and land cover, drought indices, and probability maps, respectively) are obtained from the Google Earth Engine (https://earthengine.google.com/) and available from the corresponding author, TVT, upon reasonable request.

Additional information

Funding

This research is funded by Thu Dau Mot University, Binh Duong Province, Vietnam under grant number DT.21.2-063.

Notes

1. The rainfall that is intercepted by the tree canopy and successively evaporates from the leaves.