Short-Term Forecasting of Satellite-Based Drought Indices Using Their Temporal Patterns and Numerical Model Output

Abstract

Drought forecasting is essential for effectively managing drought-related damage and providing relevant drought information to decision-makers so they can make appropriate decisions in response to drought. Although there have been great efforts in drought-forecasting research, drought forecasting on a short-term scale (up to two weeks) is still difficult. In this research, drought-forecasting models on a short-term scale (8 days) were developed considering the temporal patterns of satellite-based drought indices and numerical model outputs through the synergistic use of convolutional long short term memory (ConvLSTM) and random forest (RF) approaches over a part of East Asia. Two widely used drought indices—Scaled Drought Condition Index (SDCI) and Standardized Precipitation Index (SPI)—were used as target variables. Through the combination of temporal patterns and the upcoming weather conditions (numerical model outputs), the overall performances of drought-forecasting models (ConvLSTM and RF combined) produced competitive results in terms of r (0.90 and 0.93 for validation SDCI and SPI, respectively) and nRMSE (0.11 and 0.08 for validation of SDCI and SPI, respectively). Furthermore, our short-term drought-forecasting model can be effective regardless of drought intensification or alleviation. The proposed drought-forecasting model can be operationally used, providing useful information on upcoming drought conditions with high resolution (0.05°).

1. Introduction

Drought, one of the more extreme natural disasters observed in the world, is caused by complex mechanisms between the land surface, ocean, and atmosphere [1,2,3,4]. Since drought cannot only last for weeks, months, or even years but also develop over large spatial extents, it causes considerable problems, such as the decrease of crop yield, shortage of water, desertification, wildfires, and dust storms [5,6,7]. Many studies have reported that ongoing global warming has increased the frequency of severe drought [8,9,10]. According to Emergency Events Database (EM-DAT) provided by Centre for Research on the Epidemiology of Disasters [11], there were 33 drought events worldwide between 2008 and 2018, creating an economic loss of $18 billion. For these reasons, drought monitoring and forecasting are essential for appropriately managing drought-related damage and providing relevant drought information to decision-makers [10,12]. Drought forecasting plays a particularly vital role in risk management as a comprehensive preparation and mitigation of potential drought-caused damage in a timely manner [5,13,14].

In the United States, there are several drought-forecasting systems including the National Oceanic and Atmospheric Administration (NOAA), Climate Prediction Center’s (CPC’s), Seasonal Drought Outlook (SDO) [15,16], and the National Integrated Drought Information System’s (NIDIS’s) Drought Early Warning Systems (DEWS) (https://www.drought.gov). These systems provide monthly and seasonal drought information that integrates United States Drought Monitor (USDM) data, various satellite observations (e.g., products or reflectance of Advanced Very-High-Resolution Radiometer (AVHRR), and Moderate Resolution Imaging Spectroradiometer (MODIS)), satellite-based drought indices (e.g., Vegetation Drought Response Index (VegDRI)) and climate forecasting models (e.g., Climate Forecast System (CFS)). In South Korea, drought-forecasting information can be obtained from the National Drought Information-Analysis Center (www.drought.go.kr), which is operated by the Ministry of Environment and Korea Water Resources Corporation. The center provides present and upcoming drought conditions up to three months using the standardized precipitation index (SPI), the present conditions of water storage, and the available water capacity of the soil and the water level of dams. China officially provides the present drought conditions using ground measurements and satellite data, while they provide only the forecast of each drought factor individually such as precipitation and El Nino-Southern Oscillation (ENSO) from climate models (Weather China, http://products.weather.com.cn). However, in East Asia, there are no drought forecast systems blending in situ and satellite observations with climate model output, unlike in the United States.

There have been considerable research activities looking at drought forecasts using a number of drought-related variables—hydro-meteorological factors (e.g., precipitation, temperature, evapotranspiration, soil moisture and streamflow), drought indices (e.g., SPI, standardized precipitation evapotranspiration index (SPEI), and evaporative stress index (ESI)), and climate indices (e.g., sea surface temperature, and ENSO)—from ground stations, satellites, and reanalysis or numerical models [17,18]. Among them, drought indices have the advantage of being more readily useable for understanding drought conditions than hydro-meteorological factors [19]. This is because drought indices can reflect dryness anomalies and agricultural or hydrological impacts through using multiple drought factors. For this reason, drought indices have been used in many drought forecast studies.

In previous studies using drought indices, ground-based drought indices—SPI, SPEI, and Palmer Drought Severity Index (PDSI)—have been frequently used for forecasting drought [20,21,22,23,24]. However, those indices can cause high uncertainties over ungauged areas when calculated using spatial interpolation [25,26]. For this reason, some studies have used remote sensing, reanalysis or numerical model data in their drought forecasts. Otkin et al. [27] introduced a new index to forecast drought, named Rapid Change Index (RCI). RCI was developed to show the change in drought conditions using ESI (generated through Atmosphere-Land Exchange Invers (ALEXI) using satellite-based LST), soil moisture (obtained from North American Land Data Assimilation System (NLDAS)), and SPI (obtained from the CPC unified analysis of daily precipitation reports) across the United States. Park and Kim [28] forecasted an area of severe agricultural drought using a Landsat 8 and Shuttle Radar Topography Mission (SRTM)-based soil moisture index (SMI) and topographic data. Park et al. [29] used the time-series of satellite-based drought indices (Scaled Drought Condition Index (SDCI), Microwave Integrated Drought Index (MIDI), and Very Short-term Drought Index (VSDI)) and climate indices (i.e., Madden-Julian Oscillation (MJO)) in order to develop a drought-forecasting model for East Asia. Lorenz et al. [30] conducted drought-forecasting studies using USDM data, the anomaly of drought factors (i.e., precipitation, evapotranspiration and soil moisture), and forecasted drought factors from the Climate Forecasting System (CFS) focusing only on drought intensification with two-, four-, and six-week time periods.

Although there has been an effort to forecast drought by integrating multiple source data, drought forecasts are still challenging. This is due to the inherent complexity of drought and the spatio-temporal variability of drought-related variables associated with the global hydrologic cycle [17]. Various methods have been applied to the development of drought-forecasting models, e.g., regression [31,32,33,34], autoregressive integrated moving average (ARIMA) [35,36], and machine learning [28,29,37,38] models. However, most models have limited performance, especially with short-term (within two weeks) forecasting of drought showing lower predictive skills than long-term (e.g., seasonal) forecasts [29]. This is due to the complexity of the hydrological process related to droughts on the short-term scale (e.g., the relationship between soil moisture, temperature, and evapotranspiration after precipitation events).

Another reason for the inaccuracy of short-term forecasting is that most drought forecast studies on a short-term scale have only considered the time-series pattern of drought [27,35,39,40]. Only a few papers have used the numerical forecast model data and climatic indices to reflect the impact of upcoming weather conditions. Lorenz et al. [30] proved that the use of numerical predictors in the short-term (two weeks in this study) drought forecasting through logistic regression could increase predictive skills when the drought intensified. Park et al. [29] combined the time-series of satellite-derived drought indices and climate indices for drought forecasting at a pentad scale. Although the use of climate indices improved the drought-forecasting skills, the model performance was saturated with r ~ 0.7 regardless of the drought indices used. Therefore, it is necessary to improve drought-forecasting skills at the short-term scale considering the alleviation and intensification of drought by combining forecasted climate factors (e.g., precipitation and temperature).

In this research, we aimed to propose a drought-forecasting model on a short-term scale through the integration of numerical model outputs, topographic characteristics (i.e., climate zone, digital elevation model (DEM), and landcover), and satellite-based drought indices (i.e., Scaled Drought Condition Index (SDCI) and SPI) using convolutional long short term memory (ConvLSTM) and random forest (RF) approaches. The assumption of this study is that drought indices (SDCI and SPI for this study) perfectly represent drought. Therefore, if the drought indices are well forecasted, then the drought is also well predicted. The measurable objectives of this study were to 1) forecast drought on a short-term scale (8 days) over a part of East Asia and to 2) analyze the system’s forecasting skills when only considering temporal patterns of drought conditions and when combining numerical model outputs, topographic information, and temporal patterns of drought conditions in the short-term forecasting of drought. The reasons for the 8-day time scale were 1) satellite products are provided every eight days, and 2) drought is not a rapidly changing (e.g., daily) phenomenon. The novelty of this study can be summarized in two aspects: 1) ConvLSTM was used to develop the drought-forecasting model, which has the advantage of being able to predict time-series data by considering spatial characteristics. Because droughts have both spatial and temporal patterns, ConvLSTM can be useful in forecasting drought because it learns time-series and spatial patterns simultaneously. To our knowledge this is the first attempt to use ConvLSTM in a drought-forecasting model. 2) The proposed drought-forecasting model considers not only drought patterns but also predicted climatological factors (i.e., precipitation, temperature) in both drought intensifications and alleviations, while the recent literature has only focused on drought intensification [30].

2. Study Area and Data

2.1. Study Area

The study area is a part of East Asia (latitude: 25.17°N–45.72°N; longitude: 114.05°E–133.25°E), including east China, southeast Russia, Korea, and part of Japan (Figure 1). According to the EM-DAT (https://www.emdat.be), the study area suffered from severe drought 13 times between 2000 and 2018 with economic losses totaling over $20 million. The study area consists of diverse landcover types (i.e., water, forest, cropland, built-up, grassland, and savannas; Figure 1a) and climate zones (i.e., snow, warm temperature, and arid climates; Figure 1c). Cropland is located at low altitudes while grasslands and forests are located at relatively high altitudes. In terms of climatic characteristics (Figure 1c), the study area is generally hot and humid in summer caused by the East Asia monsoon [41], while it is dry and cold in winter. Southern and central-eastern China, Japan, and the coastlines of South Korea have warm temperate climates, while Northeast China, North Korea, and inland South Korea have snow climates.

2.2. Data

Our drought-forecasting model was developed using (1) satellite-based drought indices for documenting the historical patterns and current conditions of drought, (2) numerical forecasting model outputs for considering upcoming weather phenomena, and (3) the spatially distributed geographic characteristics of the study area.

2.2.1. Satellite-Based Drought Indices

In this study, we selected two satellite-based drought indices, SDCI and SPI, for drought forecasting over the study area. The non-vegetated area (i.e., urban area) was extracted when calculating both drought indices because drought indices are valid only over vegetated areas. The irrigation effect was not considered as the irrigation information of the study areas was not provided in detail on the spatial domain. The reasons why SDCI and SPI were selected included that (1) these two indices have been proved useful in drought monitoring in previous studies and Chinese and South Korean drought management (National Bureau of Statistics of China (http://data.stats.gov.cn); National Drought Information Portal (http://www.drought.go.kr/english/), [29,43]), and (2) they are calculated using different drought factors. SDCI is designed to incorporate multiple drought factors (i.e., vegetation health, temperature, and precipitation), while SPI uses a single factor (i.e., precipitation). Thus, they have somewhat different spatial and temporal patterns (Supplementary Figure S1), which can affect forecasting skills when using the same upcoming weather conditions. Note that this study aimed to forecast each drought index, not to compare two indices as we assumed the indices perfectly represent drought.

SDCI [26] is designed to be not only applicable to both arid and humid regions but also flexible in terms of the multiple timescales of precipitation. It combines thermal stress (temperature condition index (TCI), [44]), water stress (precipitation condition index, PCI), and vegetation stress (vegetation condition index, VCI, [45]), which can be used as drought indicators:

$$\mathrm{TCI}=\frac{{\mathrm{LST}}_{\max }-\mathrm{LST}}{{\mathrm{LST}}_{\max }{+\mathrm{LST}}_{\min }}$$

(1)

$$\mathrm{PCI}=\frac{{\mathrm{precipitation}-\mathrm{precipitation}}_{\min }}{{\mathrm{precipitation}}_{\max }{+\mathrm{precipitation}}_{\min }}$$

(2)

$$\mathrm{VCI}=\frac{{\mathrm{NDVI}-\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }{+\mathrm{NDVI}}_{\min }}$$

(3)

$$\mathrm{SDCI}=0.25\times \mathrm{TCI}+0.5\times \mathrm{PCI}+0.25\times \mathrm{VCI}$$

(4)

Similar to USDM, SDCI has six categories: no drought, abnormally dry, moderate drought, severe drought, extreme drought, and exceptional drought (

Table 1. Drought categories based on Scaled Drought Condition Index (SDCI, Rhee et al. [26]) and Standardized Precipitation Index (SPI, McKee et al. [46]).

SDCI Value (unitless)	SDCI Category	SPI Value (unitless)	SPI Category
0.0 to < 0.1	Exceptional Drought	-2.0 and less	Extreme Drought
0.1 to < 0.2	Extreme Drought	-1.99 to < -1.5	Severe Drought
0.2 to < 0.3	Severe Drought	-1.5 to < -1.0	Moderate Drought
0.3 to < 0.4	Moderate Drought	-1.0 to < 0	Mild Drought
0.4 to < 0.5	Abnormally Dry	0 or more	No Drought
0.5 to <= 1	No Drought

, Rhee et al. [26]). To obtain factors for SDCI, we used Terra MODIS products (i.e., MOD 13C1 NDVI and MOD11C2 LST, EARTHDATA (https://search.earthdata.nasa.gov) for VCI and TCI, and Tropical Rainfall Measuring Mission (TRMM) daily precipitation 3B42 (https://disc.gsfc.nasa.gov) for PCI (

Table 2. Summary of data used to develop drought-forecasting model in this study. The entire period of this study is from 2003 to 2018 and the specific period of each variable is given in the table below.

Variables (Collected Period)		Products	Spatial Resolution	Temporal Resolution
Drought indices (from 2003 to 2018)	SDCI	Terra MODIS Normalized Difference Vegetation Index (NDVI, MOD13C1)	0.05°	16 days
		Terra MODIS Land Surface Temperature (LST, MOD11C2)		8 days
	SPI	TRMM precipitation (3B42)	0.25°	daily
Numerical model (from 2015 to 2018)		GFS air temperature	0.5°	3 h (to 240 h)
		GFS precipitation
Static data		Terra MODIS landcover (MCD12C1)	0.05°	yearly
		SRTM digital elevation model (DEM)	90 m	-
		Climate zone (Kottek et al. [42])	0.5°

SPI is one of the ground-based drought indices, which is based on precipitation deficit or surplus with multiple time scales of accumulation [46]. It is calculated by a gamma probability density function using more than 30 years of precipitation data and categorized by eight classes (Table 1). It has the advantage of determining drought types, such as meteorological and agricultural drought, based on the time scales of accumulation [47,48]. However, ground measurement-based SPI does not provide spatially continuous drought information. In this study, similar to PCI, TRMM daily precipitation (3B42) was used for calculating SPI, which was applied to accumulated precipitation, 3-month (SPI3). Although SPI requires a long period (over 30 years) of data, some research has proved the capability of TRMM-based SPI for drought monitoring [48,49,50,51]. SPI was produced in 0.05 degree resolution using bilinear-resampled TRMM precipitation through SPI function [52] in Matlab 2019a.

)). PCIs with 1-, 3-, 6-, 9-month time scales were calculated at 0.05 degree for co-locating with MODIS products.

2.2.2. Numerical Model Outputs

Global Forecasts System (GFS) is a widely used weather forecast model that blends atmosphere, ocean, land/soil, and sea ice models. It provides atmospheric (e.g., air temperature, precipitation) and land-soil variables (e.g., soil moisture) for 240-hour forecasts, four times a day, with 0.5° spatial resolution for the National Centers for Environmental Prediction (NCEP, https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs) (Table 2). In this study, we used averaged air temperature and accumulated precipitation predicted for eight days. The reason for using upcoming weather conditions is that the current drought status may dramatically change depending on the upcoming weather conditions, which cannot be expected only using the historical patterns of drought indices. To minimize the differences in temperature and precipitation data between satellite products and numerical model outputs, linear fitting was conducted between the two (i.e., GFS air temperature and precipitation were converted to LST and TRMM precipitation, respectively). The mean (max and min) slopes and intercepts were 1.28 (4.5 and 0.07) and 9.62 (50 and -6) for precipitation and 1.1 (1.33 and 0.76) and 6 (71 and -91) for temperature, respectively. Although air temperature is not always highly correlated with LST, many studies have proved a relatively strong positive relationship between them [53,54]. Fitted LST and precipitation were normalized to TCI (TCI_GFS) and PCI (PCI_GFS) using the minimum and maximum values when calculating satellite-based TCI and PCI.

2.2.3. Static Data

In this study, landcover, elevation and climate zones were used as additional predictors, considering the environmental and topographic characteristics of the study area. MODIS landcover (MCD12C1, Majority land cover type 1) was used after simplifying classes through reclassification (Figure 1a). For elevation, 90 m SRTM DEM was obtained from the United States Geological Survey (USGS) elevation products site (http://eros.usgs.gov/elevation-products). It was resampled to 5 km using the mean aggregation for co-locating with MODIS products (Figure 1b).

Köppen-Geiger climate classification maps have been frequently used by researchers across a wide range of disciplines for the climatic regionalization of environmental variables [41,42]. The map was developed through rule-based methods using 50 years (from 1950 to 2000) of temperature and precipitation from reanalysis data. The map has 31 classes (e.g., equatorial, arid, warm temperate, snow, and polar) in 0.5 degrees (http://koeppen-geiger.vu-wien.ac.at/present.htm). The study region was cropped from the map produced by Kottek et al. [42] and converted to 0.05 degrees so that it matches MODIS datasets using the nearest neighbor method (Figure 1c).

3. Methodology

Figure 2 shows the process flow diagram of this study. The proposed approach is divided into two steps: step 1 uses ConvLSTM to obtain temporal patterns from historical drought conditions (i.e., SPI and SDCI) and step 2 uses RF to feed static variables (i.e., landcover, elevation and climate zone) and forecasted climate factors (i.e., temperature (TCI_GFS) and precipitation (PCI_GFS)) provided from the numerical model into the output of step 1. There are two reasons why the final drought-forecasting model combines two machine learning approaches. First, the model structure becomes complex when all independent variables are used as input variables for ConvLSTM. This requires much more memory and processing time than statistical and basic machine learning (e.g., support vector regression or RF) approaches [55]. Therefore, the temporal patterns of each drought index and spatial information were used in ConvLSTM and RF, respectively, to enable drought forecasting even in a memory-limited environment. The second reason is that ConvLSTM is a model optimized for analyzing temporal patterns such as the drought indices used as predictors in a spatial context in this study. Numerical model outputs and static data were used to help improve the model’s forecasting skills as they can provide information which was not included in the temporal patterns of the drought indices (i.e., SDCI and SPI).

During the entire study period, the ConvLSTM model was produced by real-time learning (step 1, 2003–2018 (2003–2014 for optimizing parameters)) to reflect the most recent drought condition. The RF model was applied to produce the final drought forecasts using data from 2015 to 2017, using the output of ConLSTM, static data, and numerical model outputs. The model performance was evaluated for 2018. The reason why study periods were divided into 2003–2018 and 2015–2017 is to reserve enough data for the tuning parameters of ConvLSTM and to obtain the number of reasonable samples required for the models (i.e., ConvLSTM and RF) considering the annual phenology of drought factors (i.e., temperature, vegetation, and precipitation). After obtaining drought indices (i.e., SDCI and SPI) as an 8-day interval, the ConvLSTM model for each drought index (SDCI model and ConvLSTM-SPI model for step 1) was generated by training the previous one-month (8-day x 4) of drought conditions in step 1. In step 2, the final drought-forecasting model for each drought index (SDCI model and SPI model for step 2) was developed using forecasted weather conditions, climate zone, landcover, DEM, and outputs of step 1 through RF, which reduces errors caused by training only the temporal patterns of the drought indices.

3.1. Step 1: Convolutional Long Short Term Memory (ConvLSTM)

Shi et al. [56] developed ConvLSTM (Figure 3) which combines convolutional neural networks (CNN) and long short term memory (LSTM). The two algorithms have been widely used in image classification and time-series forecasting, respectively. ConvLSTM has recently been applied to various research tasks that need to consider both time-series patterns and spatial information, such as segmentation, change detection, forecasting video frames, forecasting sea surface temperature, and air pollution research [57,58,59,60,61,62]. ConvLSTM models space-time structures through encoding spatial information, which can overcome the major limitations of LSTM, namely the loss of spatial information [63]. The structure of ConvLSTM is similar to that of LSTM, which consists of memory cells and three gates (i.e., forget, input, and output gates). The three gates play roles in maintaining or discarding memory (temporal information) [64]. The main difference between LSTM and ConvLSTM is that the internal matrix multiplications of LSTM are replaced with convolution operations in ConvLSTM. As a result, ConvLSTM can produce 2-D output while the result of LSTM is only a 1-D vector. The main equations are as follows [56]:

$${\mathrm{i}}_{\mathrm{t}}=\mathsf{\sigma }\left({\mathrm{W}}_{\mathrm{xi}}{*\mathrm{X}}_{\mathrm{t}}{+\mathrm{W}}_{\mathrm{hi}}{*\mathrm{H}}_{\mathrm{t}-1}{+\mathrm{W}}_{\mathrm{ci}}\circ {\mathrm{C}}_{\mathrm{t}-1}{+\mathrm{b}}_{\mathrm{i}}\right)$$

(5)

$${\mathrm{f}}_{\mathrm{t}}=\mathsf{\sigma }\left({\mathrm{W}}_{\mathrm{xf}}{*\mathrm{X}}_{\mathrm{t}}{+\mathrm{W}}_{\mathrm{hf}}{*\mathrm{H}}_{\mathrm{t}-1}{+\mathrm{W}}_{\mathrm{cf}}\circ {\mathrm{C}}_{\mathrm{t}-1}{+\mathrm{b}}_{\mathrm{f}}\right)6)$$

(6)

$${\mathrm{C}}_{\mathrm{t}}{=\mathrm{f}}_{\mathrm{t}}\circ {\mathrm{C}}_{\mathrm{t}-1}{+\mathrm{i}}_{\mathrm{t}}\circ \tan \mathrm{h}\left({\mathrm{W}}_{\mathrm{xc}}{*\mathrm{x}}_{\mathrm{t}}{+\mathrm{W}}_{\mathrm{hc}}{*\mathrm{H}}_{\mathrm{t}-1}{+\mathrm{b}}_{\mathrm{c}}\right)$$

(7)

$${\mathrm{O}}_{\mathrm{t}}=\mathsf{\sigma }\left({\mathrm{W}}_{\mathrm{xo}}{*\mathrm{X}}_{\mathrm{t}}{+\mathrm{W}}_{\mathrm{ho}}{*\mathrm{H}}_{\mathrm{t}-1}{+\mathrm{W}}_{\mathrm{co}}\circ {\mathrm{C}}_{\mathrm{t}}{+\mathrm{b}}_0\right)$$

(8)

$${\mathrm{H}}_{\mathrm{t}}{=\mathrm{O}}_{\mathrm{t}}\circ \tan \mathrm{h}\left({\mathrm{C}}_{\mathrm{t}}\right)$$

(9)

${\mathrm{i}}_{\mathrm{t}}$ and ${\mathrm{O}}_{\mathrm{t}}$ are input and output gates. ${\mathrm{C}}_{\mathrm{t}}$ is the memory cell that has accumulated the state information using $\mathrm{W}$(weight). When new input is entered, ${\mathrm{i}}_{\mathrm{t}}$ is activated and ${\mathrm{C}}_{\mathrm{t}}$ accumulates the state information. The past status can be forgotten when ${\mathrm{f}}_{\mathrm{t}}$ is on. ${\mathrm{O}}_{\mathrm{t}}$ can be obtained through the final state ${\mathrm{H}}_t$, which is propagated by ${\mathrm{C}}_{\mathrm{t}}$ (Figure 3). A more detailed explanation of ConvLSTM can be found in [56].

Recently, Mu et al. [60] proposed a ConvLSTM-Rolling Mechanism (ConvLSTM-RM, named “real-time learning” in the present study), which utilizes the most recent data (named “stride period” in this study) to develop a forecasting model. The RM method helped improve forecasting performance [60,65,66]. Therefore, in this study, four temporally consecutive data with an 8-day interval that stride for the recent three years were applied to forecast drought conditions in the next 8-day interval through real-time learning from 2003 to 2017. For example, to forecast drought conditions on 9 January 2018, the ConvLSTM model is updated through striding the four consecutive data with the 8-day interval from 9 January 2015 to 1 January 2018 (for three years). The stride period of three years was determined considering the computational efficiency and the impact of annual phenology. The “convlstm2d” function provided by TensorFlow Core r2.0 (https://www.tensorflow.org/api_docs/python/tf/keras/layers/ConvLSTM2D) was used in this study. After testing various combinations of parameters using data from 2003 to 2014, the ConvLSTM structure was determined to have three layers with 16, 16, and 1 filters of 3 x 3, 3 x 3, and 1 x 1 size (Figure 3). Mean square error (MSE) was used as a loss function, and the model was optimized by the adaptive moment estimation (Adam), which is a popular optimization algorithm in neural networks [67]. All experiments were run on a computer with Intel(R) Xeon(R) CPU E5-2680 v2 @ 2.80GHz and NVidia Titan Black GPU (6GB of memory). About 24 days were taken to train the ConvLSTM model in step 1 with 50 epochs and 2 batch size.

3.2. Step 2: Random Forest (RF)

RF is an ensemble approach based on classification and regression trees (CART), which overcomes the major limitations of CART, such as its sensitivity to training data configuration and the overfitting problem, by aggregating multiple independent trees [68]. RF consists of a variety of CARTs (“decision trees”) that have the same probability distributions through bootstrapping-based randomization approaches. All decision trees are aggregated using (weighted) majority voting for classification and (weighted) averaging for regression. RF provides the relative importance of independent variables, which has been widely used in previous studies, even though it is of local importance, not global [14,41,69,70,71,72]. It can be obtained as a percentage of the increased MSE for each variable using out-of-bag (OOB) data. Having a variable with a relatively large percentage of increased MSE means that it made a relatively more significant contribution to our model than to other variables.

In this study, upcoming weather conditions and terrestrial information were used as additional inputs to RF to further improve forecasting skills. The predicted drought conditions from ConvLSTM (output from step 1), static variables (i.e., landcover, climate zone, and DEM), and normalized climate factors (i.e., TCI_GFS and PCI_GFS) were used as independent variables, while the drought conditions in the next eight days were used as a dependent variable in step 2 (Figure 2). From 2015 to 2017, the calibration and validation data set were divided into 80% and 20% respectively after randomly extracting samples to ensure that landcover types, climate zones, and all ranges of drought index values (e.g., 0–1 for SDCI) were uniformly included. The RF package of R statistic software (http://www.r-project.org) was used with default settings except for the number of trees (500).

3.3. Accuracy Assessment

In order to evaluate the performances of the forecasting model for each drought index (i.e., SDCI and SPI), three statistical metrics were used: correlation coefficient (r), normalized root mean square error (nRMSE), and mean absolute scaled error (MASE; Hyndman and Koehler [73]).

$$\mathrm{r}=\frac{\mathrm{n}(\sum \mathrm{y}\hat{\mathrm{y}})-(\sum \mathrm{y})(\sum \hat{\mathrm{y}})}{\sqrt{[\mathrm{n}\sum {\mathrm{y}}^2-{(\sum \mathrm{y})}^2][\mathrm{n}\sum {\hat{\mathrm{y}}}^2-{(\sum \hat{\mathrm{y}})}^2]}}$$

(10)

$$\mathrm{RMSE}=\sqrt{\frac{\sum {(\mathrm{y}-\hat{\mathrm{y}})}^2}{\mathrm{n}}},\mathrm{nRMSE}=\frac{\mathrm{RMSE}}{{\mathrm{y}}_{\max }-{\mathrm{y}}_{\min }}$$

(11)

$$\mathrm{MASE}=\frac{1}{\mathrm{m}}\sum \frac{\left|{\mathrm{y}}_t-\widehat{{\mathrm{y}}_t}\right|}{(\sum \left|{\mathrm{y}}_t-{\mathrm{y}}_{t-1}\right|)/(n-1)}$$

(12)

n and m and the number of samples and $\mathrm{y}$ and $\widehat{\mathrm{y}}$ are the values of reference and predicted drought indices, respectively. The forecasting skills are useful when r (nRMSE) is closer to 1 (0). MASE is one of the statistical indices used to evaluate the time-series forecasting model, which has an advantage when there are very different scales (i.e., negative, positive and zero) [73]. MASE is less than 1 (MASE < 1) when the forecasting error is better than the averaged variation of the time-series data.

4. Results and Discussions

4.1. The Performance of Drought-Forecasting Model

Figure 4 depicts the model performance of the SDCI model for step 1 (ConvLSTM, Figure 4a) and step 2 (RF, Figure 4b,c). Based on the period of the three-year stride, the averaged r and nRMSE are presented from 2006 to 2017 (Figure 4a). The r and nRMSE values ranged from 0.23 to 0.92 (mean of 0.73) and from 0.07 to 0.2 (mean of 0.12), respectively. Despite the generally high accuracy, some dates had low r and high nRMSE due to sudden changes in drought conditions based on SDCI. For example, on 7 April 2006 (Day-Of-Year (DOY) 2006097) in Northeast China (Liaoning and Jilin provinces), drought conditions rapidly changed in eight days from abnormally dry to extreme drought, which was caused by little precipitation and unusually high temperatures. The region remained in normal drought conditions for about two months (DOY: 2006033–2006089), making it difficult to forecast serious drought outbreaks using temporal patterns. In other words, drought forecasting through ConvLSTM cannot reflect the impact of upcoming precipitation and temperature because ConvLSTM only considers historical patterns. Some previous studies that considered historical patterns were similarly vulnerable to sudden weather changes [24,29].

The r, nRMSE, and slope of RF calibration and validation are 0.98 and 0.90, 0.05 and 0.11, and 0.89 and 0.78, respectively (Figure 4b,c). In terms of r and nRMSE, the forecasting skill was improved when integrating upcoming weather and spatial information with ConvLSTM outputs. There was a tendency for an overestimation in exceptional and extreme drought conditions (low values of SDCI) and underestimation in no drought conditions (high values of SDCI). This was probably because 1) the samples in extreme values (0–0.1 and 0.9–1) were smaller (4% (for 0–0.1) and 5% (for 0.9–1)) than the other samples [74]; 2) random forest produces trees with reducing errors, which lead to the values trending around the mean value, especially when there are not many extreme samples [75,76], and 3) the numerical model outputs did not reflect weather well in some regions and on some dates due to low correlation with the satellite products (i.e., GPM precipitation and MODIS LST) (Supplementary Figure S2). However, compared to the validation results on Park et al. [29], who developed RF-based drought-forecasting models with the pentad interval using SDCI and climate index (MJO), our model showed larger dynamic ranges in the output, possibly due to the incorporation of numerical model outputs. Lorenz et al. [30] also demonstrated that drought-forecasting models using weather forecasts (when the drought intensified) performed better than those that only used the past and present drought conditions. The accuracy of the model without using upcoming weather conditions decreased, resulting in r of 0.85 and nRMSE of 0.13 for validation, not shown here), which supports the incorporation of upcoming weather conditions into the model.

Figure 5 shows the performance of the SPI model through step 1 (ConvLSTM, Figure 5a) and step 2 (RF, Figure 5b,c). The r and nRMSEs of the SPI ConvLSTM model ranged from 0.22 to 0.96 (mean of 0.84) and from 0.04 to 0.2 (mean of 0.09), respectively (Figure 5a). Performance varied between SPI- and SDCI-based ConvLSTM model in terms of r and nRMSE (Figure 4). This is probably because of the temporal variabilities according to type and number of drought factors when calculating drought indices. There are three cases: (1) If there is a drastic change in precipitation only, SPI rapidly changes from dry to wet or wet to dry, unlike SDCI that can have lower temporal variability caused by other factors (i.e., temperature or vegetation condition). (2) If there is an abrupt change in temperature, SDCI undergoes drastic changes. However, such a temperature change does not impact SPI because it is a function of precipitation only. (3) If both precipitation and temperature change dramatically, the result depends on the intensity of the change in each factor.

In the RF model, the r, nRMSE, and slope of RF calibration and validation were 0.98 and 0.93, 0.05 and 0.08, and 0.90 and 0.84, respectively (Figure 5b,c). In terms of r and nRMSE, both RF-based SDCI and SPI models showed better results when using both temporal patterns and numerical model outputs compared to when only temporal data were used.

Figure 6 describes the spatial and time-series distribution of the performance metrics of the SDCI and SPI model in 2018 (validation year). High correlation values (mean r of 0.62 and 0.77) for the SDCI and SPI models are shown in most areas in Figure 6a and Figure 6d, respectively. There are common distributions in the SDCI and SPI models (see r graph in middle of Figure 6a,d). Some regions (e.g., Shandong and Jiangxi province) and dates (e.g., DOY 2018193 and 2018337) showed lower correlations than other regions and dates for the SDCI model, which indicates that the numerical model did not have much impact on the improvement of the forecasting skills through RF in those specific cases. This is because there is a significant difference between upcoming precipitation and temperature from the numerical model and the precipitation and temperature in SDCI that were produced by satellite products (Supplementary Figure S2). Besides, upcoming vegetation stress was not considered as an independent variable in the RF model because of its absence in the numerical model, which could cause a decrease in correlations. The spatial distributions and time-series distributions of nRMSE (Figure 6a,e) also show generally low values (means of 0.19 and 0.17 for SDCI and SPI, respectively), which have similar patterns to r.

In Figure 6c,f, MASE shows that the forecasting ability is good (MASE < 1) when the forecasting error is less than the average time-series variation of the drought indices. Unlike the other metrics (i.e., r and nRMSE), the MASE maps and graphs have different distributions, which were caused by considering the time-series variation of each index. While most SDCI MASE values were less than 1 (mean of 0.93), SPI MASE were not (mean of 1.15). The fact that the MASE value is over 1 indicates that the error between forecasts and actual values is greater than the average of the time-series variation at each pixel. There are two reasons why SPI has relatively high MASE. One is a limitation of the MASE metric, which is sensitive to outliers [77]. Since SPI values fluctuate depending on precipitation only, the variation of the SPI is generally greater than that of SDCI. The other is that drought-forecasting skills decrease when drought conditions change rapidly. Therefore, the MASE of the SPI can be higher than SDCI even if the average of the time-series variation in the two drought indices is similar.

4.2. The Spatial Distribution of the Drought-Forecasting Model

Figure 7 depicts spatial distributions of reference SDCI (SDCI_o), outputs of forecasting model (i.e., ConvLSTM (SDCI_s1) and RF (SDCI_s2)), and forecasted climate factors (i.e., precipitation and temperature) from 1 May to 2 June in 2018. The spatial patterns of drought and no-drought conditions were well detected regardless of outputs of SDCI_s1 or SDCI_s2 (e.g., alleviation in the north-western region (i.e., Inner Mongolia) from 1 to 17 May and intensification in Korean Peninsula from 17 May to 2 June). This indicates that the spatial distribution of drought can be forecasted only using historical patterns of drought when droughts gradually intensify or alleviate [78]. However, in the case of sudden droughts, such as in the central- and north-western regions (i.e., Inner Mongolia and Henan Provinces) on DOY 2018045, the forecasting skills were limited, especially in terms of drought intensity. According to the relative importance of the independent variables used in RF, two upcoming weather data have a significant impact on the RF-based drought-forecasting model (after the ConvLSTM output). However, on that day, the quality of the TCI_GFS and PCI_GFS was lower in that 8-day interval, which lowered forecasting skills.

SDCI_s1 was slightly underestimated or overestimated when droughts alleviated or intensified, respectively, whereas SDCI_s2 reduced the differences between SDCI_s1 and SDCI_o using forecasted weather data (i.e., TCI_GFS and PCI_GFS). SDCI_s2 improved 51% of the total pixels on average (up to 68%) in terms of the value of $\left|\mathrm{forecasted}\text{-}\mathrm{actual}\right|$. In other words, the GFS data were used to improve the SDCI_s1 through the RF model. For example, the drought on 17 May was alleviated compared to the eight days before in the central-western region (e.g., Shandong and Hebei Provinces). The output of SDCI_s2 is closer to SDCI_o than SDCI_s1 (the error of 60% of the whole pixel decreased on 17 May). When compared to Lorenz et al. [30], who focused only on drought intensification, our models are applicable for both drought intensification and alleviation.

However, the improvement is relatively low when there are large gaps between the satellite data (i.e., TRMM Precipitation, MODIS LST) and the GFS data (i.e., TCIGFS and PCIGFS). For example, on 9 May, SDCI_s2 improved only about 45% of the pixels due to an underestimation of GFS precipitation in Shandong province (Supplementary Figure S3). Kumar et al. [79] compared GFS and TRMM in terms of monthly precipitation and pointed out that GFS precipitation was overestimated in south and northeast China and underestimated in central-eastern China. Another reason for the discrepancy between SDCI_o and SDCI_s2 is the degradation of forecasting skills in the numerical model as lead time increases [69].

The outputs from our forecasting model for SPI (i.e., ConvLSTM (SPI_s1) and RF (SPI_s2)) are depicted in Figure 8 with reference to SPI (SPI_o) and forecasted precipitation (PCIGFS). SPI and SDCI models showed the spatial patterns. However, some regions have different drought severity because vegetation and temperature stresses were not considered (e.g., the central-western region and Korean Peninsula on 2 June).

Similar to the results of SDCIs, the spatial patterns of drought and no-drought conditions were well detected in SPI_s1 and SPI_s2. The outputs improved when using forecasted precipitation data (PCI_GFS) (up to 62% of pixels). For example, SPI_s1 on 9 May was underestimated in the central-western regions (e.g., Jiangsu and Anhui Provinces), although the wet conditions were well detected. SPI_s2 captured wetter conditions than SPI_s1 using PCI_GFS, which was helpful to improve the forecasting skills (60% of the pixels were improved on 9 May). Another example is that the drought on 25 May was well captured when considering the dry conditions in PCI_GFS, which helped improve about 51% of pixels from SPI_s1 to SPI_s2, especially in the Korean Peninsula. In other words, if there are well-forecasted weather data, the accuracy of the drought-forecasting model can be improved. Lorenz et al. [30] also found that the weather forecasting model is responsible for improving short-term forecasting of drought. In contrast, sometimes the drought-forecasting skills were degraded when integrating forecasted data, due to the discrepancy between satellite products and numerical model outputs (e.g., Shandong Province on 9 May and Jiangxi and Fuzhou Province on 2 June) [79]. This has already been described in Lorenz et al. [30].

4.3. Novelty and Limitations

In this study, we proposed a new drought-forecasting model on a short-term scale using time-series patterns and upcoming weather conditions through ConvLSTM and RF combined. Many previous studies have developed drought-forecasting models aimed at improving the accuracy of forecasting skills (e.g., drought area and intensity) [25,26,27,28,29,30]. However, they still show limitations in forecasting skills: they have relatively simple statistical approaches (e.g., logistic regression and RF) and use historical data. To improve forecasting skills on a short-term scale we combined two models, ConvLSTM and RF, blending temporal patterns of drought and upcoming weather conditions. There are two novelties in our study: (1) the ConvLSTM approach used in this study well reflected both spatial and temporal patterns in drought forecasting. To the best of our knowledge, this is the first study in which ConvLSTM is used in the drought forecasting. (2) Our model fuses time-series patterns and upcoming weather conditions by blending approaches (i.e., ConvLSTM + RF and satellite products + numerical model outputs), which has not been tried much in previous studies. Lorenz et al. [30] also fused past drought conditions and upcoming weather conditions, but they only investigated drought intensifications. The present study, on the other hand, has examined both drought intensifications and alleviations.

However, despite the novelties of the proposed model, there are still some limitations. First, computational demand is a common problem in deep learning-based models [80]. Although ConvLSTM reflects time-series patterns well, it requires a significant computational demand in terms of memory and running time when optimizing parameters because there are 26 parameters in the Keras “convlstm2d” function. Second, machine learning approaches are generally known as black box models which give results that are hard to interpret in terms of the causal relationships between and specific importance of variables. Third, the forecasting skills were not good when there were sudden changes in drought conditions. This is because of the discrepancies between satellite products and numerical model outputs and the degradation of forecasting skills in the numerical model with increasing lead time [69].

To overcome these limitations, several plans can be made for future studies: (1) Auto-parameterization tools (e.g., Keras-tuner and AutoKeras) should be adopted for cost-effective parameterizing; (2) heatmaps should be generated to interpret the effect of each input variable on model performance [81]; (3) other machine or deep learning approaches that can reflect complex drought mechanisms should be tested to further improve drought-forecasting skills; and (4) an ensemble of various numerical models should be tested in order to reduce the gap between the satellite products and numerical model outputs.

5. Conclusions

Short-term forecasting of drought is crucial to reduce the damage to agriculture caused by drought, especially during critical crop yield development stages. Many studies have been conducted for drought forecast, but they still have limited forecasting skills. In this study, a drought-forecasting model on a short-term scale was developed using temporal patterns of drought indices and upcoming weather conditions (numerical model outputs) through the synergistic use of ConvLSTM and RF. Two satellite-based drought indices—SDCI and SPI—were selected with a short-time scale (eight days), and the GFS numerical model was used to improve drought-forecasting skills, considering upcoming weather conditions. The SDCI- and SPI-based drought-forecasting models proposed in this study (ConvLSTM and RF combined) showed competitive results in terms of r (0.90 and 0.93 for validation SDCI and SPI respectively) and RMSE (0.11 and 0.08 for validation of SDCI and SPI, respectively). Furthermore, our drought-forecasting model on a short-term scale can be applicable regardless of drought intensification or alleviation. While ConvLSTM resulted in good performance, the combined model showed better results by feeding upcoming weather conditions and topographic characteristics. The proposed drought-forecasting model can be operationally used, providing useful information on upcoming drought conditions with high resolution (0.05°).