2.1 Precipitation data

Monthly precipitation data from 59 sites located within the state of Zacatecas, Mexico (Fig. 1) were analyzed. Table I provides information on location. Data sets were kindly provided by the Comisión Nacional del Agua, the official national institution in charge of recording climatic and meteorological data in Mexico. Time series of 34 years (1976-2010) were checked for outliers or missing values.

Fig. 1 Study region of the Mexican Zacatecas state territory. 

2.2 Precipitation linear trends and anomaly time series

Precipitation time series are in general non-stationary and frequently present long-term trends. Identifying the trends and detrending the data are both of great interest and importance in data analysis. Moreover, detrending climate data is often a necessary step before meaningful spectral results can be obtained (^{Wu et al., 2007}). Linear trends of the 59 precipitation time series were estimated through a linear regression analysis, considering this simple model:

where Y _i is the i _th monthly precipitation and X _i is the i _th month.

Before carrying out the fractal analysis with the power spectrum density approach, linear trends of the series were removed following Eq. (2):

where Y _di is the i _th detrended monthly precipitation.

After removing the trends, the standardized anomalies were computed simply by subtracting the sample mean and dividing by the corresponding sample standard deviation, as pointed out by ^{Wilks (2006)}. The new time series were then referred to as “SAMP time series” and could be analyzed by using common statistical analysis and other approaches like the spectral analysis which uses several algorithms such as maximum entropy, autoregressive moving average, non-harmonic method, singular spectral analysis and the traditional Fourier analysis (^{Lana et al., 2005}).

2.3 Cluster analysis

Cluster analysis is an effective statistical tool to find homogeneous climate regions based upon observed values of meteorological parameters (^{Unal et al., 2003}; ^{Mahlstein and Knutti 2010}; ^{Karmalkar et al., 2011}). This approach groups data objects based only on information found in the data that describes the objects and their relationships. The goal is that the objects within a group be similar to one another and different from the objects in other groups (^{Tan et al., 2005}). The greater the similarity within a group and the greater the difference between groups, the better or more distinct the clustering becomes (^{Tan et al., 2005}). In this study, we used a tree clustering algorithm to join objects (i.e., SAMP) into successively larger clusters (i.e., regions under the basis of similar precipitation regimes). The 1-Pearson’s r (^{Everitt et al., 2011}) was used as linkage distance and the Ward’s method was used as linkage rule.

The joining (or tree clustering) method uses the similarities or distances between objects when forming the clusters. Similarities are a set of rules that serve as criteria for grouping or separating items. The distances can be based on a single dimension or on multiple dimensions, with each dimension representing a rule or condition for grouping objects. In the case of the 1-Pearson’s r (^{Everitt et al., 2011}) as linkage distance, low values mean that both involved elements are highly positively correlated, while two elements having low positive correlation are more distant from each other. The ^{Ward’s method (Ward, 1963)} is generally regarded as very efficient although it tends to create small-sized clusters (^{Unal et al., 2003}) and uses an analysis of variance approach to evaluate the distances between clusters. This method defines the proximity or distance between two clusters as the increase in the squared error that results when two clusters are merged (^{Tan et al., 2005}). The cluster analysis was performed in this research using Statistica v. 8.0 (^{StatSoft, 2007}).

2.4 Fractal analysis

The temporal variation of natural phenomena has been difficult to characterize and quantify, which is why ^{Mandelbrot (1982)} introduced the fractal analysis. Time series can be characterized by a non-integer dimension (fractal dimension) when treated as random walks or self-affine profiles. These series can be treated as self-affine systems because they are invariant only under anisotropic magnifications (^{Bunde and Havlin, 2013}). Self-affine systems are often characterized by roughness, which is defined as the fluctuation of the height over a length scale. For self-affine profiles, the roughness scales with the linear size of the surface by an exponent called the Hurst exponent (^{Turcotte, 1997}). This exponent, however, only gives limited information about the underlying distribution of height differences (^{Evertsz and Berkner, 1995}). The Hurst exponent as well as the fractal dimension measure how far a fractal curve is from any smooth function used to approximate it (^{Moreira et al., 1994}). There are several approaches to estimate fractal dimension for self-affine profiles but we only used the power spectrum technique because it is very sensitive and is a good exploratory tool for real data (^{Weber and Talkner, 2001}).

2.5 Power spectrum approach

Self-affine fractals are generally treated quantitatively through spectral techniques. The variation of the power spectrum P(f) with frequency f appears to follow a power-law (^{Turcotte 1997}):

The power spectrum P(f) is defined as the square of the magnitude of the Fourier transform of the SAMP. Denoting SAMP as a function of time by Z(t), we have:

where t ₀ and t ₁ are the limits of time over which the series extend. The discrete version of Eq. (4) should be used for such a standardized anomaly sampled at discrete time intervals:

The next step is to obtain a relationship between the power β and the fractal dimension D. After considering the time series Z ₁ (t) and Z ₂ (t) separated by distance r and related by

it can be observed that Z ₁ (t) has the same statistical properties as Z ₂ (t), and since Z ₂ is a properly rescaled version of Z ₁ , their power spectral densities must also be properly scaled. We can then write:

It follows that

And

where D _s denotes the fractal dimension estimated from the power spectrum and H is the Hurst exponent (see ^{Valdez-Cepeda et al., 2003a}, ^b, ²⁰⁰⁷).

Since the periodogram is a poor estimate of the power spectrum, we used the Welch method of averaging periodograms (^{Welch, 1967}). We also used the “running sum” transformation to shift the slope by a factor of +2. The Hurst exponent and the D _s (because of data-trace) had a slope between -1 and 1 on the log-log plot. Major modal contributors could then be identified in such a plot.

In order to obtain an estimate of the fractal dimension D _s , one must calculate the power spectrum P(f) (where f, f is the wavenumber, and l is the wavelength), and plot the logarithm of P(f) versus the logarithms of f. This plot should follow a straight line with a negative slope -ß if the profile is self-affine (^{Valdez-Cepeda et al., 2003a}, ^b).

The power spectrum approach involving the SAMP was carried out using the Benoit 1.3 software (^{TruSoft, 1997}). The important frequencies of SAMP were estimated through the plot yielded by power spectral density Φ _x (f) versus frequency, considering significant p ≤ 0.05 peaks.

2.6 Multivariate El Niño Southern Oscillation Index data

The ENSO is arguably the most important coupled ocean-atmosphere phenomenon causing global climate variability on inter-annual time scales (^{Wolter and Timlin, 2011}). The MEI is calculated based on the six main observed variables over the tropical Pacific (^{Wolter 1987}; ^{Wolter and Timlin 1993}, ¹⁹⁹⁸). These six variables are sea-level pressure, zonal and meridional components of the surface wind, sea surface temperature, surface air temperature, and total sky cloudiness fraction. The MEI is computed separately for each of 12 sliding bi-monthly seasons (Dec/Jan, Jan/Feb, …, and Nov/Dec). After spatially filtering the individual fields into clusters (^{Wolter 1987}), the MEI is calculated as the first un-rotated principal component (PC) of all six observed fields combined. This is accomplished by normalizing the total variance of each field first and then performing the extraction of the first PC on the covariance matrix of the combined fields (^{Wolter and Timlin, 1993}). All seasonal values are standardized with respect to each season in order to keep the MEI comparable. A MEI database with a length matching our precipitation records (1976-2010) was downloaded from the National Oceanic and Atmospheric Administration website (http://www.esrl.noaa.gov) and used to couple it with each of the standardized anomaly regional means of the monthly precipitation time series. The MEI was used because it has been related to precipitation in several regions around the world and takes into account ENSO’s seasonality. For instance, precipitation was positively and significantly correlated with the MEI during the boreal winter (Nov-Apr), mainly in northeastern and northwestern Mexico (^{Méndez-González et al., 2009}).

2.7 Wavelet coherence analysis

Wavelet analysis is a spectral method of decomposing a time series into time and frequency space. This allows the identification and analysis of dominant localized variations of power, i.e., where the variance of the time series is largest for a given frequency (^{Keener et al., 2010}). Transferring the concept of the Fourier analysis, a wavelet power spectrum (WPS) can be defined as the wavelet transformation of the autocorrelation function (^{Maraun and Kurths, 2004}). This can be implemented following the Wiener-Khinchin theorem (^{Maraun and Kurths, 2004}):

where the brackets denote expectation values and WPS _i (s) describes the power of the signal x(t) at a certain time t _i on a scale s.

As in the Fourier analysis, the univariate WPS can be extended to compare time series, for instance, x(t _i ) and y(t _i ). One can define the wavelet cross-spectrum [WCS _i (s)] as the expectation value of the product of the corresponding W _ix (s) and W _iy (s):

In contrast to the WPS, the WCS is a complex value (analogous to Fourier cross-spectra) and can be decomposed into amplitude |WCS _i (s)| and phase Φ_i(s):

The phase Φ _i (s) describes the delay between the two signals at time t _i on a scale s.

A normalized time and scale resolved measure for the relationship between the time series x(t _i ) and y(t _i ) is the wavelet coherency (WCO), defined as the amplitude of the WCS normalized to two single WPS, is the following:

A value of 1 means a linear relationship between series x(t _i ) and y(t _i ) around time t _i on a scale s. A value of zero is obtained for vanishing correlation.

Spurious peaks can occur for areas of low wavelet power due to normalization. Note the importance of estimating the expectation values of the single terms in Eq. (14). If they are not considered separately, numerator and denominator will equal each other, and the result will be identical to 1 no matter what the real WCO is. When computing WCO from observational data, one must thus smooth the numerator and denominator separately (^{Maraun and Kurths, 2004}). In the current case, WCO was performed involving the MEI series and each of the three regional SAMP under the basis of a 95% confidence level. We used the MatLab code developed by ^{Grinsted et al. (2004)} for this purpose. The code is available at http://noc.ac.uk/marine-data-products/cross-wavelet-wavelet-coherence-toolbox-matlab.

3.1 Cluster analysis

A cluster analysis was performed using the 59 SAMP time series with the aim of grouping them into successively larger groups (i.e., regions under the basis of similar precipitation regimes). The dendrogram (Fig. 2) revealed three clusters or regions (Table I) when using a linkage distance value of 1: (I) semi-arid region (11 out of 59 stations); (II) highlands region (27 out of 59 stations), and (III) canyons region (21 out of 59 stations).

Fig. 2 Dendrogram for 59 standardized anomalies of the monthly precipitation through the Zacatecas state territory using the amalgamation rule of the Ward’s method and the distance metric of 1-Pearson’s r. 

A map showing the stations and their corresponding regions is presented in Fig. 3. Those regions are distributed from NE to SW of the state territory, that is, from the driest to the wettest areas. As expected, cluster analysis results agree with the spatial distribution of the annual mean precipitation throughout the Zacatecas state territory (Fig. 4). Thus, cluster I (semi-arid region) includes weather stations with annual precipitation from < 300 to 400 mm, cluster II (highlands region) includes stations with annual precipitation from 400 to 500 mm, and cluster III (canyons region) includes stations with annual precipitation from 500 to > 700 mm.

Fig. 3 Location of the weather stations and their corresponding regions found through the cluster analysis in the Zacatecas state territory. 

Fig. 4 Distribution of annual mean precipitation along the Zacatecas state territory. 

The regions identified in this research agree well with the three floristic districts reported by ^{Balleza et al. (2005)} on the basis of the geographic distribution of 456 species of the Asteraceae family. Considering results from both studies, the semi-arid region covers the floristic district named Ia; the Highlands region agrees with the floristic district named Ib; and the canyons region covers the floristic district named II. This implies Asteraceae family species richness and their distribution could depend on the rainfall patterns as suggested by our results. In this context, ^{Balleza et al. (2005)} results suggest that the species of the Asteraceae family distributed in the Zacatecas territory can be divided into two main groups: (i) those thriving mostly in the dry environments (semi-arid region), and (ii) those thriving in the temperate, warm habitats along the highlands and canyons regions. In fact, the flowering plant family Asteraceae is one of the most influential in the diversification and evolution of many animals that heavily depend on their inflorescences to survive (^{Barreda et al., 2015}). Species of the Asteraceae family are found almost everywhere in the world except in Antarctica (^{Barreda et al., 2015}), which means that their climatic requirements are very diverse, although it seems that climatic factors restrict the phenological period of most species as pointed out by ^{Torres and Galetto (2011)}.

Rainfall patterns define the agricultural activities performed by peasants and growers of each region. Such activities represent an important economic income to most of the families settled throughout the Zacatecas territory. In the semi-arid region, goats are kept for dairy and meat production although those activities represent only low income at a family level; the major production systems are characterized by small flocks (< 100 animals), feeding based on grazing and browsing native vegetation, following regular routes and keeping the animals overnight in rudimentary shelters. The main products of these systems are kid goats (“cabritos”) and culled goats (^{Tovar-Luna, 2009}). In the highlands region, the main crop is dry beans (Phaseolus vulgaris), followed by species for forage such as corn (Zea mays) and oats (Avena sativa), and horticultural species such as peppers (Capsicum annum), onions (Allium cepa) and garlic (Allium sativum). In the canyons region, the main activity refers to the cow-calf production system followed by the production of guava (Psidium guajava), and a great variety of horticultural products.

3.2 Power spectrum analysis

A spectral analysis was performed on each regional SAMP series to estimate their self-affinity statistics (Table II): mean, median, standard deviation (Std), fractal dimension (D _s ), Hurst exponent (H) and slope (-β). All three regional SAMP series yielded straight lines on the log-log plot with -β, suggesting that P(f) ~ f ^-b (^{Valdez-Cepeda et al., 2003a}, ^b). This means that each spectrum is singular and is represented by a curve on the complex plane for the SAMP series of the case.

Spectral D _s values vary from 1.354 to 1.527, and H values change from 0.473 to 0.646. The D _s and H values of 1.354 and 0.646, respectively, suggest that the semi-arid region SAMP series tends to show persistent behavior (D _s < 1.5 and H > 0.5). The highlands region SAMP series has D _s and H values of 1.497 and 0.503, respectively, while the canyons region SAMP series has D _s and H values of 1.527 and 0.473, respectively. These results indicate that the SAMP series in both regions tend to be like the Brownian movement or random walk (D _s ≈ 1.5), as pointed out by ^{Mandelbrot (1985)}.

3.3 Important frequencies

Plots of power spectrum density [Φ _x (f)] against frequency allow us to identify the dominant frequencies on the regional SAMP series in terms of the dominant frequencies (1 year^-1) that are likely to be important to the precipitation process (Fig. 5a-c). This approach gives us some frequency components that do not consider time and length because the analysis gives us a resolution in a frequency that is determined for the window size over the analyzed time series. In other words, the results give us useful information about the frequency content of the analyzed series, but they do not indicate at which time these frequencies occur. The detected frequencies were those from confidence bands of 95% in the adjusted linear model and are likely to be important to the precipitation process. Significant peaks of the power spectrum were detected at p ≤ 0.05.

Fig. 5 Power spectrum plots of the monthly mean standardized anomalies of the monthly precipitation time series for (a) semi-arid region, (b) highlands region and (c) canyons region of Zacatecas, Mexico from January 1976 to December 2010, showing important frequencies in years. The power spectral density is given as a function of frequency for time scales of 2 to 408 months. The fitted line was used to estimate the fractal dimension (Ds). Dashed lines indicate the confidence interval at p ≤ 0.05. 

The 0.5-, 1-, 2- and 3-year frequencies are present in the power spectrums of all analyzed regional SAMP series. The highlands and canyons regions showed important frequencies of 13 and 12 years, respectively.

The 0.5- and 1-year frequencies may be linked to the seasonal cycle while the 2-year frequency could be related to the Quasi-Biennial Oscillation (^{Kane, 2009}), which is a 26-month cycle explaining solar activity and wind reversal in the lower stratosphere of the North Pole (^{Labitzke and Loon, 1989}; ^{Mendoza et al., 2001}). The 3-year frequency suggests a possible effect of a quasi-periodic event like ENSO. It is widely recognized that ENSO has an erratic cycle from 3 to 5 years (^{Weber and Talkner, 2001}), 3 to 6 years (^{Monetti et al., 2003}), 4.3 years (^{Luo et al., 2015}), or 2 to 7 years (^{Zubair, 2002}; ^{MacMynowski and Tziperman, 2008}). These discrepancies may be explained by the physical dependence of ENSO frequency on the structure of the coupled mode provided by diagnosing the relative contributions of the thermocline feedback and zonal advection feedback on ENSO evolution (^{An and Wang, 2000}). The 12- and 13-year frequencies indicate that there may be effects from the sunspot cycle, which varies from 8 to 14 years (Mendoza et al., 2001) or 11 years (^{Hu et al., 2009}) with a long-term average of 11.3 years. We then carried out a wavelet coherence analysis on each of the SAMP and the MEI time series, each with a significant 5.5-year frequency, to identify possible relationships.

3.4 Wavelet coherence analysis

Results of phase and wavelet coherence between each regional SAMP series and the Multivariate ENSO Index series are shown in Figure 6a, b, and c, which correspond to the semi-arid, highlands and canyons regions, respectively. The black figure outlines indicate areas significant to 95% confidence, that is, cones of influence; spots limited by a thick black line point out a possibly significant relation between both phenomena; and vectors indicate the phase difference between the regional SAMP and the MEI series with arrows pointing clockwise to indicate in-phase behavior and arrows pointing counter-clockwise to indicate anti-phase behavior (^{Keener et al., 2010}).

Fig. 6 Wavelet coherence and phase between standardized anomalies of the monthly precipitation (SAMP) of (a) semi-arid region, (b) highlands region, and (c) canyons region (c), and the Multivariate ENSO Index (MEI) for the Zacatecas state. Vectors indicate the phase difference between the SAMP and MEI. Thick black figure outlines indicate areas significant at a 95% confidence level using the red noise model, and the thin black line surrounds the cone of influence.