A Novel Fusion-Based Methodology for Drought Forecasting

Abstract

Accurate drought forecasting is necessary for effective agricultural and water resource management and for early risk warning. Various machine learning models have been developed for drought forecasting. This work developed and tested a fusion-based ensemble model, namely, the stacking (ST) model, that integrates extreme gradient boosting (XGBoost), random forecast (RF), and light gradient boosting machine (LightGBM) for drought forecasting. Additionally, the ST model employs the SHapley Additive exPlanations (SHAP) algorithm to interpret the relationship between variables and forecasting results. Multi-source data that encompass meteorological, vegetation, anthropogenic, landcover, climate teleconnection patterns, and topological characteristics were incorporated in the proposed ST model. The ST model forecasts the one-month lead standardized precipitation evapotranspiration index (SPEI) at a 12 month scale. The proposed ST model was applied and tested in the German federal states of Brandenburg and Berlin. The results show that the ST model outperformed the reference persistence model, XGBboost, RF, and LightGBM, achieving an average coefficient of determination (R²) value of 0.845 in each month in 2018. The spatiotemporal Moran’s I method indicates that the ST model captures non-stationarity in modeling the statistical association between predictors and the meteorological drought index and outperforms the other three models (i.e., XGBoost, RF, and LightGBM). Global sensitivity analysis indicates that the ST model is influenced by a combination of environmental variables, with the most sensitive being the preceding drought indices. The accuracy and versatility of the ST model indicate that this is a promising approach for forecasting drought and other environmental phenomena.

1. Introduction

Drought is a natural hazard and occurs in almost all regions of the world [1]. In Europe, severe droughts cause losses to agriculture, society, and ecosystems with profound impacts on crop production, the energy sector, and water supply [2,3,4]. In view of the drought threat, a suitable drought forecasting system is needed. Drought forecasting plays a key role in providing drought early warning to mitigate its impacts and to improve drought management [1].

Drought is a complex phenomenon and is among the least understood natural hazards due to its multiple causing mechanisms and contributing factors [5,6]. Drought phenomena can be categorized as meteorological, agricultural, hydrological, and socio-economic drought based on their duration and effects. Drought indices have been developed that quantify single or multiple droughts affecting variables in a standardized or threshold-based manner. The meteorological anomalies, land surface processes, and human activities interact and shape the development and evolution of drought [7]. These interactions are nonlinear and non-stationary, which complicates our understanding of drought estimation [8,9,10]. Machine learning (ML) has proven effective in dealing and solving many types of nonlinear problems. ML allows for complex mathematical calculations on mega data (e.g., large datasets such as long-time series satellite products including millions of grid points) [11,12]. The use of machine learning approaches with open-access satellite products facilitates drought forecasting from the regional to the global scale [8,13,14,15].

Multi-model fusion-based models are a promising ML method that can capture nonlinear processes. Alizadeh and Nikoo [8] argued that fusion-based ML models are more efficient than individual ML models in dealing with drought-related non-stationarities. Stacking (ST) is a hierarchical fusion framework integrating the strengths of individual estimation models [16], and it is a member of the family of ensemble learning methods (i.e., bagging, boosting, and stacking). Generally, the performance of an ensemble model exceeds that of a single learner. Compared with boosting and bagging, ST takes advantage of the characteristics of multiple base learners to achieve better predictive accuracy and stability than any of the individual models. The hierarchical fusion framework has been recently applied in several fields, such as PM 2.5 estimation [17], agriculture yield prediction [9], and landslide susceptibility mapping [18].

The ST model features a black-box structure, which does not reveal the underlying physical processes [19]. The black-box nature may complicate the quantification of model uncertainties [20]. On the other hand, a model-sensitive analysis is helpful in quantifying model uncertainties by evaluating the relative importance of input factors on model output. In practice, that means all uncertain variables, mathematical functions, and boundary conditions are allowed to vary to evaluate changes in model prediction [21,22,23].

Early drought detection is fundamental to proactive decision making and disaster preparedness. This work developed a fusion-based one-month lead meteorological drought forecast model that incorporates multi-source features from satellite images and data products. The standardized precipitation evapotranspiration index (SPEI) was applied to capture drought characteristics in the study area [24,25]. The developed ST model integrates three algorithms, namely, extreme gradient boosting (XGBoost), random forest (RF), and light gradient boosting machine (LightGBM) to forecast the SPEI at the time scale of 12 months (SPEI-12). The ST model receives climate teleconnection pattern indices, extreme climate indices, hydroclimatic conditions, vegetation coverage, and anthropogenic factors as forcing data. The sensitivity of the ST drought model is then analyzed using variance-based global sensitivity analysis (GSA). The ST model was applied and tested in the German federal states of Brandenburg and Berlin. This paper’s modeling method is useful in disaster prevention, such as meteorological drought early warning. It helps users understand the key determinants driving the meteorological drought forecast using an explainable ML algorithm.

Section 2 presents an overview of the study area and the used datasets. The ST model is described in Section 3, followed by a model sensitivity analysis in Section 4. The results are presented in Section 5. Future research directions are discussed in Section 6, followed by conclusions in Section 7.

2. Study Area and Methodology

2.1. The Study Area

This work’s methodology was applied to forecast meteorological drought severity in Berlin-Brandenburg, Germany (Figure 1). Berlin-Brandenburg is located in the North German lowlands with flat topography and consists mainly of agricultural land (49%), forestry land (35%), water bodies (3.5%), and urban areas (9%). Berlin itself consists of 18% forest, 4% agriculture, and 70% settlements and roads [20].

Berlin-Brandenburg is one of the warmest and driest regions in Germany, and thus it is particularly vulnerable to droughts. The long-term (1961–1990) average annual precipitation is less than 600 mm, and the average air temperature is near 8 °C to 10 °C [26]. In recent decades, the mean annual temperature in the region has increased by about 1 °C, exceeding the global temperature trend of 0.7 °C [27]. Previous studies indicate that evapotranspiration exceeds annual precipitation, and that the water balance in Berlin-Brandenburg is negative [26]. Regional climate simulations show a shift in precipitation patterns for the next decades with extremely dry summers and rising temperatures [28].

2.2. Datasets

Predictors include variables describing hydroclimatic, large-scale teleconnection patterns, as well as vegetation, terrain, anthropogenic, landcover, and historical drought indices. Several hydroclimatic datasets related to soil moisture, evapotranspiration, surface runoff, and terrestrial water storage were included in this study as predictors to measure the land surface response to drought conditions. Large-scale teleconnection indices were also included because drought is known to be related to large-scale atmospheric circulation and sea surface temperature conditions. Vegetation datasets were used to depict the growth conditions of vegetation in response to drought. Terrain variables provide location information to capture spatial variations in drought. Moreover, human activities and land use characteristics can impact drought occurrence and persistence through influences on water resource management, climate, and ecosystem dynamics. Finally, historical drought indices were applied to capture serial autocorrelation because drought is an inherently persistent phenomenon. These data provide a comprehensive interpretation of drought from different perspectives. Gridded data are used to reduce biases from irregularly distributed gauge stations and to account for the spatiotemporal variability of drought-related variables such as precipitation [13,29]. This study implemented gridded datasets from January 2001 to December 2018 for model training, validation, and testing to achieve available and consistent predictors. It should be noted that the calculation of the SPEI-N (N = 1, 2, …, 12) with different timescales covers the period 1989–2018, employing 30 years of data records to capture climate variability, long-term trends, and statistical significance in assessing drought conditions. All the selected datasets were of high quality and had a resolution of around 1 km to reduce bias caused by upscaling and downscaling [30,31,32,33]. Datasets that differed from the target resolution were interpolated to 1 km by bilinear interpolation. Moreover, the datasets applied in this study were intentionally restricted to those that are regularly updated so that our model can also be adapted to make predictions for future drought. The variables and metrics evaluated in this study are listed in

Table 1. Descriptions of whole datasets.

Category	Variables and Metrics	Abbreviation	Unit	Period	Explanation and Indication
Hydroclimate characteristics	Terrestrial water storage	TWS	cm	2001–2018
	V windspeed at 10 m	VW	m/s	2001–2018
	U windspeed at 10 m	UW	m/s	2001–2018
	Surface runoff	SR	m	2001–2018
	Surface latent heat flux	SLHF	J/m2	2001–2018
	Surface pressure	SP	Pa	2001–2018
	Total evaporation	EVP	m of water equivalent	2001–2018
	Humidity	RH	%	2001–2018
	Soil moisture	SM	%	2001–2018
	Mean area temperature	AT	°C	2001–2018
	Precipitation	P	mm	2001–2018
Extreme climate indices	Heavy precipitation days	R10	days/month	2001–2018	Count of days when precipitation ≥ 10 mm
	Very heavy precipitation days	R20	days/month	2001–2018	Count of days when precipitation ≥ 20 mm
	Consecutive dry days	CDD	days/month	2001–2018	Maximum consecutive days with precipitation ≤ 1 mm
	Consecutive wet days	CWD	days/month	2001–2018	Maximum consecutive wet days with precipitation > 1 mm
Large-scale teleconnection patterns	Pacific decadal oscillation	PDO	/	2001–2018
	Atlantic multidecadal oscillation	AMO	/	2001–2018
	North Atlantic oscillation	NAO	/	2001–2018
	South oscillation index	SOI	/	2001–2018
	Trans-Niño index	TNI	/	2001–2018
	Niño-based index	Niño 1 + 2	/	2001–2018
		Niño 3
		Niño 4
		Niño 3.4
Vegetation characteristics	Normalized difference vegetation index	NDVI	/	2001–2018	(NIR − R)/(NIR + R)
	Enhanced vegetation index	EVI	/	2001–2018	2.5 × (NIR − R)/(NIR + 6.0 × R − 7.5 × B + 1)
	Leaf area index	LAI	/	2001–2018
Terrain characteristic	Latitude	LAT	/	2001–2018
	Longitude	LON	/	2001–2018
	Slope	SLO	°	2001–2018
	Altitude	ALT	m	2001–2018
Anthropogenic characteristics	Night light index	NL	/	2001–2018
	Population	POP	/	2001–2018
Landcover characteristics	Shannon’s diversity index	SHDI	/	2001–2018	$-\sum _{i=1}^m(P_i\times \ln P_i$)
	Shannon’s evenness index	SHEI	/	2001–2018	$\frac{-\sum _{i=1}^m(P_i\times \ln P_i)}{\ln m}$
	Percentage of landscape	PLANDn	%	2001–2018	$\frac{\sum _{j=1}^n{a}_{ij}}{A}\times (100)$
Time characteristics	Season	S, A, W	/	2001–2018	One hot encoding
	Month-of-year	T	/	2001–2018
Historical drought index	Standardized precipitation evapotranspiration index	SPEI-Nt−1	/	2001–2018

Note: All dataset sources are listed in Table S1. NIR, R, and B represent the near-infrared, infrared, and blue laser bands, respectively. P_i denotes the proportion of the landscape occupied by patch class i. and m, a_ij, and A, denote the number of patch classes, the area (m²) of patch ij, and the total landscape area (m²), respectively. N is a timescale (N = 1, 2, 3…, 12).

.

2.2.1. Hydroclimate Datasets

Terrestrial water storage (TWS) observations were provided by the Gravity Recovery and Climate Experiment (GRACE) mission on a global scale. GRACE has been applied in drought investigations [13,34]. This study used the TWS data derived from the gridded GRACE Release-06 mascon product at the 25 km spatial resolution for the comprehensive monitoring of changes in the total volume of water stored in a region [35].

Meteorological drought is a weather-dependent phenomenon. Therefore, wind speed (v component of 10 m wind, and u component of 10 m wind), surface runoff, surface latent heat flux, surface pressure, and total evaporation from the fifth generation of the European Centre for the Medium-Range Weather Predicts (ECMWF) Atmospheric Reanalyzes of the Global Climate (ERA5 land reanalysis) dataset were also included as predictors. Furthermore, relative humidity, surface soil moisture, mean air temperature, and precipitation data were obtained from the climate data source of the German Weather Service (DWD).

2.2.2. Extreme Climate Dataset

Sun et al. [36] argued that some extreme climate indices have a close connection with droughts. Extreme indicators have been widely used to analyze extreme events in different regions of Europe [37,38]. Four climate indices were selected in this work, which are commonly used to describe drought in previous studies [39,40]. Extreme climate indices are herein calculated based on the DWD’s daily grid precipitation dataset, including heavy precipitation days (R10), very heavy precipitation days (R20), consecutive dry days (CDD), and consecutive wet days (CWD). The definition of each index is found in Karl et al. [41].

2.2.3. Large-Scale Teleconnection Patterns

Large-scale teleconnection patterns are known to be important driving factors for drought occurrence and are often used as predictors [1,42,43,44]. Linking drought occurrence to teleconnection patterns is useful for understanding how drought characteristics change in response to changes in atmospheric circulation [45]. Nine teleconnection pattern indices were selected based on previous studies in North Germany and Europe [40,46]. This work used them to represent climate anomalies (see Table 1). Detailed descriptions of those indexes can be found in Abdelkader and Yerdelen [47].

2.2.4. Vegetation and Terrain Datasets

Precipitation deficits reduce the photosynthetic capacity, which in turn changes the absorption of solar radiation in photosynthetically active wavelengths by plants [48]. Vegetation indexes (VI) are widely and successfully applied in meteorological drought predictions due to the change in greenness that reflects the influence of precipitation deficit on vegetation [49]. Three VIs were used in this study: the normalized difference vegetation index (NDVI), the enhanced vegetation index (EVI), and the leaf area index (LAI). The NDVI and EVI are derived from the MOD13A3 product, which provides monthly values at a 1 km spatial resolution. The EVI is a modified NDVI that improves sensitivity over high biomass regions and reduces the soil background effects and atmospheric influence [50]. The LAI is a critical variable in photosynthesis and precipitation interception processes [51]. It was acquired from the NOAA Climate Data Record of Advanced Very High-Resolution Radiometer Surface Reflectance (AVHRRSR).

Terrain features (latitude, longitude, slope, and altitude) are derived from the Shuttle Radar Topography Mission Digital Elevation Model.

2.2.5. Anthropogenic Datasets

Anthropogenic activities are likely to shape meteorological conditions [52]. Anthropogenically-induced climate change can exacerbate heatwaves and increase water demand, thereby aggravating droughts [53]. Human activity can be well assessed from anthropogenic light pollution [54]. The Night Light (NL) dataset and population data have been used in drought monitoring [55,56]. This study applied the NL dataset from the Harvard Dataverse to infer human activity [57]. This dataset is complemented by population distribution data (POP) from the LandScanTM High-Resolution Global Population Dataset.

2.2.6. Landcover Datasets

The land cover change influences the occurrence of precipitation and drought [58]. Land cover data involved in this study were acquired from the MCD12Q1 product with a 500 m resolution. The data comprised seventeen types of land cover (see Figure 1a). The landscape indices were used to quantify the spatial heterogeneity of land cover, geographical areas, or landscapes composed of interacting patches of vegetation or green space stands, or habitat types [59]. A total of 19 landscape indexes were derived in this work with the Fragstats V4.2.1 software representing the change characteristics of land cover [60], including Shannon’s diversity index, Shannon’s evenness index, and 17 landscape percentages.

2.2.7. Time Characteristics

Temporal information was included as a predictor to account for seasonal and trend effects. The four seasons were encoded through applying one-hot encoding [61], which represents the information in binary format, e.g., spring (March–May) as (0, 0, 0), summer (June–August) as (1, 0, 0), autumn (September–November) as (0, 1, 0), and winter (December–February) as (0, 0, 1). The month-of-year information was also used as a time characteristic.

2.2.8. The Standardized Precipitation Evapotranspiration Index (SPEI)

The SPEI quantifies the meteorological drought at multiple time scales [62]. It is a widely used standardized meteorological drought index that uses precipitation and potential evapotranspiration to determine the onset, duration, and magnitude of drought conditions with respect to normal conditions [63]. The use of the SPEI is particularly effective in drought forecasting whenever precipitation deficits and temperature anomalies have shaped past and recent drought periods, such as is the case in Germany [28,64].

The SPEI is herein calculated at the 1 km spatial scale, following the methodology described by Vicente-Serrano et al. [62]. Here, the calculation of the SPEI-N (N = 1, 2, …, 12) with different timescales is based on 30 years of data records (1989–2018) to ensure statistical robustness and reliability in quantifying the severity of drought. Its value ranges between −3 and +3 and can be categorized into several categories based on the drought severity degree [65] (see

Table 2. Drought categories based on SPEI values.

Categories	SPEI
Extremely dry	SPEI ≤ −2.0
Severely dry	−2.0 < SPEI ≤ −1.5
Moderately dry	−1.5 < SPEI ≤ −1.0
Mildly dry	−1.0 < SPEI < −0.5
Non-dry	SPEI ≥ −0.5

).

Historical drought indices were included as predictors. Autocorrelation, a commonly used measure of persistence, has been widely recognized in the literature [66,67]. The incorporation of SPEI-N_t−1 and SPEI-12_t−1 as predictors is grounded in their capability to capture the temporal patterns and autocorrelation inherent in historical drought indices. As the lag time increases, the correlation tends to decrease [68], reflecting the evolving dynamics of drought characteristics over time. This work took into account the autocorrelation of drought characteristics and introduced several timescale SPEI-N variables as predictors. The inclusion of SPEI-N_t−1 and SPEI-12_t−1 enabled us to explore distinct characteristics in detecting drought trends and periodicities associated with different timescales.

3. The Forecasting Model

The modeling framework consists of four tasks: (i) preprocessing, (ii) developing and implementing the ST model, (iii) developing and implementing the SHAP algorithm, and (iv) a sensitivity analysis (see Figure 2). The four tasks are described in detail in the following sections.

3.1. Preprocessing

Preprocessing selects the relevant variables from a pool of potential data sources. The variables without normalization are then used as input to the so-called base learners (a type of machine learning algorithm). The three best-performing ones on validation datasets are selected from the model pool, which then serves as the first layer of the ST model.

3.2. Machine Learning-Based Stacking Model

The ST model is an ensemble modeling approach that combines multiple machine learning algorithms [69], and it consists of two processing layers (Figure 2). The first layer contains several base learners. From this pool of potential base-learner algorithms, those that show good performance in the training phase are considered further. Specifically, during the training phase, which involves training and validation sets, the selection of base learners is based on the validation set. It is common to use no more than three base learners in ST models to achieve a good balance between accuracy and computational efficiency [16]. This study selected three base learners initially to construct the ST model’s first layer, and their various combinations were further verified. The three base learners, XGBoost, RF, and LightGBM, are popular models in drought and precipitation forecasting [70,71]. The XGBoost regression algorithm is a scalable ensemble-tree algorithm based on the improved gradient boosting decision tree (GBDT). The LightGBM regression algorithm achieves a high level of accuracy, and it processes faster than XGBoost. RF is an algorithm that forms part of the Bagging method in ensemble learning. An in-depth description of the XGBoost, LightGBM, and RF methods is presented by Chen and Guestrin, [72], Ke et al., [73], and Breiman, [74], respectively. The selected base learners build the first layer of the ST model. The training data are split into k folds, and k−1 folds are used to train the base learners. The first layer models are next applied to the remaining fold that was not used for the base model fitting. The test data are directly fed into the three well-trained base models. The predictions generated by three base models are then averaged, and the resulting average values serve as the input for a linear regression model (LR) to minimize the probability of overfitting during training (the second layer) [75]. Specifically, the second layer of the ST model is a meta-model that relies solely on the output results of the first layer.

3.3. Hyperparameter Optimization

Hyperparameters such as learning rate and batch size must be optimized in most ML algorithms [76]. Hyperparameter optimization is the process of selecting the best combination of hyperparameters for each independent base learner of the ST model to achieve their optimal performance on the training and validation datasets. The optimal hyperparameters are identified by combining the Opetuna and GridSearchCV approaches. Opetuna first ranks the significant hyperparameters based on the data characteristic [77]. Next, the optimal set of parameters is selected by the GridSearchCV algorithm [78]. In this study, the data were split 70% vs. 30% into a training and validation dataset (2001–2017) and a test dataset (2018) for hyperparameter optimization. The test dataset is only used to test the predictive power and generalization ability of the model. GridSearchCV conducts cross-validation on the training dataset to optimize hyperparameters, retraining each base learner model using the entire training period with the optimized parameters. Subsequently, here, we calculated the results for the validation set using the retrained base model. The optimal parameters are then used in all three selected base learners to construct the ST model.

3.4. Model Interpretation

The black-box character of many ML models hinders our capacity to understand the physical processes underlying the phenomena under investigation. Previous studies have used variable-importance diagrams to reveal the most and least important predictor variables [79]. However, this approach cannot explain the contribution of second-order effects (interaction) on a given model outcome. Dikshit and Pradhan [80] overcame this problem by implementing the SHAP method to improve the model’s interpretation with respect to spatial drought forecasting [80,81]. The SHAP method calculates the marginal contribution of each predictor to the forecast [82,83]. The average contribution of a predictor instance among all possible coalitions is given by the SHAP value. The SHAP value fairly distributes the contribution of predictors to the final prediction. The SHAP method provides multiple AI model explainers, including the Kernel Explainer, the Deep Explainer, the Tree Explainer, and the Gradient Explainer. A detailed explanation about different explainers and their corresponding plots with different uses is found in Zhang et al. [71] and Molnar [84]. The Tree Explainer was applied in this work to the outputs of the three base models. The Kernel Explainer was applied in the ST mode [85]. Three main types of interpretation can be conducted using the relevant explainer, namely, (i) the Summary plot; (ii) the Dependence plot; and (iii) the Force plot. The details of all plots can be found in Molnar [84]. The summary and force plots are two typical visualization figures focused on understanding the model’s decision process under drought conditions [80].

This study evaluated the ST model performance under different drought severities using the SHAP summary plot and force plot. The summary plot explains the relative importance of each variable on the model outcome. A SHAP force plot has three key characteristics (Figure 3). The first is the output value, which is the forecast value of a specific instance. The second is the base value: the average predicted value of the test dataset. The third is the color: variables marked in red color increase the magnitude of the model output, whereas those in blue color decrease the magnitude of the model output. The force plot shown in Figure 3 for illustrative purposes consists of a horizontal bar for four variables (SPEI-11_t−1, month-of-year (T), SPEI-12_t−1, and precipitation (P)). The force plot displayed in Figure 3 serves as a visual aid to facilitate understanding the marginal contribution of predictors to model output (i.e., forecasts) but does not explore the meaning of variables. The length of the bar represents the magnitude of the variables’ contributions to the model output. Longer bars in the SHAP force plot indicate the greater contribution of a variable to the model output. The SPEI-11_t−1 has the longest bar among the variables, followed by T. The impact of SPEI-12_t−1 and P is comparatively lower. The plot also shows the base value and the output value. The model output value (f(x) = −2.09, see Figure 3) is shown at the top of the plot, with the base value of −2.153 shown as a gray horizontal line. The horizontal axis in this visualization represents a small interval of the output and base values. The difference between the base value and the output value can be explained by the contributions of the individual variables. The red color indicates that the feature or variable increases the model output, while the blue color indicates that the feature decreases the model output. The SHAP algorithm is applied to improve the model’s transparency by determining the contribution of higher-order interactions to drought forecasting in Section 5.6.

3.5. Model Evaluation

The selection of base models from base pools and the accuracy of the SPEI prediction are assessed by the coefficient of determination (R²), the mean absolute deviation (MAE), and the root mean square error (RMSE). The best model and prediction results are comprehensive based on a higher R², lower RMSE, and MAE. The model is further compared with a persistence model (PTM) [86]. The PTM is often used as a reference for model evaluation, assuming the current SPEI at the time t will persist to the next month SPEI at the time t + 1 (see Equation (1)).

$${\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}}_{\mathrm{t}+1}={\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}}_{\mathrm{t}}$$

(1)

All SPEI values less than zero are classified as a drought in the assessment of drought occurrence [65]. The fraction of the number of grid points with SPEI < 0 in area i ($n_{0i}$) to the total number of grid cells N_t provides a measure to calculate the drought area DA (%).

$$DA=\frac{{\sum }_{i=1}^N{n}_{oi}}{N_t}$$

(2)

Multicollinearity [66] between the predictor variables is reduced by using the Akaike information criterion (AIC) and the Bayesian information criterion (BIC). The smaller the AIC and BIC values, the better the feature selection.

The prediction model’s capacity to deal with spatial and temporal non-stationarity issues is evaluated using an extended Moran’s I [9,87] which is a spatiotemporal autocorrelation index modified from the classical Moran’s I [88]. A positive value of Moran’s I represents a clustered pattern (i.e., positive correlation), a negative value represents a dispersed pattern (i.e., negative correlation), and a zero value represents a random pattern [89]. A value of the spatiotemporal Moran’s I close to 0 means that model predictions can capture spatial and temporal variations.

4. Model Sensitivity Analysis

This work applied Sobol’ method, which is a variance-based global sensitivity analysis (GSA) method to evaluate the different sources of uncertainty of predictors by decomposing the variance of the model output [23]. Contrary to more traditional SA methods, the Sobol’ method varies all input uncertainties simultaneously, which is required for nonadditive models [21,22,90].

The Sobol’ method commonly includes two measures for each parameter: the first-order index, which represents the contribution of a parameter to the response variance without interaction, and the total-order index, which represents the total contribution of a parameter to the response variance, including interaction effects. Saltelli’s extension of the Sobol’ sequence-based numerical procedures is applied in this study to generate input parameter sets for Sobol’ sensitivity analysis. More details about the GSA theory and Sobol’ method can be found in Sauter and Venema [22], and Saltelli et al. [23,91].

5. Results

Drought forecasting involves five tasks: (i) the choice of predictands (drought index); (ii) predictor selection; (iii) model selection; (iv) assessing and examining the model outcomes; and (v) model sensitivity analysis [2,92].

5.1. Predictand Selection

The SPEI-12 meteorological drought index was selected as the target variable of the ST model in this study. The selection of the target timescale of the SPEI-N index depends on the study area. Previous authors have reported that the SPEI-N of longer timescales (12 month time scales, SPEI-12) is more accurately predicted than the SPEI-N of shorter time scales [80,93].

5.2. Selection of Important Predictors

If two or more predictors are highly correlated only the predictor with the highest correlation with the predictand is considered. This is done to reduce data redundancy and model complexity. Previous drought index values were also included to exploit temporal coherence [1,94]. Furthermore, multicollinearity was tested with the AIC and BIC values (Table S2).

Table 3. Descriptions of selected variables.

Category	Selected Variables at t − 1	Abbreviation	Correlation with SPEI-12t
Hydroclimate characteristics	Surface runoff	SRt−1	0.253
	Soil moisture	SMt−1	0.365
	Precipitation	Pt−1	0.248
Extreme climate indices	Heavy precipitation days	R10t−1	0.213
	Very heavy precipitation days	R20t−1	0.158
Historical drought index	Standardized precipitation evapotranspiration index	SPEI-11t−1	0.227
		SPEI-12t−1	0.141
Time characteristics	Month-of-year	Tt−1	/

summarizes the eight selected predictors identified for the SPEI-12 from 68 predictor variables. The selection is conducted on training datasets. Eight time-series variables were selected as model predictors, which are three hydroclimatic factors, two extreme climate indices, two preceding drought indices, and the month-of-year. The selection standard of the preceding drought indices is shown in Table S3. It should be noted that SPEI-11_t−1 is not a variant of the SPEI-12_t−1 index. They describe the drought characteristics by considering different time windows. There were 7,857,530 × 8 valid training and validation data records (Years 2001–2017) available for model training once data preprocessing was implemented. Figure 4 displays the final model structure. It is a one-month lead drought forecasting model.

5.3. Model Performance and Comparison

XGBoost, RF, and LightGBM were chosen from the base learner pool to serve as the first layers of the ST model. The initial selection of these three potential base models was based on their default parameter settings without involving hyperparameter tuning. The selection standard was based on their higher average R² and smaller average RMSE (refer to Table S4) on the validation datasets. The forecast performance of the developed ST model was compared to three base models (XGBoost, RF, and LightGBM) and was also compared to the reference PTM model. The comparison was conducted for each month of the test year, 2018. The results are presented in

Table 4. Test accuracies of the drought-forecasting models with a one-month lead in year 2018, ST vs. three base models; ST vs. persistence model.

Month	XGBoost			RF			LightGBM			Stacking
	One Month Lead (t + 1)
	R2	MAE	RMSE	R2	MAE	RMSE	R2	MAE	RMSE	R2	MAE	RMSE
1	0.841	0.167	0.036	0.866	0.151	0.031	0.872	0.146	0.029	0.906	0.120	0.022
2	0.899	0.122	0.028	0.958	0.079	0.012	0.896	0.136	0.029	0.959	0.080	0.011
3	0.941	0.087	0.014	0.913	0.130	0.021	0.955	0.085	0.011	0.956	0.086	0.011
4	0.864	0.121	0.024	0.940	0.079	0.011	0.903	0.107	0.017	0.955	0.071	0.008
5	0.917	0.093	0.014	0.911	0.096	0.014	0.903	0.100	0.016	0.922	0.090	0.013
6	0.798	0.140	0.033	0.732	0.165	0.044	0.705	0.178	0.048	0.800	0.137	0.033
7	0.572	0.235	0.086	0.598	0.229	0.081	0.598	0.236	0.081	0.616	0.224	0.078
8	0.706	0.163	0.049	0.651	0.207	0.058	0.733	0.171	0.045	0.769	0.156	0.039
9	0.820	0.130	0.027	0.787	0.141	0.033	0.761	0.153	0.036	0.836	0.125	0.025
10	0.739	0.152	0.030	0.694	0.163	0.035	0.736	0.157	0.030	0.756	0.150	0.028
11	0.805	0.093	0.014	0.834	0.083	0.012	0.830	0.086	0.012	0.849	0.082	0.011
12	0.793	0.104	0.016	0.775	0.107	0.018	0.729	0.114	0.021	0.812	0.100	0.015
Ave	0.808	0.134	0.031	0.805	0.136	0.031	0.802	0.139	0.031	0.845	0.118	0.025
Std	0.097	0.040	0.019	0.113	0.048	0.021	0.102	0.042	0.019	0.098	0.042	0.019
Month	PTM						ST
	One Month Lead (t + 1)
	R2		MAE		RMSE		R2		MAE		RMSE
1	0		0.429		6.051		0.906		0.120		0.022
2	0.964		0.079		0.100		0.959		0.080		0.011
3	0.880		0.160		0.172		0.956		0.086		0.011
4	0.561		0.264		0.281		0.955		0.071		0.008
5	0.844		0.131		0.160		0.922		0.090		0.013
6	−2.423		0.713		0.752		0.800		0.137		0.033
7	−4.801		1.065		1.092		0.616		0.224		0.078
8	−2.186		0.708		0.746		0.769		0.156		0.039
9	0.471		0.267		0.287		0.836		0.125		0.025
10	−0.736		0.438		0.447		0.756		0.150		0.028
11	0.681		0.122		0.150		0.849		0.082		0.011
12	0.928		0.062		0.075		0.812		0.100		0.015
Ave	−0.401		0.370		0.859		0.845		0.118		0.025
Std	1.746		0.301		1.595		0.098		0.042		0.019

Note: In each test month, the highest R², lowest MAE, and RMSEvalues are listed in bold.

. The stacking model exhibited the best performance among all models, with the highest average R² (0.845) and the smallest average RMSE (0.025) (Table 4). In particular, compared to the PTM model, the ST model demonstrated a significant improvement in prediction accuracy. The ST model’s accuracy in summer was slightly lower compared to other seasons. This may have been due to rapid weather changes that usually occur in summer in the study region. A paired sample t-test was conducted between ST and each of the three base models and the PTM model in order to evaluate differences among the models’ results. The p-values of the t-test were below 0.05, except when compared with the PTM model (see

Table 5. Results from the paired sample t-test between the R² of ST, three base models, and the PTM model.

Model	T	p-Value
XGBoost vs. ST	−4.476	0.001
LightGBM vs. ST	−4.929	0.000
RF vs. ST	−4.250	0.001
PTM vs. ST	−2.484	0.030

), indicating that the accuracy improvement of ST was statistically significant. However, statistical metrics analyzed the variation between the forecasted and observed index values without considering the correct identification of drought and non-drought events. Therefore, the performance of the ST model, three base models, and PTM model with respect to the DA index was evaluated, with the results being listed in

Table 6. DA of the drought-forecasting models with a one-month lead in year 2018.

Month	XGBoost (%)	RF (%)	LightGBM (%)	Stacking (%)	PTM (%)	Observation (%)
1	0.013	0.016	0.016	0.008	0.503	0.752
2	0.142	0.047	0.013	0.078	0.752	2.068
3	1.595	0.836	1.510	0.689	2.068	0.661
4	0.041	0.021	0.010	0.026	0.661	1.235
5	0.555	0.503	0.150	0.849	1.235	1.003
6	19.392	16.813	10.798	20.624	1.003	25.164
7	97.011	97.075	95.495	97.271	25.164	96.211
8	99.910	99.559	99.812	99.737	89.211	99.790
9	99.997	100	100	100	99.791	99.993
10	100	100	99.995	99.995	99.993	99.995
11	99.995	99.995	99.995	99.995	99.995	100
12	100	100	100	100	100	99.995
Ave	51.554	51.239	50.650	51.606	43.365	52.239
Std	48.198	48.402	48.661	48.197	46.538	47.527

Note: In each test month, the highest DA values are listed in bold.

. The ST performed well compared with the reference PTM model, particularly in the summer season. The overall forecasting accuracy of ST was the best among all models, with an average DA (drought area) of 51.60%, and its forecasts were closest to the actual DA (52.24%).

This work resorted to a quantile–quantile plot (Q-Q plot) to assess the performance of the observed and forecasted SPEI-12 index corresponding to the ST model and the three base machine learning (ML) models across 12 test months, specifically focusing on extreme values (Figure 5). If the two distributions (reference and forecasted) being compared are similar, the blue points in the Q-Q plot will approximately lie on the identity line y = x (depicted in red). The best fitting was observed with the ST model. It is notable that the areas with extreme and severe drought were underestimated due to the limited number of extreme and severe drought samples in the training dataset. The underestimation was most pronounced for LightGBM among the three base models (Figure 5C).

5.4. The Spatial Pattern of the Predicted Drought Index

Figure 6 depicts the spatial patterns of the forecasted and observed drought indices. All models predicted the spatial structure of SPEI-12 well. However, a closer look revealed some shortcomings in all models. For example, the models tended to underestimate the intensity of severe drought events in the western and southwestern regions.

The boxplot of the differences between the forecasted and observed (referred to as prediction error) values is shown in Figure 7. The lower and upper edges of the box represent the 25th and 75th percentiles, and the horizontal line in the box denotes the median. Compared with three base learners, the median of the prediction errors for the ST model (indicated by the blue boxes) was closer to zero, suggesting reasonably accurate overall predictions. In addition, the length of the ST model’s box (interquartile range, IQR) was smaller in most test months, which indicated a consistent dispersion of prediction errors.

The spatiotemporal autocorrelation Moran’s I index, described in Section 3.5, was applied to analyze the spatiotemporal patterns of forecasting errors derived from the four models. In the ideal case, the models’ forecasting differences should be randomly distributed within the study area where the spatial non-stationarity is fully captured by the models [9,95]. This analysis aimed to investigate whether the spatial non-stationarity was fully captured in the models’ results. The average error was calculated in the 19 districts of Brandenburg and Berlin in each of the 12 test months to compare the spatiotemporal patterns and the degree of spatiotemporal aggregation. The spatial weighting matrix was generated using the inverse square distance between any two districts and is row standardized. The final spatiotemporal Moran’s I forecasting errors are listed in

Table 7. Spatiotemporal Moran’s I value of the forecasting errors calculated with ST, RF, XGBoost, and LightGBM.

Model	Spatiotemporal Moran’s I	p-Value
ST	0.0058	<0.001
RF	0.0069	<0.001
XGBoost	0.0074	<0.001
LightGBM	0.0080	<0.001

. The spatiotemporal Moran’s I values were close to 0, and the p-values were lower than 0.05, indicating that the prediction errors of four models in the 12 test months only had very slight spatiotemporal clusters. The degree of error aggregation of ST was the lowest (0.0058) among the four models.

5.5. Model Sensitivity

The ST model’s sensitivity pattern was analyzed based on the variance-based Sobol’ GSA approach described in Section 4, where the S_i indices measure the individual contributions of predictors to the output variance (ST), while S_t indices include all interaction effects. These two indices were calculated and averaged over the ST model in each testing month. Figure 8 shows Sobol’ sensitivity analysis results. It can be observed that the S_t indices were generally larger than the S_i indices, suggesting the presence of higher-order interactions. From these results, it follows that SPEI-11_t−1, SM, SPEI-12_t−1, and P may be designated as sensitive variables since the model output is mainly influenced by the uncertainty of these factors in all seasons.

Overall, the sensitivity analysis demonstrated that the performance of the ST model is most influenced by the SPEI-11_t−1 drought indices, soil moisture, SPEI-12_t−1, and precipitation as input variables.

5.6. Model Interpretation under Different Drought Severity Scenarios

The SHAP algorithm was applied to determine the contribution of higher-order interactions to drought forecasting. The meaning of the SHAP plot is described in Section 3.4. SHAP provides local interpretability by enabling us to identify the impact of each predictor on SPEI-12 prediction results for different drought conditions.

The ST model consists of a two-layer structure, where the prediction results from the first layer serve as input for the meta-model in the second layer. The SHAP algorithm for the ST model was applied to the meta-model. Figure 9 presents the results of force plots for five drought scenarios for the combined SHAP value over 12 test months: (a) overall performance, (b) extremely dry, (c) severely dry, (d) moderately dry, (e) mildly dry, and (f) non-drought. The detailed categorization of drought scenarios is listed in Table 1. Figure 9a shows that the preceding drought indices and soil moisture (SM) increased the model output, whereas the other variables slightly decreased the model output. Figure 9b–e depicts different drought conditions, and it was seen that the preceding SPEI-11_t−1 increased the model output, whereas month-of-year (T) decreased the model output. SM and precipitation (P) either increased or decreased the model output based on drought severity. The case of non-drought conditions (Figure 9f) also shows the importance of the preceding SPEI-11_t−1 in increasing the model output; T, surface runoff (SR), and P in decreasing the model output. Overall, the preceding SPEI-11_t−1 was the most influencing factor for all drought severity conditions. The final predictors of the ST model were derived from the outputs of the base learns in its first layer. Therefore, the overall predictor importance (summary plot) was based on the combined results of three base learners depicted in Figure 10A and the entire ST model (Figure 10B) for the purpose of validating the SHAP results. Although the ST model had a two-layer structure, the summary plot of the ST model was similar to the combination of the outputs from the three base learners in the first layer (see Figure 10A,B). The preceding drought indices were identified as the important features by both summary plots. In comparison with the base learners, the meta-model (second layer of the ST model) did not introduce a significant bias in the SHAP analysis.

6. Discussion

The ST model could be improved by stacking the drought classification and regression. That is, drought severity could be classified into different categories. Then, the regression model could be used for predicting sub-categories. Also, the time variable is usually included as a parallel factor, and the grid datasets themselves have potential location information. The ST model also relies on the autonomous learning ability of ML models to take into account the effect of spatial and temporal information. This approach has been demonstrated to improve the forecasting skill of ML models, but it neglects the fact that spatiotemporal information does not directly affect drought severity when applying ML models. Rather, spatiotemporal information affects drought severity through other factors attached to spatiotemporal location. Thus, a more effective representation of spatiotemporal information must be explored in drought modeling. This work applied the SHAP algorithm in the proposed two-layer structured ST model to explain the importance of variables in forecasting drought. This paper’s purpose was to validate whether the second layer, serving as the meta-model, introduces significant errors to the SHAP analysis due to the added model structure. The SHAP analysis results of the combined three-base model (first layer of the ST model) were compared with those of the entire ST model. The result demonstrated that the meta-model of the second layer did not introduce significant errors in the SHAP algorithm. It is important to exercise caution when using SHAP analysis with complex model structures (e.g., multi-layers). The stability of SHAP analysis with different model structures is a worthy topic for future research. Moreover, the structure of the ST model can be further optimized by incorporating wavelet decomposition, which is effective for capturing data interconnection in the frequency domain. This paper’s model performs short-term forecasts successfully; yet, extending the short-term (within 3 months) to medium- and long-term prediction is possible and could be explored by adding physical constraints. One potential change in future drought forecasting is to use the root zone soil moisture (RZSM) index to replace the SPEI index. The RZSM index can directly reflect regional soil moisture conditions and describe drought conditions. The wide acceptability of the ST model must be further tested with respect to a variety of drought conditions and drought indices at larger spatial scales.

In addition to the findings herein reported, we conducted a comparative analysis with the results presented by Mehr et al. [96] and with those reported by Mardian et al [97]. This paper’s model enables comprehensive monitoring of drought impacts in many areas across federal states lacking ground support at a 1 km spatial resolution in comparison with Mehr et al.’s [96] meteorological drought forecast model based on gauge stations. This paper’s ML model explains the contribution of variables in driving meteorological drought by breaking its black box structure. This paper’s ST model integrates the strengths of XGBoost, LightGBM, and RF models when compared with Mardian et al. ’s [97] individual XGBoost models. This study further explored the uncertainty errors of the SHAP algorithm introduced by different model structures.

7. Conclusions

Accurate drought forecasting is important for drought early warning and helps manage the risks associated with it. This study proposes a fusion-based ST model for one-month lead SPEI-12 drought index forecasting at 1 km spatial resolution. ML models are capable of solving complex non-linear problems but are weak in the interpretation of the underlying processes. However, the developed model can tackle the non-stationarity and underlying understanding by identifying the relationship between predictors and meteorological drought index. Specifically, 68 variables extracted from satellite products were evaluated in this study. Precipitation, soil moisture, surface runoff, heavy precipitation days, very heavy precipitation days, month-of-year time characteristics, and precedent drought indices were selected as predictor variables.

Climate teleconnection pattern indices were evaluated in this study, leading to the result that strong teleconnection is not observed between large-scale climate patterns and metrological drought in Berlin-Brandenburg. The performance of the developed ST model was compared with those of three base models (XGBoost, RF, and LightGBM) and a classic persistence model. The ST model performed best among the models, with an average R² value of 0.845 for the 12 test months. The spatiotemporal clustering analysis of prediction errors was further quantified based on the spatiotemporal Moran’s I. The error patterns analysis further showed that ST can effectively model non-stationarity in drought forecasting. The ST model can be applied to drought mapping with a relatively high spatial resolution, i.e., 1 km × 1 km.

The complex interconnection between predictors and drought is identified with the SHAP algorithm under five drought severity conditions, thus improving the understanding of the processes captured by the ST model. Moreover, model sensitivity analysis is important in practical applications (e.g., rain-fed farming systems) of the model; however, few studies include them in drought modeling. Furthermore, the model sensitivity is diagnosed with a variance-based GSA method. SHAP and GSA methods have demonstrated that the performance of the ST model is most strongly influenced by precipitation, soil moisture, and preceding drought indices, particularly SPEI-11_t−1 input variables.

This study’s methodology and results provide a valuable reference for drought early warning, assisting resource policymakers, as well as agricultural and water resource managers, in mitigating the impacts of droughts on agriculture, water supply, and ecosystems.