Ensemble prediction of regional droughts using ...

HYDROLOGICAL PROCESSES Hydrol. Process. (2013) Published online in Wiley Online Library (wileyonlinelibrary.com) DOI: 10.1002/hyp.9966

Ensemble prediction of regional droughts using climate inputs and the SVM–copula approach Poulomi Ganguli and M. Janga Reddy* Department of Civil Engineering, Indian Institute of Technology Bombay, Powai, Mumbai 400076, India

Abstract: In this study, the climate teleconnections with meteorological droughts are analysed and used to develop ensemble drought prediction models using a support vector machine (SVM)–copula approach over Western Rajasthan (India). The meteorological droughts are identified using the Standardized Precipitation Index (SPI). In the analysis of large-scale climate forcing represented by climate indices such as El Niño Southern Oscillation, Indian Ocean Dipole Mode and Atlantic Multidecadal Oscillation on regional droughts, it is found that regional droughts exhibits interannual as well as interdecadal variability. On the basis of potential teleconnections between regional droughts and climate indices, SPI-based drought forecasting models are developed with up to 3 months’ lead time. As traditional statistical forecast models are unable to capture nonlinearity and nonstationarity associated with drought forecasts, a machine learning technique, namely, support vector regression (SVR), is adopted to forecast the drought index, and the copula method is used to model the joint distribution of observed and predicted drought index. The copula-based conditional distribution of an observed drought index conditioned on predicted drought index is utilized to simulate ensembles of drought forecasts. Two variants of drought forecast models are developed, namely a single model for all the periods in a year and separate models for each of the four seasons in a year. The performance of developed models is validated for predicting drought time series for 10 years’ data. Improvement in ensemble prediction of drought indices is observed for combined seasonal model over the single model without seasonal partitions. The results show that the proposed SVM–copula approach improves the drought prediction capability and provides estimation of uncertainty associated with drought predictions. Copyright © 2013 John Wiley & Sons, Ltd. KEY WORDS

drought prediction; copulas; support vector regression; climate teleconnection; ensembles

Received 6 November 2012; Accepted 1 July 2013

INTRODUCTION Meteorological drought is defined as an extended period (a season, a year or more) of deficient rainfall relative to the statistical multiyear mean for a region (Druyan, 1996). Drought is multi-attribute in nature, characterizing several mutually correlated random variables such as severity, duration and peak. Although the magnitude of drought varies in different climatic regions of India, precipitation is the primary causal factor for droughts. In general, probabilities of droughts are higher in arid regions (e.g. western India) as compared with other parts of the country. According to recent studies, majority of districts in the northwestern part of the country, namely the states of Rajasthan, Gujarat, Jammu and Kashmir, Punjab and Haryana, have drought probabilities of more than 20% (Pai et al., 2010). Rajasthan (latitudes 23°3′N–30°12′N

*Correspondence to: Department of Civil Engineering, Indian Institute of Technology Bombay, Powai, Mumbai, 400076, India. E-mail: [email protected]

Copyright © 2013 John Wiley & Sons, Ltd.

and longitudes 69°30′E–78°17′E) is the largest state in India, covering an area of 34.22 million hectares, i.e. 10.5% of the country’s geographical area, but sharing only 1.15% of its water resources (Narain et al., 2005). Recent studies based on tree rings and other instrumental records have shown that many major droughts that occurred over different parts of the world are influenced by the atmospheric teleconnection El Nino Southern Oscillation (ENSO) phenomenon (D’Arrigo and Smerdon, 2008). On other hand, other climate indices such as the Indian Ocean Dipole (IOD) mode (Saji et al., 1999) induce an unusual rainfall distribution in surrounding regions. In the past, studies have noted that IOD and ENSO have complementarily affected the Indian summer monsoon rainfall (ISMR), and the IOD plays a key role in modulating the ISMR (Ashok et al., 2001). It was also noticed that when the correlation between ENSO and ISMR is high (low), the IOD–ISMR correlation is found to be low (high). Using an extended classification for ENSO and IOD events during 1877–2006, the study by Ummenhofer et al. (2011) noted that when El Niño events

P. GANGULI AND M. J. REDDY

coincide with positive IOD events, Indian Ocean conditions act to counter the El Niño’s drought-inducing subsidence by enhancing moisture convergence over the Indian subcontinent, which results in an average monsoon season. Other large-scale climate oscillations such as the North Atlantic Oscillation, the Pacific Decadal Oscillation (PDO) and the Atlantic Multidecadal Oscillation (AMO) act on longer timescales and sometimes interact with ENSO, which in turn may affect regional droughts (Sheffield and Wood, 2007). Özgar et al. (2009) investigated teleconnections of ENSO and PDO to droughts in Texas and found that arid regions have stronger correlations to climate anomalies as compared with subtropic humid regions. The study of Dong et al. (2006) suggested a possible link between AMO and ENSO variability – with warm AMO phase being related to weaker ENSO. During the warm phase of AMO, most part of the USA experienced less than normal precipitation including Midwest droughts in the 1930s and 1950s (Enfield et al., 2001). A significant role of AMO in causing severe droughts over the USA and Mexico was also investigated in other studies (McKee et al., 1993; Mendez and Magana, 2010). It was found that (Goswami et al., 2006) unlike that in the USA, AMO in its positive (i.e. in warm) phase increases the tropospheric temperature gradient over the entire Northern Hemisphere including Eurasia and may cause an increase in mean seasonal monsoon precipitation. From the Geophysical Fluid Dynamics Laboratory CM 2.1 climate model experiments and observations, Zhang and Delworth (2006) indicated that the warm (cold) phase of the AMO enhances (reduces) the Indian monsoon precipitation. For all India, southwest and central regions, the multidecadal wet period was found to be in phase with the positive AMO phase, and dry periods are in phase with the negative AMO phase, which indicates enhanced summer monsoon (June–September) precipitation during the warm phase of AMO (Joshy and Pandey, 2011). Although drought is a natural part of climate variability, the desire to mitigate its impact on the agriculture and water resources sectors has motivated the development of various drought prediction models. Drought forecasting is an essential tool for the implementation of appropriate drought mitigation strategies as well as in early warning and preparedness of droughts. Various drought prediction models have been developed by different researchers. Rao and Padmanabhan (1984) developed autoregressive models to forecast annual and monthly time series of Palmer Drought Severity Index (PDSI) for Iowa during 1930–1962 and for Kansas during 1887–1962. In their study, the performance of yearly drought forecast models was found to be less accurate as compared with that of the monthly forecast models. Loaiciga and Leipnik (1996) forecasted hydrological droughts using annual streamflow time series Copyright © 2013 John Wiley & Sons, Ltd.

in a coastal semiarid basin of Central California on the basis of renewal theory. Lohani and Loganathan (1997) developed a nonhomogeneous Markov chain model to forecast PDSI values for the tidewater climate division, Virginia, to develop an early drought warning system. As traditional statistical forecasting models (such as multiple linear regression) are unable to capture nonlinearity and nonstationarity associated with drought forecasts, Kim and Valdés (2003) used wavelet transformation and neural networks to capture nonlinear features for forecasting PDSI time series. Mishra et al. (2007) developed a hybrid drought forecast model combining autoregressive integrated moving average and feedforward recursive multistep artificial neural network (ANN) model in the Kangsabati River basin, India. Researchers have shown that incorporating climate state information improves the performance of long-term-lead drought forecasts. Morid et al. (2007) developed forecast models of the Effective Drought Index and SPI series for 1–12 months’ lead time using ANN and input variables as lagged values of drought indices, rainfall and climate indices – Southern Oscillation Index and North Atlantic Oscillation index in Tehran Province, Iran. An improvement in PDSI-based seasonal (September–December) drought prediction model was found by incorporating the IOD index together with Niño 3.4 Sea Surface Temperature index over Indonesia (D’Arrigo and Smerdon, 2008). Recent studies noted that the support vector machine (SVM) method, being a machine learning technique, has better capability in handling complex nonlinear relationships of hydrological problems and helps in improved forecasting of nonlinear time series (Asefa et al., 2006; Chiang and Tsai, 2012; Kalra and Ahmed, 2009). The SVM approach has greater potential to handle complexities of hydrological systems, and it may overcome some of the limitations of conventional and ANN methods such as local optima and overfitting (Shawe-Taylor and Cristianini, 2004). SVMs have a simple geometric representation, and it does not depend upon the dimensionality of the input space. The SVM has few successful applications in the hydrology and water resources field, such as for predicting streamflow (Asefa et al., 2006; Kalra and Ahmed, 2009), soil moisture estimation (Gill et al., 2006) and hydrologic drought prediction (Chiang and Tsai, 2012). From past literature, it was noted that SVM has not been applied sufficiently for drought problems, and there is a scope for testing its effectiveness for drought prediction; therefore, this study aims to develop drought prediction models using the support vector regression (SVR) method. Apart from forecasting the mean values of droughts, it is also important to specify the uncertainty associated with the forecasts. Thus, uncertainty assessment is an Hydrol. Process. (2013) DOI: 10.1002/hyp

ENSEMBLE PREDICTION OF REGIONAL DROUGHTS

multiple timescales such as 3, 6, 9 or 12 months to capture different drought states ranging from short-term, medium-term and long-term drought conditions. In this study, SPI at a 6-month timescale (SPI-6) is chosen for drought modelling, which can be useful for assessing seasonal drought conditions. The details of computation of the SPI can be found in McKee et al. (1993) and Janga Reddy and Ganguli (2012). A drought period is identified when the SPI value falls below a threshold level, which is taken as the 20th percentile value (Svoboda et al., 2002). Figure 1 illustrates drought characteristic identification using SPI. Drought length or duration (D) is taken as the number of consecutive time intervals (months) where SPI remains below the threshold level. Drought severity (S) is the cumulative values of SPI within the drought duration, i.e. Si ¼ ∑D t¼1 SPI i;t , i = 1, …, n. Drought peak (I) is the absolute value of the minimum value taken by SPI over the duration of the drought. Least square support vector regression

Support vector machines are machine learning methods that work on the basis of statistical learning theory by establishing functional dependence of the dependent variable on n-dimensional input variables. In this study, the least square SVR (LS-SVR) method is adopted for forecasting the drought index. In LS-SVR, equality type constraints are considered with a formulation in a least square sense. The solution involves solving a set of linear equations with an objective function (Suykens and Vandewalle, 1999) 1 γ n min lðw; b; ξ Þ ¼ kwk2 þ ∑ ξ 2i w;b;ξ 2 2 i¼1

2.00

1.00

Drought Events

t

0 Threshold level ti -1.00

Drought modelling using the Standardized Precipitation Index

The SPI is a normal quantile transformation applied to a fitted parametric distribution of precipitation time series (McKee et al., 1993). This SPI can characterize drought at Copyright © 2013 John Wiley & Sons, Ltd.

-3.00

te 2

1

-2.00

METHODOLOGY

(1)

3.00

SPI

important component of any prediction model, but many studies have neglected the uncertainty component and focused only on prediction of mean droughts. Only few studies have addressed this issue in the context of drought prediction: Hwang and Carbone (2009) developed autoregressive drought prediction models and quantified associated forecast uncertainty using residual sampling based on modified K-NN (where K denotes the number of neighbours and NN refers to the nearest neighbour) algorithm. Fundel et al. (2012) forecasted an 18-year (1991–2008) hydrological drought in Thur River, Switzerland, simulated from the Variable Resolution Ensemble Prediction System, a semi-distributed hydrological model employing a lead time of 1 month. But none of them have considered ensemble prediction of droughts at different lead times, incorporating climatic information. Thus, there is a greater scope for developing effective methods for predicting droughts and assessing associated forecast uncertainty. In the present study, an alternative approach is proposed for uncertainty assessment by constructing a copula-based joint probability model between observed and predicted drought index series, and the conditional distribution is used to generate ensemble of drought predictions and to specify the uncertainty in predictions. Recently, copulas have been successfully applied for probabilistic modelling and risk evaluation for several problems in the hydrology and water resources field such as drought frequency analysis (Song and Singh, 2010; Kao and Govindaraju, 2010; Sadri and Burn, 2012), uncertainty assessment pertaining to soil moisture modelling (Gao et al., 2007), uncertainty assessment in prediction of ISMR (Maity and Kumar, 2008) and downscaling of regional climate simulations (Laux et al., 2011). The present study proposes a copula-based approach for assessing the uncertainty associated with drought prediction. The main objectives of the present study are (i) to investigate the influence of climate forcing on meteorological droughts in one of the drought-prone regions in India, namely Western Rajasthan, and use the teleconnections information for predicting droughts, (ii) to develop SVRbased models to forecast droughts for lead times of 1, 2 and 3 months by utilizing its autoregressive property and lagged values of climate indices, (iii) to develop a copula-based approach for assessing the uncertainty associated with the SVR-based drought predictions.

D

Non drought duration (Dn)

Drought severity Si =

SPIi,t t=1

Drought duration (D)

Peak (I)

Time Interval (months)

Figure 1. Depiction of drought characteristics (severity, S; duration, D; and peak, I) based on Standardized Precipitation Index (SPI). Shaded regions show drought events. ti and te show initiation and termination of the drought events

Hydrol. Process. (2013) DOI: 10.1002/hyp


Subject to equality constraints, yi ¼ wT φðxi Þ þ b þ ξ i ;

i ¼ 1; …; n

(2)

where w is the weight function w ∈ Rn, b is a bias term (b ∈ R), γ is a regularization constant that controls the smoothness or flatness of the objective function and a positive real constant, ξ is a nonnegative slack variable and φ(•) is a high-dimensional feature space vector – a nonlinear function. SVR maps training vectors into a high-dimensional feature space via function φ(•). The optimization problem (Equations 1 and 2) can be solved by converting the problem into a dual-optimization problem and by using Lagrange’s multipliers. This yields the estimate of a function expressed as nSV f ðxi Þ ¼ ∑ αi k xi ; xj þ b, where k(xi, xj) represents the i¼1

kernel function, with xi and xj are the data points in the sample, and α denotes Lagrange’s multipliers. In this study, a widely used kernel function – the Gaussian radial basis function (RBF) – is used, which has the expression k(xi, xj) = exp((║xi xj║2/2σ2)), where σ > 0 is the width of the Gaussian RBF kernel. The generalization performance of SVR and its efficiency depends on the selected SVR hyperparameters. The grid search algorithm with a tenfold cross-validation is opted for the selection of SVM hyperparameters γ and σ2. Copula theory

The copula theory helps to represent the functional relationship between the n-dimensional distribution function and its univariate marginal cumulative distribution functions (CDFs). Considering a bivariate case, the two mutually correlated random variables X and Y with univariate continuous marginal distribution functions FX(x) = P(X ≤ x) and FY(y) = P(Y ≤ y), the link between the joint distribution FX,Y (x, y) and its copula C can be expressed using Sklar’s (1959) theorem as F X;Y ðx; yÞ ¼ C ½F X ðxÞ; F Y ðyÞ; θ ¼ C ðu; vÞ ∀x; y in R∈ð∞; ∞Þ

(3)

where the function C : [0, 1]2 → [0, 1] is called a twodimensional copula, with association parameter θ. In this study, three types of copulas, namely Clayton, Frank and Plackett families, are investigated for drought

modelling, and the best represented copula is used for ensemble prediction of droughts. Table I presents the expressions for copula functions and applicable range of copula parameter (θ) for modelling dependence. More details on theoretical background and properties of the copula families are available in Nelsen (2006). Parameters of the copula families are estimated using the maximum pseudolikelihood method (Genest and Favre, 2007). The procedure consists of transforming marginal variables into uniformly distributed vectors using its empirical distribution. Then, copula parameters are estimated using maximization of a pseudo-log-likelihood function. For the two-dimensional random vector X ∈ {Xi,1, Xi,2}, the procedure involves finding the ranks of the observations U ∈ {Ui,1, Ui,2} and estimating its empirical distribution using the expression U i;d ¼

1 n d ¼ 1; 2 and ∀i ¼ 1; 2; …; n ∑ I X j; d ≤X i;d n þ 1j¼1

(4) I(A) denotes an indicator function whose value is 1 if its logical argument A is true and 0 otherwise. Substituting the empirical CDFs into the bivariate copula density cθ yields log-likelihood function LU(θ). The parameter θ can be obtained by maximizing this rank-based pseudo-loglikelihood function numerically (Genest and Favre, 2007), n θ^ ¼ arg maxfLU ðθÞg ¼ arg max ∑ ln cθ U i;1 ; ; U i;2 i¼1 n Ri;1 Ri;2 ¼ arg max ∑ ln cθ ; ∀i∈f1; …; ng nþ1 nþ1 i¼1 (5) where Ri,1 and Ri,2 denote ranks of the observed data. The Cramér–von Mises distance between hypothesized copula families and the empirical copula is used for checking the adequacy of the copula models (Genest et al., 2009). The test consists of comparing the distance between rank-based empirical copula Cn(•) and the estimated parametric family of copulas Cθ(•). The Cramér–von Mises distance test statistics is given by (Genest et al., 2009) Sn ¼ ∫½0;12 nfC n ðu; vÞ C θ ðu; vÞg2 dC n ðu; vÞ n 2 ¼ ∑ C n U i;1 ; U i;2 C θ U i;1 ; U i;2

(6)

i¼1

Table I. Functional forms of copulas and their parameter range Copula family Clayton Gumbel–Hougaard Frank


Copula expression

θ∈

max([uθ + vθ 1]1/θ, 0) exp([(lnu)θ + (lnv)θ]1/θ) (1/θ)ln(1 + ((eθu 1)(eθv 1)/(eθ 1)))

[1, ∞)\{0} [1, ∞) θ ∈ R/{0}


ENSEMBLE PREDICTION OF REGIONAL DROUGHTS n ^ i;1 ≤u; U ^ i;2 ≤v ; u ; v∈½0; 1 is where C n ðu; vÞ ¼ 1n ∑ I U i¼1 the empirical copula computed from the observed data (x1, y1), …, (xn, yn). To check whether the selected model is appropriate or not, an approximate p-value for the test statistic Sn is obtained via large sample simulations using a parametric bootstrap method.

Ensemble prediction of droughts using copula. The ensembles of drought indices are generated through a copula-based approach, which involves the following steps:

1. Fitting the marginal distribution to observed SPI and predicted SPI time series using appropriate distributions, and computing CDF values v = FY(y) and u = FX(x), where Y denotes the observed SPI and X denotes the simulated SPI. 2. Fitting the joint distribution of observed and predicted time series using copula function C(u, v). 3. Selecting an appropriate copula model on the basis of analytical goodness-of-fit tests. 4. Generating random samples from the selected family of copulas using the inverse of the conditional distribution at u, i.e. v ¼ C V jU 1 ðvjuÞ, where CV│U(v│u) = ∂C(u, v)/∂u. 5. Obtaining the value of SPI (y) from the inverse probability integral transformation of v, y ¼ F 1 Y ðvÞ. Performance measures for ensemble drought prediction models

The performance of ensemble prediction models is compared using the continuous ranked probability score

(CRPS; Gneiting and Raftery, 2007) and Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) measures.

CASE STUDY AND DATA Study region

This study focused on prediction of droughts in western Rajasthan, which is one of the most drought-prone regions in India. The location map of the study region is shown in Figure 2. Western Rajasthan occupies about 57.31% of India’s total arid zone area. The climate is characterized by low, highly variable and ill-distributed rainfall, high wind speed, high evaporation losses and extremes of seasonal temperatures. Rainfall is not only low but also uncertain in the state. In Rajasthan, the southwest monsoon, which begins generally on the last week of June in the eastern parts, may last until the middle of September. Pre-monsoon showers normally begin during the middle of June, and post-monsoon precipitation occasionally occurs in October. In the winter season (i.e. during the north-eastern monsoon period – October to January, ONDJ), there are occasional little precipitation as compared with the monsoon season (June to September, JJAS). In most of the places, the highest monthly precipitation is observed during July and August (Rathore, 2005). Drought mainly occurs because of deficient rainfall during the monsoon months, which is characterized by SPI. This study initially investigates the influence of climate teleconnections (of ENSO,

Western Rajasthan

Figure 2. Location map of study region




AMO and IOD) on meteorological droughts in western Rajasthan and then uses for developing drought prediction models, so brief details of data used in this study are discussed below.

Data details

The monthly area-weighted precipitation data of nine rainfall stations in western Rajasthan meteorological subdivision (from January 1896 to December 2005) were obtained from the Indian Institute of Tropical Meteorology, Pune (http://www.tropmet.res.in). The Multivariate ENSO Index (MEI; Wolter and Timlin, 2011) is used for representing the ENSO phenomena. The data of MEI for a duration of 135 years (1896–2005) are obtained from the Climatic Prediction Center (http://www.cdc.noaa.gov/ people/klaus.wolter/MEI/). The positive values of the index represent El Niño, whereas negative values correspond to La Niña episodes. The data set of AMO is obtained for 135 years (1896–2005) from National Oceanic and Atmospheric Administration’s Earth System Research Laboratory (http://www.esrl.noaa.gov/psd/data/ timeseries/AMO/). The monthly IOD data for 46 years (1960–2005) are obtained from the Japan Agency for Marine-Earth Science and Technology website: http:// www.jamstec.go.jp/frcgc/research/d1/iod/e/iod/dipole_ mode_index.html.

RESULTS Association between ENSO and precipitation The relationship between seasonal precipitation (JJAS and ONDJ) and lagged ENSO climate indices is investigated using rank correlation analysis. Spearman’s rank correlations between standardized anomaly (i.e. deviations from mean divided by standard deviation) of JJAS precipitation and aggregated MEI series of periods January–February (JF), February-March (FM),…, July–August (JA) are computed during the period of 1896–2005. The results are presented in Table II. It was found that the JJAS season precipitation shows significant correlation with the JA season MEI as confirmed by smaller p-values of the estimate. On applying a two-tailed t-test, it is also observed that there is no significant correlation between the ONDJ season precipitation and ENSO indices. AMO–precipitation relation

Standardized anomalies of monthly precipitation data during the 135-year study period (1896–2005) are detrended and smoothed with an 11-year moving average in order to remove interannual variability from the time Copyright © 2013 John Wiley & Sons, Ltd.

series. The Spearman’s rank correlation between the smoothed time series of AMO and precipitation is found to be 0.134 with a p-value of 3.09e6, which is statistically significant at the 5% level. Figure 3 presents the decadal variability of standardized and detrended precipitation anomaly along with the decadal variability of the AMO index. It can be observed that AMO is in the positive phase (AMO+) during 1929–1961 and is in the negative phase (AMO) during 1962–1994. Analysis of droughts under ENSO and AMO phases

From the ENSO information, the climate events can be grouped into three phases: El Niño, La Niña and neutral events. The climate states are assumed to be mutually exclusive and exhaustive so that each year belongs to one particular ENSO episode. Every year from 1896 to 2005 is given an ENSO classification on the basis of the MEI index values during May to November (Wolter and Timlin, 2011). The Spearman’s rank correlations between SPI-6 and MEI time series at El Niño, La Niña and neutral phase years are found to be 0.14, 0.30 and 0.14 with corresponding p-values of 5.11e4, 2.3e8 and 0.008, respectively. All correlations are statistically significant at the 5% level, which indicates that drought in the study region is responsive to the ENSO phenomenon. Among 93 drought events identified, 43 events are classified under El Niño, 20 under La Niña and 30 under neutral episodes. Table III presents the summary statistics of droughts without and with accounting for the ENSO state. It can be observed that the droughts during the El Niño episodes are longer and severe in nature, whereas the reverse is observed for La Niña years’ droughts. Neutral phase droughts lie in between El Niño and La Niña phase droughts in terms of drought characteristics – severity, duration and peak. The number of ENSO phase droughts at AMO and AMO+ state is investigated and presented in Table IV, showing greater frequency of droughts during AMO as compared with AMO+. Table V compares two different AMO phases and the frequency of drought events under individual ENSO states in one cycle of AMO (between 1929 and 1994) comprising two contrasting 33-year periods as shown in Figure 3 for 1929–1961 (i.e. AMO+ phase) and 1962–1994 (i.e. AMO phase) and for the entire period where drought events can be classified under AMO+ and AMO phases. To test whether the frequency of droughts is modulated by the AMO and the individual ENSO state, a simple test of proportion (Hogg and Tanis, 1977) is employed. The test is performed under the assumption that the sampling distribution of the proportion of drought events in n number of ENSO years (El Niño/La Niña) within any AMO phase can be approximated by normal distribution with mean p^ and variance p^ð1 p^Þ=n. In the Hydrol. Process. (2013) DOI: 10.1002/hyp

— 0.037 (0.70) — 0.046 (0.63)


Note: Correlations that are significant at the 5% significance level are shown in bold; p-values of correlation are shown in brackets. MEI, Multivariate El Nino Southern Oscillation Index.

0.32 (