Multi-model drought predictions using temporally aggregated climate indicators

Journal of Hydrology 581 (2020) 124419

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Multi-model drought predictions using temporally aggregated climate indicators

T

⁎

Md. Mamunur Rashida,b, , Ashish Sharmaa, Fiona Johnsona a b

School of Civil and Environmental Engineering, University of New South Wales, Sydney, Australia Department of Civil, Environmental, and Construction Engineering, University of Central Florida, Orlando, FL, USA

A R T I C LE I N FO

A B S T R A C T

Keywords: Standardized precipitation index Aggregated climate indices Model combination Uncertainty

Long term skillful prediction of prolonged hydrologic anomalies is essential for proper planning to reduce the societal risk of extreme hydrological anomalies such as drought. Climate indices estimated from sea surface temperature anomalies (SSTA) of the Pacific and Indian Oceans are often used to predict monthly and seasonal rainfall in Australia and many other places around the world. This study investigates the merit of distorting the time aggregation of such indices before casting them in a predictive model. Aggregated climate indices are used to predict sustained drought and wet anomalies characterised here using a drought index (i.e. Standardize Precipitation Index, SPI) as response and the Australia as the study region of interest. The aim is to enhance the strength of relationships of drought index and climate indices (predictors) by tuning the frequency of climate predictors using an aggregation technique. Result shows that aggregated climate indices provide significant improvement in prediction of SPI over raw climate indices across Australia. As strong spatial variations in optimum aggregation window lengths are evident across Australia suggesting multiple candidate predictive models with similar accuracies, a model combination approach is also adopted. Model combination is found useful in reducing structural uncertainty and further improving the prediction efficiency. Given that the improved predictive accuracy for SSTA the current generation of climate models exhibit, the methodology developed in this study has significant implications for skillful prediction and projection of long term droughts.

1. Introduction Sustainable management of water resource systems requires skillful prediction of extended precipitation anomalies for water resource planning and allocation. Some of the important drivers of extended precipitation anomalies in Australia are related to key modes of the Pacific and Indian Oceans (Khan et al., 2015; Kirono et al., 2010; Risbey et al., 2009; Rogers and Beringer, 2017; Tozer et al., 2017; Westra and Sharma, 2010). The influence of these climate indices is strong enough to also affect long term drought (Dijk et al., 2013; Ummenhofer et al., 2009), with the climate teleconnections particularly strong at interannual and interdecadal frequencies (Özger et al., 2009). The influence of El Niño Southern Oscillation (ENSO), Pacific Decadal Oscillation (PDO), Indian Ocean Dipole (IOD) and the subtropical ridge (STR) on Australian rainfall variability are well documented (Cai and Cowan, 2013; Kirono et al., 2010; Murphy and Timbal, 2008). In Australia, typically La Niña and El Niño are associated with drier and wetter than the normal conditions respectively. For southern Australia, IOD has greater influence than ENSO during winter and

⁎

spring rainfall as a result of the strong covariation between ENSO and the IOD (Cai et al., 2011). Interannual variability of droughts over large part of Australia is linked to the co-occurrence of an El Niño and a positive IOD event whereas lower frequency variability such as multiyear drought have IOD as the dominating influence. (Risbey et al., 2009; Ummenhofer et al., 2009; Ummenhofer et al., 2011). This indicates that SSTA based climate indices have potential in predicting droughts in Australia. Among different approaches used for prediction of sustained droughts, models with large- scale climate indices is most common as their relationships with sustained droughts are reasonable (Awan and Bae, 2016; Ren et al., 2016; Vicente-Serrano et al., 2011). However, variability and plausible teleconnections between droughts and largescale climate indices are different for different temporal scales (Özger et al., 2009). For example, occurrence of long term droughts is dominated by the multi-year to decadal oscillations of climate indices (Dijk et al., 2013). Therefore, identify most useful time scales of climate indices by optimizing the relationships with droughts is crucial to improve drought prediction. Out of different modelling techniques used

Corresponding author at: Department of Civil, Environmental, and Construction Engineering, University of Central Florida, Orlando, FL, USA E-mail address: [email protected] (Md. M. Rashid).

https://doi.org/10.1016/j.jhydrol.2019.124419 Received 12 October 2018; Received in revised form 26 November 2019; Accepted 27 November 2019 Available online 02 December 2019 0022-1694/ © 2019 Elsevier B.V. All rights reserved.

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

with different predictor sets. The advantage of this approach is that an ensemble of the multiple models can be used which has been shown to provide better predictions on average than a single individual model. A second advantage of a multi-model ensemble is that it can improve our understanding of the uncertainty of the predictions of sustained drought and wet anomalies. It also circumvents the need to select a ‘best’ model from the set of alternate models (Hagedorn et al., 2005; Sharma and Chowdhury, 2011). In this study a multi-model combination approach is used to make predictions of long term drought and wet conditions obtained from alternate models.

for drought prediction, non-parametric data driven models become popular in the recent time (Anshuka et al., 2019; Fung et al., 2019). Many drought indices have been developed to represent sustained drought and wet anomalies (Mishra and Singh, 2010). Among those, Standardized Precipitation Index (SPI) is the most widely used index due to its simplicity in formulation and application. SPI can be compared across different locations due to its transformation to a standard normal distribution. SPI can be used to estimate drought and wet anomalies at different time scales. While SPI estimated at shorter time scales (e.g. 3 and 6-month) represent short-term droughts (i.e. meteorological droughts), SPI estimated at longer time scales (e.g. 12 and 36-month) is capable to capture long-term droughts (i.e. hydrological droughts) (Li et al., 2016). In Australia and many other places of the world, the droughts of most concern for water resource management occur over year to multi-year periods and are often quantified with SPI at longer time scales (Dijk et al., 2013; Sheffield and Wood, 2007). We hypothesise here that a predictive model for drought does not necessarily need to use SSTA predictors at the same time scale as the drought response of interest. Or, one may need a modulation of the SSTA predictors to a different temporal scale to assist in the formulation of the best predictive model. Such a modulation is proposed here by means of aggregating the predictor variables over variable length time scales of interest. While such aggregated predictors have been used before in formulating downscaling and other predictive models (Mehrotra and Sharma, 2006), their use in a seasonal forecasting context has not been attempted. To understand the relationships at these longer time scales, methods that emphasise or filter the low frequency components of the climate indices could be therefore useful to predict long term drought and wet conditions. Temporal aggregation is a simple way to filter the different frequency components of a time series. Wavelet transforms also often used to filter the signals at different frequency bands to use in model development (He et al., 2014; Rashid et al., 2015a; Rashid et al., 2015c; Rashid et al., 2018). Temporal aggregation of climate indices has been used to forecast seasonal rainfall (Kirono et al., 2010), standardized streamflow (Seibert et al., 2017) and to identify teleconnections with rainfall anomalies (Dijk et al., 2013; Shi et al., 2016). However the application of temporally aggregated climate indices in a modelling framework to predict long-term sustained drought and wet anomalies has not been explored. Here we investigate using filtered relationships of climate indicators to predict a drought index. The aim is to maximize the strength of relationships between the drought index and climate indices by smoothing climate indices time series using temporal aggregation. In this study, we develop a modelling framework where low frequency components of climate indices obtained through temporal aggregation were used as predictor in the K-nearest-neighbor (KNN) to predict sustained drought and wet anomalies. Directly projecting future drought and wet anomalies from GCM simulated precipitation is questionable due to the limited ability of these models to correctly simulate the low frequency variability of precipitation. While bias correction through nesting at different time scales can go some way towards correcting low frequency biases (Johnson and Sharma, 2011; Mehrotra and Sharma, 2016), there are questions about the assumed stationarity of these biases (Nahar et al., 2017). In contrast, developing a framework to project SPI from sea surface temperature (SST) is promising because of the skill of GCMs in predicting sea surface temperature (SST), especially over the tropical Indian and Pacific Ocean (Song et al., 2014; Wang et al., 2015). In addition, SST predictions may be further improved using multi modal combination (Kaiser Khan et al., 2014). Using predicted SSTs to derive the SPI hence can offer advantages when compared to the use of GCM simulated rainfall directly. Thus, this study provides a new approach that may improve future drought projections. Because the teleconnections of Australian rainfall and drought with the different climate indices varies strongly over space and time (Rashid and Beecham, 2019a; Rashid and Beecham, 2019b; Risbey et al., 2009; Tozer et al., 2017), there may be value in developing alternate models

2. Data Observed monthly gridded precipitation data from the Australian Water Availability Project (AWAP) have been used to estimate the drought response variable used in this study. The rainfall data was regridded to a spatial resolution of 0.5° using area averaging. Grid cells with zero rainfall for more than thirty percent of the data length were not considered for the analysis. These grid cells are located in central and western desert of Australia and coincide with the area identified by Jones et al. (2009) where gridded precipitation is unreliable due to the sparse gauge network. While there are several climate indices which have significant influence on Australian rainfall (Dijk et al., 2013; Tozer et al., 2017), this study uses Niño3.4 and DMI to predict SPI to simplify the presentation. Niño3.4 represents the sea surface temperature anomalies (SSTA) over the Pacific Ocean and is the leading mode of the El Niño Southern Oscillation (ENSO). DMI represents the SSTA over the Indian Ocean (Saji et al., 1999). When combined, ENSO and DMI are capable of explaining the variability of sustained rainfall anomalies across Australia (Cai et al., 2011; Fierro and Leslie, 2013; Ummenhofer et al., 2011). Niño3.4 and DMI are derived from monthly sea surface temperature values of Hadley Centre Global Ice and Sea Surface Temperature (HadISST) datasets (Rayner et al., 2003) for the period 1910 – 2010 (Table 1). 3. Methods 3.1. Standardized precipitation index (SPI) The most widely used drought index is the Standardized Precipitation Index (SPI) which can be estimated from precipitation data alone (McKee et al., 1995). The advantages the SPI are that it can be estimated for any time scale of interest and it is comparable across different locations because of the transformation to a standard Normal distribution. SPI is classically estimated by fitting a Gamma distribution to the aggregated rainfall at different time scales and then transforming this to a Normal distribution (Guttman, 1999). A range of other distributions have been recommended (Angelidis et al., 2012; Naresh Kumar et al., 2009) due to issues with fitting the Gamma distribution in areas with extremely low or high rainfall. Non-parametric distributions have also been proposed for SPI (Farahmand and AghaKouchak, 2015) to overcome the fitting issues. The parametric SPI method is sensitive to the choice of probability distribution and its parameters particularly when applied over large areas with a wide range of different climatic Table 1 SSTA climate indices considered as predictors of Australian drought in the study.

2

Climate index

Description

Region

Niño3.4 DMI

Average SSTA over 5°N to 5°S and 170°W to 120°W Indian Ocean West Pole Index (WPI) – Indian Ocean East Pole Index (EPI), where WPI and EPI is defined as: WPI: average SSTA over 10°N to 10°S and 50°E to 70°E EPI: average SSTA over 0°N to 10°S and 90°E to 1100E

Pacific Indian

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

variability compare to the raw climate indices, the reference model with raw climate indices as predictor is appropriate. The multi-model ensemble can be constructed using simple averaging of the ensemble members (Hagedorn et al., 2005; Rashid et al., 2015b) or by optimizing the weights to maximize predictive skill (Rajagopalan et al., 2002; Robertson et al., 2004) or estimating the weights conditioned to dominant predictors (Chowdhury and Sharma, 2011; Devineni and Sankarasubramanian, 2010). Weights can be static or temporally dynamic (Chowdhury and Sharma, 2009; Kaiser Khan et al., 2014). Here the multi-model combination is obtained as the weighted sum of each individual model prediction. The optimum weights for each model are estimated from the error variance-covariance matrix (Kaiser Khan et al., 2014).

conditions, the non-parametric SPI approach minimizes this issue. Because of the range of climate conditions across Australia and variety of precipitation distributions, the SPI is estimated here using the nonparametric method proposed by Farahmand and AghaKouchak (2015). 3.2. Aggregated climate indices The aim of this study is to assess whether the SPI can be better predicted using a multi- model ensemble with temporally aggregated climate indices as predictors. To isolate the most important signals present in these climate indices (CI), we have considered eighty temporal windows to aggregate the indices using a simple backwards moving average. For example, considering an aggregation window of L months (AWL), the aggregated climate index time series (ACIL) is estimated as:

3.4. Model implementation details

t−L

ACItL =

∑ CI t=t

Using the framework presented in the previous section and Fig. 1, predictions of the 12-month and 36-month SPI are made at each grid location across Australia for the period from 1910 to 2010 using the aggregated time series of Niño3.4 and DMI as the predictors. For each location and SPI time scale (12-month or 36-month), the following steps are required to select the significant predictors and formulate the four models used:

(1)

The ACI are the potential predictors for SPI and with L aggregation windows and two climate indices, the model can be formulated using a total of 2L predictor time series. The process of selecting a subset of predictors to use in the multi-model ensemble is discussed in the next section.

1. The order-1 predictors are selected from aggregated climate indices using PIC. A partial information analyses is undertaken using observed SPI as the response and all potential aggregated climate indices as predictors. The predictor with maximum PIC is chosen as the order-1 predictor for the first model (A1). A number of additional predictors [A1c] representing candidate predictors for the other three models considered, are also identified. 2. The order-1 predictor for the second model is then selected by using A1 as the response and remaining predictors in [A1c] using partial information analysis. The predictor with the lowest PIC is now chosen to be the order-1 predictor for the second model (A2). This ensures that A2 has the minimum correlation with A1, while still being strongly related to the response. It is thus assumed that using this additional predictor to form a new model will result in a predictive model that is structurally different to the model formulated earlier. 3. To select the order-1 predictor for the third model, A1 and A2 are excluded from the predictor set and a new predictor (A3) chosen from [A1c]. The partial information analysis now uses A1 as the response and A2 as a pre-existing predictor and the predictor with the lowest PIC from the remaining aggregated climate indices is then chosen as the third order-1 predictor (A3). This process is then repeated to ascertain the fourth order-1 predictor to form the predictors A1 to A4.

3.3. Predictor selection and model formulation framework One of our objectives is to develop a range of models to use in the multi-model ensemble. To maximise the effectiveness of the ensemble, it is desirable for predictors of each of the models to be as different as possible from each other (Sharma and Chowdhury, 2011). Predictors of interest are identified from the pool of predictors variables (i.e. the time series of aggregated climate indices) using a nonparametric approach based on Partial Information Correlation (PIC) (Sharma and Mehrotra, 2014). PIC is a measure of conditional dependence and similar to partial correlation. In general the PIC framework can be used to select predictors for a model by comparing the information in the predictors to the response, and is flexible enough to also consider additional preselected predictors that are presumed important from a physical perspective. Here PIC has been used to maximise the additional information provided by each predictor, equivalent to choosing new predictors with the smallest partial information correlation with a predictor that has previously been selected. A comprehensive flowchart of the proposed method for multimodal prediction of SPI is presented in Fig. 1 and further details are presented below. Let the general form of the predictive model be:

Yt ≡ Yt |[Ati , Bti]

(2)

where Yt is the response variable i.e. SPI at a specified time scale and Ati and Bti are the two predictors for the model, named order-1 and order-2 predictors respectively, both of which are aggregated climate indices for different aggregation periods. For clarity, we have restricted our study to four predictive models. Thus, At1 , At2 , At3 and At4 are the four order-1 predictors and Bt1, Bt2, Bt3 and Bt4 are the four order-2 predictors. Each of the four models is constructed using a K-nearest-neighbor (KNN) regression approach based on Partial Weights (PW) to predict SPI (Sharma and Mehrotra (2014). We have used the Nonparametric Prediction (NPRED) R Package for the analysis (Sharma et al., 2016). Additionally, a reference model is developed following the same method where raw climate indices are used as predictor instead of aggregated climate indices. The procedure of model fitting was kept the same to enable a fair comparison between two models. For example at any grid location, we have considered SPI and non-aggregated raw Niño3.4 and DMI time series to estimate PW and then predict SPI using regression based KNN model based on estimated PW. As our hypothesis is that the aggregated climate indices better represent the SPI

While we have restricted our study to four models, the number of models could be expanded by repeating step 1 to 4 but considering the remaining potential aggregated climate indices as predictors. Once the order-1 predictors have been selected for each model, the order-2 predictors are then added. This proceeds in the same way as step 3 above. For each model, the partial information analysis now uses SPI as the response, with remaining candidate predictors as the potential predictors excluding the four order-1 predictors already selected. For each model, the order-1 predictor associated is the pre-existing predictor for PIC. The predictor with the highest PIC is now chosen as the order-2 predictor for that model. Readers are referred to Sharma and Chowdhury (2011) for additional details about this approach. To construct the multi-model average, a model combination weight is estimated. The approach of Kaiser Khan et al. (2014) and Khan et al. (2017) has been used. This involves calculating the residuals for each of the four models and estimating their error covariance matrix. With four models (m = 4), the size of the error covariance matrix is 4 × 4 with 3

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

Fig. 1. Flowchart for multi-model prediction of SPI.

3.5. Assessment of model skill

the diagonals representing the error variance of each model and the offdiagonals the covariance of the errors between the models. The weights for each model are calculated using equation (3) and are static in time. −1

(

W = V' ∑

e

−1

)

−1

∑e

V

The performance of the proposed model is assessed using a cross validation approach where available 100 years of data record is divided into ten subsets consisting of 10 years long non-overlapping series. Each subset is used as a validation and the remaining (90-year period) used for calibration. A reference model is also formulated that does not consider the aggregated predictor variables used in the modelling framework described earlier. The mean square error (MSE) and Nash Sutcliffe efficiency (NSE) are used as a measure of prediction skill for all four individual models, the multi-model ensemble and the reference model. The improvement in the multi-model prediction is estimated in term of the relative differences in MSE and NSE between the multimodel and reference model.

(3)

where V is a (mX 1) dimensional column vector of ones and ∑e is the covariance matrix of prediction errors. To ensure that the multi-model prediction is unbiased, we add additional constraints that the individual model weights add up to 1 and all weights are greater than zero. The multi-model prediction (MM) at each time step is then the weighted sum of the predictions from the four models at the same time step as per Equation (4).

tMM = Y

4

∑ Wi Yt i=1

i

(4)

4. Results

t = [Y t1, Y t2, Y t3 and Y t4]' are the SPI time series predicted from where Y four models where t varies from 1 to n, the total number data points in the SPI series and W = [W1, W2, W3 and W4] are the optimum weights for each model.

4.1. Importance of temporal aggregation Our hypothesis is that by temporally-aggregating climate indices, 4

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

Fig. 2. Mean squared error (MSE) for different regression models at a typical grid location for 12 month (12 M) and 36 month (36 M) SPI. Blue line represents the MSEs obtained from models developed with aggregated climate indices for different aggregation windows and red dotted line represents the model with climate indices without aggregation (i.e. reference model). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

increases from 30 to 40 months over the north east coast, south east coast and Tasmania up to 50 to 70 months in the remainder of the country. The regional variability evident in the optimum windows suggests that not one single model to predict SPI will be successful.

their strength as potential predictors of sustained drought and wet anomalies can be improved. As a proof of this concept, simple regression models are constructed to predict SPI firstly using each climate index as a predictor. The improvements in this base model by using an aggregated climate index are then measured, using the mean square error (MSE). For the purpose of illustrating the worth of using timeaggregated predictor variables, an aggregation period (L) of 80 is used. This implies that predictors are selected from a set of 80 possible variables for both order 1 and 2 predictors for each individual model. Fig. 2 shows the MSEs obtained from the reference model and models with aggregated climate indices at different aggregation windows for 12-month and 36-month SPI at a randomly selected single location in validation. The models with aggregated climate indices provide consistently lower MSEs for each selected climate indices compared to the reference models. The blue curve in Fig. 2 demonstrates that the temporal aggregation of climate indices improves the model efficiency (i.e. decreases MSEs) up to an optimal window length when the MSE is the smallest. When aggregations longer than the optimum are used, the model error increases and in some cases becomes larger than the base model. The optimum window length varies with the climate drivers and time scale of the SPI response. There is strong regional variability in the optimal window across Australia (Fig. 3). Over eastern Australian there is a north–south gradient with longer optimum lengths in the north for Niño3.4. The length of the optimum aggregation windows over eastern Australia varies from 5 to 20 months and 30 to 50 months for the 12-month and 36-month SPI respectively. In Western Australia an east-west gradient to longer optimum windows is observed for Niño3.4 for both the 12-month and 36-month SPI. The relationship between DMI and SPI is strengthened when comparatively long windows are used (60 to 80 months) for both 12-month and 36-month SPI. Areas with these multi-year to almost decadal optimum windows include the far inland of New South Wales, South Australia and Northern Australia. Other parts of the country (e.g. north east coast, southern Australia and Tasmania) are dominated by subannual to annual variability of DMI when predicting the 12-month SPI. For longer droughts, the length of the optimum aggregation windows

4.2. Prediction of sustained drought and wet anomalies with aggregated climate indices Using the temporally aggregated climate indices as predictors, we developed multiple models using the framework described in 2.4. These models are compared to a reference model which uses unsmoothed climate indices as predictors. To illustrate the method in-depth, the results from four grid cells are first presented followed by summary results across Australia. The prediction errors (MSE) and model efficiency (NSE) for SPI from each of the four constituent models (M1, M2, M3 and M4), the multimodel ensemble (MM) and the reference model at the four selected grid locations are listed in Table 2. The order-1 and order-2 predictors for each constituent model are also listed. The value of temporally aggregating climate indices is clear, with all four constituent models outperforming the reference model at each location for both the 12month and 36-month SPI. Scatter plots of SPI predicted by the multimodel ensemble and reference model against the observed SPI along with their corresponding correlation coefficient (as shown in Figs. S1 and S2 of the supporting document), further demonstrate the superiority of using temporally aggregated climate indices to predict 12month and 36-month SPI. It should be noted that the results presented here and elsewhere in the document all represent cross-validation performances, and hence improvements in one versus the other are statistically significant. The multi-model approach allows the uncertainty in the SPI predictions originating from the choice of predictor variable to be explored. At each grid location, each of the four models performs quite differently illustrating the influence of different predictors and different aggregation windows on the SPI. This also illustrates the uncertainty in the predictive model being formulated, a scenario that is known to be benefitted through model combination. 5

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

Fig. 3. Spatial variability of length of optimum aggregation window across Australia for different climate indices and time scales of SPI (i.e. 12 M and 36 M SPI). White area indicates that data are not considered at these grid locations.

Table 2 Model performance in terms of MSE (NSE) for 12-month and 36-month SPI at four selected grid locations. ‘N’ and ‘D’ represent Niño3.4 and DMI respectively and the integers indicate the length (L) of the aggregation windows. Grid location

Model

12 month SPI

36 month SPI

Predictor

MSE (NSE)

Order-1

Order-2

Predictor

MSE (NSE)

Order-1

Order-2

G1 (35°S,150°E)

M1 M2 M3 M4 MM Reference model

N6 D3 D6 N3

D42 N12 N18 D56

0.29 0.26 0.27 0.29 0.28 0.37

(0.41) (0.18) (0.20) (0.34) (0.44) (0.10)

N56 D3 D6 N12

D76 N42 N36 N63

0.14 0.13 0.14 0.15 0.15 0.40

(0.57) (0.24) (0.31) (0.21) (0.59) (0.04)

G2 (20°S,146°E)

M1 M2 M3 M4 MM Reference model

N18 D3 D6 N24

D42 N12 N30 D30

0.15 0.16 0.16 0.16 0.17 0.40

(0.47) (0.33) (0.23) (0.46) (0.53) (0.02)

N42 N3 N6 D3

D42 N36 N56 N32

0.20 0.19 0.19 0.18 0.19 0.49

(0.64) (0.41) (0.44) (0.40) (0.66) (0.04)

G3 (14°S,132°E)

M1 M2 M3 M4 MM Reference model

N12 D30 D24 N18

D80 D76 N24 D63

0.29 0.30 0.31 0.27 0.28 0.44

(0.51) (0.16) (0.28) (0.47) (0.56) (0.01)

D80 N3 N6 D76

N56 D56 D63 N24

0.13 0.13 0.14 0.17 0.15 0.37

(0.68) (0.55) (0.56) (0.70) (0.73) (0.01)

G4 (34°S,116°E)

M1 M2 M3 M4 MM Reference model

D76 N3 N6 D80

N12 D56 D42 D18

0.25 0.17 0.17 0.25 0.20 0.37

(0.44) (0.33) (0.31) (0.39) (0.50) (0.01)

D80 N3 N6 D76

N56 D56 D63 N24

0.13 0.14 0.14 0.13 0.13 0.42

(0.74) (0.42) (0.47) (0.63) (0.75) (0.05)

6

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

Fig. 4. Density plots of observed and modelled SPI for 12 month and 36 month time scales. Blue line represents multi-model prediction with aggregated climate indices and red line represents reference model with raw climate indices. Back line represents density of observed SPI. Results for grid point G1 to G4 are represented by the plots from left to right. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

relative improvement of skill score for multi-model prediction compare to the reference model. Multi-model prediction of SPI provides improvement in skill score at varying magnitude compare to the reference model in approximately 95% grid locations. The average improvement of skill score is 23% and 28% for 12-month and 36-month SPI respectively across Australia. The improved skill for four selected locations (G1 to G4) were estimated as 5%, 11%, 17% and 22% for 12 month SPI and 12%, 16%, 28% and 41% for 36 month SPI. Extending the analysis to all of Australia, reduction of errors in the SPI predictions are observed in almost all locations (Fig. 6) for both the 12-month and 36-month SPI. The average reduction in MAE is 44% and 54% for 12-month and 36-month SPI respectively. Consistent with the results for the 4 sample locations, larger reductions are obtained for 36month SPI because as discussed earlier, the temporal aggregation of the climate indices better reflects the smoother SPI at this time scale. For example, greater than 50% reduction of error in SPI prediction is observed in 881 grid locations for 12-month SPI, whereas the count is 1705 grid locations for 36-month SPI. There is substantial spatial structure to the improvements offered by the multi-model ensemble. The largest improvements are evident over the arid inland of South Australia, Northern Australia and Western Australia. Rainfall over this region is highly persistent (Rocheta et al., 2014) which leads to more memory in SPI time series. A similar spatial pattern of improvement is observed for pdf skill score as well (Fig. 5). Thus the aggregated climate indices can represent the memory better through the smoothing of the aggregation process. This provides more skill in the predictions consistent with the results in Fig. 3. Finally, model performance of reproducing drought characteristics was assessed. First, drought events (i.e. SPI values is continuously lower than −1 at least for 2 months) were identified from the observed and modelled SPI series. Drought characteristics such as severity (sum of SPI values of each drought events), duration (i.e. average of duration of each drought events), and peak (average of largest negative values of SPI of each drought events) were estimated for each event and averaged. Drought frequency is the total number of drought events over the period 1910–2010. Fig. 7, shows the scatter plots of observed and modelled (reference model and multi-model ensemble) drought characteristics at each grid location across Australia. Results show a clear

Further improvements are obtained when model combination is applied. The multi-model prediction by combining multiple models with optimized weights reduces predictive uncertainty and improves the skill compared to any single model prediction. The multi-model approach allows for the interactions between the climate indices at a range of time scales to be incorporated into the predictions whilst still maintaining very simple model structures. The relative improvement offered by the multi-model prediction of SPI is different at the four selected grid locations which is not surprising given the regionally variable relationships between sustained rainfall anomalies and climate indices (Dijk et al., 2013; King et al., 2013; Kirono et al., 2010). The predictions from the multi-model ensemble compared to the observations and the reference model are shown in Fig. 4 where density plots are used to summarise the distribution of SPI predictions. In all four locations (Fig. S3 in the supporting document), the simple reference model consistently underestimates the frequency of droughts and wet periods with larger values of SPI (i.e SPI < −1 and SPI > 1 respectively) with the consequence that predictions of neutral conditions (i.e. 0.5 < SPI > −0.5) occurs too frequently. This is because of the inability of reference model to adequately capture the long-term variability of the observed SPI series. In contrast, the smoothing and isolation of important frequencies through the temporal aggregation means that the multi-model has much more skill in representing severe to extreme drought and wet conditions which are important for water resource management. In general, the improvements from the multimodel ensemble are larger for the 36-month SPI compared to 12-month SPI. This is because the smoothed version of the climate index time series is more valuable in representing the longer drought compared to the 12-month SPI. The multi-model estimates do not capture the full range of observed SPI values, with the probability of events with an absolute value of SPI more than 2 being underestimated. Further a probability density function (PDF) based skill score (Perkins et al., 2007) which measure the similarity between two PDFs is used to quantify the improvement in SPI prediction. This metric measures the common area between two PDFs by summing the minimum value of two distributions for different bins. Skill scores are estimated for both reference and multi-model SPI predictions against the observation at each grid locations across Australia. Fig. 5 shows the 7

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

Fig. 5. The percent improvement of PDF skill score by the aggregated climate indices compare to the raw climate indices in the SPI prediction at 12 month and 36 month time scales. White area represents no improvement or data are not considered at these grid locations.

representing drought, and in particular multi-year droughts, with strong regional variability found for the length of optimum aggregation window. The temporal aggregation also increases the number of predictors available for use with different window lengths providing different amounts of smoothing, leading to more structural uncertainty when constructing the prediction model. The PIC approach used here to select the predictors maximises the difference between predictors and provides a framework to create multi-model predictions. The multimodel ensemble clearly outperforms all of the constituent models as well as the reference model. An additional benefit is that it also allows the uncertainty in the predictions to be quantified. At each time step, the standard deviations of the four SPI predictions is calculated and used to construct prediction limits to quantify the uncertainty of the SPI. Fig. 8 shows the multi-model SPI prediction and corresponding prediction limits (95% significance level) for 12-month and 36-month

improvement when SPI is modelled with temporally-aggregated climate indices instead of raw climate indices. However, both models showed less skill for the grid locations with relatively longer and highly intense droughts (generally in the central arid region of Australia), although the multi-model ensemble still perform better than the reference model.

5. Discussion One of the important sources of uncertainty in any modelling study is the input data or predictors. This uncertainty is introduced partially by the choice of which subset of predictors to use, which can be considered as model structural uncertainty. Here a framework for identifying and using this predictor uncertainty has been presented in the context of providing improved drought predictions. Temporally aggregating climate indices was shown to maximise their skill in

Fig. 6. The percent reduction of errors (MSEs) by the aggregated climate indices compare to the raw climate indices in the SPI prediction at 12 month and 36 month time scales. White area represents no improvement or data are not considered at these grid locations. 8

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

Fig. 7. Observed and modelled (reference model and multi-model ensemble i.e. MME) drought characteristics at each grid location across Australia for 12-month (top row) and 36-month (bottom row) SPI. Red and blue circles represent each grid location across Australia. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

droughts should be considered carefully when used in decision making. This study considered four models to quantify the uncertainty in the prediction of SPI, chosen to provide a clear demonstration of the model framework. More models could be considered if needed as ascertained through a sensitivity analyses with more models which provided similar results to those presented above. Other sensitivity testing could explore the skill provided by longer aggregation windows, however data

cases at the four selected grid locations. While uncertainty slightly varies over the whole distribution of the multi-model SPI, uncertainty for the wet periods (positive SPI values) is generally higher than during the drier periods for 12 M SPI. Quantifying the uncertainty in drought predictions is useful for decision making to ensure sustainable water resource management practices. Uncertainty in the extremes is much higher compare to more neutral condition suggest that extreme

Fig. 8. Multi-model prediction of SPI and corresponding uncertainty band for 12 month and 36 month cases at selected grid locations. 9

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

12M SPI

0.23 yr

0.47yr

36M SPI

0.94 yr

1.87 yr

3.73 yr

7.4 yr

Fig. 9. Frequency bands of DMI for which correlation is maximum for 12 month and 36 month SPI across Australia. White area represents the grid locations where data are not considered.

distribution) as it is biased to the non-extreme values of the SPI series, especially when high values of the smoothing parameter (k) are used (Lall and Sharma, 1996). Fitting separate KNN models for non-extreme and extreme values of SPI would reduce the error in pdfs. Additionally, fitting a parametric probabilistic model might be useful to reduce error in pdfs, which is recommended for future work. Prediction of long term drought and wet conditions using climate indices estimated from SSTA suggests that we could consider future drought risk using SSTA information from general circulation models (GCMs). Directly projecting future drought and wet anomalies from GCM simulated precipitation is questionable due to the inability of these models to correctly simulate the low frequency variability of precipitation (Johnson and Sharma, 2011; Rocheta et al., 2014). In contrast, GCMs are reasonably reliable in predicting sea surface temperature (SST) especially over the tropical Indian and Pacific Ocean (Song et al., 2014; Wang et al., 2015). In addition, prediction of SST could be further improved by the multi-model combination (Kaiser Khan et al., 2014). Thus, this study provides a new insight into approaches that may improve future drought quantification and projections.

availability becomes limited as the aggregation window and SPI time period are increased. The spatial distribution of optimum aggregation window lengths across Australia provides a clear representation of the regional variability of teleconnections across the country. Locations with the longest aggregation windows are generally in the arid inland of Australia. Rainfall in these areas has high persistence with many dry months leading to longer spells of drought. Therefore, the estimated SPI time series are quite smoother with the highest variance at low frequencies, for example multi-year to decadal compared, to the other locations in Australia. So, it is expected that the low frequency components of climate indices would provide stronger relationships. This is prominently observed in case of DMI. Correlation of 12 month and 36 month SPI with DMI at different frequency bands as shown in Fig. 9 indicates that the frequency band of maximum correlation varies in space and reasonably resemble with the spatial distribution of length of optimum aggregation window. Raw DMI series are decomposed by discrete wavelet transform and the filtered time series corresponding to selected frequency bands are reconstructed using inverse discrete wavelet transform. This is followed by identifying which frequency band provides maximum correlation with SPI. Similarly, spatial variability of optimum aggregation windows as shown in Fig. 3, over the far inland multi-year frequency band of DMI, shows strongest relationships with 12 month SPI. Other parts of the country are dominated by inter-annual frequency bands. In case of 36 month SPI, multi-year frequency band of DMI is found as the most dominating over most parts of Australia. This suggests that aggregation of climate indices acts to enhance the strength of relationships with SPI by isolating frequency bands of interest and provides a better prediction of SPI compare to raw climate indices. In this study, we have considered SPI at 12-month and 36-month time scales because droughts of most concern for water resource management in Australia and many other parts of the world occurs over the year to multi-year periods. The reference model would be expected to perform better for shorter duration droughts as the value of aggregated climate indices would be lower in representing these events. While the aggregated model reasonably reduce the errors in SPI series compare to the reference model, the gap between the pdfs of the observed and modelled SPI is still high and further work is needed to improve it. The model has a tendency to predict neutral conditions (i.e. 0.5 < SPI > −0.5) more frequently and the extreme conditions occur less frequently compare to the observations. This might be due to the limited ability of the KNN model to capture the extremes (i.e. tail of the

6. Conclusions In this study, we investigate the potential of aggregated climate indices to enhance the prediction of prolonged drought and wet anomalies over the concurrent raw climate indices using a drought index (i.e. SPI) as response. Additionally, alternate constituent models with different subsets of predictors are developed to explore the uncertainty in the SPI predictions and a model combination approach is used to reduce the structural uncertainty. The results show that the optimum aggregation window length for varies substantially over space. This reveals the importance of using multiple constituent models with unique set of predictors to explore the extent of uncertainty. The methodology developed in this study is found useful to represent the prediction uncertainty originated from markedly different predictive models that exhibit different climatic combinations. Clear improvements are evident in the SPI prediction with aggregated climate indices over the case where concurrent raw climate indices are used as predictors. The multi-model prediction framework reduces the structural uncertainty and further improves the prediction efficiency. Relatively larger improvements are observed for more extreme drought and wet anomalies. 10

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

The methods developed in this study are can be applied to predict the response of any dynamic system that is influenced by the low frequency variability of climatic modes. One such application might be to forecast and project reservoir inflows for sustainable management and operation of large water storages.

47 (4), W04508. https://doi.org/10.1029/2010WR009272. Jones, D.A., Wang, W., Fawcett, R., 2009. High-quality spatial climate data-sets for Australia. Aust. Meteorol. Oceanogr. J. 58 (4), 233. Kaiser Khan, M.Z., Mehrotra, R., Sharma, A., Sankarasubramanian, A., 2014. Global sea surface temperature forecasts using an improved multimodel approach. J. Clim. 27 (10), 3505–3515. Khan, M.Z.K., Sharma, A., Mehrotra, R., 2017. Global seasonal precipitation forecasts using improved sea surface temperature predictions. J. Geophys. Res.: Atmos. 122 (9), 4773–4785. Khan, M.Z.K., Sharma, A., Mehrotra, R., Schepen, A., Wang, Q., 2015. Does improved SSTA prediction ensure better seasonal rainfall forecasts? Water Resour. Res. 51 (5), 3370–3383. King, A.D., Alexander, L.V., Donat, M.G., 2013. Asymmetry in the response of eastern Australia extreme rainfall to low-frequency Pacific variability. Geophys. Res. Lett. 40 (10), 2271–2277. Kirono, D.G., Chiew, F.H., Kent, D.M., 2010. Identification of best predictors for forecasting seasonal rainfall and runoff in Australia. Hydrol. Process. 24 (10), 1237–1247. Lall, U., Sharma, A., 1996. A nearest neighbor bootstrap for resampling hydrologic time series. Water Resour. Res. 32 (3), 679–693. Li, J., Zhou, S., Hu, R., 2016. Hydrological drought class transition using SPI and SRI time series by loglinear regression. Water Resour. Manage. 30 (2), 669–684. McKee, T., Doesken, N., Kleist, J., 1995. Drought monitoring with multiple time scales. In: Proceedings of the 9th Conference on Applied Climatology, American Meteorological Society Dallas, Boston, MA: 233–236. Mehrotra, R., Sharma, A., 2006. A nonparametric stochastic downscaling framework for daily rainfall at multiple locations. J. Geophys. Res.: Atmos. 111 (D15). Mehrotra, R., Sharma, A., 2016. A multivariate quantile-matching bias correction approach with auto-and cross-dependence across multiple time scales: implications for downscaling. J. Clim. 29 (10), 3519–3539. Mishra, A.K., Singh, V.P., 2010. A review of drought concepts. J. Hydrol. 391 (1), 202–216. Murphy, B.F., Timbal, B., 2008. A review of recent climate variability and climate change in southeastern Australia. Int. J. Climatol. 28 (7), 859–879. Nahar, J., Johnson, F., Sharma, A., 2017. Assessing the extent of non-stationary biases in GCMs. J. Hydrol. 549, 148–162. Naresh Kumar, M., Murthy, C., Sesha Sai, M., Roy, P., 2009. On the use of Standardized Precipitation Index (SPI) for drought intensity assessment. Meteorol. Appl. 16 (3), 381–389. Özger, M., Mishra, A.K., Singh, V.P., 2009. Low frequency drought variability associated with climate indices. J. Hydrol. 364 (1), 152–162. Perkins, S., Pitman, A., Holbrook, N., McAneney, J., 2007. Evaluation of the AR4 climate models’ simulated daily maximum temperature, minimum temperature, and precipitation over Australia using probability density functions. J. Clim. 20 (17), 4356–4376. Rajagopalan, B., Lall, U., Zebiak, S.E., 2002. Categorical climate forecasts through regularization and optimal combination of multiple GCM ensembles. Mon. Weather Rev. 130 (7), 1792–1811. Rashid, M.M., Beecham, S., 2019a. Characterization of meteorological droughts across South Australia. Meteorol. Applications. https://doi.org/10.1002/met.1783. Rashid, M.M., Beecham, S., 2019b. Development of a non-stationary Standardized Precipitation Index and its application to a South Australian climate. Sci. Total Environ. 657, 882–892. Rashid, M.M., Beecham, S., Chowdhury, R.K., 2015a. Assessment of trends in point rainfall using Continuous Wavelet Transforms. Adv. Water Resour., Elsevier 82, 1–15. Rashid, M.M., Beecham, S., Chowdhury, R.K., 2015b. Statistical downscaling of CMIP5 outputs for projecting future changes in rainfall in the Onkaparinga catchment. Sci. Total Environ. 530, 171–182. Rashid, M.M., Beecham, S., Chowdhury, R.K., 2015c. Statistical downscaling of rainfall: a non-stationary and multi-resolution approach. Theor. Appl. Climatol., Springer 1–15. https://doi.org/10.1007/s00704-015-1465-3. Rashid, M.M., Johnson, F., Sharma, A., 2018. Identifying sustained drought anomalies in hydrological records: a wavelet approach. J. Geophys. Res.: Atmos. Rayner, N., et al., 2003. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. J. Geophys. Res.: Atmos. 108 (D14). Ren, W., Wang, Y., Li, J., Feng, P., Smith, R.J., 2016. Drought forecasting in Luanhe River basin involving climatic indices. Theor. Appl. Climatol. 130 (3–4), 1133–1148. https://doi.org/10.1007/s00704-016-1952-1. Risbey, J.S., Pook, M.J., McIntosh, P.C., Wheeler, M.C., Hendon, H.H., 2009. On the remote drivers of rainfall variability in Australia. Mon. Weather Rev. 137 (10), 3233–3253. Robertson, A.W., Lall, U., Zebiak, S.E., Goddard, L., 2004. Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Mon. Weather Rev. 132 (12), 2732–2744. Rocheta, E., Sugiyanto, M., Johnson, F., Evans, J., Sharma, A., 2014. How well do general circulation models represent low-frequency rainfall variability? Water Resour. Res. 50 (3), 2108–2123. Rogers, C.D.W., Beringer, J., 2017. Describing rainfall in northern Australia using multiple climate indices. Biogeosciences 14 (3), 597. Saji, N., Goswami, B., Vinayachandran, P., Yamagata, T., 1999. A dipole mode in the tropical Indian Ocean. Nature 401 (6751), 360. Seibert, M., Merz, B., Apel, H., 2017. Seasonal forecasting of hydrological drought in the Limpopo Basin: a comparison of statistical methods. Hydrol. Earth Syst. Sci. 21 (3), 1611. Sharma, A., Chowdhury, S., 2011. Coping with model structural uncertainty in medium-

CRediT authorship contribution statement Md Mamunur Rashid: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Visualization, Writing - original draft. Ashish Sharma: Conceptualization, Methodology, Writing - review & editing. Fiona Johnson: Conceptualization, Methodology, Writing - review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This research was funded by the Australian Research Council, Australia (LP150100548), and the Crown lands & Water Division, Department of Industry, NSW, Australia. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.jhydrol.2019.124419. References Angelidis, P., Maris, F., Kotsovinos, N., Hrissanthou, V., 2012. Computation of drought index SPI with alternative distribution functions. Water Resour. Manage. 26 (9), 2453–2473. Anshuka, A., van Ogtrop, F.F., Willem Vervoort, R., 2019. Drought forecasting through statistical models using standardised precipitation index: a systematic review and meta-regression analysis. Natural Hazard. 97 (2), 955–977 doi.org/10.1007/s11069019-03665-6. Awan, J.A., Bae, D.H., 2016. Drought prediction over the East Asian monsoon region using the adaptive neuro-fuzzy inference system and the global sea surface temperature anomalies. Int. J. Climatol. 36 (15), 4767–4777. Cai, W., Cowan, T., 2013. Southeast Australia autumn rainfall reduction: A climatechange induced poleward shift of ocean-atmosphere circulation. J. Clim. 26 (23), 189–205. Cai, W., van Rensch, P., Cowan, T., Hendon, H.H., 2011. Teleconnection pathways of ENSO and the IOD and the mechanisms for impacts on Australian rainfall. J. Clim. 24 (15), 3910–3923. Chowdhury, S., Sharma, A., 2009. Long-range Nino-3.4 predictions using pairwise dynamic combinations of multiple models. J. Clim. 22 (3), 793–805. Chowdhury, S., Sharma, A., 2011. Global sea surface temperature forecasts using a pairwise dynamic combination approach. J. Clim. 24 (7), 1869–1877. Devineni, N., Sankarasubramanian, A., 2010. Improved categorical winter precipitation forecasts through multimodel combinations of coupled GCMs. Geophys. Res. Lett. 37 (24). Dijk, A.I., et al., 2013. The Millennium Drought in southeast Australia (2001–2009): Natural and human causes and implications for water resources, ecosystems, economy, and society. Water Resour. Res. 49 (2), 1040–1057. Farahmand, A., AghaKouchak, A., 2015. A generalized framework for deriving nonparametric standardized drought indicators. Adv. Water Resour. 76, 140–145. Fierro, A.O., Leslie, L.M., 2013. Links between central west Western Australian rainfall variability and large-scale climate drivers. J. Clim. 26 (7), 2222–2246. Fung, K., Huang, Y., Koo, C., Soh, Y., 2019. Drought forecasting: a review of modelling approaches 2007–2017. J. Water Climate Change. Guttman, N.B., 1999. Accepting the standardized precipitation index: a calculation algorithm. JAWRA J. Am. Water Resour. Assoc. 35 (2), 311–322. Hagedorn, R., DOBLAS-REYES, F.J., Palmer, T., 2005. The rationale behind the success of multi-model ensembles in seasonal forecasting–I. Basic concept. Tellus A 57 (3), 219–233. He, X., Guan, H., Zhang, X., Simmons, C.T., 2014. A wavelet-based multiple linear regression model for forecasting monthly rainfall. Int. J. Climatol. 34 (6), 1898–1912. https://doi.org/10.1002/joc.3809. Johnson, F., Sharma, A., 2011. Accounting for interannual variability: A comparison of options for water resources climate change impact assessments. Water Resour. Res.

11

Journal of Hydrology 581 (2020) 124419

Md. M. Rashid, et al.

seasonal rainfall variability in South Australia: potential for improving seasonal hydroclimatic forecasts. Int. J. Climatol. Ummenhofer, C.C., et al., 2009. What causes southeast Australia’s worst droughts? Geophys. Res. Lett. 36 (4). Ummenhofer, C.C., et al., 2011. Indian and Pacific Ocean influences on southeast Australian drought and soil moisture. J. Clim. 24 (5), 1313–1336. Vicente-Serrano, S.M., et al., 2011. A multiscalar global evaluation of the impact of ENSO on droughts. J. Geophys. Res.: Atmos. 116 (D20). Wang, G., Dommenget, D., Frauen, C., 2015. An evaluation of the CMIP3 and CMIP5 simulations in their skill of simulating the spatial structure of SST variability. Clim. Dyn. 44 (1–2), 95–114. Westra, S., Sharma, A., 2010. An upper limit to seasonal rainfall predictability? J. Clim. 23 (12), 3332–3351.

term hydro-climatic forecasting. Hydrol. Res. 42 (2–3), 113–127. Sharma, A., Mehrotra, R., 2014. An information theoretic alternative to model a natural system using observational information alone. Water Resour. Res. 50 (1), 650–660. Sharma, A., Mehrotra, R., Li, J., Jha, S., 2016. A programming tool for nonparametric system prediction using Partial Informational Correlation and Partial Weights. Environ. Modell. Software 83, 271–275. Sheffield, J., Wood, E.F., 2007. Characteristics of global and regional drought, 1950–2000: analysis of soil moisture data from off-line simulation of the terrestrial hydrologic cycle. J. Geophys. Res.: Atmos. 112 (D17). Shi, P., et al., 2016. Large-scale climate patterns and precipitation in an arid endorheic region: linkage and underlying mechanism. Environ. Res. Lett. 11 (4), 044006. Song, Z., Liu, H., Wang, C., Zhang, L., Qiao, F., 2014. Evaluation of the eastern equatorial Pacific SST seasonal cycle in CMIP5 models. Ocean Sci. 10 (5), 837. Tozer, C., Kiem, A., Verdon-Kidd, D., 2017. Large-scale ocean-atmospheric processes and

12