A Mechanism for Pacific Interdecadal Resonances

Pacific interdecadal variability (PIV) is an important large‐scale climate phenomenon. There is growing evidence that PIV contains three spectral resonances: a decadal (13 ± 1‐year) spectral peak, bidecadal (20 ± 5‐year) resonance, and a pentadecadal (60 ± 10‐year) resonance. Although much has been clarified about mechanisms behind PIV, there are still many open questions about the origin of these resonant modes (especially the pentadecadal mode). We describe dynamics in the Pacific basin by a toy (delayed oscillator) model that sheds light on the nature of these resonant peaks. The model suggests that unlike the bidecadal resonance, which results from local atmosphere‐ocean coupling in the extratropics, the pentadecadal and possibly also decadal resonances result from atmospheric and oceanic teleconnections between the extratropics and tropics. We show that a tiny coupling between extratropics and tropics through the ocean tunnel is sufficient to trigger the pentadecadal oscillation in the Pacific basin. Our model also explains (a) the observed three‐period locking between the bidecadal and pentadecadal modes and (b) the synchronization of anomalies in the central eastern tropics and central North Pacific with the opposite relative sign. We conclude that the role of oceanic teleconnections is probably underestimated in the current literature on PIV.


Introduction
Pacific interdecadal variability (PIV) is an important part of the Earth climate system. The uncertainties associated with understanding PIV are a serious limitation to climate prediction (Liu, 2012). In its dominant mode, the Pacific Decadal Oscillation (PDO) is defined as the leading empirical orthogonal function of monthly averaged sea surface temperature (SST) anomalies (SSTA) in the North Pacific (NP; Mantua et al., 1997). The PDO is broadly understood to result from a combination of different processes, such as stochastic wind forcing, tropics -extratropics atmospheric interaction, and long ocean memory (e.g., Newman et al., 2016;Wang et al., 2017). Even though the PDO is a complicated signal, three specific PDO (quasi)periodic modes have been identified with 13 ± 1-, 20 ± 5-, and 60 ± 10-year periods (Baumgartner et al., 1992;Bruun et al., 2017;Burroughs, 2007;Chao et al., 2000;Mantua & Hare, 2002;Mantua et al., 1997;MacDonald & Case, 2005;McCabe-Glynn et al., 2013;Minobe, 1997Minobe, , 1999Li et al., 2013;Liu, 2012;Lüdecke et al., 2013;Taguchi et al., 2007Taguchi et al., , 2012Wang et al., 2017;Yasuda, 2009). Simple structures such as (quasi)periodic oscillations are important, since they can serve as a window to help understand potentially significant physical drivers behind the PDO. This is because sustained long-term (quasi)periodic behaviour is either a result of external oscillatory forcing or an internal mode of the coupled ocean-atmosphere system in the Pacific (Minobe, 1997). The internal mode can be either an Earth system response to stochastic external forcing or a fully sustained deterministic process.
The existence of PIV long-term resonances can be demonstrated using Dominant Frequency State Analysis (DFSA; Bruun et al., 2017;Tong, 1990;Young, 2012;Young et al., 2001), which was recently developed and applied in Bruun et al. (2017) to instrumental data, such as North Pacific Index Anomaly (NPIA) and Southern Oscillation Index (SOI), as well as paleoreconstruction of sea surface temperature anomaly (SSTA) based on tree-ring data (calibrated with the Kaplan SSTA record; see Li et al., 2013). DFSA is a fully efficient eigenvalue identification approach that incorporates red noise in the estimation. It is especially suited to the task of accurately assessing the presence of small-amplitude climatic mode coefficients. An example of PIV modes identification via dominant frequency analysis applied to NPIA data set is shown in Figure 1. In the light of advanced signal analysis tools such as DFSA one can better understand why the PDO resonances sometimes remain hidden. For example, in Newman et al. (2016) the tropical and extratropical coupling is accommodated by an empirically motivated autoregressive (AR) transfer function. This is then spatially generalized using a vector AR process where low-frequency NP modes are distinguished from the higher-frequency and mostly NPIA is deseasonalized North Pacific Index (NPI, the area-weighted sea level pressure over the region 30-65 ∘ N, 160 ∘ E-140 ∘ W) calculated by subtracting the dominant (1-and 0.5-year) regular seasonal frequency states from the NPI signal. The NPIA data are annually averaged. (a) The dominant 5.8-, 14-, and 61-year frequency states (black arrows) with approximate 95% confidence intervals. The red arrow is the combined effect of the frequency states scaled to a value of 0.5. (b) The Eigenvalue magnitude of each frequency state with 95% confidence interval. The 61-year frequency state is the largest. (c) The total frequency states (red line) in the time domain. The green line is the frequency states and (the nonsignificant) trend, and blue is the DFSA model fit also including AR(1) red noise (reproduced from Bruun et al., 2017). DFSA = Dominant Frequency State Analysis; NPIA = North Pacific Index anomaly.
tropical El Niño-Southern Oscillation (ENSO) signals in the east. In that specific approach the PDO and ENSO coupling is structured empirically so that the PDO transfer function with the ENSO is a 2-year reddened memory. The limitation of that empirical AR approach is that the same transfer function then also reddens the uncorrelated noise process, inducing an autocorrelated noise structure. This may explain the debate around the existence or not of 20-and 70-year interdecadal peaks in Newman et al. (2016): the presence of autocorrelated noise can result in more uncertainty due to a less efficient and biased estimation of the actual physical transfer structure relationship. Also, the empirical form of that transfer function does not readily generalize to enable efficient estimation of the structure and form of low-frequency modes.
The Pacific interdecadal modes are an area of an intense mutual debate. In the next sections we will, for each mode individually, summarize the evidence found both observationally and in models, as well as what has been suggested in the literature about the underlying physics of PIV.
A widely used explanation for bidecadal PDO is that it originates as an ocean-atmosphere coupled mode in the NP (L. Zhang & Delworth, 2015). From the synthesis of previous knowledge Jin, 1997;Latif & Barnett, 1994, 1996Schneider et al., 2002;Zhong & Liu, 2009) one gets a better idea about the mechanism behind bidecadal PDO (L. Zhang & Delworth, 2015): The positive (negative) subsurface anomalies in western subtropics and western subpolar regions upwell to the surface triggering anomalous cyclonic (anticyclonic) winds in the western NP. The cyclonic (anticyclonic) wind anomaly cools (warms) the central NP through anomalous Ekman transport and subsequently through northward shift of the subtropic and subpolar gyres. The anomalies from subtropic and subpolar regions are advected by Kuroshio and Oyashio currents to the KOE region. The KOE anomaly then strengthens the cyclonic (anticyclonic) winds and is amplified by a westward propagating Rossby wave from the central NP (which takes 3-5 years to propagate). The anticyclonic winds lead to positive subsurface anomalies that propagate as westward Rossby waves in the subpolar and subtropic region, upwell in the subpolar and subtropic western NP, and eventually reverse the bidecadal PDO cycle (in ∼ 10 years time). The warm cycle is therefore reversed by the cold subtropic and subpolar anomalies created by a PDO cold phase a decade ago (L. Zhang & Delworth, 2015). The anomalies in the eastern NP tongue are driven by anomalous winds originating from the western half of the basin.
Furthermore an intriguing analysis of the NPI record by Minobe (1999) suggested a possible relationship between the two quasiperiodic bidecadal and pentadecadal PIV modes as follows: (a) Within a centennial time series the two modes seem to be synchronized with a relative period of three. (b) These modes follow from different mechanisms but interact with each other. (c) Within most of the twentieth century the dominant phase locking appears to correspond to a mutual phase reversal with the same relative sign (Minobe, 1999). Minobe and Jin (2004) used a prototype equation to explain how such period 3 phase-locked subharmonic resonance appears from a delayed nonlinear feedback to a periodic (i.e., 20-year) external forcing. However, the study of Minobe and Jin (2004) leaves open the question of what the exact physical mechanism is providing such feedback and what implication such mechanism has for the dynamics in the Pacific basin.
To the knowledge of authors the only suggested Earth system mechanism explaining the pentadecadal mode in the literature is via the long time memory of the subpolar NP Rossby waves (Zhong & Liu, 2009). This explanation faces several issues: (a) The modeling studies show that the subpolar Rossby wave mechanism leads to 40-to 45-year periodicities that are somewhat shorter than the 60 ± 10-year mode (Kwon & Deser, 2007;Wu et al., 2003;Zhong & Liu, 2009). (b) The SOI index, as well as tree-ring data analyses (Bruun et al., 2017;Li et al., 2013), suggests that also the pentadecadal mode is a Pacific basin-wide phenomenon. The subpolar NP and SP (NP-SP) regions have two distinct open ocean and land boundary conditions, so zonal Rossby wave propagation times are likely to differ. As a consequence we would expect very Figure 2. This figure shows schematically the tropics-extratropics exchange presented in the toy model. The scheme is only for NP, but the model assumes that similar relationships exist for the SP. The three regions mentioned in the text (KOE, central NP, and central eastern tropics) are clearly marked. The atmospheric (red) and oceanic (blue) mechanisms of the exchange are marked by the arrows. The capital letters mark the extratropics-tropics exchange pathways as follows: A = the NPA/coastal Kelvin wave-induced wind stress curl anomaly; B = the wind stress curl anomaly in the tropics triggered from the extratropics; C = the long ocean tunnel from the extratropics to the tropics (e.g., Gu & Philander, 1997); D = the coastally trapped Kelvin waves triggered from the tropics; and E = the westward propagating oceanic Rossby waves in the extratropics. NP = North Pacific; SP = South Pacific; KOE = Kuroshio-Oyashio Extension Region. different characteristics of multidecadal modes in NP and SP. (c) The subpolar NP Rossby wave mechanism does not seem to produce or explain the factor of 3 connectionbetween the bidecadal and pentadecadal modes observed in Minobe (1999). Indeed, none of the modeling studies where the pentadecadal mode occurred (Kwon & Deser, 2007;Wu et al., 2003;Zhong & Liu, 2009) successfully reproduced the period 3 mode locking observed in Minobe (1999).

Decadal Mode
In addition to those two resonances, Li et al. (2013) and Bruun et al. (2017) recently identified a third, 13 ± 1-year (decadal) statistically significant spectral peak in the Pacific SST derived from a 700-year tree ring data record. The decadal peak has also been observed in the SOI and NPIA (Bruun et al., 2017). Curiously, it has been shown by Li et al. (2013) that this periodicity cannot be convincingly linked to solar cycles (that happen at an ≈ 11-year scale) and it is most likely an independent Earth system mode (Li et al., 2013). Such decadal oscillations can be potentially explained by off-equatorial wave dynamics (i.e., White et al., 2003). Note that Bruun et al. (2017) identified in tree ring data also 180-year centennial periodicity, this one likely associated with the known 180-to 200-year de Vries solar cycle (Lüdecke et al., 2015;Wagner et al., 2001).

The Tropics-Extratropics Exchange and Our Proposal
It has been suggested (mostly analyzing tree-ring SST records) that PDO is a basin-wide phenomenon with roughly symmetric patterns across the equator (Evans et al., 2001;Garreaud & Battisti, 1999;Liu & Alexander, 2007). This supports the idea that there is a broad interaction in the whole Pacific basin. Furthermore, both observational data and models (Liu, 2012;Newman et al., 2016;Pierce et al., 2000;Tatebe et al., 2013;Yeh et al., 2015;L. Zhang & Delworth, 2015) show that PDO oscillation in the KOE and central NP is synchronized with an ENSO pattern in the central eastern tropics that has opposite sign to the KOE and central NP anomaly. Similarly, Bruun et al. (2017) found connection between SOI and NPIA data in their DFSA analysis of the low-frequency variability of ENSO. This clearly implies a form of exchange between extratropics and tropics, but it is unclear whether the origin of the PIV is in the extratropics or tropics. A number of studies Jin, 1997;Latif & Barnett, 1994, 1996Schneider et al., 2002;L. Zhang & Delworth, 2015) concluded that (at least the bidecedal) PDO originates in the extratropics and subsequently drives variability (with the opposite phase) in the central eastern tropics through anomalous wind stress curl Pierce et al., 2000). It has been also suggested (Evans et al., 2001) that the atmospheric bridge works the other way round and PDO has its origin in the tropics. Some (surgical modeling) studies (Kwon & Deser, 2007;Liu, 2002;Liu & Alexander, 2007;Wu et al., 2003) even concluded that the modes in the extratropics and tropics have both a mutually independent local origin with additional ability to reinforce each other through the atmosphere. In any case the negative correlation between the central eastern tropics and central western NP most likely happens as a result of anomalous wind stress curl forced from NP (e.g., , Pierce et al., 2000, and/or via tropical anomalies affecting the Aleutian low through Pacific-North American (PNA) teleconnection (e.g., Alexander, 1992;Liu, 2012;Liu & Alexander, 2007;Zhou et al., 2007).
In the present paper we synthesize ideas from the last couple of decades and present a simple toy model that describes the mechanism for low-frequency (period > 8-10 years) dynamics in the Pacific basin: (a) The model assumes that the local extratropics dynamics are described by the bidecadal oscillations through the atmosphere-ocean coupled mechanism summarized by L. Zhang and Delworth (2015). The mechanism has been proposed for the NP, but based on observed PDO symmetric patterns across the equator (Evans et al., 2001;Garreaud & Battisti, 1999;Liu & Alexander, 2007), we assume that a similar mechanism drives oscillations in the SP. (b) A complex mechanism for the tropics-extratropics interaction (for NP, see Figure 2) extends the extratropics anomaly to the tropics and provides feedback of tropics on extratropics. The tropics-extratropics interaction is mediated through a fast atmospheric bridge (Alexander, 1992;Kwon et al., 2011;Liu, 2012;Liu & Alexander, 2007;Newman et al., 2016;Pierce et al., 2000;Wang et al., 2012;Yeh et al., 2015;Zhou et al., 2007), slow ocean tunnel (Gu & Philander, 1997;Liu, 2012;Liu & Alexander, 2007; Journal of Geophysical Research: Oceans 10.1029/2018JC013752 Tatebe et al., 2013;R.-H. Zhang et al., 1998), and coastally trapped Kelvin waves (Newman et al., 2016;Liu, 2012;Liu & Alexander, 2007;Power & Colman, 2006).
We show that combining the mechanism for tropics-extratropics interaction with the local bidecadal mode in the extratropics produces a model that (a) adds to the bidecadal mode two more resonances: (the observed) 12-and 60-year peaks, (b) automatically produces period 3 locking between the bidecadal and pentadecadal modes as observed in Minobe (1999), and (c) synchronizes the anomalies in the central eastern tropics and central NP (or KOE) region with the opposite phase. The toy model also provides a natural explanation for the PDO as a basin-wide phenomenon and suggests that unlike the bidecadal resonance, the pentadecadal and possibly also decadal resonances originate from the interaction between extratropics and tropics. The model has therefore a remarkable ability to combine and explain all the main features observed in the periodic resonances of the PDO.

The Delayed Oscillator Model for the Pacific Basin
Let us start with the following equation describing the dynamics of the eastern tropical SST anomaly, T tr : (1) T ex is the SST anomaly in the extratropics, K O e ≥ 0 is a coupling constant for the impact of the extratropics on the tropics through the slow ocean tunnel (Gu & Philander, 1997;Liu, 2012;Liu & Alexander, 2007;Tatebe et al., 2013;R.-H. Zhang et al., 1998), d is the time it takes anomalies to propagate from the extratropics to the tropics, and K A e ≥ 0 is the coupling constant for the impact of the extratropics on the tropics through the fast atmospheric bridge Kwon et al., 2011;Newman et al., 2016;Pierce et al., 2000;Wang et al., 2012). The ENSO term contains the usual ENSO delayed oscillator dynamics (Battisti & Hirst, 1989;Bruun et al., 2017;Cane et al., 1990;Münnich et al., 1991;Neelin et al., 1998;Suarez & Schopf, 1988;Tziperman et al., 1994;Wang, 2001;Wang et al., 2017;Wang & Weisberg, 1998;White et al., 2003). The terms that couple the tropics to the extratropics can be described as follows: a. The ocean tropics to extratropics coupling term (K O e , Figure 2, C): Oceanic advection of anomalies from the extratropics to the tropics was first proposed by Gu and Philander (1997). The extratropics anomalies subduct in the NP around 140 ∘ W and 30 ∘ N and in the SP around 90 ∘ W and 30 ∘ S (Goodman et al., 2005;Gu & Philander, 1997;Liu et al., 1994) and reappear in the eastern tropics with median time delay of d = 15 years (Goodman et al., 2005;Gu & Philander, 1997). The idea of Gu and Philander (1997) obtained support from the observational data (R.-H. Zhang et al., 1998), but inconsistency with other observations was also noted (Lu & McCreary Jr. 1995;Sasaki et al., 2010;Schneider et al., 1999). There is now however growing evidence from the Southern Hemisphere that gyres propagate anomalies from the SP to the tropics and significantly affect ENSO dynamics (Giese et al., 2002;Luo & Yamagata, 2001;Sullivan et al., 2016;Tatebe et al., 2013). We include into our model the oceanic interaction proposed by Gu and Philander (1997), but to maintain reasonable consistency between our model and some of the observations (from Lu & McCreary Jr., 1995;Sasaki et al., 2010;Schneider et al., 1999), we will keep the oceanic coupling much weaker than the coupling through the atmosphere. The subduction region lies in the NP very close to the boundary between the eastern part of the central NP PDO region and the eastern NP PDO region (L. Zhang & Delworth, 2015). Based on closer inspection of the subduction region (Goodman et al., 2005) combined with the fact that central NP anomalies are much larger in amplitude than the eastern NP anomalies, we propose that the long ocean tunnel is dominated by anomalies originating in the eastern part of central NP. The T ex variable therefore describes the SST averaged through a large section of the central NP ( Figure 2). As indicated before, we assume generic latitudinal symmetry across the equator and extend this view to the SP. This will be also true about the other drivers of tropics-extratropics interaction. The mechanism will be described for NP (where it is well understood), but it will be assumed that similar interactions operate in SP. b. The atmospheric tropics to extratropics coupling term (K A e , Figure 2, B): The cold (warm) anomaly in KOE and central NP is associated with strong cyclonic (anticyclonic) winds in the NP. We follow the studies of  and Pierce et al. (2000) where it has been found that winds associated with the NP SST anomalies impact wind stress in the tropics and the slope of the thermocline. The cyclonic (anticyclonic) winds reduce (increase) the slope of the tropical thermocline introducing warm (cold) anomalies in the east and cold (warm) anomalies in the west. This can partly account for the negative correlation between the eastern tropical SST anomaly and the KOE or central NP SST anomaly (for the atmospheric teleconnection and the negative correlation see also Kwon et al., 2011;Wang et al., 2012). Since the cyclonic (anticyclonic) winds are strongly associated with cold (warm) SST in central NP Kwon et al., 2011;Newman et al., 2016;Pierce et al., 2000;Wang et al., 2012;L. Zhang & Delworth, 2015), it is reasonable to use the T ex value as a proxy for the wind stress impact on eastern tropics SST T tr .

Low-Frequency Component of Equation (1)
The ENSO term in equation (1) dominates the high-frequency interannual dynamics of T tr , whereas the extratropics to tropics coupling is expected to drive the lower-frequency interdecadal dynamics of T tr . We assume that the dynamics can be decoupled on the interannual and interdecadal scales. This essentially means two things: a. For the interdecadal (low-frequency) dynamics the ENSO term can be largely neglected. ENSO is described by a nonlinear oscillator coming from atmosphere-ocean interaction within tropics and its delayed feedback on T tr through a combination of Kelvin and Rossby waves (Battisti & Hirst, 1989;Bjerknes, 1969;Bruun et al., 2017;Cane et al., 1990;Cane & Zebiak, 1985;Jin et al., 1994;Münnich et al., 1991;Neelin et al., 1998;Suarez & Schopf, 1988;Tziperman et al., 1994Tziperman et al., , 1995Wang, 2001;Wang et al., 2017;Wang & Weisberg, 1998;White et al., 2003;Zebiak & Cane, 1987). The anomaly T tr is in the ENSO delayed oscillator models forced by the delayed feedback and damped by evaporation and heat convection processes (e.g., Battisti & Hirst, 1989;Cane et al., 1990;Münnich et al., 1991;Suarez & Schopf, 1988;Tziperman et al., 1994). These damping processes are usually described by a cubic term ∼ T 3 tr . We neglect the intra-annual oscillatory ENSO dynamics on larger (interdecadal) time scales except for the damping processes that regulate the size of the T tr anomaly. We therefore replace the ENSO term by ENSO → − T 3 tr . b. The T ex dynamics are approximated by its low-frequency component T ex →T ex (the tilde symbol will further stand for the low-frequency component of SST).
We assume thatT ex dynamics is produced by a local mechanism in extratropics described in L. Zhang and Delworth (2015) with a dominant 20-year oscillation mode. In addition,T ex is forced from the tropics mostly via atmospheric teleconnection. The impact of the tropics on the extratropics is typically seen through ENSO years as a long SST anomaly extending along the western coastline of the Americas toward high latitudes (i.e., Jong et al., 2016). There are two ways how the tropical anomalies propagate: atmospheric bridge (along the North American coastline, mostly the PNA teleconnection (Alexander, 1992;Liu, 2012;Liu & Alexander, 2007) and eastern-boundary coastally trapped Kelvin waves (Figure 2, D; Liu, 2012;Liu & Alexander, 2007;Power & Colman, 2006;Tatebe et al., 2013). Both these mechanisms propagate tropical anomalies along the eastern NP PDO tongue very fast (on the scale of months). Once the anomalies reached the higher latitudes they produce anomalous wind stress curl (Figure 2, A) triggering anomaly with the opposite sign in the central NP and SP. In the NP the mechanism is described through the PNA teleconnection impact on the Aleutian low affecting the western central part of NP (Liu, 2012;Liu & Alexander, 2007). If positive (negative) T tr anomaly propagates to the North American coastline, it supports cyclonic (anticyclonic) winds in the NP that trigger negative (positive) central NP (T ex ) anomaly through anomalous Ekman transport and southward (northward) shift of the subpolar and subtropic gyre boundary (as a part of mechanism from L. Zhang & Delworth, 2015).
Ideally one would represent the bidecadal mode dynamics in the extratropics through a complex set of equations that involve interaction between multiple regions in higher latitudes and multiple variables. However, by doing this one loses much of the insight provided by a toy model. To keep the model simple, we approximate theT ex dynamics by an algebraic (rather than differential) equation: The first term in equation (2) is the 20-year oscillation mode (produced by the dynamics within the extratropics) with amplitude T 0 ( m is the frequency of the 20 year mode). The second term is the imprint ofT tr onT ex , which is the opposite-phase oscillation pattern inT ex created through the atmospheric teleconnection by an oscillation inT ex . The K A t coupling constant describes the proportion ofT tr that is projected (with an inverse sign) intoT ex . This leads to two coupled equations (one differential and one algebraic) forT tr ,T ex : One can generalize equation (4) by including also the oceanic Rossby wave feedback ofT tr onT ex . After the T tr anomaly propagates to the extratropics (via either atmospheric bridge or coastally trapped Kelvin waves), it impacts on T ex not only through the atmosphere (Figure 2, A) but potentially also through westward propagating Rossby waves (Figure 2, E; see Power & Colman, 2006;Tatebe et al., 2013). The internal Rossby waves propagate at various latitudes with different phase speed (depending on the latitude), but effectively cross the basin to the upwelling regions in the subtropic and subpolar western NP in a decade. Then the reddened (low pass filtered) ENSO anomalies (Power & Colman, 2006;Tatebe et al., 2013) are entrained on the surface as a part of standard PDO mechanism (L. Zhang & Delworth, 2015) eventually impacting on the central NP anomaly through anomalous cyclonic (anticyclonic) winds. We estimate the effective time between the initial appearance of the anomaly in the tropics and its integrated impact (via Rossby wave pathway) on central NP to be between 10 and 15 years. The most consistent way to include the complex impact of Rossby waves at different latitudes, while respecting the internal dynamical interactions in extratropics, is to modulate the amplitude of the 20-year extratropical mode by allowing it to be a function ofT tr (t − d ′ ). The parameter d ′ (in the 10 to 15-year range) is the effective feedback time that integrates the complex processes originating from multiple Rossby wave pathways at different latitudes. If we linearly expand the amplitude as a function ofT tr (t − d ′ ), the generalized equation (4) becomes The sgn is a signature function (it returns ±1 values depending on the signature of its argument). It must be included in equation (5) to provide the correct amplitude modulation: the warm (cold) SST anomaly that arrived from the tropics meets with SST anomaly in the extratropics and warms (cools) it. If the SST anomaly in extratropics is warm (cold), the tropical anomaly further increases the amplitude of the extratropical processes, whereas if the SST anomaly in the extratropics is cold (warm) the warming suppresses the amplitude of the 20-year mode. The K O t is the contribution to the 20-year T ex mode amplitude due to the Rossby wave signal that propagates the T tr anomaly (K O t = 0 reduces equation (5) to equation (4). The model described by equations (3) and (5) contains four coupling parameters: a. K A e (atmospheric impact of extratropics on tropics): Consider an extremely simplistic model without feedbacks or internal dynamics both in the tropics and extratropics. The value of the K A e constant (and similarly K O e ) in such model corresponds to an inverse time that it takes the extratropics T ex anomaly to trigger tropics T tr anomaly with the same size than T ex . The realistic range of values for this coupling constant is in an interval around 1 year −1 (e.g., corresponds to an inverse period between few months and a couple of years). For example, Gu and Philander (1997) used for their simulation K A e = 1 year −1 . b. K O e (oceanic impact of the tropics on the extratropics): To be consistent with the observational studies from Lu and McCreary Jr. (1995), Sasaki et al. (2010) and Schneider et al. (1999), we consider the realistic range of values for this coupling constant to be at least an order of magnitude smaller than K A t . This puts the K O e value on the scale of inverse periods from the order of years to decades. The opposite was done in Gu and Philander (1997), where K O e was taken to be 5 times larger than the atmospheric coupling. c. K A t (mostly atmosphere driven impact of the tropics on the extratropics): This coupling describes the size of T tr projection intoT ex through atmospheric forcing. For example, K A t = 0.4 means that the oscillation pattern inT tr creates an oscillation pattern inT ex with 40% ofT tr amplitude. We vary the value of K A t within the range between 0 and 1. d. K O t (mostly ocean-driven impact of the tropics on the extratropics): This is the coupling parameter that appears in the extended version of the model (equation (5)). We consider the realistic contribution from this term to be somewhere between 0 and T 0 ∕4, with T 0 being the amplitude of the 20-year extratropical mode.
The last parameters that need to be specified are the amplitude of the bidecadal oscillation in central NP and SP, T 0 (equation (4)), and the damping parameter (equation (3)). Based on the results of L. Zhang and Delworth (2015), we estimate the value of T 0 to be around 0.2 ∘ C. The damping constant value was set consistently with Gu and Philander (1997) to be = 1 C −2 year −1 .

Results and Discussion
The coupling betweenT ex andT tr produces similar spectra for both tropical and extratropical SST (i.e., Figure 3). This implies that the interdecadal spectra are basin wide as suggested by multiple observational studies The spectral peaks are shown on a log scale as a function of period (in years). The dashed lines mark the odd multiples of the lowest frequency mode (these are modes corresponding to the sequence of periods given as 60/1, 60/3, 60/5, 60/7 ... years). The periodogram shows that all the resonances are an odd multiple of the lowest (1/60 year −1 ) frequency. Since the y axis is on a log scale, the time evolution of bothT tr andT ex is dominated by the 20-and 60-year oscillations. The spectra were obtained from the last 5,000 years of a 20,000-year simulation (the first 15,000 years was used to remove the transient trajectories). (Bruun et al., 2017;Evans et al., 2001;Garreaud & Battisti, 1999;Li et al., 2013;Liu & Alexander, 2007). Within the realistic part of the parameter space we observed two separate regimes: (i) where the SST dynamics has a dominant 20-year oscillation mode and (ii) where the dynamics is a combination of 20-and 60-year oscillation modes. In each case the dynamics has a periodic attractor (Hilborn, 2000) that forms a sequence of resonances defined as an odd integer multiple of the lowest frequency (1∕60 year −1 ). However, most of these resonances are very small. Figure 3 shows this characteristics for a specific combination of parameter values.
The parameter space for the two regimes was explored for the simplified model with K O t = 0 (no Rossby wave feedback from tropical SST anomaly), while keeping the T 0 and K A t parameters fixed at T 0 = 0.2 ∘ C and K A t = 1 year −1 . The results are shown in Figure 4. In particular, we find the following: Although the oceanic tunnel from the extratropics to tropics is essential to produce the 60-year resonance, the coupling through the ocean tunnel can be very small compared to the coupling through the atmosphere. In fact, the 60-year peak is produced when the K O e coupling has a tiny value (K O e = 1∕25 year −1 ) and is 50 times smaller than the coupling through the atmosphere (K A e = 1∕6 month −1 , see Figure 4). This is a very important finding as it implies that even if the coupling through ocean is very small (as concluded by some studies, e.g., Lu & McCreary Jr., 1995;Sasaki et al., 2010;Schneider et al., 1999), the ocean tunnel appears to play an essential role in the generating mechanism for the pentadecadal peak. In fact, perhaps surprisingly, Figure 4 indicates that if the ocean coupling becomes of comparable size to the coupling through the atmosphere, the 60-year resonance disappears.
We also find that if the atmospheric forcing from the tropics on the extratropics (K A t ) is small there is not enough feedback to produce the 60-year resonance (only the 20-year forced oscillation is present). In general, when the nondimensional coupling K A t approaches ≈ 0.4 (in the specific case from Figure 4 it is 0.45) the feedback kicks in and the two-mode spectra (20 and 60 years) are produced. The K A t = 0.4 value means that 40% ofT tr size is projected via atmospheric teleconnection intoT ex . The two-peak spectrum disappears once the coupling becomes too large (the cases in Figure 4 show the two-peak spectra for the interval of K A t between 0.45 and 0.68). In some small subregions of the parameter space one can have spectra that have other dominant resonances than the 20-and 60-year peaks (see Figure 4b). For Figure 4b and K O e = 1∕3 year −1 the lowest-frequency mode has a 40-year period. The spectra are generated as an integer (both even and odd) multiple of the lowest frequency mode (this case is shown in supporting information S1).
Most of the simulations were performed without the Rossby wave oceanic feedback (K O t = 0); however, it was observed that varying the K O t parameter between 0 and T 0 ∕4 = 0.05 ∘ C did not have substantial impact on the spectra. The same was true when we varied the feedback time d ′ (see equation (5) in the expected 10-to 15-year range. Similarly, the two-peak (20-and 60-year) spectra are always obtained within the 13-to 20-year range of d values (anomaly delay time through the ocean tunnel). This means that the spectra are robust over the delay time parameter range. Figure 5 shows the time series from the last 200 years of a 20,000-year run for two configurations of parameters (one with weak atmospheric feedback and the other with stronger atmospheric feedback). We see two generic features of our model: a. The atmospheric tropics-extratropics interaction in our model (Alexander, 1992;Liu, 2012;Liu & Alexander, 2007;Pierce et al., 2000) automatically produces oscillations inT tr andT ex with opposite relative phase (as required by the observations; e.g., Liu, 2012;Newman et al., 2016;Pierce et al., 2000;Tatebe et al., 2013;L. Zhang & Delworth, 2015). b. The century-long SLP record analyzed by Minobe (1999) suggested that the bidecadal and pentadecadal peaks are period 3 phase locked and change phases synchronously with the same relative sign. Our basin-scale coupling mechanism model reproduces the period 3 phase locking (between the 20-and 60-year modes). Furthermore, our model also typically reproduces the observed configuration between the 20-and 60-year peaks (synchronous phase change with the same relative sign). Minobe and Jin (2004) analyzed a prototype equation with external forcing and delayed terms and have argued that in such model the phase reversal can occur only when the delay and forcing terms are in phase. This implies that any resonance additional to the forcing period should be an odd integer multiple of the forcing period. It further implies that there must be simultaneous phase reversal between the subharmonics and the forcing mode and that the modes have the same relative sign. Furthermore, Minobe and Jin (2004) have shown that weak delayed feedback is sufficient to produce the 3 times forcing period ratio. We approach the problem from a different angle than Minobe and Jin (2004). Building on what has been learned in the last couple of decades, we investigate the exact physical mechanism behind PIV and what the role of tropics-extratropics interaction within this mechanism is. This complex mechanism has broad implications, such as opposite phase synchronization betweenT ex andT tr and the basin-wide nature of PIV. It is interesting to see that the mechanism proposed in this study, without ocean tropics feedback (K O t = 0) and for the physically most interesting case of K A t K A e > K O e , reduces to the prototype equation of Minobe and Jin (2004) with the bidecadal mode acting as forcing. Our model can then be largely seen as a detailed physical stage (and support) for the more mathematically focused analysis of Minobe and Jin (2004). This also means that the three-period subharmonics and the synchronous phase locking between the bidecadal and pentadecadal peaks can be explained as in Minobe and Jin (2004). We have also found some interesting features not previously reported by Minobe and Jin (2004). Our spectra (see Figure 3) are naturally interpreted as a sequence of odd number multiples of the lowest-frequency mode (60-year-period peak). This is particularly important, since it includes 60∕5 = 12-year period, which for some parameter values becomes significant enough to be associated with the observed (Bruun et al., 2017;Li et al., 2013) decadal (13 ± 1) resonance. Also, as noted before, for some small subregions of parameters (still falling under the prototype equation of Minobe & Jin, 2004) our model produces spectra with integer (both even and odd) multiples of a lowest-frequency mode (i.e., the case in supporting information S1 shows this for a lowest-frequency mode with a 40-year period). Within this small subregion of parameters there are also cases in which the lowest-frequency mode is not an integer multiple of the 20-year peak (see Figure 4).
It is desirable to get better understanding about why a certain range of delay times produces the three subharmonics (60-year mode) and others do not. Minobe and Jin (2004) numerically investigated the problem, and their results suggest (similarly to our study) that the three subharmonics is produced in the interval of delay times (d) somewhere between one half and three halves of the forcing period (or one sixth and one half of the three subharmonics period). Following the logic of Minobe and Jin (2004), the three subharmonics constrains the delay time because the delay term has to be in the same phase as the solution before the solution reverses phase. Although this explains why it is possible for the three subharmonics to appear in the spectra, it does not explain why it exists. We suggest that the existence of three subharmonics can be partly explained by focusing on the subinterval of d around the value given as three fourths of the forcing period. In our study this is of particular interest, as it gives the median value for the ocean tunnel propagation time d = 15 years (Goodman et al., 2005;Gu & Philander, 1997). Let us assume that due to the external forcing the state oscillates with 20-year periodicity. The delay term effectively represents signal looping back and influencing its state 15 years later. By leaving imprint on the state 15 years later it gradually influences the system's future in a sequence of integer multiples of 15 years. Because of the 20-year periodicity the signal loops back phase-shifted. The signal is phase-shifted firstly by a quarter period (5 years) and then by the half-period (10 years, which means that it arrives with the opposite phase). The three subharmonics 60-year (15 × 4) mode then corresponds to the time scale where the 20-year oscillation amplifies, since its phase matches with the looping signal.
It is left for future work to further examine the mechanism used in this paper to shed more light on the nature of the PDO periodic spectral peaks. It will be also interesting to explore the teleconnection of these large-scale Pacific periodic oscillations on the other parts of the planet. Similar (20 and 60)-year peaks have been reported in the Northern Atlantic (Longhurst, 2010;Moore et al., 2017;Olsen et al., 2012), and it is certainly possible that Pacific climate variability forces such climatic processes in the Atlantic region (Knudsen et al., 2014;Marshall et al., 2001).

Conclusions
The simple toy model presented in this study produces the pentadecadal and decadal periodic resonances observed as a part of the PDO. Unlike the bidecadal periodicity that originates from the local ocean-atmosphere coupling in the extratropics, the two other resonances originate from the exchange between the tropics and the extratropics through the atmosphere, as well as through the long ocean tunnel. Our mechanism suggests that even a tiny feedback through the ocean tunnel (50 times smaller than the coupling through the atmosphere) is sufficient to provide the time scale for the pentadecadal mode. Our model produces period 3 locking between the bidecadal and pentadecadal modes, as well as synchronized phase reversal with the same relative sign (as observed in Minobe, 1999). These last results are not surprising as our physical mechanism mathematically corresponds (in a large subrange of most relevant cases) to the prototype equation investigated by Minobe and Jin (2004). We therefore show that the type of nonlinear feedback investigated by Minobe and Jin (2004) is highly relevant in the context of realistic model of interaction within the Pacific basin. Due to the role of teleconnections in the presented model, all the model resonances are necessarily basin-wide phenomena, as found, for example, in Bruun et al. (2017) and Li et al. (2013). Even for those who remain unconvinced about the existence of some of the PIV resonances the paper provides an unambiguous logical connection between their existence and a specific mechanism of interaction within the Pacific basin. discussions on the analysis of the basin-scale ENSO system. This research was supported by the NERC National Capability in Modelling programme at the Plymouth Marine Laboratory. JB work at the University of Exeter was supported by UK research councils models 2 decisions grant (M2DPP035). The data used in the paper can be reproduced using the R code provided in the supporting information.