1 Introduction

Decadal–interdecadal climate variability affects human society by influencing marine ecosystems, fisheries, and agricultural production (e.g. Mantua et al. 1997; Chavez et al. 2003; Nigam et al. 1999). Therefore, recognizing the interdecadal variability that occurs both in the ocean and the atmosphere over large areas of the Pacific is crucial to better understand and predict changes in the climate such as global warming (Trenberth and Fasullo 2013). Among the phenomena with interdecadal timescales (> 10 years), the Pacific Decadal Oscillation (PDO) has been attracting substantial attention (Mantua et al. 1997), which highly impact the climate as well as ecosystem regime shifts (e.g. Miller and Schneider 2000; Chavez et al. 2003). The PDO is the dominant decadal-interdecadal variation in climate over the North Pacific, which is defined as the leading empirical orthogonal function (EOF) of sea-surface temperature (SST) anomalies in the Pacific north of 20° N. The associated normalized principal component is called PDO index (PDOI), providing quantitative measures to evaluate PDO. The PDO is understood as an oceanic aspect of climatic variations which consist of various physical processes, such as teleconnection from the tropics, mid-latitude air–sea interactions in the North Pacific, and the westward oceanic Rossby wave propagation (Newman et al. 2016 and references therein), with dominant temporal scales of a period of oscillation being ~ 20 years (Minobe 2000; Zhang and Delworth 2015) and ~ 50 years (Minobe 2000; Zhong and Liu 2009). Further, there is a general agreement that the Aleutian Low plays a substantial role in driving underlying oceanic variability, which could lead to PDO-like spatial SST pattern (Hasselmann 1976; Qiu et al. 2007; Liu and Alexander 2007). The Aleutian Low is the low pressure system extending above the Aleutian Islands, which is quantitatively characterized as the normalized sea level pressure (SLP) anomaly that is spatially averaged over the area 30–65° N and 160° E–140° W (North Pacific Index; NPI; Trenberth and Hurrell 1994). The NPI is known to be highly correlated with the PDOI, suggesting their close relevance. Other studies on interdecadal climate variability in the mid-latitude North Pacific center its local expression in the Gulf of Alaska, the surface air temperature of which is highly sensitive to the Aleutian Low and therefore to the PDO (Minobe 2000; Papineau 2001). Therefore, PDOI, NPI, and the surface air temperature of the Gulf of Alaska represent the most dominant decadal-interdecadal variability over the mid-latitude North Pacific.

Studies into the NPI by Minobe (1999), and the PDOI and surface air temperatures by Minobe (2000) first identified a noteworthy feature in which the climate shifts dramatically from one state to another on decadal timescales (often referred to as ‘regime shifts’) in the mid-latitude North Pacific. They further demonstrated that the three regime shifts in the 1920s, 1940s, and 1970s were accompanied by characteristic simultaneous phase reversals of 10–30 year band-passed ‘bi-decadal’ variation and 30–80 year band-passed ‘penta-decadal’ variation. On the other hand, Wilson et al. (2007) showed a contrasting feature in the surface air temperature along the coast of the Gulf of Alaska over the 1296 years from 724 to 1999, which was reconstructed from tree-ring records. This study claimed the lack of simultaneous phase reversals before the 1900s, and therefore, concluded that these characteristics are specific to the twentieth century. These two opposing viewpoints provide arguments on whether the simultaneous phase reversals in the twentieth century were by chance or of significant necessity, which remains ambiguous and controversial.

Studies on climatic regime shifts are largely constrained by the length of the available timeseries, rendering it difficult to establish a comprehensive view of the mechanism behind these opposing hypotheses on the phase-locked temporal structure of climatic regime shifts before and after the twentieth century. In the present study, the annual-mean PDOI reconstructed from tree rings over a period of 298 years (1700–1997; d’Arrigo et al. 2001) and a more reliable PDOI computed from the instrumentally observed SST over a shorter period of 119 years (1900–2018; Rayner et al. 2003) were used together in an attempt to better understand what lies behind the ‘discontinuous’ phase-locking structure in the climatic regime shifts in the North Pacific.

2 Data and methods

The 298-year (1700–1997) annual-mean PDOI (d’Arrigo et al. 2001) and the 119-year (1900–2018) annual-mean PDOI computed using the monthly SST from the Hadley Centre Sea Ice and SST dataset (HadlSST: https://www.metoffice.gov.uk/hadobs/hadisst/) (Rayner et al. 2003) are examined. The former timeseries is reconstructed using western North American tree-ring records, while the latter incorporates spatially sparse in-situ SST observations and interpolation to each grid.

Fast Fourier Transform (FFT) spectrum analysis with cosine-tapered window (two cosine lobes of 10% of the record length between which a rectangular window is interposed) is operated to the reconstructed PDOI time-series. In our analysis, we apply the Monte Carlo method, which defines significance levels with following procedures. First, 10,000 red-noise timeseries are randomly generated with the following equation:

$$y\left( t \right) = r \cdot y\left( {t - 1} \right) + \sqrt {1 - r^{2} } \cdot \varepsilon \left( t \right),$$

where \(y\left( t \right)\) is a value of red-noise timeseries at year \(t\), \(r\) a 1-year lag autocorrelation coefficient (~ 0.46 in the present case), \(\varepsilon \left( t \right)\) is a normalized random value with Gaussian distribution, which represents atmospheric white noise forcing. FFT are performed to each red-noise timeseries, and the 500-th (1000-th) largest values at each frequency are used as 95% (90%) significance level. In the present study, the frequency bands with peaks above 90% significance levels are considered significant and extracted using the bandpass filter.

2.1 Significance of phase-locks

Here we attempt to obtain quantitative measures to assess the extent at which the timeseries have a phase-locked temporal structure. To make this possible, we introduce the index, which we call ‘ratio of the phase locks’ as indicated below:

$${\text{Ratio of phase locks}} = \frac{{\text{ number of phase locks that actually took place}}}{{\text{number of maximum possible phase locks}}}.$$

The significance of the ratio of phase-locks is then evaluated using the following procedure. The phase-locks might be realized by random processes from variations with the same power spectrum density as in the reconstructed PDOI. Thus, we may consider it significant if the observed phase-locking feature appears rarely as a result of random processes. The method used for evaluating the significance in the present study consists of two stages: (1) generating pseudo-timeseries with phase randomization and (2) counting phase-locks with Monte Carlo simulation. First, the amplitudes of sine curve corresponding to the tri-decadal and the multi-decadal variations for each frequency are retrieved. A total of 1000 pseudo timeseries are then generated by taking a sum of the retrieved Fourier components with only their phases randomized (see “Appendix” for more details). The ratio of phase-locks is calculated for each tolerance from 0.0 to 6.0 years with an interval of 0.1 years, where the tolerance is defined as the extent to which the phase reversals from the same polarity to the other are considered simultaneous.

3 Results

To take advantage of the long period covered by the reconstructed PDOI, we first examine a spectrum in a frequency domain (Fig. 1a) using fast Fourier transform (FFT), which shows three 90% significant peaks centered at around 18.6 year (red vertical line in Fig. 1a), 27.8 year (blue) and 55.8 year (green) periods. Three frequency bands corresponding to these interdecadal timescale (< 1/10 cycles/year: CPY) are then extracted, including peaks with a significance level above 90%: ‘bi-decadal’ (16.52–20.48 years), ‘tri-decadal’ (23.27–32.0 years), and ‘multi-decadal’ (39.38 years and longer) variations (Fig. 1a). Minobe (1999, 2000) extracted ‘bi-decadal variation’ using a 10–30 year bandpass filter, while a more limited timescale (16.52–20.48 y) is extracted for the ‘bi-decadal variation’ in the present study. The ranges of the three frequency bands were determined according to the spectral troughs. Note that the following results are not sensitive to reasonable changes in the frequency band width extracted.

Fig. 1
figure 1

a Power spectrum density of the 298-year PDOI timeseries with 1-2-1 smoothing applied. The red, blue and green shaded boxes respectively denote the bandpass frequencies of the bi-decadal, tri-decadal and multi-decadal variations. The vertical lines in each box show 18.6 year (in red), 27.9 year (= 1.5 × 18.6 year in blue), 55.8 year (= 3 × 18.6 year in green). The red (blue) dashed line indicates 95% (90%) significance level determined by Monte Carlo method. b Reconstructed timeseries of PDOI during 1700–1997. 7-year-running mean (gray bars), 16.52–20.48 years bi-decadal (red curve), 23.27–32.0 years tri-decadal (blue) and 39.38 years-multi-decadal (green) band-passed components. The dots and the crosses denote the synchronizations of the phase reversals (zero-crossings) between the bi-decadal and the multi-decadal and between the bi-decadal and tri-decadal variations respectively. The blue- and pink-shades denote continuous synchronizations between the bi-decadal with the period of 18.6 year and the multi-decadal 3 × 18.6 year during 1770–1880 and 1920–1997 respectively with the gap during 1881–1919 (color figure online)

The band-passed timeseries of the three interdecadal bands show synchronized phase reversals as in Fig. 1b. Out of the three components, we first focus on the relationship between the bi-decadal (red line in Fig. 1b) and the multi-decadal (green line in Fig. 1b) variations. The negative to positive synchronized zero-crossings (phase reversals) were observed in around 1770, 1828, 1880, 1920, and 1975, while those with positive to negative were observed in around 1800, 1855, and 1945 (dots in Fig. 1b). This feature continues during 1780–1880 and 1920—with a skip during 1880–1920. During 1780–1880 and 1920—where the synchronous zero-crossings are evident, the multi-decadal variation retains the period of approximately 60 years (3 times that of the bi-decadal variation).

We then address the relationship between the bi-decadal and the tri-decadal (blue line in Fig. 1b) variations, where a similar feature can be seen in the zero-crossings. Since the tri-decadal variation retains a periodicity of approximately 28 years as opposed to that of bi-decadal variation (approximately 19 years, which is approximately 2/3 of 28 years), these synchronized zero-crossings can only be found every 50–60 years (~ 3 times 19 years) as negative to positive phase reversals (crosses in Fig. 1b). We call these notable synchronizations ‘phase-locked’. It is worth mentioning here that these three interdecadal (the band-passed bi-/tri-/multi-decadal) variations account for 81.7% of the total variance in the reconstructed PDOI with a running mean of 7 years.

Here we show the results of the significance test for phase-locks (see ‘Significance of phase-locks’ in the Data and methods section and “Appendix”). To quantitatively assess the ratio of phase locks as described in the Data and Methods section, we first define the possible maximum number of phase-locks as follows. For the tri-decadal variation, with a typical period of 28 years (approximately 1.5 times the period of the bi-decadal variation), the phase-locks can only occur once for every three periods of the bi-decadal variation. Considering that the bi-decadal and tri-decadal variations showed the first phase-lock around 1730 when they turned from negative to positive, the maximum number of phase-locks possible is 5 (approximately 1730, 1780, 1840, 1900, and 1950; crosses in Fig. 1b). Similar consideration was applied to the multi-decadal variation, which has a typical period of 50–60 years (approximately 3 times the period of the bi-decadal variation). The phase locking works differently in this case because phase-locks with the multi-decadal oscillation can also occur in positive to negative phase reversals, resulting in a maximum of 9 (dots in Fig. 1b). Consequently, we have 14 as a possible maximum number of total phase-locks. This process enables the definition of the ratio of the phase-locks.

Figure 2 shows that the probability of the observed phase-lock ratio is below 0.01 (0.05) for a tolerance of 1.0–2.0 (0.5–3.2) year. This suggests a remarkably significant relationship between the bi-decadal variation and the other two (tri- and multi-decadal) variations. Note that phase-locking is evident in 13 out of 14 (maximum possible) zero-crossings when the tolerance is sufficiently wide (~ 5 year), suggesting that a low phase-lock ratio with a strict tolerance (< 1 year) does not always imply the absence of a successive phase-locking temporal structure.

Fig. 2
figure 2

Diagram of statistical significance for the total number of synchronization between the bi-decadal and the tri-decadal variations and between the bi-decadal and the multi-decadal variations. The horizontal axis is the tolerance (in year), which is defined as the extent to which the phase reversals from the same polarity to the other are considered simultaneous. The bold red line indicates the ratio of number of the observed phase-locks to that of the possible maximum phase-locks. In the present case, the maximum is 14 (5 for the phase-lock with the tri-decadal variation, and 9 with the multi-decadal variation). The dashed lines show the significance level determined by Monte Carlo method, with 99%, 95%, 90% and 80% from the top. The synchronizations are counted except 30 years (former end) and 20 years (latter end) of the whole record due to uncertainty arising from window function applied to both ends

The synchronization between the bi-decadal and tri-/multi-decadal variations also reproduces the annual mean PDOI computed from the 119-year (1900–2018) monthly instrumental SST (HadlSST; Rayner et al. 2013), including the period (1998–2018) that was not covered by the reconstructed data. Figure 3a shows the sum of the three band-passed variations from the reconstructed data up until 1997 (the black curve in Fig. 3a), an extension from 1998 to 2018 (the green dashed curve in Fig. 3a), and the PDOI with 7-year running mean from HadlSST (the red curve in Fig. 3a). The correlation coefficient between the 7-year running mean from HadlSST and the sum of the three variations from 1900 to 2018 is still high at 0.85, and most importantly, the observed phase reversals in 2000 and 2010 were reproduced by the sum of the three interdecadal variations. In other words, the bi-decadal variation and the associated two variations are successful in predicting the two phase reversals that actually took place in the period that the reconstructed data did not cover.

Fig. 3
figure 3

a 7-year running mean PDO timeseries from instrumently observed Hadley Centre SST (HadlSST) during 1900–2018 (red curve). The black solid curve is the sum of the three (bi-decadal, tri-decadal, and multi-decadal) variations derived in the present study. The extension of the sum of the three variations after 1997 are drawn with the green dotted curve. b Band-passed bi-decadal (10–30 years: red curve), penta-decadal (30–80 years: blue) and the sum (10–80 years: black) timeseries of HadlSST PDO. c Same as b, but applied to the longer reconstructed PDO timeseries of 1700–1997. Black dots denote the continuously observed three synchronous zero-crossings in the twentieth century, and associated period is shaded in pink (color figure online)

We then examined the PDOI in relation to the hypotheses for climatic regime shifts given in the previous studies (Minobe 1999, 2000; Wilson et al. 2007). Although the phase-locking features are less pronounced as in Fig. 3b, the phase reversals of the 10–80 year band-passed time-series took place where the bi-decadal (10–30 year) and the penta-decadal (30–80 year) variations changed signs within a short period (< 5 year) in the 1920s, 1940s, and 1970s, which is consistent with Minobe (1999, 2000). On the other hand, applying the same bandpass filters to the reconstructed PDOI reveals that the phase-locking association that shows continuous zero-crossings between the 10–30 year and 30–80 year variations are not always observed before 1900, and are unique to the twentieth century (Fig. 3c). This is consistent with the results of Wilson et al. (2007).

Here we propose an alternative and comprehensive explanation for these two conflicting hypotheses by considering the three (bi-decadal, tri-decadal, and multi-decadal) variations as opposed to the two (10–30 years and 30–80 years) variations employed in the previous studies. The negative PDO (denoted by gray bars in Fig. 1b) regime is more conspicuous and the positive ones can hardly be recognized during the period from the late 18th to the late nineteenth century. This could be explained with the tri-decadal variation, which emphasizes the negative phases (around 1800–1830 and 1850–1880) and cancels out the positive phases (around 1770–1800 and 1830–1850) in the multi-decadal variation; simultaneous minima of the tri- and multi-decadal variations enhance the negative phases and coincident minima of the tri-decadal and maxima of the multi-decadal weaken the positive phase. In contrast, the opposite occurs in the phase relation between the tri- and the multi-decadal variations during the twentieth century, emphasizing the positive phases and cancelling the negative. The multi-decadal variation becomes more pronounced after the 1920s; the tri-decadal variation contributes less and thus does not bring the overall variation back to positive during 1940–1970, rendering both polarities more persistent and the regime shifts more remarkable. In short, the phase-locking of the bi-decadal variation with the tri-decadal variation as well as the multi-decadal variation with a gap during 1880–1920 can comprehensively explain the temporal structure of the regime shift before and after the 1900s.

4 Summary and discussion

Using 298-year long reconstructed PDOI, we found that the tri-decadal variation and the bi-decadal variation show phase-locks as well as those between the multi-decadal variation and the bi-decadal variation which was suggested in the previous works (Minobe 2000; Wilson et al. 2007). By taking into account the tri-decadal variation and its phase-locks with the bi-decadal variation, we proposed a viewpoint that comprehensively explain the continuously phase-locked temporal structure behind the climatic regime shifts, which was previously considered unique only after the 1900s (Wilson et al. 2007). We further verified that these three phase-locked variations can also be found in the more reliable 119-year instrumental PDOI.

We note that the spectral peaks that correspond with these variations have also been addressed in a number of previous studies concerning interdecadal variability, although their relevance with respect to the temporal structure resulting from synchronization among components has received little attention. Studies using other PDOI reconstructions have also found spectral peaks that may be equivalent to the bi-decadal variation, tri-decadal variation (Biondi et al. 2001; MacDonald and Case 2005; d’Arrigo and Wilson 2006), and multi-decadal variation (Biondi et al. 2001; MacDonald and Case 2005; d’Arrigo and Wilson 2006; Shen et al. 2006). Other than the PDO, the reconstructed surface air temperature of the Gulf of Alaska (Wilson et al. 2007; Wiles et al. 2014) and Nemuro, Japan (d’Arrigo et al. 2015), and the observed sea level on both a global scale (Jevrejeva et al. 2008; Chambers et al. 2012) and the local scale along the coast of Japan (Parker 2019) also demonstrate fluctuations that correspond with one or more of the three variations, supporting the reliability of the interdecadal variations observed in the present study.

The mechanism for the generation of these three variations is discussed next. First, we focus on the bi-decadal variation. A number of previous studies have highlighted mechanisms for the basin-scale variation in the North Pacific on this specific timescale (~ 20 years), with respect to both ‘internally sustained’ and ‘externally driven’ processes. Latif and Barnett (1994, 1996) were the first to provide the former perspective in the North Pacific. They argued that the positive and negative feedbacks in the local air–sea coupling and the oceanic gyre circulation sustains ~ 20 year oscillation in the North Pacific climate system. Although the mechanism for the negative feedback has been concluded unlikely by Schneider et al. (2002), many other studies on this topic have been conducted following their ‘local air–sea coupling’ framework. One well-discussed example is the work done by Zhang and Delworth (2015), which pointed out using numerical models that the self-sustained bi-decadal variation is created through the local air–sea coupling in the extratropical North Pacific climate system. The mechanism proposed in their study can be briefly addressed as follows. First, the subsurface temperature anomaly in the western subtropical/subpolar ocean reverses SST anomaly in the Kuroshio–Oyashio Extension region after transit time of 5–10 years. Then, this reversed SST anomaly induces atmospheric response that create opposite subsurface anomaly to the previous one through wind forcing over the central North Pacific and the oceanic Rossby wave adjustment. These processes complete one cycle every ~ 20 years. However, the physical mechanism of how SST anomaly affects the overlying atmosphere is still under debate. Also, Liu (2012) extensively reviewed the mechanisms of 10–20 year variability in the North Pacific in terms of both extratropical-origin (= internally sustained) and tropical-origin (= externally driven) scenario.

Furthermore, we note that there are suggestive previous publications on the effect of tidal force fluctuation on the bi-decadal variation. A previous study elucidated, using the same reconstructed PDOI as that applied in the present study, that the significant spectrum peak of the bi-decadal variation has the period of 18.6 years, suggesting that the 18.6-year lunar nodal cycle is somehow related to the bi-decadal variation in the PDO (Yasuda 2009). This phenomenon was first considered in the field of physical oceanography by Loder and Garrett (1978), where the effects that the 18.6-year cyclic fluctuation in tidal forcing have on the local SST were quantitatively demonstrated. In addition to the observational evidences for the 18.6-year tidal cycle (Cook et al. 1997; Yasuda et al. 2006; Osafune and Yasuda 2006, 2010, 2012, 2013; Wilson et al. 2007; McKinnell and Crawford 2007; Yasuda 2018), several numerical studies with air–sea coupled climate models and ocean models have also identified possible contributions to this phenomenon by inflicting periodically (18.6 years) fluctuating vertical turbulent diffusivity (e.g. Hasumi et al. 2008; Tanaka et al. 2012; Osafune et al. 2020).

The extent to which this external forcing contributes to basin-wide oceanic phenomena such as the PDO is still not well understood, provoking arguments whether it is even of any significance. Still, considering that the bi-decadal peak is sharp and close to 18.6 years as found by Yasuda (2009), we can reasonably assume that our bi-decadal variation can be relevant to the 18.6-year lunar nodal oscillation. This could be understood as a periodic pace-maker that entrains the self-sustained bi-decadal variation in the North Pacific (e.g. Zhang and Delworth 2015; Qiu et al. 2017) to its period (= 18.6 years in this case) in a nonlinear climate system.

We then focus on the other two interdecadal variations, which may be influenced significantly by the bi-decadal variation. The tri-decadal variation has a typical period of 28 years (or approximately 18.6 years times 1.5), while multi-decadal variation covers approximately 50–60 years with the period 55.8 years (= 18.6 years times 3) in between, the latter variation of which is suggested to originate in the air–sea coupling and the oceanic Rossby wave propagation in the extratropical North Pacific (Zhong and Liu 2009). Here, in addition to the proposed mechanisms of interdecadal variability that corresponds with these timescales, we introduce the phenomenon called ‘subharmonic resonance’ to suggest possible contributions of nonlinearity to the climate system to reflect the ‘phase-locked’ temporal structure. Subharmonic resonance refers to the excitement of components with integral multiple periods of that of the external forcing \(nT\) Although this idea has not often been incorporated into studies investigating oceanic decadal variability, Minobe and Jin (2004) suggested the importance of subharmonic resonance using a simplified model of tropical oceanic features. Periodic external forcing and periodic modulation were applied to the nonlinear delayed oscillator equation developed in a study of the El Niño and Southern Oscillation (ENSO) by Battisti and Hirst (1989), and as a result, demonstrated that odd-multiple resonances and even-multiple resonances are excited respectively. In terms of decadal variability in the North Pacific regarding this phenomenon, there is another important suggestion by Skákala and Bruun (2018), which used simple toy model of tropical-extratropical climatic interaction similar to that used in Minobe and Jin (2004). They showed that the bi-decadal variation discussed in Zhang and Delworth (2015) could excite variations with other periodicities in the North Pacific climate system, including one with three times longer (= 60 years) period than the original variation (= 20 years), which exhibits phase-locked temporal structure.

Taking these into account, we have good reasons to suspect that the similar mechanisms give rise to the tri-decadal and multi-decadal variations out of the bi-decadal variation. We therefore consider that the tri- and multi-decadal variations, as well as the bi-decadal variation, may be associated with the 18.6-year lunar nodal cycle.

A careful inspection of these variations, however, reveals that there are likely other drivers that contribute to variations with these timescales than the above-mentioned nonlinear process. The tri-decadal variation is typically 1.5 times longer (28 years) than that of the bi-decadal variation (19 years). Considering that subharmonic resonance only refers to excitement over periods which are multiples of full integers, we cannot exactly relate the tri-decadal variation with the bi-decadal variation via this phenomenon. For instance, other nonlinear processes may generate half-period components from multiples of the bi-decadal variation to cause the tri-decadal variation; however, verifying the occurrence of such processes requires more investigation into the generation mechanisms. Furthermore, subharmonic resonance does not fully explain the multi-decadal variation; it changed its polarity much faster during the period 1880–1920. It changed from negative to positive in the 1880s, and rapidly turned negative in the 1890s only to turn positive again in the 1920s, which is inconsistent with the temporal structure observed elsewhere (successive zero-crossings approximately every 50–60 years). This discontinuous feature in the multi-decadal variation cannot be expressed with subharmonic resonance (see “Appendix” for further discussion on this discontinuous feature and its possible effects on the power spectrum). Taking the above discussion into consideration, further validation is required to verify that the tri- and multi-decadal variations are excited by the bi-decadal variation through nonlinear processes.

In terms of the multi-decadal variation, although the reason for the discontinuity in the synchronization during 1880–1920 remains unclear, several studies suggest that climatic changes can be associated with volcanic eruptions, which generate short-term (< 10 years) cooling trends (Rampino et al. 1979; Self et al. 1981; Overpeck et al. 1997). Thus, the volcanic eruptions (for example, the eruption of Krakatau, Indonesia in 1883) might be considered a phenomenon that has driven the multi-decadal variation in a negative sense, eventually resulting in the rapid phase reversal from positive to negative in the 1910s. However, confirming these possibilities requires physical and quantitative evaluation, which is beyond the scope of the present study.