An integrated package for drought monitoring, prediction and analysis to aid drought modeling and assessment

Environmental Modelling & Software 91 (2017) 199e209

Contents lists available at ScienceDirect

Environmental Modelling & Software journal homepage: www.elsevier.com/locate/envsoft

An integrated package for drought monitoring, prediction and analysis to aid drought modeling and assessment Zengchao Hao a, *, Fanghua Hao a, Vijay P. Singh b, Wei Ouyang a, Hongguang Cheng a a b

Green Development Institute, School of Environment, Beijing Normal University, Beijing, 100875, China Department of Biological and Agricultural Engineering, Zachry Department of Civil Engineering, Texas A&M University, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 22 June 2016 Received in revised form 30 December 2016 Accepted 5 February 2017

Due to severe drought events and disastrous impacts in recent decades, substantial efforts have been devoted recently to drought monitoring, prediction and risk analysis for aiding drought preparedness plans and mitigation measures. Providing an overview of these aspects of drought research, this study presents an integrated R package and illustrates a wide range of its applications for drought modeling and assessment based on univariate and multivariate drought indices for both operational and research purposes. The package also includes statistical prediction of drought in a probabilistic manner based on multiple drought indicators, which serves as a baseline for drought prediction. The univariate and multivariate drought risk analysis of drought properties and indices is also presented. Finally, potential extensions of this package are also discussed. The package is provided freely to public to aid drought early warning and management. © 2017 Elsevier Ltd. All rights reserved.

Keywords: Drought monitoring Drought prediction Drought analysis Drought early warning Extremes

Software availability Name of software: drought Developer: Zengchao Hao Email: [email protected] Programming language: R License: This software is provided under the terms and conditions of the GNU GENERAL PUBLIC LICENSE Version 3 Availability: https://r-forge.r-project.org/projects/drought/ Cost: Free for non-commercial academic research 1. Introduction Drought is among the costliest natural hazards with disastrous impacts on society and ecosystems. Severe droughts in the past decades, such as the 2012 U.S. drought and 2010e2011 East Africa drought, led to huge economic losses or even famine around the world (Dutra et al., 2013; Smith and Katz, 2013; Hoerling et al., 2014). Reliable drought monitoring and early warning play an important role in coping with drought, which requires integrated drought monitoring, prediction and risk assessment in order to

* Corresponding author. No. 19, XinJieKouWai St., HaiDian District, Beijing 100875, China. E-mail address: [email protected] (Z. Hao). http://dx.doi.org/10.1016/j.envsoft.2017.02.008 1364-8152/© 2017 Elsevier Ltd. All rights reserved.

track the drought status, provide prediction information, and assess the risk associated with drought impacts. Drought can be classified mainly into four types, including meteorological, agricultural, hydrological, and socioeconomic drought (Dracup et al., 1980; Wilhite and Glantz, 1985), while other types of drought such as groundwater drought and ecological drought has also been used (Mishra and Singh, 2010). Scores of drought indices have been developed for the characterization of different drought types in past few decades (Heim, 2002; Mishra and Singh, 2010; Zargar et al., 2011; Hannaford et al., 2015; Hao and Singh, 2015; Bachmair et al., 2016), such as standardized precipitation index (SPI) (McKee et al., 1993) and the Palmer Drought Severity Index (PDSI) (Palmer, 1965), which have been widely used for drought characterizations. However, there is a lack of consistency on the universally accepted drought indicator, and the suitability of particular drought indices depends on a specific region, season and application. For example, for drought monitoring in the region where snowmelt has to be taken into account, such as western U.S., the Surface Water Supply Index (SWSI) (Shafer and Dezman, 1982) or other snow-related indices (Staudinger et al., 2014) would be the suitable choice. Meanwhile, drought conditions are generally related to multiple variables with deficit from various sources and an individual drought indicator may not suffice to characterize drought in different regions and seasons. This has resulted in a surge in the development of multivariate or composite

200

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209

drought indices to integrate different drought related variables or drought indices for efficient drought monitoring and early warning, such as the U.S. Drought Monitor (Svoboda et al., 2002), Vegetation Drought Response Index (VegDRI) (Brown et al., 2008), Optimal blended NLDAS drought index (OBNDI) (Xia et al., 2014), Reconnaissance Drought Index (RDI) (Tsakiris and Vangelis, 2005), Standardized Precipitation Evapotranspiration Index (SPEI) (Vicente-Serrano et al., 2010), Joint Deficit Index (JDI) (Kao and Govindaraju, 2010), Multivariate Standardized Drought Index (MSDI) (Hao and AghaKouchak, 2013), and Aggregated Drought Index (ADI) (Keyantash and Dracup, 2004). These univariate and multivariate drought indices provide useful tools for characterizing different aspects of drought from different perspectives for drought management. Early drought warning is important for drought preparedness to reduce losses, which can be achieved with drought prediction through the enhanced understanding of drought causes and evolution. In past few decades, various methods have been proposed for drought prediction, based on either dynamical or statistical methods (Mishra and Singh, 2011). Statistical methods generally rely on statistical relationships of historical records, such as regression model, time series modeling, artificial intelligence methods, probability distribution method, and Ensemble Stream€ flow Prediction (ESP) among others (Mishra and Singh, 2011; Ozger et al., 2012; Hao et al., 2016a). Dynamical drought prediction based on seasonal climate forecasts, such as NCEP Climate Forecast System Version 2 (CFSv2) (Saha et al., 2014) or North American MultiModel Ensemble (NMME) (Kirtman et al., 2014), coupled with land surface models have been increasingly used for predicting meteorological, agricultural and hydrological drought (Wood and Lettenmaier, 2006; Luo and Wood, 2007; Yuan and Wood, 2013; Yuan et al., 2013; Nijssen et al., 2014; Shukla et al., 2014). Though great strides have been made for the understanding of drought causes and the development of prediction methods, challenges still exist, such as improvement of drought prediction skill beyond 1e2 months lead time (Pozzi et al., 2013; Hoerling et al., 2014; Wood et al., 2015; Schubert et al., 2016). To facilitate drought risk analysis in operational drought management, the statistical inference of drought properties or indices is desired. Based on the run theory, a drought event can be characterized with certain properties (Yevjevich, 1967), such as duration, severity, intensity and spatial extent. Statistical analysis of these drought properties plays an important role in drought risk assessments and water resources planning and management. For example, the return period of a drought duration or severity can be estimated by fitting a suitable probability distribution, such as geometric or gamma, for frequency analysis (Kendall and Dracup, 1992; Mathier et al., 1992). Since drought properties are generally mutually correlated, traditional drought risk assessments with the univariate frequency analysis may not be sufficient and joint modeling of multiple properties is therefore required. For this purpose, multivariate approaches have been used for statistical lez and Valde s, 2003; Salas analysis of drought properties (Gonza et al., 2005; Shiau, 2006; Nadarajah, 2009; Song and Singh, 2010; Xu et al., 2015). Moreover, different types of drought may occur simultaneously (e.g., joint deficit of precipitation, soil moisture or runoff) (Beersma and Buishand, 2004; Kao and Govindaraju, 2010; Pan et al., 2013; Ma et al., 2014) and a drought is often accompanied or aggravated by high temperature, low relative humidity or other extremes (e.g., the combination of low precipitation and high temperature or high evapotranspiration) (Lyon, 2009; Hao et al., 2013; Leonard et al., 2014; Hao and Singh, 2015; Cheng et al., 2016; Mo and Lettenmaier, 2016). Hence, to facilitate the modeling of multiple drought properties or indices for risk analysis, suitable dependence modeling techniques have to be adopted to

capture various dependence structures. Operational drought management requires accurate monitoring to track drought conditions, reliable prediction to aid early warning, and objective analysis to assess risk to reduce potential impacts of drought. The objective of this study therefore is to develop an integrated drought package to meet these needs to aid drought modeling and assessment using the R software (R Development Core Team, 2015). This paper is organized as follows. An overview of this package is provided in Section 2. Section 3 illustrates the application of different components of this package, followed by the conclusions and future development in Section 4.

2. Overview of the integrated package 2.1. Computation of drought indices 2.1.1. Univariate drought index A univariate drought index has been commonly developed based on a single hydroclimatic variable (e.g., SPI based solely on precipitation), which is generally constructed using the percentile, anomaly, or the standardization in different ways to meet the need for drought characterization. The commonly used standardization scheme, based on the computation of SPI, can also be applied to other variables, such as soil moisture, streamflow, runoff, groundwater, and snow melt to derive different types of Standardized Drought Index (SDI), including the Standardized Soil moisture Index (SSI), Standardized Runoff Index (SRI), Standardized Groundwater Index (SGI) and Standardized Snow Melt and Rain Index (SMRI) (Shukla and Wood, 2008; Bloomfield and Marchant, 2013; ~ ez et al., 2014; Staudinger et al., 2014; Hao Hao et al., 2014; Nún and Singh, 2015). The development and application of SDI has attracted much attention, since the SDI meets certain desired properties of drought indicators, such as statistical consistency and comparability at different spatial scales (Heim, 2002; Steinemann et al., 2015; Hao et al., 2016b). In the package, computation of the SDI of different time scales (e.g., 3-month, 6-month) is provided, which consists of SPI as a special case, and is introduced as follows. Denote the accumulated variable (e.g., precipitation, runoff, or soil moisture) of different time scales as the random variable X with the probability density function (PDF) f(X). A suitable parametric distribution, such as gamma or log normal (McKee et al., 1993; Sheffield et al., 2004; Stagge et al., 2015), can be used to fit the accumulated precipitation, which can be expressed as:

Zx GðxÞ ¼

f ðtÞdt

(1)

0

When the aggregated monthly precipitation x ¼ 0, the cumulative distribution function in equation (1) can be expressed as:

HðxÞ ¼ q þ ð1  qÞGðxÞ

(2)

where q is the frequency of occurrence of x ¼ 0. The standard normal distribution can be used to transfer the cumulative probability (or percentile) P¼H(x) to obtain the SDI:

SDI ¼ N 1 ðPÞ

(3)

where N is the standard normal distribution function with mean zero and standard deviation unity. Except for the parametric distribution, the nonparametric method, including kernel density estimation or empirical plotting position formula (e.g., Weibull, Gringorten, Harris, Hazen, Beard,

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209

California and others) (Gringorten, 1963; Makkonen, 2006; Fuglem et al., 2013) can also be used for the computation of drought indices. For example, the Weibull and Gringorten plotting position formulae can be expressed as:

Weibull

P¼

i nþ1

Gringorten

P¼

i  0:44 n þ 0:12

(4)

where n is the length of observations; i is the rank of observed values. The probability P can then be substituted into equation (3) to compute the SDI. 2.1.2. Multivariate drought index In the past decade, there have been extensive efforts in developing the multivariate drought index (MDI) that integrate a suite of drought related variables or indices for the multivariate or composite drought characterization (Hao and Singh, 2015). For a suite of variables Z1, Z2, …, Zn, the MDI can be expressed with a function F as:

MDI ¼ FðZ1 ; Z2 ; :::; Zn Þ

(5)

where Z1, Z2, …, Zn can be drought related variables or drought indices; the function F can be a linear combination, joint distribution, Principal Component Analysis (PCA) or others. The joint distribution is commonly used to derive the MDI, since it enables statistical characterizations of the joint or conditional behavior of two or more random variables. In this case, function F in equation (5) can be a parametric or nonparametric distribution in the multivariate case to model various dependence structures of multiple hydroclimatic variables (Hao and Singh, 2016). In the past decade, various transformations of the joint probability distribution have been employed to develop the MDI (Beersma and Buishand, 2004; Kao and Govindaraju, 2010; Ma et al., 2014; Huang et al., 2016; Hao et al., 2016b). In the following, the derivation of multivariate or composite drought indices based on the joint distribution is introduced. 2.2. Joint distribution and Kendall distribution function Denoting two random variables X and Y representing hydroclimatic variables (e.g., precipitation, soil moisture and runoff) or indices (e.g., PDSI), the joint distribution F(X,Y) of the two random variables can be expressed as:

Fðx; yÞ ¼ PðX  x; Y  yÞ

(6)

Equation (6) defines the probability of two variables lower than certain values, which represent the case when both variables lower than certain values is of interest (e.g., both precipitation and soil moisture is low). Other cases based on the combination of two variables such as P(X  x or Y  y) can also be defined based on the joint distribution function. Though a variety of multivariate distribution functions are available for modeling the dependence structure in equation (6), recent decades have witnessed a surge in applications of the copula (Nelsen, 2006; Joe, 2014) in constructing the joint distribution, which is advantageous in that the construction of joint distributions to model the dependence is independent of marginal distribution types. Assuming U¼FX(x) and V¼FY(y) as marginal distributions of the random vector (X, Y), the joint distribution of the random vector can be expressed with the copula C as:

Fðx; yÞ ¼ PðX  x; Y  yÞ ¼ Cðu; vÞ

(7)

There are a variety of copula families, such as empirical copula,

201

meta-elliptical copula, Archimedean copula (e.g., Gumbel, Clayton, Frank), extreme-value copula, and vine copula (Genest and Favre, 2007; Kurowicka and Joe, 2011; Hao and Singh, 2016). For example, the Gumbel copula Cg(u,v) in the bivariate case can be expressed as:

i1=q o n h Cg ðu; vÞ ¼ exp  ð  loguÞq þ ð  log vÞq

(8)

where q is the parameter of the Gumbel copula. However, in many cases parametric forms of the joint distribution are not available and nonparametric forms are therefore resorted to. Let x(i) and x(j) (1  i, j  n) be order statistics of the samples of size n. The empirical copula of two random variables X and Y can be derived as (Nelsen et al., 2003):

C

0



i j ; n n

 ¼

  # x  xðiÞ ; y  yðjÞ n

¼

n1 n

(9)

where n1 is the number of samples (xk, yk) satisfying (xk  x(i) and yk  y(j)), 1  k  n. The Kendall distribution function is a commonly used concept with respect to the joint distribution. The copula C(u,v) with margins U and V gives the probability measure with C(u,v) ¼ p, which corresponds to infinite combinations of u and v leading to the same probability p (e.g., marginal pairs (u1,v1) and (u2,v2) may lead to the same probability p). In the bivariate case, the Kendall distribution function Kc defines a measure of the set [(u,v)2[0,1]2j C(u,v)q]. Specifically, for a continuous random vector (X, Y) with marginals U and V, the Kendall distribution function K(p) can be defined as (Genest et al., 2007):

KðpÞ ¼ PðCðu; vÞ  pÞ

(10)

where p is a probability level. The analytical form of the Kendall distribution function is available for certain copulas, such as Archimedean family, but does not exist for all copulas (in which case simulation approach can be , 2009; Salvadori et al., 2011). resorted to) (McNeil and Neslehova For example, the Kendall distribution function for the Gumbel copula can be expressed as:

  log p KðpÞ ¼ p 1 

(11)

q

Since analytical expressions of Kendall distribution may not exist for many copulas, the empirical Kendall distribution function can be used, which can be expressed as (Nelsen et al., 2003):

KC ; ðpÞ ¼

n2 n

(12) 0

where n2 is the number of samples satisfying C’ (i/n, j/n)p (C is the empirical copula). Equation (12) can be used for the construction of the Kendall distribution function in a nonparametric way. 2.3. Derivation of multivariate drought index The development of the MDI based on the joint distribution function can be summarized into two cases, either based on the joint probability or the Kendall distribution function. The joint probability (or its transformation) of two variables X and Y characterizes the joint behavior and can be regarded as a measure of the drought condition. For example, the joint probability P(X  x, Y  y) ¼ p can be employed as a bivariate drought indicator (Beersma and Buishand, 2004; Ramadas and Govindaraju, 2016). A

202

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209

low value of p (say p ¼ 0.2) indicates an overall drought condition, while a high value of p (say 0.8) implies a wet condition. Apart from the joint probability p (regarded as the indicator of the joint deficit), it would be of particular interest to know the probability of the random vector with C(u,v)q, which can also be used to derive the index using the Kendall distribution (Kao and Govindaraju, 2010). Similar to the SDI, it is generally convenient to express the drought indicator in a standardization manner for statistical consistency across space and time and thus the standardization has been used to derive the MDI either based on the joint probability (Hao and AghaKouchak, 2013) or the Kendall distribution (Kao and Govindaraju, 2010). To summarize, the development of a multivariate or composite drought index using the joint distribution (denoted as the Joint Drought Severity Index, JDSI) can be achieved through the joint probability or Kendall distribution of multiple hydroclimatic variables or indices. The JDSI in these two cases (denoted as JDSIJ and JDSIK) can be expressed in a standardized manner as:

JDSIJ ¼ N 1 ½PðX  x; Y  yÞ 

(13)

JDSIK ¼ N 1 ½PðX  x; Y  yÞ  p  ¼ N 1 ½K ðpÞ 

(14)

The development of JDSI in equations (13) and (14) provides a general framework to define a multitude of MDIs. For example, the recently developed Multivariate Standardized Drought Index (MSDI) (Hao et al., 2014) and Joint Deficit Index (JDI) (Kao and Govindaraju, 2010) or others are all special cases of JDSIJ or JDSIK. One issue related to the JDSI is that joint probability of the multiple variable are generally not uniform within [0 1] (and the MDI is not normally distributed) (Hao et al., 2016b). Thus the transformation can be employed to rescale the MDI to be normally distributed, for which the normal quantile transformation (NQT) can be used:

NQTðZÞ ¼ N1 ½GðZÞ

(15)

where Z is a random variable, which can be the MDI; G is a parametric or nonparametric distribution. This procedure is essentially the percentile approach, which has been commonly used for the classification of drought severity into categories for drought management (Svoboda et al., 2002; Steinemann, 2003; Goodrich and Ellis, 2006; Mo and Lettenmaier, 2014), in the multivariate case. 2.4. Statistical drought prediction 2.4.1. Baseline drought prediction In this package, the baseline drought prediction method proposed by Lyon et al., (2012) is employed, which is similar to the concept of the Ensemble Streamflow Prediction (ESP) method (Twedt et al., 1977; Day, 1985). By incorporating the prior state of precipitation (initial condition) and seasonal cycle of the climatological precipitation (climatic condition), the drought prediction can be achieved, based on the inherent persistence of drought indicator SPI. Since this ESP based prediction method has been shown to hold potential for drought prediction through comparison with the dynamical prediction (Quan et al., 2012; Mo and Lyon, 2015), we introduce this method for drought prediction in both the univariate case based on SPI (Lyon et al., 2012) (or other SDIs) and the multivariate case based on MDIs (Hao et al., 2014, 2016c). Assume monthly precipitation (or other variables) of n-year historical records is known and the prediction of the SPI of the smonth time scale (e.g., s ¼ 6 represents 6 month accumulation) for the target month in the year nþ1 is required. In the following, APp,q stands for the precipitation of the month q (q ¼ 1,2, …, 12) in the

year p (p ¼ 1,2, …, nþ1). For the s month accumulation of precipitation ending with month i, the accumulated precipitation APt,i of the year t can be expressed as:

APt;i ¼ Pt;isþ1 þ ::: þ Pt;i1 þ Pt;i

(16)

The accumulated precipitation of l (l < s) month ahead in the future can be expressed as:

APt;iþl ¼ Pt;iþlsþ1 þ ::: þ Pt;iþl1 þ Pt;iþl

(17)

The accumulated precipitation APt,iþ1 can also be expressed as (Lyon et al., 2012):

APt;iþl ¼ Ct;il þ Dt;il

(18)

where Ct,il is the “common” term of the current observation that can be expressed as:

Ct;il ¼ Pt;iþlsþ1 þ ::: þ Pt;i1 þ Pt;i

(19)

and Dt,il is the “disjoint” term in the future that can be expressed as:

Dt;il ¼ Pt;iþ1 þ ::: þ Pt;iþl

(20)

The common term Cil can be regarded as the initial condition, which is already known and provides the initial state of the accumulated precipitation for the prediction of SPI in the target month i þ l. The disjoint term Dil can be regarded as the climatic information, which is unknown and provides the climatology information of precipitation in the target month i þ l. For the dynamical prediction, the prediction of the disjoint term Dt,il of precipitation is achieved from seasonal climate forecast (Yoon et al., 2012; Yuan and Wood, 2013). For the statistical prediction based on ESP, the prediction of disjoint term Dt,il is achieved by resampling from the historical record (from the 1st year to the year t-1), assuming all historical records will occur with the same likelihood in the target month. For example, when the 1- month lead prediction of the 6 month SPI for August in the year t is performed (APt,8), the common term Ct,il includes the precipitation of the first 5 months (i.e., Pt,3, …, Pt,6, Pt,7) from March to July of the current year t and the disjoint term Dt,il includes monthly precipitation in August Pt,8, which is predicted as an ensemble of historical monthly observations in August (i.e., P1,8, …, Pt-1,8). Here we take the two month lead prediction (l ¼ 2) of SPI, based on the 6-month accumulation (s ¼ 6), as an example to illustrate the prediction. Based on n-year historical records, the 2-month lead prediction of the accumulated precipitation of the target month iþ2 in the year t ¼ nþ1 can then be expressed from equation (17) as:

APnþ1;iþ2 ¼ Pnþ1;i3 þ ::: þ Pnþ1;i þ Pnþ1;iþ1 þ Pnþ1;iþ2

(21)

As stated before, the disjoint terms Pnþ1,iþ1 and Pnþ1,iþ2 need to be predicted for the prediction of APnþ1,iþ2, which can be achieved through resampling from the precipitation in the target month i þ 1 and i þ 2 in n-year historical records. Specifically, an ensemble of n members of APnþ1,iþ2 can be obtained by replacing Pnþ1,iþ1 and Pnþ1,iþ2 in equation (21) with each element of the historical record (P1,iþ1, …, Pn,iþ1) and (P1,iþ2, …, Pn,iþ2), which can be expressed as:

AP ðjÞ nþ1;iþ2 ¼ Pnþ1;i3 þ ::: þ Pnþ1;i þ Pj;iþ1 þ Pj;iþ2 ¼ 1; 2; :::; n

j (22)

where AP(j)nþ1,iþ2 is the jth member of the predicted ensemble of APnþ1,iþ2. Accordingly, the 2-month lead prediction of the accumulated

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209

precipitation (APnþ1,iþ2) of the target month in the year nþ1 is an ensemble of n members (i.e., AP(1)nþ1,iþ2, …, AP(n) nþ1, iþ2). Denote Z a random variable of the 6-month accumulated precipitation ending with month i þ 2. The jth member AP(j)nþ1,iþ2 (j ¼ 1,2, …, n) can be combined with the accumulated precipitation in the n-year historical record to form a sample zj¼(AP1,iþ2, AP2,iþ2, …, APn,iþ2, AP(j)nþ1,iþ2). The SPI can be computed, based on the sample zj, which can be expressed as:

SPI ¼ N 1 ½FZ ðzÞ

(23)

where Fz is the distribution function of the sample. The predicted jth member of SPI of the target month can then be obtained from equation (23), which is the one corresponding to the accumulated value AP(j)nþ1,iþ2. Totally n member of SPI prediction can be obtained based on n-year historical records. 2.4.2. Multivariate drought prediction The drought prediction with multivariate drought indices has also been highlighted in recent years, since the multivariate prediction from multiple indices may provide a complimentary prediction skill to aid drought early warning (Hao et al., 2014; Behrangi et al., 2015; Mo and Lyon, 2015; Hao et al., 2016c). The prediction of MDI can be obtained with ESP by sampling multiple variables of the same year from historical records (termed as Multivariate Ensemble Streamflow Prediction, MESP), therefore preserving the cross dependence among different variables (Hao et al., 2016c). The predicted accumulated variables in equation (22) or SDI in equation (23) can be substituted into the function F in equation (5) to obtain the multivariate drought prediction. Similar to the prediction of the univariate SDI, these procedures also result in an ensemble of the predicted MDI, which can be used for probabilistic drought prediction. 2.4.3. Severity and probability of drought prediction Based on the prediction of n ensemble members of drought indices (DI) (e.g., SDI or MDI), the ensemble mean or median can be used as the predicted drought severity. In other words, the prediction of drought severity (DIp) can be expressed as:

DIp ¼ meanðDI1 ; :::; DIn Þ

or

DIp ¼ medianðDI1 ; :::; DIn Þ

PðDI < DI0 Þ ¼

#ðDIi < DI0 Þ i ¼ 1; 2; …; n nþ1

interest. For a random variable X with the cumulative distribution function F(x) in the univariate setting (e.g., drought duration D, severity S or spatial extent E), the return period T of the variable exceeding a certain threshold xp can be defined as (Singh et al., 2007):

m

m

(25)

The drought severity and the probability of prediction in equations (24) and (25) can be used to display prediction results. 2.5. Drought risk analysis 2.5.1. Univariate frequency analysis The traditional risk analysis of drought properties in the univariate case is generally based on the concept of return period or return level by fitting a probability distribution to the variable of

m

   ¼ ¼ Tp ¼  P X > xp 1  P X  xp 1  F xp

(26)

where m is the average or expected inter-arrival time of drought events (for the case of annual maxima, m equals 1 year); p is a critical probability level or threshold; xp is the design return level (or quantile) corresponding to the given return period Tp (or probability p); and P (X  x) is the exceedance probability. The return period for the non-exceedance probability P (X  x) can be defined in a similar way. 2.5.2. Multivariate frequency analysis Since multiple drought properties (e.g., duration and severity) or variables (e.g., precipitation and soil moisture) may be mutually correlated, the multivariate frequency analysis (MFA) is also desired. Extensive efforts have been devoted to generalize the concept of return period from the univariate case to the multivariate case (Shiau, 2006; Genest et al., 2007; Salvadori et al., 2011, 2013; Hao and Singh, 2013). For these purposes, the joint distribution is generally constructed to obtain the joint and conditional return periods for multivariate risk analysis. In a bivariate or multivariate setting, different definitions of the €ler et al., return period exist in hydrology (Salvadori et al., 2011; Gra 2013), either based on the joint probability distribution function or Kendall distribution function. The first type of return period can be classified as the AND (both D and S are exceeded) and OR (either D or S are exceeded) joint return period. For a continuous random vector (D, S) with the joint distribution expressed with a copula C(u,v), the joint return period for both the AND and OR cases can be expressed as (Shiau, 2003, 2006):

TAND ¼

m PðD  d

and

S  sÞ

¼

m 1  FD ðdÞ  FS ðsÞ þ Cðu; vÞ (27)

(24)

In addition, a parametric or nonparametric probability distribution function can be fitted to the predicted ensemble members of drought indices for probabilistic drought forecasting and assessing prediction uncertainty. For example, the nonparametric kernel density estimator can be used for estimating the PDFs of the ensemble members (Hao et al., 2016c), which is employed in this study. In addition, the probability of drought indices lower than certain values can also be obtained to assess the risk of drought. For example, the probability of predicted drought index DI lower than certain threshold DI0 (e.g., 0.8) can be expressed with the Weibull plotting position formula as:

203

TOR ¼

m PðD  d

or

S  sÞ

¼

m 1  Cðu; vÞ

(28)

where U and V are the marginal distributions of random vector (D, S); where FD(d) and FS(s) are marginal distributions of duration and severity, respectively; C(u,v) is the joint probability of the drought duration and severity. In addition, the joint return period can also be derived based on the Kendall distribution function (the secondary return period), which can be denoted as Kendall case. The return period based on Kendall distribution K(p) can be expressed as (Genest et al., 2007; Salvadori et al., 2013):

Tk ¼

m 1  KðpÞ

¼

m 1  PðCðu; vÞ  pÞ

(29)

Notice that the joint return period and joint drought index are closely correlated (compare equations (13) and (14) with equations (28) and (29)), which are all based on either the joint probability or Kendall distribution. The conditional frequency analysis of drought duration and severity refers to the risk of D (or S), conditioned on a certain level of S (or D). Specifically, the conditional return period of drought duration given drought severity S exceeding a certain value (denoted as TDjSs) can be defined as (Shiau, 2003, 2006):

204

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209

TDjSs ¼

m ½1  FS ðsÞ½1  FD ðdÞ  FS ðsÞ þ Cðu; vÞ

(30)

In a broader sense, multivariate drought frequency analysis refers to a variety of cases for drought risk analysis, for which the variables in equations (27)e(30) may refer to different drought indices or types (SPI and SRI, or meteorological and hydrological  ska, 2014; Hao drought) (Wong et al., 2013; Tokarczyk and Szalin and Singh, 2015). The drought condition may be related to the simultaneous deficit (or anomaly) in different variables, such as precipitation, temperature, soil moisture and streamflow, leading to different types of concurrent or integrated drought (or concurrent extremes). As such, the joint modeling of multiple variables or indices (e.g., SPI, PDSI, or temperature) would be important for the characterization of drought risk. For example, the occurrence of hydrological drought (denoted as V) conditioned on the meteorological drought (denoted as U) may be of particular interest for drought analysis. In this case, the conditional probability of V given U ¼ u0 can be used for the analysis, which can be expressed with a copula C as (Zhang and Singh, 2007):

PðV  vjU ¼ u0 Þ ¼

 vCðU; VÞ U ¼ u0 vU 

(31)

Based on equation (31), the conditional probability and return period can also be derived (Yue and Rasmussen, 2002; Zhang and Singh, 2007), which can be employed to assess the conditional risk of a specific type of drought given the condition of other types of drought or extremes.

3. Application The monthly precipitation, soil moisture and runoff data from 1932 to 2011 in Climate Division 3 in Texas, USA (obtained from the Climate Prediction Center (CPC)) are used to show how this drought package can be applied for drought monitoring, prediction and analysis. These data are used to compute drought indices, including the Standardized Precipitation Index (SPI), Standardized Soil Moisture Index (SSI), and Standardized Runoff Index (SRI) of 6month time scale (McKee et al., 1993; Shukla and Wood, 2008; Hao et al., 2014), representing the meteorological drought, agricultural drought and hydrological drought, respectively. The multivariate drought indices, JDSIJ and JDSIK are used to combine drought information from monthly precipitation (or SPI) and runoff (or SRI) for the integrated drought monitoring. These drought indices can then be used for the drought prediction based on the ESP concept and for drought frequency analysis. The framework of the package is shown in Fig. 1. In R, the “drought” package can be obtained from the R-Forge website and installed as:

Fig. 1. The framework of the package for drought monitoring, prediction and analysis. FD(d) and FS(s) are marginal distributions of duration and severity, respectively; F(d,s,…) is the distribution of drought properties; G(SDI1, SDI2, …) is the joint distribution of drought indices.

SDIðX; ts; distÞ where X is the monthly hydroloclimatic variables such as precipitation, soil moisture, runoff, groundwater, or snow; ts is the time scale of the drought index; dist is a character string to specify the empirical or parametric distribution function used to compute the probability. In the current version, the option of dist can be the Gringorten or Weibull plotting position formula (“EmpGrin”, “EmpWeib”) or the parametric Gamma and Weibull distributions (“Gamma” and “Lognormal”). The 6-month SPI, SSI and SRI for the period 1999e2012 were shown in Fig. 2 (a). It can be seen that severe droughts occurred during the period 1999e2001, 2005e2007 and 2010e2011, during which the drought indices were below the threshold value 0.8 (corresponding to D1 drought from the U.S. Drought Monitor (Svoboda et al., 2002)). For example, a severe drought event stroke Texas in 2011, causing huge losses to various sectors, including agriculture (Nielsen-Gammon, 2012). This drought event was indicated by these three drought indices. For example, for June 2011, the values of SPI, SSI and SRI were as low as 1.58, 0.82 and1.06, respectively. The SSI and SRI captured the onset of these severe drought events with lags of a few months to the SPI. The main reason was that drought generally originated from the precipitation deficit or meteorological drought. For example, in March 2011, the severe drought was captured by SPI (1.31), while no drought had been shown from SSI (0.22) and SRI (0.39).

install:packagesð“drought”; repos ¼ “http : ==R­Forge:R­project:org”Þ

3.1. Drought monitoring Three drought indices, SPI, SSI, and SRI, were computed, based on the period from 1932 to 2011 to monitor the meteorological, agricultural and hydrological drought, respectively. These standardized drought indices can be easily computed with the function in the R package as follows:

The multivariate drought index can also be derived with the package. The JDSI is used as the example to illustrate an application of the multivariate drought index to monitor the drought condition from multiple variables or indices. From this package, the JDSI in equations (13) and (14) can be computed as:

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209

focus on implementing this method to illustrate the application for drought predictions with the drought event in 2006 in Texas as a case study. The statistical prediction with the ESP method for the univariate and multivariate drought index can be achieved with the package though the function ESPPred by:

Drought Index

2 (a) 1 0 -1 -2

SPI SSI SRI -3 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 4 Drought Index

(b) 2 0

-4

JDSIJ

JDSIK

1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011

Fig. 2. (a) Monitoring drought condition with drought indices SPI, SSI and SRI; (b) Monitoring drought condition with multivariate drought indices JDSIJ and JDSIK from monthly precipitation and runoff.

JDSIðX; Y; ts; typeÞ where X and Y are two hydroclimatic variables; ts is the time scale of the drought index; type is used to specify the method used to define the JDSI (i.e., type ¼ 1 and 2 represents the use of joint probability and Kendall distribution in equations (13) and (14), respectively). In the current version, only the multivariate drought index based on the empirical methods introduced in equations (9) and (12) are used. Different types of parametric copulas, which can be obtained from the R package “copula” (Yan, 2007), can also be used to derive JDSI. For the illustration purpose, the JDSIJ and JDSIK, based on the joint probability and the Kendall distribution function of the precipitation and runoff are shown in Fig. 2(b), representing the combined meteorological and hydrological drought. Similar to the drought condition from the univariate drought index, the multivariate drought index also captured severe drought events during 1999e2001, 2005e2007, and 2010e2011. However, drought severity, along with drought onset and recovery, differs from the univariate drought index. For example, for the 2011 Texas drought, the JDSIJ for June was around 1.80, which was lower than the univariate drought index, indicating more severe drought. This is understandable, since JDSIJ is defined as the inverse of the probability of both precipitation and runoff being lower than certain values. Thus, the combined drought condition would be more severe than that indicated by the individual indicator or variable (i.e., P(X  x,Y  y)P(X  x)). From Fig. 2(b), the severity of JDSIJ is generally lower than that from JDSIK, both of which show the historical drought events. A detailed comparison of different drought indices is beyond the scope of this study and will be carried out in the future.

ESPPred ðX; Y; L; m; tsÞ where X and Y are the vector of monthly precipitation and runoff; L is the lead time; m is the end month for the prediction; ts is the time scale. The outputs of this function include the prediction of SPI, SRI, and JSDIJ. The predictive PDF of ensemble members of the prediction can be estimated with parametric distributions or an empirical method to facilitate the probabilistic drought prediction (Hao et al., 2016c). For illustration purposes, the predictive PDF of SPI for September 2006 of 1e5 months lead time is estimated with the kernel density estimation and shown in Fig. 3. From Fig. 3, the predicted SPI (1.25) for one month lead corresponding to the highest density is generally close to the observed SPI (1.41), indicating the 1-month lead prediction of the SPI is close to the observed drought severity. The main reason of the satisfactory performance is that for 1month lead prediction based on SPI6, the accumulated precipitation of 5 months from observations was used in computing the predicted SPI6 of the target month. The predictive probability of the SPI lower than the threshold 0.8 was obtained as 0.88, which indicates high likelihood of drought occurrence in September 2006. With the increase of lead time, the variability of the predicted ensemble members is wider than that for the 1-month prediction, which indicates the degradation of drought prediction. For multivariate drought prediction, the prediction of JDSIJ, based on the precipitation and runoff, is used as an example to illustrate the application of the prediction method. Similar to the univariate case, the probabilistic prediction of JDSIJ can be displayed through the PDF. For example, the PDF of predicted JDSIJ for 1e3 month lead prediction for September 2006 is shown in Fig. 4(a). With the increase of lead time, the predictive PDF tends to become flat, implying a larger variability of prediction. For illustrative purposes, the 1e2 month lead prediction of JDSIJ for the period from June to November 2006 is shown in Fig. 4(b). Generally, the median of prediction is close to that of observation (e.g.,

1.4 1-mon. 2-mon. 3-mon. 4-mon. 5-mon.

1.2

1

lead lead lead lead lead

0.8 PDF

3

-2

205

0.6

0.4

3.2. Drought prediction 0.2

The drought prediction with the ESP method introduced in section 2.4 is employed here for the drought prediction based on the univariate drought index SPI (Lyon et al., 2012) and multivariate drought index JDSIJ (Hao et al., 2014, 2016c) based on monthly precipitation and runoff in Climate Division 3 in Texas. Since the prediction performance has been well assessed in the literature, we

0 -4

-3

-2

-1

0

1

2

3

SPI Fig. 3. The PDF of predicted SPI of different lead time based on the ESP method (The observed value is denoted as star).

206

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209

3 1-mon. lead 2-mon. lead 3-mon. lead

(a) PDF

2 1 0 -3

-2.5

-2

-1.5

-1

-0.5

0

0.5

JDSIJ 0

(b)

1-mon. lead 2-mon. lead Observation

JDSIJ

-0.5 -1 -1.5 -2 Jun.

Jul.

Aug.

Sep.

Oct.

Nov.

Fig. 4. The PDF of predicted JDSIJ based on precipitation (or SPI) and runoff (or SRI) for 1e3 month lead prediction of for the period September 2006 (a) and 1e2 month lead prediction of JDSIJ from June to November 2006 (b).

November), though large discrepancies exist for the prediction of certain periods (e.g., July). For example, the median of JDSIJ of 1 and 2 month lead for November 2006 was 1.36 and 1.24, which is close to the observed JDSIJ -1.45. 3.3. Drought analysis The univarite and multivariate frequency analysis of drought properties can be achieved with the function DSFreq as:

of drought duration and severity. The inference function for marginal (IFM) method was used to estimate the copula parameter, in which parameter estimation of the marginal distribution can be split from that of the joint distribution (Joe, 1997). Three copulas from the Archimedian family, namely Clayton, Frank, and Gumbel, were selected to construct the joint distribution for comparison. Based on the AIC, the Gumbel copula was finally selected to model the joint distribution of drought duration and severity. The joint return period of drought duration and severity exceeding certain values was computed from equation (27), as shown in Fig. 6(a), along with the observed pairs of duration and severity. Based on Fig. 6, the joint return period from the copula method for the 2005e2006 drought was around 20 years. The conditional return period of the drought duration conditional on the severity (s  1, 2 and 3) was computed from equation (30), as shown in Fig. 6(b). For example, the conditional return period of the drought duration D  12 months conditioned on the drought duration S  3 was around 34 years. As stated before, the multivariate drought frequency analysis can also refer to the concurrence (or the combination) of different drought types. In this section, the risk of joint deficit of precipitation and runoff, representing the concurrent meteorological and hydrological droughts, was analyzed. The minimum value of SPI and SRI for each year was extracted and the joint distribution was constructed, based on copula (the Gumbel copula was selected based on AIC). The conditional probability of hydrological drought conditioned on meteorological drought can then be evaluated based on equation (31). Specifically, for the Gumbel copula, the conditional probability of SRI given SPI can be expressed explicitly as (Venter, 2002; Zhang and Singh, 2007; Aas et al., 2009):

PðV  vjU ¼ u0 Þ ¼ Cg ðu0 ; vÞu0 1 ð  log u0 Þq1

h

(32)

DSFreqðX; Y; ELÞ where X and Y can be the drought properties such as duration, severity and spatial extent (or drought indices); EL is the expected inter-arrival time. The outputs of this function include the univariate and multivariate return period defined in equation (26)e(29). 3.3.1. Univariate frequency analysis For the frequency analysis in this section, the drought severity during a drought event with duration of n months for a specific drought indicator (using SRI as an example) is defined as Si¼1, … n(0.8-SRIi), where 0.8 is the threshold to define drought. In this section, the SRI of 6 month time scale was used for univariate frequency analysis, from which drought duration and severity were extracted (with threshold 0.8). For example, for drought during 2005e2006, a drought duration of 14 months was obtained. To model the marginal distributions of drought duration and severity, four distributions, including exponential, gamma, Weibull, and lognormal, were selected as the candidate distributions. The exponential and gamma distribution were selected to model drought duration and severity, respectively, based on the Akaike information criterion (AIC) (Akaike, 1974). The cumulative probability of drought duration and severity from the theoretical distribution is shown in Fig. 5(a, b), which is compared with that from the Gringorten plotting position formula. Generally, the theoretical distribution was satisfactory in modeling drought duration and severity. The return periods of drought duration and severity were then computed, as shown in Fig. 5(c, d). For example, for drought duration and severity during the 2005e2006 drought, the return periods were around 18 and 17 years, respectively. 3.3.2. Multivariate frequency analysis The copula method was used to construct the joint distribution

where U and V are the marginal probability of SPI and SRI, respectively. Based on equation (32), the conditional probability distribution of hydrological drought given the meteorological drought can be estimated, as shown in Fig. 7(a). For example, given the meteorological drought index SPI ¼ 0.8, the conditional probability of the SRI lower than 0.5 was relatively high. The conditional return period of hydrological drought given the meteorological drought (e.g., SPI ¼ 0.8) was derived accordingly as the reciprocal of the probability in equation (31), as shown in Fig. 7(b). For example, the conditional return period for SRI equal to or lower than 0.8 was around 2.5 years. The conditional probability and return period of different drought properties and types provide useful information for drought management.

4. Conclusions and future developments An integrated R package “drought” is introduced for drought monitoring, prediction and analysis, which is available from the website: http://r-forge.r-project.org/projects/drought/. It allows users to compute univariate and multivariate drought indices, perform statistical drought prediction and also conduct risk analysis of drought. Both the theoretical background and practical applications of each component are presented to illustrate the implementation of this package. For the monitoring component, the computation of standardized drought index (SDI) is introduced, based on a similar manner of the SPI computation. In addition, the development of the Joint Drought Severity Index (JDSI) is also introduced, based on the joint probability and the Kendall distribution function of multiple hydroclimatic variables, which provide a general framework for characterizing the concurrent drought (or

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209

1

0.6 0.4 0.2 0

0

10 20 30 Duration (month)

0.4

0

40

0

5

10 15 20 Severity

25

100 Return Period (year)

Return Period (year)

0.6

0.2

100 80 (c) 60 40 20 0

(b)

0.8 Probability

0.8 Probability

1

(a)

207

0

10 20 Duration (months)

30

80 (d) 60 40 20 0

0

5

10 Severity

15

20

Fig. 5. Comparisons of fitted distributions of drought duration and severity with the empirical probability estimated from Gringorten plotting position formula (a, b) and the corresponding univariate return periods (c,d).

Fig. 6. Joint return periods of drought duration and severity (a) and conditional return period of drought duration given severity (b).

other concurrent extremes). The statistical drought prediction of univariate and multivariate drought indicators based on the ESP concept is also presented, which provides probabilistic drought predictions for early warning. At last, risk analysis of drought properties and types is also illustrated for univariate and multivariate (including conditional) frequency analysis. Due to the limitation of space, we only cover basic properties of drought modeling and assessment in this paper. There are potential extensions for the future development of this package. Currently, only commonly used drought indices are provided and other indices will be updated in the future for comparison purposes. In addition, the characterization of multiple drought indices for developing drought indices and risk analysis is mainly introduced

Fig. 7. Conditional probability (a) and conditional return period (b) of hydrological drought SRI given meteorological drought (SPI ¼ 0.8).

in the bivariate case. The extension to high dimensions will be added to the package based on multivariate analysis methods. For example, due to its advantage in flexible dependence modeling, the vine copula can be used for developing integrated drought indices to integrate multiple drought indices (Hao and Singh, 2015) and for statistical drought prediction. Moreover, other drought prediction methods and analysis tools (e.g., trend analysis and breakpoint analysis) will be embedded in this package in the future to aid drought analysis. The current package mainly provides the drought modeling and assessment at the local or grid scale. It can be extended straightforward to couple with certain data processing packages such “netcdf4” (Pierce, 2015) for the drought management at regional scales. For example, the Phase 2 of the North American Land Data

208

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209

Assimilation System (NLDAS-2) (Ek et al., 2011) provides land surface datasets in the Continental United States (CONUS) from 1979 to present with spatial resolution of 0.125 , which are available in the netCDF (network Common Data Form) format. The package can then be used for the drought monitoring, prediction and analysis in the CONUS with the NLDAS-2 products after suitable data processing (Xia et al., 2014; Hao et al., 2016a). It is expected that the package will help with drought practitioners in the related drought modeling for research and operational purposes. Acknowledgments This work is supported by National Natural Science Foundation of China (NSFC, No. 41601014) and Youth Scholars Program of Beijing Normal University (Grant No. 2015NT02). The R package copula is used for the copula modeling in this study. The authors are grateful to the Editors and reviewers for the valuable comments and suggestions. References Aas, K., Czado, C., Frigessi, A., et al., 2009. Pair-copula constructions of multiple dependence. Insur. Math. Econ. 44 (2), 182e198. Akaike, H., 1974. A new look at the statistical model identification. IEEE Trans. Autom. Contr. 19 (6), 716e723. Bachmair, S., Stahl, K., Collins, K., et al., 2016. Drought indicators revisited: the need for a wider consideration of environment and society. Wiley Interdiscip. Rev. Water 3 (4), 516e536. Beersma, J.J., Buishand, T.A., 2004. Joint probability of precipitation and discharge deficits in The Netherlands. Water Resour. Res. 40 (12), W12508. Behrangi, A., Nguyen, H., Granger, S., 2015. Probabilistic seasonal prediction of meteorological drought using the bootstrap and multivariate information. J. Appl. Meteor. Climatol. 54 (7), 1510e1522. Bloomfield, J., Marchant, B., 2013. Analysis of groundwater drought building on the standardised precipitation index approach. Hydrol. Earth Syst. Sci. 17, 4769e4787. Brown, J.F., Wardlow, B.D., Tadesse, T., et al., 2008. The Vegetation Drought Response Index (VegDRI): a new integrated approach for monitoring drought stress in vegetation. Gisci. Remote Sens. 45 (1), 16e46. Cheng, L., Hoerling, M., AghaKouchak, A., et al., 2016. How has human-induced climate change affected California drought risk? J. Clim. 29 (1), 111e120. Day, G.N., 1985. Extended streamflow forecasting using NWSRFS. J. Water Res. Plan. Man. 111 (2), 157e170. Dracup, J., Lee, K., Paulson Jr., E., 1980. On the definition of droughts. Water Resour. Res. 16 (2), 297e302. Dutra, E., Magnusson, L., Wetterhall, F., et al., 2013. The 2010e2011 drought in the Horn of Africa in ECMWF reanalysis and seasonal forecast products. Int. J. Climatol. 33 (7), 1720e1729. Ek, M.B., Xia, Y., Wood, E., et al., 2011. North American Land Data Assimilation System Phase 2 (NLDAS-2): development and applications. GEWEX News 21 (2), 6e7. Fuglem, M., Parr, G., Jordaan, I., 2013. Plotting positions for fitting distributions and extreme value analysis. Can. J. Civ. Eng. 40 (2), 130e139. Genest, C., Favre, A.-C., 2007. Everything you always wanted to know about copula modeling but were afraid to ask. J. Hydrol. Eng. 12 (4), 347e368. liveau, J., et al., 2007. Metaelliptical copulas and their use in Genest, C., Favre, A.C., Be frequency analysis of multivariate hydrological data. Water Resour. Res. 43, W09401 doi:09410.01029/02006WR005275. lez, J., Valde s, J., 2003. Bivariate drought recurrence analysis using tree ring Gonza reconstructions. J. Hydrol. Eng. 8 (5), 247e258. Goodrich, G.B., Ellis, A.W., 2006. Climatological drought in Arizona: an analysis of indicators for guiding the governor's drought task force. Prof. Geogr. 58 (4), 460e469. €ler, B., van den Berg, M., Vandenberghe, S., et al., 2013. Multivariate return peGra riods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation. Hydrol. Earth Syst. Sci. 17 (4), 1281e1296. Gringorten, I.I., 1963. A plotting rule for extreme probability paper. J. Geophys. Res. 68 (3), 813e814. Hannaford, J., Acreman, M., Stahl, K., et al., 2015. Enhancing drought monitoring and early warning by linking indicators to impacts. In: Joaquin Andreu, A.S., Paredes-Arquiola, Javier, Haro-Monteagudo, David, van Lanen, Henny (Eds.), Drought: Research and Science-policy Interfacing. CRC Press, pp. 287e292. Hao, Z., AghaKouchak, A., 2013. Multivariate Standardized Drought Index: a parametric approach for drought analysis. Adv. Water Resour. 57, 12e18. Hao, Z., AghaKouchak, A., Nakhjiri, N., et al., 2014. Global integrated drought monitoring and prediction system. Sci. Data 1, 140001. Hao, Z., AghaKouchak, A., Phillips, T.J., 2013. Changes in concurrent monthly precipitation and temperature extremes. Environ. Res. Lett. 8 (3), 034014.

Hao, Z., Hao, F., Singh, V.P., 2016c. A general framework for the multivariate multiindex drought prediction based on multivariate Ensemble Streamflow Predictions (ESP). J. Hydrol. 539, 1e10. Hao, Z., Hao, F., Singh, V.P., et al., 2016b. A theoretical drought classification method for the multivariate drought index based on distribution properties of standardized drought indices. Adv. Water Resour. 92, 240e247. Hao, Z., Hao, F., Xia, Y., et al., 2016a. A statistical method for categorical drought prediction based on NLDAS-2. J. Appl. Meteor. Climatol. 55 (4), 1049e1061. Hao, Z., Singh, V.P., 2013. Entropy-based method for bivariate drought analysis. J. Hydrol. Eng. 18 (7), 780e786. Hao, Z., Singh, V.P., 2015. Drought characterization from a multivariate perspective: a review. J. Hydrol. 527, 668e678. Hao, Z., Singh, V.P., 2016. Review of dependence modeling in hydrology and water resources. Prog. Phys. Geogr. 40 (4), 549e578. Heim, R.R., 2002. A review of twentieth-century drought indices used in the United States. Bull. Amer. Meteor. Soc. 83 (8), 1149e1166. Hoerling, M., Eischeid, J., Kumar, A., et al., 2014. Causes and predictability of the 2012 great plains drought. Bull. Amer. Meteor. Soc. 95 (2), 269e282. Huang, S., Huang, Q., Leng, G., et al., 2016. A hybrid index for characterizing drought based on a nonparametric kernel estimator. J. Appl. Meteor. Climatol. 55 (6), 1377e1389. Joe, H., 1997. Multivariate Models and Dependence Concepts. Chapman & Hall, London. Joe, H., 2014. Dependence Modeling with Copulas. CRC Press, Florida. Kao, S.C., Govindaraju, R.S., 2010. A copula-based joint deficit index for droughts. J. Hydrol. 380 (1e2), 121e134. Kendall, D., Dracup, J., 1992. On the generation of drought events using an alternating renewal-reward model. Stoch. Hydrol. Hydraul. 6 (1), 55e68. Keyantash, J.A., Dracup, J.A., 2004. An aggregate drought index: assessing drought severity based on fluctuations in the hydrologic cycle and surface water storage. Water Resour. Res. 40 (9), W09304. Kirtman, B., Min, D., Infanti, J., et al., 2014. The North american multimodel ensemble: Phase-1 seasonal-to-interannual prediction; Phase-2 toward developing intraseasonal prediction. Bull. Amer. Meteor. Soc. 95 (4), 585e601. Kurowicka, D., Joe, H., 2011. Dependence Modeling: Vine Copula Handbook. World Scientific, Singapore. Leonard, M., Westra, S., Phatak, A., et al., 2014. A compound event framework for understanding extreme impacts. Wiley Interdiscip. Rev. 5 (1), 113e128. Luo, L., Wood, E.F., 2007. Monitoring and predicting the 2007 US drought. Geophys. Res. Lett. 34 (22), L22702. Lyon, B., 2009. Southern Africa summer drought and heat waves: observations and coupled model behavior. J. Clim. 22 (22), 6033e6046. Lyon, B., Bell, M.A., Tippett, M.K., et al., 2012. Baseline probabilities for the seasonal prediction of meteorological drought. J. Clim. Appl. Meteor. 51 (7), 1222e1237. Ma, M., Ren, L., Singh, V.P., et al., 2014. New variants of the Palmer drought scheme capable of integrated utility. J. Hydrol. 519, 1108e1119. Makkonen, L., 2006. Plotting positions in extreme value analysis. J. Appl. Meteor. Climatol. 45 (2), 334e340. e, B., et al., 1992. The use of geometric and gammaMathier, L., Perreault, L., Bobe related distributions for frequency analysis of water deficit. Stoch. Hydrol. Hydraul. 6 (4), 239e254. McKee, T.B., Doesken, N.J., Kleist, J., 1993. The relationship of drought frequency and duration to time scales. In: Eighth Conference on Applied Climatology. Am. Meteorol. Soc., Anaheim, CA. , J., 2009. Multivariate archimedean copulas, d-monotone McNeil, A.J., Neslehova functions and [₁-norm symmetric distributions. Ann. Stat. 3059e3097. Mishra, A.K., Singh, V.P., 2010. A review of drought concepts. J. Hydrol. 391 (1e2), 202e216. Mishra, A.K., Singh, V.P., 2011. Drought modeling-A review. J. Hydrol. 403 (1e2), 157e175. Mo, K.C., Lettenmaier, D.P., 2014. Objective drought classification using multiple land surface models. J. Hydrometeorol. 15 (3), 990e1010. Mo, K.C., Lettenmaier, D.P., 2016. Precipitation deficit flash droughts over the United States. J. Hydrometeorol. 17 (4), 1169e1184. Mo, K.C., Lyon, B., 2015. Global meteorological drought prediction using the North American Multi-Model Ensemble. J. Hydrometeorol. 16 (3), 1409e1424. Nadarajah, S., 2009. A bivariate pareto model for drought. Stoch. Env. Res. Risk. A. 23 (6), 811e822. Nelsen, R.B., 2006. An Introduction to Copulas. Springer, New York. Nelsen, R.B., Quesada-Molina, J.J., Rodríguez-Lallena, J.A., et al., 2003. Kendall distribution functions. Stat. Probab. Lett. 65 (3), 263e268. Nielsen-Gammon, J.W., 2012. The 2011 Texas drought. Tex. Water J. 3 (1), 59e95. Nijssen, B., Shukla, S., Lin, C.-Y., et al., 2014. A prototype Global Drought Information System based on multiple land surface models. J. Hydrometeorol. 15, 1661e1676. ~ ez, J., Rivera, D., Oyarzún, R., et al., 2014. On the use of Standardized Drought Nún Indices under decadal climate variability: critical assessment and drought policy implications. J. Hydrol. 517, 458e470. € Ozger, M., Mishra, A.K., Singh, V.P., 2012. Long lead time drought forecasting using a wavelet and fuzzy logic combination model: a case study in Texas. J. Hydrometeorol. 13 (1), 284e297. Palmer, W., 1965. Meteorological Drought. Research Paper No. 45. U.S. Weather Bureau, Washington, D.C. Pan, M., Yuan, X., Wood, E.F., 2013. A probabilistic framework for assessing drought recovery. Geophys. Res. Lett. 40 (14), 3637e3642.

Z. Hao et al. / Environmental Modelling & Software 91 (2017) 199e209 Pierce, D., 2015. ncdf4: Interface to Unidata netCDF (Version 4 or Earlier) Format Data Files. R package. http://CRAN.R-project.org/package¼ncdf4. Pozzi, W., Sheffield, J., Stefanski, R., et al., 2013. Towards global drought early warning capability: expanding international cooperation for the development of a framework for global drought monitoring and forecasting. Bull. Amer. Meteor. Soc. 94 (6), 776e785. Quan, X.-W., Hoerling, M.P., Lyon, B., et al., 2012. Prospects for dynamical prediction of meteorological drought. J. Clim. Appl. Meteor. 51 (7), 1238e1252. R Development Core Team, 2015. R: a Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. http:// www.R-project.org/. Ramadas, M., Govindaraju, R., 2016. Choice of hydrologic variables for probabilistic drought classification. J. Irrig. Drain. Eng. 142 (3), 05015013. Saha, S., Moorthi, S., Wu, X., et al., 2014. The NCEP climate forecast system version 2. J. Clim. 27 (6), 2185e2208. Salas, J., Fu, C., Cancelliere, A., et al., 2005. Characterizing the severity and risk of drought in the Poudre River, Colorado. J. Water Res. Plan. Man. 131 (5), 383e393. Salvadori, G., De Michele, C., Durante, F., 2011. On the return period and design in a multivariate framework. Hydrol. Earth Syst. Sci. 15 (11), 3293e3305. Salvadori, G., Durante, F., Michele, C., 2013. Multivariate return period calculation via survival functions. Water Resour. Res. 49 (4), 2308e2311. Schubert, S.D., Stewart, R.E., Wang, H., et al., 2016. Global meteorological drought: a synthesis of current understanding with a focus on SST drivers of precipitation deficits. J. Clim. 29 (11), 3989e4019. Shafer, B., Dezman, L., 1982. Development of a Surface Water Supply Index (SWSI) to assess the severity of drought conditions in snowpack runoff areas. In: Proceedings of the Western Snow Conference. Sheffield, J., Goteti, G., Wen, F., et al., 2004. A simulated soil moisture based drought analysis for the United States. J. Geophys. Res-Atmos. 109, D24108. Shiau, J., 2003. Return period of bivariate distributed extreme hydrological events. Stoch. Env. Res. Risk. A. 17 (1), 42e57. Shiau, J., 2006. Fitting drought duration and severity with two-dimensional copulas. Water Resour. Manag. 20 (5), 795e815. Shukla, S., McNally, A., Husak, G., et al., 2014. A seasonal agricultural drought forecast system for food-insecure regions of East Africa. Hydrol. Earth Syst. Sci. 18 (10), 3907e3921. Shukla, S., Wood, A.W., 2008. Use of a standardized runoff index for characterizing hydrologic drought. Geophys. Res. Lett. 35 (2), L02405. Singh, V.P., Jain, S.K., Tyagi, A., 2007. Risk and Reliability Analysis: a Handbook for Civil and Environmental Engineers. ASCE Press, Virginia. Smith, A.B., Katz, R.W., 2013. US billion-dollar weather and climate disasters: data sources, trends, accuracy and biases. Nat. Hazards 67 (2), 387e410. Song, S., Singh, V.P., 2010. Meta-elliptical copulas for drought frequency analysis of periodic hydrologic data. Stoch. Env. Res. Risk. A. 24 (3), 425e444. Stagge, J.H., Tallaksen, L.M., Gudmundsson, L., et al., 2015. Candidate distributions for climatological drought indices (SPI and SPEI). Int. J. Climatol. 35 (13), 4027e4040. Staudinger, M., Stahl, K., Seibert, J., 2014. A drought index accounting for snow. Water Resour. Res. 50 (10), 7861e7872. Steinemann, A., 2003. Drought indicators and triggers: a stochastic approach to

209

evaluation. J. Am. Water Resour. Assoc. 39 (5), 1217e1233. Steinemann, A., Iacobellis, S.F., Cayan, D.R., 2015. Developing and evaluating drought indicators for decision-making. J. Hydrometeorol. 16 (4), 1793e1803. Svoboda, M., LeComte, D., Hayes, M., et al., 2002. The drought monitor. Bull. Amer. Meteor. Soc. 83 (8), 1181e1190.  ska, W., 2014. Combined analysis of precipitation and water Tokarczyk, T., Szalin deficit for drought hazard assessment. Hydrol. Sci. Journal/journal Des. Sci. Hydrol. 59 (9), 1675e1689. Tsakiris, G., Vangelis, H., 2005. Establishing a drought index incorporating evapotranspiration. Eur. Water. 9 (10), 3e11. Twedt, T.M., Schaake Jr., J.C., Peck, E.L., 1977. National Weather Service extended streamflow prediction [USA]. In: Proceedings Western Snow Conference (Albuquerque, NM). Venter, G.G., 2002. Tails of copulas. In: Proceedings of the Casualty Actuarial Society. pez-Moreno, J.I., 2010. A multiscalar drought Vicente-Serrano, S.M., Beguería, S., Lo index sensitive to global warming: the standardized precipitation evapotranspiration index. J. Clim. 23 (7), 1696e1718. Wilhite, D.A., Glantz, M.H., 1985. Understanding: the drought phenomenon: the role of definitions. Water Int. 10 (3), 111e120. Wong, G., Van Lanen, H., Torfs, P., 2013. Probabilistic analysis of hydrological drought characteristics using meteorological drought. Hydrol. Sci. J. 58 (2), 253e270. Wood, A.W., Lettenmaier, D.P., 2006. A test bed for new seasonal hydrologic forecasting approaches in the western United States. Bull. Amer. Meteor. Soc. 87 (12), 1699e1712. Wood, E.F., Schubert, S.D., Wood, A.W., et al., 2015. Prospects for advancing drought understanding, monitoring and prediction. J. Hydrometeorol. 16 (4), 1636e1657. Xia, Y., Ek, M.B., Peters-Lidard, C.D., et al., 2014. Application of USDM statistics in NLDAS-2: optimal blended NLDAS drought index over the continental United States. J. Geophys. Res-Atmos. 119 (6), 2947e2965. Xu, K., Yang, D., Xu, X., et al., 2015. Copula based drought frequency analysis considering the spatio-temporal variability in Southwest China. J. Hydrol. 527, 630e640. Yan, J., 2007. Enjoy the joy of copulas: with a package copula. J. Stat. Softw. 21 (i04), 251e256. Yevjevich, V., 1967. Objective Approach to Definitions and Investigations of Continental Hydrologi Droughts. Hydrology Paper 23. Colorado State U, Fort Collins. Yoon, J.-H., Mo, K., Wood, E.F., 2012. Dynamic-model-based seasonal prediction of meteorological drought over the contiguous United States. J. Hydrometeorol. 13 (2), 463e482. Yuan, X., Wood, E.F., 2013. Multimodel seasonal forecasting of global drought onset. Geophys. Res. Lett. 40 (18), 4900e4905. Yuan, X., Wood, E.F., Chaney, N.W., et al., 2013. Probabilistic seasonal forecasting of African drought by dynamical models. J. Hydrometeorol. 14 (6). Yue, S., Rasmussen, P., 2002. Bivariate frequency analysis: discussion of some useful concepts in hydrological application. Hydrol. Process. 16 (14), 2881e2898. Zargar, A., Sadiq, R., Naser, B., et al., 2011. A review of drought indices. Environ. Rev. 19, 333e349. Zhang, L., Singh, V.P., 2007. Bivariate rainfall frequency distributions using Archimedean copulas. J. Hydrol. 332 (1), 93e109.