The data used are the monthly PDO index and the SST for 1961– 2010, which are published on the website of the National Oceanic and Atmospheric Administration (NOAA). In this paper, the SST data have a resolution ratio of 2^° × 2^° , selecting 20^°   N– 60^°   N, 140^°   E– 120^°   W; i.e., the covering range of the PDO model. In the actual calculation, the years from 1981 to 2010 are taken as the reference climate state, and the anomaly of original observation data is used as a reference climate state for a calculation sequence.

The variance and autocorrelation coefficient of the sequence are calculated through moving in the research in order to find the early warning signal of an abrupt change of the PDO index sequence. The specific method is shown in Fig.  1.

Fig.  1.

	Figure Option ViewDownloadNew Window
	Fig.  1. We calculated the variance and autocorrelation coefficient by sliding the window.

a) We calculated variance by sliding the window, L₁, L₂, L₃, … , L_n … denote windows of the same length, S₁, S₂, S₃, … , S_n … represent the variances of the corresponding windows, L is the total length of the sequence, and M_T is the sliding step; b) we calculated the autocorrelation coefficient by sliding the window, L₁(L₁₂), L₂(L₂₂), L₃(L₃₂), … , L_n(L_n2)… represent windows of the same length, α ₁ denotes the autocorrelation coefficients of L₁ and the L₁₂, α ₂ denotes the autocorrelation coefficient of the L₂, while L₂₂, … the α _n refer to the autocorrelation coefficients of the L_n and L_n2, L_T represents the lag time, and L and M_T have the same meaning as those in Fig.  1(a).

(iii) Relationship between the critical slowing down and the increased autocorrelation coefficient

A critical slowing down will tend to increase both the autocorrelation coefficient and the variance of the fluctuations in a stochastically forced system when the system approaches to a bifurcation at the threshold of the control parameter.^{[21, 22]} It is assumed that there is a repeated disturbance of the state variable after each period Δ t (additive noise). Between disturbances, returning to equilibrium is approximately exponential with a certain recovery speed λ , which can be described with an auto-regressive model^{[21, 22, 33]}

where x_n is the deviation of the system variable from an equilibrium state, ε _n is a random term with standard normal distribution, i.e., white noise in the system, and s represents the standard deviation. If λ and Δ t are independent of x_n, then this model can also be written as a 1-order auto-regressive (AR(1)) process:

where the autocorrelation c = e^{λ Δ t} is 0 for white noise and close to 1 for red noise.

The analysis of variance to AR(1) process leads to

Generally speaking, in the process where a system tends to its critical point, the recovery rate from small-amplitude disturbance will become increasingly slower.^{[33– 36]} When the system is closer to the critical point, the recovery rate λ tends to 0, and the autocorrelation item c approaches to 1. Thus, the increased autocorrelation coefficient can be considered as an early warning signal that the system is approaching a critical point. In the calculation of this study, s represents the variance of the whole sequence, which is a constant value, while the Var value in Eq.  (6) changes with the window size and moving step.

(iv) Lyapunov exponent

The Lyapunov exponent^{[37, 38]} is an important physical quantity that is used to describe whether a system is in a regular motion or a chaotic motion. When at least one of the Lyapunov exponents of a system is positive, the system is considered to be chaotic. When all of the Lyapunov exponents of a system are negative, the system is considered to be convergent. Therefore, the LLE of a system is often calculated in order to determine whether the system is in a regular motion or a chaotic motion. The currently used methods for calculating the largest Lyapunov exponent of a chaotic time sequence includes the following:^[39] the Nicolis method, the Jacobian method, the Wolf method (phase-space reconstruction method), the P-norm method, and the small-data method. Of these methods, the most common and widely used are the Wolf and small-data methods. This study adopts the Wolf method to calculate the LLE of the sequence. The specific operation steps of the Wolf method are described as follows.

The method of extracting the Lyapunov exponent from the time sequence of a single variable is based on the phase-space reconstruction of the time sequence. Wolf et al.^[37] made an estimation of the Lyapunov exponent directly based on the evolution of the phase trajectory, phase plane, and phase volume, and so on. These methods are collectively referred to as the Wolf method, which is widely used to study the chaos and chaos time sequence predictions based on the Lyapunov exponent.

Assume that the chaos time sequence is x₁, x₂, … , x_k, … , x_N; the embedded dimension is m; and the time delay is τ . The reconstructed phase space will then be

In the actual calculation of this study, the embedded dimension m is 3 and the time delay τ is 12. Since the data used in this study are the monthly data, it is defined that the data of the current month has the smallest difference from the data of the corresponding month of the previous and following years. Take the initial point Y (t₀), and assume that its distance with the nearest neighbor point Y₀ (t₀) is L₀. Then trace the time evolution of these two points until the t₁ time, when their spacing exceeds a certain value ε > 0, . Next, retain Y (t₁). Then find another point Y₁ (t₁) near Y (t₁) to make , the included angle with respect to it will be as small as possible. Finally, continue the above procedure until Y(t) reaches the final point N of time sequence. At this point, trace the total iterative times of the evolution process, and then the largest Lyapunov exponent will be as follows:

In this study, the LLE of the PDO index sequence at different time intervals is calculated through the Wolf method, which is used to determine the chaotic characteristics of the PDO index sequence at different times. The monthly data of the PDO index sequence during the period from 1951 to 2010 are selected, which totals 720. The LLE during the period of 1951 to 2010 is calculated to be 1.23 by using the Wolf method, and it is defined as the total LLE. Then, the LLE at each time interval (every five years) is further calculated and compared with the total LLE. The method of calculating the LLE at each time interval is as follows. First, a PDO time sequence (raw data minus the average and then divided by the mean square deviation) is standardized. Then, a random function is used to generate 1000 random arrays with a length of 60, followed by the averaging of these 1000 random arrays to obtain a random array with a length of 60. Then, the standardized treatment is realized. The calculated random number is used to gradually replace the value of the PDO index sequence at each corresponding time interval (gradually replace the values at the time intervals of 1951 to 1955, 1956 to 1960, … , 2001 to 2005, 2006 to 2010). The Wolf method is then utilized to calculate the LLE of the new sequence. The LLE of the PDO sequence is obtained at each time interval through the described replacements. The LLE reflects the chaotic characteristics of the system. It is determined that a larger exponent value will lead to stronger chaos characteristics.

Figure  2 shows the information of the PDO index sequence, PDO has obvious inter-decadal variations. The research results^{[7, 13, 40, 41]} show that the cold and warm phase positions of the PDO experienced an abrupt transition (abrupt change) in 1976. The movement average trend line in Fig.  2 also illustrates that the phase position of the PDO before 1976 was at a cold phase position, and immediately changed to a warm phase position after 1976. Therefore, the PDO experienced a transition from cold to warm phase positions in 1976. This study mainly discusses the early warning signal of this abrupt change in 1976, retains the sequence before the abrupt change in 1976, and calculates the variance of the sequence, such as the length of parameter L, which is 312 months (monthly data within 1951 to 1976) before the abrupt change, and the autocorrelation coefficient in order to determine the early warning signal of an abrupt change, as shown in Fig.  1.

Figure  3 shows the early warning signals of the PDO index with different windows and moving step lengths. Figure  1 shows the variances of the moving calculation sequence and the autocorrelation coefficient, as well as the selection of the window moving and step length. For the purposes of this study, L₁, L₂, L₃, … , L_n … , respectively, include 60 months (5 years) and 120 months (10 years). Similarly, L₁(L₁₂), L₂(L₂₂), L₃(L₃₂), … , L_n(L_n2)… also, respectively, include five years and 10 years. The moving step length M_T is taken to be 1 month, while the lag time L_T is taken to be 1 month. Figures  3(a) and 3(b) show the variance signals of the PDO index sequence. The moving step length (M_T as shown in Fig.  1) indicates that the sequence with the selected window size would slide backwards by a fixed step length in order to obtain a new sequence, and then the variance of this new sequence could be calculated. The arrows in Figs.  3(a) and 3(b) show that the variance started to increase gradually around the year 1970 (figure  3(a) shows 1970, while figure  3(b) shows 1969). It can thus be known from the above theoretical analysis that the critical slowing down causes the lower internal change rate of the system, and the system’ s state at any time is very similar to its previous state. Therefore, the autocorrelation coefficient tended to be 1, and the variance according to Eq.  (6) tended to be infinity. In principle, the critical slowing down reduces the system’ s ability to track the fluctuation and has an opposite effect on the variance. The phenomenon of critical slowing down, for example, the variance increase occurs when the climate system approaches to its critical point, can be regarded as an early warning signal for an abrupt change of the climate system, or as an early warning signal for the PDO’ s abrupt change which occurred around 1970. Therefore, it can be seen that the interval between the occurrence time of the early warning signal for the PDO’ s abrupt change and the occurrence time of an abrupt change is six years. This is a good indicator for climate predictions of a future abrupt change.

Fig.  2.

	Figure Option ViewDownloadNew Window
	Fig.  2. Curves of PDO index sequence changing with time (The bars represents the value of the PDO index: red for positive, blue for negative; and, a black curve for the trend information extracted from moving average).

Fig.  3.

	Figure Option ViewDownloadNew Window
	Fig.  3. Critical slowing down phenomena in different window lengths and different sliding steps of PDO index, showing the signal of variances in (a) 5-year and (b) 10-year window lengths, and the signals of autocorrelation coefficient in (c) 5-year and (d) 10-year window lengths.

Figures  3(c) and 3(d) show the autocorrelation coefficient detection of the PDO index sequence. It may be noted that the lag time (LT as shown in Fig.  1) and the moving steps (MTs as shown in Fig.  1) have different meanings in Figs.  3(c) and Fig.  3(d). The lag time refers to the original sequence, where the selected window size will be lagged for a selected step length to obtain another sequence with the same length, and to calculate the correlation (which is its own lag correlation) between the obtained sequence and the previous sequence. However, the moving step is the same as that of the variance signal. The arrows in Figs.  3(c) and 3(d) show that the autocorrelation coefficient starts to increase gradually (figure  3(c) shows 1967, while figure  3(d) shows 1970), the critical slowing down causes the lower internal change rate of the system, and the system’ s state at any time is very similar to its previous state. Therefore, the autocorrelation coefficient will tend to be 1. The phenomenon of the critical slowing down, such as the variance increase, increase of the autocorrelation coefficient, etc. when the climate system approaches to its critical point shows that the climate system will experience an abrupt change. Therefore, we can conclude that the early warning signal occurred in, approximately, 1970. We then determine that the time of early warning signal, which is discovered through the variance and autocorrelation coefficient, is the same. This further demonstrates the feasibility of finding the early warning signal of an abrupt change on the basis of the phenomenon of the critical slowing down.

Following the comparison of Figs.  3(a) and 3(b) with Figs.  3(c) and 3(d), it is easily determined that the different window sizes and moving steps would have an influence on the stability of the results. When the data volume is fixed, the larger window size and longer moving step mean that result is more stable. It is also more convenient to find the early warning signal of the abrupt climate change. References  [25], [26], and [42] contain a detailed description regarding the influences of the window size and moving step on the stability results and the relevant research on the identification of the spurious signal. For example, when the data volume is fixed, the larger window size and longer moving step give a stabler result. Different lag steps with the same window also influence the result, and the selection of the window and lag step makes the signal disappear, and also influences the signal stability.

The preceding sections of this study test the early warning signal of the PDO abrupt change based on the phenomenon of the critical slowing down. However, the chaos characteristics and source of the early warning signal of the system approaching to its critical point are still a problem that urgently needs to be solved. The Lyapunov exponent is an important physical quantity which reflects whether the system is in a regular motion or a chaotic motion. Therefore, the LLE at each time interval should be calculated in order to investigate the intrinsic nature of the system, and to study its chaos characteristics.

Figure  4 illustrates the LLEs at each time interval and the total LLE of the PDO. It can be seen from Fig.  4 that the total LLE is 1.23 > 0, which indicates that the PDO system at this time interval is chaotic. It is found that a larger LLE value will make the chaotic characteristics of the system stronger. Only the LLEs within the time intervals of 1951– 1955, 1971– 1975, and 1996– 2000 are smaller than the LLE of the total sequence. We can obtain the specific calculation method of the LLE at each time interval by using the above theory. The random array replaces the original sequence at each time interval, and makes the LLE smaller, which means that the chaotic characteristics of the original sequence at these three time intervals are stronger. For example, the chaotic characteristics at the time intervals of 1951– 1955, 1971– 1975, and 1996– 2000 are stronger. The abrupt climate change refers to a phenomenon^[1] where the climate crosses through a critical threshold and transits from one stable state (or a stable and continuous change trend) to another stable state (or a stable and continuous change trend). Through a comparison of the several time intervals between before and after 1976, in combination with the analysis of Figs.  2 and 4, the LLE during the period from 1956– 1970 is found to be positive. That is to say, the chaotic characteristics of the original sequence during this time interval is weaker. Similarly, the chaotic characteristics of the original sequence is stronger at the period of 1971– 1975, and weaker at the period of 1976– 1995. This means that the original sequence enters into a chaotic state from 1971 to 1975. Therefore, the phenomenon of a critical slowing down could be detected within this time interval, as could the phenomenon of early warning signal. Similarly, the chaotic characteristics of the original sequence are also stronger during the periods of 1951– 1955 and 1996– 2000. By comparing with the PDO index sequence shown in Fig.  2, it can be determined that the PDO transits from a cold phase to a warm phase position in, approximately, 1960 and 2000. Therefore, the chaotic characteristics of the original sequences at intervals of 1951– 1955 and 1996– 2000 are stronger. The conclusion that the PDO experiences an abrupt change during the years of 1960 and 2000 has been demonstrated in previous studies.^{[1, 7, 13, 43]} However, no consensus has been reached for these two abrupt changes, as has been reached for the abrupt change in 1976, due to the consideration of the data length and the duration of time which has passed. In the same way, the early warning signal of the PDO’ s abrupt change during 2000 can be obtained by calculating the variance and autocorrelation coefficient of the sequence. However, the early warning signal of the abrupt change in 1960 cannot be found by using the above method, due to problems with the data length and the moving window (figure omitted).

Fig.  4.

	Figure Option ViewDownloadNew Window
	Fig.  4. The LLEs of PDO index in each time period and total time series (The bars represent the LLEs of the PDO index in each time period, thick black dotted line is for the LLE of the total series).

Based on the above factors, the LLE can successfully reflect the inherent chaotic characteristics of the PDO system and the changes of the chaotic characteristics shown at the different time intervals, thus providing the theoretical basis for the detectability of an early warning signal. In addition, the LLE further validates the feasibility of regarding the phenomenon of a critical slowing down as an early warning signal of an abrupt change.

The PDO is a type of long-term existing Pacific climate change model. The PDO index of the United States National Climatic Data Center is expanded and reconstructed on the basis of the SST field of the NOAA, and the empirical orthogonal function (EOF) decomposition is conducted after removing the global warming trend in the SST field in the north of Pacific at 20^°   N. Its first principal component is defined as the PDO index.^{[27, 44]} Therefore, the early warning signals of the PDO’ s abrupt change should be further studied from the aspect of the background field of SST. Figure  5 shows the relationship between the phenomenon of the critical slowing down caused by the PDO’ s abrupt change and the SST change rate in the corresponding areas.

Fig.  5.

	Figure Option ViewDownloadNew Window
	Fig.  5. Relation among the SST, increased autocorrelation, and increased variance of PDO index.

It can be concluded that the early warning signals of these two phenomena of critical slowing down (increase of autocorrelation coefficient and increase of the variance) before an abrupt change of the PDO occur at the same time. The increase of the autocorrelation coefficient of the PDO index means that the transition of the PDO between the positive and negative phase positions is more regular. For example, the spatial distribution-type change of SST in the area is more regular. Therefore, the change of SST at each lattice point of the area is also more regular, and the autocorrelation coefficient of SST change at each lattice point is also increased. The autocorrelation coefficient probability distribution diagram is given for the SST at each time interval to represent the autocorrelation coefficient of the sequence at each lattice point of the space. The characteristics of the increase of the autocorrelation coefficient of SST indicate a left advertence in the autocorrelation coefficient probability distribution diagram, and its peak value appears to be large.

In the same way, the increase of the PDO index variance indicates the increased change amplitude of the PDO phase position. For example, it transits from an obvious positive (negative) phase position to an obvious negative (positive) phase position. This also indicates the increased change amplitude of the spatial distribution pattern of the SST in the area, or the increased SST change amplitude at each of the lattice points of the area. Therefore, the variance of the SST change at each lattice point also increases, while the variance difference of the SST in the corresponding area decreases (the change rate of positive and negative phase position of the SST was synchronous). This is also reflected by the distribution diagram for the SST variance coefficient probability rate at each time interval, and its probability distribution rate indicates a centralized probability distribution and a narrow value range.

Fig.  6.

	Figure Option ViewDownloadNew Window
	Fig.  6. The Distributions of the autocorrelation coefficient of the SST.

Figure  6 shows the spatial distribution diagrams for the autocorrelation coefficient (the lag time of the autocorrelation coefficient of the SST was − 1 in this study; i.e., lagging for one month, which corresponds to the autocorrelation coefficient calculation of the PDO) of the SST at time intervals of 1966– 1970, 1971– 1975, 1976– 1980, and 1981– 1985. This study focuses on the early warning signals of the PDO’ s abrupt change before and after 1976, thus only these four time intervals are selected for analyses. Similarly, the spatial distribution diagram could also be made for the autocorrelation coefficient at each time interval. It is known from the analysis in Fig.  5 that the autocorrelation coefficient of the SST in the corresponding area also increases if the autocorrelation coefficient of the PDO index augments. It may be concluded from Fig.  6 that the autocorrelation coefficient of the SST at the time interval from 1971 to 1975 is larger than that at the other time intervals. In other words, there is a critical slowing down phenomenon when the autocorrelation coefficient is increased, but such a trend is not obvious. Therefore, in order to explain this problem, further analysis should be conducted by using the probability distribution diagram for the autocorrelation coefficient at each time interval, as shown in Fig.  7. In Fig.  5, if the autocorrelation coefficient of the PDO index sequence is increased, then the autocorrelation coefficient of the corresponding SST is also increased. Moreover, the probability distribution diagram for the autocorrelation coefficient of SST shows a left advertence, and its peak value is large. Comparisons among these characteristics and the creation of the contrastive analysis on the probability distribution diagram for the autocorrelation coefficient at each time interval are shown in Fig.  7 (the autocorrelation coefficient probability distribution in the total time interval of 1951– 2010 is close to the normal distribution). It can be easily determined that the characteristics at the time intervals of 1971– 1975 and 1996– 2000 comply with the judgment basis, which illustrates that in the SST during 1971– 1975 there appears a phenomenon of a critical slowing down. Therefore, these results demonstrate that the reliability and reasonability of an early warning signal for an abrupt change in 1976 of the PDO system can be found in 1970 from the aspect of the background of SST. In the same way, the phenomenon of a critical slowing down occurs between 1996 and 2000, and this can also be regarded as the early warning signal of an abrupt change^{[7, 43]} which occurs around the year 2000. From the above discussion, the phenomenon of a critical slowing down, such as the increase of the autocorrelation coefficient occurring before the PDO’ s abrupt change, can be regarded as an early warning signal of its abrupt change. Such a phenomenon shows a good correlation with the SST in its corresponding area. The analysis on the autocorrelation coefficient from the background field of the SST also demonstrates the reliability and reasonability of the determined early warning signals.

Fig.  7.

	Figure Option ViewDownloadNew Window
	Fig.  7. Probability distributions of the autocorrelation coefficient of the SST in each time period.

Like the analysis of autocorrelation coefficient, the relationship between the early warning signal of the PDO’ s abrupt change which is based on the increased variance and the variance of SST in the corresponding area, will be discussed below from the point of view of variance. Figure  8 shows the distribution diagrams for the variance of SST in the four time intervals. In accordance with the analysis shown in Fig.  5, it is known that if the PDO index sequence exhibits the phenomenon of a critical slowing down, such as a variance increase, then the SST in the corresponding area of the PDO will show an increased change rate, and the variance differences of SST in the corresponding area will be reduced. Using the comparative analysis of Fig.  8, it can be found that the variance difference of the SSTs in the area from 1971 to 1975 is smaller, and the SST model in the area is relatively stable. Figure  9 shows the probability distribution diagrams for the variance of SST at each time interval. In accordance with the aforementioned theory and the analysis shown in Fig.  5, it is known that if in the PDO system there occurs a phenomenon of a critical slowing down, such as a variance increase, then the probability distribution of the variance of the SST in the corresponding area will show a concentrated distribution (ensuring that the change rates of the positive and negative phase positions of the SST are synchronous), which means that the value range of its probability distribution diagram is relatively narrow. In accordance with this judgment basis, and by making a comparative analysis on the distribution characteristics at each time interval, as shown in Fig.  9, we can easily find that the distribution characteristics at the time intervals of 1966– 1970 and 1971– 1975 meet the requirements. This indicates that in the SST there occurs a phenomenon of a critical slowing down at this time, and it also

Fig.  8.

	Figure Option ViewDownloadNew Window
	Fig.  8. Distributions of the variance coefficient of the SST.

Fig.  9.

	Figure Option ViewDownloadNew Window
	Fig.  9. Probability distributions of the variance coefficient of the SST in each time period.

The above-mentioned research is performed by using the annual and monthly continuous SST data. Although the anomaly treatment has been made for the sequence, its sequence still has the inter-decadal, inter-annual, and seasonal change trends and it cannot reflect the inter-monthly autocorrelation and variance characteristics of SST. Therefore, the monthly change trend of the SST is obtained through a low-pass filtering treatment, before calculating the autocorrelation coefficient and variance coefficient of SST. In this study, the sequence without a low-pass filtering treatment is also calculated (figure omitted), and the result is very similar to that given earlier in this paper, which means that the change of SST is dominated by a high-frequency seasonal change signal. Therefore, the phenomenon of a critical slowing down which occurs before the PDO’ s abrupt change can be regarded as the early warning signal of its abrupt change, and this phenomenon has a good correlation with the SST in its corresponding area. The analysis regarding the autocorrelation coefficient from the background field of SST also demonstrates the reliability and reasonability of the determined early warning signal.