## Abstract

The spectral analysis of waterhead measurements on a network of piezometers has allowed the analysis of the spatio-temporal variability of data in the Vega de Granada aquifer (southern Spain). This analysis has detected several long-term cycles in the piezometric time series for the period 1966–2001. A complete analysis, using different spectral methods of those time series, has detected four different cycles: a decadal period (11 years), a 3.2-year period, an annual period and a semi-annual period. The annual cycle is ubiquitous, as a reflection of the annual hydrological cycle. The decadal cycle can be traced to hydrological variations involving a climatological regime induced by sunspot activity (11-year cycle). The piezometers located close to streamflows reaching the aquifer or along the main river, in areas where the main river recharges the aquifer, present a statistically significant decadal cycle. The explanation for the presence and spatial distribution of the 3.2-year cycle is similar to the decadal cycle, though it is identified in fewer piezometers.

In a general sense, spatial and temporal changes in groundwater levels in an aquifer are conditioned by both anthropogenic and natural factors. The former are related, mainly, to the effect of pumping, artificial recharge, etc. The latter have a number of causes, usually involving recharge by rain infiltration or superficial runoff, although other factors may be considered, such as changes due to atmospheric pressure or, in the case of coastal aquifers, the evolution of tides (Tison 1956; Nilsson 1968; Carr 1971; Fetter 1994). Frequently, series of the waterhead observed in control points over a long time show cyclic behaviour ranging from short-term (i.e. hours, days) to long-term variations (i.e. semester, annual or decadal cycles). However, our ability to study such cycles is affected, firstly, by the real presence of such cycles in the hydrodynamic behaviour of the aquifer system, and secondly, if present, with the availability of a record of measurements long enough to allow the cycles to be detected statistically. This aspect concerning experimental data series is important because long-term cycles in short records may have the appearance of a trend, which is obviously spurious.

Spectral analysis is a powerful statistical technique widely used for analysing the presence and statistical significance of cycles in time series. Given *N* experimental data of waterhead recorded in a piezometer at an equally spaced time sampling interval Δ, the range of frequencies that can be investigated goes from the Nyquist frequency 1/2Δ to the Rayleigh frequency 1/*N*Δ. Additionally, frequencies close to the Nyquist one will be affected by aliasing if the sampling interval Δ is not short enough to represent the cycle with the highest frequency affecting the waterhead evolution. Also, frequencies close to the Rayleigh one are likely to be affected by red noise (i.e. noise in the low frequencies). The origin of red noise may be the presence of a trend, long-period cycles or random components. Although spectral techniques can be used to ameliorate the previous inconveniences, the effective range of frequencies that can be investigated with confidence will be an interval shorter than the one defined by the Nyquist and Rayleigh frequencies. In our case study, for example, with the sampling interval Δ being a month and *N* ranging from 200 to 350, so that the length of the time series is from 20 to 30 years, the identifiable cycles will range from 6 months to 120 months (10 years), i.e. the cycles occur at least twice in the time series. Such are the limitations imposed by the sampling interval and the length of the records.

A careful spectral analysis determined the statistical significance of spectral peaks in the power spectrum of the waterhead time series for 53 piezometers irregularly distributed throughout the Vega de Granada alluvial aquifer (SE Spain).

## Hydrogeology of the Vega de Granada aquifer

The Vega de Granada is an important Mediterranean aquifer located in an alluvial plain surrounded by mountains in southern Spain, within the province of Granada (Fig. 1). The aquifer presents a multilayer structure, with the superposition of sedimentary materials showing a varied range of permeabilities. The surface of the aquifer has a smooth topography between the altitudes of 530 and 760 m above sea level and it is the receptor of a drainage basin of 2900 km^{2}. The annual mean rainfall over the aquifer is 450 mm, although it can reach 1000 mm above some points of the drainage basin such as the Sierra Nevada range. The aquifer has a surface of around 200 km^{2}. Its Quaternary alluvial sediments reach a thickness of 250 m in the middle and diminish towards the northern and southern borders to 50 m. The sediments are mainly gravel, sands, silts and clay, with frequent spatial changes of lithofacies, resulting in a multilayer unconfined aquifer. There is an increase of the silt and clay fraction from the central axis of the aquifer (defined by the Genil river) towards the northern and southern ends (PNUD/FAO 1972). Figure 2 shows a cross-section along the main river, where lithofacies and permeabilities were interpreted from drilling data. As seen, the changes in permeability in the aquifer are considerable, which may explain the spatial variations in the spectra calculated; however, the role of the distance to sources and to the aquifer borders would appear to be more important (Luque-Espinar 2001).

The mean transmissivity of the aquifer is around 4000 m^{2}/day with a wide range of variation from 40 000 m^{2}/day in some central sectors, to 100 m^{2}/day at the border when clay materials are more frequent. The mean effective porosity is estimated at 6%, with most values ranging from 1 to 10% (PNUD/FAO 1972). The mean groundwater flow direction is from east to west, with the steepest gradients in the NE and eastern sectors (Fig. 1).

The main inputs into the aquifer come from the infiltration of superficial runoff and return from irrigation water, plus infiltration of rainfall water. An important part of the superficial runoff comes from snow melting in the Sierra Nevada range (with altitudes up to 3000 m). Total runoff water is estimated at 400 hm^{3}/year, and the renewable aquifer resources are estimated at 180 to 230 hm^{3}/year (PNUD/FAO 1972; ITGE 1989). The output from the aquifer is mainly natural drainage towards rivers, while well pumping represents around 30 hm^{3}/year (ITGE 1989).

## Waterhead data

The original data used in this study come from the control network of the aquifer, carried out by the Instituto Geológico y Minero de España (IGME) and Confederación Hidrográfica del Guadalquivir (Ministry of Environment). The piezometric series expand up to 35 years with the oldest data from 1968. The number of piezometric series used in this work was 53, scattered all over the aquifer. Observations were made on a monthly basis for most parts, and sometimes every two months. Thus, the series average nearly 150 measurements, which is sufficient for studies of this nature. The data are measured almost monthly, but not on the same day of the month. The sampling with a constant step of one month does not introduce any considerable distortion of the spectral content. The database was checked for errors and, afterwards, the time series were sampled with a constant step of one month using linear interpolation. Figure 3 shows the time series at four piezometers, located in different sectors of the aquifer. These series can be considered representative of the water level variations (Luque-Espinar 2001). Piezometer A (Fig. 3) is situated at the western end, in a sector of the aquifer where fine-grained sediments predominate. The piezometric time series presents a cyclical aspect of long periods with minor piezometric variations; no monthly fluctuations are seen. Piezometer B (Fig. 3) is near the eastern border, in the main recharge area. In this area, gravel-size sediments are abundant, and permeability values are very high. Evolution over time shows long cycles with monthly fluctuations. Piezometer C (Fig. 3) is located in the central part of the aquifer, near the main river, which contains water all year round. The high permeability values are due to the abundance of gravel and sand. The series is characterized by distinctive annual cycles within the long-range variations. Finally, piezometer D (Fig. 3) is located between A and C, in the aquifer's discharge area. The sand-sized detritic fraction predominates, with intercalating gravel and lime parcels. Monthly oscillations are very pronounced, while long-term ones are scarcely observed.

## Spectral analysis

The statistical technique used to find cyclic components in a time series is known as spectral analysis (Jenkins & Watts 1968; Yevjevich 1972; Bras & Rodríguez-Iturbe 1985). The signal component represents the structured part of the time series, made up of a small number of embedded periodicities or cycles repeated over a long time. The noise is a random component; it may be white noise, but more often will be red noise. A time series can be represented by a finite number of measurements. In the present case, a piezometric time series is represented by a succession of waterhead data measured at regular or irregular time intervals. When a cycling component is added to another cyclic component of a longer period than the length of the time series, it will give an apparent trend. This, together with possible real trends and other factors, gives rise to noise in the low frequencies, known as red noise.

Harmonic analysis is another name used to denote the estimation of cyclic components in the time series. The time series is supposed to be a linear combination of sinusoidal functions of known periods but of unknown amplitude and phase. The modulus of the amplitude is related to the variance of the time series, explained by the oscillation at each frequency. The representation of the square of the modulus versus frequency is known as the power spectrum. There are a number of methods that can be used to infer the power spectrum: the periodogram (Papoulis 1984), the Blackman–Tukey approach (Blackman & Tukey 1958), maximum entropy (Burg 1972), and Thomson multitaper (Thomson 1982), among others. Each methodology has advantages and disadvantages, for which reason a good strategy is to use various methods and compare the results. This was done with the time series of the piezometric head; nevertheless we can affirm that the indirect method of Blackman–Tukey is a robust approach which gives acceptable results with our data sets. As an example of this comparative analysis, Figure 4 shows the results obtained with the methods of (a) maximum entropy; (b) Thomson multitaper; and (c) Blackman–Tukey. The results obtained were similar, yet we consider the Blackman–Tukey approach to be more robust and therefore more adequate for the analysis of time series.

The power spectrum is calculated from the covariance function by:
where *Ŝ*(ω) is estimated power spectrum for frequency ω, *Ĉ*(*k*) is estimated covariance function for the *k*th lag and λ(*k*) is weighting function, known as a lag-window, which is used to give less weight to the covariance estimates as the lag increases. For large lags, the estimated covariance function is less reliable. The lag-window used was the Tukey window:
where *M* is maximum number of lags for the covariance function used in the spectral estimation. The maximum number of lags is *N*−1, with *N* being the number of experimental data; however, with large values of *M* a great number of peaks will be seen in the estimated power spectrum, most representing spurious cycles. On the other hand, if *M* is very small, significant cycles would not be seen in the estimated power spectrum. For this reason we used a number of *M*=*N*/2 in order to resolve peaks, and a value of *M*=*N*/4 in order to find out which are the most significant peaks.

In addition to using a small value for *N*, confidence levels were estimated for the inferred power spectrum. Our approach consists of fitting a background power spectrum with no cyclic component, but rather a smooth continuous spectrum, which is done by fitting the spectrum of an autoregressive process of order one, i.e. AR(1). The parameter of this process is estimated from the experimental data. Then, we take into account the known result for the one-side confidence band of power spectrum estimator:
where P is probability operator, *Ŝ*(ω) is power spectrum estimate for frequency ω, *S*(ω) is underlying power spectrum for frequency ω, υ is number of degrees of freedom – for the Blackman–Tukey estimate with a Tukey lag-window, the number of degrees of freedom is 2.67*N*/*M*, is the α quantile of a chi-squared distribution with υ degrees of freedom and α is significance level. For this study, we established confidence levels of 90%, 95% and 99%.

## Results

The harmonic analysis of the time series of waterhead variations on the network of piezometers detects the presence of four distinctive periodicities: a decadal cycle, a cycle of 3.2 years, an annual cycle and a semi-annual cycle. The decadal cycle is related to the climatic 11-year cycle, in turn related to sunspot activity in the context of the North Atlantic Oscillation (NAO). The 3.2-year cycles could be related to the climatic cycle of the NAO (it has been recognized that the NAO has an influence on climate on a cyclic basis known as the quasibiennial oscillation, which ranges between two and four years). The annual cycle is related to the hydrological annual cycle, and similarly the half-year cycle is related to two precipitation seasons in a single year. Table 1 gives the percentages presented by the cycles according to the established confidence intervals.

The 11-year cycle is apparent in most of the piezometric series (85%), marking one of the main features of the temporal behaviour of the aquifer piezometric level, together with the annual cycle. The 3.2-year cycle, present in one-third of the series, also has a sparse distribution throughout the aquifer. The annual cycle appears in all the piezometric series, as expected. Finally, the semi-annual cycle is weakly represented, in just 21% of the series. This might be related to a precipitation cycle (in this area there are two precipitation seasons in a single year), which would have a minor influence on the recharging of the aquifer.

Figure 5 shows the typical power spectra found in the harmonic analysis. Not all cycles are equal in their presence and intensity in the power spectra of the different piezometers. As shown in Table 1, while decadal and annual cycles are detected in most of the piezometers, the 3.2-year and half-year cycles are detected only in a set of them. According to our probabilistic estimation, there is also a difference in the statistical significance of every cycle from one piezometer to the next and, consequently, there is a spatial variability in their importance in different areas of the aquifer.

The procedure used to assess this spatial variability was to give a code to every cycle at each piezometer according to its statistical significance in the estimated power spectrum of each piezometer. The categorical code used was:

the cycle is not distinguishable in the power spectrum;

the cycle is distinguishable but is not statistically significant at the 90% level;

the cycle is statistically significant at the 90% level but not at the 95% level;

the cycle is statistically significant at the 95% level but not at the 99% level;

the cycle is statistically significant at the 99% level.

The previous categorization, although arbitrary, allows us to highlight differences in the behaviour of cycles when the categories are represented on maps. Figure 6a–d present the spatial distribution of the statistical significance of peaks for the decadal, 3.2-year, annual and semi-annual cycles, respectively.

Figure 6a presents the spatial distribution of the statistical significance of the decadal cycle. The most significant presence of this cycle (category 4) can be seen in the borders of the aquifer and in relation to rivers and other watercourses that represent input of water from an important drainage network. We believe this is because climatic variations of around 11 years, related to the sunspot cycle, may only be seen in the variability of rainfall when integrated in an area such as the drainage basin. This implies a variability in the amount of water that enters into the aquifer from surface drainage networks; it would be detectable near the mouth of those tributary channels, but its effect is damped in the interior of the aquifer or farther away from the surface watercourse. The lesser significance of the 11-year cycle in the central area can be attributed to the small variation of the water level in this zone, where the signs of this cycle would be reduced. In the eastern sector, in contrast, the fluctuations in the piezometric level are pronounced, and decreases or increases accumulate, forming waves over an approximately 11-year period, well reflected in the spectral analysis. Moreover, for this cycle we found that the piezometric series were highly parallel to graphs of the NAO index (Hurrell 1995) in that the more positive the NAO index, the lower the piezometric levels. In this sense, it is generally accepted that the NAO marks climatic behaviour at these latitudes (Hurrel 1995; Qian *et al.* 2000; Rodrigo *et al.* 2000), though some authors relate this cycle to sunspot activity (Eddy 1976; Reid 1993).

Most researchers point to the NAO as the cause of the three-year cycle (Pozo-Vázquez *et al.* 2000), although some studies conclude that the behaviour of the rainfall in this region could be influenced by the El Niño South Oscillation (ENSO) (Rodo *et al.* 1997). Thus, the statistically significant presence of this cycle would be produced in a set of piezometers located in sectors where the vertical permeability in the proximities of surface currents is high. There is also a preferential location at points near the aquifer borders, or by torrential watersheds and more permeable detritic stretches. The NAO cycle is a subtle one that may be difficult to detect because its influence may be small in comparison with the decadal and annual cycles. We believe this cycle is exhibited in the form of minor pulses on the piezometric levels that will essentially coincide with stormy episodes. The lack of statistical significance or the absence of this cycle from other piezometers of the aquifer could be attributed mainly to a lower permeability of the sector where the piezometer is located, and to the recharge being more localized.

Figure 6c reflects how the annual cycle is highly significant (category 4) in most of the piezometers. This ubiquitous presence is logical because it is related to the yearly hydrological cycle, which affects the whole aquifer. From Figure 6d it may be seen how the semi-annual cycle is detected only at some piezometers. While there is no simple explanation for its spatial distribution, it might be related to some other climatic cycle (in this area there are two precipitation seasons in a single year).

## Conclusions

Most studies of the temporal variability of piezometric levels that appear in the literature refer to short-term cyclic variations, with periods ranging from hours to days or weeks, such as the variations deriving from tidal waves that may be observed in coastal aquifers. Some series presenting long-term fluctuations have been studied qualitatively, and then compared with rainfall data, with which there is not necessarily a parallel (Hanson & Dettinger 2005).

The most relevant aspects of the present study are the long piezometric series involved and the spectral methodological focus, as there are few bibliographic precedents. Detecting longer cycles requires a long series of records, preferably on a network of monitoring points: there may be a spatial variability in the hydrodynamic parameters of the aquifer affecting the preservation of cycles in some places, while the cycles are unobservable at other places. There were four different cycles detected in this study (not all being detected at all piezometers): decadal, 3.2-year, annual and semi-annual. The decadal cycle is perhaps the most interesting one. It is the cycle with the longest range that is related to sunspot activity, yet it is not expected in waterhead evolution. It is seen more clearly at the piezometers in close proximity to drainage channels at the border of the aquifer or along the main river entering the aquifer. This is the case because climatic variations of rainfall in relation to sunspot activity are amplified by the effect of the drainage basin, which will produce a runoff series well correlated with the climatology. Thus, it has the same effect on piezometers close to where those drainage channels enter the aquifer. The presence of the three-year cycle would have a similar general explanation, yet because the rainfall regime and its spatial distribution are different, variations in the piezometric level would be seen only in the more permeable sectors of the aquifer borders. In this way, the relationship between the surface drainage network and the aquifer recharge would be established as a reflection of the climate. The existence of cycles in the series can therefore be traced to climatic variations and the hydrodynamic characteristics of the aquifer itself, as the signal received by the piezometric levels is modulated by the functioning of the aquifer. Local or distorting phenomena that may be caused by stormy episodes are not manifest in the piezometric levels, which would clearly reflect only the predominant cycles of the regional climate.

In conclusion, we underline the importance of long-term monitoring of piezometric levels. Processing this information with the methodology put forth here constitutes, in our opinion, a sound reference for the study of climatic changes.

## Acknowledgments

We are grateful for the financial support given by the Spanish MEC (Project BTE2002-00152) and Junta de Andalucía (Group RNM122). The third author is a Ramón and Cajal Grant Holder from the Ministry of Education and Science (MEC).

- © The Geological Society of London 2008