Using GRACE Satellite Gravimetry for Assessing Large-Scale Hydrologic Extremes

Global assessment of the spatiotemporal variability in terrestrial total water storage anomalies (TWSA) in response to hydrologic extremes is critical for water resources management. Using TWSA derived from the gravity recovery and climate experiment (GRACE) satellites, this study systematically assessed the skill of the TWSA-climatology (TC) approach and breakpoint (BP) detection method for identifying large-scale hydrologic extremes. The TC approach calculates standardized anomalies by using the mean and standard deviation of the GRACE TWSA corresponding to each month. In the BP detection method, the empirical mode decomposition (EMD) is first applied to identify the mean return period of TWSA extremes, and then a statistical procedure is used to identify the actual occurrence times of abrupt changes (i.e., BPs) in TWSA. Both detection methods were demonstrated on basin-averaged TWSA time series for the world’s 35 largest river basins. A nonlinear event coincidence analysis measure was applied to cross-examine abrupt changes detected by these methods with those detected by the Standardized Precipitation Index (SPI). Results show that our EMD-assisted BP procedure is a promising tool for identifying hydrologic extremes using GRACE TWSA data. Abrupt changes detected by the BP method coincide well with those of the SPI anomalies and with documented hydrologic extreme events. Event timings obtained by the TC method were ambiguous for a number of river basins studied, probably because the GRACE data length is too short to derive long-term climatology at this time. The BP approach demonstrates a robust wet-dry anomaly detection capability, which will be important for applications with the upcoming GRACE Follow-On mission.


Introduction
Assessment of hydrologic extremes is valuable for water resources management, hazard preparedness, and food security [1,2].Identifying occurrences of hydrologic extremes has become critically important in light of a significantly larger number of recent extreme events [3][4][5] and the expected more widespread emergence of such events in the future [1,6].Among all possible indicators for quantifying hydrologic extremes [7][8][9], the terrestrial total water storage (TWS), which is a holistic and unique measure of the status of the Earth's water resources, has attracted significant attention in the recent decade [10][11][12][13][14][15].TWS includes snow, surface water, soil moisture, and groundwater stored at all levels, thus representing the upper bound of all usable water resources.Spatiotemporal variations of a region's TWS in response to hydrologic extremes may partly reflect the vulnerability or resilience of that region, a source of information that can be especially valuable for understanding the direct impact of hydrologic extremes at the regional scale.Systematic use of TWS as an indicator for extreme events has been rare in the past, mainly because of lack of reliable in situ observations [13,16,17].Since its launch in 2002, the gravity recovery and climate experiment (GRACE) satellite mission has enabled remote sensing of TWS anomalies (TWSA) at regional to continental scales, and routinely provided the much-needed global tracking of water storage changes until late 2017 [18,19].Its successor, the GRACE Follow-On mission, is scheduled to be launched in early 2018.
As the temporal length of GRACE records increases, a number of recent studies have looked into the potential merit of GRACE satellite gravimetry for drought monitoring [18,21,[35][36][37][38][39][40][41] and long lead-time flood prediction [11,12,42].A potential strength of TWS-based drought indicators is that they inherently account for storage capacity, including deep soil moisture and groundwater, thereby encompassing the impact of water demands and potentially providing more holistic and realistic diagnoses of drought conditions than commonly used indices (e.g., standardized precipitation index (SPI) and palmer drought severity index (PDSI)) at regional scales [18,35].Built on the same premises, Reager et al. [12] proposed to use the GRACE TWSA as an indicator of saturation excess, thus helping to assess the predisposition of a river basin to floods with 5-11 months of lead time.
Many existing studies applied the climatology derived from GRACE, or TWS-climatology (TC), to identify anomalous wet/dry events in TWS.Frappart et al. [43] used a standardized TWSA index to understand how interannual variability of rainfall impacted the land water storage in the Amazon River basin during 2003-2010.Thomas et al. [41] introduced a GRACE-based water storage deficit approach for hydrological drought characterization, in which water storage deficits were calculated as the negative residuals after removing monthly TC, and drought duration was calculated as the number of months with continuous deficits.Humphrey et al. [39] identified global drought events using the notion of water storage deficits similar to that of Thomas et al. [41].Kusche et al. [44] quantified statistical distributions of recurrence levels of annual high and low water storage from GRACE, with respect to TC.A potential limitation of the TC is related to the relatively short duration of GRACE data (<15 years), while long-term records (>30 years) are generally needed to reliably characterize climatology.Furthermore, TC-based approach may lead to a systematic bias if there is a trend (e.g., a drying or wetting signal) in the observations, an issue discussed in numerous papers on the notion of nonstationarity [45][46][47].
When used for TWS extreme detection, the climatology-based TC method is essentially a statistical Z-score method for anomaly detection.Like many other climate time series, TWS time series are nonlinear, nonstationary and tend to vary at multiple temporal scales, making filtering and detrending of TWS a nontrivial task.Detection of changes in nonstationary time series has been widely studied, and a large number of algorithms exist, such as the Mann-Kendall test [48,49], change point analysis [50,51], and numerous machine learning methods [52].In this study, we used an abrupt change detection algorithm, the breakpoint (BP) method, to detect points in time when abrupt changes occur.The underlying assumption is that hydrologic extremes are generally accompanied by abrupt changes in TWS time series; therefore, the onset, magnitudes, and duration of these abrupt changes may reflect the nature of extremes.
The particular BP detection algorithm used in this study originated from vegetation phenological change detection [53] and near real-time ecosystem disturbances [54].Its application to water storage dynamics is limited.Maeda et al. [40] applied the BP analysis to studying disruption of hydroecological equilibrium in the southwest Amazon River basin.They identified two abrupt changes (BPs) in TWSA that corresponded to known Amazon droughts in 2005 and 2010; however, when they performed BP analysis on monthly precipitation, they were not able to see any abrupt changes (in Section 3.3 we will show that the BP analysis should be performed on the cumulative precipitation residuals instead).
Several challenges may be perceived when using BP analysis to identify TWS extremes.First, the extreme wet events, typically of high intensity and short duration, are more identifiable than dry events, which have slow onset time and longer durations.Second, TWS extremes are often mixed with long-term and seasonal trends.So far, the predominant method used for detrending TWS time series in the GRACE research community is fitting and removing a trend model consisting of a long-term linear trend and seasonal cycles (modeled using Fourier harmonics) by using linear regression [22].Such a method implicitly assumes that the underlying process has regular seasonal cycles, which is not always true for all basins.Thus, the detrending process may either amplify or weaken the actual extreme event signals.Third, the choice of return period, which is a parameter required for many extreme event analysis procedures (including the BP detection algorithm used in this study), is not easy for short duration records.In this study, we used empirical mode decomposition (EMD) as an exploratory tool to help guide the BP detection process and tackle some of the aforementioned challenges.
EMD, a nonlinear, nonstationary time series decomposition method proposed by Huang et al. [55], decomposes a time series adaptively into a hierarchy of nonstationary oscillating modes, each containing information at a specific time scale.The resulting modes, called the intrinsic mode functions (IMF), separate multiple scales embedded in a nonstationary time series naturally without any a priori criterion [56], thus avoiding subjectivity in traditional harmonics-based decomposition.EMD has been used extensively in climate series analysis, financial event analysis, and transport planning [56][57][58][59].
A main finding of this study is that the EMD, through its scale separation process, can reveal the TWS extremes through one of the IMFs, thus making it a promising exploratory tool for assisting event extraction, interpretation, and BP detection.
The main objectives of this work were two-fold: (1) assessing the performance of TC-based and EMD-assisted BP detection methods for identifying abrupt changes in TWS, and evaluating the results both quantitatively and qualitatively using historical climate records, and (2) examining the synchronicity of observed TWSA with precipitation forcing (Unless otherwise specified, the BP method implies EMD-assisted BP method for the rest of this discussion).Specifically, the study was designed to address the following questions: (1) how results of BP and TC methods differ from each other; (2) how abrupt changes in TWSA are related to antecedent abrupt changes in precipitation anomalies; and (3) is there linkage among TWSA anomalies across different basins?To help answer the latter two questions, an event-based framework was adopted to quantify the coincidence rates of precipitation and TWSA extremes.So far, a systematic study is lacking on the performance of abrupt-change detection methods in extracting TWSA extremes.In anticipation of the GRACE Follow-On mission, such a global comparative study is necessary to gain insights about the general applicability and strength of both TC and BP methods, thereby supporting the meaningful use of TWSA as an indicator of hydrologic extremes.
This paper is organized as follows.Section 2 describes the GRACE satellite gravimetry and the precipitation data used.For satellite gravimetry, we used both the recently released mass concentration (mascon) solutions and gridded spherical harmonic (SH) products.For precipitation, we considered global land data assimilation system (GLDAS) Version 2.1 forcing data [60] and tropical rainfall measuring mission (TRMM) 3B43 V7 product.Section 3 presents formulations of our EMD-assisted BP detection method and the principles of event coincidence analysis.Section 4 shows results obtained for the world's 35 largest river basins.The EMD process was applied to identify temporal scales for BP analyses and event extraction.Events extracted using the TC and BP methods were then carefully examined and compared to those obtained using SPI.Section 5 provides a discussion on the results.The main findings and conclusions are summarized in the last section.Because of the nature of this work, the use of a large number of acronyms is inevitable.Thus, a list of acronyms is provided at the end of the text for reference.

GRACE Total Water Storage Data
Our study period spans from April 2002 through December 2015.Gridded 1 × 1 • monthly GRACE SH products (RL05) from The University of Texas Center for Space Research (CSR) and NASA Jet Propulsion Laboratory (JPL) were downloaded from the JPL Tellus site (http://grace.jpl.nasa.gov).A multiplicative scaling factor dataset distributed with the gridded products is commonly applied to restore signal attenuation caused by destriping and smoothing [61].However, the multiplicative scaling factors, which are derived from an LSM, are themselves uncertain due to uncertainties related to either unrepresented processes in the LSM (e.g., irrigation and lateral water redistribution) or missing storage compartments (e.g., deep groundwater storage) [62].
Recently, several GRACE mascon solutions have been released.In general, these mascon solutions implement noise reduction during processing steps, which is more optimal than the SH approach that only performs noise reduction and signal restoration as part of the post-processing.The regional mascon solutions allow the land-ocean boundaries to be delineated, greatly reducing leakage at this interface relative to the global SH processing.JPL-M uses near-global geophysical models, including GLDAS models, to constrain the solutions during processing [63].Instead of using external geophysical models, the CSR mascon (CSR-M) solution is derived using a regularization matrix that is temporally variable and represents geophysical properties for a particular data measurement interval while retaining the long-term geophysical characteristics of GRACE signals [64].
Recently, Scanlon et al. [65] compared seasonal signals and long-term trends derived using gridded SH products (CSR) and two mascon solutions (JPL-M and CSR-M) for 176 global basins.They showed that the performance of mascon solutions is equivalent or superior to SH solutions.Thus, in this study we mainly used CSR-M, which is distributed as a 0.5 × 0.5 • gridded product and can be downloaded from CSR's web site, http://www.csr.utexas.edu/grace.

Precipitation Data
Monthly TRMM datasets (TMPA/3B43 V7, 0.25 • ) were downloaded from NASA Mirador data portal (https://mirador.gsfc.nasa.gov).The TMPA/3B43 V7 product was generated by combining information from passive microwave and precipitation radar onboard TRMM satellite, the Visible and InfraRed Scanner onboard the Special Sensor Microwave Imager, and gauge observations wherever available [66].The spatial coverage of TRMM 3B43 V7 is 50 • S-50 • N and temporal coverage is from January 1998 to the present.The TRMM satellite mission ended in 2015, but the TMPA product will continue to be produced through early 2018 (https://pmm.nasa.gov/TRMM).
A number of high-latitude basins that we would like to study are outside the TRMM coverage area.Thus, we also extracted precipitation forcing data from the global land data assimilation system (GLDAS) V2.1 monthly product (1 × 1 • ), which covers the globe from 60 • S to 90 • N and temporally from January 2000 onward.GLDAS data were also retrieved from the NASA Mirador data portal.A major issue found with the previous version of GLDAS (V1) is that its precipitation forcing had a lower bias than the Global Precipitation Climatology Project (GPCP) dataset [67].In the recently released GLDAS V2.1, the GPCP 1 • daily dataset and an updated disaggregation routine were used to prepare GPCP fields.In Supporting Information Part A (SI. A), we compared monthly time series of GLDAS V2.1 precipitation data to TRMM data over selected river basins.We showed that the two products gave very similar results with a high Nash-Sutcliffe efficiency (≥0.95).On the other hand, the match of precipitation forcing between the previous version of GLDAS (V1) with TRMM is poor (see Figure S1).Thus, the GLDAS V2.1 precipitation forcing was selected for use in this study.
SPI time series was calculated from GLDAS V2.1 precipitation using the R package SPEI [68], with a time scale of three months (i.e., SPI3).The SPI values generally range from −3 (extreme dry) to +3 (extreme wet).We used an SPI threshold of 1 (−1) to identify wet (dry) events, which represent moderate events in terms of severity.SPI that is calculated using either TRMM or GLDAS V2.1 forcing may also be affected by the limited duration of these products.A previous study evaluated the performance of TMPA products in identifying meteorological droughts by using ten years of TRMM record (2000-2009) [69]; they showed that dry events could be identified spatially and temporally using the TMPA research precipitation datasets (3B42 V6 and 3B42 V7) and results were comparable with other non-TMPA, satellite-gauge combined precipitation products.

Methods
The terrestrial water budget equation is written as where P, ET, and Q denote precipitation, evapotranspiration, and surface runoff, respectively.In this work, only the coincidence between TWS and P events are quantified and compared.In the following, we first describe the EMD and BP algorithms and their coupling, and then describe a procedure for event coincidence analysis.

Empirical Mode Decomposition
EMD may be regarded as a sifting process that adaptively extracts zero-mean, oscillatory modes (i.e., IMFs) from a time series.The advantage of EMD, compared to the more commonly used wavelets, is that it is a data-driven algorithm that decomposes a signal in a natural way without prior knowledge about the signal of interest [56].The original EMD procedure is summarized in the following pseudo algorithm [55]: Set y 0 (t) = TWSA(t), and set i = 1, where i is the index of IMF to be extracted 2.
Extract the i-th Identify all local extrema in l k−1 (t) c.
Generate the upper and lower envelopes, e max (t) and e min (t), of l k−1 (t) by cubic spline interpolation d.
Calculate the mean value of upper and lower envelopes as f. Check whether l k (t) satisfies properties of IMF, namely, the number of extrema and the number of zero crossings must either equal or differ at most by one, and at any point the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero.If the properties are satisfied, define l k (t) as the i-th IMF c i (t); otherwise, set k = k + 1 and go to Step 2b.

Define the current residual y
Repeat Steps 2-3 until either the residual y i (t) becomes a monotonic function, or the number of zero crossings and extrema is the same as that of the successive sifting step.
After applying the EMD, the original time series is decomposed into the sum of IMFs and a residual term where N is the number of IMFs; and r N (t) denotes the final monotonic residual term from EMD, which is equivalent to the trend term.Some other interesting properties of the IMF are: they are nearly orthogonal (i.e., different IMFs show information from different frequencies) and normally distributed; the intervals between extrema are not equal, but the mean period of each IMF, defined as the mean interval between two adjacent maxima (or minima) in that IMF, is almost twice of the previous one [56].
We will use the notion of mean period extensively in the following discussion.
In practice, EMD may be affected by mode mixing, which refers to a situation where a single IMF consists of a mixture of oscillations of dramatically different temporal scales [56].Mode mixing may affect the uniqueness of decomposition using EMD.Wu and Huang [56] proposed a more robust, ensemble EMD method to alleviate the mode mixing issue and improve stability of EMD.In the ensemble EMD, a white noise series is added to the original data before applying EMD.The procedure is repeated for an ensemble of white noise realizations.The ensemble mean is then used as the final results.We used the MATLAB ensemble EMD toolbox [70].The number of realizations used was 300 and the standard deviation of white noise was set to 0.01.

Breaking Point Detection
We adopted an abrupt change detection algorithm, breaks for additive seasonal and trend (BFAST), originally developed to conduct phenological change detection [71].The starting point of BFAST is the following additive time series model [71], where y t = TWSA(t), T t and S t are trend and seasonal components, e t represents residual term not included in T t and S t , and n is the total number of months in study period.BFAST attempts to split the seasonal and trend components into multiple segments, and the endpoints of those segments (except for the two boundaries) are BPs.To do this, an iterative procedure is used to estimate T t and S t and their abrupt changes.During each iteration, the trend component is split into piecewise linear segments such that the trend model within each segment is in which a i and b i are the intercept and slope associated with the trend T t,i of the i-th segment, and m is the number of BPs.Similarly, the seasonal component of each segment is parameterized using a harmonic model where β j and φ j are the amplitude and phase of the seasonal model, ω j are angular frequencies, and T p,j are the corresponding periods.These steps are iterated until no further changes can be made in the number and positions of BPs.The ordinary least squares moving sum test is used to test whether one or more BPs occur at a given significance level (e.g., 0.05), and the Bayesian Information Criterion is used to determine the optimal number of BPs.More implementation details can be found in Verbesselt et al. [53].For the seasonal model, we assume it includes both annual and semiannual harmonic cycles, namely, k = 2, and T p = 1 year and 0.5 year in Equation (7).Unlike the conventional TWSA decomposition that assumes a single long-term linear trend, BFAST assumes piecewise linear trends for locating abrupt changes.However, BFAST still uses a harmonic seasonal model (Equation ( 7)) which may shadow the true effect of small-magnitude hydrologic extremes.
The main required user input of BFAST is the minimum interval (denoted as h) between two adjacent BPs.Using too large a value will miss some BPs, while using too small a value will give too many false positives.In this study, we first identified the most relevant IMF from EMD, and then used the corresponding mean period as a cue to set h.The subtle difference between an IMF's mean period and h is worthy of more explanation.By construction, EMD forces an IMF to have the same number of zero-crossings and extrema, but the signal amplitude is not constrained.In fact, the effect of extrema is usually amplified in one of the IMFs.Thus, by visually comparing all the IMFs with known historical extremes for a basin, we can identify this "extreme indictor" IMF easily.On the other hand, the minimum interval h is a user-set parameter to help BFAST to avoid missing those back-to-back extreme events (e.g., the Amazon 2009 flood followed by its 2010 extreme drought).If the value of h is greater than the event intervals, the effect of all events in between will be averaged (i.e., shadowed).In our EMD-assisted BP analysis, a reasonable choice for h is one half the indicator IMF's mean period (i.e., distance between adjacent maximum and minimum values) such that all significant extrema can be identified as BPs.
Alternatively, one may either perform a conventional frequency analysis of extreme events for each basin based solely on historical records, or run BFAST using different h values.Both are more labor intensive and even infeasible when the number of basins is large or when historical information is limited.Therefore, our EMD-assisted BP procedure provides a more rigorous and objective way for extreme event analysis.We emphasize that only partial knowledge of extremes at a basin is required.This is because all IMFs have distinct frequency patterns, and knowledge of the timing of two or more extremes is usually sufficient for choosing the right extreme indicator IMF.
Depending on the sign of fitted trends immediately before and after a BP, four types of BPs (abrupt changes) may be possible, namely, +/+, −/−, and +/−, and −/+.Specifically, we will use the following BP classification convention in the Results section: Type 1, negative break in a positive trend (positive trends before and after the break); Type 2, positive break in a negative trend (negative trends before and after the break); Type 3, increase to decrease trend reversal (+/− trends); and Type 4, decrease to increase trend reversal (−/+ trends).Physically speaking, the first two BP types may be considered abrupt changes that do not have a long-term impact on a basin (e.g., flash drought or flood), while the last two represent changes that tend to reverse the pattern of a basin's interannual dynamics.The R package, BFAST [71], was used.

Event Coincidence Analysis
Quantification of the concurrence of anomalous events in TWSA and SPI, which are independently derived, gives a means for cross-examining events detected by TC-and BP-based methods.After an event detection method is validated, the calculated concurrence rate between extreme TWS and precipitation events then helps to infer whether climate forcing is the dominant driver of TWS changes in a river basin.
In this work, an event coincidence analysis (ECA) method was applied to quantify potential causal relationships between SPI anomalies and TWSA extreme events.Unlike the linear correlation analysis, ECA is nonparametric and only requires knowledge of the timing of events.Given two binary event time series, ES A and ES B , which can be events of either the same or different types, ECA counts how often events in those two series co-occur, within a user-defined time delay τ and a sliding window ∆T.The time delay τ is used to account for possible fixed shifts between event series, while the sliding window ∆T is used to accommodate irregularity in event concurrence times.For the purpose of this study, we were mainly interested in a so-called precursor coincidence rate defined as [72]: where N X (X = A, B) is the number of events and t X are event occurrence times; H is the Heaviside step function, and I [0,∆T] (•) is an indicator function that is equal to 1 when the argument falls within the interval [0, ∆T] and 0 otherwise.Thus, if ES A represents TWSA event series and ES B represents SPI anomaly series, then r p measures the fraction of TWSA events that are preceded by at least one precipitation event.Similarly, if ES A and ES B represent TWSA event series extracted from two different basins, then r p is a measure of how TWSA events in the two basins may co-occur through some common climate forcing mechanism (e.g., El Niño Southern Oscillation or ENSO) that affect both basins.
In our analyses, the delay parameter τ was set to zero, while the width of sliding window in Equation ( 8) was determined by considering correlation strength between precipitation and TWSA.To do this, we decomposed each term on the right hand side of Equation ( 1) into constant, linear, and time-variable residual components, which leads to the following equation after rearranging terms [34]: where subscripts 0 and 1 indicate constant and linear terms, respectively, the asterisk (*) represents the residual term, C is an integration constant, and integration of the right-hand side equation is from the beginning of study period t s to t.After fitting and removing a quadratic trend in the form of a + bt + ct 2 from GRACE TWSA using linear regression, we get an estimate of the left hand side of Equation ( 9).Similarly, after removing a linear trend from monthly P time series and integrating the residuals, we can obtain an estimate of the contribution of P to the integral on the right hand side of Equation (9).In essence, the data processing technique described in Equation ( 9) allows a correlation analysis between cumulative P residuals (now a state variable) and TWSA.Using the technique, Crowley et al. [34] qualitatively compared GRACE-derived TWSA with cumulative P residuals for the Congo basin during 2002-2006 and found good agreements.
In this study, we applied the aforementioned data processing technique to quantify lagged correlation between cumulative P residuals and TWSA (ET and Q were not considered here).The lag at which the correlation is maximum, in turn, was referenced when setting the width of sliding window.For instance, if the correlation lag is 8 months, then the width of ∆T is set to the same length to capture the potential causal relationship between P and TWSA, at least in an average sense.
The original precursor coincidence rate given in Equation (8) does not differentiate between dry (negative anomaly) and wet (positive anomaly) events.Thus, for TC-based ECA, we introduced a weighted coincidence rate, r p , that combines event coincidence rates calculated separately for each type of event where superscripts + and − indicate wet and dry events, and N represents the number of each type of events.By definition, the coincidence rate ranges from 0 (no coincidence between two event series) to 1 (complete coincidence).In this work, we used the R package CoinCalc for ECA [72].

Results
The 35 largest (by area) river basins from the total runoff integrating pathway (TRIP) database [73] (see Figure 1) were selected to demonstrate the TC-based and BP detection methods.Climates of the selected basins, which are classified using the mean annual aridity index [74], are humid (H) (22 basins), semi-humid (SH) (2 basins), semi-arid (SA) (10 basins), and arid (A) (1 basin) (4th and 5th columns in Table 1).The percentage of irrigated area in each basin is listed in the last column of Table 1.
The 35 largest (by area) river basins from the total runoff integrating pathway (TRIP) database [73] (see Figure 1) were selected to demonstrate the TC-based and BP detection methods.Climates of the selected basins, which are classified using the mean annual aridity index [74], are humid (H) (22 basins), semi-humid (SH) (2 basins), semi-arid (SA) (10 basins), and arid (A) (1 basin) (4th and 5th columns in Table 1).The percentage of irrigated area in each basin is listed in the last column of Table 1.

EMD Results
EMD was first applied to the basin-averaged TWSA time series obtained using the CSR-M GRACE product.In Figure 2, results for the Amazon River basin are shown.The EMD analysis shows an upward linear trend over the study period (bottom panel).The frequencies of 5 IMFs decrease from the top to bottom.IMF1 is the dominant mode in signal amplitude, consisting of a sinusoid varying at annual cycle (mean period = 12 months).IMF2 is at approximately biennial scale (mean period = 28 months).IMF3 represents a longer-term, interannual oscillation (mean period = 56 months).IMF4 and IMF5 are similar, both representing small-amplitude, interdecadal oscillations; thus, they should be combined according to the guidance in [70] (p.12).In the following discussion, we put the relevant references after all identified events for ease of interpretation and confirmation.We also consulted the online international disasters database, EM-DATA [75], for flood and drought events.
The extrema of IMF2 correspond well to the known hydrologic extremes occurred at the Amazon River basin, including the extreme drought of 2005 [76] and 2010 [40,43], and the extreme floods in 2009 [77], austral summer in 2011-2012 [78], and austral summer of 2013-2014 [79].The less discussed wet year 2006 and dry year 2007 [43] are also discernable.Thus, in this case IMF2 is the "extreme indicator" IMF.Note that the index of extreme indicator IMF may vary for other basins, depending on the mean return period.
Amazon River basin is known for its strong cyclic patterns.So, can we find similar indicator IMFs for other basins?

EMD Results
EMD was first applied to the basin-averaged TWSA time series obtained using the CSR-M GRACE product.In Figure 2, results for the Amazon River basin are shown.The EMD analysis shows an upward linear trend over the study period (bottom panel).The frequencies of 5 IMFs decrease from the top to bottom.IMF1 is the dominant mode in signal amplitude, consisting of a sinusoid varying at annual cycle (mean period = 12 months).IMF2 is at approximately biennial scale (mean period = 28 months).IMF3 represents a longer-term, interannual oscillation (mean period = 56 months).IMF4 and IMF5 are similar, both representing small-amplitude, interdecadal oscillations; thus, they should be combined according to the guidance in [70] (p.12).In the following discussion, we put the relevant references after all identified events for ease of interpretation and confirmation.We also consulted the online international disasters database, EM-DATA [75], for flood and drought events.
The extrema of IMF2 correspond well to the known hydrologic extremes occurred at the Amazon River basin, including the extreme drought of 2005 [76] and 2010 [40,43], and the extreme floods in 2009 [77], austral summer in 2011-2012 [78], and austral summer of 2013-2014 [79].The less discussed wet year 2006 and dry year 2007 [43] are also discernable.Thus, in this case IMF2 is the "extreme indicator" IMF.Note that the index of extreme indicator IMF may vary for other basins, depending on the mean return period.
Amazon River basin is known for its strong cyclic patterns.So, can we find similar indicator IMFs for other basins?Figure 3a-d show the indicator IMFs identified for four additional large river basins, namely, Mississippi (H), Nile (SA), Zambezi (SH), and Murray Darling (SA).Different from the Amazon River basin, all basins show strong nonstationarity during the study period (solid TWSA lines in Figure 3).
In the Mississippi River basin, the indicator IMF (dashed line) identifies the extreme U.S. midwest flood in 2011 [12], the extreme drought in 2012 [41]; also, the dry event in 2002-2003, the hurricane-related wet year 2005, wet event in the beginning of 2008, wet year 2013-2014 can be observed [75,80].
For the Nile River basin, the indicator IMF reveals major dry events in 2003,2005,2006,2009,2011, some of which were related to ENSO; the wet events in 2003, 2005, 2007 were also documented [75,81].
Remote Sens. 2017, 9, x FOR PEER REVIEW 11 of 25 Figure 3a-d show the indicator IMFs identified for four additional large river basins, namely, Mississippi (H), Nile (SA), Zambezi (SH), and Murray Darling (SA).Different from the Amazon River basin, all basins show strong nonstationarity during the study period (solid TWSA lines in Figure 3).
In the Mississippi River basin, the indicator IMF (dashed line) identifies the extreme U.S. midwest flood in 2011 [12], the extreme drought in 2012 [41]; also, the dry event in 2002-2003, the hurricane-related wet year 2005, wet event in the beginning of 2008, wet year 2013-2014 can be observed [75,80].
For the Nile River basin, the indicator IMF reveals major dry events in 2003, 2005, 2006, 2009, 2011, some of which were related to ENSO; the wet events in 2003, 2005, 2007 were also documented [75,81].In the case of Zambezi River basin, which is known for its strong annual storage variations [82], the shape of indicator IMF is more interesting and is dominated by a long extreme drought occurred in eastern Africa in 2003-2006 (separated by a wet event in 2004), followed by a sudden increase in 2007 [83]; afterward, the region encountered flood events almost annually between 2007 and 2014 [75].The Murray Darling basin in Australia experienced a prolonged drought period [84] that ended in 2010; the recent wet events in 2012 and 2014-2015, and dry events in 2011, end of 2013, and 2015 may also be seen [75].In this case, IMF2 reveals all major TWSA extremes occurred at the basin during the study period.
In the case of Zambezi River basin, which is known for its strong annual storage variations [82], the shape of indicator IMF is more interesting and is dominated by a long extreme drought occurred in eastern Africa in 2003-2006 (separated by a wet event in 2004), followed by a sudden increase in 2007 [83]; afterward, the region encountered flood events almost annually between 2007 and 2014 [75].The Murray Darling basin in Australia experienced a prolonged drought period [84] that ended in 2010; the recent wet events in 2012 and 2014-2015, and dry events in 2011, end of 2013, and 2015 may also be seen [75].
The full set of IMFs for each basin are plotted in SI Part B, which suggest that in all basins except Mississippi a nonlinear EMD trend exists that appeared to be relatively flat between 2002 and 2009/2010 and then started to increase after that.
The mean periods for all basins are listed in Table 1, which generally range from two to three years.This information was then used to derive the h values required for the following BP analyses.
Remote Sens. 2017, 9, x FOR PEER REVIEW 12 of 25 The full set of IMFs for each basin are plotted in SI Part B, which suggest that in all basins except Mississippi a nonlinear EMD trend exists that appeared to be relatively flat between 2002 and 2009/2010 and then started to increase after that.
The mean periods for all basins are listed in Table 1, which generally range from two to three years.This information was then used to derive the h values required for the following BP analyses.

TC and BP Results
The BP and TC-based methods were performed using the same basin-averaged TWSA time series.Here we showcase the results for four major river basins, Amazon, Mississippi, Nile, and Zambezi.
Figure 4a-d show BPs detected in each of the four basins.On the left axis of each subplot, the BFAST trend (dark green solid line) is overlaid on the original TWSA data (dark gray solid line), together with BP times (vertical dotted line) and the associated 95% confidence intervals (shaded green boxes near the top of each subplot).For comparison, monthly SPI3 series is plotted on the right axis (filled line plots) and its ±1 thresholds are shown with horizontal, dotted lines in gray.

TC and BP Results
The BP and TC-based methods were performed using the same basin-averaged TWSA time series.Here we showcase the results for four major river basins, Amazon, Mississippi, Nile, and Zambezi.
Figure 4a-d show BPs detected in each of the four basins.On the left axis of each subplot, the BFAST trend (dark green solid line) is overlaid on the original TWSA data (dark gray solid line), together with BP times (vertical dotted line) and the associated 95% confidence intervals (shaded green boxes near the top of each subplot).For comparison, monthly SPI3 series is plotted on the right axis (filled line plots) and its ±1 thresholds are shown with horizontal, dotted lines in gray.The BP detection procedure seeks to optimally divide the whole time series into multiple short periods, separated by abrupt changes (BPs).In our study, we found that all BPs are related to abrupt changes in the trend.For the Amazon River basin, the major drought in 2005 was part of a declining trend that started in late 2004 and was only interrupted by the wet year in 2006.Similarly, the extreme drought in 2010 was part of a sharp declining trend that started since the end of extreme Amazon flood in 2009.Similar to the extreme indicator IMFs resulting from the EMD, BPs are related to local extrema in a time series-the number of BPs detected agrees with the number of extrema in the extreme indicator IMF.Unlike the indicator IMFs, however, BP detection reveals not only abrupt changes, but also the actual trends in between.The magnitudes of these short-term trends reflect the rate of decline or recovery.In this case, the transition from the extreme wet 2009 to extreme dry 2010 makes the trend in 2009-2010 the steepest decline.In many cases, the abrupt changes may also be accompanied by sharp jumps in the trend to adapt to a new short-term regime.For instance, the positive jump at the end of 2005 drought, the positive jump related to 2009 flood, and the negative jump related to the dry event in 2013.
Using the BP classification scheme mentioned in Section 3, we can see examples of all four types of abrupt changes in Figure 4a.Specifically, the 2003-2005 dry period was interrupted by a small-magnitude wet event in 2004, making it a Type 2 event preceded and succeeded by negative trends.The 2010-2014 wet period was interrupted by a dry event in 2012, making the BP a Type 1 event that was both preceded and succeeded by positive trends.The 2009 flood was a Type 3 event (+/−), while the 2010 extreme drought was a Type 4 event (−/+).
For the Mississippi River basin, notable BPs are related to a wet event in 2004, drought period in 2005-2006, extreme flood in 2011, and a severe drought that started in 2012.Like in the case of Amazon River basin during 2009-2010, the short period 2010-2012 included a massive flood followed by a drought event.The 2012 Central U.S. flash drought impacted major agricultural areas across the region and had a severity comparable to historical extreme droughts during the 1930s and 1950s [85,86].The 2012 BP is thus accompanied by a sharp drop in TWSA magnitude.A major difference between Figure 4b and its corresponding IMF (Figure 3a) is in the period 2008-2010, where the IMF showed a wet event in 2008 and two dry events in 2007 and 2009, while the BP results suggested a positive trend in the period.This phenomenon may be attributed to the mixing of signals by the harmonic seasonal component assumed in BFAST.
Both the Nile River basin and Zambezi River basin are characterized by high-frequency precipitation and TWSA variations.In Nile River basin, major BPs corresponded to 2003, 2007, 2010 wet events, and 2005, 2009, 2011 dry events.In Zambezi River basin, a severe drought occurred in 2005, followed by a 3-year wet period until 2008, and then a downward trend continued after 2010.In these cases, the BP algorithm detected the major events, but also missed some smaller-magnitude, short-term events, for example, the sharp downward spike in late 2005 in Nile River basin.We note that the indicator IMF for the Nile River basin (Figure 3b) detected the 2005 event.Similarly, the indicator IMF for the Zambezi River basin (Figure 3c) also identified the short-term events happened in the beginning of 2012 and in late 2003.
Figure 5a-d provide results obtained using the TC-based method.The normalized TWSA index values were smoothed using the 3-month moving average.Overall, the TC method seems to capture some SPI extreme events well, including their durations, for example, the 2005, 2010 droughts and 2009 flood in the Amazon River basin; the 2009 flood and 2012 midwest drought in the Mississippi River basin; and the beginning of extreme dry period in the Zambezi River basin in 2005.However, some event durations were significantly exaggerated, as compared to SPI.For example, the Mississippi 2005-2009 drought period and the 2005-2008 Zambezi drought.
It may be argued that TWS takes longer to recover to its normal conditions than SPI [41]; however, differences between TWSA extremes and rainfall anomalies need to be clarified and quantified, at least in terms of the event timings.To address this latter issue, we applied the ECA method to quantify concurrence between SPI anomalies and TWSA by considering the inherent lags identified during the following correlation analysis.

ECA Results
Lagged correlations between the cumulative rainfall residuals and CSR-M TWSA were quantified using the method described under Section 3. The lags for all 35 basins are already shown in Figure 1 (normalized by 12 months), which are binned into different color groups according to the maximum correlation.The actual numerical values are reported under the 7th and 8th columns in Table 1.The lagged correlation between cumulative P residuals and TWSA ranges from 0.19 Left axis: standardized TWSA smoothed using 3-month moving average (dark green line with filled dots); right axis: SPI (filled gray areas) and ±1 thresholds (horizontal dotted lines).TWSA anomalous events were extracted using the ±1 thresholds.

ECA Results
Lagged correlations between the cumulative rainfall residuals and CSR-M TWSA were quantified using the method described under Section 3. The lags for all 35 basins are already shown in Figure 1 (normalized by 12 months), which are binned into different color groups according to the maximum correlation.The actual numerical values are reported under the 7th and 8th columns in Table 1.The lagged correlation between cumulative P residuals and TWSA ranges from 0.19 (Colorado River basin) to 0.96 (Chari basin), and all correlation values reported here are significant with a p-value less than 0.02.About half of the 35 basins show no lags between cumulative P residuals and TWSA.The maximum lag is up to nine months (Volga River basin).In general, non-zero lags were observed for north hemisphere, high-latitude basins, while zero lags were identified for basins in tropical climate regions.Human water management activities may be responsible for the relatively long lag (seven months) and low correlation obtained for the Colorado River basin.
For event extraction using the TC-based method, we used ±1 as thresholds.Similarly, the rainfall anomalies were also obtained by using thresholds ±1 on SPI, as mentioned in Section 2. While SPI cutoff values are linked to different levels of event severity, no such standard seems to exist for TWSA extreme events.Previously, Humphrey et al. [39] used an absolute water storage deficit value of 30 mm to extract drought events, while Thomas et al. [41] defined droughts as events when water storage deficits last three months or longer.Here, the thresholds for the normalized TWSA were chosen so that TC methods would include similar events as SPI.
For the BP method, times of abrupt changes (BPs) were used to form event series.The time lags found in the aforementioned correlation analysis were used to define the width of sliding window between TWSA and SPI for both TC and BP methods in ECA calculations.The final results of ECA are reported in Table 1.
The ECA results reported in Table 1 (9th column) suggest that almost all basin-scale abrupt changes (BPs) are perfectly related to some antecedent SPI events.Only Congo and Nile river basins have values less than 1.In contrast, events extracted using the climatology-based TC method gave much lower coincidence rates (10th column in Table 1) over many river basins.Relatively low TC coincidence rates are generally found in north hemisphere, high-latitude basins, such as in St. Lawrence (0.55) and Yukon (0.52); in heavily managed basins, such as Colorado basin (0.58); or in semiarid areas, such as Niger (0.56).These results indicate that BP detection method achieved better accuracy in uncovering extreme events, as relative to SPI events, than the TC method.
At last, we quantified the pairwise, coincidence rates of TWSA extreme events for all 35 river basins.For each river basin, event coincidence rates were calculated between that basin and all other basins.The event series were formed by using results of the BP detection method.The results were used to build a coincidence matrix as shown in Figure 6.The meaning of this BP coincidence matrix is similar to the correlation-based TWSA connectivity map, which measures how TWS in different basins tend to co-vary [87,88].The difference is that a nonparametric ECA measure is used in our analyses.
Figure 6 suggests that geographically close basins (e.g., Yenisei/Ob/Lena and Niger/Chari) generally show high coincidence rates.Basins that are affected by similar climatic circulation forcing (e.g., Orinoco/Parana, Indus/Orange) also show high coincidence in BP events.Other cases, such as the relatively high coincidence rate between the high-latitude Kolyma and Mississippi basin is intriguing.Similar patterns could be observed from the connectivity map in Sun et al. (Figure 2, [88]).Interestingly, the Orinoco River basin that is north of the Amazon River basin seems to have strong linkage with a large number of river basins.Many other basins, such as the Huanghe River basin in North China Plain, Yukon, Brahmaputra, and Murray-Darling are relatively disconnected with other basins.
Our basin-scale analyses show that the BP method provides a more robust tool for identifying hydrologic extremes.A consistent causal relationship (in the ECA sense) was found between BP and SPI events for river basins considered here.On the other hand, the event coincidence rates derived using the TC method are mixed due to the limited length of GRACE records.

Discussion
In this study, two non-parametric methods, EMD and ECA were used to understand dynamics of TWSA extremes.Overall, our EMD exploratory analysis results suggest that the extreme indicator IMF observed for the Amazon River basin (Figure 2), which exhibits strong seasonality, is not fortuitous.In fact, we were able to identify and subsequently corroborate an appropriate extreme indicator IMF for all major basins considered here (Table 1), indicating the usefulness of EMD as a preliminary analysis tool for isolating the inherent temporal scale related to TWSA extremes.As far as we know, only one previous study [44] had tried to quantify return periods of TWSA events by using GRACE climatology and the traditional frequency analysis (e.g., deriving one-in-five-year probability of anomalously high TWSA with ~10 years of data).A main strength of the EMD analysis is that it is entirely data driven and is not limited by the stationarity assumption.The observation that the mean return period of TWSA extremes ranges from 15 to 42 months for all basins studied (Table 1) is an interesting finding and warrants further studies for predictive analyses.
The event coincidence metric provides a simple measure for quantifying causal relationships between events.The event coincidence between P and TWSA events was found to be 1.0 for most of

Discussion
In this study, two non-parametric methods, EMD and ECA were used to understand dynamics of TWSA extremes.Overall, our EMD exploratory analysis results suggest that the extreme indicator IMF observed for the Amazon River basin (Figure 2), which exhibits strong seasonality, is not fortuitous.In fact, we were able to identify and subsequently corroborate an appropriate extreme indicator IMF for all major basins considered here (Table 1), indicating the usefulness of EMD as a preliminary analysis tool for isolating the inherent temporal scale related to TWSA extremes.As far as we know, only one previous study [44] had tried to quantify return periods of TWSA events by using GRACE climatology and the traditional frequency analysis (e.g., deriving one-in-five-year probability of anomalously high TWSA with ~10 years of data).A main strength of the EMD analysis is that it is entirely data driven and is not limited by the stationarity assumption.The observation that the mean return period of TWSA extremes ranges from 15 to 42 months for all basins studied (Table 1) is an interesting finding and warrants further studies for predictive analyses.
The event coincidence metric provides a simple measure for quantifying causal relationships between events.The event coincidence between P and TWSA events was found to be 1.0 for most of the river basins, regardless of the climate regime, suggesting the dominant role of precipitation in driving TWS dynamics [87].
As mentioned in Section 2.1, TWSA may be affected by uncertainty in data processing and inherent measurement errors.A standard approach in the research community has been comparing multiple GRACE products.Thus, we repeated the same basin-scale BP analyses done in Section 4.2, but using the scaled SH product from CSR (CSR-SF, SF is short for scaled factor) and the JPL-M mascon solution.The results are compared in Figure 7.
Remote Sens. 2017, 9, x FOR PEER REVIEW 18 of 25 the river basins, regardless of the climate regime, suggesting the dominant role of precipitation in driving TWS dynamics [87].
As mentioned in Section 2.1, TWSA may be affected by uncertainty in data processing and inherent measurement errors.A standard approach in the research community has been comparing multiple GRACE products.Thus, we repeated the same basin-scale BP analyses done in Section 4.2, but using the scaled SH product from CSR (CSR-SF, SF is short for scaled factor) and the JPL-M mascon solution.The results are compared in Figure 7.  Overall, Figure 7 suggests that abrupt changes detected using the three different GRACE products are generally in agreement with each other.For the Amazon River basin, all three products gave very similar results.Slight differences were noticed only around the 2010 extreme drought and the 2014 wet event.For the Mississippi River basin, the main difference is that both JPL-M and CSR-SF missed the onset of the severe 2012 drought and only gave a BP towards early 2013.CSR-SF was the only product that detected the short-term meteorological drought happened in early 2008.On the other hand, CSR-SF did not detect the BP that happened in late 2014 that corresponded to a dry event.
For the Nile River basin, the three GRACE products showed more discrepancies than in other cases, which may be attributed to the differences in the magnitudes of TWSA time series-the semiarid nature of Nile River basin and strong annual variability makes it more difficult to track TWS changes, leading to larger uncertainty in various GRACE products.For Zambezi River basin, CSR-M and JPL-M show similar BPs, but CSR-SF was not sensitive to many small-scale events between 2006 and 2014.
In this study, we mainly focused on analyzing the basin-averaged TWSA time series.The same EMD-assisted BP analysis procedure may be applied at the intra-basin scale to investigate the spatial distribution of TWSA extremes.As an illustration, Figure 8  Overall, Figure 7 suggests that abrupt changes detected using the three different GRACE products are generally in agreement with each other.For the Amazon River basin, all three products gave very similar results.Slight differences were noticed only around the 2010 extreme drought and the 2014 wet event.For the Mississippi River basin, the main difference is that both JPL-M and CSR-SF missed the onset of the severe 2012 drought and only gave a BP towards early 2013.CSR-SF was the only product that detected the short-term meteorological drought happened in early 2008.On the other hand, CSR-SF did not detect the BP that happened in late 2014 that corresponded to a dry event.
For the Nile River basin, the three GRACE products showed more discrepancies than in other cases, which may be attributed to the differences in the magnitudes of TWSA time series-the semiarid nature of Nile River basin and strong annual variability makes it more difficult to track TWS changes, leading to larger uncertainty in various GRACE products.For Zambezi River basin, CSR-M and JPL-M show similar BPs, but CSR-SF was not sensitive to many small-scale events between 2006 and 2014.
In this study, we mainly focused on analyzing the basin-averaged TWSA time series.The same EMD-assisted BP analysis procedure may be applied at the intra-basin scale to investigate the spatial distribution of TWSA extremes.As an illustration, Figure 8 shows the maximum value (in absolute magnitude) of all BPs detected at each grid block for the Mississippi River basin.The negative maximum BPs in the southeast part (Arkansas-White-River and Lower Mississippi) of the basin are related to droughts in 2005, whereas the positive maximum BPs observed in the northwest part are related to a wet season in the spring of 2013 that ended the 2012 flash drought at the region.

Conclusions
A comparative study was performed on the skill of climatology-based (TC) and breakpoint (BP) analyses for detecting hydrologic extremes over the world's 35 largest river basins by using GRACE TWSA products.The empirical mode decomposition (EMD) was used as a preliminary step before BP analyses to uncover the mean return period of extreme events for each basin.The TC method requires access to stable long-term climatology and uses a pre-defined threshold to identify anomalous events, whereas the EMD-assisted BP method iteratively identifies the optimal number, time, magnitude, and type of BPs in a time series with little subjectivity.Our comparison results show that the BP method explains the occurrences of hydrologic extremes better in all basins examined, and correlate well with precipitation anomalies.

Conclusions
A comparative study was performed on the skill of climatology-based (TC) and breakpoint (BP) analyses for detecting hydrologic extremes over the world's 35 largest river basins by using GRACE TWSA products.The empirical mode decomposition (EMD) was used as a preliminary step before BP analyses to uncover the mean return period of extreme events for each basin.The TC method requires access to stable long-term climatology and uses a pre-defined threshold to identify anomalous events, whereas the EMD-assisted BP method iteratively identifies the optimal number, time, magnitude, and type of BPs in a time series with little subjectivity.Our comparison results show that the BP method explains the occurrences of hydrologic extremes better in all basins examined, and correlate well with precipitation anomalies.
The potential merit of TWSA over forcing-based index (SPI) or simple water-budget based (e.g., PDSI) indices for monitoring large-scale, hydrologic extremes has been demonstrated in previous studies [18,35,41].This study further demonstrates that EMD-assisted BP detection method provides an alternative and more robust tool for detecting hydrologic extremes without being restricted by the limited GRACE climatology (14 years for the study period) and the everlasting stationarity question.
Currently GRACE products are limited by their spatial and temporal resolution.Recent studies demonstrate the feasibility of downscaling GRACE data spatially and temporally using data assimilation and statistical approaches, potentially extending the applicability of BP analysis to smaller regions [89][90][91][92] and shorter timescales [93,94].In the future, these new data processing techniques can be combined with EMD-assisted BP analysis presented in this study to detect the onset of extreme events.

Figure 1 .
Figure 1.Map of the 35 river basins, numbered in decreasing basin area and colored by correlation between basin-averaged cumulative P residuals and total water storage anomalies (TWSA).Pie charts indicate lags (clockwise) between cumulative P residuals and TWSA, which are normalized by 12. Zero-lags are shown as clear circles.Basin numbers correspond to basin IDs in Table1.The actual correlation values and lags are reported under P-TWS Corr and P-TWS Lag columns in Table1.

Figure 1 .
Figure 1.Map of the 35 river basins, numbered in decreasing basin area and colored by correlation between basin-averaged cumulative P residuals and total water storage anomalies (TWSA).Pie charts indicate lags (clockwise) between cumulative P residuals and TWSA, which are normalized by 12. Zero-lags are shown as clear circles.Basin numbers correspond to basin IDs in Table1.The actual correlation values and lags are reported under P-TWS Corr and P-TWS Lag columns in Table1.

Figure 2 .
Figure 2. Empirical mode decomposition (EMD) analysis of Amazon River basin TWSA.Top panel shows the raw data and bottom panel shows the EMD trend.All 5 IMFs are shown in the middle panels.In this case, IMF2 reveals all major TWSA extremes occurred at the basin during the study period.

Figure 2 .
Figure 2. Empirical mode decomposition (EMD) analysis of Amazon River basin TWSA.Top panel shows the raw data and bottom panel shows the EMD trend.All 5 IMFs are shown in the middle panels.In this case, IMF2 reveals all major TWSA extremes occurred at the basin during the study period.

Figure 6 .
Figure 6.Heat map of BP coincidence among 35 river basins.The number of BPs for each basin is labeled on the diagonal axis.Each pixel on heat map represents pairwise basin event coincidence rate.Basin pairs having coincidence rate <0.2 are shown as blanks.

Figure 6 .
Figure 6.Heat map of BP coincidence among 35 river basins.The number of BPs for each basin is labeled on the diagonal axis.Each pixel on heat map represents pairwise basin event coincidence rate.Basin pairs having coincidence rate <0.2 are shown as blanks.
shows the maximum value (in absolute magnitude) of all BPs detected at each grid block for the Mississippi River basin.The negative maximum BPs in the southeast part (Arkansas-White-River and Lower Mississippi) of the basin are related to droughts in 2005, whereas the positive maximum BPs observed in the northwest part are related to a wet season in the spring of 2013 that ended the 2012 flash drought at the region.Remote Sens. 2017, 9, x FOR PEER REVIEW 19 of 25

Figure 8 .
Figure 8. Maximum BP (in absolute magnitude) obtained for each grid block in Mississippi River basin.

Figure 8 .
Figure 8. Maximum BP (in absolute magnitude) obtained for each grid block in Mississippi River basin.

Table 1 .
The actual correlation values and lags are reported under P-TWS Corr and P-TWS Lag columns in

Table 1 .
The actual correlation values and lags are reported under P-TWS Corr and P-TWS Lag columns in

Table 1 .
List of 35 global basins examined in this study.AI represents aridity index, CoinBP and CoinTC represent event coincidence rates calculated using BP and TC methods, respectively, P-TWSA lag is in months.Under the Climate column: H-human, SH-semihumid, SA-semiarid, and A-arid.Irrigation area is given as percentage of basin areas.