Pacific Decadal Oscillation Climate Variability and Temporal Pattern of Winter Flows in Northwestern North America (original) (raw)

1. Introduction

It is important to explore the linkages between large-scale climate variability and regional-scale hydroclimatic processes because interannual and decadal-scale climate variability is instrumental in planning regional water resources, which are of significant ecological, cultural, and economic value. In some of the previous studies, the variability of winter precipitation and temperature in northwestern North America (NWNA) was associated with the Pacific decadal oscillation (PDO), in addition to other large-scale ocean and atmospheric phenomenon (e.g., Stahl et al. 2006; Woo and Thorn 2008). Also, it has been shown that the Arctic Oscillation and PDO have teleconnections with river flows of the Hudson Bay in Canada (Déry and Wood 2004) and with the flows of southeastern Alaskan rivers (Neal et al. 2002). Quite recently, the influence of PDO on winter and other seasonal streamflows and their temporal pattern of occurrences have been studied in an attempt to explore changing climate–streamflow interactions (e.g., Woo and Thorn 2008; Burn 2008; Brabets and Walvoord 2009). Also, in an earlier study by Chen et al. (2004), it was speculated that the relation between climate variability and groundwater levels in a major aquifer in Manitoba is influenced by PDO variability. Along similar lines, the current study concentrates on the influence of PDO on winter (January–March) flows and changing climate–winter flow interactions in NWNA so as to advance our understanding of such mechanisms. During the January–March period, many of the streams run under iced conditions and hence winter flows are believed to be largely groundwater fed.

The PDO is a pattern of Pacific Ocean sea surface temperature (SST) anomalies, noted about 12 years ago by Mantua et al. (1997) when studying the influence of climate variability on fish production in Alaska. It is believed that the NWNA is a PDO-sensitive region, wherein patterns of climatic fields such as precipitation and temperature and hydrological response (i.e., streamflows) are influenced by this climate variability. The PDO signal is stronger for the winter and therefore it would influence winter and spring flows more compared to other seasonal flows. It has been suggested that the PDO pattern fluctuates between positive (warm) and negative (cold) phases on an interdecadal time scale of about 20–30 yr (Mantua et al. 1997). However, the mechanism by which this climate pattern persists for several years in the same phase has not yet been identified. For the winter period over NWNA, a positive (negative) PDO phase corresponds to periods of high (low) temperature and low (high) precipitation. The relationship of PDO and climatic fields is relatively straightforward compared to its relationship with streamflows because of the influence of basin characteristics such as geology, land forms, vegetation, etc. The former relationship could also get complicated because of local effects such as altitude and distance to large water bodies (Stahl et al. 2006; Woo and Thorn 2008).

It is shown in this paper that the serial structure of PDO indices strongly resembles that of a stochastic process with long-term persistence (LTP). The idea of the LTP mechanism was envisioned more than 50 years ago by Hurst (1951) when he was designing water storage facilities for the Aswan Dam. Since then, it has been shown by numerous researchers that the LTP mechanism is present in many hydroclimatic time series, for example, streamflows (Montanari et al. 1997; Koutsoyiannis 2006), rainfall (Montanari et al. 1996), temperature (e.g., Syroka and Toumi 2001; Gil-Alana 2003; Koutsoyiannis and Montanari 2007) and wind speed (Govindan and Kantz 2004). Despite these developments, the LTP mechanism has not received overwhelming popularity. However, its use has been extended to other scientific areas for stochastic modeling purposes such as network traffic modeling (e.g., Liu et al. 2006), medical sciences (e.g., Leite et al. 2007), among others. As far as the analysis of trends in hydroclimatic time series is concerned, the LTP mechanism considerably affects the statistical significance of trends, for example, see Cohn and Lins (2005), Kallache et al. (2005), Hamed (2008), and Khaliq et al. (2009a) in this regard. Climatic fields (i.e., temperature and precipitation) are the main drivers of the hydrological response and, because of their physical link with the PDO phenomenon (which possesses characteristics of an LTP-like stochastic process), it would be plausible to investigate trends in streamflows with the LTP assumption. However, there is a possibility that the intensity of LTP in streamflows could become modified because of basin characteristics. Because of this reason and the influence of other unidentified factors, the serial structure of streamflows could also resemble that of a stochastic process with short-term persistence (STP) or no persistence at all. On the basis of these arguments, it would be sensible to investigate temporal changes in winter flows with multiple stochastic assumptions rather than assuming merely independence (IND), which is the simplest assumption.

The objectives of this paper are twofold. On one hand, the influence of PDO climate variability on winter flows is investigated, using a much larger database compared to previous studies wherein only smaller subsets of the same database were explored, in an attempt to develop more detailed PDO–streamflow linkages. On the other hand, the influence of various stochastic assumptions (i.e., IND, or STP, or LTP) on the changing patterns of winter flows is investigated. The influence of such assumptions, particularly those of type STP and LTP, on trend estimates in streamflows has not been investigated in earlier studies in NWNA and hence the significance of trends reported may have been overstated. As it is difficult to investigate LTP-like structure, particularly in short samples, that are common in northern regions, a significant portion of this paper is devoted to the investigation of trends in the presence of STP-like serial structure by introducing a suite of semiparametric methods.

The remainder of this paper is organized as follows. Section 2 describes the study area and various databases used in the study followed by a description of the analysis methods in section 3. Section 4 contains results of the study wherein the influence of PDO on winter flow regimes is explored thoroughly to elucidate PDO–streamflow linkages followed by an analysis of the influence of various stochastic assumptions on the changing patterns of winter flows in NWNA. Finally, section 5 provides a discussion and a summary of the results obtained and their implications.

2. Study area and data

NWNA includes Alaska, Yukon and Northwest Territories, British Columbia, and Alberta. Several climate zones exist in this region. Hare and Thomas (1979) subdivided the region climatically into the Pacific, Cordillera, Prairies, Boreal, and Arctic (Fig. 1). The western Cordillera with chains of lofty mountains, plateaus, and valleys dominates the region. East of the Cordillera lies the interior plains and the Canadian Shield. Typically, spring melt alone or spring melt with rain generates high flows (as well as spring floods) that are several orders of magnitude larger than the winter flows. Spring high flows are followed by a decline in flow, which is revived occasionally by summer/fall rainstorms.

Various databases used in this study are described below along with a short description of the studies wherein smaller subsets of the same databases were explored earlier. Woo and Thorn (2008) considered 110 streamflow stations, but temporal trends were estimated only for a smaller subset of these stations that had at least 35 yr of record for the period 1960–99; four of the six stations explored in Neal et al. (2002) were also considered in this study, and nearly all stations represent natural flow regulation. Brabets and Walvoord (2009) analyzed annual, winter, April–September monthly, and fall flow regimes for only 21 stations located in upper, middle, and lower parts of the Yukon River basin for the period 1944–2005. Annual, summer, and winter low flow regimes of 62 pristine river basins included in the Canadian reference hydrometric basin network (Brimley et al. 1999; Harvey et al. 1999) and located in Yukon and Northwest Territories, British Columbia, and Alberta were investigated in Khaliq et al. (2008) within the framework of a Canada-wide study. Though there were some similarities between these studies with respect to the chosen study areas, it is important to mention that only some of the stations were found common in these studies. Therefore, for the present analyses, a much larger dataset of 179 stations, consisting of a majority of the stations explored earlier in separate studies, was adopted to develop more detailed spatial analysis. A list of these stations is provided in Table 1, and detailed characteristics of these stations can be found in the references quoted above. The last year of data was fixed at 2007. Data for the Alaskan stations were obtained online (at http://waterdata.usgs.gov/nwis) and for the Canadian stations (at http://www.wsc.ec.gc.ca/rhc-wsc/). As the focus is on the January–March period during which many of the streams run under iced conditions, it is very common to find missing or zero values in instrumental records. Therefore, for a year (January–March) to be included in the analysis, at most 10 missing values are allowed. This criterion resulted in 163 stations that are used in this paper for studying PDO–winter flow interactions. The effect of missing values and that of leap years is taken into account by studying seasonal average daily flow rather than the accumulation of daily flows for the January–March period. For the analysis of temporal trends, a station was retained in the analysis if at most 5 yr were found missing. This additional criterion resulted in 72 stations with record lengths ≥50 yr.

Monthly values of the PDO indices for the period 1900–2007 were obtained from the University of Washington and the National Oceanic and Atmospheric Administration’s Joint Institute for the Study of the Atmosphere and Ocean (JISAO 2009). This PDO index is based on monthly SST anomalies in the Pacific Ocean north of 20°N. Annual and seasonal mean values of the indices were obtained by simply averaging the corresponding monthly values.

3. Methodology

a. Methods for detection of change points in time series of PDO indices

It has been suggested in the literature that the PDO shifts phases from positive to negative (and vice versa) on an interdecadal time scale. When such shifts had occurred over time is the sole purpose of change-point analysis included in this paper. A change point is a “point in time” that divides a time series into two independent epochs with significantly different mean values. In this paper, a change point is taken as the first year of the new regime (or epoch) in a manner similar to the work of Elsner et al. (2004) and Khaliq et al. (2007). Two parametric and two nonparametric tests are used to identify the most probable time points when a PDO regime shift, from negative (positive) to positive (negative), had occurred during the period 1900–2007. Parametric tests include cumulative departure plots combined with a Student’s t test (Dahman and Hall 1990) and the Worsley likelihood ratio test (Worsley 1988), and nonparametric tests include the Pettitt test (Pettitt 1979) and a repeated rank-sum test (Siegel and Castellan 1988; Lanzante 1996).

b. Methods for correlation analysis and PDO–streamflow interactions

To investigate the PDO–winter flow interactions and cross-correlation among various stations chosen for the study, both parametric [i.e., the Pearson’s product moment correlation (PPMC)] and nonparametric [i.e., the Spearman’s rank order correlation (SROC)] measures of linear association are employed. These two measures of correlation will help address the sensitivity issues related to underlying distributional assumptions. Implementation details of these two methods can be found in Walpole and Myers (1985). In addition to correlation analyses, the PDO–streamflow linkages are investigated using plots of (i) averaged percentage differences and (ii) averaged standardized differences from the period-of-record mean. The results of such plots would help develop generalized interpretations of the influence of positive and negative phases of the PDO on winter flows. Compared to the first set of plots, the second set of plots would help address the influence of second moment of winter flows on differences from the period-of-record mean.

c. Methods for analysis of serial structure of time series of PDO indices and winter flows

In addition to first-order autocorrelation analysis, which is often investigated by assuming a Markovian process of order one, the stochastic structure of PDO indices, for various monthly combinations, and that of the winter flows is analyzed using three different measures of the Hurst exponent (commonly denoted by the symbol H): (i) rescaled adjusted range statistic (RARS) (Mielniczuk and Wojdyłło 2007), (ii) aggregated standard deviation (ASD) (Koutsoyiannis 2003, 2006), and (iii) detrended fluctuation analysis (DFA) (Peng et al. 1994; Kantelhardt et al. 2001). In the literature (e.g., Kurnaz 2004), it has been argued that the DFA method, compared to the other two, is superior because it takes into account the effects of local trends, which could bias the estimates of the Hurst exponent. Together, the results of these three methods would help to overcome such biases. Values of the exponent H are interpreted as follows: the range 0.5 < H < 1 (0 < H < 0.5) corresponds to a persistent (antipersistent) process and the value H = 0.5 corresponds to an independent process, which is valid in an asymptotic sense. To take into account the sample size uncertainties on the estimates of the Hurst exponent, Monte Carlo (MC) simulation–based confidence intervals were developed in a manner as in Khaliq et al. (2008).

According to the findings of a review paper on trends (Khaliq et al. 2009b), there is no universally accepted approach for investigation of trends in time series of hydroclimatic observations, so various methods have been adopted for this purpose, including parametric and nonparametric ones. However, in a large majority of the studies on trends, nonparametric methods have been used, particularly the Mann–Kendall (MK) trend test (Kendall 1975). The same test is retained in this paper as well; a description of this test is provided in the appendix. The MK test was originally designed for applications where the time series data were assumed to be IND and identically distributed. For such time series data, the observation at time point t is not dependent on observations recorded at previous time points (t − 1), (t − 2), … Performance of the MK test for IND times series is well documented (e.g., Helsel and Hirsch 1992), and its use is generally well accepted among practitioners. However, if the assumption of IND cannot be reliably verified and dependence does exist in a time series, then the performance of the MK test is seriously affected (e.g., von Storch 1995; Hamed and Rao 1998; Yue and Wang 2004; Khaliq et al. 2009a); that is, for a positively (negatively) autocorrelated time series, the MK test would suggest a significant trend more (less) often compared to an IND series. To address these issues, semiparametric methods that are based on the joint use of the MK testing procedure and time series modeling and simulation concepts are presented. Such approaches have been lacking in the literature on trends. For the analysis of temporal changes in time series of winter flows, it is assumed that the time series consists of a deterministic trend component At and a pure stochastic component Bt; that is,

i1525-7541-11-4-917-e1

i1525-7541-11-4-917-e1

This formulation implies that the effects of both deterministic and stochastic components can be separated. The deterministic component could be assumed linear; that is, At = a + bt, which had been generally the case for the majority of the studies on trends. However, this restriction is ad hoc and can be relaxed. The stochastic component Bt could have a STP- or LTP-like dependence structure.

When the MK test is applied to an observed time series Yt, the test statistic _S_obs reflects the effects of both deterministic and stochastic components At and Bt, respectively. If one can develop a distribution of the statistic S by simulating a large number of time series of the stochastic component Bt and calculating _S_sim for each of these simulated time series, then the significance of the observed statistic _S_obs can readily be obtained from the simulated distribution of S. If _S_obs lies in the tails of the simulated distribution of S, then the deterministic trend component is more likely to be statistically significant. It should be noted that the distribution of _S_sim reflects the effects of stochastic behavior of data and, hence, the influence of random variations is accounted for.

For this category of dependence, it is often assumed that the autocorrelation structure of observational records conforms to an autoregressive moving average (ARMA) process, that is, ARMA (p, q) in which p is the number of autoregressive terms and q the number of moving average terms, like short-term correlation structure. In the majority of the studies on trends, the ARMA(p, q) structure is simplified further by assuming p = 1 and q = 0. Under this assumption, commonly known as AR(1), the standardized stochastic component Bts = (Bt − μB)/σB in which μB and σB are the mean and standard deviation of the stochastic component, can be written as

i1525-7541-11-4-917-e2

i1525-7541-11-4-917-e2

where ɛ_t_ is an uncorrelated noise process with variance _σ_ɛ2 and r_1 is the lag-1 autocorrelation coefficient. When the observational records conform to a normal marginal distribution, the process ɛ_t is often modeled as having mean 0 and variance _σ_ɛ2. Under nonnormal situations, the type of distribution of the noise process must be taken into account. Since the exact distribution of the noise process is unknown, some convenient procedures need to be adopted to simulate a sample from the stochastic component, which should resemble the observed counterpart in the statistical sense. The following MC simulation–based approaches are considered:

  1. The noise terms are bootstrapped to generate random components, which in turn are used to generate realizations of the stochastic component. This approach is referred to as MCBS in the rest of the paper.
  2. The ɛ_t_ terms are modeled as a polynomial function of the standard normal variable z, as explained below:
    i1525-7541-11-4-917-e3
    i1525-7541-11-4-917-e3
    The parameters υ and u can be estimated by ordinary or weighted least squares, taking the parameter g equal to the value of ɛ_t_ at z = 0. The latter constraint ensures continuity of the relationship. The weighted least squares approach helps to avoid the influence of so-called outliers on the estimated parameters. The degree of polynomials was decided depending upon the closeness of fit as measured by least squares. This approach is referred to as MCPOL.
  3. The stochastic component Bt is log transformed before obtaining the standardized stochastic component Bts. The log transformation is often used to overcome nonnormality. In this case, the distribution of ɛ_t_ would closely resemble the standard normal distribution; this approach is referred to as MCLOG.
  4. The ɛ_t_ terms are assumed to follow generalized extreme value (GEV) and five-parameter Wakeby distributions. The former (latter) approach is referred to as MCGEV (MCWAK).

In addition to these four semiparametric methods, two modified MK (MMK) tests, based on the effective sample size approach, one proposed by Hamed and Rao (1998) and the other by Yue and Wang (2004), are also implemented for comparison purposes. These approaches are referred to as MMK1 and MMK2, respectively. The modified prewhitening (MPW) approach suggested by Wang and Swail (2001) is also used. This approach is based on the joint use of the MK testing procedure and the Sen (1968) nonparametric slope estimator in an iterative mode to address the influence of lag-1 autocorrelation on trend and vice versa. For implementation details of these approaches, the reader is referred to the respective references.

For this category of dependence, the modified MK scaling test developed by Hamed (2008) is used. This method is referred to as MMKS, and a detailed step-by-step procedure can be found in Hamed (2008). The MC simulation and fractional autoregressive integrated moving average (FARIMA) modeling approach, with normal marginal distribution (Cohn and Lins 2005) and nonnormal marginal distribution (Khaliq et al. 2009a), can also be implemented. For the sake of brevity, the results of the latter set of methods are not included in this paper. However, the overall characteristics of their results were similar to those of the MMKS test.

3) Interaction of deterministic and stochastic components

In practice, it is difficult to estimate parameters of the deterministic and stochastic components independently since both interact. A common approach had been to remove the deterministic trend component from observations by using the Sen slope estimator and estimating parameters of the stochastic component from trend-free observations (e.g., Hamed and Rao 1998; Yue and Wang 2004). In this paper, parameters of the deterministic and stochastic components were estimated jointly using nonlinear optimization procedures by maximizing an objective criterion, such as the likelihood function of the assumed stochastic process. Additional details for such optimization-based methods can be found in, for example, Cohn and Lins (2005). In the case of an AR(1) process, this technique can be readily modified and implemented by assuming a linear deterministic trend component and jointly estimating the trend-independent lag-1 autocorrelation and the magnitude of trend in an iterative mode using simple root-finding procedures (Press et al. 1989): (i) assume a magnitude of trend and remove it from all observations; (ii) estimate the lag-1 autocorrelation from trend-independent values from step (i); and (iii) repeat steps (i) and (ii) until the difference between lag-1 autocorrelations of current and the previous iteration converges to a chosen very small value, for example, 0.001.

4. Results

a. Change points in time series of PDO indices

Time series of mean monthly values of the December–March PDO indices for the period 1900–2007 are plotted in Fig. 2a. This time series was subjected to four change-point detection tests, which suggest that 1942 and 1976 are the end points of two different regimes. Thus, this time series can be partitioned into three segments/regimes: 1900–42, 1943–76, and 1977–2007. Because of very limited winter flow records during 1900–42, this study concentrates on the second (i.e., 1943–76) and third (i.e., 1977–2007) segments only. It is important to mention here that a time series of mean monthly values of the October–March PDO indices was used in the work of Woo and Thorn (2008) compared to that of the December–March period used in the present study. For the October–March time series, the results of change-point analysis were the same as reported above for the December–March time series. Thus, it makes little difference using either of these two time series.

The PDO indices for the period 1943–76 are predominantly negative and those for the period 1977–2007 are predominantly positive. According to the published literature, the former period is referred to as the negative (or cold) phase and the latter as the positive (or warm) phase of the PDO. However, these multidecade epochs contain intervals of up to few years in which the PDO indices were reversed, for example, the positive PDO indices for the period 1958–1961 and strong negative PDO indices for the period 1989–91. It is also important to mention here that, at the time of writing this paper, it was not possible to find any credible published evidence concerning entrance of the PDO into the negative phase by year 2007. Hence, given this uncertainty, it is assumed that the positive phase persists until the year 2007. It is shown later in the paper that the time series of PDO indices is strongly persistent. However, the change-point detection tests were applied under the assumption of IND, which was required for these tests. If persistence is exploited fully, then it is possible that the statistical significance of identified change points may not be as strong as it is with the IND assumption. Cumulative departures from the overall mean for the monthly and annual time series of PDO indices are shown in Fig. 2b. This figure demonstrates that the regime shifts are not present only in time series of winter PDO indices, they are also present in monthly and annual time series; this result further strengthens the observation that shifts occurred around 1942 and 1976 over the history of PDO indices. A review of relevant literature led us to conclude that there are more variations for the exact location of the former shift compared to the latter. This point was also noticed when the same four change-point detection tests were applied to monthly and various averaged bimonthly, trimonthly, and so on, time series of PDO indices.

b. Effect of PDO variability on winter flows

First, the cross-correlation of winter flows observed at stations studied is presented followed by analyses of PDO–winter flow correlations and average percentage differences and standardized differences from the period-of-record mean. Of these analyses, the first analysis will show the extent of cross-correlation present among various stations in NWNA. By definition, correlation measures the degree of linear relationship between two variables. A large positive (negative) value indicates positively (negatively) correlated data; absolute magnitude of the correlation coefficient measures the strength of the relationship. Distributions of the PPMC and SROC, considering all possible pairs of stations, are shown in Fig. 3a. Both shapes of the distributions and proportions of significant cross-correlations, for the overall period and for the periods corresponding to positive and negative phases of the PDO, indicate that NWNA is dominated by positive cross correlations. However, the possibility of negative cross-correlations cannot be ruled out. The proportions of significant cross-correlations are larger for the overall period and for the positive phase of the PDO compared to the negative phase. This could be because fewer stations, including stations with short records, are available for the negative phase period compared to that for the positive phase. However, longer datasets would certainly have led to more robust results. The results of a similar analysis for the PDO–winter flow correlations are presented in Fig. 3b. The proportion of stations with significant positive (negative) correlations varies between 3% (1%) and 18% (2%). Hence, similar observations as for the station cross-correlations can be made from this figure as well; that is, the majority of stations are positively correlated with PDO but the opposite behavior cannot be ruled out.

Average percentage differences and standardized differences from the period-of-record mean are shown in Figs. 4a and 4b. The original sign (i.e., positive and negative) of the PDO indices is retained in these figures regardless of the long-term positive and negative signs of the PDO phases for 1943–76 and 1977–2007. Both figures suggest that one can expect below (above) average winter flows during the negative (positive) phase of the PDO. The relationship of PDO–winter flows is statistically significant at a level as large as 0.006. However, the relationship shown in the latter plot is slightly stronger than the former. It is discussed above that winter flows for some of the stations are negatively correlated with PDO; therefore, it was tempting to study the influence of these stations on the PDO–winter flow associations. By excluding stations, which are negatively correlated with the PDO, the PDO–winter flow association, shown in Fig. 4, strengthens further; that is, the relationship becomes significant at a level of 0.0004. Similar plots were also developed separately for the group of stations for which the winter flows are negatively correlated with the PDO indices. This group was formed if either the PPMC or SROC was found negative: 27% of stations studied fall in this group. The results for PDO–winter flow associations for this group of stations are shown in Fig. 5. This figure demonstrates a very weak relationship between PDO indices and winter flows for the overall period and for the individual periods corresponding to negative and positive phases, suggesting that (i) winter flows at these stations were influenced by factors other than the PDO variability and (ii) it is very likely that inclusion of these stations in the analysis would contaminate the true PDO–winter flow relationship. Of the 27% of stations negatively correlated with the PDO, 75% are located in British Columbia and Alberta with the majority in the Montane Cordillera ecozone, while the remaining 25% are located in the Yukon (one station) and Northwest (six stations) Territories and Alaska (two stations). The spatial distribution of correlations for the majority of these stations was found very similar to that shown in the work of Woo and Thorn (2008) on a monthly time step.

c. LTP structure of PDO indices

A detailed analysis of the serial structure of the monthly, annual, trimonthly, six-monthly, and December–March mean monthly values of the PDO indices is carried out by estimating the Hurst exponent using the RARS, ASD and DFA methods. The results of this analysis are shown in Fig. 6. For each of the three methods, estimates of the Hurst exponent were obtained from the slopes of the linear relationships fitted by a least squares technique. For the DFA method, the effects of local trends of first-degree polynomial type were considered only. Though possible, the effects of trends of higher degree polynomials were not analyzed. Estimates of the 90% and 95% confidence intervals were obtained by developing distributions of the Hurst exponent, in a similar manner as for the PDO indices, from 10 000 Normal (0, 1) samples and estimating 2.5th, 5th, 95th, and 97.5th percentile values from the resulting distributions. The results of all three methods strongly suggest that the serial structure of PDO indices resembles that of a stochastic process with LTP. This observation has serious implications for the analysis of trends in time series of hydroclimatic variables in NWNA. If we believe that climatic fields (i.e., temperature and precipitation) are influenced by the PDO, which possesses characteristics of a LTP-like stochastic process, then it is also plausible to speculate that the LTP mechanism could be present in time series of winter flows through the “PDO–climatic fields” physical link, though its intensity could become modified because of basin characteristics such as geology, landforms, vegetation, and others.

d. LTP structure of winter flows

A similar analysis as presented above for the PDO indices was carried out for winter flows observed at 16 selected stations with relatively longer records ranging from 59 to 65 yr. The results of this analysis, shown in Fig. 7, demonstrate that the presence or absence of LTP depends, to some extent, on the choice of the method used for estimating the Hurst exponent. According to the results of the RARS, ASD, and DFA methods, estimated values of the Hurst exponent are >0.5 for 16, 10, and 14 stations (out of 16). Thus, if only the H > 0.5 criterion is used, then the majority of the 16 stations may have weak to strong LTP. However, if sample size uncertainty is taken into account by considering approximate 90% confidence intervals, then there are only 7 (according to the RARS), 6 (according to the ASD), and 7 (according to the DFA) stations (out of 16) where LTP could be present. Compared to the very strong possibility of LTP in PDO indices, the results of winter flows for the selected 16 stations portray a moderate picture only. Thus, it is reasonable to suspect that the strength of LTP in winter flows may have been modified owing to basin characteristics or other unknown factors. Given this assertion, it would be plausible to study trends in winter flows with multiple stochastic assumptions (i.e., IND, STP, and LTP) to reveal a more detailed picture of trends in winter flows and that is what has been presented in the section to follow.

There are 72 stations that satisfy the “station inclusion criterion” for the analysis of trends. Statistical significance of trends in time series of winter flows observed at these stations was evaluated under the following three categories.

  1. The trends and their statistical significance were studied with the IND hypothesis, that is, by assuming that the temporal structure of winter flows does not correspond to that of a STP- or LTP-like stochastic process and hence the MK test was applied without any modification.
  2. A deterministic linear trend and a STP-like stochastic structure of type AR(1) were assumed. The AR(1) structure was assumed only if the trend-independent value of the first autocorrelation coefficient was greater than and equal to 0.10. This is an arbitrary but reasonable threshold to take into account the effect of first-order autocorrelation on trend significance. A threshold of 0.05 has been used in previous studies on trends (e.g., Wang and Swail 2001). For the implementation of effective sample-size-based modified MK test (i.e., MMK2), an iteratively estimated value of _r_1 was used. For the MMK1 test, an iteratively estimated trend was removed before estimating _r_1 from ranks of the observations. For semiparametric methods, based on the time series modeling and simulation approach, 10 000 simulated samples were used to estimate significance of the observed MK test statistic S.
  3. First, a deterministic linear trend was assumed and that was optimized in the presence of a LTP-like stochastic structure of type FARIMA (0, d, 0). Second, a regression-type linear trend was estimated without reference to the stochastic structure of the time series. Third, a nonparametric “Sen slope”-type linear trend was estimated also without reference to the stochastic structure of the time series. The optimized or the linear-model-type or the Sen slope–type deterministic component was removed, one at a time, before estimating trend significance with the MKS method. These methods are referred to as MKS-FM, MKS-LM, and MKS-Sen, respectively.

Combined results of the signs of the MK test statistic S, Sen slope estimator, slope of the linear regression model, iteratively estimated slope, and optimized slope indicate that at least (most) 17% (22%) stations show a negative trend and 78% (86%) stations show a positive trend in winter flows. Hence, it is possible to conclude that NWNA is dominated by positive trends. The interaction of deterministic and stochastic components is studied in Fig. 8, wherein the values of _r_1 estimated from raw and trend-removed data for three different approaches are compared. The results of this figure demonstrate that a deterministic trend component, if present, would bias the _r_1 estimate and hence it is important to estimate both deterministic and stochastic components in a joint manner to overcome their influence on each other.

The total number of stations with significant trends, noted at four levels of significance (i.e., 0.001, 0.01, 0.05, and 0.1), is shown in Fig. 9. Stations where trends are significant at the 0.001 (0.1) level can be seen as having very strong (very moderate to very strong) trends. Similarly, stations where trends are significant at the 0.01 (0.05) level can be seen as having strong (moderate) to very strong trends. The effect of various structural assumptions on trend significance is obvious from the results shown in Fig. 9. Seven stations show very strong trends with the IND hypothesis and that number reduces to zero with the STP–LTP joint hypothesis. Similarly, at the 5% significance level, the MK test suggests a significant trend for 28 stations compared to an average number of 23 by the STP-based tests alone and 11 by the STP/LTP-based combined tests. Though there are slight differences between the results of the STP-based methods, the number of stations with significant trends was found almost similar except for the MMK2 test, which seems out of sync. The modification proposed by Yue and Wang (2004) for this test, based on the original work of Bayley and Hammersley (1946), was questioned earlier in the work of Hamed and Rao (1998), who proposed the MMK1 test. After a comparison with the STP-based semiparametric tests, it seems that the MMK1 test provides more reliable results than the MMK2 test. The results of the semiparametric tests are mutually comparable. In addition to considering the effect of lag-1 autocorrelation on trend significance, these tests allow one to approximately take into account the effect of the underlying true distribution of observational records, which is unknown in almost all cases.

Field significance of trends in NWNA was assessed using the false discovery rate approach (Wilks 2006; Khaliq et al. 2009b). For the implementation of this method, 5% station/local and field significance levels were used. Trends in winter flows in NWNA were found field significant for the MK, MCBS, MCPOL, MCLOG, MCGEV, MCWAK, MMK1, and MPW testing methods only. The maximum (i.e., 13) number of discoveries, which is an indirect indication of the strength of the signal, was found for the MK testing procedure and the minimum (i.e., 4) for the MCLOG procedure. When the joint STP–LTP hypothesis was invoked, trends were not found field significant, even at 10% level. The results of the MMK2 testing method were found similar to those for the case of joint STP–LTP hypothesis. In summary, results of the analysis reveal considerable influence of various stochastic assumptions about the serial structure of time series on trend significance. This conclusion suggests that it is important to explore fully and to take into account the effect of serial structure of time series for a detailed investigation of trends.

5. Discussion and concluding remarks

It is shown in this study that mean monthly winter (January–March) flows have statistically significant associations with December–March mean monthly PDO indices. This observation is in line with the common belief that NWNA is a PDO-sensitive region. Though this observation is valid across NWNA in a general manner, discrepancies have been noticed for some stations, which have negative associations with the PDO indices; these stations tend to deteriorate the true PDO–winter flow relationship. Most probably, winter flows at these stations and at some other stations where the PDO signal is weak are affected by factors other than the PDO variability, and hence one could speculate that local climatic factors (such as altitude, topography, distance from large water bodies, among others) may have contributed to such behavior. In addition to other investigators, this point was recently raised in the work of Woo and Thorn (2008), which seems to be a valid argument; that is, areas with complex topography can complicate the relation between the regional climate signal and streamflow in specific locations. According to Glass (1999), glacier-fed rivers have much different runoff patterns than nonglacier-fed rivers, even if their basins contain as little as 5% glacier ice. This factor may also influence PDO–winter flow relationships. Additional discussion on the physical processes that might support a cause-and-effect between PDO and winter flows can be seen in Neal et al. (2002) and Woo and Thorn (2008) and references therein.

Also Woo and Thorn argue that trends induced due to decadal shifts of the PDO phenomenon are more important compared to period-of-record linear trends, which are largely driven by climate change concerns and have been the focus of numerous studies over the past several years. It is shown in this paper that the serial structure of PDO indices conforms to a stochastic process with very strong LTP. Trends induced due to long excursions of the mean level of the process are regular features of a stochastic process with LTP (Koutsoyiannis 2003, 2006). Thus, it is still a reasonable hypothesis to investigate period-of-record linear/nonlinear trends as a deterministic component superimposed by a stochastic component of type STP or LTP, as demonstrated in this paper and also in the work of Cohn and Lins (2005). Therefore, the influence of transient trends such as the ones induced owing to decadal shifts of the PDO phenomenon can be taken into account when evaluating the statistical significance of period-of-record trends driven by climate change concerns, provided appropriate stochastic assumptions about the serial structure of the time series be invoked, for example, the LTP-like hypothesis. However, the problem one faces is that it is extremely difficult to investigate the presence of LTP, particularly in short, noise-corrupted, time series of observational records. Hence, there is a possibility that one would misdiagnose LTP as STP or no persistence at all. An earlier work of Lettenmaier and Burges (1978) along similar lines is an extremely useful reference. Also, Brabets and Walvoord (2009) discuss the suitability of analyzing period-of-record linear trends in the presence of decadal shifts due to the PDO phenomenon.

Another important point, which was raised recently in the work of Chen and Grasby (2009), is the starting point of a time series selected for investigation of trends. If a fluctuating phenomenon like the PDO is believed to be the source of local climate variability, then it is desirable to have the starting points of all time series, included in a trend investigation study, as close as possible. Time series that span different lengths of the PDO cycle may result in misleading sign (i.e., positive or negative) of the period-of-record trends; particularly, this aspect will be of considerable importance if the trends are believed to be driven by climate change concerns. To avoid the possibility of such errors, trends in the present study were investigated for nearly a full cycle of the PDO phenomenon. A similar practice in the future would help produce robust conclusions.

During the winter period from January–March, when most of the streams are ice-covered and terrestrial surface layers are frozen, groundwater is assumed to be the major input to many rivers in NWNA (Walvoord and Striegl 2007). According to the results of trend analysis reported in the present study for 72 stations, which have at least 50 years of observations, winter flows are increasing at about 82% of stations. This means that the groundwater contribution is generally increasing over time. It is suggested in the work of Brabets and Walvoord (2009) that such an increasing pattern could be the result of a warming effect, particularly for the Yukon River basin, which is a part of the current study as well. However, the trends and their significance obtained with the STP–LTP joint hypothesis for the Alaskan stations alone were found to be field significant at a level of about 25%. Thus, the current results are not found to be as strong as suggested earlier. Compared to this, similar results obtained with the STP–LTP joint hypothesis for the whole of NWNA were found field significant at about the 13% significance level. For both STP- and IND-based hypotheses, the results of trend analysis were found field significant at a level less than 5%. These observations clearly demonstrate the influence of various stochastic assumptions, about the serial structure of time series, on the significance of trends.

In summary, the contribution of this paper would be a reassessment of the effects of the PDO phenomenon on winter flows in NWNA using a much larger dataset than the ones previously explored and an improved understanding and interpretation of temporal changes in winter flow regimes. This target is achieved by adopting and introducing new tests of trend analysis that allow incorporating the effects of various serial structures on trend significance. Though the current study concentrates on temporal changes in the mean level of the underlying observation-generating mechanism, it was felt that an equally important area of investigation is the temporal change in variability of winter flows, which is found to be increasing over time at numerous stations. The PDO is not the only climate pattern that influences streamflows in NWNA, the Arctic Oscillation and other climate patterns, such as the El Niño–Southern Oscillation, also may influence streamflows but are not studied in this paper. Such possibilities along with the analysis of other seasonal flow regimes will be explored in future studies. Lastly, it is worrisome that the trend significance is sensitive to the serial structure of time series; this raises serious concerns about the usefulness of univariate trend analysis approaches. Perhaps, a multivariate setting, where various variables involved interact with each other, may be a better alternative to address the issue of change detection in time series of instrumental records.

Acknowledgments

The financial support provided by the Canadian Treasury Board to the Adaptation and Impacts Research Section of the Atmospheric Science and Technology Directorate of Environment Canada is gratefully acknowledged. The authors thank Harry Lins from the U.S. Geological Survey and an anonymous referee for their helpful comments.

REFERENCES

APPENDIX

The Mann–Kendall Test

The Mann–Kendall test searches for a trend in a time series without specifying whether the trend is linear or nonlinear. It is based on the test statistic S, defined as

i1525-7541-11-4-917-ea1

i1525-7541-11-4-917-ea1

where yi and yj are the sequential data values, N is the total number of data values in the time series, and

i1525-7541-11-4-917-ea2

i1525-7541-11-4-917-ea2

A positive (negative) value of S indicates an upward (downward) trend. For N ≥ 8, Mann (1945) and Kendall (1975) have documented that the statistic S is approximately normally distributed with the mean, E[_S_] = 0, and variance

i1525-7541-11-4-917-ea3

i1525-7541-11-4-917-ea3

where ti is the number of ties of extent i (i.e., the number of data in the tied group) and n is the number of tied groups. The standardized test statistic

i1525-7541-11-4-917-ea4

i1525-7541-11-4-917-ea4

follows the standard normal distribution. At α significance level, the null hypothesis of no trend is rejected if the absolute value of Z is greater than the theoretical value Z_1−_α/2.

Fig. 1.

Fig. 1.

Fig. 1.

Map of the study area and various climate zones according to Hare and Thomas (1979), figure adapted from Woo and Thorn (2008), where the location of 110 stations used in their study is shown.

Citation: Journal of Hydrometeorology 11, 4; 10.1175/2010JHM1254.1

Fig. 2.

Fig. 2.

Fig. 2.

(a) Time series of December–March mean monthly PDO indices for the period 1900–2007, along with negative (cold)- and positive (warm)-phase partitions obtained with the help of statistical tests of change-point detection. (b) Plots of cumulative departures from the period-of-record mean for the monthly, annual, and December–March averaged times series of PDO indices. Linear trends for each partition are shown in (a) and ending years of each partition (vertical lines with year labels) are shown in both (a) and (b).

Citation: Journal of Hydrometeorology 11, 4; 10.1175/2010JHM1254.1

Fig. 3.

Fig. 3.

Fig. 3.

Frequency distributions of the Pearson product moment correlation (PPMC) and Spearman rank order correlation (SROC) for the entire period 1943–2007 and the periods corresponding to cold and warm phases of the PDO (i.e., 1943–76 and 1977–2007). Results for the station cross-correlations are shown in (a) and those for the PDO–winter flows are shown in (b). Proportion of significant (at 5% level) positive and negative (in brackets) correlations are shown inside each panel.

Citation: Journal of Hydrometeorology 11, 4; 10.1175/2010JHM1254.1

Fig. 4.

Fig. 4.

Fig. 4.

(a) Average percentage differences and (b) standardized differences from the period-of-record mean are shown: (left) time series plots of these differences and (right) the relationships of these differences with the PDO indices. The relationships are significant at the 5% significance level.

Citation: Journal of Hydrometeorology 11, 4; 10.1175/2010JHM1254.1

Fig. 5.

Fig. 5.

Fig. 5.

As in Fig. 4, but for the group of stations negatively correlated with the PDO indices. These relationships are strongly nonsignificant.

Citation: Journal of Hydrometeorology 11, 4; 10.1175/2010JHM1254.1

Fig. 6.

Fig. 6.

Fig. 6.

Estimates of the Hurst exponent (symbols) using the rescaled adjusted range statistic (RARS), aggragated standard deviation (ASD), and detrended fluctuation analysis (DFA) methods for (a) monthly and annual average, (b) trimonthly average, (c) six-monthly average, and (d) December–March average values of the PDO indices for the period 1900–2007, along with 90% and 95% confidence intervals (small and large horizontal bars joined with vertical lines), obtained using Monte Carlo simulation technique using Normal (0, 1) samples of the same length as the PDO indices.

Citation: Journal of Hydrometeorology 11, 4; 10.1175/2010JHM1254.1

Fig. 7.

Fig. 7.

Fig. 7.

Estimates of the Hurst exponent (circles) using the RARS, ASD, and DFA methods for 16 selected stations (station names are shown along the x axis) with long-term winter flows, along with 90% and 95% confidence intervals (i.e., small and large horizontal bars), obtained using Monte Carlo simulation techniques.

Citation: Journal of Hydrometeorology 11, 4; 10.1175/2010JHM1254.1

Fig. 8.

Fig. 8.

Fig. 8.

Comparison of lag-1 autocorrelation estimated from raw and trend-removed data (OPT: lag-1 autocorrelation estimated by maximizing the likelihood function of an AR(1) process, ITE: lag-1 autocorrelation estimated using an iterative procedure described in the text, MPW: lag-1 autocorrelation estimated using the MPW approach).

Citation: Journal of Hydrometeorology 11, 4; 10.1175/2010JHM1254.1

Fig. 9.

Fig. 9.

Fig. 9.

Number of stations with (a) very strong, (b) strong to very strong, (c) moderate to very strong, and (d) very moderate to very strong trends in time series of winter flows observed in NWNA.

Citation: Journal of Hydrometeorology 11, 4; 10.1175/2010JHM1254.1

Table 1.

List of NWNA hydrometric stations used in the study. Stations (72) marked with an asterisk have record lengths ≥50 yr and were used in the analysis of trends. Stations marked with a tilde were discarded because of a large number of missing values.

Table 1.

Table 1.