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

: Drought forecasting is essential for e ﬀ ectively 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 e ﬀ orts in drought-forecasting research, drought forecasting on a short-term scale (up to two weeks) is still di ﬃ cult. 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 e ﬀ ective regardless of drought intensiﬁcation or alleviation. The proposed drought-forecasting model can be operationally used, providing useful information on upcoming drought conditions with high resolution (0.05 ◦ ).


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 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].

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.

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.

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

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.

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.
Remote Sens. 2020, 12, 3499 5 of 21 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:  Table 2)). PCIs with 1-, 3-, 6-, 9-month time scales were calculated at 0.05 degree for co-locating with MODIS products. 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.

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.

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). 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). 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.  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.

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 Remote Sens. 2020, 12, 3499 8 of 21 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]: 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]: O t = σ(W xo *X t +W ho *H t-1 +W co ∘C t +b 0 ) i t and O t are input and output gates. is the memory cell that has accumulated the state information using W (weight). When new input is entered, i t is activated and C t accumulates the state information. The past status can be forgotten when f t is on. O t can be obtained through the final state H , which is propagated by C t (Figure 3). A more detailed explanation of ConvLSTM can be found in [56].  i t and O t are input and output gates. C t is the memory cell that has accumulated the state information using W (weight). When new input is entered, i t is activated and C t accumulates the state information. The past status can be forgotten when f t is on. O t can be obtained through the final state H t , which is propagated by C 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.

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).

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]).
n and m and the number of samples and y andŷ 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. , 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 Performance of Drought-Forecasting Model
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.   Figures 5b and 5c). 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 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.    Figure 6a and d, respectively. There are common distributions in the SDCI and SPI models (see r graph in middle of Figures 6a and 6d). 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 (Figures 6a and 6e) 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 and 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 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 Figures 6a and 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. 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 timeseries 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. Figure 6. Spatial distribution and time-series patterns of r, nRMSE, and MASE from random forest model results for SDCI and SPI. The vivid red in the six maps indicates areas of relatively high errors (a-f). The green and purple in three time-series graphs indicate the time-series patterns of SDCI and SPI, respectively. 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 Figure 6. Spatial distribution and time-series patterns of r, nRMSE, and MASE from random forest model results for SDCI and SPI. The vivid red in the six maps indicates areas of relatively high errors (a-f). The green and purple in three time-series graphs indicate the time-series patterns of SDCI and SPI, respectively.

The Spatial Distribution of the Drought-Forecasting Model
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. 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 |forecasted-actual|. 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.

The Spatial Distribution of the Drought-Forecasting Model
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).
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 TCIGFS and PCIGFS was lower in that 8-day interval, which lowered forecasting skills. Figure 7. Spatial distribution of forecasted SDCI from 1 May to 2 June. The vivid red and blue present dry and wet conditions, respectively, in SDCI maps (i.e., reference SDCI (SDCI_o), SDCI from step 1 (SDCI_s1) and SDCI from step 2 (SDCI_s2). The light red and blue present dry and wet conditions, Figure 7. Spatial distribution of forecasted SDCI from 1 May to 2 June. The vivid red and blue present dry and wet conditions, respectively, in SDCI maps (i.e., reference SDCI (SDCI_o), SDCI from step 1 (SDCI_s1) and SDCI from step 2 (SDCI_s2). The light red and blue present dry and wet conditions, respectively, caused by temperature (TCI GFS ) and precipitation (PCI GFS ) from Global Forecasts System (GFS).
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].
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 21 Figure 8. Spatial distribution of forecasted SPI from 1 May to 2 June. The vivid red and blue present dry and wet conditions, respectively, in the SPI maps (i.e., reference SPI (SPI_o), SPI from step1 (SPI_s1) and SDCI from step2 (SPI_s2)). The light red and blue present dry and wet conditions, respectively, using precipitation (PCIGFS) from Global Forecasts System (GFS).
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 Figure 8. Spatial distribution of forecasted SPI from 1 May to 2 June. The vivid red and blue present dry and wet conditions, respectively, in the SPI maps (i.e., reference SPI (SPI_o), SPI from step1 (SPI_s1) and SDCI from step2 (SPI_s2)). The light red and blue present dry and wet conditions, respectively, using precipitation (PCI GFS ) from Global Forecasts System (GFS).

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.

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 • ).

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/12/21/3499/s1, Figure S1. Spatial distributions of r between both drought indices (SDCI and SPI) over 2018; Figure S2. Spatial distributions of r and nRMSE between satellite products and linear-fitted numerical model outputs; Figure S3 Drought conditions for SDCI on DOY 2018129 in Shandong province.
Author Contributions: S.P. carried out the manuscript writing and contributed to the data analysis and research design. J.I. supervised this study, contributed to the research design and manuscript writing, and served as the corresponding author. D.H. contributed to data processing and analysis. J.R. contributed to the discussion of the results and manuscript writing. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.