Applying machine learning for drought prediction in a perfect model framework using data from a large ensemble of climate simulations

. There is strong scientiﬁc and social interest to understand the factors leading to extreme events in order to improve the management of risks associated with hazards like droughts. In this study, artiﬁcial neural networks are applied to predict the occurrence of a drought in two contrasting European domains, Munich and Lisbon, with a lead time of one month. The approach takes into account a list of 28 atmospheric and soil variables as input parameters from a single-model initial condition large 5 ensemble (CRCM5-LE). The data was produced the context of the ClimEx project by Ouranos with the Canadian Regional Climate Model (CRCM5) driven by 50 members of the Canadian Earth System Model (CanESM2). Drought occurrence is deﬁned using the Standardized Precipitation Index. The best performing machine learning algorithms manage to obtain a correct classiﬁcation of drought or no drought for a lead time of one month for around 55-60 % of the events of each class for both domains. Explainable AI methods like SHapley Additive exPlanations (SHAP) are applied to gain a better understanding 10 of the trained algorithms. Variables like the North Atlantic Oscillation Index and air pressure one month before the event prove to be of high importance for the prediction. The study shows that seasonality has a high inﬂuence on goodness of drought prediction, especially for the Lisbon domain.


Introduction
Droughts remain to be one of the most dangerous hazards, having a serious and large-scale impact on environment, society and economy.Recent events like the Summer 2018 drought in huge parts of Central Europe led to severe forest fires and crop failures.The damage was estimated to several hundred millions euros solely in Germany (Federal Ministry of Food and Agriculture, 2018).Moreover the effect of global warming leads to major changes in the earth's climate system, having a direct influence on the frequency and severity of extreme events like droughts (Spinoni et al., 2016).An increase in frequency of drought occurrence is a major threat for current and future generations, and comprehensive knowledge on the phenomenon of drought is needed in order to take action early and to prevent humanitarian catastrophes.This goes in conjunction with drought prediction.Precise drought prediction would enable to mitigate the dangers connected to drought occurrences, such that e.g.stakeholders would be able to store the maximal possible amount of water in the endangered regions.This would help are projected to keep rising under RCP4.5 as well as RCP 8.5.(Spinoni et al., 2017).Observational data offers only a limited field for drought investigation as it can be seen from the following approximation.Systematical weather observations started in 1781 by Societas Meteorologica Palatina (Kington, 1980).In this study SPI1<-1 is used as a threshold for drought occurrence.
It corresponds to the 15% driest months (John Keyantash, 2018) and can be estimated by a total amount of 430 observed events until the year 2020 (Eq.1).
50 yr/member • 50 members • 12 months/yr • 15% = 4500 events (2) This is a difference of an order of magnitude.The more data is available the better the predictions that can be derived by a drought predicting machine learning model and the more can be learned about drought formation.According to von Trentini et al. (2020) precipitation in summer and winter derived from the European gridded data set (E-OBS) does fall to a high percentage into the range produced by CRCM5-LE for the historic period.Therefore, the CRCM5-LE proves applicable to this study and its larger amount of extreme events can be used as input to the machine learning algorithms.In this study a variety of ANNs are trained.Best performing models are investigated to using explainable AI methods to understand the results.
While no comparable study exists for the Munich domain, Santos et al. (2014) performed a drought prediction based on SPI6 for Portugal for the months April, May and June using the following input variables: sea surface temperatures (JFM), NAO (DJFM) and cumulative precipitation (NDJFM for SP I6 April , DJFM for SP I6 M ay , JFM for SP I6 June ).Best results were achieved for the prediction of SPI6 for April with a correlation coefficient of 0.98.SPI6 for May and June referred to a correlation coefficient of 0.78 and 0.77 respectively.
2 Data and Methods

Datasets
To investigate the predictability of droughts data from the single-model initial condition large ensemble (SMILE) consisting of 50 members, the Canadian Regional Climate Model 5 Large Ensemble (CRCM5-LE) is used.The data was produced within the scope of the ClimEx Project (Leduc et al. (2019), www.climex-project.org).The CRCM5-LE was generated by dynamical downscaling of the data provided by the 50-member initial condition Canadian Earth System Model 2 using the Canadian Regional Climate Model 5 (Martynov et al., 2013).The data has a resolution of 0.11 deg (12 km) and is produced for the  AMO is calculated by spatial averaging over the 0−60 • N, 0−80 • W area of the anomaly of sea surface temperature (E Trenberth, 2011).Additionally the 10-year running mean of AMO is calculated as an input variable, as it is widely used in various studies and was shown to be correlated with precipitation (Enfield et al., 2001).To limit the computation time a preselection of variables for the input data is performed.In order to eliminate redundant variables Pearsons R between the CRCM5 variables for the chosen domains is calculated.Pearsons R (ρ X,Y ) is a measure of linear correlation between two variables X and Y. ρ X,Y equals 1 if the correlation is total positive, 0 if there is no linear correlation and -1 if the correlation is total negative (Guyon and Elisseeff, 2003).For two samples x and y the Pearsons R is defined in the following way: (3) The bar refers to the average over the index i (Guyon and Elisseeff, 2003).ρ is calculated for all possible permutations of the 41 input variables.The ones correlating to a high degree are examined and a threshold of 0.95 is chosen.In Tab. 2 a list of sorted out variables and the corresponding values of Pearson R is given.The high correlation values can be explained by a physical relationship between the variables: e.g. the total evaporation (evspsbl) is almost the same as evaporation from land (evspsblland), as there are no relevant water bodies in the chosen domains.Out of the full list of 42 variables 14 are sorted out as being redundant.

Standardized Precipitation Index
The Standardized Precipitation Index (SPI) is a precipitation based index introduced by McKee et al. (1993).For the calculation of SPI a continuous monthly precipitation dataset is used.The index can be calculated on different timescales: typically, it is 1, 3, 6, 12 or 24 months.As a first step the precipitation values are accumulated for the needed timescale.The resulting dataset is fitted to a Gamma distribution for each month separately and then transformed to a normal distribution, such that the mean SPI is zero.The SPI value for a given precipitation is then the number of standard deviations from normal.Because of  the normalization SPI is especially useful to represent wetter and drier climates, as well as to account for differences among seasons.Here SPI1 is calculated for Lisbon and Munich each using the data from 1955-2005 from all members as reference.

Machine learning
This study investigates drought predictability applying the technique of supervised machine learning for this purpose.Machine learning is a promising tool for the analysis of complex and data-rich phenomena as droughts (McGovern et al., 2017).The python package Keras, a high-level neural network package, is used for the design of the machine learning models (Chollet et al., 2015), as it allows to design neural networks in an easy way by adding layers.Three crucial elements are needed to perform drought prediction by supervised machine learning: input data, a target variable to be predicted and a computation pipeline, which includes the machine learning algorithm.
The data from the years 1957 -1999 is used as training data, the years 2000-2005 are used for the testing purpose.Each of the time periods is available 50 times as we are dealing with an ensemble of 50 members.This results in 2150 model years for training and 250 years for testing.A small fraction of the training data is used for the validation of the machine learning algorithms.The target variable chosen for the prediction of droughts is SPI1.Two classes for the prediction are identified in the following way: SP I1 < −1 is defined as an event and is initialized with 1, SP I1 > −1 is initialized with 0 and corresponds to a non-drought event.The lead time of one month is chosen for the prediction.After the feature selection 28 variables originating directly from the CRCM5-LE dataset are used as input.In addition to those the teleconnection indices NAO, SCA, EA, EA/WR, AMO and AMO10 are used as input.
To predict e.g. a drought or non-drought in April of 1980, the data for twelve months before the event is used as input, this is NAO and other teleconnection and atmospheric variables for the period April 1989 -March 1980; for a prediction of an event in May 1980, May 1989-April 1980 is used as input.We perform a time series prediction with no limitation on special months or seasons to be inspected.
For this analysis we use a supervised machine learning algorithm, an Artificial Neural Network (ANN).ANNs are algorithms which design is inspired by the architecture of the human brain with its neurons (Russell and Norvig, 2009).; they both consists of connected nodes.A link between the node i and the node j serves to propagate the activation a i from i to j.To each connection a numeric weight w i,j is assigned.The output of the node is computed by: (4) (Russell and Norvig, 2009, p. 728).The activation function defines the output of the node.In order to have stable learners with confident predictions a function with a soft threshold is recommended (Russell and Norvig, 2009).In this study the following three activation functions are used: Sigmoid, Rectified Linear Unit (ReLU), Exponential Linear Unit (ELU).Sigmoid activation is especially useful for the output layer (Russell and Norvig, 2009), while ReLU and ELU both have the property of allowing very fast optimization (Maas, 2013) .
Sigmoid function, also called logistic function, is defined in the following way: (5) (Russell and Norvig, 2009).This function has an output between 0 and 1.This can be interpreted as a probability of belonging to the class 1.One of the main disadvantages of the sigmoid activation function is the vanishing gradient problem: at higher, almost saturated layers with values of 1 or -1, the gradients become nearly 0 resulting in a slow optimization convergence (Russell and Norvig, 2009, p. 726).
ReLU refers to Rectified Linear Unit and shows better performance when dealing with the vanishing gradient problem (Maas, 2013).ReLU is defined in the following way: ELU refers to the Exponential Linear Unit and was introduced by Clevert et al. (2016).Clevert et al. (2016) claim that in experiments the ELU activation led to faster learning and significantly better generalization performance than ReLU and sigmoid activation.The function is defined as: α controls the value to which an ELU saturates for negative inputs.Per default the value is set to 1 such that the function saturates at -1.
Two kinds of layers are used in this study: Dense and Dropout.Dense refers to a regular fully connected neural network layer.Dropout refers to a layer which is randomly setting a fraction of inputs to zero at each update.This technique is used to prevent overfitting and therefore improving the performance of the algorithm (Chollet et al., 2015).The first part of the study concentrates on the methodological search for the best performing algorithms.A pipeline to search for the best performing architecture, value for L2-Regularization and loss function is built up.
The model performance is evaluated using Accuracy and F1-score (Sasaki, 2007).The latter one is especially useful when training on datasets with an imbalanced class distribution, as it is in the case of our dataset.Accuracy is defined in the following way: F1-score is a harmonic measure between precision and recall.Precision is the amount of true positives with respect to the amount of positively classified data.Recall is the amount of true positives with respect to the total number of positives in the data.F1-score is defined in the following way: Due to the class imbalance within the dataset we require that the accuracy on each class is at least 50%.In that case given the distribution of the test dataset of 1803 non-drought events to 387 droughts for Lisbon and 1848 non-drought events to 352 drought events for Munich a marginal F1-score of 0.26 for Lisbon and 0.24 for Munich is given.
The second part of the study analyzes the best performing algorithms (one for Lisbon domain, one for Munich domain) by applying explainable AI methods.SHAP (SHapley Additive exPlanations) is a state of the art method for interpretation of machine learning models, which was inspired by game theory (Lundberg and Lee, 2017).It estimates for each input feature an average marginal contribution to the prediction of the result and therefore allows a comparison of the contributions among different features.In addition to that the difference in predictability among the seasons is calculated and compared to gain a better understanding on the influence of seasonal weather patterns.

Results
This study consists of two parts: the first parts deals with a systematical search for the best performing setup of the ANN model for the two domains of interest: Munich and Lisbon.A repeated training is conducted by varying the values of parameters like the architecture of the hidden layers, L2-Regularization and the loss function.In the second part of the analysis the best performing models for the two domains are analyzed using explainable AI methods.

Model training results
For the design of the ANN it is crucial to perform fine tuning of the model parameters to find the optimal setup.An architecture has to have enough layers and neurons to capture the complexity of the dataset (Goodfellow et al., 2016).In order to find In such way a network is searched which captures the given complexity of the dataset.This is reached with an algorithm consisting of at least five layers.Additionally two dropout layers, which are setting a specified number of nodes to zero in a random way, are introduced in order to fight overfitting.

L2-Regularization
L2-Regularization is a broadly applied method to prevent overfitting on the training data (Bishop, 2007).The main idea behind regularization is to add a penalty term to the loss function, which will punish the classifier for complexity and force some of the weights to zero (Russell and Norvig, 2009).In case of L2-Regularization the punishing term is proportional to L2-norm of the weight vector.The weight of the punishing term λ determines the relative importance of the regularization.
The results of the training with different values of λ for L2-Regularization are shown in Tab. 3. Training results are displayed in this particular case as the regularization is introduced to prevent overfitting.Generally the performance on the test dataset is more important and will be inspected in following experiments.If λ is set to zero the regularization term vanishes.Especially in those cases the overfitting is high.For Lisbon overall higher performance could be seen for values of λ around 0.01, 0.001 or 0.0001.Models that are trained on the Munich dataset perform better with the value of λ of 0.001.Since the performance of the model on the F1-score has a higher importance for an imbalanced dataset than the pure accuracy the value of 0.001 is chosen for the following ANN model training.3. Results of ANN training for different values for λ for L2-Regularization.λ of 0.001 (highlighted in grey) is chosen for both domains for following training, since the performance of the model on the F1-score has a higher importance for an imbalanced dataset than the pure accuracy.

Loss function
As a next step the influence of the different loss functions on the model performance is investigated.Loss function is a function to evaluate how well a specific algorithm manages to fit the training data (Janocha and Czarnecki, 2017).It is an important part of the optimization function which has a direct influence on the updating of the weights of the ANN (Russell and Norvig, 2009).In addition to overall accuracy and F1-metric, the accuracies on the non-drought and drought classes on the test dataset are displayed.The results are shown in Tab. 4. Binary cross-entropy, mean absolute error and hinge loss functions show the best performance for the Munich domain.In contrast to that for the Lisbon domain only the mean absolute error loss function has an accuracy of higher than 0.5.Also in the case of the Munich domain mean absolute error shows a higher performance on the F1-score.Therefore mean absolute error is used for further analysis.

Model architecture
Lastly the models are trained on both domains using different architectures.Table 5 is displaying the model training results on the test dataset.The column "architecture" refers to the number of neurons in each Dense (De) layer separated by the *-sign.
For Dropout (Dr) layers the fraction of weights which are randomly set to zero is given.The model architecture consists of overall seven layers.For example the architecture for the model in the first line of Lisbon it shows an Accuracy of at least 0.5 on both classes and a higher performance on the F1-score.the nominal values of the variable contributions.The nominal Shapely values are normed and recalculated to a percentage of contribution to the prediction, e.g. the NAO1 value explains roughly 2.3% of the prediction for the Lisbon domain.
We see that for both domains the contribution to the prediction is broadly distributed among the many input variables.
Between Lisbon and Munich Shapely values show a distinct difference in the nominal values of the feature contributions: values for Lisbon are about 6 times higher than those for Munich (e. g. the contribution of NAO1 for Munich is around 0.3% and for Lisbon around 1.9%).
For the Lisbon domain, the variables with a higher impact are sea level pressure (psl), surface pressure (ps) and NAO one month before the event.The first two variables are strongly autocorrelated for the Lisbon domain due to its location at the sea.The strong influence of ps/psl and NAO shows the influence of the atmospheric pressure system on drought formation in Lisbon.It is also striking that the influence of the local pressure seems to be higher than the influence of NAO.The next two variables for the Lisbon domain with the strongest contribution to the prediction are Northward Near-Surface Wind (vas) and Evaporation (evspsbl).The latter variable has a very direct influence on the formation of drought given that if evaporation is getting lower, also the probability of formation of rain clouds decreases (Sheffield and Wood, 2011).The contribution of vas to drought formation in Lisbon needs to be further studied.For the Munich domain the highest influence is found for NAO1, psl1 , EAWR5 and ps1.The results indicate that NAO is the most influential drought predictor for Munich.Additionally the contribution of EAWR5 and SCA5 on the Munich domain cannot be neglected as they are found within the top five predictors.
A further investigation of this relationship is of interest for the understanding of drought formation in Munich.

Seasonality
In order to evaluate the influence of seasonality on the prediction the performance of the model is calculated separately for the four seasons.Since the distribution between the drought and non-drought classes is different among the seasons (e.g.range of 17% to 19% of drought events for the Lisbon domain) a rescaling of the number of drought and non-drought events is performed to ensure comparability among the results.To compare the performance a precision recall plot is used (Saito and Rehmsmeier, 2015).Recall and precision are calculated for each of the four seasons (MAM, JJA, SON, DJF) and for the two half years (MAMJJA and SONDJF) using the estimated scaling factors.Results of the calculation are shown in Fig. 5.
The dotted line is marking the line under which the classifier shows no skill.The line is defined as a proportion of drought events against overall amount of events (Saito and Rehmsmeier, 2015).For the Lisbon domain it becomes evident that the model performance is very different across seasons: higher precision of around 0.23 can be found during the winter half year.
However for the spring season and summer half year the recall rises, while precision goes down.For the Munich classifier the results for the different seasons are closer together in terms of recall.It shows a worse performance for the winter months (DJF), while fall, spring and summer show an overall better model performance.This is an indication that for the Munich domain better drought predictability is possible in spring, fall and summer.
An additional analysis is conducted to calculate the Shapely values separately for the four season and the two domains in order to understand the influence of the different variables on the prediction.The results of the analysis can be seen in Fig. 6 and 7.The results for the Lisbon domain show that NAO1 is the strongest predictor in winter and spring season, while the contribution of pressure on drought predictability is higher in fall, followed by NAO1.On the contrary for the summer season NAO1 is not among the top 10 predictors, but other teleconnection indices like EAWR5, NAO7 and SCA7.Those

Discussion and conclusion
Drought is a multiscale phenomenon and its formation and evolution is different to every climatology and season.In this study, we i) explored the possibilities of using the data provided by CRCM5-LE to predict droughts using ANN and ii) applied explainable AI methods to gain a better understanding of the results.A drought event is defined as a SPI1 less than -1 at the given  where five layers are fully connected and two layers are Dropout layers, obtain best results.Best performing models obtain accuracies of 57% for the Lisbon domain and 55% for the Munich domain.
The precision of the prediction in both cases is rather moderate, as a high percentage of data is misclassified.For Lisbon, classifier precision remains at around 22 %.This means that one out of four predicted drought events is an actual drought.
For the Munich case, this ratio is even lower and amounts to 18 %.However, the models provide an important basis for the development of future drought predicting models and offer a fruitful ground for the investigation of influence of single input variables during different seasons on drought formation.
Compared to the study by Santos et al. (2014), which investigated drought predictability in Portugal, the weak prediction accuracies of our study are not surprising.In Santos et al. (2014) SPI6 for April, May and June is predicted, however precipitation amounts for the months until March were also given as input.As SPI6 is calculated using the sum of 6 months precipitation, the model is receiving over the half of the information it needs for the calculation of the value.As no similar studies exist for the Munich domain, no comparison can be performed.
The second half of the study concentrates on the analysis of the obtained algorithms using explainable AI methods.Among the strongest predictors for the domains are NAO, psl and ps one month before the event.This underlines the importance of the atmospheric system on the drought formation.For the model trained for the Lisbon domain the variables of Northward Near-Surface Wind (vas) and Evaporation (evspsbl) followed.For the Munich domain, EAWR and SCA five month before the This study indicates that seasonality is a crucial factor for drought predictions.Precision and recall of the prediction is getting lower in summer for the Lisbon domain and for winter for Munich domain.Moreover while for Munich domain the spread of precision and recall across the seasons is rather low, huge differences are found for Lisbon domain: the trained model obtained higher recall and lower precision for spring and higher precision and lower recall for fall when comparing to the baseline of all data.The results show that for the Lisbon domain NAO1 is the strongest predictor in winter and spring season, while the contribution of pressure on drought predictability is higher in fall, followed by the contribution of NAO1.For Munich domain NAO1 is found to have one of the highest contributions for spring, summer and fall, while it could not be found among the ten strongest predictors for winter.
Further investigations are of interest for scientific research on both objectives.In terms of drought prediction, further research is possible within the same setting.The field of AI is evolving rapidly, showing new algorithms, methods and frameworks, such that there is a high potential for finding better suited algorithms (Hao, 2019).One of the main limitations of this study remains that an application of the obtained framework on observation data is not possible, due to the fact that observational data lacks a multitude of variables which are used as input in this study e.g.Heat Fluxes, radiation, etc. However the results obtained by shapely value calculation are of high importance for the choice of variables for a development of a future model which potentially could be applied to observational data.Given the high Shapely importance of NAO for drought prediction, other large scale variables, such as atmospheric blocking, can be added to the input variables.Moreover, the application to new domains is of interest to investigate the regionality of drought prediction possibilities.Explainable AI methods offer an important approach to improve the current limitations of machine learning models; their application is of high importance in the field of physical geography since it enables providing a physical interpretation to statistical results.
years 1950-2099 for a European and an eastern North America domain.For the years 1950-2005 the historical greenhouse gas concentrations and aerosol emissions are being used.Starting from 2005, the model introduces the RCP8.5 (IPCC, 2013) forcing scenario.A total of 42 atmospheric variables is available in a temporal resolution of one to three hours.They are used on monthly basis as input to the machine learning algorithms.The list of variables is provided in Tab. 1.

Figure 2 .
Figure 2. First four leading eigenfunctions of the mean sea level pressure in CanESM2.Percentage of variance the mode explains is given on top of the figures.

Figure 3 .
Figure 3. Learning curve for two chosen fitting examples: algorithm complexity insufficient (top) and overfitting (bottom).
3. Dense layer with 1000 neurons 4. Dropout Layer randomly setting 50% of weights to zero 5. Dense layer with 500 neurons 6. Dense layer with 100 neurons 7. Dense layer with 5 neurons We require the accuracy on both classes individually to be higher than 0.5 and search for an F1-score as high as possible.In case of the Lisbon domain three trained models are satisfying the criterion of at least 50% accuracy on each class: the model in the first, in the fourth and in the last row.Best performance in terms of F1-score is obtained for the last model with the following architecture: 5000*0.5*4000*0.5*1000*500*100.For the Munich domain only the first and the fourth models are satisfying the criterion of at least 50% accuracy on each class.For further analyses the first model is chosen, as it shows the highest F1-score.The following model architecture is used for the Munich domain: 4000*0.5*1000*0.5*500*100*5.In the next step those models are analyzed using explainable AI methods.3.2Explainable AI methods for the analysis of best performing algorithms3.2.1 Shapely valuesFor the Munich and Lisbon domain Shapely values are calculated using the results of the best performing models on the test dataset.For the calculation each of the 12 month used as input to the predicting algorithm for each variable is considered individually, resulting in 28 atmospheric variables *12 + 6 teleconnection indices * 12 = 408 variables.The number behind the variable name refers to the number of months before the event (NAO1 -NAO value one month before the predicted event).The results are shown in Fig.4.Since the calculation of Shapely values is computationally expensive they are calculated 5 times on a subset of 500 data points.The error bars displayed in black on the plot indicate that the uncertainties are smaller than Performance of models for the Lisbon and Munich domains for different variations of architecture on the test dataset.Acc nd refers to the accuracy on the non-drought class and Acc d to the accuracy of the drought class.Sixth model architecture is chosen for following analysis for Lisbon and first for Munich (highlighted in grey), due to an Accuracy of at least 0.5 on both classes and a higher performance on the F1-score.

Figure 4 .
Figure 4. Mean Shapely values normalized to the contribution to the prediction for the top 15 variables with the highest importance for Lisbon (left) and Munich (right) on the test dataset.The number behind the variable name refers to the number of months before the event (NAO1 -NAO value one month before the predicted event).The results indicate that for the Lisbon domain psl1 and ps1 are the most influential drought predictors, for Munich this is NAO1.

Figure 5 .
Figure 5.The effect of seasonality on precision and recall for Lisbon (blue) and Munich (green).The results indicate that for the Munich and Lisbon domain better drought predictability is possible in spring, fall and summer.
site.The first half of the study deals with the systematic search for best performing models.For the Lisbon domain the model with L2-Regularization of 0.001, mean absolute error as loss function and the architecture 5000*0.5*4000*0.5*1000*500*100,where five layers are fully connected and two layers are Dropout layers, obtain the best results.For the Munich domain the model with L2-Regularization of 0.001, mean absolute error as loss function and the architecture 4000*0.5*1000*0.5*500*100*5,

Figure 6 .
Figure 6.Shapely values for Lisbon calculated separately for the four seasons and sorted by the maximum contribution in DJF (top left), MAM (top right), JJA (bottom left) and SON (bottom right) for the test dataset.evspsbl abbreviated as evp.

Figure 7 .
Figure 7. Shapely values for Munich calculated separately for the four seasons and sorted by the maximum contribution in DJF (top left), MAM (top right), JJA (bottom left) and SON (bottom right) for the test dataset.sfcWindmax abbreviated as sfcWm.

Table 2 .
List of sorted out variables

Table 4 .
Performance of the model for different loss functions on the test dataset.Acc nd refers to the accuracy on the non-drought class and Acc d to the accuracy of the drought class.Mean absolute error (highlighted in grey) is chosen for following analysis, since for Munich and