Secular trend of global drought since 1950

Drought is a recurring extreme climate event over most parts of the world featured by long duration and low predictability. The secular trend of drought is of particular interest for investigators in agriculture, climate change and sustainability domains. In this study, we applied the ensemble empirical mode decomposition (EEMD) method and analyzed the spatio-temporal characteristics of the secular trends of meteorological drought over global land surface during the period 1950–2015 using a self-calibrating Palmer Drought Severity Index (PDSI) product. We found that there were 25.98% PDSI samples had turning point, namely the shift of trend, in the corresponding secular trend series; the probability distribution of the turning points position (period) extracted by EEMD closely follows a normal distribution with mean value at Nov. 1981. We showed that there is large discrepancy in the secular trend types extracted by EEMD and Mann–Kendall test, and exemplified the risk of using a monotonic trend to capture the changes of the intrinsic secular trend of PDSI series. We suggested that there was an accelerated drying trend over global land surface as a whole, but large areas with wetting trend existed in the meantime, especially at the high latitudes in the northern hemisphere. Additionally, we found that the PDSI secular trend change rate exhibits a multidecadal variability of about 50 years or so and it implies a potential relationship with periodic variations of the oceanic and atmospheric current. We showed that the secular trend of PDSI series extracted by EEMD could provide more detailed spatio-temporal characteristics, featured by the shifts of trend and nonlinear property of the secular trend, of global drought than that of the non-parametric or linear regression methods. The secular trend of PDSI could present more insights about the transition and progress of wetting/drying trend over global land surface at multidecadal scale.


Introduction
Drought is one of the most serious natural disasters characterized by abnormal reduction of precipitation relative to local normal precipitation and extended from months to years from a meteorological perspective (Mishra and Singh 2010, Dai 2011b, Trenberth et al 2013. The droughts in literature could be categorized as three physical types: meteorological, agricultural, and hydrological, in addition, the socioeconomic drought if the societal water resource demand exceeds some kind of natural supply associated with meteorological, agricultural, and hydrological drought (American Meteorological Society 1997, Keyantash and Dracup 2002). The meteorological drought is defined as shortage of precipitation at local scale over a period of time and is closely correlated with hydrological and agricultural droughts (Lu et al 2017, Penagallardo et al 2019, and it is one of the most extensively studied kind of drought due to the availability of relative long historical gauge observations, analysis data sets, and associated meteorological variables (e.g. temperature and wind) (Chen et al 2002, Huffman et al 2009, Harris et al 2014. In this study, we focused on the spatio-temporal characteristics of meteorological drought at global scale. Drought is a recurring extreme climate event over most parts of the world even in wet and humid areas and featured by its long duration and low predictability (Dai 2011b). During the past decades, there has been an increased interest in studying drought in the context of climate change at decadal to century time scale, as climate change has profound impacts on the spatio-temporal distributions of key meteorological factors correlated with drought, e.g. precipitation, surface vapor pressure deficit, and temperature (Seager et al 2010, Trenberth 2011, Ji et al 2014. For example, numerous studies have tried to correlate regional droughts with the oscillations of sea surface temperature, e.g. the El Niño/Southern Oscillation (ENSO), by teleconnections (Palmer and Branković 1989, Hong and Kalnay 2000, Shabbar and Skinner 2004, Trenberth et al 2013. Another serious concern is about the widespread drying trend over land which could be partially attributed to the sustained global warming, and it is estimated that dry areas has increased by about 1.74% of the global land area per decade from 1950 to 2008 (Hansen et al 2010, Dai 2011a, 2013, Ji et al 2014. Analysis based on the Coupled Model Intercomparison Project phase 5 (CMIP5) showed that dry land area will increase by 23% and 11% relative to the 1961-1990 baseline under representative concentration pathways (RCPs) RCP8.5 and RCP4.5 respectively, by the end of this century (Huang et al 2016); similar results could be found in Feng and Fu (2013).
The secular trend of time series climate data, including drought, is of particular interest for investigators in climate change domain. The most commonly used trend analysis methods, e.g. the simple linear trend, nonlinear regression analysis or more complicated trend analysis methods like Fourier-or wavelet-based filtering, are under the premise of linear and stationary assumption (Wu et al 2007, 2011, Yin et al 2017. For drought, it is not unexpected that many studies using single linear or non-parametric trend to describe the long-term characteristics of drought (e.g. Burke et al 2006, Dai 2013, Ficklin et al 2015, but this kind of trend may mask some opposite trends or turning points in the original time series data while this kind of information might be crucial to a particular system. Additionally, the precipitation and drought might violate the linear assumption and the nonstationarity might be common as regional climates systematically change (Chou et al 2012, Zhang and Cong 2014, Apurv and Cai 2019. As suggested by Wu et al (2007), the aforementioned trend extraction methods subjectively fit data with a prior determined functional forms, while there is no physical bases to support the data to follow a specific functional form. Also, the linearly trend extraction methods are sensitive to short-term fluctuations and abrupt changes, and the components of a long-period oscillation may be confused with shortterm trends (Verbesselt et Verbesselt et al 2010, Forkel et al 2013 implicitly assume that the change rate abruptly converts at the turning or breaking point, nevertheless the drought generally is a gradual process so that the linear trend change detection methods may not be justified. As an alternative, the ensemble empirical mode decomposition (EEMD) is a temporal local analysis method suitable for determining the intrinsic trend from nonstationary and nonlinear processes (Wu et al 2007, Wu andHuang 2009). The EEMD is a noise-assisted data analysis method based on the EMD method, which adaptively decomposes the data into a series different intrinsic modes of oscillations (IMFs) and a residual component, originally developed by Huang et al (1998). The residual component of EEMD is monotonic or contains at most one extreme and represents the secular trend of the time series process (Wu et al 2011, Ji et al 2014. The EEMD method has been an appealing and insightful means in revealing the secular trends of underlying environmental and physical patterns in various applications, e.g. climate change (Qian et al 2009, Ji et al 2014, remote sensing (Yin et al 2017, Pan et al 2018, and meteorology (Franzke 2012).
The objective of this study is to quantitatively investigate the spatio-temporal characteristics about the secular trend, including its turning point, of global drought in recent decades since 1950. The secular trend of a drought index product is analyzed per pixel globally using the EEMD. We applied the commonly used non-parametric methods, e.g. Mann-Kendal and Pettitt's tests, to obtain the trend and breaking point period (month) and to compare with the corresponding results from EEMD. We analyzed the long term spatio-temporal characteristics of global drought based on the secular trend. The paper is arranged as follows. Section 2 introduces the data and methods. The spatio-temporal characteristics of the secular trend of global drought are presented in section 3, followed by the discussions in section 4 and summary in section 5.

Data
For spatio-temporal assessment of drought, many quantitative drought indices designed for different types of drought have been developed and the Palmer Drought Severity Index (PDSI) is the most commonly used one for the meteorological drought (Heim 2002, Dai 2011a). The PDSI is a relative metric and was originally formulated by Palmer (1965) which incorporates antecedent precipitation and current moisture supply (precipitation), and moisture demand (potential evapotranspiration, PET) into a hydrological accounting system within a two-layer bucket soil model for soil moisture calculation (Heim 2002). The standardized PDSI ranges from −10 (dry) to 10 (wet) with value less/greater than −3/3 representing severe to extreme drought/very to extremely wet. To improve the comparability of PDSI caused by the empirically derived climatic characteristic and duration factors in the original PDSI, Wells et al (2004) proposed a self-calibrating PDSI (sc_PDSI) by automatically calculating aforementioned parameters with local historical climatic data. Another widely applied improvement to Palmer's PDSI is to replace the Thornthwaite equation (Thornthwaite 1948) with more sophisticated Penman-Monteith equation (Thornthwaite Allen et al 1998) to get a more realistic PET (Dai 2011a). The PDSI and sc_PDSI have been widely applied to investigate the spatio-temporal characteristics of drought, and they are more suitable to assess the effect of global warming on drought and wet spells than other statistically based drought indices as supported by the Palmer's water balance model, which incorporates both precipitation and surface air temperature (Mishra and Singh 2010, Dai 2011a, 2013, Sheffield et al 2012. However, it is worth noting that quantitative interpretation of the dryness/wetness for a given PDSI value cannot be separated from local mean climate conditions (Dai et al 2004).
We used a monthly sc_PDSI_pm (self-calibrating PDSI with potential evapotranspiration calculated by the Penman Monteith equation) product (Dai et al 2004, Dai 2011a over global land areas in this study. The spatial resolution of this data is 2.5 • × 2.5 • , and its temporal range is from Jan. 1850 to Dec. 2014 currently. To calculate sc_PDSI_pm, as stated by Dai (2011a), the monthly surface air temperature was derived from CRUTEM3 anomalies and CRU climatology; the monthly precipitation data was produced by merging datasets from Dai et al (1997) for 1850-1947, Chen et al (2002) for 1948-1978, and the Global Precipitation Climatology Project (GPCP) (Huffman et al 2009(Huffman et al ) since 1979. The additional surface variables, e.g. specific humidity, wind speed, air pressure, and net solar radiation, since 1948 were taken from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis. However, it is noted that Dai (2011a) simply used the long-term mean values of these surface variables before 1948 due to the scarcity of data and this could cause the sc_PDSI_pm product before 1948 has limited variations. To avoid influence of the uncertainties in the forcing data, we limited our study period in 1950-2015 and excluded the grids with missing data. There are 2433 qualified grids covering most of the land areas globally during the study period. The excluded grids are about 11.78% of all the land surface PDSI grids and mainly concentrated in the northern Africa around the Sahara desert and high latitude areas in the northern hemisphere, i.e. either extreme arid or permafrost zones where meteorological data are not readily available. It might be assumed that the excluded data would not have much influence on the results.

EEMD
The EEMD is designed to overcome the so-called mode mixing problem in EMD (Wu and Huang 2009) by taking the mean of an ensemble of trials, i.e. adding white noise of finite amplitude to the original signal x (t), as the true IMF component if the trials are enough for the noise in each trail to cancel out in the ensemble mean, and it is not sensitive to noise and addition of new data on the original time series data.
The EEMD decomposes the data x (t) into a series IMFs, c j , with decreasing frequencies and a residual trend component, r n , as follows: where n is the total number of IMFs and is close to log 2 N, with N being the length of x (t), and r n is the residual component of the data. The inherent fluctuations in the data series have been removed by EEMD and no more IMFs could be extracted from r n , except the secular trend information (Wu et al 2007, Wu andHuang 2009). In this study, the ensemble number of EEMD is set to 400, and the amplitude of the white noise is set to 0.2 standard deviation of the corresponding data at each grid.

Non-parametric methods
In this study, we also used two non-parametric methods, i.e. Mann-Kendal and Pettitt's tests, and the Breaks For Additive Seasonal and Trend (BFAST) to verify (1) the extra information contained in the secular trend identified by EEMD while omitted by the traditional monotonic or linear trend methods, and to compare (2) the turning point positions extracted by EEMD with one classical change-point detection method, i.e. the Pettitt's test, and BFAST, a generic change detection approach. The rank-based Mann-Kendall (MK) trend test (Mann 1945, Kendall 1975) is a widely used nonparametric method to detect the significance of trends for non-normally distributed data, which is common in hydrological and climate time series, and is less sensitive to outliers (Yue et al 2002a). The null hypothesis of MK test is that the time series has no trend, while the alternative hypothesis is that a monotonic (upward or downward) trend exists. As the distribution of sc_PDSI is non-Gaussian in many cases (Dai 2011a), we suggested that MK test might be more suitable than the linear regression method to detect the significance of the monotonic trends of sc_PDSI at each cell. It is noted that the positive/negative serial correlation, a common phenomenon in hydroclimatic time series, might increase/decrease the magnitude of trend detected by MK even in the absence of a trend, so we applied the pre-whitening technique to remove serial correlation, especially the lag 1 autoregressive (AR(1)) process in the original time series before EEMD analysis (von Storch 1995, Yue andWang 2002b).
The Pettitt's test (Pettitt 1979) is a rank-based test method commonly used to detect the non-parametric single abrupt change in univariate time series data set, especially in hydrological or climate series (Wijngaard et al 2003, Rougé et al 2013. The Pettitt's test identifies the significance and timing of the major change in a time series, where the change means a shift in the mean. We compared the temporal difference between the positions of the change-points detected by the Pettitt's test and the turning points in the secular trends extracted by EEMD. We also applied the generic trajectory based change detection method BFAST (Verbesselt et al 2010) to detect the turning points in PDSI series. The BFAST decomposes the time series into trend, seasonal and noise components, and then follows an iterative algorithm to detect the turning points, including the magnitude and direction of the change (Verbesselt et al 2010).

Statistical analysis of trends and turning points
According to the MK test, there were 40.93% and 59.07% PDSI samples followed upward and downward trends (including statistically insignificant samples), respectively; in contrast, the monotonic increase and decrease ratios of the secular trend of r n in EEMD were 38.51% and 35.51%, respectively (table 1). From the perspective of the secular trend of PDSI extracted by EEMD, only 21.09% (monotonic increase) and 26.39% (monotonic decrease) samples' trend types coincided with MK test (including statistically insignificant samples in MK test). Except the large discrepancy shown in the monotonic ones, total 25.98% samples (in which 12.17% decrease-increase and 13.81% increase-decrease) had secular trend with turning point in EEMD, but this information was totally masked in the MK test.
There were 632, 2266, and 2222 turning points (the PDSI sample size is 2433) for EEMD, Pettitt's test, andBFAST in 1955-2010, respectively; the numbers of turning points falling outside of 1955-2010 were 25, 32, and 164 for EEMD, Pettitt's test, and BFAST accordingly. In Pettitt's test, the statistically insignificant turning points were excluded; while in BFAST, only one major turning point was kept for a particular PDSI series. It showed that the Pettitt's test and BFAST tend to find a turning point for each PDSI series and, compare with EEMD, are much sensitive to the fluctuation of data.
The medians of EEMD, Pettitt's test and BFAST are 375 (March 1981), 353 (May 1979), and 325 (January 1977), respectively. The mean and standard deviation of EEMD are 382.5 (Nov. 1981) and 133.8 (11.2 years), 367.7 (Aug. 1980) and 159.5 (13.3 years) for Pettitt's test, and 349.2 (Jan. 1979) and 196.1 (16.3 years) for BFAST, respectively (figure 1). All of the medians and means of the turning points of the three methods distribute around 1980. But compared with BFAST and Pettitt's test, the turning points of EEMD had the smallest interquartile range and standard deviation. Actually, the probability density of the turning points extracted by EEMD in figure 1 is more appealing. The histogram of the turning points of EEMD and a fitted normal distribution is shown in figure 2. Further, we tested the normality of the turning point distributions of the three methods using the Kolmogorov-Smirnov test (K-S test), and the results showed that EEMD had a p-value of 0.2487, while the p-values of BFAST and Pettitt's test were far less than 0.05, so only the turning points of EEMD was not significantly different from the normal distribution (a p-value greater than 0.05 means that the null hypothesis, i.e. the data is normally distributed, is not rejected in the K-S test).

Trend of global monthly averaged PDSI
All the valid PDSI samples were averaged per month, and the EEMD was applied to decompose the series and get the corresponding secular trend (figure 3). Both of the secular trend and linear regression showed that the averaged PDSI decreased continuously from 1950 to 2015; the decreased PDSI values, −0.383 and −0.390, were very close for the secular trend and linear regression, respectively. However, the secular trend of the averaged PDSI showed that the change rate accelerated throughout, as shown by the 1st derivative of this secular trend series (figure 3(b)). The change rate of PDSI in the secular trend almost quadrupled by 2015 compared with that of 1950.

Spatio-temporal features of PDSI secular trend
The valid PDSI samples were classified into four types according to the corresponding PDSI secular trend types. Figure 4 showed the spatial distribution of each type. The PDSI grids with monotonic increase or decrease secular trend types occupied about 74% of the total samples (table 1). The monotonic increase ones mainly clustered in the central and eastern North America, western Canada, South America, Europe, central Asia, Russian Far East, southern Asia, southern Africa, and Midwest Australia. The monotonic decrease ones mainly concentrated in the southern Canada, north of the South America, central and northern Africa, part of Europe (e.g. Germany and Turkey), eastern Australia, and especially, the eastern part of China and Mongolia.
The PDSI grids with decrease-increase or increase-decrease secular trend types occupied about 26% of the total samples (table 1). Figure 5 showed the spatio-temporal distribution of PDSI secular     trend types with turning point, i.e. the decreaseincrease and increase-decrease ones respectively. The mean timings of the turning points for the decrease-increase and increase-decrease secular trend types were about 1983 and 1980, respectively. Except the sporadic distributed ones, the decrease-increase samples mainly distributed in the northern Africa, northeast China, eastern Russia, and high latitude of northern America ( figure 5(a)); while the increasedecrease samples mainly distributed in the North America, Ukraine, and Romania, Western Australia, and some parts of South America and central Africa ( figure 5(b)). In North America, the samples with increasedecrease type of PDSI secular trend were obviously more than the decrease-increase ones, especially in the continental United States, showing a prevailing drying trend over this region mainly since 1980s. Similarly, the increase-decrease and decrease-increase types in South America mainly distributed in the south and north part of this continent, respectively. In Africa, the PDSI samples with decrease-increase secular trend type were dominant and mainly distributed in the south of the Sahara. In Eurasia, the decrease-increase PDSI samples mainly distributed in the northeast China and Russia; the increase-decrease ones showed no obvious clustering trend, except in Ukraine and Romania. The PDSI samples in Australia could be roughly divided into two parts, i.e. west and east; the secular trend in the west part of Australia was dominated by the monotonic increase type as shown in figure 4, but some PDSI samples showed declining trend since around 2000s.
The zonal average of PDSI secular trend, including its relative change, along the latitudinal direction is used to quantitatively visualize the temporal evolution of PDSI globally (figure 6). It showed that the averaged PDSI secular trend declined generally at most of the latitudinal bands, but in areas around 60 • N and above two bands with increased PDSI secular trend could be identified ( figure 6(a)). The relative change of the averaged PDSI secular trend showed the magnitude of PDSI secular trend change in the latitudinal direction with additional spatio-temporal pattern exposed ( figure 6(b)). It revealed several latitudinal bands with relative high PDSI change rates, for example, the areas (roughly) between 60 • -75 • N, 0 • -20 • N, and 40 • -55 • S; the areas around 60 • N and ∼70 • -75 • N showed obvious wetting trend, and the drying trend in 0 • -20 • N was most obvious.  Figures 7 and 8 showed the spatio-temporal change rate of PDSI secular trend in different ways. In figure 7, the accumulated changes of PDSI secular trend since 1950 was shown in every 10 years (15 years from 2000 to 2015). It showed that the drying trend (shades of blue) gradually emerged and mainly concentrated in Canada, north of South America, north of Africa, northeast China, north of Mongolia, and eastern Russia; and this drying trend was more highlighted since 1980s (figures 7(c)-(f)). The continental United States generally turned wet (yellow green colors) but the trend was not persisted especially since 2000s (figure 7(f)). The contrast of wetting and drying in Midwest and eastern Australia became more apparent since 1980s (figure 7(c)). In contrast, the wetting trend was most evident in the north of North America, central Asia, and eastern Russia, especially after 2000s (figure 7(f)). Figure 8 showed the change rate of PDSI secular trend every 10 years (15 years from 2000 to 2015), and it provided an opportunity to compare the change  1960, 1970, 1980, 1990, 2000, and 2015 minus 1950 for each grid, respectively. rate of PDSI secular trend in different period. Compared with the periods of 1970-1980 and 1980-1990, the changes of PDSI secular trend in the other periods were more obvious. In 1950-1960 (figure 8(a)), the relative high increase of PDSI secular trend appeared in the continental United States, north Canada, south of South America, east Europe and central Asia, south Africa, and Midwest Australia; while the relative large of PDSI secular trend decrease could be found in north Africa, central and east Russia, and northeast China; but this trend, both increase and decrease of PDSI secular trend, gradually weakened till 1990s.

PDSI secular trend changes
Since 1990s (figures 8(e)-(f)), the PDSI secular trend change intensified again globally. But the spatial pattern of changes during this period looks much like the inversion of 1950-1970 (figures 8(a) and (b)). For example, the dominant increased PDSI secular trend in the North and South Americas, east Europe, and central Asia in 1950-1970turned into decreased trend in 1990 while the prevalent decreased PDSI secular trend in the North Africa, central and east Russia, andnortheast China in 1950-1970 largely turned into increased trend since 1990s. It is also noted that the areas with increased PDSI secular trend in Alaska of the United States had largely expanded since 1990s compared with that of 1950s.

Extraction of PDSI secular trend and turning point: a methodology perspective
For trend analysis, the EEMD is unique in that the secular trend is not extracted directly from the original data but the residual component of data with no more IMFs existed so it is robust to short term fluctuations, and consequently, the turning point, if existed, is more close to reflect the long-term shift in the trend. Compared with EEMD, the MK test can only classify the trends of PDSI as upward or downward statistically. The fundamental difference between these two methods is that the MK test is a statistically method while EEMD is an analytical one. Counting both of the statistically significant and insignificant samples, about 49% samples in MK test with upward trend showed different trends in EEMD; for the downward trend samples, the ratio was about 55% (table 1). This study is not intended to verify which kind of 'trend' is more representative to a data itself, but this discrepancy exemplified the risk of using a monotonic trend to summarize a data series, especially for data with nonlinear and nonstationary characteristics in which monotonic trends might not suffice to capture the changes of the intrinsic secular trend, and this is also true for the linear regression method (Wu et al 2007).  1960-1950, (b) 1970-1960, (c) 1980-1970, (d) 1990-1980, (e) 2000-1990, and (f) 2015-2000. Note that the change rate is referred as the change of PDSI secular trend during a specified period divided by the length of that period.
The turning point analysis was limited in 1955-2010, i.e. turning point falling outside this temporal range will not be considered. This artificial setting is to reduce the influence of short-lived process at the early or late sampling period for methods like BFAST (De Jong et al 2013), although the number of turning points falling outside of 1955-2010 for all the three methods only had negligible effects to our results. The quantity disparity of turning points among EEMD, Pettitt's test and BFAST showed that the turning point defined in EEMD is fundamentally different with the turning point or more specifically, the abrupt change, detected by Pettitt's test and BFAST. The turning point in EEMD is extracted from the secular trend component of a series, and it means that the original data fluctuations have been removed before the extremum, or turning point, could be identified (Wu et al 2007). Literally, the extraction of turning point by EEMD could be categorized as a global method; while the Pettitt's test and BFAST could be regarded as local methods, in which both methods apply varied algorithms to segment the series at the turning point(s). For example, the change trajectories at both sides of a specific turning point (abrupt change) may follow the same direction in Pettitt's test and BFAST, but this kind of point will very likely to be omitted in EEMD because the secular trend of the series might be a monotonic function locally. Considering that the Pettitt's test and BFAST methods find turning points for almost every PDSI sample in this study, it implies that the physical meanings of the retrieved turning points might be in doubt and this could be partially reflected in the probability distribution of the turning points shown in figure 1.
The means of the turning points in EEMD, Pettitt's test and BFAST were all around 1980, but the corresponding probability distributions (figure 1) showed that the turning points of EEMD is more statistically meaningful and pretty close to follow a normal distribution (figure 2). It is noted that the mean time of turning points in EEMD is coincidence with the warming trend over global land surface around 1980 (Vose et al 2005, Ren et al 2013, Ji et al 2014, and it is also suggested that recent surface warming could have enhanced evaporative demand over land and coupled back to the drying since the 1980s (Dai 2011b). Additionally, it is also suggested that the 1982/83 ENSO induced global precipitation anomalies might be a main driver for the rapid expansion of dry land (Dai and Wigley 2000). As the computation of PDSI heavily relies on the temperature and precipitation, it is reasonable that global scale climate anomaly caused trend shifts in PDSI should statistically reflect this kind of transition in turn.

PDSI secular trend and turning point
The secular trend of the averaged PDSI series for all the valid samples showed an accelerated drying trend over land surface since 1950 ( figure 3). This pattern is broadly consistent with that of Dai (2011b). Recent studies also found continued expansions of drylands (defined as regions where the ratio of annual precipitation to annual potential evapotranspiration less than 0.65) since 1950 s with the supports of global climate models in CMIP5 and observational data from the Climate Prediction Center (e.g. Fu 2013, Huang et al 2016). The secular trend of the averaged PDSI series also suggested that, although the averaged PDSI decreased very close as derived from the linear regression and EEMD, the decreasing rate of this drying trend was accelerating and this trend might expect to continue in the near future, at least from current data point of view.
A number of modelling efforts also suggested potential future increase of drought globally. For example, Feng and Fu (2013) and Huang et al (2016) predicted that global drylands will increase by 23% and 11% (or 10%) relative to 1961-1990 baseline respectively, projected under the Representative Concentration Pathways (RCPs) RCP8.5 and RCP4.5 with an ensemble mean of CMIP5 models, by the end of this century. Sheffield and Wood (2008) found that soil moisture decreased and drought occurrence monotonically increased globally; more specifically, the spatial extent of severe soil moisture deficits and frequency of short-term (4-6-month duration) droughts since the mid-20th century will double, while the long-term droughts (>12 months) will become three times more common, with eight IPCC 4th Assessment Report (AR4) models under the Special Report on Emissions Scenarios (SRES) B1, A1B and A2 future climate scenarios and simulations up to the end of 21st century. Dai (2011b) also found widespread drought till the end of 21st century with drought indices, i.e. PDSI_pm and sc_PDSI_pm, derived from multi-model ensemble-mean monthly data from 22 coupled climate models participated in the IPCC AR4 under the SRES A1B. And it is noted that the drought occurrence will increase even under a reduced greenhouse gas emission pathway because of the already accumulated greenhouse gases and the thermal inertia of the oceans (Wigley 2005, Sheffield andWood 2008).
The PDSIs with monotonic secular trends occupied 74% of all the samples, in which the ratio of the monotonic increase samples (38.5%) was very close to the monotonic decrease ones (35.5%) (table 1). Dai (2011b) did a linear regression with sc_PDSI_pm data from 1950 to 2008 to analyze the trend of PDSI globally, and it is intriguing to compare with the secular trends obtained by EEMD, especially the inconsistences of the broad pattern of PDSI trend between Dai's work and this study. Dai's results (figure 7(c) in 2011b) suggested that Alaska and northern Canada (except the Arctic areas), southern Europe, East and South Asia, eastern Australia, and much of Africa were drying, while the United States, western China and most of Western Australia were becoming wetting. In this study, however, we noted many PDSI samples with decrease-increase secular trend form appeared in Alaska, northern Africa, and East Asia (including eastern Russia); in the western United States, a number of PDSI samples showed increase-decrease secular trend form; in northern South America, large part of the PDSI samples showed wetting trend, and it is similar for the southern Africa; in southern Europe, the drying PDSI samples mixed with the wetting ones; the regions around Ukraine were dominated by PDSI samples with increase-decrease secular trend; in Eurasia, the wetting PDSI samples were more than the drying ones; and the northeast China was becoming wetting (figure 4).
The spatio-temporal distribution of the secular trend types (figure 4), together with the turning point, presents some interesting points that might be masked by the linear regression method. The spatiotemporal distribution of the turning points showed obvious aggregation tendency and could provide hints about the long-term changes of the regional climate systems (figure 5). For example, the western United States was becoming drier partially (Gutzler andRobbins 2011, Jiang et al 2019); the northeast China was turning wet while southern China was becoming drier (Zhai et al 2005); and part of the Sahel region was becoming drier despite the overall greenup trend in that area (Hickler et al 2005).

Spatio-temporal patterns of PDSI secular trend
The zonal average of PDSI secular trend showed that most of the land surface were drying on average over time (figure 6); although large areas land surface with wetting trend existed in the meantime, its amplitude was relative weak compared with the drying counterpart ( figure 7). It is noted that the secular trend of the global monthly averaged PDSI series also showed an accelerated drying trend (figure 3).
Two bands of wetting trend were observed in the high latitudes in the northern hemisphere, in which the band around ∼70 • -75 • N emerged in 1950s, earlier than the band around 60 • N, but it seems that the extension of the tongue-like wetting trend in the band around 60 • N was strengthening; in contrast, the northern tropical band was dominated by a relative strong drying trend. The amplitude of this globalized drying intensified since 1950 and some spatial features could be inferred from the accumulated PDSI secular trend change maps (figure 7). The drying in the northern tropical area was dominated by the southern Sahel region, while the drying trends in the middle and high latitudes in the northern hemisphere were mainly determined by the eastern Asia and the northern North America. It is also notable that some wetting hotspots largely enhanced after 2000s (figure 7(f)), for example the eastern Russia and central Asia. Because the land surface area above 40 • S in the southern hemisphere is relative small, the spatio-temporal patterns of PDSI secular trend in those areas will not be discussed here.
The maps of the PDSI secular trend change rate per decade provide an interesting perspective about the evolution of drought globally (figure 8). The most highlighted feature is that it seems a periodic turnover of wetting/drying hotspots existed at global scale. The broad pattern of the PDSI secular trend change rate per decade in 1960-1950 (figure 8(a)) is much like the inversion of that of 2015-2000 (figure 8(f)), while the middle periods provide transition to finish this process in about 50 years or so. This multidecadal variability implies that the spatio-temporal patterns of global drought might have relationship with the periodic variations of the oceanic and atmospheric current ranging from 50 to 70 years (Delworth and Mann 2000, Mantua and Hare 2002, Semenov et al 2010.

Limitations and prospects
The PDSI, as well as the Standardized Precipitation Index (SPI) (Guttman 1998) and the Standardized Precipitation-Evapotranspiration Index (SPIE) (Vicente-Serrano et al 2010), are widely used indices to characterize meteorological drought, and it is not out of expectation that different drought indices may characterize droughts in varied ways. For example, there is discrepancy in the spatial distribution of the meteorological drought trend in this study and Spinoni et al (2019), in which the global scale SPEI and SPI dataset were analyzed, and this discrepancy obviously originates from both the trend extraction methods applied and the essential differences in SPI/SPEI and sc_PDSI_pm (Trenberth et al 2013). The results of this study have some implications about the historical secular trend of global meteorological drought from a new methodological perspective of EEMD, but attention should be taken when extrapolating the related findings regard of the limitations of the specific PDSI dataset (Trenberth et al 2013).
The primary sources of uncertainties stem from temporal and spatial inhomogeneities of the global fields of the forcing data, e.g. precipitation, solar radiation, surface humidity, and wind speed, required to calculate the PET in PDSI and the corresponding calibration period (Dai and Zhao 2017). In addition, the ET estimation may also contain substantial uncertainties (Wang and Dickinson 2012). However, despite the uncertainties in meteorological forcing data, the sc_PDSI_pm dataset used in this study had showed significant correlation with observed soil moisture, streamflow, and water storage changes over large areas globally, including the cold regions (Dai 2013). As argued by Trenberth et al (2013), in which discrepancies in global drought trend separately obtained by Dai (2013) and Sheffield et al (2012) using different PDSI_pm products were thoroughly analyzed, the authors suggested that Dai's sc_PDSI_pm dataset adopted a more reasonable base period to calibrate the PDSI moisture categories, and used a more reliable precipitation data set (Dai and Zhao 2017).
The primary sources of uncertainties stem from temporal and spatial inhomogeneities of the global fields of the forcing data, e.g. precipitation, solar radiation, surface humidity, and wind speed, required to calculate the PET in PDSI and the corresponding calibration period (Dai and Zhao 2017). In addition, the ET estimation may also contain substantial uncertainties (Wang and Dickinson 2012). However, despite the uncertainties in meteorological forcing data, the sc_PDSI_pm dataset used in this study had showed significant correlation with observed soil moisture, streamflow, and water storage changes over large areas globally, including the cold regions (Dai 2013). As argued by Trenberth et al (2013), in which discrepancies in global drought trend separately obtained by Dai (2013) and Sheffield et al (2012) using different PDSI_pm products were thoroughly analyzed, the authors suggested that Dai's sc_PDSI_pm dataset adopted a more reasonable base period to calibrate the PDSI moisture categories, and used a more reliable precipitation data set (Dai and Zhao 2017).
Although the PDSI model contains uncertainties and limitations (Alley 1984), it is still a fairly comprehensive index of relative drought by synthesizing precipitation, ET, runoff, and soil properties (Trenberth et al 2013). To diminish the discrepancies in different studies and the uncertainties in PDSI products, high quality/consistent meteorological forcing data, including the accurate estimation of ET (Wang and Dickinson 2012), is still of critical priority and is the basis for comparison between different studies at regional or global scales. Further, as suggested by Sheffield et al (2012), an alternative solution to above mentioned problem is to use physically realistic hydrological modelling driving by in situ and satellitebased data, although longer historical estimation cannot be established.

Conclusions
This study applied EEMD, a method suitable for nonstationary and nonlinear processes, analyzed the spatio-temporal characteristics of the secular trends of a monthly sc_PDSI_pm product per pixel over global land surface from 1950 to 2015. We compared the EEMD with two non-parametric methods, i.e. Mann-Kendal and Pettitt's tests, and BFAST in extracting the trend and turning point of PDSI series. We found that there is large discrepancy in the secular trend extracted by EEMD and the trend types detected by the MK. The probability distributions of the turning points position (period) in the secular trend extracted by EEMD is much close to a normal distribution than that of the Pettitt's test and BFAST, and is physically meaningful if taking global scale climate anomaly background into consideration. There was an accelerated drying trend over global land surface as a whole, but large areas with wetting trend featured by relative weak amplitude still existed in the meantime, especially at the high latitudes in the northern hemisphere. We demonstrated that the secular trend of PDSI series extracted by EEMD could provide more detailed spatio-temporal characteristics, represented by the shifts of trend and nonlinear property of the secular trend, of global drought than that of the nonparametric or linear regression methods. The adaptively extracted secular trend of PDSI by EEMD could present more insights about the transition and progress of wetting and drying trend over global land surface at multidecadal scale. Future studies could further improve the present analysis by using long-term land surface observational or model simulated data, e.g. soil moisture and streamflow data, to validate the drought trends derived from PDSI, and it would be of particular interest to study the response of vegetation growth to drying/wetting at regional or global scales using satellite-derived vegetation growth data.