Dramatic Fluctuations of Devils Lake, North Dakota:

Climate Connections and Forecasts

Connely K. Baldwin and Upmanu Lall

Utah Water Research Laboratory, Utah State University, Logan, UT 84322-8200

Introduction

The recent (1992-date) record rise in the level of the Devils Lake, North Dakota, has led to a number of questions as to the nature of regional and global climate variability, and the utility of existing methods for forecasting lake levels and assessing the associated flood risk. A purpose of the work presented here was to explore the connection of the Devils Lake volumetric fluctuations to interannual and longer regional and global climate fluctuations, and to test the performance of recently proposed time series forecasting methods. Wiche and Vecchia (1995), and Osborne (1998) provide background information on the lake, and prior forecasting and climate analysis. Key trends in hydroclimatic variables in the Devils Lake region are first identified and discussed in the context of large-scale hydroclimate variations. Hypotheses as to operative climatic mechanisms that have led to the recent rise in the lake level are developed from this analysis. Two types of long-range lake forecasts are then considered. A forecast of lake levels for the near future (1-5 years or an inter-annual period) is developed for assessing the potential of continued flooding and associated needs for disaster relief. Second, since closed basin lakes typically exhibit long memory, procedures to estimate conditional probabilities of lake levels for extended horizons (e.g., over a 30-year flood control project, or inter-decadal periods) given current conditions were explored. Nonlinear time series analysis methods using the historical volumes of Devils Lake and selected climate indices as predictors were used to develop the inter-annual forecasts as conditional means of expected future volumes. A variety of time series modeling approaches were explored for the inter-decadal forecasts. Results are presented here for a linear, Bayesian autoregressive time series model that incorporates model and parameter uncertainty. We conclude that direct applications of existing time series analysis methods are not well suited for the development of long-range probabilistic forecasts of Devils Lake. The recent trends of Devils Lake are consistent with large-scale changes seen elsewhere. However, whether these changes are part of the natural long-term variability of climate or represent a changed climate due to human influence remains inconclusive. Consequently, while we are able to relate the recent Devils Lake trend to causative hydroclimatic factors, we are unable to confidently predict the long-term future levels of the lake. Only qualitative remarks are offered to characterize the uncertainty associated with using the past as a guide to the future of Devils Lake.

Trends in Local, Regional and Hemispheric Hydroclimatic Variables

Like most closed basin lakes, the Devils Lake exhibits dramatic fluctuations (Figure 1) over decadal and longer periods, that derive from climatic fluctuations. There is a limited literature (e.g., Lall and Mann, 1995, Mann et al, 1995, Pusc, 1993, Wiche et al., 1986) diagnosing the climatic causes of such fluctuations. The contributing drainage area of such lakes varies with climate state. The chain lakes above Devils Lake and other depressions store water during dry periods, but contribute runoff during protracted wet periods. This change in drainage area may be a key explanation for the dramatic changes in the Lakes volume, subsequent to moderate changes in the climate signal. Notable increases occur in 1950s, 1970, 1980s, and the late 1990s. These periods are also important at the regional streamflow stations.

It is useful to first look at the annual cycle of monthly changes in the Devils Lake volume (Figure 2), to motivate the search for trends in monthly precipitation, temperature, Sea Level Pressure (SLP), and cloudiness. The influence of ice cover in the winter months decreases the magnitude of changes observed. On average, the lake volume increases in spring (April through June), due to snowmelt runoff. Decreases occur in the summer months, with the greatest decrease usually occurring in August. However, increases or decreases can occur in any calendar month. At first glance, one would suspect that the winter/spring precipitation and summer/fall temperature (and hence atmospheric circulation patterns), are most important for diagnosing the changes in the volume of the lake. Monthly trends and the base climatology of these variables are reviewed next.

Figure 1. (a) Historical, standardized time series of the Great Salt Lake, Utah and Devils Lake, North Dakota. Note the dramatic fluctuations of both N. American closed basin lakes, that are often in phase over the period of common record, while the pattern for the amplitudes of positive or negative excursions is quite different. The “X” shows the approximate DL volume around 1840. Note that the annual DL lake level data prior to 1941 has been stochastically disaggregated to the monthly time-scale. (b) Annual change in volume of the two lakes show much weaker correspondence.

Figure 2. Climatology of monthly lake volume changes. Based on monthly lake volume data from October 1941 to May 1998.

Figure 3. Location of precipitation stations.

The climatology of monthly precipitation amounts for selected stations (Figure 3) is shown in (Figure 4). The data for these stations was obtained from ftp://ftp.ncdc.noaa.gov/pub/data/ushcn/daily/. Note that precipitation peaks in June, and that September/October precipitation can occasionally be as large as the average June precipitation. The range can be dramatic (e.g., 0.5 to 12 for June at Bottineau), even though the total average amounts reflect an arid climate. Winter precipitation is on average much smaller than the summer/fall precipitation in the region. It is consequently remarkable that the largest annual volume increase for the lake is seen in April. The subtle role of subsurface hydrologic processes and evaporation in modulating the lake inflows and hence the volume is indicated.

Monthly precipitation data were used to investigate the variability in precipitation. Data were obtained from the U.S. HCN web site (ftp://ftp.ncdc.noaa.gov/pub/data/ushcn/) and the FILNET version of the data, are used. The accompanying literature notes that the FILNET data has been adjusted for the time of observation bias, maximum and minimum temperature system (MMTS) bias, station moves and changes bias, and contains estimated values for missing/outlier data. These adjustments provide a long, serially complete data set useful. An examination of these records reveals largely consistent precipitation trends in the recent period. Month-by-month trends for the Langdon station are presented in Figure 5. Note that no major positive trends in the winter/spring (Dec-Apr) precipitation are evident, while May through August and October/November exhibit a positive trend in the recent period. Given that the summer and fall are usually periods of lake volume decrease, these precipitation trends are significant in that they may imply a situation where the lakes annual cycle of decrease may be reversed. The increased regional wetness in the summer may also lead to a higher contributing drainage area and runoff, and to a lower potential evaporation associated with increased local humidity. Recently, Karl et al. (1998) investigated the secular trends in both amount and intensity of precipitation in the U.S. with relatively long data sets (1910-1996). They found positive trends of heavy rainfall in the region covering the nearby upper Mississippi River basin in all seasons, except winter, where a negative trend was found. The largest trend reported for the region was for summer, consistent with the finding of Angel and Huff (1995). These analyses are consistent with our observations from the local station data.

An examination of the 1896 to 1996 monthly temperature record at the Langdon station (Figure 6) shows no unprecedented changes over the recent period of concern. January, April, May and July temperatures have reversed earlier warming trends during the 1990s. The period of decreasing and extremely low lake levels may be associated with lower summer precipitation and warmer summers. The cooler and wetter summer/fall regime is associated with the recent rise of the Devils Lake.

A comparison between the July and October precipitation trends and the local Sea Level Pressure (SLP) trends is instructive. From Figure 7, we see that for both months the increasing precipitation trend is accompanied by a decrease in the mean monthly SLP for the month. Similar trends in the local SLP are evident for June and August. This observation sets the stage for an examination of trends in the larger-scale precipitation, river flow and atmospheric circulation fields.

Figure 4. Climatology of precipitation at selected stations. The boxplot includes the median (white line), the 25th and 75th percentiles (the box limits), the 10th and 90th percentiles (whiskers), and the outliers (lines outside whiskers). The abbreviations correspond to stations shown in Figure 3.

Figure 5. Raw time series (light solid) and 20-year Loess smooth (heavy solid) of monthly precipitation at the Langdon Experiment Station. Note the May through August and October trends.

Figure 6. Raw time series (light solid) and 20-year Loess smooth (heavy solid) of monthly temperature at the Langdon Experiment Station.

Figure 7. Sea level pressure (SLP) anomaly and Langdon precipitation time series plots for July and October. The heavy line is a 20-year Loess smooth of the raw data (light solid). Note the recent large decrease in SLP in both months and the associated increase in precipitation. The summer SLP values in 1992-98 correspond to a 6-8 m average reduction in the height of the regional 700 mb geopotential surface.

Figure 8. Annual average flow for each water year of the period of record (1874 to 1996) of the Mississippi River at the Clinton, IA gaging station. The raw flows (solid) and two Loess smooths using a 10-year span (dotted) and 60-year span (dashed) are shown.

Baldwin and Lall (1999) note changes in the seasonal and long-term mean flow of the Upper Mississippi River (Figure 8). Specifically, they note changes in the seasonality of flow, precipitation, and atmospheric circulation in the region that translate into a secondary October peak and an enhanced summer peak for the flow at Clinton, IA in the recent period. The general trends in the flow of the Upper Mississippi River are consistent with those of the Devils Lake. Specifically, both the Devils Lake volume and the heaviest smooth of the Upper Mississippi River flows (dashed line in Figure 8) change direction around 1940. There are indications that this change-point is a regional- to continental-scale phenomenon.

An increase in mean annual cloud amounts of 27% occurred at Bismarck, North Dakota between 1930-1950 (Figure 9a). A Bayesian change-point analysis (Lee and Heghinian, 1977) indicates that the most likely year the regime changed is 1940. This coincides with the beginning of a general increase in Devils Lake volume over the next 60 years. One explanation may be that the increase in cloudiness reduces potential evaporation by decreasing incoming solar radiation. Since evaporation is the only way for water to escape a closed-basin lake such as the Devils Lake, such a change could have a dramatic impact. The climatology of cloudiness for the Devils Lake area (Figure 9b) indicates that the months of July and August have the highest proportion of clear days, months that correspond to decreases in the Devils Lake volume on average (Figure 2). Hence, increases in cloudiness in these months would likely decrease evaporative losses from the Devils Lake.

This change in cloudiness is not a local phenomenon. A similar increase occurred in many cities in the western U.S. between 1937-1942 (Baldwin et al, 1999; Lee and Heghinian, 1977). Thus, it appears that the climatic factors that have an impact on Devils Lake are at least regional and approach continental scale.

Figure 9. (a) Time series and 5-year moving average of mean annual cloudiness at Bismark, ND (Steurer and Karl, no date); (b) Climatology of cloudiness at Devils Lake, ND (Jensen, no date).

A continental perspective on the spatial pattern of the 1992-98 precipitation anomalies for July and October is presented in Figure 10, using the U.S. Climate Division data. The July rainfall trend extends over the Upper Mississippi river basin and into the Pacific Northwest, and the extreme Northeastern part of the country. By contrast the southern part of the country has tended to be drier than normal in the 1992-8 period in July. The picture for October reveals that the Devils Lake region is near the core of a significant national pattern of regional positive and negative precipitation anomalies. As a contrast, the precipitation anomaly patterns for a set of dry years (1940-45) for Devils Lake are shown in Figure 11. Note the opposite, but qualitatively similar spatial structure of the anomaly fields, especially for October. Opposite phases of the same climate mechanisms may thus be at play. Similar plots (using the same years) for the winter seasons did not reveal any large, structured precipitation anomalies across the country.

Figure 10. Composites of standardized U.S. precipitation anomalies for 1992-1998 July (a) and October (b). Standard deviations are computed with respect to the 1950-1995 mean at each climate division. The color codes are in units of standard deviation.

Figure 11. Composites of standardized U.S. precipitation anomalies for 1940-1945 July (a) and October (b). Standard deviations are computed with respect to the 1950-1995 mean at each climate division. The color codes are in units of standard deviation.

We now examine the spatial structure of N. Hemisphere summer and fall atmospheric circulation fields through composited anomalies of extratropical SLP for the same years as the precipitation composites presented in Figures 10 and 11. The average summer pressure fields (Figure 12a) are marked by climatological highs in the Pacific and Atlantic Oceans near 35N, and continental low pressure centers over Asia, Southwest U.S., and Eastern Canada. Closed basin lakes such as the Devils Lake and the Great Salt Lake occur in regions where the long-term potential evaporation rate exceeds the precipitation rate. As seen from Figures 12 and 13, such basins correspond to areas with a ridge in atmospheric pressure on average. Persistent changes in large-scale pressure patterns that weaken this ridge bring storm systems to the area. The summer circulation fields are generally weaker than the winter/fall fields (compare Figures 12a and 13a). By October, continental high pressure centers (East Asia, central N. America) become established, and low pressure centers in the northern oceans become more pronounced leading to a wave pattern as the extratropical jet stream gains strength and moves towards the equator. The ridge in the Devils Lake region is pronounced in the average picture reflecting the low probability of large-scale storms that can penetrate the region and lead to precipitation. The 1940-45 October (Figure 13b) is marked by an expanded and stronger central Pacific high pressure system, a more pronounced Aleutian low, and an enhanced continental ridge, and a weaker low East of Canada. These combine to keep storms out of the Devils Lake region. By contrast, the 1992-97 wet period (Figure 13c) is marked by a weakening of the October ridge over the Devils Lake area and an associated broad weakening of the continental and coastal high pressure systems around N. America, and of the low pressure systems East of Canada and in the Aleutians. An intensification and displacement of the East Asian high, and the North Asian low-pressure systems is also noted. These asymmetric changes of the average pressure systems in the regions that control the jet-stream dynamics may allow localized circulations such as the summertime low-level jet to remain established in October, and lead to the spatial structure of the rainfall anomalies (Figure 11) over N. America. Similar changes are evident in the summer analysis in Figure 12.

Figure 12. (a) Mean summer (June, July, August) SLP (mb) for the Northern Hemisphere, north of 20N; (b) The SLP anomaly field for the 1940-45 dry period; and (c) the anomaly field for the 1992-97 wet period. The approximate location of Devils Lake is marked by the *. The Great Salt Lake is located at 42N and 112W.

Figure 13. (a) Mean October SLP (mb) for the Northern Hemisphere, north of 20N; (b) The SLP anomaly field for the 1940-45 dry period; and (c) the anomaly field for the 1992-97 wet period. The approximate location of Devils Lake is marked by the *.

Connections to Low Frequency Climate Indices

An explanation of the recent rise of Devils Lake and a forecast of its future volumes consequently requires the characterization of operative multi-year climate regimes that may affect regional summer and fall precipitation. Three quasi-oscillatory climate regimes have recently been identified as the primary determinants of interannual to interdecadal climate variability in the Northern Hemisphere. These are the El Nino Southern Oscillation (ENSO, Battisti and Sarachick, 1995) with interannual (3-7 year) scales of characteristic variability, the Pacific Decadal Oscillation (PDO, Mantua et al., 1997) with interdecadal (16-20-year) characteristics, and the North Atlantic Oscillation (NAO, Kushnir, 1994) with decadal (8-10-years) variability. These three modes are often termed quasi-oscillatory, since they have recurrent phases with preferred time scales rather than an exactly periodic behavior like the seasonal cycle. The ENSO, PDO, and NAO are hypothesized to be the outcome of bi-directional, multi-scale interactions between slowly varying ocean conditions (temperature, height and velocities) and rapidly varying atmospheric conditions (winds and pressure systems). They reflect persistent ocean and atmospheric conditions in the tropical Pacific (ENSO), the extratropical Pacific (PDO) and the Atlantic (NAO) Oceans that can rapidly switch and lead to hemispheric impacts on atmospheric pressure and storm migration patterns. ENSO is by far the best studied and understood of these modes. Interactions between these climate modes and their teleconnections to continental climate have been noted, but are not yet well understood. These are areas of active diagnostic research using data and models. In the current context, the interaction between these modes in producing the atmospheric conditions associated with Devils Lake variations is of interest.

ENSO has been rather active lately, with 1997-8 being one of the strongest El Nino events on record, and an extended El Nino event persisted from 1991-95. Several “indices” that provide a measure of the ENSO state are available. There is also some discussion (Trenberth and Hoar, 1995, Rajagopalan et al, 1995) as to whether the frequency of ENSO events is changing and whether it is related to anthropogenic climate change. Since this debate is inconclusive, we shall strive to use the indices related to these oscillations as diagnostic tools for understanding oceanic teleconnections to Devils Lake. Here, we have used the Sea Surface Temperature (SST) in the NINO3 region (5N-5S, 150W-90W) in the Eastern tropical Pacific. A positive value of this index corresponds to El Nino conditions, and a negative value to La Nina conditions. The PDO is represented by the PDO Index develop by Nathan Mantua (Mantua et al., 1997). It corresponds to the leading spatial pattern from a Principal Component Analysis of N. Pacific (30N to 70N) SSTs. An NAO Index is also available, as the difference of normalized sea level pressures (SLP) between Lisbon, Portugal and Stykkisholmur, Iceland. The time series of the NINO3, PDO, and NAO indices have statistically significant cross-correlations at various lags, but have rather different time scales of variation based on spectral analysis.

The time-frequency variation in the three indices, and in the monthly change in the volume of Devils Lake and the Great Salt Lake is examined (Figure 14) through a wavelet analysis (Torrence and Compo, 1998) of their monthly time series. The spectra of these series are generally quite different in their overall character with different frequency bands emphasized. Note the intermittent nature of the narrow band oscillations (indicated by each global wavelet analysis and the corresponding time-frequency plot) in the series, over the period of record. The characteristic time scales of NINO3 are interannual, of the NAO, decadal, and the PDO interannual to interdecadal. The two lakes are more similar to each other than to the indicators, as would be expected. One expects the lakes to have a general long memory behavior. However, they are more organized in frequency bands than say the PDO, which despite preferred interannual and multi- decadal attributes is redder in

Figure 14. Wavelet spectra for NINO3, PDO, NAO, and monthly volume change of Devils Lake, and Great Salt Lake. The full record spectrum is shown to the right. The dashed red line in each such plot is the 95% significance level for red noise. The wavelet spectrum shows the corresponding variation in spectral power for each period as a function of time. For Devils Lake, most of the power at all frequencies is concentrated in the last decade, emphasizing the unusual nature of this period. Considerable time variation in the time-frequency structure as well as common time-frequency structures across the series are also notable.

character. The Devils Lake annual cycle is not nearly as clear as for the Great Salt Lake reflecting somewhat different basin/climate dynamics. The wavelet spectrum for the Great Salt Lake shows a much greater time-frequency commonality with the climate index spectra, than does the Devils Lake. The recent period is clearly the most anomalous for the Devils Lake record.

The role of the three low-frequency climate patterns in generating such an anomaly is explored through an analysis of the correlation between the three climate indices and summer and October continental precipitation in Figure 15. While the spatial patterns of correlation of each index with continental precipitation differ for each season, particularly for the summer, they (a) show strong spatially coherent response structures for each index, and (b) have statistically significant correlations with the general Devils Lake region that are consistent with the increased regional precipitation for these seasons shown in Figures 5 and 7. Note that NINO3 has been positive on average over this period with a protracted positive anomaly over 1991-95 and the large positive anomaly through 1997-8, and is positively correlated with the Devils Lake region precipitation for summer. Likewise PDO has been in its positive phase and is positively correlated while the NAO, which is negatively correlated with precipitation in the region, has been in its negative phase. The correlation patterns of the three indices with the summer 700 mb geopotential height are shown in Figure 16. The correlation of the Devils Lake region’s pressure with NINO3 is not significant. However, the correlations of the PDO and NAO with the regional pressure surface are consistent with the corresponding indications for precipitation. Thus it appears that the conditions in the Pacific and the Atlantic Ocean may jointly influence the atmospheric regime that leads to anomalous summer precipitation in the Devils Lake region. Indeed correlations between the indices and the atmospheric pressure over the U.S. for June through August are consistent with patterns that would steer storms to the Devils Lake area.Similar, but weaker correlations are noted for October.

The connections of the climate index teleconnections to Devils Lake volume changes were further explored through multi-taper spectral analysis. The spectral coherence (correlation at specific frequency bands) between the monthly indices and monthly changes in lake volume was estimated using the multi-taper method of spectral analysis (Thomson, 1982, Mann and Park, 1995). The analysis was performed for all three climate indices and the product of all possible pair-wise combinations (e.g., NAO*PDO) of these indices as a crude approximation to the physical interaction of the climate modes. Significance was assessed in terms of the spectra of the individual series correlated and the coherence at a given frequency. The results indicate that there are three dominant modes of coherent variability, with periods of 12-25 years, 3-5 years, and 2 years. PDO is coherent with the lake volume change with periods of 12-25 and 2 years, while NINO3, NAO, and the product of NINO3 and NAO are coherent with the lake with periods of 3-5 years. Note that the time-frequency analysis using wavelet methods shows considerable time variation in the frequency structure of the Devils Lake monthly volume changes and the climate indices, except for the recent period. Thus it is likely that these teleconnections are largely episodic. Our hypothesis is that these spectral peaks do not necessarily correspond to exactly periodic behavior and coupling at these frequencies. Rather, they represent organization in a space-time dynamical process that leads to preferred, narrow band, spectral signatures.

The significance of the interaction terms in the coherency analysis motivates the search for smaller-scale patterns of SST that may be more directly related to DL volume. This is accomplished by computing the lag cross-correlation between changes in DL volume and gridded SST values in the Pacific and Atlantic Oceans. The patterns observed in the Pacific Ocean were very similar to observed patterns already represented in other indices (PDO, ENSO). However, in the Atlantic Ocean five areas of SST (Figure 17) were identified as potential predictors. Preliminary forecasts of DL volume using all of these SST areas as predictors, as well as the indices for the large-scale climate patterns, consistently chose one climate index (PDO) and one SST area (SEC), as detailed in a later section. A lag cross-correlation analysis indicates that both ENSO and NAO are significantly correlated with SSTs in the SEC area.

Since the Devils Lake and basin hydrology act as integrators of the larger-scale atmospheric circulation and of regional precipitation, it is also useful to look at the time series of the Devils Lake volume in conjunction with running (cumulative) sums of the departures of these indices from their long-term mean. The perspective here is that the large-scale climate states suggested by these indices translate into preferred, concurrent local precipitation and temperature signals, that are modified in some way by the basin hydrology, and then added up over time as lake volumes. Thus, the lake volume can be thought of as a random walk through time where the increments of the random walk are determined in some way by the underlying climate states. Since we dont

Figure 15. Correlation of U.S. summer (June, July, August) rainfall with (a) NINO3, (b) PDO and (c) NAO, and of October precipitation with (d) NINO3, (e) PDO and (f) NAO Indices for the period 1958-1998. Correlations larger than 0.32 in absolute magnitude are significant at the 95% level. The summer correlations are stronger than those for October, with all 3 indices significantly correlated to the Devils Lake region in the summer, and only the NINO3 significant in October.

Figure 16. Areal averages of sea surface temperature in the Atlantic Ocean that are correlated with changes in Devils Lake volume. The acronyms indicate the approximate location of the area: South-East coast (SEC), Gulf of Mexico (GM), West African coast (WAC), European coast (EUR), and Labrador (LAB).

know precisely how these influences can be described, it is useful to compare the cumulative tendencies of the presumed causative factors (selected climate indices) with the resulting random walk (the Devils Lake volume). This comparison is presented below in Figure 18. Recall that the slope of a cumulative sum curve represents the local mean of the process. In this respect, the NINO3 and the SEC index have recently been at rather different levels than their long-term mean (which is 0), and the PDO has recovered from a 1942-1976 period of an anomalously negative mean. If we assume that these indices are useful predictors of concurrent changes in monthly Devils Lake volume, then these rather long-term changes in state are important for understanding the potential changes in the long-run probabilities of lake volume.

Figure 18. Cumulative trends (obtained through a running sum of anomalies taken about the long-term mean) of NINO3, PDO and SEC climate indices compared to the Devils Lake volume. Note the turning points in the climate series marked by the solid vertical lines around 1942 and 1976. The cumulative SEC and NINO3 indices are at extreme anomalous levels, like the Devils Lake volume at the end of the period.

In summary, evidence for correlation between summer precipitation and atmospheric pressure regimes in the Devils Lake region and indices of quasi-oscillatory climate modes related to the Pacific and Atlantic Oceans was found. The superposition of specific phases of these modes of natural climate variability may then be responsible for the recent, anomalous summer/fall wetness and rise of the lake. Since the indices used to diagnose these climate modes are themselves correlated, one has to be careful in making causal statements about their individual or joint role. Further, the analyses presented thus far consider concurrent variability of the climate indices and the regional hydroclimatic variable of interest. It is not clear that these 0-lag correlations are particularly useful for long lead forecasting of Devils Lake volumes. They merely establish that slowly varying sea surface temperature conditions that are known to be associated with interannual to interdecadal atmospheric circulation pattern changes are relevant for understanding the fluctuations of Devils Lake. No conclusive long run demonstrations of the ability to statistically or deterministically forecast these climate indices have so far emerged. Thus, at present, it is important to recognize that it is possible to establish such links and to think about what the low-frequency nature of the underlying climate state implies for operational and planning models. Clearly, one would need to develop paradigms for both the nature of dependence of the hydrologic variable (Devils Lake) on the climate state, and for the time variation of the underlying climate state. Formal investigations of this sort were not pursued here. The forecasting methods described in the next two sections approach the problem stated here empirically, in ways that have been useful in other related contexts.

Interannual Forecasts Using Nonlinear Time Series Analysis Methods

A nonlinear time series modeling framework that has its roots in the literature (see Diks, 1999; Kantz and Schreiber, 1997; Cutler and Kaplan, 1997; Abarbanel, 1995) on state space reconstruction of dynamical systems from time series data is considered here. The essence of this approach is that attributes of a multivariate dynamical system can be reconstructed from time series of one or more state variables through embedding, An embedding is defined by constructing a multivariate state space using appropriately time lagged copies of the observed time series. For instance if we had observations of only the Devils Lake volume, (denoted as xt, t=1n), then Takens embedding theorem indicates that the dynamics of the lake volume can be reconstructed by forming the embedding Dt defined as (xt, xt-;#61556;, xt-2;#61556;,. xt-(d-1);#61472;;#61556;), where ;#61556; is an appropriate sampling frequency, d is the total number of lags considered and d (2p+1), and p is the number of state variables that interact to produce the fluctuations of the lake volume series, xt. Forecasts of future values of the state variable are then possible if the following relationship can be identified from the available data:

(1)

where T is the lead-time for the forecast, and f(.) is a function that describes the dynamical relationship between the past states and the future state.

In practice, we have a finite amount of data, and appropriate values of ;#61556; and d, as well as the form of f(.) are not known a priori. If the function f(.) is linear in its arguments, this modeling framework is similar to traditional autoregressive time series models. Forecasts from such systems tend to dissipate to the long-term mean value of the time series from whatever condition they are in, unless special efforts are made to include oscillatory or other long-memory features. As was discussed earlier, one would be hard pressed to explain the fluctuations of Devils Lake or other similar processes under the traditional, linear, stationary system paradigm. However, if f(.) is nonlinear, there is a possibility of regime and oscillatory dynamics as well as chaotic dynamics. The system may stay in a particular regime for some time and then switch out to a different regime of behavior. Chaos or loss of predictability is often associated with the regime transitions, as exemplified in the famous work of Lorenz (1963). Oscillations that are the outcomes of positive and negative feedbacks within the system may occur within regimes or across regimes at a variety of time scales. These are the aspects of the climate system that were of interest in the preceding section.

Lall et al (1995) developed a forecasting model for the Great Salt Lake (GSL) using the above ideas, where a nonparametric, spline regression methodology was used to estimate f(.), and statistical criteria were used to choose ;#61556;, d, and the subset of lagged coordinates used in building the model. They noted that certain regime transitions (e.g., the start of the 1983 rise of the lake) of the Great Salt Lake were not predictable even a few months in advance, but in general one could expect relatively accurate forecasts 1-4 years into the future, even during the extreme rise and fall of the lake. Data on the Great Salt Lake had been reconstructed back to 1847 for the analysis. This allowed some of the extreme fluctuations in the 19th century that are similar to the 1980’s GSL fluctuations to be represented in the data set available for model building. Unfortunately, for Devils Lake, the record could only be reconstructed back to 1905, limiting the ability to reconstruct the dynamics associated with the extreme recent fluctuations. Consequently, in our work here, we have used an extension (see Moon, 1995; Ames (1998)) of the Lall et al (1995) algorithm, that allows for the reconstruction of the dynamics of a target variable (xt) using time series of selected climate indicators. The general forecasting model is represented as :

(2)

where y1, y2,..ym, refer to m potential auxiliary predictors (e.g., climate indicators), with associated sampling frequencies ;#61556;;#61489;;#61484;;#61472;;#61556;;#61490;;#61484;;#61472;;#61486;;#61486;;#61486;;#61556;m, and embedding dimensions d1, d2, dm, and et is an error process that includes components due to measurement error and due to approximation error in estimating f(.). The approximation error may result from under-specification of the true state space (useful predictors are missing), or from limitations of the numerical scheme used to fit f(.).

Forecasts using expressions like those in equations (1) or (2) can be produced using an iterated or direct approach. For the 1-step iterated approach, the next value of the time series (T=1 in the equations) is forecast. This is then used as a known value at time step t+1, and the model is re-applied to forecast the next time step (t+2). This process is repeated T times until the desired lead-time for the forecast is reached. Only the existing data is used for fitting the 1-step ahead forecasting function function f(.). The estimated values of xt+1, xt+2 etc., are used only to compute new iterates and not to re-fit the function f(.). The iterated approach is consequently similar to the traditional autoregressive modeling approach. In the direct method, separate 1-step, 2-step, T-step ahead models are fit to the data, and are directly applied to generate 1, 2, T-step ahead forecasts. The two forecasting methods can provide different results depending on the relative signal-to-noise ratio (relative magnitude of the variance of the error term et, to the variance of xt), and on local variations in predictability that depend on the nonlinearity of the underlying f(.). Both methods of forecasting were evaluated in a cross-validated testing mode with the Devils Lake data. The models were fit to selected portions of the data and tested on the remainder.

For any candidate set of predictors in a particular fitting exercise, the predictors retained in the model as well as the complexity of the model (e.g., number and placement of knots for the regression spline) is selected using traditional statistical criteria (e.g., Generalized Cross Validation or GCV, and the Schwarz Criteria, SC). The use of logarithmic and square root transforms of the Devils Lake volume were also explored in the model building process. Suitably chosen predictors with different intrinsic time scales of fluctuations (e.g., interannual to decadal for ocean temperatures and seasonal for local precipitation) can potentially be used to reconstruct the short and long run dynamics of Devils Lake. A variety of numerical algorithms (e.g., spline regression, locally weighted polynomial regression, and neural networks) were explored for estimating f(.). Multivariate, adaptive regression splines (Friedman, 1991) encoded in a Windows 95 application (Ames, 1998) that focuses on time series model building and forecasting were used in the work reported here.

Pre-screening of potential predictors:

One interesting implication of Takens theorem is that if multiple, lagged variables are used to reconstruct the state space of a dynamical system, there is a potential for significant coordinate redundancy, since it is conceptually possible to reconstruct the state space by lagging any one of the state variables. The statistical criteria we used to select predictors seek to limit such redundancies and their effects on the forecast scheme. However, as the number of potential predictors and hence choices for the statistical criteria increases, there is increasing potential for model mis-specification. Consequently, it is important to pre-screen the potential predictor set, before a model such as in (2) is fit to the data. We used the correlative analyses described in the preceding section and fit a set of candidate models with different subsets of potential predictors. Efforts were made to include potential predictors in each candidate set that span a set of intrinsic time scales of variability.

After the appropriate transformations are applied, a series of trial forecasts are performed to identify the most important predictors. The predictors considered were the PDO, NAO, and NINO3 indices, the five SST areas indicated in Figure 1, and the monthly precipitation anomalies for climate division 3 of North Dakota (PCP). The cumulative sum of the predictor was also considered as a predictor (Corradi, 1995). Several MARS models were fit using different combinations of these predictors at several different lead-times and from different starting dates in the historical record. The most important predictors identified are the PDO index and the SEC area of SST. The PCP predictor is important in certain cases. Often models that used just the time history of Devils Lake performed as well as those that used climate predictors.

Validation forecasts and outlook for future Devils Lake levels

Once these predictors were identified, iterated and direct forecasts at four different lead-times were made to explore the predictability of the DL volume series in the historical record. Comparisons for direct forecasts are shown in Figure 19 for: (a) the period of relatively steady lake volume beginning in 1981, (b) the 1987 transition to decreasing lake volume, (c) the increasing lake volume from 1997-1999, and (d) a blind outlook for future lake levels. We consider four different lead-times: 12, 18, 24, and 36 months. The potential predictors provided to the model for these forecasts were the past volumes of Devils Lake (with a square root transform), PDO, SEC, and N. Dakota Climate division (3) precipitation. Cumulative sums of anomalies from

Figure 19. Forecasts of Devils Lake volume, converted to levels, starting at various times for different lead times as shown are given in the legend. These are all direct forecasts, based on models fit using only data prior to the start date of the forecast. For a 36-month lead forecast starting in January 1981, only data up to January 1978 (i.e. 36 months prior) is used for model fitting. The subsequent months data is used to generate the forecast with the fitted model. Thus, at the end of the 36 month forecast, data from December 1980 would be used. MARS chooses different predictors from the candidate set. There was considerable variation in the predictors used for models fit at different times, and for different leads. For example, the 36 month ahead forecast model fit for data up to May 1999, uses Devils Lake (t-36), Cum. Sum SEC (t-36, t-72, t-84), Cum. Sum PDO (t-84) and N. Dakota Climate Division Precipitation (t-84) as predictors.

the long-term average were used for all predictors except Devils Lake. Considerable variation was noted in the fitted models for different lead times and for different times. One finding was that sometimes using the precipitation (PCP) as a predictor can improve the longer lead-time forecasts. The 24- and 36-month lead forecasts mentioned above used PCP as a predictor. The 24-month forecast mentioned above without PCP failed to predict the decrease after July 1997 (Figure 19c). Using PCP as a predictor also improved the 36-month lead forecast starting from January 1997. The forecast made without PCP was too extreme, over-predicting the volume for January 1999 by 207,000 acre-feet. However, when PCP was included as a predictor the accuracy improved, predicting the correct volume within 12,000 acre-feet (0.6%).

The outlook on future Devils Lake levels (up to mid 2002) indicate decreasing lake levels. This is a rather surprising forecast given the recent increases. However, the outlook agrees with long-term forecasts of the Upper Mississippi River (UMR) streamflow, which indicate that 2000 will be near-normal, but that 2001 will be below-normal (Baldwin, 1998, also available at: http://pub.uwrl.usu.edu/cbald/dlake/cmon-new2-fcast-updated.gif). We compare these forecasts with those of the UMR since they have been demonstrated to be more accurate than those of Devils Lake directly, and provide a useful validation.

Bayesian time series forecasting methodology

Wiche and Vecchia (1995) presented an effective implementation of established, stationary, hydrologic time series analysis methods for the analysis of Devils Lake volumes. Unfortunately, such methods, can have a hard time reproducing features such as the recent rise of the Devils Lake, even with parameter uncertainty considered. In addressing our second objective for forecasting that entailed the determination of long-run lake volume probabilities conditional on current state, we considered several alternatives to the Wiche and Vecchia work. These included Fractionally Integrated Autoregressive Moving Average models (ARFIMA), and a Bayesian autoregressive modeling approach (ARCOMP) that considers uncertainty in both the model parameters and coefficients. ARCOMP, is a relatively new approach due to Huerta and West (1999), that admits certain long memory and quasi-periodic sub-processes. A direct application of the ARFIMA model to the monthly Devils Lake volume (log transform) time series leads to the selection by AIC of a (AR=5, d=0.32, MA=1) model, with model coefficients (AR: 1.83, -0.97,0.08, -0.004,0.068; MA: 0.73). The forecasts from this model tended to the mean of the series, and were unsuccessful for the 1990s. The ARCOMP procedure is described below.

Consider, a univariate autoregressive process of order p, AR(p):

(3)

where the ;#61542;I are autoregressive coefficients for lag i, and ;#61541;t is a an independent, noise process.

This process can also be written as:

(4)

where B is the backshift operator; ;#61537;j are the roots of the characteristic polynomial associated with the AR(p) process, including the (R=p-2C) real roots rj, for j=2C+1,p, and the 2C complex roots, , for j=12C, corresponding to quasi-periodic processes with frequency ;#61559;j; ztj and atj are latent processes corresponding to the complex and real roots respectively; and bj, dj, and ej, are some real constants.

Huerta and West note that (a) state space models can be written as ARMA(p’,q’) models, that can in turn be approximated as high order AR(p) models, and (b) as per (4) one can include a certain number of quasi-periodic components determined by the order of complex roots admitted (C). These observations are interesting, because they allow one to investigate low frequency trends and quasi-periodic behavior in univariate time series such as the Devils Lake volumes where these features are of interest.

In classical linear, autoregressive modeling, the order, p, of the model is selected and fixed at some level p*, using a criteria of best fit, such as the Akaike Information Criteria (AIC). Uncertainty of the AR coefficients and its impacts on simulations are then assessed within this framework. Huerta and West take a rather different approach. They assume user specified upper bounds on C and R. For instance, we could assume that the upper bound on C is 5 based on an assumption that potentially an annual cycle, two quasi-periodic periodic components with periodicities in the 3 to 5 year range related to ENSO, a quasi-periodic component with a decadal time scale related to NAO, and a quasi-periodic component at inter-decadal scales related to the PDO, were the main components of the climate system that are likely to be seen in the Devils Lake fluctuations. The actual frequencies ;#61559;j and their amplitudes rj for these components are not specified. Upper and lower bounds (typically 2;;#61548;;n/2, where ;#61548;=;#61490;;#61552;;#61487;;#61559;) on the frequencies are specified, and the amplitudes are bounded in absolute value by 1. The number of real roots R, could be chosen as a suitably large number, say 10. This would imply an upper bound on p of 20 (2C+R). In the absence of any information as to the number of admissible C, R, components, one may consider larger values for the upper bounds. Instead of seeking the “optimal” values for the model order and the associated model coefficients, Huerta and West use a Bayesian approach in which a prior probability distribution (typically approximately uniform, that admits mass on unit and zero roots) is specified on the values of the model parameters, and the data is used to then develop a posterior distribution for these parameters. Admitting zero roots allows for consideration of model order uncertainty, while admitting unit roots allows nonstationary components. A Markov Chain Monte Carlo (MCMC) approach is used to identify superior model coefficients. At the end of the MCMC simulations, contingent on the upper bounds for C and R, we have posterior probability distributions associated with each of the amplitudes rj for both the complex and the real roots. Thus, one could diagnose roots for which the posterior probability mass is significantly away from 0, and hence there is useful information. The corresponding frequencies ;#61559;j for such complex roots can then be highlighted. Probabilistic model forecasts for xt+T consequently encode the posterior probability distributions for each of the admissible parameter values.

The ARCOMP algorithm was applied to the Devils Lake monthly volume data and up to 10,000 simulations were generated from different starting points using the posterior probability density functions of the AR components. The Devils Lake data was transformed using a square root transform (as also done by Wiche and Vecchia, 1995) prior to the ARCOMP analysis. Exceedence probabilities of different DL levels were computed from the simulations from the end of the model-fitting period to enable comparisons with the USGS model forecasts for 2000 and beyond. Split sample analyses were used to assess the performance of ARCOMP on historical lake fluctuations using only prior period data. Various specifications of C and R, for the same maximal order p were considered. The forecasts from these combinations were generally quite comparable. The main difference is in the interpretation of the different models. When complex roots are allowed, the modes of the posterior probability distributions are predominantly at periodicities of 80, 36, 26, 18 and 12 months. Typically, the posterior distribution for the number of real processes (R) has a strong mode at 2, with some mass on 3 and 4. This suggests a fairly low order autoregressive model with some quasi-periodic components. However, the posterior modes for the dominant real roots, indicate a highly nonstationary model, particularly if data from the 1990’s is included.

Split sample results are presented here (Figure 20a) for forecasts for the model calibrated to data prior to mid-1994. One notes, that within 2 years (i.e., 1996) from the start of the forecast the Devils Lake trajectory has already gone above the 0.2% exceedence level threshold indicated from ARCOMP simulations. Note also, that the median forecast from ARCOMP is the maintenance of the 1994 volume. This feature persists for any forecast started after 1994, including one from 1999 conditions. Somewhat different forecasts are realized for the case where only real roots are admitted. The lake now tends to its average state in the long run, as with the Wiche and Vecchia forecasts, but the extreme exceedence probabilities still indicate that the lake will spill over for the entire future projection. Given the posterior density estimates of the AR model structure and coefficients, and the forecasting results, we have concluded that the primary information that can be gleaned from this analysis is that a linear autoregressive structure (even allowing for high order models to approximate quasi-periodic and low frequency trends as exhibited by some state space and ARMA models) is likely inappropriate for describing the Devils Lake monthly volume process. The indicated uncertainty bands are very wide and perturbations in the assumptions translate into quite different median forecasts. The nonstationarity of the system is emphasized by the near unit root solutions to the parameter vectors. The Wiche and Vecchia decomposition of the lake volume data into an annual mean and difference process is a clever decomposition of the data that may deal with some aspects of the nonstationary data better than the direct application. One can think that the ARCOMP procedure essentially generalizes the Wiche and Vecchia procedure by considering model and parameter uncertainty in an ARMA modeling framework. Given the limited success of the Wiche and Vecchia procedure, and the indications from the direct application of ARCOMP to the raw DL data, we conclude that there is a need for developing new modeling tools that can represent such apparently nonstationary processes from relatively short data sets. Interestingly, split sample experiments with ARCOMP using data prior to the recent rise (e.g., up to 1980), were better able to forecast the next several years. This suggests that the 1990’s regime shift from prior conditions is not adequately modeled by a linear model.

Figure 20. ARCOMP probabilistic forecasts of lake level exceedence based on 1,000 to 10,000 simulations from models using monthly Devils Lake volume data up to (a) 1994, (b) 1999, and (c) 1999 considering only real roots.

Summary

The main results from our investigations are:

1)The recent trend in the Devils Lake volume is likely a consequence of changes in the seasonality of annual rainfall, and may be determined to a great degree by increases in summer and fall precipitation, that are associated with corresponding changes in the atmospheric circulation for those seasons, that are manifested as decreases in the regional atmospheric pressure. The seasonality of the factors at play in the recent period is different than for the Great Salt Lake and the Devils Lake in the 1980s. Winter/spring/fall precipitation changes were a more dominant influence in that period, but for a shorter duration.

2)These changes in atmospheric circulation and regional precipitation have large spatial structure and are not likely due to increases in local convection and moisture recycling related to local conditions. There is some evidence that a combination of factors related to Pacific and Atlantic ocean-atmosphere oscillations is important. The summer-fall precipitation in the larger region that exhibits the consistent precipitation anomaly structure for the wet and dry periods for Devils Lake is influenced to some extent by features such as the night-time low-level jet, the Southwest monsoon, and northern frontal systems that bring ocean moisture to the region. These transient features are not directly reflected in the monthly atmospheric data that we analyzed, and had access to. Consequently, the correlative analyses used with climate indices and atmospheric pressure time series are inferential and diagnostic, rather than causal. It may be useful to pursue more direct investigations to better pin down the climate mechanisms responsible.

3)The spectral signature of the series analyzed reveals that while there are well separated, narrow band interannual, and interdecadal oscillatory components shared between the Devils Lake and the climate indices, their expression is rather time dependent and the recent record of the Devils Lake is manifested as a singularity in the system where the dominant frequencies of interest, both for the lake and for the climate indices (NINO3, PDO, and NAO) are concurrently at anomalous levels, their interaction (i.e., cross-ocean factors) is important in determining the local precipitation and lake response. This combination of factors and the lakes state does not have an analog in the 1905-1999 record.

4)The connection of the Devils Lake fluctuations to the climate indices raises hopes of climate based forecasts of the lake and associated regional hydroclimatology. However, given that no successful means of forecasting these indices for the long run yet exists, despite their clear-cut low-frequency character, makes the task difficult. An additional complication is that the Devils Lake basin hydrology and lake dynamics appear to be acting as a nonlinear, notch filter on the climate system, in that until some thresholds of wetness (unknown but definable) are crossed the lake volume squelches the climate signal. The changes in seasonality of the precipitation forcing may play a critical role in determining the nature and magnitude of this threshold behavior. Conversely, after the threshold is crossed, the climate signal is amplified. This is similar to but more dramatic than what happens with the Great Salt Lake.

5)A question that has been brought up in the climatic context of Devils Lake has been the possibility that a changed climate due to increased carbon dioxide (CO2) in the atmosphere may be responsible for the changes in the precipitation and in the lake volume. Such questions are invariably difficult to answer given the limitations of numerical models of the Earths climate and the limited time history over which such assessments can be done. We did not directly pursue investigations to investigate such an attribution. However, given the longer, paleoclimatic context for the region and for other lakes such as the Great Salt Lake, it is evident that the type of conditions being currently experienced have occurred in the past (see for instance the marker (X) in Figure 1) prior to our notion of anthropogenic climate change. Consequently, such questions can be answered in a useful way only through investigation of climatic mechanisms, i.e. modes of the ocean-atmosphere system, that would lead to anomalous moisture transport to the region, and to investigate whether the frequency of such modes is likely to undergo changes over time, in particular due to anthropogenic forcing.If changes in the frequency of such events are indicated, then the relative risk of such occurrences is likely to increase. We noted that the regime residence time and regularity/duration of switching of the low-frequency climate conditions indicated by the climate indices used have varied quite a bit over the historical period. Whether such variations occur in the natural climate or whether they are forced by greenhouse effects is difficult to diagnose, given that current coupled ocean-atmosphere models do not adequately reproduce these low frequency modes. However there are indications from several such models of the increased incidence of El Nino like conditions under a warming scenario, which may in turn translate into positive summer/fall precipitation in the region as indicated by the correlations identified here. However, the models are unable to define the nature of the PDO/NAO variations that have longer time scales and may be just as important for the region. Indeed, the persistent nature of the current event would likely be linked to the more slowly varying ocean states (PDO, NAO) than the tropical Pacific (El Nino).

6)Various methods of near-term and long-term forecasting were tested. In particular, the nonlinear time series methods based on state space reconstruction using Multivariate Adaptive Regression Splines (MARS) that have been very successful for Great Salt Lake forecasts, were used for the 1-4 year ahead forecasts. While for certain combinations of parameters, and at certain times, these models produced forecasts that looked quite promising, there was significant variation in the actual predictors selected and the complexity of the fitted models over the combinations explored. Further, many of the fitted models were unstable, in that small perturbations led to dramatically varying and unrealistic forecasts. The indicated model confidence limits, accounting for the degrees of freedom were often very wide. Our confidence in these forecasts is consequently limited. In this respect the performance was quite different from our applications to the Great Salt Lake volume data. Various differences in the character of the Devils Lake and Great Salt Lake volume fluctuations were highlighted earlier. The longer Great Salt Lake record actually included analogs for the extreme 1980s rise and fall for the lake, and consequently, the forecasts for that period were more robust. The Great Salt Lake forecasts for the subsequent period (to date) have been remarkably accurate even as far as 4 years ahead.

7)The longer run forecasts based on ARFIMA or ARCOMP type of models suffered from similar problems. Confidence limits were very wide, particularly as the 1990s data was included in these models, and the models indicated solutions on the non-stationary boundary. This is hardly surprising in light of the earlier results, and the anomalous nature of the recent period relative to the rest of the record. Spectral forecasting methods that may have been able to account for such nonstationarities were considered, but their direct application to the problem was not pursued. We were discouraged by the pronounced amplitude variation of the signals for Devils Lake over the period of record, and were concerned that robust forecasts from these methods would not be likely.

8)Given the historical perspective of Devils Lake and the Great Salt Lake, one can say that while there is considerable uncertainty associated with numerical forecasts of Devils Lakes future levels, such conditions or excursions for closed basin lakes are not uncommon. They reflect the complex interaction between basin scale and hemispheric scale threshold processes. Key characteristics of such excursions are that, once the threshold is crossed in either direction, the lake will rise or fall rapidly. The rate and amount of such changes are generally far larger than what strategies for human adaptation, such as pumping or diversions, can handle. Thus, flood plain zoning or other measures that moves critical facilities out of such areas are desirable. A more intensive regional/continental analysis of climatic factors and their land surface projection using model, historical and paleo information may be useful for improving the perception and estimates of risk.

References

Abarbanel, H. D. I., 1995. Analysis of observed chaotic data, Springer Verlag, New York, 272 pp.

Ames, D.P., 1998. Seasonal to interannual streamflow forecasts using nonlinear time series methods and climate information, Utah State University, M.S. Thesis.

Angel, J. R., and F. A. Huff, 1997. Changes in heavy rainfall in midwestern United States, J. Water Resour. Plng. and Mgmt., 123, 246-249.

Baldwin, C.K., 1999. Dynamics of streamflow generation: The upper Mississippi River, Utah State University, M.S. Thesis, 147 p.

Baldwin, C.K., and U. Lall, 1999. Seasonality of Streamflow: the Upper Mississippi River, Water Resource Research, 35(4) p. 1143-1154.

Baldwin, C.K., U. Lall, and F.H. Wagner, 1999. Trend and Variability Analyses for the Rocky Mountain/Great Basin Regional Assessment of Climate Variability and Change, Poster presented at AGU Spring meeting, Boston, MA. June1-4. Available online at: http://pub.uwrl.usu.edu/cbald/trends-poster.pdf

Battisti, D. S., and E. S. Sarachik, 1995. U.S. national report to International Union of Geodesy and Geophysics: Contributions in oceanography, American Geophysical Union, Reviews of Geophysics, 33, 1367-1376.

Corradi, V., 1995. Nonlinear transformations of integrated time series: A reconsideration, J. Time Series Anal., 16, 6, 539-549.

Cutler, C., and D. T. Kaplan, 1997. Nonlinear Dynamics and Time Series, American Mathematical Society, Providence, R.I., 252 pp.

Diks, C., 1999. Nonlinear Time Series Analysis,World Scientific, Singapore, 209 pp.

Freidman, J.H., 1991. Multivariate adaptive regression splines, Annals of Statistics, 19, 1-142.

Huerta, G. and M. West, 1999. Priors and component structures in autoregressive time series models. Journal of the Royal Statistical Society (Series B), 61, 881-899.

Jensen, Ray E. No Date. Climate of North Dakota. National Weather Service, North Dakota State University, Fargo, North Dakota. Jamestown, ND: Northern Prairie Wildlife Research Center Home Page.http://www.npwrc.org/resource/othrdata/climate/climate.htm (Version 02APR98).

Kantz, H., and T. Schreiber, 1997, Nonlinear Time Series Analysis, Cambridge University Press, Cambridge, U.K., 304 pp.

Karl, T. R., R. W. Knight, and N. Plummer, 1998. Secular trends of precipitation amount, frequency, and intensity in the United States, Bull. Amer. Meteor. Soc., 79, 231-241.

Kushnir Y., 1994. Interdecadal variations in North Atlantic sea surface temperature and associated atmospheric conditions, J.Clim., 7, 141-157.

Lall, U., T. Sangoyomi, and H.D.I. Abarbanel, 1995. Nonlinear Dynamics of the Great Salt Lake: Nonparametric Short-Term Forecasting, Water Resour. Res., Vol. 32, No. 4, p. 975-985.

Lee and Heghinian, 1977. A shift of the mean level in a sequence of independent normal random variables: A Bayesian approach, Technometrics 19(4), 503-506.

Lorenz, E. N., 1963. Deterministic nonperiodic flow, J. Atmos. Sci., 20, 130-141.

Mantua, N. J., S. R. Hare, Y. Zhang, J. M. Wallace, and R. C. Francis, 1997. A Pacific interdecadal climate oscillation with impacts on salmon production, Bull. Amer. Meteor. Soc., 78, 1069-1079.

Moon, Y-I., 1995. Large-scale Atmosphere Variability and the Great Salt Lake, Ph.D. dissertation, Utah State Univ., Logan, Utah.

Osborne, L., 1998. Variables Affecting Climate in the Devils Lake Basin, Regional Weather Information Center, John D.

Miscellaneous