Global trend analysis of the MODIS drought severity index

Recently, Mu et al. (2013) compiled an open access database of a remotely sensed global drought severity index (DSI) based on MODIS (Moderate Resolution Imaging Spectroradiometer) satellite measurements covering a continuous period of 12 years. The highest spatial resolution is 0.05× 0.05 in the geographic band between 60 S and 80 N latitudes (more than 4.9 million locations over land). Here we present a global trend analysis of these satellitebased DSI time series in order to identify geographic locations where either positive or negative trends are statistically significant. Our main result is that 17.34 % of land areas exhibit significant trends of drying or wetting, and these sites constitute geographically connected regions. Since a DSI value conveys local characterization at a given site, we argue that usual field significance tests cannot provide more information about the observations than the presented analysis. The relatively short period of 12 years hinders linking the trends to global climate change; however, we think that the observations might be related to slow (decadal) modes of natural climate variability or anthropogenic impacts.


Introduction
Severe droughts or floods are threatening events for both ecosystems and human society.There are several indices used widely for drought assessment integrating large amounts of data (precipitation, snowpack, stream flow, etc.).The best known index is probably the Palmer drought severity index (PDSI) (Palmer, 1968;Alley, 1984) determined by monthly water supply (precipitation), water outputs (evapo-ration and runoff), and preceding soil water status.New variants of the original approach have emerged in order to overcome some limitations of the Palmer model (Alley, 1984;Keyantas and Dracup, 2002), such as the self-calibrating PDSI by Wells et al. (2004) or PDSI incorporating improved formulations for potential evapotranspiration (Heim Jr., 2002).Remote sensing data from the Moderate Resolution Imaging Spectroradiometer (MODIS) combined with NCEP (National Centers for Environmental Prediction) reanalysis meteorological records and statistical procedures together have supported the development of an evaporative drought index (EDI) by Yao et al. (2010Yao et al. ( , 2014) ) at 4 km spatial and 1 month temporal resolutions.Nevertheless, the development and improvement of drought indices are incomplete tasks and numerous challenges remain for the future (Vicente-Serrano et al., 2011;Trenberth et al., 2014).
To the best of our knowledge, the most comprehensive and longest PDSI trend analysis has been provided by Dai et al. (2004).A monthly PDSI data set from 1870 to 2002 has been derived using historical precipitation and temperature data for global land areas on a grid of 2.5 • × 2.5 • .An empirical orthogonal function (EOF) analysis resulted in a linear trend in the twentieth century, with drying over northern and southern Africa, the Middle East, Mongolia, and eastern Australia, and moistening over the United States, Argentina, and parts of Eurasia (Dai et al., 2004).A follow-up study by Dai (2011) compared the original and three other variants of PDSI records, but the main conclusion remained the same: warming in the second half of the last century is responsible for much of the drying trend over several land areas (Dai, 2013) change may not cause droughts but it is expected that when droughts occur they are likely to set in quicker and be more intense (Trenberth et al., 2014).However, similarly to the open questions on an optimal definition of a drought index, debates on the trends are also not entirely closed (Sheffield et al., 2012;Damberg and AghaKouchak, 2013;Spinoni et al., 2013;Trenberth et al., 2014).
Here we report on a global trend analysis of the remotely sensed DSI time series constructed by Mu et al. (2013).Special attention is paid to testing the statistical significance of local linear trends by data shuffling (see Sect. 2).The main result is that 17.34 % of the land area exhibits significant trends of either signs (12.03 % drying and 5.31 % wetting) and most of these locations form large, geographically connected areas.We emphasize that the usual field significance tests (Benjamini and Hochberg, 1995;Douglas et al., 2000;Ventura et al., 2004;Wilks, 2006;Renard et al., 2008) cannot give more reliable estimates, because a DSI value as defined provides a fully local characterization, and the same numerical value can be related to very different local circumstances.The relatively short period of 12 years hinders linking the trends to global climate change; however, we think that the observations might reveal a slow (decadal) mode of natural climate variability.Correlations with other atmospheric and oceanic variables are found at various (statistically insignificant) levels; therefore, at the moment we cannot prove any causal relationship or propose a solid explanation of the observations.

Data and methods
In order to better exploit the strengths of continuous satellite observations, Mu et al. (2013) recently developed a remotely sensed global drought severity index (DSI) and compiled an open access database spanning 12 years, between 2000 and 2011, at a temporal resolution of 8 days.The highest spatial resolution is around 5 km (0.05 • × 0.05 • ) with an almost global coverage.The definition of DSI incorporates the normalized difference vegetation index (MOD13 NDVI product), together with the evapotranspiration and potential evapotranspiration ratio data (MOD16 ET / PET product).
The NDVI vegetation index measures the fraction of photosynthetically active radiation absorbed by plant canopies, basically the amount of vegetation present on the ground (Zhou et al., 2001;Huete et al., 2002;Fensholt and Proud, 2012).By design, the dimensionless NDVI is a transformation of the near-infrared/red spectral reflectance ratio, and it varies between −1.0 and +1.0.Relatively large negative values occur when the red reflectance exceeds the near-infrared one corresponding to water surfaces, values around zero are characteristic for barren areas, while positive values span from grasslands to midlatitude forests and to tropical jungles.
Terrestrial evapotranspiration (ET) is the sum of evaporation and plant transpiration from land surface to the atmo-sphere.The computation algorithm of MOD16 ET is based on the theory of the Penman-Monteith energy balance by using remote sensing inputs of the leaf area index, land cover, albedo, and enhanced vegetation index, as well as meteorological parameter values of radiation, air temperature, pressure, and humidity (Mu et al., 2007(Mu et al., , 2011)).Over a sufficiently long time period, ET is less than or equal to precipitation for most vegetated geographic locations, apart from sites where the irrigation or subsurface water supply may shift the balance.Tropical forests have the highest ET values, dry areas and areas with short growing seasons exhibit low ET, while values for temperate and boreal forests lie between the two extremes.Potential evapotranspiration (PET) is the amount of water that would be evaporated and transpired in case of hypothetically infinite water availability, incorporating the energy available for evaporation and the ability of the lower atmosphere to transport evaporated moisture away from the land surface.ET and PET are expressed in terms of depth of water (mm) for a given time period (day, week, month, etc.); thus, the ratio ET / PET has a value of 1.0 where evapotranspiration fully satisfies potential conditions, and declines toward zero where the surface dries.
A given DSI value is obtained by standardization of the sum of previously and separately standardized NDVI index and ET / PET ratio (Mu et al., 2013): Temporal mean values * and standard deviations σ * are determined over the available time period for each grid point separately.It is an important detail that DSI is derived using ET / PET without NDVI during the classified dormant season, because of greater noise in the non-growing season NDVI signal (Mu et al., 2013).Permanently non-vegetated locations such as deserts, high mountains, extended lakes, or large cities cannot provide useful input for DSI data; therefore, such grid points are filtered out by a quality assessment procedure (Fensholt and Proud, 2012).Note that zero DSI values are not equivalent with missing data, the standardization procedure defined by Eqs. ( 1) and ( 2 8 days, apart from the necessary cuts at the end of each year. Example time series and linear trends are shown in Fig. 1 for three nearby locations (at the same latitude) in Argentina, where significant negative trends are identified (see below).Note that the time series exhibit apparently weak seasonal and annual variations, in spite of the fact that the climate of the province of Córdoba is humid subtropical with four marked seasons.This is mostly because the ET / PET term in DSI (see Eq. 1) has a feeble seasonal variability in many places (Lafleur et al., 2005); moreover, NDVI is incorporated only during the classified snow-free growing seasons as noticed before.Statistical significance of slopes of linear fits is verified by the standard permutation test (Manly, 2007).Since most of the DSI signals exhibit marked persistence on timescales of weeks or even months (see Fig. 1), the basic unit of data shuffling was one whole calendar year.We cut a given record into 12 pieces and built a test set from randomly shuffled and glued years.The mean slope and standard deviation (σ ) were determined, and we accepted a fitted slope of a measured record to be significant when its distance from zero was larger than 2σ of its own test set.Figure 2 illustrates that a test set size of 100 samples provides essentially the same statistics as 100 000 random samples; however, for the sake of minimizing errors we fixed the test set size as 1000 samples.Obviously, the larger the sample size the closer the histogram of obtained slopes to a pure Gaussian shape (not shown here); however, the mean and variance do not show detectable sensitivity to the size of the test sets (Fig. 2).

Results and discussion
Condensation of the multitude of information into a single number such as DSI has the drawback that a given numerical value does not carry any explicit meaning.Two zeros at two distant grid points convey the message only that the situations at both sites are compatible with the local 12-year mean state of the very sites considering the functional status of the ecosystems.Nonzero DSI values of distant locations   The histogram of all fitted slopes for standardized DSI records is shown in Fig. 3 (blue bars), note the logarithmic vertical scale.The shape is clearly not Gaussian with a mean value of −0.00875 and standard deviation of 0.04971 year −1 .Statistically significant local slopes are obtained for 852 373 sites (17.34 %) at 2σ level, the numbers for 2.5σ and 3σ limits are 269 900 (5.49%) and 16 321 (0.33 %), respectively.Negative (drying) trends (12.03 %) have a mean slope of −0.05466 ± 0.04535 year −1 , while significant positive slopes (5.31 %) are around 0.02892 ± 0.04685 year −1 .Spatial correlations certainly bias these numerical values; however, we argue that a proper interpretation of DSI must be entirely based on local information, so field significance tests (Benjamini and Hochberg, 1995;Douglas et al., 2000;Ventura et al., 2004;Wilks, 2006;Renard et al., 2008) cannot be applied here for quantitatively incomparable records.
The main result of the present analysis is illustrated in Fig. 4. Note that reddish and blueish colouration indicates sites of significant DSI trends, and the zero level is not white (the latter is used to identify missing locations).There are several geographically connected areas exhibiting "drying" (South America, middle Asia or sub-equatorial Africa) or "wetting" (middle and northern Africa, Indian Peninsula or eastern Spain) tendencies.
As for an interpretation of the observed trends, a further complication arises from the composite aspect of DSI defined by Eq. (1).NDVI and ET / PET are far from being independent characteristics; however, they are not strictly coupled.Therefore, a given trend can be the result of a dominating tendency of changes in one of the two terms or in both of them.Peculiar situations can occur when the two partial trends cancel each other out, and this is not a purely theoretical possibility.Fortunately, Fensholt and Proud (2012) per-formed a linear trend analysis of monthly MODIS NDVI data over the period 2000-2010.They found significant slopes for 11.8 % of the pixels (the same spatial resolution as for DSI) on a global scale: 5.4 % characterized by positive trends and 6.3 % with negative trends.The comparison of Fig. 4 and Fig. 2a of Fensholt and Proud (2012) reveals a high level of correspondence; however, DSI is definitely a more sensitive indicator of negative trends.The most important difference already mentioned is that Fensholt and Proud (2012) evaluated NDVI data in the growing seasons, while DSI provides continuous records by means of ET / PET ratio.For this reason, the areas indicated by blue rectangles in Fig. 4 can be considered as exhibiting ET / PET dominated trends.The majority of such grid points indicating statistically significant negative DSI tendencies are located in the middle of the Eurasian continent, in Southeast Asia, in the northern part of South America, and in the eastern part of middle Africa.Three smaller regions denoted by orange rectangles in Fig. 4 (northern Mexico, South Africa, and northern Australia) show significant NDVI trends; however, DSI records do not obey similar (or exhibit much weaker) tendencies.These are candidate areas where NDVI and ET / PET terms cancel each other out, but further local studies are necessary to get a clear interpretation.
In order to demonstrate the power of high resolution mapping, we illustrate enlargements of South America (Figs. 5,  6) and India (Fig. 7).In both cases, we show examples where significant trends at isolated locations have a plausible explanation, and they are not observational errors.
In Fig. 6, an area at the border of Uruguay and Brazil is depicted (Rivera Department), where the satellite picture of the yellow rectangle clearly indicates intense agricultural activity (see the straight cuts between forested area and crop fields).Both Fig. 4 and the NDVI analysis of Fensholt and Proud (2012) distinguish the very spot with a significant positive trend inside an extended area of decreasing vegetation greenness signal (Fig. 5).Essentially the same observation is formulated by Barbosa et al. (2015), who studied spatial patterns of standardized difference vegetation index (SDVI) and standardized precipitation index (SPI) for the interval 1998-2012 in South America.(Monthly SDVI is equivalent with the first term in Eq. ( 1), precipitation data are treated analogously.)The map in Fig. 6a of Barbosa et al. (2015) is consistent with Fig. 5; moreover, regions of definitely decreasing rainfall trends largely overlap with the red areas in Fig. 5.The main conclusion of Barbosa et al. (2015) is that in 46 % of the study area, significant decreasing or increasing greenness tendencies cannot be linked to changes in rainfall over time, indicating human impacts or the influence of other climatic factors, such as temperature.
A similarly large area, where DSI trends have the opposite sign is mapped in Fig. 7.Note that the tendencies here are also NDVI dominated (see Fig. 4b); however, the number of pixels with statistically significant slopes is definitely larger than in Fensholt and Proud (2012).The same greening process is identified by de Jong et al. ( 2011) for a longer period of 1981-2006 in the Global Inventory Modeling and Mapping Studies (GIMMS) NDVI data set (Bai et al., 2008).The yellow circle in Fig. 7 (see also the Inset) identifies the pixels around the Indira Sagar Reservoir, which is constructed as the key project of a large multipurpose river basin development on the Narmada River.Full-scale energy production started in 2005, just in the middle of the observed period.It is most likely that the gradual filling up of the reservoir resulted in a decreasing vegetated area and thus the DSI signal reflects a decreasing trend where degrading NDVI contribution in Eq. ( 1) is the main factor.
The fact that positive DSI trends are significant in a much larger area than NDVI greening in India indicates that ET / PET contribution can be equally important.Since the Indian Peninsula is a "water-limited" area considering long- term evaporation (McVicar et al., 2012), precipitation tendencies might play a key role in interpreting DSI shifts.Kothyari and Singh (1996) studied long-term time series of summer monsoon rainfall and identified decadal departures above and below the long time average alternatively for three consecutive decades.Singh and Sontakke (2002) reported on an increase in extreme rainfall events over northwest India during the summer monsoon and a decline of the number of rainy days along east coastal stations in the past decades, resulting in a westward shift in rainfall activities.Similarly, Murumkar and Arya (2014) demonstrated by means of wavelet analysis that prominent annual rainfall periods exist ranging from 2 to 8 years at all the studied stations after the 1960s.Large-scale spatial and temporal correlations between the trends of rainfall and temperature are found by Subash and Sikka (2013), without a direct relationship between increasing rainfall and increasing maximum temperature of monthly or seasonal patterns over meteorological subdivisions of India.As for the particular area, even glaciers can be listed as candidate explanatory factors, since they influence runoff into lowland rivers and recharge riverfed aquifers (Bolch et al., 2012).
A satisfactory interpretation of DSI trends certainly requires several explanatory variables, since the index itself is a rather complex quantity.Figure 8 illustrates example time series of standard atmospheric parameters (daily mean temperature, precipitation and relative humidity) for two weather stations in eastern Spain (Tortosa and Zaragoza) seemingly embedded in a rather large wetting region (data from Klein Tank et al., 2002).There is no sign of any trends in the time series during the study period, even when we check the previous decades (not plotted here) (Vicente-Serrano, 2007).This semiarid region is water-limited too (McVicar et al., 2012); therefore, precipitation trends should be correlated with DSI shifts, especially in Tortosa.
The relatively short period of 12 years is not enough to connect the results with global climate change.We think that the observed significant DSI trends over extended geographic areas might be related to a decadal mode of the natural climate variability or extended anthropogenic impacts, or a combination thereof.Apart from time-span limitations, recent studies on drought trends (Sheffield et al., 2012;Dai, 2011Dai, , 2013) ) lead to somewhat controversial conclusions, as noticed by Trenberth et al. (2014).One major issue in determining reliable long-term trends in drought due to climate change is to separate the effects of natural variability, especially El Niño-Southern Oscillation (ENSO) (Trenberth et al., 2014).During El Niño events, the main rainfall systems in the tropics move eastward over the tropical Pacific Ocean leaving weakened monsoons behind (Panda and Kumar, 2014; Barreiro et al., 2014).In the La Niña phase, dry areas are more common in places where it is wet during El Niño events.Indeed, as Miralles et al. (2014) pointed out in a recent study, ENSO dominated the multi-decadal variability of terrestrial evaporation at the global scale.Their main conclusion is that the recent decadal decline in global average continental evaporation is not the consequence of a persistent reorganization of the terrestrial water cycle, rather it is an indication of transitions to El Niño conditions (Miralles et al., 2014).

Conclusions
The objective of the present work was to perform a global linear trend analysis of the remotely sensed drought severity index (DSI) compiled by Mu et al. (2013).The methodology was very similar to other studies on vegetation indices (NDVI or SDVI) (de Jong et al., 2011;Fensholt and Proud, 2012;Barbosa et al., 2015) with the main difference that our significance test is based on the direct permutation method.Detailed comparisons mostly with the results of Fensholt and Proud (2012) indicate that DSI performs somewhat better in detecting significant negative trends.This can be a consequence of the definition of DSI, where the ET / PET ratio provides a continuous contribution to the signal, while NDVI participates only in the growing season.On the other hand, this definition complicates the interpretation of observed trends, because the two terms are not functionally related; therefore, the separation of contributions is not trivial.We have illustrated the power of high resolution mapping by zooming in on restricted regions and providing reasonable explanations of why local trends can have opposite signs to the surrounding extended area.
Work is in progress in three directions in order to find a better explanation of the observed trends.Firstly, it is a plausible goal to repeat the analysis separately for the two terms (standardized NDVI and ET / PET ratio) to identify precisely the role of these factors.Secondly, it is reasonable to compare DSI with the many existing drought indices.This is a demanding task, mostly because the validation of the various time series certainly requires direct comparisons with field observations.Thirdly, decadal trends over extended geographic areas call for an explanation related to global climate change, especially when the subject has such societal implication as a drought severity index.This is a highly nontrivial problem too, because the separation of natural climate variability from unambiguous climate shifts is hindered by the length and reliability of available data (Trenberth et al., 2014).Nevertheless cross-correlation analyses with relevant atmospheric variables are necessary to begin the procedure.Candidate indices are El Niño-Southern Oscillation, Northern Annular Mode/North Atlantic Oscillation, Southern Annular Mode, sea surface temperature (SST) anomalies, sea ice cover (SIC), Atlantic Multidecadal Variability, etc.

Figure 2 .
Figure 2. Significance analysis of fitted slopes by using test sets of randomly shuffled whole year DSI records.The ensemble means and standard deviations are plotted for 10 sites along latitude 30.025 • S, evenly spaced by 0.05 • westward starting from 63.775 • W. Test set sizes are 100 (black circles), 1000 (red squares), 10 000 (orange diamonds), and 100 000 (blue stars).Maroon crosses indicate fitted slopes (year −1 ) for the original measured time series.

Figure 3 .
Figure 3. Histogram of all fitted slopes (blue) over the continents, where DSI data are available (4 914 440 locations).Its is clearly not a Gaussian shape (note the logarithmic vertical scale), the global mean value is −0.00875 ± 0.04971 year −1 .Orange bars denote the histogram of mean slopes computed for 1000 randomly shuffled test sets for each geographic location.Red bars indicate significant slopes at 2σ or higher level.

Figure 4 .
Figure 4. Geographic distribution of sites where DSI trends are significant at 2σ or higher levels.Linear trend slopes in units of year −1 are color coded.Blue rectangles denote areas where significant DSI trends are ET / PET dominated (based on a comparison with the map of Fig. 2a in Fensholt and Proud, 2012).Orange rectangles identify regions where NDVI trends are significant (see Fig. 2a in Fensholt and Proud, 2012), but DSI trends are missing or very weak.

Figure 5 .Figure 6 .
Figure 5. Enlargement view of South America at the highest spatial resolution of 0.05 • × 0.05 • .Trends in units of year −1 are color coded.

Figure 7 .
Figure 7. Enlargement view of the Indian Peninsula at the highest spatial resolution of 0.05 • × 0.05 • .Trends in units of year −1 are color coded.Yellow circle locate the Indira Sagar Reservoir, approx.22.17 • N, 76.65 • E (see the inset).

Published by Copernicus Publications on behalf of the European Geosciences Union. P. I. Orvos et al.: Global trends of drought severity index
. Increased heating itself from global climate