Interannual hydroclimatic variability and its influence on winter nutrient loadings over the Southeast United States

It is well established in the hydroclimatic literature that the interannual variability in seasonal streamflow could be partially explained using climatic precursors such as tropical sea surface temperature (SST) conditions. Similarly, it is widely known that streamflow is the most important predictor in estimating nutrient loadings and the associated concentration. The intent of this study is to bridge these two findings so that nutrient loadings could be predicted using season-ahead climate forecasts forced with forecasted SSTs. By selecting 18 relatively undeveloped basins in the Southeast US (SEUS), we relate winter (January-February-March, JFM) precipitation forecasts that influence the JFM streamflow over the basin to develop winter forecasts of nutrient loadings. For this purpose, we consider two different types of low-dimensional statistical models to predict 3-month ahead nutrient loadings based on retrospective climate forecasts. Split sample validation of the predictive models shows that 18–45 % of interannual variability in observed winter nutrient loadings could be predicted even before the beginning of the season for at least 8 stations. Stations that have very high coefficient of determination ( > 0.8) in predicting the observed water quality network (WQN) loadings during JFM exhibit significant skill in predicting seasonal total nitrogen (TN) loadings using climate forecasts. Incorporating antecedent flow conditions (December flow) as an additional predictor did not increase the explained variance in these stations, but substantially reduced the root-mean-square error (RMSE) in the predicted loadings. Relating the dominant mode of winter nutrient loadings over 18 stations clearly illustrates the association with El Ni ño Southern Oscillation (ENSO) conditions. Potential utility of these season-ahead nutrient predictions in developing proactive and adaptive nutrient management strategies is also discussed.


Introduction
Concerted efforts to improve national water quality conditions resulted in the enactment of the 1972 Clean Water Act with section (303) d requiring the states and territories to list impaired water bodies and develop total maximum daily loads (TMDLs) for these waters.Despite these efforts and frequent updates to TMDLs, the US Environmental Protection Agency's recent update reveals that nutrients affect 20 % of impaired and 12 % of the assessed river miles (EPA, http://www.epa.gov/owow/tmdl/,2006).The increase in aquatic nutrients might result from population growth as well as from increased fertilizer application (Meybeck, 1982;Vitousek et al., 1997).However, natural variability associated with weather (e.g., hurricanes) and climatic events (e.g., El Niño) could also induce significant increase in nutrient concentrations beyond critical levels (Chen et al., 2007) even if the basin is not experiencing any pressure from urban development or changes in agricultural practice.Thus, for these virgin basins, we have the opportunity to estimate the seasonal nutrient loadings due to the potential changes in runoff during the season.
The National Research Council (NRC, 2002(NRC, , 2001) ) has emphasized that a detailed understanding of various sources of uncertainties, including the role of climate change and climate variability, is required for improving water quality prediction in natural systems.One of the dominant and wellunderstood modes of global climatic variability is El Niño Southern Oscillation (ENSO) that has a periodicity of 3-7 yr and exhibits anomalous warm/cold sea surface temperature (SST) conditions in the equatorial Pacific, thereby modulating the climate particularly in the tropics and sub-tropics (Ropelewski and Halpert, 1987).Considerable research now exists on the recurrence and regime structure of ENSO and its teleconnections to rainfall/streamflow, and their potential predictability of interannual hydroclimatic variability over the United States (Ropelewski and Halpert, 1987;Dettinger and Diaz, 2000;Devineni and Sankarasubramanian, 2010).It is also well known that instream nutrient concentration and loadings primarily depend on streamflow variability (Borsuk et al., 2004;Paerl et al., 2006;Lin et al., 2007) and antecedent flow conditions (Vecchia, 2003;Alexander and Smith, 2006).Recent studies on the relationship between coastal water quality conditions and SST conditions also show that there is a strong association between climatic modes and concentrations of phosphorous (Childers et al., 2006), aquatic vegetation (Cho and Poirrier, 2005), and chlorophyll and phytoplankton levels (Arhonditsis et al., 2004).However, systematic research in associating climatic variability to instream nutrient variability and utilizing that linkage to estimate season-ahead nutrient loadings is very limited.
Most of the studies on estimating instream nutrient concentrations have focused primarily on predicting the average annual concentrations using runoff and various basin attributes (Smith et al., 1998(Smith et al., , 2003;;Mueller and Spahr, 2006;Mueller et al., 1997).Studies have also recommended approaches to predict daily and seasonal loadings and concentration nutrients using streamflow and their time of observation (Cohn et al., 1992;Runkel et al., 2004).However, these nutrient models rely on the observed information (e.g., streamflow) during that season, which has limited utility in developing season-ahead estimates of nutrients.Findings from the hydroclimatic literature clearly show that interannual variability in streamflow can be predicted by developing low-dimensional models contingent on SST conditions (Devineni et al., 2008) as well as using precipitation forecasts from general circulation models (GCMs) (Sankarasubramanian et al., 2008).Similarly, water quality literature emphasizes that streamflow is the most important descriptor in explaining nutrient variability (Cohn et al., 1992;Runkel et al., 2004;Cohn, 2005).To our knowledge, this is the first effort that associates the interannual variability in all of the above noted three variables -climate, streamflow and total nitrogen (TN) -to develop TN forecasts over a region.The purpose is to understand the controls that are required for developing skillful seasonal nutrient forecasts and also to assess how the skill in hydroclimatic predictions translates into skill in nutrient forecasts over the regional scale.For this purpose, we consider two low-dimensional models that consider season-ahead precipitation forecasts and streamflow conditions as predictors for developing season-ahead estimates of winter nitrogen loadings over the SEUS.
The manuscript is organized as follows: a brief description of precipitation forecasts, streamflow and water quality databases employed in the study is first provided in Sect. 2. Following that, Sect. 3 provides the details of the low-dimensional statistical models and skill measures utilized in developing and evaluating the season-ahead nitrogen loadings forecasts.In Sect.4, we present results from the winter nutrient forecasts developed using the lowdimensional models.Next, we discuss the potential implications of the findings in the context of developing adaptive water quality management plans.Finally, in Sect.6, we summarize the findings and conclusions from the study.

Data sources
In this section, we discuss various hydroclimatic and water quality databases employed for associating climate forecasts with the nutrient loadings over the Southeast US.We consider 18 watersheds (Fig. 1) for understanding the role of hydroclimate in influencing interannual nutrient variability.Previous studies have shown that winter precipitation and streamflow over the Southeast US are heavily influenced by the ENSO variations (Ropelewski and Halpert, 1987;Devineni and Sankarasubramanian, 2010).The selected 18 watersheds span over seven states and the streamflow with drainage area ranging from 136 km 2 to 44 547 km 2 (Table 1).Given that the selected watersheds are minimally impacted by anthropogenic influences, we hypothesize that the interannual variability in winter nutrients could be explained by the precipitation variability as well as by the antecedent flow conditions.For this purpose, we assemble hydroclimatic and water quality databases for developing season-ahead nutrient forecasts over these 18 watersheds.

HCDN streamflow database
Given that the intent of the study is to associate interannual variability in winter nutrient loadings to climatic variability, we focus our analysis on 18 undeveloped basins over the Southeast United States (SEUS) from the Hydro-Climatic Data Network (HCDN) database (Slack et al., 1993).Daily streamflow records in the HCDN basins are purported to be relatively free of anthropogenic influences such as upstream storage and groundwater pumping, and the accuracy ratings of these records are at least "good" according to United States Geological Survey (USGS) standards.The HCDN database contains the mean daily discharge for about 1600 sites across the continental United States with an average length of 48 yr. Figure 1 shows the location of 18 HCDN stations, and Table 1 provides the list of the 18 stations considered in this study along with their drainage areas.Since the streamflow data (Q) in the HCDN database are available only up to 1988, we have extended it up to 2009 based on the USGS historical daily streamflow database.

WQN Water Quality Network database
USGS provides national and regional descriptions of stream water quality conditions in Water Quality monitoring Network (WQN) across the nation (Alexander et al., 1998)   monitoring networks from both large watersheds (National Stream Quality Accounting Network, NASQAN) and minimally developed watersheds (Hydrologic Benchmark Network, HBN).We employ the observed daily concentrations of total nitrogen (TN) available for these 18 stations from the NASQAN.Observed streamflow during the time of sampling is also available as part of the WQN database.The available water quality data vary from 10-30 yr depending on the measured water quality variable and station.By ensuring the selected watersheds are from HCDN basins, we basically ensure that both the streamflow and water quality data are minimally affected by anthropogenic influences.For additional details about WQN, see Alexander et al. (1998).The selected 18 HCDN stations have observed TN concentrations for 12-22 yr (Table 1).However, the number of samples for each station ranges from 54-152 daily observations with an average of 5-7 observations per year.

Simulated nutrients database
Though nutrient data in the WQN database are available for 12-23 yr over 18 watersheds (see Table 2), their samplings are intermittent.Using the daily observation over this period, we first obtain continuous daily nutrients for the observed period using the LOADing ESTimation (LOADEST) program developed by USGS (Runkel et al., 2004).LOADEST is a statistical model that estimates daily loadings based on the observed daily streamflow and the centered time (dtime) of the year of the observation (Runkel et al., 2004).
ln(L j ) = a0 + a1 ln(Q j ) + a2 ln Q 2 j + a3 sin(2π dtime) +a4 cos(2π dtime) + εj ... (1) where L j denotes the observed daily loadings from the WQN database with j denoting the day of observation; Q j is the observed daily flow and dtime is the centered time which is a function of the observation's number of days (from 1 January) in the calendar year; a0-a4 denote the model coefficients, and εj is the estimated residual for the model.The expression dtime is centered to avoid multi-collinearity, and dtime also represents the seasonality in loadings pattern.For a detailed expression on dtime, see Cohn et al. (1992).
The LOADEST model allows the user to select the bestfitting regression model from 11 predefined regression models using the Akaike information criterion (AIC) (Akaike, 1974).Five regression models that include a linear time trend are not appropriate, since we are employing observed streamflow to estimate simulated loadings beyond the observed period.Therefore, the simulated nutrient loadings based on the remaining regression models (i.e., model forms: 1, 2, 4 and 6 as defined in Runkel et al., 2004) in the LOADEST program do not have any time trend.Equation (1) represents the model form 6. Model form 1 (Eq. 1) considers only the first two (three) terms on the right-hand side (RHS) of Eq. ( 1), whereas model form 3 considers all the terms except the third term in the RHS of Eq. ( 1).For further details on model forms, see Runkel et al. (2004).Table 2 shows the goodness of fit statistics (coefficient of determination (R 2 ) and AIC) in predicting the observed daily loadings in the WQN database (Table 1) and the coefficients of the best fitting regression model for TN for the selected 18 stations.From Table 2, we infer that R 2 ranges from 0.83-0.97indicating good fit of the observed daily loadings over 18 stations.
In this study, we primarily focus on developing nutrient forecasts for the winter season.Given that streamflow peaks in winter over the Southeast US, it has been shown that loadings are at their peak during the same season (Muller and Spahr, 2006).From the perspective of winter forecasts too, this is a season that has significant skill in predicting observed precipitation by the ECHAM4.5 model (Devineni and Sankarasubramanian, 2010).Given that the selected 18 basins are virgin watersheds whose flows are minimally impacted by anthropogenic influences, the primary source of total nitrogen is from non-point that includes fall foliage and post-harvest agricultural lands (Muller and Spahr, 2006).Thus, increased flow during the winter season primarily carries the TN loadings from fall foliage and agricultural lands.Given the focus on the JFM season, we also report the ability of the LOADEST model in predicting the observed JFM nutrients in Table 2. From R 2 (LOADEST) and RMSE during JFM, we clearly see that the performance of the LOADEST model in predicting JFM nutrients is poor in stations 5, 6, 13, 16 and 18.This potentially will have impact in developing nutrient forecasting model for the sites with low R 2 in JFM.
To relate retrospective climate forecasts (discussed in the next section), we estimated the daily TN loadings from 1957-2009 using the observed streamflow data available from the extended HCDN database and the best fitting regression model (in Table 2).The simulated daily loadings obtained from the LOADEST model over the period 1957-2009 are aggregated during JFM to develop winter loadings (L t ) of TN where t represents the year.Given that the estimated loadings are based on adjusted maximum likelihood estimation (AMLE) procedure in the LOADEST model, the simulated daily and the aggregated winter loadings are statistically unbiased (Cohn, 2005).We also computed R 2 (LOADEST) and the root-mean-square error (RMSE (LOADEST) ) for the simulated winter TN loadings obtained from the LOADEST model.For additional details on computing errors in seasonal predictions, see Cohn (2005).Further, to ensure that there is no trend in the winter loadings, we performed Mann-Kendall test.At 1 % significance level, null-hypothesis with Kendall's tau being not equal to zero was rejected in all of the sites for TN.We also performed regional Mann-Kendall test to account for spatial correlation among the 18 stations (Douglas et al., 2000).The p-value for TN is 4 % indicating no trend at the regional level.Our study will consider the simulated winter TN loadings (L t ) available during 1957-2009 for relating the interannual hydroclimatic variability to nutrient variability over 18 stations in the SEUS.

Climate forecasts database
Seasonal climate forecasts are typically developed either using atmospheric GCMs (AGCMs) or using coupled GCMs (CGCMs).In the case of former, it is a two-tiered system, in which SSTs are forecasted first using a statistical/dynamical model and then they are forced with AGCMs.In CGCMs, since ocean and atmospheric models are coupled, they are run in a continuous mode.Recent studies clearly show that AGCMs are more skillful than CGCMs (Goddard et al., 2003).Further, Devineni and Sankarasubramanian (2010) show that ECHAM4.5 precipitation forecasts explain 25-36 % of the variability in observed precipitation over the Southeast US.For this study, we utilize the retrospective winter precipitation forecasts from ECHAM4.5 general circulation model forced with constructed analogue SSTs (http://iridl.ldeo.columbia.edu/SOURCES/.IRI/.FD/ .ECHAM4p5/.Forecast/ca sst/.ensemble24/.MONTHLY/ .prec/,International Research Institute of Climate and Society (IRI) data library) (Li and Goddard, 2005).Retrospective precipitation forecasts from ECHAM4.5 are available for 5 months in advance for every month beginning January 1957.To force the ECHAM4.5 with SST forecasts, retrospective monthly SST forecasts were developed based on the observed SST conditions in that month based on the constructed analogue approach.For additional details on forcing ECHAM4.5 using constructed analogue SST forecasts, see Li and Goddard (2005).
Figure 1 also shows the locations of 56 grid points of precipitation forecasts from ECHAM 4.5 along with their latitude and longitude over SEUS.For this study, we utilize only the forecasted mean (which is obtained by computing the average of 24 ensembles) of winter retrospective precipitation forecasts issued in the beginning of January for developing 3month ahead retrospective nutrient forecasts over the period 1957-2007.

Low-dimensional models development and performance validation metrics
Given that winter streamflow over the SEUS is predominantly rainfall driven with limited snow accumulation, we hypothesize that precipitation is the primary driver in controlling the JFM loadings.To verify this, we correlate simulated JFM loadings with both observed precipitation (Fig. 2a) and principal components of the forecasted precipitation from ECHAM4.5 (Table 3).Principal components basically signify the reduced dimensions of the gridded precipitation forecasts.Given that the precipitation forecasts from the selected grid points are spatially correlated, it is better to obtain low/reduced dimensions of the forecasts using principal component analysis.In this study, we only employ Spearman rank correlation for performing all correlation analyses.Similarly, the computed rank correlation was checked for statistical significance (i.e., 1.96/(n − 3) 0.5 at 95 % confidence interval, where n denotes the number of data points used in calculating the correlation).Thus, the computed correlation in Fig. 2a needs to be greater than 0.29 (n = 50) to indicate statistically significant relationship between the observed precipitation and simulated loadings.From Fig. 2a, we infer that the correlation between observed precipitation and simulated loadings is statistically significant and greater than 0.55 for all the basins.Given this dependency, we first identify relevant grid points (Table 3) of JFM precipitation forecasts that have statistically significant correlation with JFM-observed precipitation for each watershed.Nearest grid points that are significantly correlated to each watershed (Fig. 1) are selected.The variance explained by the first principal component (PC1) of the precipitation forecasts from these grid points is around 74-95 % indicating the strong spatial correlation among the gridded forecasts.Further, from Table 3, we also infer that rank correlations between PC1 of precipitation forecasts with streamflow, and seasonal loadings of TN, are statistically significant for all stations (> 0.29) with the only exception being station 18.The primary reason for such low correlation for station 18 is due to the poor coefficient of determination from the LOADEST model (Table 2) in predicting observed WQN data.Given that the PC1 of precipitation forecasts explain almost the same amount of observed variability in precipitation, streamflow and TN (Table 3 and Fig. 2a), it is logical to develop TN forecasts conditioned on precipitation forecasts.This is expected since the selected basins are virgin watersheds with non-point loadings being the primary source of TN.Given that PCs of precipitation forecasts are significantly correlated with the observed precipitation, streamflow and TN loadings, predictive models of TN loadings can be a powerful tool to obtain season-ahead nutrient forecasts.

Source of climatic information influencing the winter TN variability
To understand the source of climate information that modulates the TN variability over the SEUS, we performed principal component analysis on the simulated loadings (L t ) of TN over 18 stations.The first component approximately explains 59 % of total variability in TN loadings over 18 stations.It is well known in the hydroclimatic literature that ENSO is one of the important climatic conditions that influence the winter precipitation, temperature and streamflow over the SEUS (Ropelewski and Halpert, 1987).Figure 2b shows the correlation between the first component of JFM TN loadings over 18 stations and JFM Nino3.4 -an index used to denote ENSO conditions by averaging the SSTs (Kaplan et al., 1998) over the tropical Pacific (170 • W-120 • W; 5 • S-5 • N).From Fig. 2b, we infer that roughly 36 % of the variability in the first principal component of nutrient loadings over SEUS could be explained purely based on ENSO conditions.ENSO plays an important role on the winter climate of the US since its peak activity typically coincides during December-February.In fact, the precipitation forecasts from ECHAM4.5 incorporate the forecasts of tropical SST conditions (i.e., Nino3.4 region), which are obtained from constructed analogue SST forecasts for forcing the ECHAM4.5.Thus, ENSO is one of the sources of climatic variability that primarily influence both JFM hydroclimatic and nutrient variability over the SEUS.Based on the information provided in Fig. 2 and Table 3, we understand that there is scope for using the low-dimensional components of precipitation forecasts for developing season-ahead forecasts of TN loadings over 18 selected stations.

Low-dimensional models
Given that our interest is primarily in understanding how large-scale hydroclimatic information could be utilized for seasonal nutrient predictions over the SEUS, we consider two low-dimensional models: principal component regression (PCR) and canonical correlation analysis (CCA).Lowdimensional models reduce the correlated predictors and predictands so that a subspace of uncorrelated predictors and predictands could be used for regression model development (Tippet et al., 2003;Sankarasubramanian et al., 2008).Further, these low-dimensional models also recalibrate the GCM forecasts so that any bias in predicting the JFM long-term mean of the nutrients is removed based on the regression model.Brief description of the low-dimensional models is provided next.

Principal component regression (PCR):
PCR, which is otherwise known as model output statistics (MOS) (Wilks, 1995), eliminates systematic errors and biases in GCM fields and also recalibrates the principal components (PCs) of GCM fields to predict the hydroclimatic variable of interest using regression analyses.The predictand could be streamflow (Q t ) or loadings (L t ) over a watershed.Since the gridded precipitation forecasts over a given region are spatially correlated, employing precipitation forecasts from multiple grid points as predictors would raise multicollinearity issues in developing the regression.To avoid this, we employ PCR based on Eq. ( 2): where L t denotes the estimate of daily average TN loadings during the JFM season in year t, PC k t denotes the k-th PCs from the retained K PCs of precipitation forecasts and βs denote the regression coefficients, whose estimates are obtained by minimizing the sum of squares of error.We employ stepwise regression to select K PCs out of the rotated grid points of precipitation (given in Table 3) for developing the PCR model.From Table 4, we infer that most of the stations (except stations 8, 10 and 11) require only up to the first four principal components for developing the PCR model.

Canonical correlation analysis (CCA):
In PCR, we develop separate regression models for each site.Given that the predictands, the winter loadings, across the basins are also spatially correlated, one could utilize that information to develop a reduced set of regression models.CCA is a low-dimensional regression modeling framework that utilizes inter-site correlations to relate many (multiple predictands)-to-many (multiple predictors).Consider winter loadings available from m sites represented by L T = (L 1 , L 2 ,..., L m ) (dimension: nXm), whose corresponding p grid points of precipitation forecasts (p > m) are represented as X T = (X 1 , X 2 ,...,X p ) (dimension: nXp) where "T" denotes the transpose.Then, canonical correlation analysis finds a linear combination of the p predictors, L * = b T L, that maximally correlates with the linear combination of m predictands (X * = a T X).Mathematically, the canonical correlation is obtained by choosing the vectors a and b that maximize the relationship where denotes the variance-covariance matrix between the two matrices in the subscript.For a detailed mathematical treatment of CCA, see Wilks (1995).Number of components from m predictands and p predictors to be retained for the regression is decided based on step-wise regression.Squared values of canonical correlation represent the percentage of variance explained in each predictand by the predictors under that dimension.Thus, CCA allows us to develop a reduced set of models that can be used to predict loadings for each site based on the precipitation forecasts.
Before performing CCA, we first group the basins based on k-means clustering (Hartigan and Wong, 1979) so that CCA could be performed on each cluster.Based on clustering, four groups were identified (Table 5) with the sites having the highest loadings placed under group 1 and the lowest average loadings placed under group 4. Separate CCA was performed for each group.For instance, CCA on group 1 is performed on loadings from two sites (m = 2) and the corresponding grid points of precipitation forecasts for the two sites (# 13 and #18) from Table 3 are combined (p = 20) as predictors.The skill in predicting the winter loadings for each station is evaluated based on two different skill scores, which are discussed next.

Skill scores for nutrients forecasts
To evaluate the skill in predicting the interannual variability in winter TN loadings using climate forecasts, we consider two error metrics -coefficient of determination (R 2 ) and root-mean-square error (RMSE) per unit area of the watershed (A).These metrics need to account for both sources of errors: error in predicting the observed JFM nutrient loadings from the WQN database using the LOADEST program (R 2 (LOADEST) , RMSE (LOADEST) ) (see Table 2) and the error in predicting simulated JFM nutrients from LODEST based on the low-dimensional model (R 2 (PCR/CCA) , RMSE (PCR/CCA) ).Since these two models are developed independently, R 2 and RMSE in predicting winter nutrient loadings using climate information could be expressed as follows: For each station, R 2 (PCR/CCA) and RMSE (PCR/CCA) were computed based on the estimated TN loadings from the lowdimensional models and the simulated winter TN loadings for the period 1957-2006.Thus, we compute skill measures using Eqs.( 3) and ( 4) to quantify our ability to explain the interannual variability in winter TN loadings using precipitation forecasts from ECHAM4.5.

Results and analyses
To ensure that the skill in forecasting winter nutrients is reliable, we evaluate the low-dimensional models based on two different types of validation: leave-X out cross-validation (LCV) and split-sample validation (SSV).Both these methods are commonly adapted in forecasting literature for validating the model (Wilks, 1995).

TN loadings forecasts based on PCR models
For validating the PCR models under LCV, the methodology suggested by Towler et al. (2009) is modified to evaluate the skill of the model over 51 yr .The LCV steps for PCR models are described as follows: (i) 10 % of the data (5 yr) are randomly removed along with the year for which the prediction is desired; (ii) a PCR model is developed using the remaining 45 yr of loadings (L t ) and retained PCs; (iii) the developed model is then used to predict the left-out year, and (iv) steps (i) to (iii) are repeated to develop prediction for each year, and skill measures (R 2 and RMSE) were computed based on the 51 yr of predicted data.This entire procedure (i)-(iv) is repeated 100 times and a box plot of R 2 (Fig. 3) and the median of RMSE (in Table 4) are presented.
Figure 3 shows the box plot of R 2 under LCV for 18 stations.Under LCV, we compute R 2 based on the predicted loadings for 51 yr.Hence, R 2 needs to be higher than 0.08 (correlation > 0.29) to demonstrate statistically significant skill in predicting season-ahead nutrient loadings.However, more than 12 stations exhibit R 2 greater than 0.16 over 100 trials of LCV.From Fig. 3, 16 stations show statistically significant skill for TN.The developed PCR model under LCV explains more than 10 % of interannual variability in TN loadings in all the 100 different fittings (Fig. 3) except stations 6 and 18.For the rest of the 16 sites, the correlation between the predicted nutrient loadings obtained using climate forecasts and the loadings simulated from LOADEST using   the WQN database is greater than 0.29, which is statistically significant for the 51 yr of data considered.The forecasted TN loadings show statistically insignificant relationship with observed TN loadings that the correlation coefficients are 0.28 and 0.16 for stations 6 and 18 respectively.Poor goodness of fit (see Table 2, R 2 (LOADEST) ) from the LOADEST model is the primary reason behind the poor performance of these two stations during the winter season.Further, station 18 shows poor correlation between the principal components of precipitation forecasts and JFM loadings (Table 3).Another possible reason for such poor prediction by LOADEST model in those two stations is the limited number of years of data availability (see Table 1) with station 6 (18) WQN observations spanning 14 (12) years having a total of 56 (57) daily samples.The median RMSE (Table 4) computed under LCV also shows that error in predicting the observed WQN loadings during the winter season is lesser than 1 kg day −1 km −2 for most of the stations.
Under split-sample validation (SSV), PCR models are developed using L t and PCs available over the calibration period (PCR: 1957(PCR: -1986) and skill measures are computed in predicting L t during the validation period .Hence, R 2 needs to be higher than 0.21 (correlation > 0.46 for 21 yr of data) to demonstrate statistically significant skill in predicting season-ahead nutrient loadings.Based on this, Fig. 4 indicates that 11 stations (2-4, 7-11 and 13-15) show significant skill in predicting TN loadings.Stations 6 and 18 perform poorly because of the limited number of years of WQN data which results in very low R 2 of the LOADEST model.Apart from these two stations, stations 1, 5, 12, 16 and 17 also show insignificant skill (R 2 < 0.21).Stations 5 and 16 perform poorly due to the poor skill of the LOAD-EST model during JFM (see Table 2, R 2 (LOADEST) ).Similarly, stations 1, 12 and 17 perform well under LCV, but the skill is statistically insignificant under SSV, even though the simulated TN (Fig. 2) loadings exhibit significant correlation with both observed precipitation (Table 2) and forecasted  6, 8, 9, 10, 11, 15, 17 895.95precipitation (Table 3).It is important to discuss the importance from the perspective of developing real-time nutrient forecasts.LCV basically exploits all the available data in the future years to predict the TN loadings for the left-out year.Thus, LCV shows that there is potential in developing seasonal nutrient forecasts for sites 1, 12 and 17.However, under SSV, we cannot guarantee that skill in developing real time (without using future information) forecasts, since the trained model using the data from 1957-1986, is not capable of developing a statistically significant forecast for the period 1987-2007.But, as we collect more data in the future, we may be able to develop statistically significant forecasts for these stations.Thus, LCV shows the potential skill in developing the forecast, whereas SSV shows the demonstrable skill in developing real-time forecasts.Thus, based on two different validation methods, we understand that 11 stations (2-4, 7-11 and 13-15) exhibit statistically significant skill in predicting the observed WQN loadings using the PCR model developed separately for each site.Next, we evaluate the ability to predict the loadings in these stations under a different low-dimensional model -canonical correlation analyses -that utilizes the spatial correlation in the TN loadings to develop a predictive model.

TN loadings forecasts based on canonical correlation analyses
Four different CCAs were performed on each group (listed in Table 5), and the developed models were evaluated under LCV and SSV.For LCV, we simply perform leave-5 out cross-validation instead of repeated fitting of the model (as described in Sect.4.1 for the PCR model).Under leave-5 out cross-validation, we randomly leave out five predictands and predictors along with the year for which the prediction is desired for each station under a given group (in Table 5) and CCA model is developed using the rest of the 46 yr of data.The developed CCA model was employed to predict the year for which the prediction is desired.This procedure was repeated for all the years of observation under a given group to develop the CCA model estimated loadings.The R 2 CCA and RMSE CCA were estimated between the simulated winter loadings from the LOADEST model and predicted loadings from the CCA model for each station under all the four groups.R 2 CCA and RMSE CCA were further adjusted  according to Eqs. (3) (Fig. 5) and Eq. ( 4) (Table 4) to account for the errors in the LOADEST model in predicting the WQN database.From Under SSV (Fig. 6 and Table 4), we compute R 2 based on the predicted loadings during 1987-2007 using the CCA model developed over the period 1957-1986.From Fig. 6, CCA model did not exhibit any skill in predicting the winter loadings in stations 5, 6, 11-13, 16-18.Comparing this with the PCR model performance, the CCA model performs similarly with the exception being very low R 2 at station 11.For station 11, PCR (R 2 = 0.47) performs significantly better than the CCA model (R 2 = 0.14).One possible reason for such poor performance of CCA model is that station 11 has low correlation with the rest of the sites under group #4 which has 8 stations.Under SSV, R 2 of the CCA model for the rest of the stations is almost similar to that of R 2 of the PCR model.However, the RMSE CCA is consistently higher than the RMSE PCR (Table 4).This implies that the conditional bias (overprediction and underprediction) of the CCA model is much higher.One possible reason for such increased conditional bias under CCA is due to increased heteroscedasticity in the observed loadings under a given group.However,   the ability of the CCA model to explain the observed variance in loadings is almost comparable to that of the PCR model, indicating that the source of interannual variability in winter nutrients is the same across the region.To summarize, using ECHAM4.5 precipitation forecasts alone, we infer both low-dimensional models demonstrate significant ability in predicting the observed winter TN loadings in nine coastal stations (#2-4, 7-10 and 14-15) based on two different validation methods.

Discussion
Analyses presented in Figs.2-6 show that interannual variability in nutrient loadings could be predicted well before the beginning of the season contingent on the climate forecasts.
Though instream loadings primarily depend on streamflow and precipitation variability during the season, antecedent moisture/flow conditions also play a critical role in influencing the nutrient loadings from the watershed (Vecchia, 2003;Alexander and Smith, 2006).At seasonal time scales, antecedent flow conditions could be considered as the surrogate for basin storage or initial conditions in influencing the streamflow variability.To understand the role of antecedent storage conditions, we consider the observed December streamflow at each station as an additional predictor along with the gridded precipitation forecasts (Table 3) to develop nutrient forecasts for each station.Forecasts of TN loadings were developed using both PCR model (Fig. 7a) and CCA model (Fig. 7b), and the modified R 2 and RMSE (Table 6) are computed (Eqs.3 and 4) based on SSV by evaluating the model over the period 1987-2007. Comparing . Comparing    statistically significant skill by adding streamflow as an additional predictor.However, by using streamflow as an additional predictor, R 2 of the CCA model substantially improved over all stations, which indicates the importance of incorporating basin information in spatial-dimension reduction.Further, RMSE of both PCR and CCA models is substantially reduced by adding the observed December streamflow as an additional predictor.This implies that antecedent storage/flow conditions are very critical in reducing the conditional bias in developing season-ahead TN forecasts resulting in reduced over/underprediction compared to the models developed using the precipitation forecasts alone.Thus, from a process control perspective, given the good skill in the reconstructed seasonal nutrient loadings, the interannual variability in nutrient loadings could be partially explained based on climatic variability.But, to obtain improved prediction (i.e., RMSE), it is important to incorporate both climatic variability and antecedent storage conditions in developing season-ahead nutrient forecasts.

Role of basin characteristics
We also compared the R 2 (Fig. 7a) and RMSE (Table 6) to basin characteristics such as drainage area (figure not shown).But, we did not find any significant relationship between the skills of the nutrient forecasts to basin characteristics.It is difficult to associate the spatial variability in the skill with the basin characteristics just using 18 stations.Further, at interannual time scales, we understand that climatic signals are the primary source of variability (Fig. 2b) followed by basin storage conditions (Fig. 7).Performing these analyses over the continental scale may provide additional information on the role of basin characteristics in influencing the skill in predicting nutrients.

Validating nutrient forecasts with observed precipitation
To understand how well the forecasts predict the extreme seasons of loadings over 9 stations (#2-4, 7-10 and 14-15), we compare the predicted loadings with the observed wet and dry precipitation seasons.For this purpose, we computed the 33rd and 67th percentiles of observed JFM precipitation and JFM-simulated loadings based on the available record over the period 1957-2007.Based on this, if the observed precipitation in a given year falls below (above) 33rd (67th) percentile of precipitation, then that particular year is denoted as below-normal (above-normal).If the observation falls between 33rd and 67th climatological percentiles, then it is denoted as normal year.Similarly, based on the 33rd and 67th climatological percentiles of TN loadings, we also grouped the forecasted nutrients into three categories: below-normal, normal and above-normal.Based on these two tercile categories available for precipitation and forecasted nutrients, we computed the false alarm rate for the forecasted nutrients estimated using forecasted precipitation and streamflow as predictors (Fig. 7a) based on the PCR model over the period 1987-2007.The false alarm ratio gives the probability of predicting above-normal (below-normal) nutrient year given the observed precipitation for that year falls under below-normal (above-normal).If this ratio is very high, then it indicates the forecasting model issues false alarms too often contrary to the observed category.Based on Table 7, the probability of issuing false-alarm is around 4-5 % for all the 9 stations.It is important to note that this false alarm rate also incorporates the errors that would have resulted from the forecasted precipitation.Results from Table 7 show that the forecasted nutrients modulate well with the observed extreme seasons in precipitation over the validation period 1987-2007 and could be utilized for issuing categorical nutrient forecasts.

Potential for application -water quality trading
Perhaps the most important utility of the season-ahead forecasts of nutrient loadings is in promoting water quality trading.Some of the successful water quality trading programs in the country (e.g., Tar-Pamlico River Basin and Neuse River Basin in NC) typically allow trading nutrient loadings across different point sources as well as with non-point sources (e.g., farmers participating in the voluntary nutrient reduction program) through the basin-level trading association so that the seasonal/annual load caps are always met from the basin.
Research on climate forecasts and water allocation clearly shows that probabilistic streamflow forecasts could be effectively utilized to specify the failure probability of reservoir  releases as well as in ensuring the end of season target storage conditions being met with high probability (Sankarasubramanian et al., 2009).Similarly, in the context of seasonal water quality management, the developed forecasts of loadings could be used to estimate the probability of violation of target loadings for the upcoming season (Towler et al., 2010;Borsuk et al., 2002).One could also develop an optimal nutrient loading model such that the probability of violating the total loadings from multiple sources is within the acceptable level.Thus, utilizing season-ahead forecasts of nutrient loadings and updating them throughout the season provide an opportunity to develop adaptive nutrient control strategies that ensure target nutrient loadings and desired concentration.

Alternate forecasting models and nutrient forecasts for urbanized watersheds
Thus, the intent of this study is to understand how well climate and basin storage conditions control the development of skillful forecasts of TN loadings and to evaluate the performance of two low-dimensional models in issuing seasonahead nutrient forecasts utilizing climate forecasts and basin storage conditions.In principle, the analyses provided here could also be extended with other sophisticated statistical models including nonparametric and Bayesian hierarchical models to estimate the entire conditional distribution of loadings.Similarly, one could also develop nutrient forecasts by forcing the mechanistic water quality model with forecasted streamflow and water temperature, which in turn could be obtained based on dynamical downscaling (Leung et al., 1999) or statistical downscaling (Devineni et al., 2008) based on the climate forecasts.
Given that the developed models are statistical models, it requires observed nutrients and streamflow information for application in other basins.However, CCA model could be employed for prediction in ungauged basins within the cluster since it employs both spatial and temporal dimension reduction.However, to apply this approach for urban watersheds with significant anthropogenic influences, the developed model could be used to predict the upstream loadings for the undeveloped segment of the watershed.However, to estimate the loadings leaving the fully developed watersheds, one needs to know the loadings from the point sources so that net loadings leaving the watershed could be estimated in total.We are currently working on extending this approach with weather forecasts, so that point and non-point sources could be controlled in both virgin basins as well as in basins whose loadings are significantly impacted by anthropogenic influences.

Summary and conclusions
The study primarily focused on understanding the process controls in estimating winter nutrient loadings by considering 18 HCDN watersheds over the SEUS.Given the discontinuous observed daily TN loadings, the study reconstructed simulated TN loadings using the LOADEST model for the winter season.The ability to predict these simulated loadings was validated with two low-dimensional models that utilize winter precipitation forecasts and pre-season flow conditions.
Out of 18 stations, a total of 9 stations (#2-4, 7-10 and 14-15) exhibited statistically significant skill in predicting the observed winter nutrient loadings under both lowdimensional models based on two different validation methods.However, the reported skill in predicting the TN loadings accounts for both error from the LOADEST model as well as the error from the low-dimensional models.Findings from the study could be summarized as the following "controls" that influence the skill in predicting seasonal TN loadings: stations that have very high R 2 (LOADEST) (> 0.8) in predicting the observed WQN loadings during the winter (Table 2) exhibit significant skill in loadings.This high R 2 from the LOADEST model could be considered as a potential criterion for developing nutrient forecasts if the basin's hydroclimatology exhibits significant association with climatic signals.Incorporating antecedent flow conditions (December flow) as an additional predictor did not increase the explained variance in these stations, but substantially reduced the RMSE in the predicted loadings.Understanding the source of climatic variability that controls the TN variability revealed that Nino3.4,an index denoting ENSO conditions over the tropical Pacific, accounted for 36 % of the observed spatial variability in the TN loadings over the SEUS.Given that using climate forecasts have been very beneficial in improving reservoir management over seasonal time scale (Sankarasubramanian et al., 2009), we argue the need to develop nutrient loading forecasts conditioned on climate forecasts.Our future work will utilize these seasonal nutrient forecasts in developing adaptive water management plans over the SEUS.

Figure 1 :
Figure 1: Location of the 18 Hydro-Climatic Data Network (HCDN) stations along with the considered grid points of precipitation forecasts over the Southeast United States (SEUS).

Fig. 1 .
Fig. 1.Location of the 18 Hydro-Climatic Data Network (HCDN) stations along with the considered grid points of precipitation forecasts over the Southeast United States (SEUS).

2Figure 2 :Fig. 2 .
Figure 2: Rank correlation between the simulated total nitrogen (TN) loadings from the LOADEST model and observed precipitation (a) and the association between El Nino Southern Oscillation (ENSO) index -Nino3.4-and the first principal component (PC1) of TN loadings over the selected 18 stations. 3

Figure 3 :
Figure 3: Box-plot of R 2 (based on equation (3)) of principal component regression (PCR) model predicted TN loadings obtained using PC's of forecasted precipitation under leave-one-out cross validation (LCV).

Fig. 3 .
Fig. 3. Box plot of R 2 (based on Eq. 3) of principal component regression (PCR) model predicted TN loadings obtained using PCs of forecasted precipitation under leave-one-out cross validation (LCV).

Figure 4 :
Figure 4: Modified R 2 (based on equation (3)) of PCR model predicted TN loadings obtained using PC's of forecasted precipitation under split sample validation (SSV).

Fig. 4 .
Fig. 4. Modified R 2 (based on Eq. 3) of PCR model predicted TN loadings obtained using PCs of forecasted precipitation under split sample validation (SSV).
Fig. 5, stations 5, 6 and 18 do not exhibit statistically significant correlation in predicting the loadings from the WQN loadings.As discussed under PCR model, stations 5, 6 and 18 did not perform well because of the limited number of years of WQN data and low R 2 of the LOADEST model during the winter season.The rest of the sites exhibited statistically significant relationships in explaining the observed variability in the WQN loadings.For stations 8 and 9, CCA model explains 30-40 % (correlation 0.55-0.63) of the observed winter variability in TN within the WQN database.With regard to RMSE, the performance of CCA model and PCR model is almost similar at most of the stations having an error less than 1 kg day −1 km −2 . 5

Figure 5 :
Figure 5: Modified R 2 (based on equation (3)) of canonical correlation analysis (CCA) model predicted TN loadings obtained using PC's of forecasted precipitation under LCV.

Fig. 5 .
Fig. 5. Modified R 2 (based on Eq. 3) of canonical correlation analysis (CCA) model predicted TN loadings obtained using PCs of forecasted precipitation under LCV.
Fig. 7a and b with Fig. 4 (PCR model) and 6 (CCA model), we infer that only station #13 (under CCA model) has resulted in

Figure 6 :
Figure 6: Modified R 2 (based on equation (3)) of CCA model predicted TN loadings obtained using PC's of forecasted precipitation under SSV.

Fig. 6 .
Fig. 6.Modified R 2 (based on Eq. 3) of CCA model predicted TN loadings obtained using PCs of forecasted precipitation under SSV.

Figure 7 :Fig. 7 .
Figure 7: R 2 (based on equation (3)) of PCR(a) and CCA(b) model predicted TN loadings obtained using both PC's of forecasted precipitation and December streamflow under SSV.

Table 1 .
Baseline information for 18 selected stations showing the number of years of observed daily records of TN available in the Water-Quality Monitoring Network (WQN) database.Values in the parentheses under number of years column show the total number of daily observations available for each station.

Table 2 .
Performance of LOADEST model in predicting the observed TN loadings from the WQN database.Models with linear time components (Model No.: 3, 5, 7-9) are not considered.The skill in predicting JFM WQN loadings is separately given in the last two columns.

Table 3 .
Rank correlation between observed winter streamflow, TN loadings with the first principal component of the winter precipitation forecasts for the 18 selected stations.Locations of grid points indicated in the Table are shown in Fig.1.

Table 4 .
Skill, expressed as root-mean-square error (RMSE) (based on Eq. 4), in predicting winter TN loadings using climate forecasts.Tablealsogives the number of principal components considered and the percentage variance explained by them for the total grid points selected (given in Table3) for each station.

Table 5 .
Grouping of 18 selected stations based on k-means clustering.

Table 6 .
RMSE (Eq.4)of forecasted TN loadings based on PCR and CCA models that consider ECHAM4.5 precipitation forecasts and December streamflow as predictors under SSV.

Table 7 .
The false alarm rate of nutrient forecasts developed using precipitation forecasts and December streamflow.