LongerLead WaterSupply Forecasts Statistical Forecasting with Optimal Climate

Longer-Lead Water-Supply Forecasts - Statistical Forecasting with Optimal Climate Predictor Selection Hamid Moradkhani Seasonal to two year forecasting workshop for the Colorado Basin March 21 -22 Department of Civil and Environmental Engineering 1

Background and Motivation Ø Interest in water supply forecasting has grown prominently in the western US due to population growth and increasing demands for water Ø In the Western US, the NWRFC and the NRCS jointly issue seasonal water supply outlook (WSO) forecasts of naturalized or unimpaired flow Ø Successful management of the West’s water supply is necessary to provide an uninterrupted and dependable water supply to meet all the needs Ø One important aspect of successfully managing the West’s supply of water is accurate and reliable forecasts of seasonal streamflow volumes Ø Longer-lead forecasts are useful to water managers and decision makers but difficult to make due to the uncertainty in future winter and spring climate conditions and the lack of snowpack information

Developing a Forecast Model l Model selection: Ø Statistical regression-based models Ø A regression model consists of l Dependent variable l Predictor variable(s) Ø A regression model established a linear relationship between the predictor variable(s) and the dependent variable (predictand)

Developing a Forecast Model l The dependent variable is the total volumetric flow over a particular period at a specific point in a basin. Total Volumetric Flow = Sum(Apr to Sept Volumes) Total volume is what we want to predict

Study Basins

Yakima River Basin l l The Yakima River is located in central Washington in the Yakima River Basin and is approximately 215 miles in length. The Yakima Basin drains approximately 6, 150 square miles of area. The basin is bordered on the west by the Cascade Mountains, on the north by the Wenatchee Mountains, on the east by the Columbia, and on the south by the Simcoe Mountains. The climate in most of the basin is dry, with a mean annual precipitation over the entire basin of 27 inches. Most of the precipitation in the basin falls during the winter months as snow in the mountains.

Sprague River Basin l The Sprague River is located in southwestern Oregon in the Upper Klamath River Basin and is approximately 75 miles in length. l The Sprague River drains an area east of the Cascade Mountains that is approximately 1, 600 square miles in area. l l The climate in most of the basin is much drier than that of western Oregon, and has more extreme temperatures, especially in the winter months. Snowfall accounts for 30 percent of the annual precipitation in the valleys and as much as 50 percent in the mountains.

Rogue River Basin l l The Rogue River is located in southwestern Oregon in the Rogue River Basin and is approximately 220 miles in length. The Rogue River drains an area between the Cascade Mountains and the Pacific Ocean that is approximately 5, 160 square miles. The climate of southwestern Oregon is cool and wet in the winter and among the hottest and driest in the western Cascades in the summer. Mean annual precipitation in the headwaters of the basin range from 20 to 30 inches, and 80 to 100 inches near the Oregon coast.

Study Basins Rogue River Basin location: Southwestern Oregon Basin size: 1, 600 square miles Basin location: Southwestern Oregon Basin size: 5, 160 square miles Streamflow (1000 -AF) MAP (in) Sprague River Basin Month Basin location: Central Washington Basin size: 6, 150 square miles MAP (in) Streamflow (1000 -AF) Yakima River Basin Month

Developing a Forecast Model l l Our objective is to establish a linear relationship between Ø several predictor variables (x), and Ø a predictand (y) Given a set of data that consists of n observations and j predictor variables the problem becomes finding the regression function Each of the j predictor variables has its own regression parameters, bj, and regression constant, b 0. The best possible estimates of the regression parameters are the ones which minimize the sum of the squared errors

Identifying Potential Predictors l Explore the relationship between streamflow and: Ø Snow water equivalent Ø Precipitation for past months Ø Streamflow for past months l Streamflow in the West is the result of the Ø Accumulation of seasonal snowpack over the winter months Ø Melting of this snowpack over the spring and summer

Identifying Potential Predictors Total Fall and Winter Precipitation = l l Fall and winter precipitation can also provide information about spring runoff This gives us information about the soil moisture state in the basin Sum(Oct to Forecast Issue Date)

Identifying Potential Predictors l l Past months streamflow also provides useful information about future streamflow volumes. Combining these three predictors would results in the following: Forecast Issue Date SWE Antecedent Precipitation Forecast Window Antecedent Streamflow Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr - Sep

Identifying Potential Predictors l l Recent research has shown that forecast accuracy can be improved by including large-scale climate information into forecasts models. Some useful climate teleconnection indices commonly used for streamflow forecasting in the western US: Ø Ø Ø Southern Oscillation Index (SOI) Pacific Decadal Oscillation (PDO) Multivariate El Nino Southern Oscillation Index (MEI) Pacific North American Index (PNA) Trans-Nino Index (TNI)

Identifying Potential Predictors l l l So how do we incorporate large-scale climate information into seasonal forecasts… Goal: Establish a relationship between climate information and spring runoff. Analysis: Investigate how climate information from the previous year relates to spring runoff of the forecast year.

Identifying Potential Predictors l Correlation analysis can be performed by taking monthly or 3 month aggregations of climate data and correlating with spring runoff.

Identifying Potential Predictors l Combining SWE, climate information from past months, precipitation from past months, and streamflow from past months Forecast Issue Date SWE 3 -Month Climate Data Antecedent Precipitation Forecast Window Antecedent Streamflow Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr - Sep

Identifying Potential Predictors l We have a set of potential predictors that capture some of the physical driving processes of streamflow, now how do we use them in a regression model? l Looking at the data we would find that our predictors are intercorrelated.

Intercorrelated Predictor Variables l l l If the predictor variables have minimal correlation among themselves, then application of the OLS regression methodology works as intended. If, on the other hand, the predictor variables are mutually correlated, then the predictor variables contain redundant information, and this condition leads to parameter estimates that are biased. Several remediation techniques are available for dealing with datasets that would otherwise give bias parameter estimations. Ø Ø Principal component regression Partial least squares regression Z-score regression Independent component regression

Principal Component Regression l l Principal component regression (PCR) combines principal component analysis (PCA) with multiple linear regression. PCA is a technique that creates new variables from the original predictor variables that can be scored and ordered. Principal components: Ø are composites of many of the original predictor variables Ø they represent a smaller number of variables that maintain the important patterns in the data. Principal components are generated from linear combinations of the predictor variables.

Principal Component Regression l Principal components are ordered Ø The first principal component accounts for as much of the variability in the data set as possible Ø Each succeeding component accounts for as much of the remaining variability as possible. l Most importantly, the principal components created from PCA are uncorrelated, which allows them to be used in multiple linear regression.

Principal Component Regression l PCR Steps: 1. Start by standardizing the predictor dataset – subtract the sample mean and divide by the sample standard deviation. 2. Compute the sample variance-covariance matrix of the standardized dataset. 3. Calculate the eigenvectors and eigenvalues of the variancecovariance matrix. 4. Calculate the principal component time series by multiplying each predictor variable times the corresponding eigenvector. 5. Use the principal components obtained from step 4 in linear regression model.

Principal Component Regression l l l PCA produces the same number of principal components as there are predictor variables. Some measure must be used to select the number of principal components to retain in the regression model. We investigated two methods: Ø The regression parameters are tested using a standard t-test and a user specified critical t-value Ø The Prediction REsidual Sum of Squares (PRESS) statistic is calculated using a leave one out cross-validation procedure

Principal Component Regression l l l Principal components are added into the regression model oneby-one The addition of each principal component is checked for statistical significance using the standard t-test. The principal components are retained in the regression model as long as the regression parameters are statistically significant (t-value > 1. 6).

Principal Component Regression l The PRESS statistic is calculated for each of the extracted factors j using: PRESS(j) = ∑(qobs – qsim(j))2 l Where Ø qobs observed streamflow value Ø qsim (j) predicted streamflow value for the extracted factor j l The PRESS statistic with the minimum value is then used to determine the number of extracted components to keep in the final model.

Partial Least Squares Regression l l l Partial least squares regression (PLSR) is based on the principal components of both the predictor (X) and the dependent variable (Y) PLSR is different than PCR in that PLSR searches for a set of components that explains the maximum covariance between X and Y, where PCR concentrates on the variance of X only PLSR decomposes X and Y into a score matrix (Sx, Sy) times a loading matrix (Lx, Ly) and a residual matrix (Ex, Ey): X = Sx * Lx + Ex Y = Sy * Ly + Ey l This is referred to as the outer relations

Partial Least Squares Regression l PLSR seeks to minimize Ey while maintaining the correlation between X and Y by an inner relations: X = Sx * Lx + Ex (1) Y = Sy * Ly + Ey (2) Sy = D * Sx + Ei (3) Outer relations Inner relations Where, D = diagonal correlation matrix between X and Y Ei = error term l Inserting equation (3) into equation (2) gives a predictive model for Y Y = D * Sx * Ly + Eyy where Eyy is to be minimized

Partial Least Squares Regression l l PLSR also produces the same number of components as there are predictor variables. In this study the PRESS statistic is used for determining the optimal number of components to retain PRESS(j) = ∑(qobs – qsim(j))2 Where Ø qobs observed streamflow value Ø qsim (j) predicted streamflow value for the extracted factor j

Development of Climate Predictors using ICA l l l This research investigated the use of Independent Component Analysis (ICA) for the decomposition of large-scale climate data in order to see if statistically significant climate predictors could be extracted. The idea is that the oceanic and atmospheric systems receive contributions from many sources and that the observations represent a linear mixture of independent signals. Before ICA was implemented a quick correlation analysis was performed using the NCEP Reanalysis data provided by the Physical Sciences Division NOAA/OAR/ESRL (http: //www. cdc. noaa. gov/).

Climate Signals Correlated With Spring Runoff Volumes Label Pattern PNA Pacific North American Index AMO Atlantic Multidecadal Oscillation (unsmoothed) EPNP Eastern Pacific Oscillation EA Index WP Western Pacific Index PE Polar/ Index NAO Oscillation SCAND Index SOI Southern Oscillation Index SOLAR FLUX Solar Flux (10. 7 cm) PDO Pacific Decadal Index OLR Outgoing Long Wave Radiation Equator QB 30 Quasi-Biennial Oscillation (30 mb zonal wind) SLP Sea Level Pressure QB 50 Quasi-Biennial Oscillation (50 mb zonal wind) SLP Sea Level Pressure GLAAM Globally Integrated Angular Momentum 850 MB TRADE WINDS 850 mb Trade Wind Index MEI Multivariate ENSO index 1000 MB ZWIND NCEP 1000 mb Zonal Wind NINO 1+2 Extreme Eastern Tropical Pacific SST KAPLAN SST Surface Temperature NINO 3 Eastern Tropical Pacific SST NCEP SST Surface Temperature NINO 3. 4 Eastern Central Tropical Pacific SST NCEP SUR PRESS NCEP Surface Pressure NINO 4 Central Tropical Pacific SST NCEP SLP NCEP Sea Level Pressure ONI Oceanic Nino Index NCEP AIR TEMP NCEP Surface Air Temperature TNI Trans-Nino Index NCEP 1000 MB GEO NCEP 1000 mb Geopotential Height WHWP warm pool NCEP REL HUMID NCEP Surface Relative Humidity TNA Tropical Index NCEP PRCP RATE NCEP Surface Precipitation Rate TSA Tropical Index NCEP PRCP WATER NCEP Surface Precipitable Water

Climate Signal Correlation With Apr-Sept Streamflow Volume Yakima River Basin Sea Level Pressure Geopotential Height Surface Air Temperature Relative Humidity Rogue River Basin Surface Pressure Geopotential Height Surface Air Temperature Precipitable Water

Independent Component Analysis l l ICA is a method that separates mutually independent signals from observation data The ICA model: X = A * S X = observed data A = some unknown mixing matrix S = independent components l Assumptions: Ø Observed signals are linear mixtures of unobserved independent signals Ø Independent signals are non-Gaussian and mutually independent

Independent Component Analysis What is ICA? “Independent component analysis (ICA) is a method for finding underlying factors or components from multivariate (multi-dimensional) statistical data. What distinguishes ICA from other methods is that it looks for components that are both Statistically Independent and Non. Gaussian. ” The simple “Cocktail Party” Problem Mixing matrix A x 1 s 1 Observations Sources x 2 s 2 n sources, m observations x 1 = a 11 s 1 + a 12 s 2 x 2 = a 21 s 1 + a 22 s 2 x = As

Independent versus Uncorrelated Two variables of y 1 and y 2 are independent if and only if: P(y 1, y 2) = P 1(y 1)P 2(y 2) l Uncorrelatedness is a weaker form of independence. l When two variables are independent they are also uncorrelated however the inverse is not true.

The independent components must be non-Gaussian for ICA to be possible. Maximizing the non-Gaussianity of WTx results in the independent components.

Measures of Non-Gaussianity • We need to have a quantitative measure of non-Gaussianity for ICA Estimation. • Kurtosis : gauss = 0 (sensitive to outliers) • Entropy : Gauss = largest • Neg-entropy : Gauss = 0 (difficult to estimate) Approximations where v is a standard Gaussian random variable and :

ICA vs PCA üPrinciple Component Analysis (PCA) translates a set of possibly correlated variables into a set of values of uncorrelated variables which is called principal components. üIndependent Component Analysis (ICA) goes beyond and finds independent variables. Therefore some preprocessing steps would make the problem of ICA estimation simpler: § Centering § Whitening

ICA Preprocessing l l l Centering: is achieved by subtracting the mean vector (m) from the mixed signal “x”. After estimating the mixing matrix with centered data, one can add the mean vector of s (obtained from A-1) back to the centered estimates of s. Whitening: to linearly transform the observed vector “x” so that its components are uncorrelated and their variances equal unity.

Independent Component Analysis Mixing Matrix [A] Mixed Signal [X] Original Signals [So]

Independent Component Analysis ICA Mixed Signal [X] Estimated Signals [S]

Independent Component Analysis Estimated Signals [S] Original Signals [So]

Seasonal Climate Signal Correlation with Rogue Apr-Sept Streamflow Volume NCEP/NCAR Reanalysis

Seasonal Climate Signal Correlation with Sprague Apr-Sept Streamflow Volume NCEP/NCAR Reanalysis

Development of Climate Predictors using ICA l ICA climate predictor identification steps: 1. Select one of the large-scale climate data sets to analyze. • The climate data should consist of monthly values and will have dimensions n x m where n is the number of historical observations and m is the number of months. • The number of historical observations, n, used in this study was 30 (1979 -2008).

Development of Climate Predictors using ICA

Development of Climate Predictors using ICA 2. Perform the following steps i times, where i = 1 to m: • Use PCA to find i principal components for the selected climate data set. • Whiten the data using the eigenvectors and eigenvalues from PCA above. • Use ICA algorithm to estimate i independent components. In this study it is assumed that the number of signals is equal to the number of principal components. • Retain each of the estimated i independent components. The number of components retained, r, is given by r = m(m+1)/2. 3. Calculate the Pearson’s correlation coefficient between each r independent climate signal and spring runoff. Retain the climate signal that has the highest correlation coefficient.

Development of Climate Predictors using ICA 4. Repeat steps 2 through 3 k times, where k = 50 for this study. • Repeating steps 2 through 3 is necessary because the ICA program implements a stochastic algorithm and the ICA decomposition is only unique up to sign, scaling and permutation. • This bootstrapping technique is a way of ordering and selecting only those signals that are linearly correlated with spring runoff. The result of step 5 will be a matrix of dimensions n x k. 5. Perform PCA on the matrix n x k from step 4 above, keeping only the first principal component. 6. Whiten the data using the eigenvectors and eigenvalues from step 5. 7. Perform ICA using the whitened data from step 6 to estimate one independent component. The result from this step is a n x 1 column vector which becomes the climate predictor variable in the multiple linear regression model.

ICA climate predictor selection for the Yakima River Basin • • New predictors have strong correlation with spring runoff ICA procedure found 9 new predictors to include in regression model

ICA climate predictor selection for the Rogue River Basin • • New predictors have strong correlation with spring runoff ICA procedure found 9 new predictors to include in regression model

ICA climate predictor selection for the Sprague River Basin • • New predictors have strong correlation with spring runoff ICA procedure found 7 new predictors to include in regression model

Sprague River Basin Data type and Sources Data Type Site No Site ID Agency Elev. Lat Long Period of Record SNOTEL FISH LAKE SWE 478 21 B 04 S NRCS 3371' 47º 32' -121º 4' 1984 -2008 SNOTEL OLALLIE MEADOWS SWE 672 21 B 55 S NRCS 3700' 47º 22' -121º 26' 1984 -2008 SNOTEL SASSE RIDGE SWE 734 21 B 51 S NRCS 4200' 47º 23' -121º 3' 1984 -2008 SNOTEL STAMPEDE PASS SWE 788 21 B 10 S NRCS 3860' 47º 16' -121º 20' 1984 -2008 SNOW TUNNEL AVENUE SWE N/A 21 B 08 NRCS N/A 47º 26' -121º 31' 1984 -2008 SNOTEL FISH LAKE PRCP 478 21 B 04 S NRCS 3371' 47º 32' -121º 4' 1984 -2008 SNOTEL OLALLIE MEADOWS PRCP 672 21 B 55 S NRCS 3700' 47º 22' -121º 26' 1984 -2008 SNOTEL SASSE RIDGE PRCP 734 21 B 51 S NRCS 4200' 47º 23' -121º 3' 1984 -2008 SNOTEL STAMPEDE PASS PRCP 788 21 B 10 S NRCS 3860' 47º 16' -121º 20' 1984 -2008 YAKIMA RIVER AT CLE ELUM, WA STRM 12479500 N/A USGS 1902' 47º 11' -120º 56' 1984 -2008 ICA GLAAM SIGNAL ICA NA NA NOAA NA 90 N 90 S -180 W 180 E 1979 -2008 ICA QBO 30 SIGNAL ICA NA NA NOAA NA 60 N 60 S 90 W -73 E 1979 -2008 ICA NCEP SST SIGNAL ICA NA NA NOAA NA 60 N 30 S 200 W 240 E 1979 -2008 ICA NCEP SLP SIGNAL ICA NA NA NOAA NA 60 N 30 S 140 W 240 E 1979 -2008 ICA NCEP AIR TEMP SIGNAL ICA NA NA NOAA NA 60 N 30 S 200 W 240 E 1979 -2008 ICA NCEP 1000 MB GEO SIGNAL ICA NA NA NOAA NA 60 N 30 S 140 W 240 E 1979 -2008 ICA NCEP PRCP WATER SIGNAL ICA NA NA NOAA NA 5 N 5 S 150 W 160 E 1979 -2008 Data Description

Sprague River Basin Predictor Selection MONTHLY PREDICTORS USED IN ICA CLIMATE MODEL PREDICTORS AVAILABLE Sept Oct Nov Dec Jan Feb Mar Apr SNOTEL TAYLOR BUTTE SWE SNOTEL SILVER CREEK SWE SNOTEL SUMMER RIM SWE SNOTEL QUARTZ MOUNTAIN SWE SNOTEL STRAWBERRY SWE SNOTEL TAYLOR BUTTE PRCP SNOTEL SILVER CREEK PRCP SNOTEL SUMMER RIM PRCP SNOTEL QUARTZ MOUNTAIN PRCP SNOTEL STRAWBERRY PRCP SPRAGUE RIVER ANTECEDENT STRMFLW ICA EPNP SIGNAL ICA GLAAM SIGNAL ICA AMO SIGNAL ICA TNI SIGNAL ICA NCEP 1000 MB ZWIND SIGNAL ICA NCEP PRCP WATER SIGNAL ICA NCEP PRCP RATE SIGNAL

Predictor Identification & Forecast Model Development 53

Volumetric Streamflow Forecast For Sep. 2008 -Apr. 2009 Rogue River Basin 54

Volumetric Streamflow Forecast For Sep. 2008 -Apr. 2009 Yakima River Basin

Model Evaluation and Performance measures Ø The model skill (coefficient of determination and root mean square error) Ø A leave-one-out (LOO) cross-validation procedure Ø The forecast skill of each cross-validation model is evaluated using the Linear Error in Probability Space (LEPS) score: L = 3*(1 -|Pf – P 0| + Pf 2 – Pf + P 02 – P 0) – 1 Ø Benchmark Efficiency:

Yakima River Basin Error Statistics FORECAST ISSUE DATE Sept 1 st Oct 1 st Nov 1 st Dec 1 st Jan 1 st Feb 1 st Mar 1 st Apr 1 st 1 1 2 1 1 1 Calibration RMSE 122. 34 108. 85 88. 32 85. 75 74. 51 75. 4 75. 97 64. 96 Cross-Validation RMSE 137. 93 122. 32 94. 79 95. 16 82. 07 81. 25 84. 14 73. 21 Calibration R 2 0. 78 0. 85 0. 86 0. 9 0. 89 0. 92 Jackknife Cross-Validation R 2 0. 64 0. 72 0. 83 0. 87 0. 88 0. 87 0. 9 --- --- 0. 84 0. 8 0. 79 0. 72 0. 59 0. 62 0. 69 0. 72 0. 77 0. 78 0. 73 0. 76 No. of Comps Benchmark Efficiency Linear Error in Prob. Space

Model Forecast Skill while using Climate Predictors Rogue River Basin BE BE Yakima River Basin Forecast Issue Date LEPS Forecast Issue Date

Thanks! Question/s?