Characteristics of the study area

The state of Zacatecas is located between the coordinates 21º 01’ 00” and 25º 07’ 00” N and 100º 43’ 00” and 104º 22’ 00” W. Its area is 75 284 km ². The main climate types are arid and semiarid (73%) and temperate subhumid (17%). Mean annual temperature is 17ºC and average annual high and low temperatures are 30ºC (May) and 3ºC (January), respectively. The mean annual precipitation in Zacatecas is 510 mm, with 300 mm in the north and 860 mm in the south. Most of the rainfall occurs in summer (75%), from June to September. The climate of Zacatecas limits agriculture, the main economic activity to which 1.74 million hectares is dedicated; 88% of this area is rainfed.

Climatological information

To identify the drought events, as well as to determine their duration, magnitude, intensities, and frequencies, this study used daily precipitation and temperature data series from 9 weather stations in the state of Zacatecas (Figure 1), taken from the CICLOM database (^{Conagua, 2014}) for the period from January 1, 1961, to December 31, 2012. The geographic descriptors (coordinates and altitude) of the weather stations and the mean annual values corresponding to precipitation, P, and temperature, T, are presented in Table 1.

Figure 1 Weather stations selected for drought analysis. 

Selection and quality control of the climatological information

A key stage in the analysis of hydroclimatological time series, as in the case of drought analysis, is the data selection and quality control. Selection of the climatological information used is primarily based on four criteria: (1) that the weather station is active, to be able to evaluate the droughts during the last few years; (2) that each weather station has precipitation and temperature records covering at least 40 years; (3) that the percentage of missing data is less than or equal to 15% in order to maintain data representativeness; and (4) that the stations are distributed throughout the study zone so as to have good spatial representation of the drought.

For drought analysis, a preliminary quality control of the precipitation and temperature data series was carried out, checking for homogeneity of the climatological data series and filling missing data. To check for homogeneity of precipitation and temperature data, parametric and non-parametric statistical tests were applied to the annual series. The tests selected were the Student’s t-test, the Cramer test (^{Mirza, 1997}), the Standard Normal Homogeneity Test (SNHT), and the Buishand test (^{Hänsel, Mederios, Matschullat, Petta, & Mendoça-Silva, 2016}).

To select the method for filling in missing data, several methods were tested: polynomial interpolation, arithmetic average, inverse distance, and normal ratio (^{Tabios III & Salas, 1985}; ^{Dingman, 2002}). The normal ratio had the best results (Equation (1)):

where p^m is the missing datum to be estimated for the day, month, or year in the station of analysis m; g is the neighboring stations; ng is the number of neighboring stations considered in the analysis; p-m and p-g are the annual mean precipitation at the station m with missing information and at the neighboring stations g, respectively; and pg is the precipitation registered at the neighboring stations on the day, month, or year in which the datum at station m was missing. This same method was also applied to the temperature data.

Selection of the scales of drought analysis (k)

The time scale, k, plays an important role in the study of drought characteristics. Short time scales are relevant because they are representative of meteorological and agricultural droughts. A time scale of 3 months gives a seasonal estimation of precipitation that reflects the moisture conditions on the short- and medium-term. A 6-month scale estimates precipitation on a medium-term scale and can be associated with reservoir levels and anomalous flows, while a time scale longer than 12 months reflects long-term patterns in precipitation that are likely related to flows, reservoir levels, and even groundwater levels, and permits the identification of hydrological drought.

This study selected scales of 3, 6, and 12 months for the SPI, RDI, and SPEI indices (SPI-3, SPI-6, SPI-12, RDI-3, etc.), with the aim of obtaining short-, medium-, and relatively long-term perspectives of drought characteristics.

Calculation of potential evapotranspiration (ETP)

To calculate the RDI and SPEI indices, potential evapotranspiration was estimated using the ^{Thornthwaite (1948)} method, as proposed by ^{Tsakiris and Vangelis (2005)} and ^{Vicente-Serrano et al. (2010)}. This method of calculating ETP (Equation (2)) is one of the simplest, since it only requires mean monthly precipitation data and the latitude of the site where ETP is to be estimated.

where ETP is the monthly potential evapotranspiration (mm/month); T is the monthly mean temperature (ºC), and I is the heat index (Equation (3)), which is calculated as the sum of 12 monthly index values.

where NDM is the number of days of the month and N is the maximum number of hours of sunlight (Equation (6)).

and ws is the hourly sunrise angle (Equation (7)).

In which φ is the latitude (radians) of the weather station and δ is the solar declination (radians) calculated as:

where J is the average Julian day of the month.

The Standardized Precipitation Index (SPI)

The SPI uses precipitation as the only variable for the analysis, and estimates whether, in a given region and period, there is a deficit or excess of precipitation relative to normal conditions (^{Hayes, Svoboda, Wilhite, & Vanyarkho, 1999}). The calculation of the SPI is based on long-term precipitation records. ^{Guttman (1994)} points out that preferably 50 to 60 years or more of records are needed, although other researchers (^{Wu, Hayes, Wilhite, & Svoboda, 2005}; ^{Li et al., 2014}) recommend at least 30 years. The SPI is calculated by fitting an appropriate probability distribution function (pdf), generally the gamma distribution function, to cumulative monthly precipitation data for each time scale, k, and site of interest. The SPI values are obtained by transforming the gamma distribution value into normal standard values. The gamma probability density function is defined by:

where a is the shape parameter; b is the scale parameter; x is the cumulative monthly precipitation (x>0), in mm, and Гa is the gamma function. Parameters a and b can be estimated using several methods, and their values by the maximum likelihood method are (^{Haan, 1977}):

where A is an auxiliary variable that is defined by:

where n' is the number of precipitation data xi different from zero and x- is the mean value of the xi data.

Integrating the gamma density function (Equation (9)), relative to x, we obtain the cumulative distribution function Gx

Considering t=x/b, Equation (13) is then reduced to the incomplete gamma function.

Given that the gamma function is not defined for x=0, and the precipitation series can have zero values, then the cumulative probability can be calculated as (^{Angelidis, Maris, Kotsovinos, & Hrissanthou, 2012}):

where Hx is a function of the cumulative probability; q is the probability of having precipitation values equal to zero; and Gx is the cumulative probability of the incomplete gamma function.

If m is the number of zeros in the precipitation series, then q=m/n, where n is the total number of data in the precipitation series, according to the scale of analysis (k). When the series does not have precipitation values equal to zero, then q=0, and thus Hx=Gx, and when x=0, then H0=q.

Finally, to define the SPI value, the cumulative probability function, Hx, is transformed into the standard normal variable, Z, which has a mean of zero and a variance of one (^{Edwards & McKee, 1997}), using the approach proposed by ^{Zelen and Severo (1965)}:

where c0=2.515517, c1=0.802853, c2=0.010328, d1=1.432788, d2=0.189269, and d3=0.001308.

A drought event occurs when, at any time scale, SPI is continuously negative and reaches a value of -1.0 or less. It ends when the SPI becomes positive (^{McKee et al., 1993}; ^{Vrochidou & Tsanis, 2012}). The SPI values can be easily compared simultaneously on the spatial and temporal dimensions (^{Lopez-Bustins, Pascual, Pla, & Retana, 2013}). Based on the range of SPI values, drought events are classified from mild to extreme (Table 2).

The Reconnaissance Drought Index (RDI)

The RDI can be calculated at different time scales and its interpretation is similar to that of the SPI (Table 2) (^{Vicente-Serrano, Van der Schrier, Beguería, Azorin-Molina, & López-Moreno, 2015}; ^{Xu et al., 2015}). One of the most notable theoretical limitations of this drought index is that it is not valid when the ETP value is equal to zero, a very common situation in cold regions during winter (^{Vicente-Serrano et al., 2015}). RDI is expressed in three ways: the initial value (αk), normalized RDI (RDIn), and the standardized RDI (RDIst). The initial value (αk) is presented in aggregated form using a monthly time scale and can be calculated on a monthly, seasonal, or annual basis (^{Asadi-Zarch, Malekinezhad, Mobin, Dastorani, & Kousari, 2011}). Its values are calculated with the following equation:

where Pij and ETPij are the precipitation and potential evapotranspiration for month j of year i, respectively, in mm; Na is the number of years of available information; and k is the scale of analysis at which the RDI is calculated. The normalized and standardized forms of the RDI are (^{Tsakiris, Pangalou, & Vangelis, 2007}):

where δ-k(i) is the arithmetic mean of δk(i) and yk(i) are the natural logarithms of δki, whose mean and standard deviation are y-k and σ^yk, respectively.

Different studies with data from different places and different time scales have found that, in all cases, the values of δk followed either a gamma or a lognormal distribution, although in most of the cases the gamma distribution was more successful ^{(Kousari et al., 2014}; ^{Cai, Zhang, Yao, & Chen, 2015}).

^{Tigkas (2008)} and ^{Asadi-Zarch et al. (2011)} have shown that the RDIst can be better calculated by fitting the gamma distribution to the distribution of δk frequencies, since this method tends to solve the problem of calculating RDIst for short time scales, which can include values of zero precipitation (δk=0), for which Equation (20) is not defined.

In our study, the RDI was obtained by fitting the gamma distribution to δk values, in a way similar to SPI, using Equations (10) through (17), and substituting the values of δk instead of x.

The Standardized Precipitation-Evapotranspiration Index

The main limitation of the SPI is that it is based only on precipitation and ignores other variables that affect atmospheric water demand, such as temperature, wind speed, solar radiation, and vapor pressure deficit (^{McEvoy, Huntington, Abatzoglou, & Edwards, 2012}). To overcome this limitation, ^{Vicente-Serrano et al. (2010)} developed the SPEI, which combines the PDSI’s sensitivity to changes in the evaporation demand, which is caused by temperature fluctuations and tendencies, with the SPI’s simplicity of calculation and multi-scale nature (^{Vicente-Serrano et al., 2010}; ^{Banimahd & Khalili, 2013}; ^{Li et al., 2014}).

The SPEI is based on the monthly climatic water balance (precipitation minus potential evapotranspiration, called series D), which is calculated at different time scales and is mathematically similar to the SPI, but includes the role of potential evapotranspiration (^{López-Moreno et al., 2013}). Its calculation follows an approach similar to the one used to calculate the SPI, but with a three-parameter log-logistic probability distribution instead of the two-parameter gamma distribution (Figure 2). The SPI drought classification can also be used to evaluate the SPEI (Table 2).

As a first step to calculate the SPEI, the difference between precipitation and potential evapotranspiration for month j of year i (series Dij) is calculated, in mm. The series Dij is then added in the time scale (k) of interest to obtain the series Dkij, in millimeters.

where Pij and ETPij have the same meaning as in Equation (18).

To calculate the SPEI, a three-parameter probability distribution needs to be used, since for a two-parameter distribution, the variable x has a lower limit of zero 0<x<∞, whereas for the three-parameter distribution, x can take values in the range (γ<x<∞), where γ is the origin parameter of the distribution. For this reason, ^{Vicente-Serrano et al. (2010)} proposed using the log-logistic distribution, since the Dk series can have negative values. Its probability density function is:

where α, β, and γ are the shape, scale, and origin parameters, respectively, for the Dk values in the range of γ>Dk<∞.

The log-logistic distribution parameters can be obtained using different methods, of which the L-moments method is the most robust and the easiest to apply (^{Ahmad, Sinclair, & Werritti, 1988}; ^{Singh, Guo, & Yu, 1993}).

where Γ1+1/β is the gamma function of 1+1/β and ws is the moment of high probability of s s=0,1,2:

where xi is the ordered random sample x1<x2,⋯,<xn of the Dk values and n is the sample size.

The probability distribution of the Dk series, according to the log-logistic distribution, is given by (^{Vicente-Serrano et al., 2010}):

Finally, to obtain the SPEI values, Fx values are transformed into the standard normal variable, Z, substituting Fx for Hx in Equations (16) and (17).

Comparison of drought indices

The SPI, RDI, and SPEI indices follow similar methodologies and their values have the same statistical significance; therefore, they are comparable. To compare the drought indices, we used the Pearson (r) and Spearman (rs) coefficients.

Mean drought characteristics

Table 4 presents a summary of the number of drought events and the mean duration of drought periods for short-, medium-, and long-term time scales, determined with the three drought indices evaluated. The mean number of drought episodes decreased with longer time scales.

The average number of drought events at a short-term time scale varied from 36 - 48 months for SPI-3, from 41 - 49 for RDI-3, and from 38 - 52 months for SPEI-3. The mean duration of the drought events ranged from 2.3 to 2.9, 2.2 to 2.7, and 2.0 to 3.2 months for SPI-3, RDI-3, and SPEI-3, respectively. At 70% of the weather stations, the drought event of largest magnitude reached the severe drought category with SPI-3 and RDI-3, while with SPEI-3, it reached the category of extreme drought at 90% of the weather stations. The longest duration of drought periods registered with SPI-3 and RDI-3 ranged from 5 to 8 months, and from 4 to 7 months for SPEI-3.

Over the last two decades, the most intense drought events occurred at the 3-month scale. According to SPI-3 and RDI-3, these events were from February to May of 1999 and from January to May of 2011, and according to SPEI-3 they occurred from May to July of 2011 and from May to June of 1998.

For the 6-month scale, the drought events detected by the three indices were relatively similar. Average droughts ranged from 28 to 41, 26 to 39, and 27 to 39 months for SPI-6, RDI-6, and SPEI-6, respectively. The mean duration of these events ranged from 2.6 to 3.1 months for SPI-6, from 2.7 to 3.2 months for RDI-6, and from 2.9 to 4.1 months for SPEI-6. With SPI-6 and RDI-6, the most intense events reached the category of extreme drought at all the stations; with SPEI-6 only 20% reached this category and the rest of the stations registered severe drought. According to SPEI-6, the longest durations of dry periods were very similar among the stations, ranging from 6 to 9 months, while according to SPI-6 and RDI-6, the range was greater, from 4 to 15 months.

According to SPI-6, RDI-6, and SPEI-6, the northeastern part of Zacatecas had major droughts in the 1970s and 1980s, in terms of magnitude and intensity. During the 1990s and the first decade of 2000, major droughts occurred in the central and southern parts of the state.

At the 6-month scale, RDI-6 detected the longest drought period (15 months), which occurred in the northern region of Zacatecas (station 32052). The longest drought detected by SPI-6 and RDI-6 was at station 32001 (northern Zacatecas), lasting 13 months, from June 1987 to June 1988. SPEI-6 detected a drought for this same period and station, but in two subperiods, one 6 months long (from July to December 1987) and another lasting 3 months (May to July 1988). The longest drought detected with SPEI-6 was 9 months at stations 32032 and 32054, from March to November of 2011, in southern and west-central Zacatecas, respectively. In this same period, SPI-6 and RDI-6 detected a period of drought of only 4 months, from March to June of 2011.

Regarding the 12-month time scale (medium- and long-term), the average number of drought events was similar for the three indices, ranging from 10 to 21 for SPI-12, from 11 to 21 for RDI-12, and from 13 to 21 for SPEI-12. The mean duration of the drought events ranged from 4.8 to 8.9 months for SPI-12, from 4.2 to 9.4 months for RDI-12, and from 4.6 to 9.1 months for SPEI-12. At this scale, with SPI-12 and RDI-12 the most intense events reached the category of extreme drought at 70% of the stations, and 30% severe drought, while with SPEI-12 only 10% of the stations showed extreme drought and the remaining 90% reported severe drought. The longest dry periods were similar with SPI-12 and RDI-12, and ranged from 16 to 40 and from 14 to 40 months, respectively, while the periods were shorter with SPEI-12, ranging from 12 to 27 months.

During the 52 years of climatological records analyzed in the state of Zacatecas, the SPI detected the most intense droughts on the annual time scale, coinciding with the results obtained by ^{Campos-Aranda (2014}, ²⁰¹⁵), who analyzed droughts with SPI and RDI in central Zacatecas. With the 3- and 6-month time scales, the most intense droughts were detected by SPEI and RDI, respectively.

The most intense droughts at the three time scales occurred between 1998 and 2011. RDI and SPEI also identified extreme droughts in 1974, 1979, 1982, and 1989. Figure 4 shows the spatial distribution of the SPI, RDI, and SPEI for the month with the greatest drought at the three time scales analyzed. With the 3-month time scale, SPI-3 and RDI-3 identified April 2011 as the month with the highest drought intensity. These indices have a spatial behavior that is practically equal (Figures 4a and 4b). Most of Zacatecas (82%) suffered severe drought and only 8% extreme drought. The most intense drought detected by SPEI-3 occurred in June 2011, with extreme drought in 59% of the state, mainly in the central part (Figure 4c). At the 6-month scale, April 2011 was also the month of greatest drought intensity according to SPI and RDI; extreme drought occurred in practically the entire state (Figures 4d and 4e). The SPEI-6 identified July 2011 as the month with the greatest drought; most of the territory of Zacatecas underwent extreme drought (49%) and severe drought (39%). Extreme drought occurred mainly in central and northeastern Zacatecas. At the 12-month scale, September 2011 was identified by the three indices as the month with the highest drought intensity. SPI-12 and RDI-12 had similar spatial behavior, identifying around 55% of Zacatecas with extreme drought (Figures 4g and 4h), mainly in the south and west of the state. SPEI-12 identified extreme drought in 33% of Zacatecas, in the same regions as did SPI-12 and RDI-12, but in a smaller area.

Drought frequency

Drought frequency was analyzed using the different drought classes (Table 2) related to the values of the drought indices (SPI, RDI, SPEI). At the 3-month scale, for the drought indices analyzed, the average time (percentage of average frequency) of dry and wet periods was 17.1 and 16.9%, respectively. Normal precipitation conditions were present 66.0% of the time. SPI-3 and RDI-3 showed a higher frequency of moderate droughts (12.5 and 12.2%, respectively) than SPEI-3, while the frequency of severe droughts was lower for SPI-3 and RDI-3 (4.1 and 4.8%, respectively) than for SPEI-3. Of the three indices, SPI-3 and RDI-3 identified a lower number of extreme droughts. Figure 5a shows the frequency percentage of dry and wet periods for SPI-3, RDI-3, and SPEI-3. With SPI-3 and RDI-3, only the west-central region of the study area (station 32054) had extreme droughts (14 and 15 events, respectively), while with SPEI-3, all the stations had extreme drought events, ranging from 3 to 16 events at stations 32052 and 32032, in northern and southern Zacatecas, respectively.

Figure 4 Spatial distribution of the SPI, RDI, and SPEI indices for the month with the most intense drought at scales of 3, 6, and 12 months. 

Figure 5 Frequency distribution of the values of the SPI, RDI, and SPEI indices: a) 3-month scale; b) 6-month scale; and c) 12-month scale for the period 1961-2012. 

At the 6-month scale, according to the three multiscale indices, normal precipitation conditions were present an average of 67.4% of the time. The frequency of dry events was 16.7% and that of wet events was 15.9%. SPEI-6 resulted in a higher frequency of moderate droughts (11.0%) than SPI-6 and RDI-6 (8.7 and 7.7%, respectively). The same occurred with severe droughts (Figure 5b). However, SPI-6 and RDI-6 presented a higher frequency of extreme droughts than SPEI-6: 3.4 and 3.8% respectively. With SPI-6 and RDI-6, 26 and 29 extreme droughts in central Zacatecas were detected (stations 32018 and 32030), while with SPEI-6 the highest number of extreme droughts detected was 15, which occurred in the north (station 32001). Finally, at the 12-month scale, the average frequency percentage of normal, wet, and dry events was 66.6, 16.9, and 16.5%, respectively. The frequency of droughts with SPI-12, RDI-12, and SPEI-12 was 16.3, 15.8, and 17.4%, respectively. The behavior and values of the frequency of the different classes of droughts for SPI-12, RDI-12, and SPEI-12 were similar to those detected by SPI-6, RDI-6, and SPEI-6. And SPEI-12 presented a higher frequency of moderate and severe droughts than SPI-12 and RDI-12. The opposite occurred for extreme droughts, as seen in Figure 5c. With SPI-12 and RDI-12, 23 and 22 extreme drought events were detected at station 32032 (southern Zacatecas), while with SPEI-12, the highest number of extreme drought events was 17, and occurred at station 32001 (northern Zacatecas).

Correlation of the drought indices

Considering the three time scales and the nine stations analyzed, the SPI and RDI were the mostly highly correlated indices, with correlation coefficients (r) ranging from 0.920 to 0.996. SPI and SPEI had the lowest correlation coefficients, which ranged from 0.336 to 0.982. Figure 6 shows the general trend of the relationships between SPI vs RDI, SPI vs SPEI, and RDI vs SPEI at the 3-, 6-, and 12-month time scales, with dispersion diagrams for the series of indices at El Sauz station (32018), located in the central region of Zacatecas.

At the 3-month scale, the correlation coefficients (r) for SPI-3 vs RDI-3 ranged from 0.920 to 0.969, while SPI-3 vs SPEI-3 and RDI-3 vs SPEI-3 ranged from 0.336 to 0.804 and from 0.583 to 0.867, respectively. In general, the correlations for the 6-month scale were slightly higher than for the 3-month scale, ranging from 0.943 to 0.964 for SPI-6 vs RDI-6, from 0.410 to 0.834 for SPI-6 vs SPEI-6, and from 0.643 to 0.906 for RDI-6 vs SPEI-6. At the 12-month scale, the correlation coefficients ranged from 0.952 to 0.996 for SPI-12 vs RDI-12, from 0.926 to 0.982 for SPI-12 vs SPEI-12, and from 0.952 to 0.993 for RDI-12 vs SPEI-12. For the scales and stations analyzed, the smaller the scale, the lower the correlation between indices (Table 5). At the 3- and 6-month scales, the lowest correlations were in the north of the state (station 32052) and the highest in the south (station 32032).

The Spearman correlation coefficients (Table 5) were similar to those of the Pearson correlation coefficients. As for coefficient r, at the 3- and 6-month scales the SPI vs RDI indices had higher values of rs, while SPI vs SPEI had lower values. For the 12-month scale, the RDI-12 vs SPEI-12 and the SPI-12 vs RDI-12 had the highest rs coefficients. The similarity between the r and rs coefficients increased as the scale increased.

Figure 6 Dispersion diagrams for the SPI vs RDI, SPI vs SPEI, and RDI vs SPEI series for the El Sauz station (32018), period 1961-2012. 

In this study, SPI and RDI were more suitable for characterizing meteorological drought. According to ^{Gocic and Trajkovic (2014)} and ^{Xu et al. (2015)}, RDI and SPI are more suitable indices than SPEI for characterizing drought in arid and semiarid regions. Our results coincide.

Moreover, selection of the probability distribution function for the analysis of drought indices may affect the results, underestimating or overestimating drought events. Therefore, it is helpful to use a frequency analysis to determine the distribution function that best fits the precipitation, precipitation/evapotranspiration, and precipitation - evapotranspiration series when evaluating the SPI, RDI, and SPEI indices, respectively.