Incorporating Pacific Ocean climate information to enhance the tree-ring-based streamflow reconstruction skill

Streams and rivers are crucial sources of surface water. Rising population, changing climate, and increasing water demands may adversely affect the sustainability of these water resources (Dawadi & Ahmad 2013). Unpredictability of the changing climate may lead to exceedingly high (floods) and low (droughts) streamflow conditions that can cause significant environmental impact (Choubin et al. 2014; Zhang et al. 2016; Amiri & Mesgari 2018; Saifullah et al. 2019). Better planning with regard to water resources can help mitigate the severity of these anthropogenic influences; however, failing to do so may result in adverse socio-economic impacts to society (Ahmad & Simonovic 2000, 2001; Chen et al. 2017). A study conducted by Smith & Katz (2013) for the United States estimated that $210.1 billion and $85.1 billion were expended for drought and non-tropical flood events, respectively, for the years 1980–2011.

Informed and sustainable strategies for water management are dependent upon the knowledge of the past hydrologic variability and the factors that influence it (Zhang et al. 2014; Amiri & Mesgari 2016, 2019; Amiri et al. 2017; Nguyen & Galelli 2018; Martínez-Sifuentes et al. 2020). However, the typical instrumental record lengths do not extend beyond the last century (Tamaddun et al. 2016). Hence, streamflow reconstruction becomes a vital tool for understanding the long-term regional hydrologic variability (Woodhouse 2001; Carrier et al. 2016). Tree-ring chronologies (TRCs) are typically utilized for hydrologic (streamflow, snowpack, and precipitation) reconstructions. This is because trees contain climatic information within their rings, specific to the environment in which the tree has grown (Coulthard et al. 2020). Moisture and temperature receptive tree-rings bring valuable information to the reconstruction technique about past climates, by lengthening the hydrologic data to several centuries, constrained only by the length of the TRC. State-of-the-art streamflow reconstructions using TRC were performed by Stockton and Jacoby in 1976 for the Upper Colorado River (Stockton & Jacoby 1976). The reconstruction skill has since been improved upon because of updated methodologies and the addition of new and revised tree-ring datasets (Meko et al. 2001; Woodhouse & Lukas 2006; Anderson et al. 2019; Chen et al. 2019).

Despite the many advantages of using TRC for hydrologic reconstructions, there are also certain limitations. Growth of the tree is influenced by climatic, biological, and environmental factors (Schweingruber 2012). Climatic factors that influence the tree-ring formation include precipitation, availability and duration of solar radiation, temperature variations, and relative humidity (Wise & Dannenberg 2019). Thus, a lagged response of the tree to such climatic factors due to some biological adaptations may result into a partial or negligible effect of these factors on the tree-ring formation. Furthermore, tree-ring growth is affected by localized factors including site elevation, slope effects, and the presence of other nearby trees (Paulsen et al. 2000). The absence of favorable conditions may disrupt the internal processes of cell division and expansion, and transport and storage allocation (Jarvis & Linder 2000). The presence or absence of certain stresses may lead to negligible growth of tree-rings or the development of more than one annual ring. This brings into question the basic assumption of dendrochronology that pertains to the growth rate of one annual tree-ring, thus leading to uncertainty or noise in TRC (Fritts 2012).

Sea surface temperatures (SSTs) and climate indices have been utilized to study and understand hydrologic variability of continental United States (Mccabe & Dettinger 2002; Corringham & Cayan 2019). The SST and atmospheric pressures influence the streamflow volumes of the continental United States, and many climate changes occur due to the teleconnections between the atmosphere and the ocean (Shabbar & Skinner 2004). Some of the primary climate indices of SST that describe the interannual, decadal-to-multidecadal patterns, and changes in atmospheric pressure (Arctic Oscillation, Atlantic Multidecadal Oscillation, El Niño–Southern Oscillation (ENSO), North Atlantic Oscillation, Pacific Decadal Oscillation (PDO), Pacific North American Pattern, and Southern Oscillation Index (SOI)) describe SST for a predetermined region, which may introduce spatial bias. However, this spatial bias can be removed by utilizing SSTs for the entire Pacific or Atlantic Ocean available in the form of reconstructed SST (Evans et al. 2002; Smith & Reynolds 2004). To overcome the limitations of TRC, incorporating the oceanic and atmospheric climate information can potentially lead to the improvement of the tree-ring-based streamflow reconstruction methodology (Bukhary et al. 2014, 2015; Patskoski et al. 2015; Mukhopadhyay et al. 2018). Mukhopadhyay et al. (2018) successfully reconstructed streamflows for 301 stations across the United States by incorporating the indices of PDO and ENSO along with TRC as predictors in a regression model.

The teleconnections between the hydrologic variable and the SST can be explained by different methods such as Barnet and Preisendorfer, combined principal component analysis (PCA), canonical correlation analysis, second-field PCA, and singular value decomposition (SVD). Bretherton et al. (1992) compared these statistical techniques and determined SVD application to be the simplest and easiest overall and for identifying the leading modes of variability. The statistical methodology of SVD was applied by Shabbar & Skinner (2004) to find coupled regions of winter global SST and Palmer drought severity index (PDSI) for Canadian regions. Over 80% of squared covariance was described by the first three SVD-coupled modes between the SST and the PDSI, while the influence of ENSO and PDO, revealed in the second and third SVD modes, showed about 48% of the squared covariance (Shabbar & Skinner 2004). Aziz et al. (2010) applied SVD for evaluating primary climate indices of ENSO, PDO, and SST influencing the Upper Colorado River Basin snowpack in the United States. Pacific Ocean SST was found to be the driving climate variable. The results indicated the square covariance for the first SVD mode to be in the range of 42–78%, while the first, second, and third modes together described a total variance of 82–92%.

Different modeling approaches exist to reconstruct streamflows such as artificial neural networks, support vector machines, and regression models (Ghumman et al. 2018; Modaresi et al. 2018). Wang et al. (2006) analyzed four rivers located in Asia, Europe, and North America, and determined that nonlinearity is apparent for a smaller timescale (e.g. daily streamflow series) compared with a larger timescale such as annual streamflow series, which is linear. Stepwise linear regression (SLR) methodology is a commonly used method for hydrologic reconstruction (Woodhouse 2001; Margolis et al. 2011; Lara et al. 2015; Patskoski et al. 2015). Woodhouse (2001) used the technique of SLR to reconstruct the annual discharge of Middle Boulder Creek in the Colorado Front Range, extending from 1703 to 1987. The study used TRC as a predictor variable and explained a variance of 70% in streamflow. Margolis et al. (2011) employed the SLR methodology to reconstruct the streamflow of upper Santa Fe River using tree-rings and explained a variance of 62% for the reconstruction length of 1592–2007 and a variance of 50% in streamflow for the reconstruction length of 1305–2007. Variance explained by TRC only, for streamflow reconstruction, may be attributed to the region's moisture supply; however, the variance explained by the SST and the climate indices may be attributed to the moisture supply due to exogenous climatic factors (Patskoski et al. 2015; Mukhopadhyay et al. 2018).

Regional rivers are a crucial source of water supply for the southwest. Changing climate may have played a hand in the recurring hydrologic droughts affecting the region (Kalra et al. 2017; Bukhary et al. 2018, 2020), especially the California basin, as shown by a study by California Department of Water Resources (CADWR 2020). Particularly, the Sacramento River Basin (SRB) and the San Joaquin River Basin (JRB) have been affected by the recurring droughts. The available streamflow records are usually constrained to the past hundred years, and hence may not be able to provide adequate knowledge about the length, intensity, and frequency of historical droughts, which may limit the long-term planning of water resources.

To understand the long-term hydrologic variability, streamflow reconstruction can be effectively employed for SRB and JRB. Studies have shown that incorporating SST along with TRC (Anderson et al. 2012a; Bukhary et al. 2014; Patskoski et al. 2015; Mukhopadhyay et al. 2018), and SOI and PDO along with TRC (Graumlich et al. 2003; Anderson et al. 2012b; Bukhary et al. 2015), resulted in improved hydrologic reconstructions for western United States. Studies have determined that the PDO influences the snowmelt (Dettinger & Cayan 1995; Tamaddun et al. 2017) and rainfall (Fierro 2014; Liu et al. 2016) in California, whereas Sagarika et al. (2016) identified the SST of the El Niño region to have a positive relationship with streamflows in California. Fierro (2014) determined SOI to influence California rainfall. A research gap exists in this area to assess the effectiveness of using climate indices and SST for the improvement of the TRC-based reconstruction of streamflow in the region. The contribution of the current study is in addressing this research gap. This study is an extension of the efforts for improving the traditional streamflow reconstruction methodology that only uses TRC as a predictor variable, by incorporating the Pacific Ocean SST, PDO, and SOI as predictor variables, along with TRC, in an SLR model. Two watersheds in northern California, i.e. SRB and JRB, were used in the analysis. Reconstructions were generated for the duration of 1801–1980 with an overlap period of 1933–1980. Moreover, the SVD method was utilized to determine the coupled regions of Pacific Ocean SST tele-connected with the streamflow.

The study area map for the SRB (70,474 km²) and the JRB (41,129 km²), including the locations of the TRCs and the streamflow stations, is shown in Figure 1. The largest river of SRB and California is Sacramento River located in the Lower Sacramento Basin (Domagalski et al. 1998), whereas the largest river of JRB is San Joaquin River. The lower SRB and the San Joaquin Basin are situated in the semi-arid central Valley of California. The northern and eastern portions of SRB and southern portions of JRB are more arid. The Sacramento River originates from the Modoc Plateau, Klamath Mountains, and Mount Shasta, whereas the headwater of San Joaquin River is located in the Sierra Nevada snow mountains. Sacramento River's largest tributary is the Pit River. Other major tributaries of Sacramento River include the American, Feather, and Yuba Rivers, where some of the gauges used in this study are also located. Major tributaries of JRB are the Tuolumne, Stanislaus, and Merced Rivers. Other main tributaries include the Chowchilla, Cosumnes, Fresno, and Mokelumne Rivers. Both SRB and JRB are characterized by hot/dry summers (mean daytime temperature ≈ 90°F) and cool/damp winters (precipitation occurs between mid-autumn to mid-spring). Mean annual temperatures for the basins have risen by 2°F for the past century (Van Lienden et al. 2014). Sacramento and San Joaquin rivers help irrigate the arid lands of California. Both rivers combine at the Sacramento-San Joaquin Delta before ultimately falling into the Pacific Ocean.

This section informs about the data sources of streamflow gauges, PDO, SOI, SST, and TRC data utilized in the current study. Table 1 lists the sources, lengths, and website links for the datasets employed for this study.

A total of 11 streamflow gauges were used for reconstruction, of which six stations were located in the SRB, while the other five gauges were situated in the JRB (Figure 1; Tables 1 and 2). The historic full natural flow datasets for two gauges in the JRB and four gauges in the SRB were downloaded from the website of CADWR. Unimpaired flow series (free from anthropogenic influences) for two gauges in the SRB and three gauges in the JRB, which are also a part of United States Geologic Survey's (USGS) Hydro-Climatic Data Network 2009 (Lins 2012), were downloaded from the website of USGS.

The location of the TRC dataset is shown in Figure 1. Eighty-one TRCs employed in the current work were downloaded from the International Tree-Ring Data Bank (ITRDB) website (Table 1). Downloaded TRC, designated as ‘sitecodeR.crn’ at the ITRDB website, is raw ring widths standardized and further processed for the removal of autocorrelation. This is because high autocorrelation or linear codependency exists among the TRCs. However, linear regression is based on the assumption that the independent variables are not significantly correlated. If autocorrelation is not removed from the TRC predictors, this can lead to computation of statistically inaccurate and unstable regression coefficients, and imprecise rejection of variables (Hidalgo et al. 2000).

Data sources for both PDO and SOI are shown in Table 1. The PDO is a standardized decadal-scale variability indicator of North Pacific Ocean SST (Mantua & Hare 2002) and represents the leading year-round pattern of SST anomalies of that region obtained by the empirical orthogonal function analysis (Mantua et al. 1997). The downloaded PDO dataset is a 530-year reconstruction developed from a drought/flood index based on historical records of eastern China. Studies have shown that this particular PDO dataset along with TRC resulted in improved hydrologic reconstructions for western United States (Graumlich et al. 2003; Anderson et al. 2012b; Bukhary et al. 2015) and thus was selected as a predictor for this study. The SOI represents the large-scale air pressure variations, during El Niño and La Niña phenomena, occurring between the eastern and western tropical Pacific. It is computed by measuring mean monthly differences in sea-level air pressure between the regions of Tahiti and Darwin, Australia. A dataset for the SOI used in the current study, based on the work of Ropelewski & Jones (1987), was comprised of monthly values that were averaged to create yearly values to be used in this study.

SST data were obtained from two sources: Evans et al. (2002) SST dataset and Smith & Reynolds (2004) SST dataset (Table 1). Smith & Reynolds (2004) used an instrument-based comprehensive ocean–atmosphere dataset to reconstruct global monthly SST having a spatial coverage of 88.0°N to 88.0°S latitude and 0.0°E to 358.0°E longitude. Each cell has a resolution of 2° by 2°. This study used 3,331 SST data cells extending from 30.0°S to 68.0°N latitude and 100.0°E to 270.0°W longitude. The dataset used in the current study is the extended reconstruction SSTs version 3b (ERSST.v3b) (Table 1).

Apart from Smith & Reynolds (2004) SST, the current study also utilized the SST dataset provided by Evans et al. (2002), which can be categorized as Evans-Instrument-based SST (Evans-Instrument SST) starting 1856 and Evans-Coral-based SST (Evans-Coral SST) beginning 1800. The global monthly reconstructed SST dataset spans over the latitude of 60.0°S to 65.0°N and the longitude of 110.0°E to 65.0°W. The resolution of each grid cell is 5° by 5°. The downloaded Evans-Instrument SST dataset comprised yearly (April–March) mean anomalies. Evans-Coral SST is based on the 13 yearly mean coral-derived oxygen isotope time series. For both Evans-Instrument SST and Evans-Coral SST, this study used 655 SST data cells extending from the latitude of 30.0°S to 68.0°N and the longitude of 100.0°E to 270.0°W.

The length of the Smith & Reynolds (2004) SST dataset is comparatively shorter than Evans SST data, while the Evans et al. (2002) SST dataset has a lower resolution grid. Hence, because of the higher grid resolution (and hence greater accuracy) compared with Evans et al. (2002) SST, Smith & Reynolds (2004) SST was chosen to find the linkages between Pacific Ocean SST sectors and streamflow stations. SST sectors tele-connected with streamflow gauges were isolated using the SVD method. Evans SST was employed during forecasting only, and the same teleconnections were assumed for the dataset as that of Smith & Reynolds (2004) SST. Therefore, spatial coordinates of Smith and Reynolds SST sectors that were identified by the application of SVD were likewise used for Evans et al. (2002) SST.

This section explains the reconstruction approach utilized in this study that consists of the technique of SVD (adopted from Bretherton et al. 1992, section 4.1) and SLR (adopted from Anderson et al. 2012a, section 4.2). Matlab 2013a was used as the modeling platform.

The singular values are arranged in a way that the first-mode SVD or first singular value is larger than the second-mode SVD or second singular value and so forth. Then, the squared covariance fraction (SCF) was calculated for the eigenvectors, which is squared eigenvectors normalized to unity (Bretherton et al. 1992) and generates squared covariance for each mode expressed as a fraction or percentage. The first SVD mode explains the maximum covariance. The normalized square covariance (NSC) was also computed, which is squared singular values added and then divided by the product of the number of streamflow stations and SST data cells (Bretherton et al. 1992).

SST data were contained in the left singular matrix, and the streamflow data were contained in the right singular matrix. The first columns of both the matrices were projected onto their left into their corresponding standardized anomalies matrix, resulting in their corresponding temporal expansion series (TES). Correlation values were determined for every streamflow station and every SST grid cell by correlating the first TES of gage data with the matrix of the SST and by correlating the first TES of SST with the matrix of the gage data. Correlation values indicated the coupled regions of the SST and the streamflow data. The SST regions that showed correlations with streamflow stations, at confidence level 90% or above, were retained as predictors for the application of regression procedure.

PCA was applied on the observed streamflow station datasets to create a representative regional streamflow for each basin: SRB only, JRB only, and combined SRB + JRB. The method of PCA reduces large datasets while retaining vital information. Principal components having eigenvalues >1.0 were retained. Pearson's correlation method was used for prescreening of TRCs. Pearson's correlation method was used to correlate regional streamflow with 81 TRCs (p < 0.05, at confidence level 90% or above) for a 20-year moving window for the temporally stable relationship (Timilsena & Piechota 2008; Anderson et al. 2012a). Some of the TRC were located close to streamflow stations but were not selected in prescreening. This may occur due to different reasons already outlined in the literature review, as limitations of TRC.4.2.

In this study, the predictors of TRC, PDO, SOI, Smith and Reynolds SST, Evans-Instrument SST, and Evans-Coral SST were utilized in a stepwise linear regression technique which utilizes the technique of forward selection and backward elimination. Stepwise regression analysis does multiple regressions to add (forward selection at cut-off p-value = 0.05) or remove (backward selection/elimination at p-value = 0.1) the predictor variables based on the strong or weak correlation, respectively, in a stepwise fashion. The selected independent predictor variables are used to predict the dependent predicted variable. For this analysis, reconstructions were developed by the application of SLR on the datasets of TRC, PDO, SOI, and SST. TRCs were used together with oceanic–atmospheric predictors to generate (a) streamflow reconstruction using TRC only, (b) streamflow reconstruction using TRC, PDO, and SOI, and (c) streamflow reconstruction using TRC and Pacific Ocean SST (Smith and Reynolds SST, Evans-Inst SST, and Evans-Coral SST). This resulted in five models of SLR for streamflow reconstruction, which are as follows:

• SLR-1: Traditional TRC only

• SLR-2: TRC, PDO and SOI

• SLR-3: TRC and Smith and Reynolds SST

• SLR-4: TRC and Evans-Instrument SST

• SLR-5: TRC and Evans-Coral SST

The five SLR models were applied to develop reconstructions that were representative of the three regions, i.e. SRB, JRB, and Combined SRB + JRB. The SLR model was verified using calibration R². The increased number of potential predictors may result in the over-fitting of the SLR model (Rencher & Pun 1980). In this study, the validation of the SLR model was accomplished by using the leave-one-out cross-validation technique (Michaelsen 1987; Woodhouse & Lukas 2006), Durbin–Watson (D–W) statistics method (Lara et al. 2015; Anderson et al. 2019), Nash–Sutcliffe Efficiency (NSE; Mccuen et al. 2006; Brocca et al. 2013), RMSE-observations standard deviation ratio (RSR; Mair & Fares 2010; Carrier et al. 2013), scatter plots, and standard error (SE) (Woodhouse & Lukas 2006; Anderson et al. 2012a). Leave-one-out cross-validation is employed for assessing the predictive performance of the model for an independent or unknown data values. One data point is removed temporarily from the training data, and that value is predicted by the model by utilizing the remaining data points. The approach is replicated for all the data points in the dataset. Cross-validation R² or predicted R² assists in identifying if the model over-fits, i.e. not generalized to new or unknown datasets. SE of regressions is an indicator of the precision of the model fit to the known data. It is found by the estimation of the average of length (absolute values) between the data values and the regression line. Low SEs indicate good regression skill of the model. The D–W statistic test investigates the autocorrelation in the residuals that may emerge from the regression analysis. D–W statistics always range between 0 and 4 (Durbin & Watson 1951). The presence of autocorrelation signifies that the important predictors are not included. Non-autocorrelation is indicated by the value 2. A value approaching zero shows positive autocorrelation, whereas negative autocorrelation is revealed by a value moving towards the value 4. RSR is a fraction obtained by dividing the root mean square error with the standard deviation of the observed values (Singh et al. 2004), while NSE (range infinity to 1) indicates how well the model predicts (Nash & Sutcliffe 1970). The study by Moriasi et al. (2007) assessed the performance of a model by using RSR and NSE error statistics for a monthly time step and reported that RSR in the range of 0–0.5 indicated a ‘very good’ model and RSR > 0.7 indicated an ‘unsatisfactory’ model; NSE = 1 indicated a 1:1 match for the observed and simulated values, NSE ≥ 0.75 indicated a ‘very good’ model, and NSE ≤ 0.5 indicated ‘unsatisfactory’ (Moriasi et al. 2007). The criteria established by Moriasi et al. (2007) were also applied to the current study to assess the performance of the SLR models.

The following section describes the SVD application results (section ‘Application of SVD’), results of the five regression models using TRC along with oceanic–atmospheric predictors, and the application of the statistical methods to evaluate the regression model performance for SRB, JRB, and combined SRB + JRB (section ‘Application of SLR’).

The SVD approach was applied to the Smith and Reynolds Pacific Ocean annual mean SSTs. Results for the application of SVD are shown in Figure 2. For SRB, first-mode SCF explained a variance of 97%. The second-mode value was 2.5%, and the third-mode value was 0.3%. The first-mode SCF value for the JRB was 92%, whereas the second- and third-mode values were 6.2 and 1.2%, respectively. For combined SRB + JRB, 91.7% of variance was explained by the first-mode SCF. The second- and third-mode SCF values were 3.9 and 3.3%, respectively. The sum of the first three modes of SCF values such as 99.79% for SRB, 99.27% for JRB, and 98.87% for combined JRB + SRB indicates strong linkage between the streamflow stations and Pacific Ocean SST. NSC values are 3.2, 3.7, and 3.8% for SRB, JRB, and combined SRB + JRB, respectively. Four, five, and three significant regions were identified through the application of SVD for SRB, JRB, and the combined SRB + JRB, respectively (Figure 2).

For the current work, the technique of SVD was applied to identify the teleconnections between the Pacific Ocean SST and streamflow stations. As highlighted by Figure 2, the significant regions were found in the northern, equatorial, and central Pacific Ocean. A close inspection of the SVD-SST map shows that some of the significant regions are located in the Niño 3 (bounded by 90.0° to 150.0° west and 5.0° north to 5.0° south), Niño 4 (located in 150.0° west to 160.0° east and 5.0° North to 5.0° South), and Niño 3.4 (120.0° west to 170.0° west and 5.0° north to 5.0° south). SST anomalies in these regions are indicative of ENSO conditions. Pacific SST shows significant linkages with the streamflow of western United States (Sagarika et al. 2015). Studies have also demonstrated the ENSO influences over streamflow in western United States (Mccabe & Dettinger 2002; Corringham & Cayan 2019). Pacific SST anomalies or ENSO-induced changes affect atmospheric circulation and movement of jet streams, impacting weather patterns and precipitation over the entire North American continent (Eichler & Higgins 2006; Bhandari et al. 2018). Jet streams travel in a wave motion consisting of high air pressure crests and low air pressure troughs influencing the hydrologic parameters of precipitation and thus streamflow. Regions such as California that are generally equator-ward of the polar jet stream, a transient southward excursion of the wintertime jet, should typically result in more precipitation due to added trough passages (Soukup et al. 2009). Pacific jet stream is also a contributing factor of the wintertime precipitation in California depending on its strength and if the jet is more north or more south than normal. The winter season contributes to over 80% of the precipitation in California, thus significantly affecting the regional streamflows (Wahl et al. 2019).

For SRB, prescreening of TRC, based on Pearson's correlation method, resulted in the selection of 20 out of 81 TRCs that showed positive and significant correlations (p < 0.05) with observed streamflow. The application of SLR resulted in successful streamflow reconstructions for the six stations in SRB (Table 3 and Figure 3). For all stations in SRB, at least one SLR model resulted in high-quality reconstructions because of the addition of oceanic–atmospheric predictors. Variance was calculated for the years 1929–1980 (Table 3). Models of SLR-3 and SLR-5 explained a variance of 75 and 72% compared with 69% explained by TRC only, for station AMF. For the FTO station, Evans-Coral SST explained a variance of 77% compared with 74% explained by TRC only. For the SBB station, SOI was dropped and PDO was retained as a predictor, which improved R₂ from 0.75 (when using TRC only) to 0.78 (when using TRC along with PDO). For the YRS station, using the traditional TRC predictor described a variance of 72%. Incorporating Smith and Reynolds SST improved the variance to 76%. Similarly, for USGS station 11383500, inclusion of Smith and Reynolds SST improved the variance to 73% compared with when using TRC (variance = 70%) by itself. For station 11381500, variance was improved by incorporating Evans-Coral SST (variance = 74%) and Smith and Reynolds SST (variance = 77%), compared with using TRC (variance = 71%) only. The oceanic–atmospheric predictors, not selected, can be attributed to the reason that they contained climate information that was already captured by TRC. D–W statistic values for all the stations lie in the range of 1.5–2.22, most of them close to 2, indicating the presence of small or negligible autocorrelation. The best D–W statistic results are shown by the YRS station where the value ranges between 2 and 2.1. The worst D–W statistic results lie in the range of 1.5–1.72 for the SBB station. For all the stations, difference between predicted R² and R² is low, indicating that the model did not overfit values. The low values of SE indicated that the model had a good regression skill. NSE and RSR calculations (Table 3) indicated that the performance of SLR models for SRB ranged between good and very good. Scatter plots with a 45° bisector line (Figure 4) displayed a good positive correlation between observed and reconstructed streamflow.

Streamflow reconstructions for five stations in the JRB were accomplished successfully (Figure 5 and Table 4). For each station analyzed in the JRB, inclusion of at least one SST predictor showed improved results compared with using the base case of TRC only (Table 4). SOI and PDO were not retained as predictors for any of the stations in the JRB, because they did not provide any additional hydrologic variance, already described by TRC. For the CSN station, Evans-Instrument SST along with TRC showed a variance of 66% which was an improvement compared with when using TRC only (variance 62%). Inclusion of Evans-Coral SST demonstrated R² of 0.91, which was an improvement over the TRC-only model (R² = 0.86) for the SJF station. For USGS gage # 11274500, inclusion of Smith and Reynolds SST resulted in slight improvement of R² = 0.85 compared with the base case of TRC only (R² = 0.84). Incorporating Evans-Coral SST resulted in an improved variance of 85% for USGS gage # 11266500, and 60% for USGS gage # 11315000, over the TRC-only model (variance 74 and 51%, respectively). The oceanic–atmospheric variables not retained as predictors in the regression analysis did not incorporate any additional information in the SLR model not already contained in the TRC. D–W statistic results were in the range of 1.82–2.51, indicating slight positive and negative autocorrelation or non-autocorrelation. Results show low SE values. Estimates of RSR and NSE helped analyze the performance of the SLR models (Table 4). Good positive correlation between observed and reconstructed streamflow was demonstrated by the scatter plots (Figure 6). USGS gage # 11315000 and CSN station showed satisfactory model performance, whereas USGS gage # 11266500, 11274500, and SJF station showed very good model performance based on the criteria provided by Moriasi et al. (2007).

For the combined case of SRB + JRB, SLR models for six of the 11 stations did not retain oceanic–atmospheric variables as predictors for the reconstruction of streamflow. The six stations were AMF, YRS, CSN, SJF, and USGS gage # 11274500 and 11266500. The other five stations showed improved results due to the incorporation of SST. For the FTO station, inclusion of both Smith and Reynolds SST and Evans-Coral SST improved the variance to 77%, whereas the variance explained by the TRC predictor was 75%. For the SBB station, the value of R² was 0.69 for TRC. Inclusion of Evans-Coral SST and Smith and Reynolds SST increased the R² value to 0.72 and 0.79, respectively. Smith and Reynolds SST also improved the variance to 75 and 76% for USGS gage # 11383500 and 11381500, respectively, compared with the TRC-explained variance of 72%. Results indicate small SE values. D–W statistics ranged between 1.32 and 2.55. SLR Models 1, 2, and 4 for station SBB showed the lowest D–W statistics of 1.32. D–W statistics for USGS gage # 11315000 ranged between 2.33 and 2.55. The best D–W statistic of 1.96 was displayed by the results of CSN station.

The results showed that inclusion of the climate variables of PDO, SOI, and SST together with TRC enhanced the methodology of streamflow reconstruction compared with when employing TRC by itself for SRB, JRB, and combined SRB + JRB. Overall, incorporating Smith and Reynolds SST (SLR-3) and Evans-Coral SST (SLR-5) as predictors demonstrated better results compared with Evans-Instrument SST, SOI, and PDO. The current study is yet another illustration of the connections between SST, atmosphere, and streamflow volume (Mccabe & Dettinger 2002; Anderson et al. 2012a; Patskoski et al. 2015). The application of SLR resulted in successful reconstructions for SRB, JRB, and combined SRB + JRB. The results are consistent with the existing work conducted in this region (Meko et al. 2001; Meko & Woodhouse 2005). Meko et al. (2001) performed tree-ring-based reconstructions for four streamflow gauges in the SRB. The study explained the variance of 64, 69, and 81% in streamflow for the yearly durations of 869–1099, 1100–1629, and 1630–1905, respectively. Based on the results of the current study, compared with when only using TRC, incorporating information related to climate signals displayed an enhanced reconstruction skill for the regions of SRB and JRB as indicated in Tables 3 and 4, Figures 3 and 5.

Similar results were generated by other studies in western United States (Graumlich et al. 2003; Anderson et al. 2012a). Graumlich et al. (2003) performed Yellowstone River reconstruction (in western United States) by using the traditional reconstruction predictor of TRC together with SOI and PDO in a simple regression model. Results showed that using the predictor variable of TRC only in a regression model explained a variance of 42% in streamflow. Incorporating the Pacific Ocean climate indices of PDO and SOI as predictors along with TRC improved the variance to 59%. Using the climate indices limited the length of Yellowstone River reconstruction to the twentieth century because of the shorter length of the SOI and PDO datasets compared with the longer-spanned TRC-based reconstructions (1706–1977). However, incorporating climate indices improved the skill of the reconstruction methodology (Graumlich et al. 2003). Similar results were displayed for the SBB station in the SRB using the SLR model in the current study. Although the SOI predictor was not retained, inclusion of the PDO predictor along with TRC improved the variance to 78% from 75% when using TRC only. None of the other stations in SRB, JRB, or combined SRB + JRB retained SOI or PDO as predictors. Snowpack reconstructions for the Upper Green River Basin in western United States were performed by Anderson et al. (2012a); the traditional reconstruction methodology using TRC only was improved by the addition of SOI, PDO, and SST. The base case of using TRC only explained a variance of 34% in snowpack. The addition of SOI together with TRC explained an improved variance of 57%. Incorporating SST data (Smith and Reynolds SST and Evans SST) further improved the reconstruction skill by explaining variance in the range of 60–65% (Anderson et al. 2012a). These studies show that climate variables of PDO, SOI, and SST influence the volume of streamflow in western United States. The current study exhibited similar results when using SST as a predictor variable in combination with TRC, and improved reconstructions were displayed for both SRB and JRB. Incorporating SST explained the variance in the range of 73–77% for SRB compared with 69–75% when using TRC only; and for JRB variance using SST as a predictor was in the range of 60–91% compared with 51–86% when using TRC only.

The occurrence of drought is typical to arid and semi-arid regions. Among various indicators, streamflow is one of the most commonly used drought-monitoring indicators. Hence, results generated by the current study may also be indicative of the strong links between oceanic indices and drought (Kingston et al. 2013). California has a history of recurring droughts. According to CADWR (2020), based on precipitation records of California State, the period spanning between October 2012 and October 2014 was marked as the driest consecutive 3-year period. Observed instrumental records offer a relatively small window in the long-term hydro-climatic variability of the region. However, using paleoclimate records for streamflow reconstruction may provide useful insights into the near long-term climatic variability of the region as shown by the current study (CADWR 2020).

In the current study, streamflow reconstruction was performed for the regions of SRB, JRB, and combined SRB + JRB by utilizing an SLR model and incorporating the climate indices of PDO and SOI and the Pacific Ocean SST (Smith and Reynolds SST and Evans SST), together with the traditional TRC as predictors. The application of the methodology resulted in successful reconstructions that are indicative of the long-term hydro-climatic variability of the region. The improved reconstruction methodology has expanded the streamflow records back to 1800. The results using the five SLR models, each using a different predictor variable, showed that reconstruction methodology is improved upon by incorporating the oceanic–atmospheric predictors compared with when using the conventional predictor of TRC by itself. It was concluded that the traditional predictor of TRC may not be a true indicator of climatic influences. Incorporating oceanic–atmospheric predictors of SST, PDO, and SOI together with TRC improved the reconstruction skill of the model. The application of the SVD method isolated the coupled Pacific Ocean SST sectors and thus improved the streamflow reconstruction methodology. Incorporating Smith and Reynolds SST (SLR-3) and Evans-Coral SST (SLR-5) as predictors along with TRC demonstrated better results compared with Evans-Instrument SST (SLR-4), SOI and PDO (SLR-2) along with TRC, and TRC only (SLR-1).

Although the study satisfactorily provided streamflow reconstructions using climate information, there were certain limitations associated with the current work. Due to the availability of climate data of shorter time span compared with TRC, the length of the reconstructions was curtailed to 1800, while the TRC may extend back to several centuries. Still the shorter reconstructions were helpful in depicting a more accurate picture of near long-term regional hydro-climatic variability. Overall, this study helped in greater understanding of the past regional hydrologic variability. Improvements in reconstruction methodology may help water managers to develop efficient and informed strategies for future water availability. Future work may focus on improving the reconstruction methodology by adopting several stochastic or data-driven type modeling approaches. Furthermore, the methodology can be easily transferred to other watersheds.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

All relevant data are included in the paper or its Supplementary Information.