Download Frequency Content and Mechanisms in Hydrologic Time Series: Spectral Analysis and more Lecture notes Statistics in PDF only on Docsity!
HYDROLOGICAL PROCESSES
SCIENTIFIC BRIEFING
Hydrol. Process. 16 , 565–574 (2002) DOI: 10.1002/hyp.
Practical applications of spectral analysis to hydrologic
time series
Sean W. Fleming, 1 *
A. Marsh Lavenue, 1†
Alaa H. Aly 1‡^ and
Alison Adams 2
(^1) Waterstone Environmental
Hydrology and Engineering, Inc., 1650 38th St., Suite 201E, Boulder, CO 80301, USA (^2) Tampa Bay Water, 2535 Landmark Drive, Suite 211, Clearwater, FL 33761-3930, USA
*Correspondence to: S. W. Fleming, Department of Earth and Ocean Sciences, University of British Columbia, Geophysics Building, 129- Main Mall, Vancouver, BC V6T 1Z4, Canada. E-mail: fleming [email protected] † (^) Present address: Accentera Inc., 936 Poplar Place, Boulder, CO 80304 USA. ‡ (^) Present address: 1825 Spencer St.,
Longmont, CO 80501, USA.
Abstract Fourier-transform-based spectral analysis and filtering techniques, although potentially very useful, have seen little practical application in hydrology. We provide an overview of the Fourier transform and spectral analysis and present examples of how these methods may be applied to practical hydrologic problems: determination of the frequency content of a time series; inference of the physical mechanisms responsible for this frequency content; and evaluation of the performance of a process-based simulation model used for water resource management. In all cases, the methods performed well and were reasonably straightforward to implement, highlighting their general utility. Copyright 2002 John Wiley & Sons, Ltd.
Introduction
The Fourier transform and spectral analysis have seen sophisticated use in surface and subsurface hydrology. This includes evaluation of highly nonlinear (i.e. chaotic and/or fractal) behaviour (e.g. Pasternack, 1999; Kirchner et al ., 2000), geostatistical applications (e.g. Robin et al ., 1993; Gelhar, 1993), climate change investigations (e.g. Kite, 1989; Lall and Mann, 1995), and reformulation of various hydrological inverse problems in the frequency domain (e.g. Duffy and Gelhar, 1985; Duffy and Al-Hassan, 1988). However, there is little discussion in the literature of more pragmatic hydrologic applications of spectral analysis, and this is not yet a standard method in the skill set of the typical practicing water resource scientist or engineer. The overall purpose of our paper is to bridge this gap between theory and practice. We first provide a synthesis of the more salient practical aspects of Fourier-transform-based spectral analysis, including its implementation in the popular Microsoft Excel spreadsheet. We then present illustrative examples of the application of such methods to practical hydrological problems. Taken together, the background and examples should provide a useful template or starting point for employing the techniques in other applied problems. The work presented here arose from evaluation of the performance of ISGW, a regional-scale integrated surface water/groundwater simulation model developed and used as a water resource management tool in the Tampa Bay, Florida, area (SDI, 1997). Three issues were: (1) what is the frequency content of observed discharge time series? (2) why are these cycles present in stream flow? (3) does the ISGW model honour the observed frequency content? Spectral analysis and filtering were applied to observed and ISGW-simulated hydrologic time series to address these
Received 8 May 2001 Copyright 2002 John Wiley & Sons, Ltd. (^) 565 Accepted 19 September 2001
S. W. FLEMING ET AL.
questions. Note that we present neither an overall evaluation of the ISGW model nor results of all spec- tral analyses performed. Rather, we draw from this experience examples to illustrate some of the poten- tial practical hydrological uses of spectral analysis.
Method
A number of tools are available to infer structure in time series. We used the Fourier transform because: (1) it is a very well-established method extensively used, for example, in electrical and biomedical engi- neering, geophysics, astronomy, atmospheric science, and oceanography; (2) unlike some forecast-oriented statistical methods, it does not require educated gues- ses as to the frequency content of the signal, thus per- mitting effective empirical determination of the full spectrum; (3) the results are physically intuitive; (4) it can be directly used as the basis for other analysis and processing techniques, such as frequency-domain fil- tering; and (5) a wide variety of software is available to implement the method. A large number of texts thoroughly address Fourier-transform-based spectral analysis and filtering from various perspectives (e.g. Dobrin, 1976; Press et al ., 1992; Emery and Thom- son, 1998; Smith, 1998; Bracewell, 2000); Peters (1998), Bates (1998), and Peters and Bates (1998), chapters in a single book, provide a particularly good initial introduction to the subject. A summary of the more salient aspects for application to practical hydro- logic problems is provided here. Although it is hoped that this brief and accessible tutorial is sufficiently comprehensive and rigorous for many purposes, we strongly encourage the reader to consult the above texts for a more solid grounding in the subject. The forward Fourier transform takes a time-domain signal, gt (i.e. the hydrologic time series under consideration, such as stream discharge measurements at a particular location over time) and transforms it into a frequency-domain signal, Gf, referred to as the Fourier transform of gt; f is frequency and t is time (specifically, elapsed time from the start of the data record analysed, which is taken to be t D 0). The forward Fourier transform is given by:
Gf D
gt e^2 ift^ dt 1
where i D
p 1 and all other quantities are as descri- bed previously. It is helpful to think of gt and
Gf as simply two different ways of representing the same function (e.g. Press et al ., 1992). Alterna- tively, the Fourier transform may be visualized as the decomposition of a time series into sine waves of varying amplitude, phase, and period; summing together sine waves with the characteristics described by the Fourier analysis would result in the origi- nal time series. This physical interpretation is, of course, complicated by the fact that the Fourier trans- form is complex-valued. Also note that, unlike the simpler Fourier series, the Fourier transform can be effectively applied to aperiodic functions. The inverse Fourier transform, which transforms the frequency- domain signal back to the time domain, is:
gt D
Gf e 2 ift^ df 2
There are a variety of equivalent expressions for the forward-inverse Fourier transform pair. For example, the form presented in Equations (1) and (2) corre- sponds to frequency units equal to the inverse of the time unit of the input signal, gt; if t is in months (e.g. a monthly stream discharge record), then f is in month^1. Slightly different expressions are appro- priate if alternate frequency units are used (e.g. Press et al ., 1992):
G ω D
gt eiωt^ dt and
gt D
G ω e iωt^ dt 3
where ω is angular frequency (i.e. radians/unit time); ω D 2 f. Although perhaps more prevalent in text- books, Equation (3) is generally more cumbersome than Equations (1) and (2) for application to practi- cal hydrologic problems. Also note that all the above transformations are one-dimensional insofar as the input and output signals are each a function of only one independent variable (t and f respectively). This is appropriate to hydrologic time series analysis; two- and three-dimensional forms exist for other applica- tions (e.g. Telford, 1990). For a real-valued input signal gt, the Fourier transform Gf is complex and may be presented either in rectangular notation:
Gf D Rf C iIf 4
S. W. FLEMING ET AL.
approach in terms of software availability. Potential limitations of the FFT include requirements that t is constant and that N D 2 M, where M is any inte- ger. The first requirement is made by most PSD algorithms; an exception is the Lomb–Scargle peri- odogram [Lomb, 1976; Scargle, 1982; see Press et al. (1992) for a summary]. If the input time series is not of length 2M, the second requirement may be accom- modated by truncating the data set; an alternative is to pad it with zero values to length 2M, provided that the time series is statistically stationary and has a zero mean. In a manner similar to removing the mean, subtract- ing a linear trend from the data effectively removes extremely low frequency components. Although not strictly necessary, without this the transformed sig- nal may be dominated by these typically physically uninteresting components. Power leakage across fre- quencies in the discrete transformed signal may be addressed by padding the time series with zero val- ues to give finer frequency resolution, tapering of the time series, and application of moving averages to the signal (‘data windowing’; e.g. Press et al ., 1992; Bates, 1998). Note that windowing involves a trade- off between leakage reduction and decreased resolu- tion (e.g. Bates, 1998; Smith, 1998), and can also involve a loss of information that may be significant for short data sets (e.g. Burroughs, 1992). If window- ing is applied, it is wise also to determine the PSD without this procedure to evaluate its effect upon a particular data set. We used the Microsoft Excel^ data analysis pack- age for this study, which implements Equations (1) and (2) using the FFT. Given the widespread use and availability of this software, we briefly describe its FFT capability here. This function is accessed under the Tools/Data Analysis/Fourier Analysis menus; if the Data Analysis menu is not available, the Analy- sis ToolPak must be selected under Tools/Add-Ins. If not present there, it must be loaded from the orig- inal Microsoft Office^ CD. t and gt should be present in the spreadsheet in adjacent columns. When prompted for the input field, provide only the col- umn containing gt. The output consists of a single column containing only Gf, in rectangular nota- tion [see Equation (4)]. It is in typical wrap-around format; the column of corresponding f values may easily be manually calculated in the spreadsheet [refer to Figure 12Ð 2 Ð 2 Ðb of Press et al. (1992)]. The output
must also be parsed into separate columns of Rf and If to perform subsequent manipulations, such as the PSD calculation. Regardless of the software used, a good initial check of the procedure is to ver- ify that the amplitude spectrum is symmetrical about f D 0 over the interval fc to fc. Illustrative examples of results from one discharge and one precipitation station in the Tampa Bay area are provided in the following section. Trends were determined using linear regression analysis and removed from the data prior to transformation. Time series were truncated to length 2M, with removal occurring at the start of the available data. Data win- dowing was not found to be necessary. PSD calcu- lations and filtering procedures were manually per- formed in the spreadsheet.
Results
Spectral analysis of observed discharge and precipitation Discharge and precipitation data from the Blackwater Creek and CNR-T3 gauges respectively are given in Figure 1 (refer to caption for site descriptions). The power spectrum of observed stream discharge is pro- vided in Figure 2. As N D 256 and t D 1 month, fc is 0Ð5 month^1 and the PSD frequency inter- val is 0Ð003 906 month^1. The PSD shows peaks at frequencies of 0Ð1680 month^1 , 0Ð0820 month^1 , and 0 Ð0273 month^1 or periods of approximately 6 months, 12 months, and 36 months respectively. Smaller peaks at f D 0 Ð3320, 0Ð1016, and 0Ð 0547 month^1 (periods of about 3, 9, and 18 months) are difficult to discern from ambient noise. Assuming that there is relatively little interference with discharge through engineered structures such as large-scale reservoirs and wellfields, it is reasonable to suggest that, in Florida, variation in precipitation dominates the discharge time series. If so, the fre- quency contents of the two signals should be sim- ilar. This hypothesis was tested by determining the PSD of the precipitation record at CNR-T3 (Figure 2). The very close correspondence between the frequency peaks (particularly at 3, 6, and 12 months) for two such physically distinct hydrologic variables strongly suggests that the periodicities identified in both are genuine, and confirms that precipitation is the dom- inant physical control over periodicities in discharge within the study area.
SCIENTIFIC BRIEFING
0
5
10
15
20
25
discharge [m
3 /s]
Jan-71 Jan-76 Jan-81 Jan-86 Jan-91 Jan-
0
100
200
300
400 precipitation [mm/month]
ISGW-simulated discharge observed discharge precipitation
Figure 1. Observed and ISGW model-simulated discharge at Blackwater Creek near Knights, a stream gauge located about 35 km ENE of Tampa, FL, and observed precipitation at rain gauge CNR-T3, located about 2Ð5 km ENE of the stream gauge. This stream gauge lies at the bottom of a well-defined sub-basin, rain gauges (including CNR-T3) are present within the watershed, and there are no large unnatural hydrologic disturbances within the watershed, such as major dams and reservoirs (J. Hughes, personal communication, 2000). The watershed area to the discharge gauge is approximately 290 km^2 (SDI, 1997). Precipitation data were available at daily intervals and summed to provide a monthly time series to facilitate comparison with the monthly discharge data. Note that, in this figure, the discharge data are monthly means, whereas the rainfall data are monthly totals. Discharge and precipitation data were continuous and available from 1971 through 1999
By way of method validation and further physical interpretation, plots of normalized average monthly discharge and precipitation at these gauges are provided in Figure 3. Annual peaks in discharge and precipitation occur in late summer, with secondary peaks in early spring. The primary peak is significantly greater than the secondary one, and these maxima are about 6 months apart. This accounts for the 6 and 12 month periodicities in the power spectra. Also note that though both the discharge and precipitation values show a very
strong annual maximum, the approximately 6 month seasonal variation is much more clearly defined in discharge, which shows strong, narrow peaks in March and September, than in the rainfall record, which exhibits a broader temporal distribution of rainfall over the year. This observation is reflected by the power spectra: in the discharge PSD the 6 and 12 month bands contain nearly equal power, whereas in the precipitation PSD the 6 month band contains only about half the power of the 12 month band.
SCIENTIFIC BRIEFING
0
1
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
normalized average monthly value [-]
discharge precipitation
Figure 3. Independently normalized average monthly rainfall and discharge at CNR-T3 and Blackwater Creek gauges over the period 1971 – 99, illustrating the physical basis for the presence and characteristics of the 6 and 12 month periodicities evident in the precipitation and stream flow PSDs
upon river flow. Another possibility is that atypically high precipitation may cause operators of small-scale, private reservoirs within this rural watershed to open the floodgates, suddenly increasing discharge by an inordinate amount. Signal power in the 3 month band is low, and it is difficult to correlate this frequency component unam- biguously to the monthly average rainfall and river flow data. Nonetheless, the distinct presence of this periodicity in both the discharge and rainfall PSDs leads us tentatively to regard it as genuine (perhaps reflecting subtle seasonal variations not clearly appar- ent from monthly mean values obtained by averaging over decades), with the caveat that it could simply represent noise. In either case, the low signal power indicates that the 3 month periodicity makes a small contribution to the overall precipitation and discharge
signals. Finally, the low signal power exhibited by the 9 and 18 month periodicities, the lack of a cor- relation of these frequencies to the raw time series or average monthly values, the absence of a physi- cal mechanism clearly responsible for the frequencies, and the fact that the periodicities do not appear with equal strength in the discharge and rainfall PSDs (the 18 month periodicity is essentially nonexistent in the precipitation power spectrum) combine to sug- gest that these frequencies cannot, on the basis of this analysis, be taken as genuinely present in the data. In this simple example, the main practical benefit of the spectral analysis is to reveal structure directly in the time series that could have been easily overlooked otherwise, or (equivalently) to provide independent, semi-quantitative confirmation of certain properties
S. W. FLEMING ET AL.
that might have been inferred using other approaches. These include both annual and semi-annual hydrolog- ical cycles, the significance of major tropical storms, and confirmation and better description of the rela- tionship between discharge and precipitation. An important addendum is to recall that, clearly, the data set duration and temporal discretization cho- sen must correspond to the goals of the analysis. For detailed examination of rainfall-runoff relation- ships, such as the basin response to storm events, a shorter duration of high-density (daily or, ideally, hourly) data would be appropriate, whereas exami- nation of climatic cycles would require decades of coarser (monthly or yearly) data. The data duration and density used in the medium time-scale examples provided here are adequate for their exploratory pur- poses.
Simulation model validation
A potentially highly valuable approach to assess- ing the reliability of a calibrated simulation model involves determination and comparison of the period- icities present in observed and model-simulated time series. Standard time-domain measures of model per- formance, such as root mean square error (RMSE) between measured and predicted discharge, stage, or head, are excellent overall metrics of calibration and prediction quality; moreover, a perfect time-domain match implies a perfect frequency-domain match. In practice, however, matches are imperfect and some maximum time-domain RMSE is defined, sometimes rather arbitrarily, as acceptable for a given purpose. This issue may be particularly significant when, as in the case of the ISGW model considered here, a regional-scale coupled groundwater–surface water study is undertaken and acceptable time-domain cal- ibration quality is based upon the ability to match several state variables over decades, and over a large geographic area. From a practical perspective, then, such standard performance measures typically do not speak directly to the ability of a model to replicate or forecast the structure of an observed hydrologic time series. The normalized power spectrum of simulated dis- charge was determined and is graphically compared with that of observed discharge in Figure 4. The PSDs are similar, confirming that the ISGW model ade- quately captures the periodic structure of the observed
time series at a temporal resolution of 1 month. Such a comparison essentially represents one iteration in the frequency-domain inversion procedure described, for example, by Duffy and Gelhar (1985) and Duffy and Al-Hassan (1988). Here, we perform the com- parison once as a means of evaluating a ‘completed’ model, since such replication of the overall observed frequency content may be viewed as a minimum requirement for acceptance of a calibrated model. This aspect of the application is novel and may be of particular use, for example, in regulatory or litigation support contexts.
Spectral filtering A number of advanced hydrological applications of the Fourier transform and spectral analysis were listed in the introductory section. A relatively straightfor- ward but novel and potentially powerful extension, which may be simply implemented in a spreadsheet program, is frequency-domain digital signal process- ing. For example, by applying a band-pass filter to the frequency-domain data prior to the inverse transform, which eliminates all identified periodic- ities but one, the contribution made by that fre- quency to the time-domain signal may be assessed in a highly intuitive manner. Similarly, all the peri- odicities that were clearly identified using the PSD may be passed through the filter, but those fre- quency bands that lie in the region between the peaks in the power spectrum may be omitted. Pre- liminary work with the Blackwater Creek discharge record suggests that this latter procedure easily yields a measure of the overall relative contributions of periodic versus aperiodic (i.e. stochastic, or chaotic or other aperiodic causal) components to the time series. Note that manual reconstruction of constituent sinusoidal components on the basis of FFT results (e.g. Peters, 1998) might be used to perform similar functions. Geophysical signal processing and elec- trical engineering texts, such as Dobrin (1976) and Smith (1998), provide good introductory explana- tions of how to extend the concepts and proce- dures outlined previously in this paper to include filtering.
Conclusions
Fourier-transform-based spectral analysis techniques have been previously used in sophisticated studies,
S. W. FLEMING ET AL.
Brillinger DR. 1993. The digital rainbow: some history and applica- tions of numerical spectrum analysis Canadian Journal of Statistics 21 (1): 1–19.
Burroughs WJ. 1992. Weather Cycles: Real or imaginary? Cam- bridge University Press: Cambridge, UK; 201 pp. Dobrin MB. 1976. Introduction to Geophysical Prospecting. Mc- Graw-Hill: New York, NY; 630 pp. Duffy CJ, Al-Hassan S. 1988. The time and frequency response of tracer experiments. Journal of Hydrology 97 : 59–73.
Duffy CJ, Gelhar LW. 1985. A frequency domain approach to water quality modeling in groundwater: theory. Water Resources Research 21 (8): 1175– 1184.
Emery WJ, Thomson RE. 1998. Data Analysis Methods in Physical Oceanography. Pergamon Press: London, UK; 650 pp. Gelhar LW. 1993. Stochastic Subsurface Hydrology. Prentice Hall: Englewood, NJ; 390 pp. Kirchner JW, Feng X, Neal C. 2000. Fractal stream chemistry and its implications for contaminant transport in catchments. Nature 403 : 524 – 526.
Kite G. 1989. Use of time series analysis to detect climate change. Journal of Hydrology 111 : 259– 279. Lall U, Mann M. 1995. The Great Salt Lake: a barometer of low- frequency climatic variability. Water Resources Research 31 (10): 2503 – 2515. Lomb NR. 1976. Least-squares frequency analysis of unequally spaced data. Astrophysics and Space Science 39 : 447–462.
Mitchell Jr JM, Stockton CW, Meko DM. 1979. Evidence of a 22- year rhythm of drought in the western United States related to the Hale solar cycle since the 17th century. In Solar–Terrestrial
Influences on Weather and Climate , McCormac BM, Seliga TA (eds). D. Reidel: Dordrecht, Holland; 346 pp. Pasternack GB. 1999. Does the river run wild? Assessing chaos in hydrological systems. Advances in Water Resources 23 : 253– 260. Pelletier JD. 1997. Fractal models in geology. PhD Dissertation. Cornell University: Ithaca, NY; 348 pp. Peters TM. 1998. Introduction to the Fourier transform. In The Fourier Transform in Biomedical Engineering , Peters TM, Willi- ams J (eds). Birkhauser: Boston, MA; 1–24. Peters TM, Bates JHT. 1998. The discrete Fourier transform and the fast Fourier transform. In The Fourier Transform in Biomedical Engineering , Peters TM, Williams J (eds). Birkhauser: Boston, MA; 175 – 194. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. 1992. Numerical Recipes in Fortran 77: the Art of Scientific Computing , 2nd edn. Cambridge University Press: Cambridge, UK; 933 pp. Robin MJL, Gutjahr AL, Sudicky EA, Wilson JL. 1993. Cross- correlated random field generation with the direct Fourier transform method. Water Resources Research 29 (7): 2385– 2397. Scargle JD. 1982. Studies in astronomical time series analysis. II. Statistical aspects of spectral analysis of unevenly spaced data. Astrophysics Journal 263 : 835– 853. SDI Environmental Services. 1997. Water resource evaluation and integrated hydrologic model of the central northern Tampa Bay region , SDI project WCF-690. SDI: Tampa Bay, FL; 175 pp. Smith SW. 1998. The Scientist and Engineer’s Guide to Digital Signal Processing. California Technical Publishing: San Diego, CA; 640 pp. Telford WM, Geldart LP, Sheriff RE. 1990. Applied Geophysics , 2nd edn. Cambridge University Press: Cambridge, UK; 770 pp.
