[report] Nonlinear Dynamics Of The Great Salt Lake: Short Term Forecasting And Probabilistic Cost Estimation

E441, Fall 2008

Nonlinear Dynamics of the Great Salt Lake: Short Term Forecasting and Probabilistic Cost Estimation Cameron W. Bracken Department of Environmental Resources Engineering, Humboldt State University

A methodology has been presented which uses ideas from chaos theory to reconstruct and subsequently forecast the water surface level in the Great Salt Lake. Locally weighted polynomial (LWP) models are used to model the phase space. This methodology includes the generation of ensemble forecasts from a suite of models with different parameter combinations. The generalized cross validation statistic is minimized by varying the embedding dimension d E , characteristic time lag T, the local polynomial degree p and, neighborhood size α. Techniques from attractor reconstruction theory are used to select appropriate search ranges for the parameters. The utility of this method is shown in its ability to accurately blind forecast within the expected forecast time horizon (1 year). In the blind forecast beginning in 1985, the method accurately blind forecasts the observed peak. The ability to generate ensemble forecasts enables the generation of probabilistic cost estimates at every future time step. This type of information enables policy makers to asses risk and make informed decisions. Abstract.

et al., 1979]. Of particular interest is the sharp rise in water level from 1983–1987. Short and mid term forecasts are important for property and business owners near the GSL. In 1986 the Great Salt Lake Pumping project was initiated, for which accurate forecasts would have been particularly useful [UDNR, 2007]. By producing ensemble forecasts, a probabilistic estimate of monetary damage may be made. Sangoyomi [1993] compiled a biweekly record of GSL volume (ac-ft) beginning in 1848. Original measured data is stage (feet above mean sea level) from which volume is calculated via a stage-volume relationship such as the the one given in [James et al., 1979]. The variable used here will be standardized stage which is identical to the standardized volume (Figure 3).

1. Introduction Theoretically, any natural system could be modeled deterministically with physical equations derived from first principals. In practice, this is not possible due to lack of data, lack of precision and accuracy in the data itself, lack of computing resources and possibly a lack of physical equations which adequately describe the system. Additionally, Lorenz [1963] described a type of system, known as dynamical or chaotic, in which long term prediction is inherently impossible. In these types of systems, a simple deterministic relationships may appear as randomness. In fact, a high order dynamical system appears as complete randomness [Abarbanel et al., 1993]. Traditional statistical methods explain a time series as a realization of a random process resulting from a system with many degree of freedom [Regonda et al., 2005]. An alternate explanation is that the system governed may be governed by low order deterministic chaos. In general a chaos is (1)“aperiodic long-term behavior in a (2)deterministic system that exhibits (3)sensitive dependence on initial conditions” [Strogatz , 1994]. A consequence of (3) is that long term forecasting of these systems is inherently impossible. Short and mid term forecasting is still possible and may be improved upon by incorporating concepts from chaos theory. The Great Salt Lake (GSL) is a terminal or closed basin lake located in north western Utah between 40◦ –42◦ N and 110◦ –112◦ W (Figure 1). A terminal lake is one in which no surface outflow occurs. That is, the only hydrologic losses occur from evaporation and subsurface flow. The GSL is so named because of its high salt content which is typical of terminal lakes. Terminal lakes with large drainage areas (such as the GSL) effectively filter out high frequency atmospheric fluctuations and tend to be strongly influenced by long term atmospheric and climatic fluctuations [Arbarbanel and Lall , 1996]. This gives confidence in the to the theory that the system may be described a low order dynamical system. The Great Salt Lake stage has fluctuated wildly in the past 150 years of historical record. Sharp rises and drops in lake stage are of particular interest to humans because of potential for property damage or loss of business [James

2. Literature Review Chaos or dynamics in a deterministic system was first described by Lorenz [1963] where the famous Lorenz equations were derived. Aside from revealing that fairly simple systems can exhibit chaotic behavior, this discovery revealed the possibility that long term prediction of chaotic physical systems may be inherently impossible. The field of dynamics has evolved significantly finding broad applications in many fields including Chemistry, Biology, Medicine and Engineering [Strogatz , 1994]. One remarkable discovery in the field of dynamics is the technique of attractor reconstruction. Lorenz himself described the technique as the most surprising development in nonlinear dynamics [Strogatz , 1994]. Embedding theory, also known as attractor reconstruction, describes how the dynamics of a system can be recovered by samples taken from a single state (phase) space coordinate Takens [1981]. This remarkable technique has become a valuable tool in time series modeling [Abarbanel et al., 1993]. These methods involve reconstructing the chaotic attractor in dimension d with some characteristic time lag T. Models are constructed of the d dimensional phase (state) space which describe how a system evolves from an initial point x. The form of the model is for I time steps in the future is xt+ I = f ( x I , xt− I , xt−2I , ..., xt+(d−1) I );

f : Rd → R

(1)

where f may be local or global. In fact if I = 1 and f is a linear function then the model is a linear autoregressive

Engineering 441, Fall 2008.

1

2

Bracken: Nonlinear Dynamics of the Great Salt Lake

(AR) model [Lall et al., 1996]. The model can be thought of as a nonlinear autoregressive process. Attractor reconstruction was first used by Packard et al. [1980] to reconstruct the time derivatives of a time series, st = st0 +nτs where τs is the sampling interval and n = 1, ..., N. The theory behind this technique was first described mathematically by Takens [1981]. The overview of this method provided by Abarbanel et al. [1993] is strongly recommended. Attractor reconstruction has been used in forecasting of natural systems [Bordignon and Lisi, 2000; Islam and Sivakumar , 2002; Coulibaly and Baldwinb, 2005], particularly with success in the case of the GSL [Sangoyomi, 1993; Abarbanel et al., 1995; Arbarbanel and Lall , 1996; Regonda et al., 2005; Lall et al., 2006]. Arbarbanel and Lall [1996] used a local linear model to reconstruct the phase space. Their single step forecasts were amazingly accurate though

Idaho

Great Nevada

Salt Lake

Wyo Utah

no confidence intervals on the forecast were given. The decrease in forecast skill was apparent using 10 step forecast, that is the forecasted value at every point is blind forecasted 10 days earlier. Regonda et al. [2005] focused on blind forecasts only. Arbarbanel and Lall [1996] calculated liapunov exponents to show that the forecast time horizon is about 1 year. Lall et al. [1996] modeled the phase space with multivariate adaptive regression splines and was able to show significant forecast skill over traditional linear AR models; AR lags of up to 73 were used. In many cases when using statistical models with parameters Θ, one is faced with competing models that, based on some criteria C(Θ), are very similar. This criteria may be a measure of goodness of fit or parsimony (Coefficient of determination, Mallows C p , Akaike information criterion) or predictive risk (Cross validation, generalized cross validation [Craven and Wahba, 1979]). Traditional regression would choose a single the “optimal” parameter set Θ∗ ignoring any other models. The multi-model method claims that a sufficiently different parameter set, say Θ0 with a similar criteria value, C(Θ0 ), contains valuable information about the system that is lost if it is not included in some overall or ensemble model. The multi-model method been used with success in the forecasting skill in many hydrologic systems [Krishnamurti et al., 1999, 2000; Rajagopalan et al., 2002; Grantz et al., 2005; Hagedorn et al., 2005; Regonda et al., 2005, 2006]. Here we use the forecast algorithm of Regonda et al. [2005] which uses a locally weighted polynomial (LWP) model. LWP models have been applied successfully in many fields of hydrologic modeling including ensemble stream flow forecasting and simulation [Grantz et al., 2005; Regonda et al., 2006; Prairie, 2006; Opitz-Stapleton et al., 2007] rainfall and climate forecasting [Krishnamurti et al., 2000; Singhrattna et al., 2005] and the GSL itself [Regonda et al., 2005; Lall et al., 2006]. One of the benefits of local regression is its ability to capture arbitrary local nonlinearities which are expected to be prevalent in this setting.

3. Methodology N

Outline of Great Salt Lake Basin

Figure 1. Great Salt Lake and drainage basin (Adapted from James et al. [1979]).

This section contains the details of the methods outlined in Sections 1 and 2. Consider a dynamical system with an attractor dimension d A from which the time series sn = st0 +nτs has been sampled at the interval τs . According to the theory, the dynamics of a system can be reconstructed by creating a d E dimensional vector of lagged or embedded time series. y(n) = [st0 +nτs , st0 +(n+T )τs , st0 +(n+2T )τs , ..., st0 +(n+(dE −1)T )τs ]

(2)

Stage (ft. above MSL) 4195 4205

= [sn , sn+T , sn+2T , ..., sn+(dE −1)T ]

1880

1920

1960

Figure 2. Standardized Great Salt Lake Stage.

2000

where T is an unknown integer and n = 1, ..., N, known as the characteristic lag Abarbanel et al. [1993]. Viewing these embedded time series as the coordinates of a dynamical system is the basis of attractor reconstruction. The embedding theory of Takens [1981] guarantees that a d A dimensional attractor can be reconstructed (unfolded) unambiguously from s(n) in a dimension d E > 2d A . The parameter d E is said to be the embedding dimension of s(n). That is any point xi in the phase space lies on one or zero trajectories of the dynamical system. The notion that no two trajectories can cross in phase space is an important result from chaos theory [Strogatz , 1994]. Reconstructing the phase space in this fashion the evolution of the system to be examined from any initial condition x0 . In practice there are actually a large number of embedded time series that depend on how much data is available prior

3

Bracken: Nonlinear Dynamics of the Great Salt Lake the time of the forecast. If we want to forecast from time step I onward then there are n E embedded time series where n E = I − (d E − 1) T − 2. First, placing each embedded time series in a vector sn+T = [sn , sn+T , sn+2T , ..., sn+(dE −1)T ]. All the embedded time series are gathered into an n E × d E matrix S where

and Wahba, 1979] and is indirectly dependent on all of the parameters of our model. The GCV is 1 N

N

∑ ei2

i =1 (6) 2 (1 − m N) where Θ = [d E , T, α, p] is a vector of the parameters, ei is   s1+ T the model residual and m is the number of parameters.   s2+ T Once constructed for a given parameter set Θ, S and r   S= = .. do not change though they will be different for any param  . eter combination of which there are quite a lot. Following s I −1−(dE −1)T −1 the multi-model method all models with GCV values within   10% of the top model are included and used in the ensems1 s1+ T ··· s1+(dE −1)T   ble forecast. This suite of models is then used to obtain s2 s2+ T ··· s2+(dE −1)T   an ensemble of predictions. This ensemble of predictions  . .. .. .. ..   also gives an idea of the forecast uncertainty, i.e. the wider   . . . . variation in ensemble members, the more uncertain the pres I −2−(dE −1) s I −2−(dE −1)T +T · · · s I −2−(dE −1)T +(dE −1)T diction. (3) The model for predicting K time steps into the future has the form The bottom rightmost element has been written out explics I +K = f (z( I )) + ε I (7) itly but simplifies to s( I − 2), that is two time steps before where K = 0 is considered the first future point and z( I ) the first forecast time. Generally each i, j entry of P is given is called the feature vector which is simply the most recent by embedded time series ending at I − 1 where Si,j+1 = si+ jT (4)   s I −(dE −1)T −1 where i = 1, ...n E and j = 0, 1, ..., d E − 1. In addition to  s I −(d −2)T −1  recording the conditions in the phase space in the the maE   z( I ) =  (8) . .. trix S the next successive lagged value of the time series is   . recorded in a vector r of length n E where s I −1   s1+(dE −1)T +1   For every new time step into the future z( I ) is reconstructed s2+(dE −1)T +1     r= (5) by including the latest forecasted point at the end of the ob..  served time series and repeating as if the observed series   . was one point longer. Forecasts are assumed to become sucs I −2−(dE −1)T +(dE −1)T +1 cessively worse for larger K’s as more forecasted points are The last entry in r simplifies to s I −1 which is the last time included. By preforming the forecast with each model in before the start the forecast. The vector r is the key to re- the ensemble, an estimate of the forecast uncertainty can be lating a state of the system in the phase space to a value of obtained by observing the spread in the ensemble members. the time series. A pseudo code for constructing the phase space calibra- 3.2. Probabilistic Damage Estimates tion model looks like: Rising water levels cause damage to businesses and propDO i from 0 to n E − 1 erty which reside on the GSL lakefront. In the case of the GSL pumping project costs are related to increased pumpDO j from 1 to d E ing. From such information a stage versus cost relationship Si,j = si+ j∗T can be developed. Here a hypothetical relationship is used to r i = s d E ∗ T + i +1 END DO END DO For even further clarification a numeric example from the GSL. Say the forecast time step is I = 3000 and sn have data starts from n = 1 and d E = 4 and T = 15. There are n E = I − 2 − (d E = 1) T = 3000 − 2 − (4 − 1) · 15 = 2953 embedded time series! It should be mentioned that this method is computationally intensive, especially with large historical records.

Space

As mentioned in Section 2 local polynomial models (LWP) are used in this analysis to fit the phase space. These models are have to a local polynomial degree p and are fitted to k nearest neighbors in the phase space where k = αN, α ∈ (0, 1] and N is the number of data points. If α = p = 1 the local regression collapses to multiple linear regression and thus is a generalization which is is more flexible. The local regression is preformed with the LOCFIT package [Loader , 1999]. The generalized cross validation (GCV) statistic is used here as a measure of predictive risk [Craven

Cost (Million $) 2.345

Phase

2.335

3.1. Locally Weighted Polynomial Model and Ensemble Generation

GCV (Θ) =

4180

4190 4200 Stage (ft. above MSL)

Figure 3. Hypothetical stage-cost curve.

4210

4

Bracken: Nonlinear Dynamics of the Great Salt Lake

AMI ( T ) −

N

∑

n =1

P(sn , sn+T ) log2



P(sn , sn+ T ) P(sn ) P(sn+ T )

 (9)

2

3

The forecast algorithm as a whole proceeds as follows: 1. Using the appropriate methods (ACF, AMI, FNN) obtain a range of values for d E and T 2. Construct S and r for a particular combination of α, p, de , and T. 3. Fit a LWP model to S and r and calculate the GCV value for the particular fit. 4. Repeat steps 1-2 for all combinations of α, p, de , and T. 5. Collect the parameter values for all models within 10% of the model with the lowest GCV value.

-2

where P(·) denotes probability, P(·, ·) denotes joint probability and AMI ( T ) ≥ 0 [Abarbanel et al., 1993]. Abarbanel et al. [1993] suggests that the first local minimum of AMI ( T ) gives a good estimate of T.

3.5. Forecast Algorithm

1

Until now the lag T has been an unknown integer, though the choice of its value is important. One of the benefits of the multi-model ensemble method is that the pressure is taken off of choosing a single best value; the choice of an acceptable range is nonetheless important. If T is chosen too small, each point in phase space will be essentially the same as the last and if T is chosen too large then information about the system will be lost. Two common methods are using the linear autocorrelation function (ACF) and the average mutual information criteria. A somewhat naive approach to determining T is to examine is linear dependence of s(n) and s(n + T ) form the ACF. One would examine the ACF as a function of T and find the first zero crossing. Unfortunately this method does not take into account the nonlinear dependence of the variables and leads to a gross overestimation of T when used in the traditional way [Abarbanel et al., 1993]. Fortunately the ACF may still be used to estimate T by finding its first local minimum [Arbarbanel and Lall , 1996]. The average mutual information (AMI) criteria in general is used to describe the amount of information one gains on average about an observation ai from an observation bk where ai and bk are measurements from systems A and B respectively. In terms of the time series s(n) the AMI criteria describes the about of information gained about sn+T having observed sn . The AMI is

e(t + 2T ) 0

3.3. Methods for Choosing Time Delay

of the linear part of the correlation integral curve ceases to change as the dimension increases. This corresponds the the attractor being completely unfolded.

-1

demonstrate the utility of ensemble forecasts for generating cost estimates.

-2

(10)

therefore the dimension is the slope of the line log C (e) = a + d log e.

1 e(t)

2

3

0.9

1.0

Figure 4. Reconstructed GSL attractor viewed in R3 (i.e. d E = 3), unfolded with T = 4.

0.5

C (e) ∝ ed ,

0

ACF 0.7 0.8

Embedding theory guarantees that an embedding dimension d E > 2d A is sufficient to unfold the attractor but not necessary. In practice a d E ≤ 2d A may sufficiently unfold the attractor. A reasons for not simply choosing a large embedding dimension are (1) that computational costs become larger and (2) relationships between variables in higher dimensions may simply be noise in observed data [Abarbanel et al., 1993]. Two common techniques for determining the embedding dimension are the false nearest neighbor (FNN) technique and the saturation of system invariants or the correlation integral method. The false nearest neighbor (FNN) method attempts by geometrical criteria to determine weather two points in dimension d are neighbors due to the projection of the attractor into dimension that is insufficiently low. The correlation integral is another geometric method which starts by placing a ball of radius e enclosing and average fraction C (e) of the points in dimension d. The number of points enclosed is expected to grow like

-1

0.6

3.4. Methods for Choosing Time Delay

3 1 2 0 -2 -1 e(t + T )

(11)

In practice this linearity only holds for some middle region of the range of e because the ball eventually encloses all of the points [Strogatz , 1994]. The dimension in which the slope

0

50

100 T

Figure 5. Autocorrelation function.

150

200

5

Bracken: Nonlinear Dynamics of the Great Salt Lake 6. Form an ensemble by making a forecasts at time step I from each of the models for K time steps into the future.

4. Application

5. Results Blind forecasts beginning in 1979, before the sharp rise in the GSL water level and for 1985, after 2 years of dramatic water level rise and just before the implementation of the GSL pumping project. In addition single step forecasts are presented beginning at both times. And an example of a probabilistic cost estimate.

1.0

AMI ( T ) 1.5

2.0

2.5

Fraction of false nearest neighbors 0.6 0.4 0.5 0.7

The GSL water level time series was compiled by Sangoyomi [1993]. The analysis is done on the standardized time series with a mean 4201 feet and a standard deviation 4.4 feet. A reconstructed attractor with d E = 3 is shown in figure 4. The autocorrelation function (ACF) is shown in figure 5. The first local minimum occurs at T = 15. As a side note the criteria of a linear ACF would result in a time lag T > 400! The Average mutual information (AMI) function is shown in figure 6. The first local minimum occurs at T = 16 which is in good agreement with the lag suggested form the ACF. This leads to a search range for T of 13–19. The correlation integrals for the GSL time series are shown in Figure 7. The slope of linear part is independent

of dimension at around d E = 4. Figure 8 shows the fraction of the false nearest neighbors as a function of embedding dimension. The FNN method suggests an embedding dimension d E = 4 which agree with the observation of the correlation integrals. This prompts the search for d E in the range 3–5.

0.5

1 0

50

100 T

150

200

2

4 3 5 Embedding dimension

6

7

Figure 8. Correlation Integrals for d E = 1, ..., 7.

Stage (ft. above MSL) 4196 4197 4198 4199

0.6

1.0

4200

Figure 6. Average Mutual information criteria.

C (e) 0.4

dE = 1 dE = 2 dE = 3 dE = 4

4195

0.2

dE = 5 dE = 6 dE = 7 0.5

1.0

1979 2.0

5.0 e

Figure 7. Correlation Integrals for d E = 1, ..., 7.

10.0

1980

1981

1982

1983

Figure 9. Blind forecast beginning January 1st, 1979. The dashed line is the median forecast and the dotted lines are the 5th and 95th percentiles of the forecast. Circles represent measured data.

6

Bracken: Nonlinear Dynamics of the Great Salt Lake For this forecast 25 models were selected within 10% of the top model. Surprisingly enough the parameter ranges were identical to those selected for the 1979 blind forecast. These ranges agree with those reported by Regonda et al. [2006]. The one step forecast does exceptionally well again for this time interval (Figure 13). Figure 12 is an example of the probabilistic cost forecast that can be produced at any time future time step. Herein lies the utility in an ensemble forecast. This type of information allows policy makers to asses risk and make decisions accordingly.

6. Conclusion A methodology has been presented which uses ideas from chaos theory to reconstruct and subsequently forecast the water surface level in the Great Salt Lake. Locally weighted

0e+00

4196

Stage (ft. above MSL) 4198 4199 4197

Probability Density 2e-05 4e-05 6e-05 8e-05

4200

Recall the estimated forecast time horizon was 1 year. The blind forecast in Figure 9 starting in 1979 verifies this. The blind forecast captures the observations well until after almost exactly a year the observed data begins to veer off from the forecast. For this forecast 27 models were selected within 10% of the top model. The parameter ranges for embedding dimension was 4–5, time lag 14–15 (units of 15 days), local polynomial degree 1–2, and neighborhood size .1–.6. The one step forecast shown in Figure 10 does amazingly well, though the utility of a 15 day forecast is limited. The Blind forecast shown in Figure 11 beginning in 1985 does well at capturing the peak and not overshooting. The most important piece of information at that point in time was weather the lake would continue to rise at its current rate or start to decline. The medina of the ensemble forecasts would have forecasted that the peak lake level to go no higher than it actually did over two years in advance.

5020 1979

1980

1981

1982

1983

5120

Figure 12. 1985 probability density function at 50 steps ahead.

4202

4206

4204

Stage (ft. above MSL) 4210 4206 4208

Stage (ft. above MSL) 4207 4208 4209 4210

4212

4211

Figure 10. Same as figure 9 but one step forecast beginning January 1st, 1979.

5040 5060 5080 5100 Cost (Ten Thousand $)

1985

1986

1987

1988

1989

Figure 11. Same as figure 9 but blind forecast beginning January 1st, 1985.

1985

1986

1987

1988

1989

Figure 13. Same as figure 10 but one step forecast beginning January 1st, 1985.

Bracken: Nonlinear Dynamics of the Great Salt Lake polynomial models are used to model the phase space. This methodology includes the generation of ensemble forecasts from a suite of models with different parameter combinations. The generalized cross validation statistic is minimized by varying the embedding dimension d E , characteristic time lag T, the local polynomial degree p and, neighborhood size α. Techniques from attractor reconstruction theory are used to select appropriate search ranges for the parameters. The utility of this method is shown in its ability to accurately blind forecast within the expected forecast time horizon (1 year). In the blind forecast beginning in 1985, the method accurately blind forecasts the observed peak. The ability to generate ensemble forecasts enables the generation of probabilistic cost estimates at every future time step. This type of information enables policy makers to asses risk and make informed decisions. One obvious downside of this type of model is its inability to generate long term series. Though a portion of the phase space and attractor within it are reconstructed, the models do not cover the entire phase space (granted that the reconstructed portion is of the most interest). As a result when generating long term series, the model fails as soon as a point is generated where there is no data. So, while this method is useful for short term forecasting and decision making, its uses in long term planning are limited. That is, unless a method to correct for the long term instability in the model is devised. On that note one way might be to use the median forecast of all the ensemble members as the update to all the models in the next time step. This would admittedly destroy the benefit of producing an ensemble forecast but could enable the generation of long term series by stabilizing the forecasts. Another extension might be to weight the forecasts of each model according to their relative GCV score.

References Abarbanel, H. D. I., R. Brown, J. J. Sidorowich, and L. S. Tsimring (1993), The analysis of observed chaotic data in physical systems, Reviews of Modern Physics, 65(4), 1331–1390. Abarbanel, H. D. I., U. Lall, Y.-I. Moon, M. E. Mann, and T. Sangoyomi (1995), Nonlinear dynamics and the Great Salt Sake: a predictable indicator of regional climate, Energy, 21, 655–665. Arbarbanel, H. D. I., and U. Lall (1996), Nonlinear dynamics of the Great Salt Lake: system identification and prediction, Climate Dynamics, 12, 287—297. Bordignon, S., and F. Lisi (2000), Nonlinear analysis and prediction of river ow time series, Environometrics, 11, 473–477. Coulibaly, P., and C. K. Baldwinb (2005), Nonstationary hydrological time series forecasting using nonlinear dynamic methods, Journal of Hydrology, 307, 164–174. Craven, P., and G. Wahba (1979), Smoothing noisy data with spline functions: estimating the correct degree of smoothing by the method of generalized cross-validation, Numerical Mathematics, 31, 377–403. Grantz, K., B. Rajagopalan, M. Clark, and E. Zagona (2005), A technique for incorporating large-scale climate information in basin-scale ensemble streamflow forecasts, Water Resources Research, 41. Hagedorn, R., F. J. Doblas-Reyes, and T. N. Palmer (2005), The rationale behind the success of multi-model ensembles in seasonal forecasting. Part I: Basic concept, Tellus, Ser. A, 57, 219–233.

7

Islam, M. N., and B. Sivakumar (2002), Characterization and ˘ dynamics: a nonlinear dynamical view, prediction of runo¨ıˇ nA Advances in Water Resources, 25, 179–190. James, L. D., D. S. Bowles, W. R. James, and R. V. Canfield (1979), Estimation of Water Surface Elevation Probabilites and Associated Damages for the Great Salt Lake, Tech. rep., Utah Water Researsh Laboratory, College of Engineering, Utah State University. Krishnamurti, T. N., C. M. Kishtawal, T. E. LaRow, D. R. Bachiochi, Z. Zhang, C. E. Williford, S. Gadgil, and S. Surendran (1999), Improved weather and seasonal climate forecasts from multi-model superensemble, Science, 285, 1548–1550. Krishnamurti, T. N., C. M. Kishtawal, Z. Zhang, T. E. LaRow, D. R. Bachiochi, C. E. Williford, S. Gadgil, and S. Surendran (2000), Multimodel ensemble forecasts for weather and seasonal climate, Journal of Climatology, 13, 4196–4216. Lall, U., T. Sangoyomi, and H. D. I. Abarbanel (1996), Nonlinear dynamics of the Great Salt Lake: Nonparametric short-term forecasting, Water Resources Research, 32, 975–985. Lall, U., Y.-I. Moon, H.-H. Kwon, and K. Bosworth (2006), Locally weighted polynomial regression: Parameter choice and application to forecasts of the Great Salt Lake, Water Resources Research, 42, 11 pp. Loader, C. (1999), Local Regression and Liklihood, Springer. Lorenz, E. N. (1963), Deterministic nonperiodic flow, Journal of the Atmospheric Sciences, 20, 130–141. Opitz-Stapleton, S., S. Gangopadhyay, and B. Rajagopalan (2007), Generating streamflow forecasts for the yakima river basin using large-scale climate predictors, Journal of Hydrology, 341, 131–143. Packard, N. H., J. P. Crutchfield, J. D. Farmer, and R. D. Shaw (1980), Geometry from a time series, Physics Review Letters, 45, 712–716. Prairie, J. (2006), Stochastic nonparametric framework for basin wide streamflow and salinity modeling: Application for the colorado river basin, Ph.D. thesis, University of Colorado. Rajagopalan, B., U. Lall, and S. Zebiak (2002), Optimal categorical climate forecasts through multiple GCM ensemble combination and regularization, Monthly Weather Review, 130, 1792–1811. Regonda, S. K., B. Rajagopalan, U. Lall, M. Clark, and Y.-I. Moon (2005), Local polynomial method for ensemble forecast of time series, Nonlinear Processes in Geophysics, 12, 1–10. Regonda, S. K., B. Rajagopalan, M. Clark, and E. Zagona (2006), A multimodel ensemble forecast framework: Application to spring seasonal flows in the gunnison river basin, Water Resources Research, 42. Sangoyomi, T. (1993), Climatic variability and dynamics of Great Salt Lake hydrology, Ph.D. thesis, Department of Civil and Environmental Engineering, Utah State Univ. Singhrattna, N., B. Rajagopalan, M. Clark, , and K. K. Kumar (2005), Forecasting thailand summer monsoon rainfall, International Journal of Climatology, 25, 649–664. Strogatz, S. H. (1994), Nonlinear Dynamics and Chaos, Westview Press. Takens, F. (1981), Detecting strange attractors in turbulence., Dynamical Systems and Turbulence, Lecture Notes in Mathematics, 898, 366–381. UDNR (2007), Great salt lake pumping project, Tech. rep., Utah Division of Natural Resources. C. W. Bracken, Environmental Resources Engineering, Humboldt State University, Arcata, CA, USA. ([email protected])