Indo-Pacific variability on seasonal to multidecadal timescales. Part II: Multiscale atmosphere-ocean linkages

The coupled atmosphere-ocean variability of the Indo-Pacific on interannual to multidecadal timescales is investigated in a millennial control run of CCSM4 and in observations using a family of modes recovered in Part~I of this work from unprocessed SST data through nonlinear Laplacian spectral analysis (NLSA). It is found that ENSO and combination modes of ENSO with the annual cycle exhibit a seasonally synchronized southward shift of equatorial surface zonal winds and thermocline adjustment consistent with terminating El Nino and La Nina events. The surface wind patterns associated with these modes also generate teleconnections between the Pacific and Indian Oceans, leading to a pattern of SST anomalies characteristic of the Indian Ocean dipole. Fundamental and combination modes representing the tropospheric biennial oscillation (TBO) in CCSM4 are also found to be consistent with mechanisms for seasonally synchronized biennial variability of the Asian-Australian monsoon and Walker circulation. On longer timescales, the leading multidecadal pattern recovered by NLSA from Indo-Pacific SST data in CCSM4, referred to as west Pacific multidecadal mode (WPMM), is found to significantly modulate ENSO and TBO activity with periods of negative SST anomalies in the western tropical Pacific favoring stronger ENSO and TBO variability. Physically, this behavior is attributed to the fact that cold WPMM phases feature anomalous decadal westerlies in the tropical central Pacific and easterlies in the tropical Indian Ocean, as well as an anomalously deep thermocline in the eastern Pacific cold tongue. Moreover, despite the relatively low SST variance explained by this mode, the WPMM is found to correlate significantly with decadal precipitation over Australia in CCSM4.


Introduction
The Indo-Pacific domain is an arena where prominent modes of coupled atmosphere-ocean variability act on seasonal to multidecadal time scales. Notably, it is the region where El Niño-Southern Oscillation (ENSO), the dominant mode of interannual climate variability (Bjerknes 1969;Sarachik and Cane 2010), takes place. ENSO is well known to exhibit a phase synchronization with the seasonal cycle (Rasmusson and Carpenter 1982), with El Niño and La Niña events typically developing in the boreal summer, peaking during the following winter, and dissipating in the ensuing spring-a feature not fully explained theoretically despite decades of research. Recent studies (Stein et al. 2011(Stein et al. , 2014McGregor et al. 2012;Stuecker et al. 2013Stuecker et al. , 2015aRen et al. 2016) attribute this phase synchronization to quadratic nonlinearities in the coupled atmosphere-ocean system producing a cascade of frequencies known as ENSO tones, whose corresponding temporal and spatial patterns are called ENSO combination modes. In particular, these studies propose that the ENSO combination modes at the sum and difference of the dominant ENSO frequency and the annual cycle frequency play Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/ JCLI-D-17-0031.s1. an important role in restoring El Niño or La Niña conditions to neutral conditions through nonlinear interactions that project significantly on the corresponding surface atmospheric circulation and thermocline depth fluctuations. Another prominent ENSO feature under active investigation, as linear theories do not explain it well, is an asymmetry between El Niño and La Niña phases, manifested, for example, in positively skewed SST anomaly distributions (Burgers and Stephenson 1999;Kang and Kug 2002;An and Jin 2004) and longer average duration of La Niña conditions (Kessler 2002;Okumura and Deser 2010;DiNezio and Deser 2014).
In addition to ENSO and ENSO combination modes, other modes of interannual Indo-Pacific climate variability that have received significant attention are the Indian Ocean dipole (IOD; Saji et al. 1999;Webster et al. 1999) and the tropospheric biennial oscillation (TBO; Meehl 1987Meehl , 1993Meehl , 1994Meehl , 1997Li and Zhang 2002;Li et al. 2006). The IOD is generally defined as an opposite-sign SST anomaly pattern developing in the eastern and western Indian Ocean, whereas the TBO is associated with biennial variability of the Asian-Australian monsoon (Li 2010). Seasonal locking is a prominent aspect of both the IOD and the TBO that, as in the case of ENSO, is not fully understood despite being observed and studied intensively. The existence and significance of linkages between ENSO, the IOD, and the TBO is another topic of active research (Goswami 1995;Annamalai et al. 2003;Li et al. 2006;Dommenget 2011;Zhao and Nigam 2015).
The Indo-Pacific domain is also a region where pronounced decadal variability patterns arise, with the Pacific decadal oscillation (PDO; Mantua et al. 1997) and the interdecadal Pacific oscillation (IPO; Power et al. 1999) arguably being the most frequently studied modes. The domains of definition of these modes overlap to some degree (the PDO is defined over the North Pacific, whereas the IPO is defined from low-pass-filtered SST for the whole Pacific basin), and they are often considered to be a manifestation of the same phenomenon. The physical nature of these modes, and more broadly low-frequency Indo-Pacific fluctuations, has been investigated intensively over the last 20 years, and while a number of hypotheses have emerged, consensus on their validity is limited. Some studies find that the dominant decadal modes retrieved using classical methods such as empirical orthogonal function (EOF) analysis are predominantly ENSO residuals after low-pass filtering, in particular, residuals due to El Niño-La Niña asymmetries resulting from nonlinear dynamics (Rodgers et al. 2004;Vimont 2005;Sun and Yu 2009;Wittenberg et al. 2014). For instance, Vimont (2005) constructs a decadal ENSO-like signal through a linear combination of three interannual modes with no contribution of a decadal mode and demonstrates that a low-frequency residual may result from different spatiotemporal characteristics of faster components, which in turn can be associated with the ENSO precursor, peak, and ''leftover.'' Moreover, Rodgers et al. (2004) find that some of the dominant patterns of decadal SST and thermocline depth variability, as extracted via EOF analysis of low-passfiltered and meridionally averaged data over the equatorial Pacific, bear strong similarities to corresponding composite patterns of El Niño-La Niña asymmetries. In light of this finding, they argue that the dominant patterns recovered from decadal low-pass-filtered data are a consequence of ENSO asymmetries.
At the same time, many studies propose that the lowfrequency Indo-Pacific SST variability in both nature and coupled climate models requires a more complex explanation than solely spatial characteristics of lowpass-filtered ENSO and suggest independent mechanisms for multidecadal fluctuations of the tropical climate. Some associate decadal variability with one or more stochastic processes (Cane et al. 1995;Thompson and Battisti 2001), tropical/ENSO-like dynamics (Knutson and Manabe 1998;Yukimoto et al. 2000), or remote forcing from the midlatitudes (Kleeman et al. 1999;Vimont et al. 2003;Solomon et al. 2003). Frequently, decadal and multidecadal anomalies are not associated with a single physical mode, but are considered as multiple phenomena acting simultaneously Newman et al. 2016). Alternatively, they are viewed as fluctuations in the mean background state (Wang and Ropelewski 1995). Such changes may correspond to global warming (Timmermann et al. 1999), the Last Glacial Maximum , climate shifts (Ye and Hsieh 2006), or quasi-decadal variability (Hasegawa and Hanawa 2006), the latter most often attributed to the IPO (Salinger et al. 2001;Imada and Kimoto 2009). Decadal modulations of interannual modes such as ENSO, the TBO, and the IOD have also been studied (e.g., Arblaster 2011, 2012). Arguably, progress in elucidating the relation between decadal and interannual variability is hampered by the fact that the extraction of decadal modes through conventional data analysis techniques such as EOF analysis often requires ad hoc preprocessing of the data (e.g., low-pass filtering), compounding the physical interpretation of the results.
Irrespective of the physical nature of Indo-Pacific multidecadal variability and the statistical origin of the corresponding modes, its climatic impacts have been documented in numerous studies. Power et al. (1999) find that ENSO's impact on Australian precipitation depends strongly on the IPO phase, and da Silva et al. (2011) report an IPO modulation of ENSO's impact on South American precipitation. Dai (2013) has found a strong correlation of the IPO with anomalously dry and wet decades occurring over the western and central United States in the twentieth century. Moreover, multidecadal droughts in the Indian subcontinent and the southwestern United States have also been attributed to the IPO (Meehl and Hu 2006), and the rapid sea level rise (fall) observed over the last 17 years over the western (eastern) tropical Pacific (called east-west dipole by Meyssignac et al. 2012) has been associated with the PDO (e.g., Zhang and Church 2012). Meehl et al. (2013), England et al. (2015), and others have associated the decades of global warming hiatus observed during the twentieth century with negative phases of the IPO. At the same time, several studies (Solomon and Newman 2012;Seager et al. 2015) also stress the need to consider modes of Indo-Pacific internal variability other than the IPO (or trends associated with anthropogenic forcings) to build a more complete explanation of the ever-expanding observational record. In particular, the existence of unknown low-frequency modes of Indo-Pacific SST has been proposed by studies on recent sea altimetry (Palanisamy et al. 2015), warming hiatus (England et al. 2014), and African paleodroughts (Tierney et al. 2013).
In Slawinska and Giannakis (2017, hereafter Part I), we recovered a family of spatiotemporal modes of Indo-Pacific SST variability in models and observations using a recently developed eigendecomposition technique called nonlinear Laplacian spectral analysis (NLSA; Giannakis and Majda 2011, 2012a,b, 2013, 2014. A key feature of this method is that it can simultaneously recover modes spanning multiple time scales without prefiltering the input data. In Part I, we used NLSA to recover a family of Indo-Pacific SST modes from a 1300-yr integration of the Community Climate System Model, version 4 (CCSM4; Gent et al. 2011) and an 800-yr integration of the Geophysical Fluid Dynamics Laboratory Climate Model, version 3 (GFDL CM3; Griffies et al. 2011), as well as from monthly averaged HadISST data (Rayner et al. 2003) for the industrial and satellite eras. This family consists of modes representing (i) the annual cycle and its harmonics, (ii) the fundamental component of ENSO and its associated combination modes with the annual cycle (denoted ENSO-A modes), (iii) the TBO and its associated TBO-A combination modes, (iv) the IPO, and (v) a lowfrequency mode featuring a prominent cluster of SST anomalies in the western equatorial Pacific referred to as the west Pacific multidecadal mode (WPMM).
Compared to modes obtained via classical EOF-based approaches, the mode family identified in Part I has minimal risk of mixing intrinsic decadal variability with residuals associated with low-pass-filtered interannual modes (since NLSA operates on unprocessed data). Moreover, unlike EOF analysis, which mixes distinct ENSO combination frequencies into single modes (McGregor et al. 2012;Stuecker et al. 2013), the NLSAbased ENSO-A modes resolve individual combination tones and have the theoretically correct mode multiplicities. In particular, the leading four such modes form two oscillatory pairs with frequency peaks at either the sum or the difference between the fundamental ENSO frequency and the annual cycle frequency. In Part I, we saw that the SST anomalies due to the ENSO-A modes are especially prominent in the region west of the Sunda Strait employed in the definition of IOD indices (Saji et al. 1999), indicating that these modes may be useful IOD predictors.
The seasonal, ENSO, and ENSO-A modes recovered via NLSA were found to be in good qualitative agreement among the CCSM4, GFDL CM3, and HadISST datasets (the latter, for both industrial and satellite eras) and were also verified to be robust against variations of NLSA parameters and spatial domain, as well as measurement noise. The TBO, TBO-A, and WPMM (whose corresponding SST anomaly amplitudes are typically an order of magnitude smaller than ENSO) were best recovered in CCSM4, possibly because the CCSM4 dataset was the longest studied, although representations of the TBO (but not its combination modes) and modes resembling the WPMM were also found in GFDL CM3 and industrial-era HadISST. IPO modes were identified in CCSM4, GFDL CM3, and industrial-era HadISST, although in all cases the recovered patterns exhibited some degree of mixing with higher-frequency (e.g., interannual) time scales. The poorer quality of the TBO and decadal modes from GFDL CM3 and industrialera HadISST was at least partly caused by the shorter time span of these datasets compared to CCSM4, since a similar quality degradation was observed in experiments with subsets of the CCSM4 dataset. We were not able to recover TBO and decadal modes from the satellite-era HadISST dataset.
In this paper, we employ the mode family from Part I to study various aspects of coupled atmosphere-ocean variability in the Indo-Pacific domain on seasonal to multidecadal time scales. In addition, we study modes derived from Twentieth Century Reanalysis, version 2c (20CRv2c; Compo et al. 2011), SST data, which also yielded TBO, TBO-A, and WPMM of reasonably high quality. First, we examine the phase synchronization of ENSO to the seasonal cycle using composites and spatiotemporal reconstructions of SST, surface winds, and thermocline depth associated with the ENSO and ENSO-A modes and use these modes to study El Niño-La Niña asymmetries. We also study the reconstructed SST, surface wind, and precipitation patterns associated with the TBO and TBO-A modes and find that these patterns are consistent with the biennial variability of the Walker circulation and Asian-Australian monsoon expected for the TBO. Moreover, we examine the relationships between the WPMM and the ENSO and TBO modes and find that the WPMM has sign-dependent modulating relationships with ENSO and the TBO, enhancing the amplitude of the latter during its phase with negative western Pacific SST anomalies. Analogous modulation relationships are found to hold with decadal precipitation over Australia.
This paper is organized is as follows. In section 2, we describe the model and observational datasets studied in this work and give a brief overview of the NLSA calculations performed along with our procedure for building composites. In section 3, we present composites and spatiotemporal reconstructions of SST, surface winds, and thermocline depth associated with the ENSO and ENSO-A mode families and discuss El Niño-La Niña asymmetries and Pacific-Indian Ocean couplings associated with ENSO's spatiotemporal evolution. Section 4 presents and discusses analogous results for the TBO and its combination modes. In section 5, we present the corresponding composites and reconstructions for the WPMM and discuss its physical role in ENSO decadal variability and decadal precipitation over Australia. We state our conclusions in section 6. Details on the SST modes derived from 20CRv2c are included in appendix A. Appendix B contains spatiotemporal reconstructions of individual events created from satelliteera HadISST data. Movies illustrating the dynamical evolution of the NLSA modes are provided as supplementary material. A MATLAB code for NLSA algorithms, including example scripts for analyzing CCSM4 data, is available for download at the first author's personal website (http://cims.nyu.edu/;dimitris).

Dataset description and summary of analyses
We study the same Indo-Pacific domain as in Part I (ocean and atmosphere grid points in the longitudelatitude box of 288E-708W and 608S-208N), and analyze monthly averaged data from the same 1300-yr control integration of CCSM4 (CCSM 2010;Gent et al. 2011;Deser et al. 2012). In addition to the SST field studied in Part I, we examine surface zonal and meridional winds, precipitation, and maximum mixed layer depth (XMXL) output on the respective native atmosphere and ocean grids. In this control integration, both atmosphere and ocean grids have a 18 nominal resolution. We use XMXL as a proxy for thermocline depth.
We also study monthly averaged observational and reanalysis data for the industrial and satellite eras, that is, the periods January 1851-December 2012 and January 1983-December 2009, respectively; the latter period was also studied in Part I. For the industrial-era analysis, we use 20CRv2c data (NOAA 2014;Compo et al. 2011) for SST, surface winds, and precipitation, and we combine this data with Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCC) Historical Ocean Reanalysis (CHOR_RL; CMCC 2016; Yang et al. 2017) data for the 208C isotherm z T20 . We treat the latter as a proxy for thermocline depth. For the satellite-era analysis, we use the same monthly averaged HadISST data (Met Office Hadley Centre 2013;Rayner et al. 2003) as in Part I, and we also employ monthly averaged MERRA (NASA GMAO 2014;Rienecker et al. 2011) and CMCC Global Ocean Physical Reanalysis System (C-GLORS) data (CMCC 2014;Storto et al. 2016) for surface winds and z T20 , respectively. The 20CRv2c, CHOR_RL, HadISST, MERRA, and C-GLORS data are all output on uniform longitude-latitude grids of 28 3 28, 1 /28 3 1 /28, 18 3 18, 1 /28 3 2 /38, and 1 /48 3 1 /48 resolution, respectively. For conciseness, we collectively refer to the 20CRv2c, CHOR_RL, HadISST, MERRA, and C-GLORS data as observational data. Moreover, when convenient, we collectively refer to XMXL and z T20 as thermocline depth data.
Throughout this paper, we work with the NLSA eigenfunctions derived from SST data from CCSM4 and satellite-era HadISST, which are discussed in sections 4 and 7 of Part I, respectively. In addition, we employ NLSA eigenfunctions computed from SST data from 20CRv2c; these eigenfunctions were not discussed in Part I, but they include ENSO, ENSO-A, TBO, TBO-A, WPMM, and IPO modes as in CCSM4. Analogous modes were also identified from industrial-era HadISST data in Part I, but the WPMM and TBO modes from 20CRv2c have higher quality than their HadISST counterparts. Eigenfunction time series and frequency spectra for the 20CRv2c-derived modes are shown in appendix A. For completeness, in appendix B we show spatiotemporal reconstructions of individual El Niño and La Niña events based on the corresponding satelliteera HadISST modes. Note that we do not use surface wind, thermocline depth, or precipitation data in the computation of the NLSA eigenfunctions.
As discussed in Part I, an important parameter in NLSA is the embedding window length Dt, which controls time-scale separation in the recovered modes. In the CCSM4 analysis, we use a 20-yr embedding window, which allows us to recover annual cycle, ENSO, ENSO-A, TBO, TBO-A, IPO, and west Pacific multidecadal modes with high fidelity, as well as a family of secondary ENSO modes (and their associated annual cycle combination modes) with frequencies adjacent to the fundamental ENSO frequency and a more complex temporal evolution. In the observational analyses we must content ourselves with smaller embedding windows; in what follows, we use Dt 5 8 and 4 yr, respectively, for the 20CRv2c and HadISST analyses. As discussed in Part I, section 5a (see, in particular, Fig. 9 therein), at these smaller embedding window lengths, the ENSO and ENSO-A modes do not split into fundamental and secondary, and as a result have a more broadband character.
In what follows, all spatiotemporal reconstructions and explained variances are computed as described in section 2c of Part I. As in Part I, we use the term ''nonperiodic explained variance'' to refer to the variance explained by a subset of NLSA modes relative to the raw data, or the reconstructed data using all of our identified modes, after removal of the seasonal cycle. Using the reconstructed fields, we also create composites of SST, surface wind, mixed layer depth, and precipitation associated with the NLSA eigenfunctions. To construct these composites, we select a reference physical field of interest averaged over a region representative of the family of modes under study. Optionally, for seasonally locked patterns, we also select a reference month (e.g., in the case of the ENSO and ENSO-A eigenfunctions, we use averaged SST over the Niño-3.4 region with January as the reference month). We then identify all time stamps in the reconstructed data for which the averaged field exceeds a preselected signed threshold (e.g., one standard deviation of the averaged field over the region of interest), subject to the additional condition that the calendar month of the selected time stamps is the same as the reference month if a reference month is used. Composites are then created for every field of interest by averaging over the identified time stamps and over those time stamps shifted by monthly leads or lags in a given range. For example, in the case of ENSO and ENSO-A composites, we examine 4-month leads and lags around January as in Deser et al. (2012). We label these composites using a superscript to denote lead-lag in years relative to the reference month; for example, in the case of the ENSO composites just mentioned, April 0 would stand for the April of the year of the reference January.
Note that it is important to identify the time stamps for compositing using reconstructed physical fields as opposed to NLSA eigenfunctions, since (as with other techniques that employ time-lagged embedding) the temporal correspondence between the phase of an oscillatory pattern in eigenfunction space and the manifestation of this pattern in physical space is not known a priori. It should also be noted that, because of the high temporal coherence of the NLSA eigenfunctions, the spatial patterns in our composites are qualitatively similar to individual event reconstructions, such as those shown in appendix B. Such reconstructions can be also directly visualized in the movies in the supplementary material.
3. Atmosphere-ocean covariability of ENSO and ENSO combination modes a. Surface circulation and thermocline patterns It is well known that interannual SST variability in the Pacific Ocean is strongly coupled with surface atmospheric circulation and thermocline variability (Alexander et al. 2002;Lengaigne et al. 2006;Vecchi 2006;McGregor et al. 2012;Stuecker et al. 2013Stuecker et al. , 2015a. In this section, we analyze the nature of these couplings in the context of the ENSO and ENSO-A combination modes recovered from the model and observational datasets described in section 2. In particular, we examine the evolution of SST, surface winds, and thermocline depth associated with this family of modes. We create composites via the approach described in section 2 using January SST anomalies reconstructed jointly from the ENSO and ENSO-A modes and averaged over the Niño-3.4 box to identify significant El Niño and La Niña events. The (signed) threshold used to identify significant events is equal to one standard deviation of the averaged SST anomalies. Identified events are subsequently employed to build lead-lag composites every 4 months for the period starting in November 22 and ending in July 0 .
Here, we discuss primarily El Niño and La Niña composites created using SST data from CCSM4 and 20CRv2c as inputs; in the latter case, the resulting thermocline depth composites are based on CHOR_RL z T20 data. In addition, the evolution of individual events, shown in movies 1 and 3 in the supplementary material for CCSM4 and 20CRv2c, respectively, is also instructive. For instance, in movie 1, ENSO activity is strong in the CCSM4 simulation period January 1174-December 1178 (and the WPMM is in its cold phase; see section 5). Meanwhile, movie 3 contains the strong 1997/ 98 El Niño and the subsequent La Niña of 1999/2000. Individual events reconstructed from satellite-era HadISST are discussed in appendix B and visualized in movie 4 in the supplementary material. Note again that, because of the high temporal coherence of the ENSO and ENSO-A modes from NLSA, individual significant ENSO events bear strong similarities to the composites.
We first consider the El Niño composites from CCSM4. As shown in Fig. 1, positive SST anomalies begin building up in the eastern Pacific in the prior February-April (FMA 21 ). During that period, anomalous westerly winds develop in the western tropical Pacific, and a positive XMXL anomaly develops in the eastern equatorial Pacific, where it continues to increase in the ensuing months. In movie 1 in the supplementary material, it can be seen that the eastern Pacific positive XMXL anomaly develops as a traveling anomaly originating in the western Pacific. The eastern Pacific SST anomalies also increase in the months following FMA 21 , until El Niño peaks in November 21 -December 21 , reaching around 1-K SST anomalies in the centraleastern equatorial Pacific. Around the time of the El Niño peak, the anomalous western Pacific westerlies undergo a pronounced southward shift to about 108S, which is accompanied by a return of the eastern equatorial Pacific thermocline depth to near-climatological levels by July 0 -August 0 . The corresponding CCSM4-derived La Niña composites (see Fig. 2) behave to a large extent as sign-reversed analogs of the El Niño composites described above.
More detailed aspects of the thermocline response associated with the southward shift of meridional winds is captured in reconstructions (movie 1 in the supplementary material) are based on the ENSO-A modes alone. For instance, it can be seen that as the anomalous westerlies shift to the south of the equator in February-March 1176 after a strong La Niña event, a positive XMXL anomaly is generated in the central equatorial Pacific that subsequently propagates to the west, reaching the Maritime Continent in August-September 1176, where it contributes to the return of the previously negative XMXL anomaly there (associated with the fundamental ENSO mode) to near-climatological levels. That ENSO-A anomaly is then reflected back and begins propagating toward the east, reaching the eastern Pacific in March-April 1177-approximately during the time when an El Niño event begins to develop. To create these composites, significant El Niño events were selected using a 1s threshold of reconstructed January SST anomalies averaged over the Niño-3.4 region. The composite evolution based on these events is plotted every 4 months starting in November 22 14 months prior to the reference January and ending six months later in July 0 . Statistical significance of the composite patterns was assessed via a z test with significance level 0.05. Gridpoint anomalies not meeting this test were masked in white. The evolution shown here illustrates the (a)-(i) development, (j)-(l) peak, and (m)-(r) dissipation of El Niño events. Notice the anomalous westerly winds in the western equatorial Pacific in (e),(h),(k) and west-east anomalous thermocline depth gradient in (f),(i),(l) as El Niño conditions develop in March 21 -November 21 , followed by a southward shift of zonal winds in (n) and return of climatological thermocline conditions in (o),(r) as the event decays in February 0 -July 0 . This process can also be visualized in movie 1 in the supplementary material.
Overall, the atmosphere-ocean process described above is in good agreement with the predictions of ENSO combination mode theory (Stein et al. 2011(Stein et al. , 2014McGregor et al. 2012;Stuecker et al. 2013Stuecker et al. , 2015a. In particular, using intermediate-complexity atmosphere and ocean models, McGregor et al. (2012) propose a mechanism whereby the SST anomalies from a developed ENSO event induce a Gill-type Rossby-Kelvin wave in the lower troposphere, which leads to an ageostrophic boundary layer flow. The latter is stronger in the South Pacific than in the North Pacific due to weaker climatological wind speeds and associated Ekman pumping in DJF and MAM, giving rise to a seasonally locked southward shift of meridional winds. Moreover, using a reduced ocean model, they propose that this southward shift of meridional winds causes two distinct oceanic responses involving generation of equatorial Kelvin waves or equatorial heat recharge/discharge that are responsible for terminating ENSO events.
This seasonally locked southward shift of meridional winds preceding ENSO termination is well captured in Figs. 1 and 2 and movie 1 in the supplementary material. As discussed in Part I, a representation of the southward zonal wind shift can also be obtained through the two leading EOFs of atmospheric circulation fields (McGregor et al. 2012;Stuecker et al. 2013), but this representation mixes the three distinct ENSO and ENSO-A frequencies into a single principal component pair. In contrast, the NLSA eigenfunctions recovered from SST isolate the components of ENSO and ENSO-A variability through three distinct pairs of modes evolving at the correct theoretical frequencies. This more detailed and theoretically consistent decomposition should be useful for future studies of ENSO termination.
Next, we consider the SST, surface wind, and z T20 composites associated with the ENSO and ENSO-A modes extracted from the 20CRv2c and CHOR_RL datasets. As shown in the El Niño composites in Fig. 3 (and in the spatiotemporal reconstructions in movie 3 in the supplementary material), SST anomalies preceding the El Niño peak phase in December 21 -January 0 begin developing in the equatorial eastern Pacific in December 22 through January 21 -February 21 , and continue growing until the event's peak. The buildup of SST anomalies is accompanied by anomalous westerlies in the western Pacific and the establishment of an anomalous shallow (deep) thermocline in the western (eastern) tropical Pacific. The anomalous zonal winds migrate to the south (in a process that begins ;2 months prior to the event peak), reaching a maximum southerly latitude of about 108S in March 0 , at which point they start to diminish and the SST and z T20 anomalies in the eastern Pacific decay to climatological levels. ENSO-neutral conditions are established by September 0 (not shown). The La Niña composites (Fig. 4) exhibit a sign-reversed pattern, but are otherwise broadly similar to the El Niño composites described above. As discussed in appendix B and visualized in movie 4, the El Niño and La Niña events reconstructed from satellite-era HadISST data (shown in Figs. B1 and B2, respectively) are broadly consistent with their 20CRv2c counterparts.
In general, the composites based on CCSM4 and observations are in good agreement over the equatorial Pacific where the majority of ENSO activity takes place-this is to be expected since CCSM4 is known to simulate realistic ENSO variability (although with notable biases in certain aspects, such as a 30% amplitude overestimate; Deser et al. 2012). Still, important differences remain, including a westward displacement of the anomalous divergence associated with the ENSO response of the Walker cell in CCSM4 relative to 20CRv2c and HadISST. This region can be identified as the center of a prominent dipole of zonal wind anomalies developing over the western Pacific warm pool and the Maritime Continent during the growth and peak phases of El Niño and La Niña events. In CCSM4, the center of that dipole is located at the longitudes of Sumatra (;1008E), but in 20CRv2c and HadISST it is located over New Guinea (i.e., 308-408 to the east); compare, for example, the November 21 composites in Fig. 1 (CCSM4) with Figs. 3 (20CRv2c) and B1 (HadISST). Consistent with this discrepancy is a sign discrepancy of the SST anomalies in the seas north of Sumatra. There, during peaking El Niño events, the SST anomalies are negative in the CCSM4 composites, but positive in the 20CRv2c composites and HadISST reconstructions. An analogous but opposite-sign discrepancy takes place during La Niña events. These discrepancies between CCSM4 and HadISST can also be seen in the composites in Figs. 9-12 in Deser et al. (2012), which are based on Niño-3.4 indices.

b. Pacific-Indian Ocean couplings
Besides the tropical Pacific Ocean, the ENSO and ENSO-A modes have strong surface atmospheric and thermocline response in the Indian and South Pacific Oceans. As shown in Fig. 1 and movie 1 in the supplementary material, coincident with peaking El Niño events is the formation of a prominent anticyclonic surface circulation in the southeastern Indian Ocean off the western coast of Australia (see, e.g., the November 21 -March 0 composites). These patterns, along with the corresponding thermocline anomalies, are consistent with the development of (i) positive SST anomalies in the western Indian Ocean (which are especially pronounced off Africa's east coast) via downwelling and advection of warm water from the Maritime Continent and (ii) negative SST anomalies in the eastern Indian Ocean via upwelling and advection of cold water from the southern Indian Ocean. This dipole SST pattern peaks in February 0 and decays together with the main El Niño SST anomalies in the ensuing boreal spring. As noted in Part I, such an SST pattern is characteristic of the IOD and explains approximately 35% of the raw nonperiodic variance in the Indian Ocean regions used to define IOD indices (Saji et al. 1999). During La Niña events ( Fig. 2 and movie 1 in the supplementary material), the southeastern Indian Ocean anticyclone is replaced by a cyclone (centered at the same location), and the positive SST anomalies in the western Indian Ocean are replaced by positive SST anomalies in the eastern Indian Ocean, especially in the region west of the Sunda Strait.
In the case of 20CRv2c data, significant Pacific-Indian Ocean couplings also take place during El Niño (Fig. 3) and La Niña ( Fig. 4) events (see also movie 3 in the supplementary material), although with some notable differences compared to the CCSM4 data. In particular, as stated in section 3a, the anticyclonic surface circulation over the southeastern Indian Ocean is shifted to the east in 20CRv2 and its zonal winds tend to be weaker in its southern branch than in its northern one (compare, e.g., the November 21 composites in Figs. 1 and 3). Anomalies in Indian Ocean thermocline are also shifted eastward (toward the east-northeast of Madagascar). Nevertheless, the dipole-like zonal pattern that they develop is consistent with upwelling/downwelling due to zonal advection in response to the atmospheric circulation anomalies described above. Also consistent with the eastward shift of surface wind anomalies is that in the months during and after El Niño and La Niña peaks, the SST field develops a monopole anomaly pattern as opposed to the dipole pattern in CCSM4. In fact, weak SST anomaly patterns with a sign reversal over the Indian Ocean can be solely seen during the months preceding El Niño and La Niña peaks, for example, July-November 1997 and 1999 in Figs. B1 and B2, respectively. A similar discrepancy between the Indian Ocean SST anomalies in CCSM4 and the El Niño and La Niña events reconstructed from the satellite-era HadISST data can be seen in Figs. B1 and B2, respectively. Overall, FIG. 4. As in Fig. 3, but for La Niña. This process can be more directly visualized in movie 3 in the supplementary material.
these discrepancies indicate that there may be significant biases of Indian Ocean SST variability in this model.

c. El Niño-La Niña asymmetries
In this section, we demonstrate the skill of the NLSAderived ENSO modes in capturing the higher-order statistics of ENSO, and in particular the positive SST and zonal wind skewness associated with asymmetries between El Niño and La Niña events (Burgers and Stephenson 1999;Kang and Kug 2002;Kessler 2002;An and Jin 2004;Okumura and Deser 2010;DiNezio and Deser 2014). In what follows, we study these statistics by means of skewness maps for reconstructed SST and surface zonal wind fields (Figs. 5 and 6, respectively) constructed from various combinations of ENSO modes derived from CCSM4, 20CRv2c, and industrial-era HadISST data. We also examine the corresponding histograms computed from averaged reconstructed anomalies over ENSO active regions (Figs. 7 and 8). In the case of SST, the region for averaging is the eastern half of the Niño-3 box (i.e., 58S-58N, 1208-908W); in the case of surface winds, we average over the Niño-4 box to capture western Pacific atmospheric variability important to the ENSO life cycle. Table 1 summarizes the statistics of these anomalies, as well as the standard deviation, kurtosis, and relative occurrence of El Niño versus La Niña conditions, as described below. Throughout, we assess the statistical significance of these results using standard errors computed via bootstrapping.
To construct the maps in Figs. 5 and 6, we first compute the third moment m ENSO,i 5 (y ENSO,i 2 y ENSO,i ) 3 of the reconstructed anomalies y ENSO,i associated with a selected ENSO mode family at every grid point i, and then plot the quantity g ENSO,i 5 sgn(m ENSO,i )jm ENSO,i j 1/3 (here, y stands for either SST or surface zonal wind anomalies; see section 2c in Part I for additional details on reconstructed fields). Note that g ENSO,i has the same physical dimension as y ENSO,i and the same sign as m ENSO,i . Hereafter, we refer to this quantity as skewness. An analogous quantity to g ENSO,i was used to quantify skewness by DiNezio and Deser (2014). Note that for the surface wind reconstructions from the HadISST modes, we used only the portion of the timespan of these modes overlapping with the MERRA reanalysis dataset (see section 2).  Table 1). Grid points with jg ENSO,i j , « ENSO,i were masked in white.
The histograms in Figs. 7 and 8 were constructed by binning ENSO reconstructed SST and surface zonal wind anomalies averaged over the ENSO active regions described above. For these spatially averaged anomalies, we also compute normalized skewness coefficients G ENSO 5 (y ENSO 2 y ENSO ) 3 /s 3 ENSO , where y ENSO denotes spatially averaged anomalies, and s ENSO denotes the corresponding standard deviations. Note that the calculation of G ENSO is numerically well conditioned since s ENSO is always statistically significantly greater than zero (which would not be the case for certain gridpoint standard deviations). To assess the probability of occurrence of large deviations from the mean, we also compute the kurtosis k ENSO 5 (y ENSO 2 y ENSO ) 4 /s 4 ENSO . In addition, we compute the relative occurrence frequency of El Niño and La Niña conditions as the ratio between the total time that the Niño-3 reconstructed SST anomaly is positive and the total time that it is negative. In the case of the zonal wind reconstructions, we compute the analogous ratio between Niño-4 westerly and easterly conditions. The values of these statistics are listed in Table 1. 1) CCSM4 DATA WITH A 20-YR EMBEDDING WINDOW As discussed in section 5a of Part I, depending on the length of the lag embedding window, besides the fundamental ENSO and ENSO-A modes the NLSA spectrum may also contain secondary ENSO modes with frequencies adjacent to the main ENSO band. These modes form degenerate oscillatory pairs, and appear as the embedding window Dt in NLSA is increased from interannual to decadal lengths. As stated in section 2, the embedding window used to compute modes from the CCSM4 dataset is Dt 5 20 yr, and a number of clear secondary modes are present in the spectrum; see Fig. 9 in Part I. This behavior is consistent with NLSA progressively resolving a broadband spectrum (here, due to ENSO) into distinct modes as Dt grows (Berry et al. 2013;Das and Giannakis 2017, manuscript submitted to Nonlinearity). In particular, the primary ENSO modes capture a relatively narrowband, regular ENSO signal with prominent modulations on decadal time scales, whereas the secondary ENSO modes exhibit a less coherent behavior with spectral power concentrated around the fundamental ENSO frequency and weaker low-frequency amplitude modulations.
The time series and frequency spectrum of a secondary ENSO mode is displayed in Fig. 9, where the fundamental ENSO mode is also shown in comparison. The corresponding spatial patterns (not shown here) exhibit ENSO-like SST anomalies in the Pacific Ocean. It should be noted that the secondary ENSO modes appear in the NLSA spectrum together with their associated combination modes with the annual cycle, and these modes were also used to compute the skewness maps and histograms in Figs. 5-8 involving secondary ENSO modes. Figure 5c shows a skewness map derived from the fundamental and leading four pairs of secondary ENSO modes together with the associated ENSO combination modes recovered from CCSM4 using a 20-yr embedding window. There, skewness is positive over the eastern Pacific cold tongue region where the majority of ENSO activity occurs, and takes negative values in the western Pacific. Moreover, the Indian Ocean features a dipolar skewness pattern with predominantly positive (negative) values in the western (eastern) parts of the basin. These features are qualitatively consistent with skewness patterns constructed from detrended observational SSTs (see, e.g., Fig. 3a in Burgers and Stephenson 1999). The positively skewed SST statistics over ENSO active regions are also evident in the distribution of the reconstructed SST anomalies over the eastern half of the Niño-3 box, shown in Fig. 7c. More quantitatively, we find that the skewness and normalized skewness of this variable are g ENSO 5 0.76 K and G ENSO 5 0.33, respectively, both of which are statistically significantly different from zero (see Table 1). This value for g ENSO is higher than the 0.19 K value reported by DiNezio and Deser (2014) for CCSM4; this difference could be due to issues such as the choice of averaging domain and bandpass filtering of the data. Our value for G ENSO is consistent with reported Niño-3 skewness values of CMIP5 models (e.g., Fig. 6a in Weller and Cai 2013).
As shown in Table 1, the averaged SST anomalies have a weakly platykurtic distribution with kurtosis k ENSO 5 2.90, that is, more frequent, but less extreme, departures from the mean than Gaussian random variables with k ENSO 5 3. Moreover, the rate of occurrence of El Niño versus La Niña conditions shows a 52/48 preference for La Niña; this is consistent with known behavior in CCSM4 (Deser et al. 2012;DiNezio and Deser 2014), as well as observations (Kessler 2002;Okumura and Deser 2010).
As perhaps one might expect, using fewer secondary ENSO modes or a longer embedding window leads to lower skewness and variance values from those shown in Figs. 5c and 7c. Thus, to assess the contributions of the individual modes in the ENSO family used for reconstruction, we examine skewness maps and averaged SST anomaly histograms based on either the primary ENSO and ENSO-A modes (Figs. 5a and 7a), or the secondary ENSO and ENSO-A modes alone (Figs. 5b and 7b). Interestingly, the skewness maps and histograms based on the primary modes exhibit no statistically significant skewness, and the corresponding rate of occurrence of El Niño and La Niña conditions is very FIG. 7. Histograms of SST anomalies dSST due to ENSO modes averaged over the eastern half of the Niño-3 region. The ENSO modes used in each panel are the same as the corresponding panels in Fig. 5. The histograms were computed by binning averaged SST anomalies over 51 bins of uniform width 3/25 K with centers dSST 225 , . . . , dSST 25 arranged symmetrically about zero. Red bars show the corresponding normalized bin counts r 225 , . . . , r 25 (range of values displayed on left vertical axes). Blue bars show the differences r i 2 r 2i for i 2 f1, . . . , 25g between the bin counts at negative and positive bins (range of values displayed on right vertical axes) to highlight asymmetries. nearly equal. At the same time, the SST reconstructions based on the primary modes have a leptokurtic distribution (k ENSO 5 3.93), which is consistent with an oscillator evolving coherently in a typical amplitude state and undergoing comparatively rare excursions to large/ small amplitude values (see Fig. 9b). In contrast, the reconstructions based on the secondary modes have a statistically significant positive skewness (G ENSO 5 0.30), as well a near-Gaussian kurtosis value (k ENSO 5 3.00) characteristic of a more noisy signal. These reconstructions also show a 53/47 preference for occurrence of La Niña conditions over El Niño conditions. As is evident in Fig. 7, the standard deviation of the reconstructed anomalies with the secondary ENSO modes included is considerably higher than that from the fundamental modes alone; specifically, in the former case we have s ENSO 5 0.69 K, whereas in the latter case that value drops to 0.31 K (see Table 1). Overall, these results indicate that NLSA with long embedding windows separates ENSO variability into a temporally coherent FIG. 8. As in Fig. 7, but for ENSO reconstructed surface zonal wind anomalies du averaged over the Niño-4 region. In this case, the histogram bin widths are 1/5 m s 21 .
TABLE 1. Summary of statistics of reconstructed ENSO SST and surface meridional wind u anomalies averaged over the eastern half of the Niño-3 region for SST and the full Niño-4 region for u, showing standard deviation s ENSO , skewness g ENSO , normalized skewness G ENSO , kurtosis k ENSO , and relative occurrence frequency of El Niño and La Niña conditions for SST anomalies and westerly and easterly u anomalies. The ENSO modes used for reconstruction are the same as in Figs. 7 and 8 To explore the physical origin of this behavior, we now examine the corresponding skewness maps (Fig. 6) and Niño-4 histograms (Fig. 8) for reconstructed zonal surface winds. In particular, in both models and observations, it has been established that the integrated impact of highfrequency atmospheric dynamics such as westerly wind bursts can disturb the regularity of ENSO life cycle, especially its El Niño phase (Fedorov 2002;Fedorov et al. 2003;Rong et al. 2011), leading to positively skewed ENSO SST statistics. While our monthly averaged data preclude us from directly identifying short-lived events such as westerly wind bursts, it is nevertheless interesting to examine whether our reconstructed surface circulation fields have at least qualitative features that are consistent with the SST results discussed above and can be explained as signatures of high-frequency atmospheric events influencing ENSO variability. Indeed, as can be seen in Fig. 6c, the skewness maps based on surface zonal winds reconstructed using both fundamental and secondary ENSO modes exhibit positive values in central equatorial Pacific regions associated with anomalous westerlies such as westerly wind burst and other high-frequency atmospheric phenomena impacting ENSO, but the reconstructions based on the fundamental modes alone have negligible skewness. Similarly, there is a marked increase in skewness between the histograms in Figs. 8a and 8c, which were computed using the fundamental ENSO modes alone and the union of fundamental and secondary ENSO modes, respectively. Overall, these results are consistent with the hypothesis that the secondary ENSO modes capture the aggregate effect of high-frequency atmospheric variability with a positive skewness (westerly winds) and an associated positive SST skewness (strong El Niño events).

2) OBSERVATIONAL AND CCSM4 DATA WITH A 4-YR EMBEDDING WINDOW
In the case of the 20CRv2c data, our NLSA modes were computed using an 8-yr embedding window (see appendix A), and as a result they exhibit a less clear distinction between primary and secondary ENSO modes. Here, we have computed skewness maps and histograms using the primary ENSO and ENSO-A modes together with a single family of secondary modes analogous to those shown in Fig. 9. As shown Fig. 5e, the collection of these modes recovers successfully the basic features seen in skewness maps from raw observational SST data, having positive and negative values over the eastern and western tropical Pacific, respectively (in the eastern Pacific we have g ENSO 5 0.25 K and G ENSO 5 0.07; see Table 1). In separate calculations, we have verified that the skewness maps obtained from the primary ENSO modes from 20CRv2c already exhibit statistically significant positive skewness in the eastern Pacific (in contrast to the insignificant skewness of the primary modes from CCSM4 with the 20-yr embedding window). Qualitatively similar results (Fig. 5f), although with somewhat higher skewness values in the eastern equatorial Pacific (g ENSO 5 0.31 K and G ENSO 5 0.25), are also obtained from industrial-era HadISST data with the 4-yr embedding window used in Part I.
The skewness maps from 20CRv2c and industrial-era HadISST and MERRA are also fairly consistent with respect to the zonal surface winds (Figs. 6e,f). Namely, they both feature a skewness dipole centered over New Guinea (with positive values in the western Pacific and negative values in the Indian Ocean), although in this case 20CRv2c exhibits larger skewness values (g ENSO 5 0.32 K versus 0.17 K in HadISST and MERRA, although the latter dataset actually has larger normalized skewness values; see Table 1). The skewness maps in Figs. 6e,f have more significant differences in the subtropics and southern midlatitudes; for example, in the southeastern Pacific off the South American coast where 20CRv2c shows positive skewness, whereas the skewness from HadISST and MERRA is weaker and generally negative there.
The skewness results from the observational data have broadly consistent patterns, although different magnitudes, compared with those obtained from CCSM4 using the 20-yr embedding window and both the primary and secondary ENSO modes (Figs. 5c and 6c). The results from the observational datasets are also in reasonably good agreement with the skewness maps obtained from CCSM4 using a 4-yr embedding window (Figs. 5d and 6d). In particular, in the CCSM4 analysis with Dt 5 4 yr we only use primary ENSO modes, and the fact that the resulting skewness maps agree in pattern with the CCSM4 analysis with Dt 5 20 yr using both the primary and secondary ENSO modes (but not the primary modes alone) is indicative of the statement made earlier and in Part I that as the embedding window length increases, the recovered ENSO modes split into primary and secondary modes, capturing distinct aspects of ENSO variability. It should be noted that the CCSM4 reconstructions performed with Dt 5 20 yr and both the primary and secondary ENSO modes lead to significantly higher skewness values than their Dt 5 4 yr counterparts, but that difference is likely at least partly caused by the fairly large number of secondary modes used for reconstruction in the former case. It is also worthwhile to note that the primary ENSO modes from HadISST produce significantly higher skewness and kurtosis values than the primary modes from CCSM4 with Dt 5 4 yr (see Table  1), despite these analyses having the same embedding window length.
As perhaps expected from the discussion in section 3b, a region where the CCSM4 skewness maps (for both embedding window lengths used here) exhibit strong differences from the observational maps is the Indian Ocean, where CCSM4 displays a dipole pattern of SST skewness and a strong negative skewness of surface zonal winds, whereas both 20CRv2c and HadISST show a monopole pattern of weakly positive SST skewness and a weaker negative skewness of surface zonal winds.

3) DISCUSSION
Before closing this section, we discuss our results in the context of the findings of other studies on ENSO dynamics and skewness. As mentioned above, the coupling of high-frequency atmospheric dynamics to ENSO has been a topic of significant interest in both model and observational studies. For example, Christensen et al. (2017) report climate simulations with stochastic parameterization that improve the statistics of the ''atmospheric noise'' associated with westerly wind bursts, leading to amplification of El Niño-La Niña asymmetries. At the same time, many aspects of the atmosphere-ocean coupling impacting ENSO asymmetrically remain undetermined. A number of studies have emphasized the role of the background SST state (Gebbie et al. 2007;Tziperman and Yu 2007;Jin et al. 2007;Levine et al. 2016;Christensen et al. 2017) and the associated nonlinear atmospheric response through organized convection (Neale et al. 2008;Slawinska et al. 2014Slawinska et al. , 2015 in generating atmospheric disturbances that act on shorter time scales and impact ENSO asymmetrically. Other than convection, proposed mechanisms producing skewed ENSO statistics include nonlinear dynamical heating (Jin et al. 2003), nonlinear thermocline feedback (DiNezio and Deser 2014), and nonlinear wind stress response to SST (Kang and Kug 2002).
Given that NLSA is capable of extracting and separating primary ENSO modes that correspond to narrowband quasiperiodic parts of ENSO variability from secondary modes that can be associated with more broadband, irregular signals and ENSO's positive skewness (at least in CCSM4), we believe that the method can provide a useful framework for improved characterization of nonlinear ENSO dynamics, and in particular the details of the relevant atmosphere-ocean coupling. Such a framework may also contribute to reduction of ENSO biases in current-generation climate models (Zhang and Sun 2014), as well as potential improvements in ENSO predictability (Astudillo et al. 2017;Petrova et al. 2017), where the infamous ''spring barrier'' has been attributed by some studies to ENSO's coupling to highfrequency atmospheric disturbances and their nonlinear interactions with the annual cycle (Fedorov et al. 2003;Lopez and Kirtman 2014;Levine and McPhaden 2015).

Atmosphere-ocean covariability of the TBO and TBO combination modes
In this section, we study SST, surface wind, and precipitation patterns associated with the TBO and TBO-A modes. As discussed in Part I (in particular, section 5), the skill of NLSA in extracting low-variance patterns such as the TBO decreases by shortening the time span of the data or by adding observational noise. In particular, although a TBO mode was successfully recovered from industrial-era HadISST data, we were not able to identify one in the satellite-era HadISST dataset. Moreover, combination modes between the TBO and the annual cycle, representing potentially important aspects of this pattern such as its seasonal locking, were not captured by NLSA in either the industrial-era or the satellite-era HadISST data. Given that, and the fact that HadISST does not provide an atmospheric circulation and precipitation product, here we do not study the TBO with this dataset further. However, in the case of the 20CRv2c data, NLSA recovers both the TBO and an associated TBO-A combination mode with the annual cycle (the latter with some mixing between the TBO and TBO-A frequencies). The basic temporal properties of these modes are described in appendix B. Figure A1 shows the corresponding eigenfunction time series and frequency spectra.
In what follows, we study the NLSA-derived TBO and TBO-A modes from CCSM4 and 20CRv2c by means of SST, surface circulation, and precipitation composites. These composites were created as described in section 2, using a one standard deviation (1s) threshold of July precipitation anomalies, averaged over a Western Ghats domain (88-228N, 728-788E) to identify significant events. The lead-lag period covered is January 0 -July 11 . The composites from CCSM4 and 20CRv2c are displayed in Figs. 10 and 11, respectively. In addition, reconstructions of specific TBO events in CCSM4 and 20CRv2c can be seen in movies 2 (e.g., simulation period December 1169-December 1171) and 3 [e.g., the year 1987 (1988) associated with a weak (strong) Indian monsoon] in the supplementary material, respectively. Note that in Figs. 10 and 11 we focus on an Indian and western Pacific Ocean subdomain to bring out spatial detail in the regions where the TBO process is predominantly active. Reconstructions over the full Indo-Pacific domain studied here can be found in movies 2 and 3 in the supplementary material.
We first discuss the composites from CCSM4. As shown in Figs. 10a-c, January 0 corresponds to an anomalously weak Australian monsoon phase, featuring horizontal surface divergence (subsidence) and negative SST anomalies over the western Pacific, together with negative precipitation anomalies stretching south from the Maritime Continent to northern Australia. This warm pool divergence drives anomalous easterlies over the eastern and central Indian Ocean and anomalous westerlies over the central tropical Pacific. The anomalous easterlies correspond to a weakened western Walker cell, enabling the development of positive SST anomalies in the western Indian Ocean and FIG. 10. Composites of (a),(d),(g),( j) SST; (b),(e),(h),(k) surface wind; and (c),(f),(i),(l) precipitation anomalies associated with the TBO and TBO-A mode families extracted from CCSM4 SST data via NLSA. The composites are plotted every 6 months using July as the reference month, and span one cycle of the TBO. Statistical significance of the composite patterns was assessed via a z test as in Fig. 1. The phases shown here feature (a)-(c) an anomalously weak Australian monsoon in January 0 , (d)-(f) an anomalously strong Indian monsoon in July 0 , (g)-(i) an anomalously strong Australian monsoon in January 11 , and ( j)-(l) an anomalously weak Indian monsoon in July 11 . This process can also be visualized in movie 2 in the supplementary material.
intensification of deep convection (as evidenced by positive precipitation anomalies in east Africa). Meehl and Arblaster (2002) associate this pattern with an atmospheric Rossby wave response over Asia and positive temperature anomalies over the Indian subcontinent, preconditioning for a strong Indian monsoon in the upcoming boreal summer. This is consistent with the July 0 composites that exhibit an anomalously strong Indian monsoon (see positive precipitation anomalies over the Indian subcontinent and Southeast Asia in Fig. 10f). During that time, the western Walker cell is anomalously strong (Fig. 10e), and anomalous precipitation occurs over the tropical eastern Indian Ocean and western Pacific driven by positive SST anomalies there (Fig. 10d). Meehl and Arblaster (2002) attribute the development of these anomalies to thermocline adjustment resulting from eastward-propagating Kelvin waves generated by the anomalous circulation in the preceding winter and spring months.
The positive SST anomalies in the waters off Australia and the Maritime Continent persist over the ensuing boreal fall and winter (Fig. 10g). As a result, as precipitation shifts southeastward from India and Southeast Asia in response to the annual cycle of insolation, it becomes anomalously strong. This results in an anomalously strong Australian monsoon in January 11 (Fig. 10i), completing the first half of the TBO cycle. In January 11 -January 12 , the mechanism described above is effectively sign-reversed, and an anomalously weak Indian monsoon takes place (Figs. 10j-l), followed by a weak Australian monsoon in the ensuing boreal winter, completing the TBO cycle.
Next, we turn to the TBO composites from 20CRv2c in Fig. 11. There, a biennial oscillation featuring a weak Australian monsoon in January 0 (Figs. 11a-c), followed by a strong Indian monsoon in July 0 (Figs. 11d-f), a strong Australian monsoon in January 1 (Figs. 11g-i), and a weak Australian monsoon in July 1 (Figs. 11j-l) is clearly evident. This pattern is broadly consistent with that in the CCSM4-based composites in Fig. 10, but at the same the two sets of composites have a number of significant differences despite the fact that the 20CRv2c composites have generally weaker statistical significance than their CCSM4 counterparts. Over the Indian Ocean, the main difference between the CCSM4 and 20CRv2c composites is that the former exhibit more dipolar SST anomaly patterns than the latter (see movies 2 and 3 in the supplementary material). This discrepancy is reminiscent of that described in section 3 in the case of ENSO, and is again accompanied by an apparent westward shift of anomalous surface winds associated with the Walker cell in CCSM4. Over the tropical Pacific, the differences are more pronounced, as SST anomalies and the associated surface circulation response are 1808 out of phase between CCSM4 and 20CRv2c (cf. movies 2 and 3 in the supplementary material). FIG. 11. As in Fig. 10, but for the TBO mode families extracted from 20CRv2c SST data. Another visualization of this process can be found in movie 3 in the supplementary material.
In summary, despite notable differences, we find that in both CCSM4 and 20CRv2c the spatiotemporal SST, precipitation, and surface wind patterns represented by the TBO and TBO-A modes are broadly consistent with the mechanism for biennial monsoon variability proposed by Meehl and Arblaster (2002). It is interesting to note that as with ENSO, both the fundamental and combination (TBO-A) modes are important in correctly representing this seasonally locked process, but due to the weak (;0.1 K) SST anomalies involved, these modes are challenging to capture from unprocessed (other than monthly averaged) observational data. To our knowledge, TBO-A modes have not been previously identified in analyses of either model or observational data.

Decadal-interannual interactions and climatic impacts a. Amplitude modulation of ENSO and the TBO
In this section, we study the modulation relationships between the WPMM and other interannual and decadal NLSA modes recovered from model and observational data. Figure 12 displays the temporal pattern corresponding to the WPMM, together with the amplitude envelopes of the fundamental ENSO modes, the fundamental TBO modes, and the IPO mode recovered from the CCSM4 data (see section 4 in Part I). There, one can clearly notice that the envelopes of the fundamental ENSO and TBO modes correlate negatively with the temporal pattern of the WPMM; that is, the interannual modes are strengthened during negative phases of WPMM with negative SST anomalies in the western tropical Pacific. The corresponding pattern correlation (PC) coefficients r are 20.63 and 20.52 (both at numerically zero p values), respectively, for ENSO and the TBO. On the other hand, the correlation between the amplitude of the IPO mode and the WPMM, r 5 0.11 ( p 5 0.2), is not significantly different from zero. A similar lack of statistically significant correlations is observed between the IPO-like modes and the amplitudes of the primary ENSO and TBO modes.
In the case of the modes recovered from 800 years of GFDL CM3 SST data (see section 6 in Part I), we also observe amplitude modulation relationships (not shown here) between the WPMM and the interannual modes, but these relationships are somewhat different than in CCSM4. In particular, we find that in GFDL CM3 the negative phase of the WPMM correlates with increased ENSO and TBO amplitude, but during the positive WPMM phase there is no significant ENSO and TBO 2) in (c). Here, p values were computed using the two-sided Student's t test with (1300 yr) n ENSO 2 2 5 323 degrees of freedom, representing the number of ''independent'' ENSO events at a frequency n ENSO ' 0.25 yr 21 over the 1300-yr CCSM4 dataset. The notation p ' 0 means that p is numerically equal to zero for this number of degrees of freedom. amplitude suppression as observed in CCSM4. As a result, the corresponding time-averaged correlation coefficients are about a factor of 2 smaller than the CCSM4 results in Fig. 12. Another contributing factor in this decrease of correlation is that the temporal patterns of the WPMM and TBO in GFDL CM3 are noisier than in CCSM4. In a separate analysis, we computed modes from an 800-yr portion of the CCSM4 data and observed that the WPMM becomes corrupted by higher-frequency noise, as does the WPMM in GFDL CM3 (see Fig. 11 in Part I). Nevertheless, despite this quality degradation, the correlation between the WPMM and ENSO amplitude (r 5 20.59 with p ' 0) is still higher than in GFDL CM3. This suggests that there is a qualitative difference in the decadal variability between the two models and the corresponding linkages with the interannual modes. Indeed, it is known that in CCSM4 Pacific decadal variability has significant correlations with ENSO (Deser et al. 2012), whereas in GFDL CM3 such relationships are known to be weaker ).
Next, we examine the corresponding correlations in the 20CRv2c analysis. Figure 13 shows the WPMM and IPO time series together with the ENSO and TBO amplitudes from 20CRv2c. Evidently, the results are noisier in this case than in CCSM4, but we still detect a statistically significant correlation r 5 20.53 (p 5 0.001) between the ENSO amplitude and the WPMM. On the other hand, the PC coefficient between the WPMM and the TBO amplitude (r 5 20.20) has the same sign as in the CCSM4 analysis but is not statistically significant (p 5 0.2). The PC coefficient between the WPMM and the IPO amplitude is r 5 20.25 (p 5 0.1), which is higher than the r 5 20.11 value found in CCSM4 (see Fig. 12). As discussed in appendix A, after decadal low-pass filtering, the NLSA IPO mode becomes highly correlated with the IPO mode from EOF analysis of low-pass-filtered Indo-Pacific SST data-the PC coefficient of this ''denoised'' IPO mode with the ENSO amplitude is r 5 20.38 (p 5 0.015), which is again higher than the corresponding value in CCSM4.
In summary, we see that, among the 20CRv2c modes identified here, the WPMM is the decadal mode with the highest correlation to the ENSO amplitude. Of course, it is possible that the diminished correlation compared to CCSM4 is due to intrinsic (i.e., independent of the data analysis algorithm and time series length) differences in the relationships between decadal and interannual variability in CCSM4 and in nature. Nevertheless, the PC results from 20CRv2c provide at least some observational evidence that the WPMM exhibits modulating relationships with ENSO.
To physically account for the amplitude modulation relationships between the WPMM and the interannual modes, we examine its associated reconstructed SST, surface wind, mixed layer depth, and precipitation patterns. As with ENSO and the TBO, we study these patterns through composites, displayed in Figs. 14 and 15 for FIG. 13. As in Fig. 12, but for the WPMM, ENSO, TBO, and IPO modes recovered from the 20CRv2c dataset. The correlation coefficients r between the plotted time series in each panel (and the corresponding p values) are r 5 20.53 ( p 5 0.001) in (a), r 5 20.20 ( p 5 0.2) in (b), and r 5 20.25 ( p 5 0.1) in (c). These p values were computed using the two-sided Student's t test with (162 yr) n ENSO 2 2 ' 38 degrees of freedom, representing the number of ''independent'' ENSO events at a frequency n ENSO ' 0.25 yr 21 over the 162-yr 20CRv2c dataset. the cold WPMM phase from CCSM4 and 20CRv2c, respectively, with negative SST anomalies in the western equatorial Pacific. Unlike our composites for the seasonally locked ENSO and TBO, in this case we did not use a reference month and instead averaged all snapshots in the WPMM reconstructed data where the average SST anomalies in the box 58S-58N, 1508-1708E were less than 21s. The composites corresponding to the positive WPMM phase are qualitatively similar, up to sign reversals, to those for the negative phase, so we do not discuss them here.
In the CCSM4 composites (Fig. 14), it is evident that during the WPMM's cold phase, the decadal cluster of negative SST anomalies in the western Pacific is collocated with surface wind divergence, leading to anomalous westerlies in the central and eastern Pacific and anomalous easterlies in the western Pacific region north of the Maritime Continent. The thermocline depth associated with the cold WPMM phase features strong positive anomalies in the western Pacific region collocated with the anomalous easterlies, but also exhibits appreciable negative and positive anomalies in the central and eastern Pacific, respectively. The 20CRv2c composites in Fig. 15 are in fairly good qualitative agreement with this behavior. There, the cold SST anomaly in the western Pacific is also collocated with anomalous surface wind divergence, although the strength of the surface circulation pattern in the vicinity of that divergence region, and particularly in the western equatorial Pacific, is weaker. Another notable difference is that the positive SST anomalies in the eastern tropical Pacific in the WPMM composites from 20CRv2c are stronger than in CCSM4, leading to a more dipolar SST anomaly pattern over the tropical Pacific. Moreover, while the WPMM thermocline depth patterns from 20CRv2c exhibit negative anomalies in the central Pacific and positive anomalies in the western and eastern Pacific as in CCSM4, the western Pacific anomalies are substantially weaker. It should be noted that the z T20 data used in Fig. 15 come from a different reanalysis product (CHOR_RL) than 20CRv2c, which carries a risk of reconstruction biases, particularly for low-variance patterns such as the WPMM. Nevertheless, despite these differences, in both CCSM4 and 20CRv2c the WPMM is associated with a flatter meridional thermocline profile in the tropical Pacific, and such conditions have been shown to correlate with periods of stronger ENSO variability (Kirtman and Schopf FIG. 14. Composites of (a) SST, (b) XMXL, (c) surface winds, and (d) precipitation over Australia based on the negative (cold) WPMM phase as recovered from CCSM4 data via NLSA. Statistical significance of the composite patterns was assessed via a z test as in Fig. 1. The anomalous westerlies in the central Pacific in (c) and the anomalously shallow (deep) thermocline in the central (eastern) equatorial Pacific in (b) create a background known to correlate with strong ENSO activity. Similarly, the prominent cyclonic circulation in the Indian Ocean and surface wind convergence to the north of Australia in (c), creates an environment favoring strong Australian monsoons, thus strengthening the TBO. 1998; Kleeman et al. 1999;Fedorov and Philander 2000;Rodgers et al. 2004), as observed above. Other studies (Gehne et al. 2014;Karamperidou et al. 2014) indicate that western and central Pacific subsurface and surface oceanic conditions may be a significant source of ENSO predictability.
Over the Indian Ocean, the WPMM composites from CCSM4 and 20CRv2c have more pronounced differences, but nevertheless qualitative similarities can still be found. In CCSM4, the cold WPMM phase features a prominent cyclonic circulation over the Indian Ocean (centered at ;208S, ;908E) that strengthens the western Walker cell and advects warm, moist air from the Maritime Continent toward northern Australia (see Fig. 14c). These conditions are favorable to the Australian monsoon and thus may be related to the observed TBO strengthening during cold WPMM phases. This circulation pattern is also consistent with the formation of a dipole of SST anomalies in the Indian Ocean featuring positive (negative) SST anomalies in the eastern (western) part of the basin (Fig. 14a). In 20CRv2c (Fig. 15c), the cyclonic circulation pattern in the Indian Ocean is weakened, but still visible. Moreover, the dipole of Indian Ocean SST anomalies in CCSM4 is replaced by a monopole pattern (Fig. 15a). The more pronounced differences between the WPMM patterns from CCSM4 and 20CRv2c over the Indian Ocean compared to the Pacific Ocean may contribute to the weakness of the correlation between the WPMM and the TBO amplitude in the latter dataset.

b. Linkages with Australian hydroclimate on seasonal to multidecadal time scales
It is well known that ENSO bears significant influences on many climate modes of variability on seasonal to decadal time scales (Bjerknes 1969;Alexander et al. 2002;Deser et al. 2010). Given that, and the pronounced amplitude modulations of ENSO by the WPMM discussed in section 5a, it appears plausible that the WPMM exhibits linkages to a number of ENSOdependent phenomena. In this section, we demonstrate such a linkage using Australian hydroclimate as an example. Here, we study 1300 years of precipitation as simulated by CCSM4 and integrated over the Australian continent; we denote this variable by P Australia . We also study the corresponding precipitation data from 20CRv2c, with the difference that in this case we replace averages over the Australian continent by averages over rectangular longitude-latitude regions, as we did not have access to a land mask for this dataset.
In agreement with previous studies (Schubert et al. 2016), ENSO explains a significant amount of this region's hydroclimatic fluctuations on seasonal time scales. One way of assessing this contribution, which we will use throughout this section, is through the ratio of the L 1 norms (i.e., the time averages of absolute values of anomalies) of reconstructed and raw precipitation anomaly time series. We refer to this ratio as the ''relative fluctuation strength'' associated with a given mode as it measures the relative average FIG. 15. As in Fig. 14, but for the WPMM recovered from 20CRv2c data. amplitude of fluctuations as opposed to the energy or variance of these fluctuations (captured by the L 2 norm). In CCSM4, the relative strength of P Australia fluctuations associated with the ENSO and ENSO-A modes is approximately 13%, and this number doubles to almost 27% if secondary ENSO and ENSO-A modes are included. The convective and large-scale contributions to this value are about the same (52% in the case of the latter). However, only approximately 6% of the P Australia fluctuations results from ENSO combination modes [this compares well with IOD impact as studied by, e.g., Maher and Sherwood (2014)], split as approximately 3.5% and approximately 15% between the convective and large-scale components, respectively. On average, the monthly anomalies of P Australia associated with ENSO are almost 6 times larger than the ones due to the WPMM. In 20CRv2c, the relative strength of the P Australia fluctuations associated with ENSO and ENSO-A modes is 7% (note that in this case we do not have available a decomposition into convective and large-scale contributions). As in CCSM4, the ENSO modes carry significantly stronger fluctuations than the WPMM-in this case, 2.5 times larger.
Because of the decadal amplitude modulations of the ENSO and ENSO combination modes, the corresponding reconstructed Australian precipitation patterns also exhibit decadal modulations. For example, in CCSM4, ENSO and the corresponding reconstructed precipitation patterns during simulation years 1165-85 (whose corresponding SST patterns are shown in movie 2 in Part I) are stronger than during the century that preceded this period. These two periods correspond to opposite WPMM signs, in agreement with the results of section 5a. In fact, the monthly P Australia anomalies associated with ENSO are, on average, approximately 50% larger during the negative WPMM's phase than during its positive phase. As discussed in section 5a, the correlation between WPMM and ENSO amplitude is weaker in 20CRv2c than in CCSM4, but nevertheless negative (positive) WPMM periods with high (low) ENSO-reconstructed precipitation can still be found.
ENSO's impact on the Australian climate and its decadal variability has been extensively studied in the literature. In particular, the study by Power et al. (1999) introducing the IPO examines precipitation (among other variables) and finds that some decades are characterized by significant ENSO-related variability of the Australian hydroclimate [and thus good performance of forecasting schemes that employ the Southern Oscillation index (SOI)], while others do not exhibit such a relation (with significantly lower skill of SOI predictors).
Notably, they associate these decadal modulations with slow changes of Indo-Pacific SST, and based on that finding derive the index of Indo-Pacific low-frequency variability known as the IPO. Yet, even though Power et al. (1999) establish that the PC coefficient between the IPO index and low-pass-filtered Australian precipitation (based on twentieth-century measurements over the whole continent) can be as high as 0.8, forecasting scores conditioned on the particular phase of the IPO mode do not improve, leaving open questions on the underlying physical mechanism.
In the case of the CCSM4 dataset studied here, the correlation coefficient between the 13-yr running average of P Australia and the IPO index, the latter obtained through the first EOF of the 13-yr running mean of CCSM4 Indo-Pacific SST data (after linear detrending to remove model drift) as in Power et al. (1999), is 0.45 ( p value equal to 1 3 10 27 using the two-sided Student's t test with 1263 degrees of freedom as in Fig. 16). This value is comparable to the 20.5 correlation coefficient between the SOI and low-pass-filtered P Australia reported in the model-based study of Arblaster et al. (2002). Both of these results are indicative of the fact that that in climate models the correlation between the IPO and decadal Australian precipitation tends to be weaker than the value observed in nature.
We now examine the corresponding relationships between decadal precipitation over Australia and ENSO or the WPMM in CCSM4. Following the approach undertaken by Power et al. (1999) and Arblaster et al. (2002), we begin by examining the 11-yr running average of (i) the raw P Australia variable and (ii) the reconstructed P Australia from a collection of primary and secondary NLSA ENSO modes recovered using a 20-yr embedding window (see section 3c) and their associated annual cycle combination modes. Low-pass filtering leads to partial cancellation of precipitation anomalies contributed by different ENSO phases, with residuals associated with El Niño-La Niña asymmetries (captured at least partially by these ENSO modes) contributing to decadal variability of both the raw and ENSOreconstructed P Australia . We find that the relative strength of the low-pass-filtered P Australia fluctuations associated with the ENSO modes is approximately 17.5%, which is lower by about 10% than in the case of the monthly signal. Conditioned on the positive and negative WPMM phases, this number becomes approximately 15% and approximately 20%, respectively. At the same time, we find that the relative fluctuation strength associated with the WPMM is over 45%. In other words, although the WPMM explains only about 2.5% of nonperiodic Indo-Pacific SST variance, its associated strength of decadal, spatially averaged Australian precipitation fluctuations in CCSM4 is 2.5 times stronger than that attributed to the ENSO modes.
Next, we turn to the linkages between decadal P Australia and the WPMM (still in CCSM4), which we examine by means of a precipitation field composite for the negative WPMM phase over Australia (Fig. 14d) and the time series of the WPMM and low-pass-filtered P Australia shown in Fig. 16a. There, it can be seen that the negative phase of the WPMM (i.e., negative SST anomalies over the western tropical Pacific and positive SST anomalies to the south of Australia) are associated with anomalously wet periods that can last almost a century. On the other hand, the positive phase (not shown here) results in prolonged droughts. As with the WPMM temporal pattern, the low-pass-filtered P Australia exhibits pronounced fluctuations on multidecadal time scales, and the correlation coefficient between the two time series is 20.62 ( p ' 0).
Spatially, the negative WPMM phase is associated with precipitation enhancement over the whole Australian continent, with over 53% related to convective rainfall and the remaining 47% attributed to largescale precipitation. This convective intensification occurs predominantly over northern Australia in a dynamically consistent manner with the associated surface atmospheric winds entraining anomalously warm and moist air, which is driven primarily by subsiding circulation and divergent air motion centered over the anomalously cold western Pacific SST. Moreover, additional enhancement of precipitation can be attributed to thermodynamic factors such as positive surface temperature anomalies over Australia, with warm SST surrounding the whole continent. These patterns are sign-reversed (and remain consistent) during the positive WPMM phase.
We now examine the relationship between the WPMM and decadal Australian precipitation in 20CRv2c. Figures 15d and 16b show the analogous precipitation composites and time series to Figs. 14d and 16a, respectively. As in CCSM4, the cold phase of the WPMM recovered from 20CRv2c is seen to correlate with enhanced precipitation over Australia and the surrounding seas to the north of the continent, although in contrast to CCSM4 the precipitation anomalies south of about 358S are generally negative. The corresponding SST (Fig. 14a) and surface circulation (Fig. 15c) patterns are broadly consistent with Australian precipitation enhancement as in CCSM4. In 20CRv2c, the PC coefficient between the WPMM and low-pass-filtered precipitation, averaged over the longitude-latitude box of 1158-1508E FIG. 16. Eigenfunction time series of the WPMM (blue lines) and 11-yr running mean of spatially averaged Australian precipitation (red lines, arbitrary units) for (a) the CCSM4 and (b) the 20CRv2c datasets. In (a) the spatial averaging region consists of the CCSM4 land grid points over the Australian continent. The PC coefficient of the two time series is equal to 20.62 at p value numerically equal to zero. This p value was computed using the twosided Student's t test with (1300 yr 3 0.1 yr 21 ) 2 2 5 128 degrees of freedom, representing the number of ''independent'' decadal precipitation events over the 1300-yr CCSM4 dataset. In (b) the precipitation averaging domain is the longitude-latitude rectangle of 1058-1508E and 208-108S (box 2 in the main text). The PC coefficient of the two time series is equal to r 5 20.53 at p 5 0.001. This p value was computed using the two-sided t Student's test with the same number of degrees of freedom as in Fig. 13; in particular, we did not assign a number of degrees of freedom based on the number of decades in the 20CRv2c dataset as done in the CCSM4 analysis because of the noisier nature of the WPMM time series in 20CRv2c. and 358-108S (hereafter, box 1) encompassing the Australian continent is r 5 20.45 ( p 5 0.004; see the caption to Fig. 16 for degrees of freedom). This relationship is diminished compared to CCSM4, but part of the reduction can be attributed to differences in the spatial precipitation patterns associated with the WPMM. In particular, modifying the averaging region to the more northerly domain 208-108S, 1058-1508E (box 2) allows the PC coefficient to become 20.53 ( p 5 0.001). If, in addition, the WPMM time series is low-pass filtered by decadal running averaging to remove high-frequency noise (see appendix A), these values increase to r 5 20.60 (p 5 0.01 using 15 degrees of freedom to account for low-pass filtering) and r 5 20.67 (p 5 0.003).
Besides pattern correlation, we examine the relative precipitation fluctuation strength associated with the WPMM. For decadal precipitation averaged over boxes 1 and 2, this is equal to 23% and 28%, respectively. By comparison, the relative fluctuation strength associated with the IPO mode (after low-pass filtering to remove interannual variability; see appendix A) is 8.5% and 3.8%, respectively. Interestingly, the PC coefficients between the low-pass-filtered IPO and averaged decadal precipitation over boxes 1 and 2 (in both cases, r 5 0.62 and p 5 0.04) can actually be higher than the corresponding values for the low-pass-filtered WPMM (depending on the averaging region). These results are again indicative of possibly stronger relations between the IPO and Australian precipitation observed in nature than in models.
In summary, in this section we have examined linkages between Indo-Pacific SST and Australian hydroclimate that act across seasonal to decadal time scales. Despite the known deficiencies of global climate models in simulating precipitation Maher and Sherwood 2016), we have seen that the modulating relationships between the WPMM and decadal Australian precipitation are in fairly good qualitative agreement between CCSM4 and 20CRv2c. Of course, as has been pointed out in a number of studies (e.g., Arblaster et al. 2002), the presence of decadal-scale correlations in observational data should be treated with caution as it may be of limited statistical significance. Such a risk was also highlighted in the sensitivity analysis performed in Part I, where we saw that different portions of the 1300-yr dataset of comparable span to industrial-and satellite-era datasets may or may not lead to detectable decadal modes (at least via NLSA operating on unprocessed data). Nevertheless, the physical consistency of the atmospheric and oceanic patterns associated with the WPMM and the level of agreement between the results from the model and observational data studied in this work provide motivation for more detailed mechanistic studies on this mode and its linkages with other aspects of decadal climate variability in models and observations.

Concluding remarks
In this paper, we have employed a family of modes identified via NLSA in Part I of this work, as well as a set of modes recovered here from 20CRv2c SST data, to study various aspects of coupled atmosphere-ocean variability in the Indo-Pacific sector of model and observational datasets on seasonal to decadal time scales. One of our main objectives has been to study the SST, surface wind, and thermocline depth patterns associated with a family of ENSO modes and combination modes between ENSO and the annual cycle (denoted ENSO-A modes), which were recovered via NLSA from monthly averaged, but otherwise unprocessed, Indo-Pacific SST data. We found that in both model data (CCSM4) and observational data (HadISST and 20CRv2c SST for mode extraction; MERRA, C-GLORS, and CHOR_ RL for surface wind and thermocline depth reconstruction) these patterns are in good agreement with mechanisms proposed for the phase synchronization between ENSO and the seasonal cycle (Stein et al. 2011(Stein et al. , 2014Stuecker et al. 2013Stuecker et al. , 2015aMcGregor et al. 2012;Ren et al. 2016), whereby quadratic nonlinearities in the governing equations for the coupled atmosphere-ocean system lead to a seasonally-dependent southward shift of zonal winds following the boreal winter peak of El Niño (La Niña) events, generating an upwelling oceanic response restoring the ENSO state to neutral or La Niña (El Niño) conditions. In previous data-driven studies, the ENSO-A modes have mainly been identified through EOF analysis of atmospheric circulation variables in the tropical Pacific (e.g., McGregor et al. 2012), mixing together the combination tones with the fundamental ENSO frequency into individual modes. Here, we have shown that physically consistent circulation, SST, and thermocline patterns, having the theoretically expected frequencies and mode degeneracies, can be recovered from SST data and on a larger Indo-Pacific domain.
We also showed that the surface wind patterns associated with the ENSO and ENSO-A modes produce a coupling between the Indian Ocean and Pacific SST, with SST patterns resembling the IOD and explaining up to 35% of the raw nonperiodic variance in regions employed in IOD index definition. These results indicate that in both CCSM4 and observations, the ENSO and ENSO-A modes are good predictors for IOD variability. However, unlike the Pacific sector, where the ENSO patterns from CCSM4 and observations were found to be in good qualitative agreement, more significant differences were present in the Indian Ocean. In particular, the reconstructed SST anomalies from CCSM4 exhibit a dipolar pattern during peaking and decaying El Niño and La Niña events, as opposed to the observational reconstructions where dipolar SST anomaly patterns are present in the months leading to El Niño and La Niña events, and become replaced by monopolar patterns during the peak and decay phases. We attributed this discrepancy to a longitudinal difference in the location of the center of the anomalous divergence (convergence) of the Walker cell due to El Niño (La Niña) between CCSM4 and the 20CRv2c and HadISST/C-GLORS data. Given that this discrepancy is consistent among the two observational datasets studied here, it is with some probability the outcome of a bias in CCSM4.
Another aspect of ENSO that we have studied, oftentimes considered a hallmark of nonlinear dynamics (Burgers and Stephenson 1999;An and Jin 2004;DiNezio and Deser 2014), is statistical skewness due to El Niño-La Niña asymmetries. We showed that in CCSM4, Indo-Pacific SST anomalies reconstructed jointly using the leading ENSO and ENSO-A modes from NLSA have negligible skewness in the eastern Pacific cold tongue region, but including a family of secondary ENSO modes present in the NLSA spectrum along with their associated annual-cycle combination modes leads to positively skewed SST anomalies. These secondary ENSO modes have a more rapidly decorrelating, stochastic character and were found be associated with positively skewed surface zonal wind statistics in the tropical western Pacific. Thus, our results are consistent with a number of studies highlighting the role of high-frequency fluctuations (e.g., westerly wind bursts) in producing El Niño-La Niña asymmetries (Fedorov et al. 2003;Rong et al. 2011;Christensen et al. 2017). In the 20CRv2c and HadISST analyses, we used a shorter delay-embedding window (a parameter in NLSA controlling the time-scale separation of the recovered modes) due to the smaller number of available samples, so our results did not exhibit a distinction between primary and secondary ENSO modes as clear as in CCSM4; instead, the leading ENSO and ENSO-A modes directly produced positively skewed statistics in the tropical eastern Pacific. Assuming that the CCSM4 analysis with a longer embedding window is representative of the behavior in nature, our results suggest that ENSO is decomposable into a set of primary ENSO modes with high temporal coherence (and likely predictability) and a set of less predictable secondary modes capturing positively skewed statistics due to El Niño-La Niña asymmetries.
Besides ENSO, we studied the SST, surface wind, and precipitation patterns associated with a family of biennial modes and biennial-annual combination modes representing the TBO in CCSM4 and 20CRv2c. We showed that despite the small SST anomalies involved (about an order of magnitude smaller than typical ENSO amplitudes), these patterns are consistent with coupled atmosphere-ocean mechanisms for biennial variability of the Asian-Australian monsoon , whereby an anomalously weak Australian monsoon in boreal winter tends to be followed by anomalously strong Indian and Australian monsoons in the ensuing boreal summer and winter. At the same time, we found that the TBO patterns from CCSM4 and 20CRv2c exhibit significant discrepancies in the Pacific Ocean, even at the level of anomaly signs. Despite these discrepancies, the facts that our TBO modes were recovered without bandpass filtering the input data, and that they exhibit physically consistent precipitation, SST, and surface wind patterns, support the existence of the TBO as a physically meaningful process. Furthermore, while the TBO is traditionally defined using precipitation variables, our results demonstrate that it can be inferred from SST data alone, although our analysis does not address the causality relationships between it and other modes of climate variability such as ENSO (Goswami 1995;Meehl and Arblaster 2002;Li et al. 2006;Meehl and Arblaster 2011). In Part I, we were less successful in identifying coherent TBO patterns in shorter (satellite era) HadISST data, corroborating similar difficulties faced by other studies (e.g., Stuecker et al. 2015b), but the encouraging results from CCSM4 and industrial-era observational data suggest expanding this study to other datasets and/or including additional physical variables that would improve the detectability of TBO patterns. Notice that, to our knowledge, it is one of the first times that TBO combination modes have been detected from either model or observational data.
On the decadal time scale, we studied amplitude modulation relationships between the WPMM (the dominant mode of Indo-Pacific multidecadal variability recovered by NLSA from CCSM4 data in Part I and 20CRv2c data in this paper, featuring a prominent cluster of SST anomalies in the west-central tropical Pacific) and ENSO, the TBO, and decadal precipitation over Australia. In CCSM4, we found that the WPMM exhibits a significant anticorrelation with the decadal envelope of ENSO and the TBO (at the 20.63 and 20.52 level, respectively), with cold (warm) WPMM phases corresponding to decadal periods of enhanced (suppressed) ENSO and TBO activity. In particular, we showed that cold WPMM periods are associated with anomalous westerlies in the central equatorial Pacific and anomalously flat zonal thermocline profile in the equatorial Pacific, leading to a background configuration that previous studies have shown to favor ENSO activity and predictability (Kirtman and Schopf 1998;Kleeman et al. 1999;Fedorov and Philander 2000;Rodgers et al. 2004). We also demonstrated that cold WPMM phases can favor Australian monsoon activity through anomalous westerlies in the Indian Ocean strengthening the western Walker Cell, which may explain the enhanced TBO activity during those periods. In 20CRv2c, similar modulating relationships were found to take place between the WPMM and ENSO amplitude, although at a weaker (but still statically significant), 20.53 level. The correlation between the WPMM and the TBO amplitude (PC coefficient equal to 20.20) was found to have the same sign as that in CCSM4, but was not deemed to be statistically significant. Moreover, we showed that the WPMM anticorrelates (at PC coefficients as strong as 20.67) with decadal precipitation over Australia in both CCSM4 and 20CRv2c, and we justified this correlation on the basis of advection of warm moist air over the Australian continent by the atmospheric circulation patterns associated with the WPMM. Thus, this mode may be a good predictor of prolonged droughts and/or anomalously wet periods there. Overall, our results allude to the WPMM being an important pattern of Indo-Pacific variability that could potentially be useful in conjunction with the IPO to characterize the decadal behavior of the climate.
In summary, the family of NLSA modes identified in this work provides an attractive framework to characterize numerous aspects of the Indo-Pacific climate system. Potential future applications of these modes include data-driven predictive studies (e.g., via an extension of NLSA for analog prediction; Comeau et al. 2017), as well as analyses of paleoclimatic records as a way of detecting decadal and multidecadal modes in nature (Gergis et al. 2012;Tierney et al. 2015; Vance  (Ghil et al. 2002;Thomson 1982) using a time-bandwidth product p 5 2 and K 5 2p 2 1 5 4 Slepian tapers. The effective frequency resolution is Dn 5 p/(Ndt) ' 1/1857 yr 21 , where N 5 1857 is the number of samples in the eigenfunction time series from the 20CRv2c data. Dashed lines in the power spectral density plots indicate 95% confidence intervals computed using an F test with (2, 2K 2 2) degrees of freedom. Palmer et al. 2015;Dobrowolski et al. 2016;Dixon et al. 2017;Christiansen and Ljungqvist 2017). From a data analysis standpoint, this work has demonstrated that appropriately designed methods taking dynamics into account can potentially overcome some of the limitations of classical linear techniques in revealing coherent climatic processes across multiple time scales, but our study is certainly only an initial step in this direction. Further assessment and refinement of the patterns identified here, potentially via newly developed Koopman operator techniques Slawinska et al. 2017), should be a fruitful avenue for future work. from the Center for Prototype Climate Modeling at NYU Abu Dhabi and NSF EAGER Grant 1551489. D. Giannakis was supported by ONR Grant N00014-14-0150, ONR MURI Grant 25-74200-F7112, and NSF Grant DMS-1521775. This research was partially carried out on the high-performance computing resources at New York University Abu Dhabi. We thank Andrea Storto, Malte Stuecker, and Sulian Thual for stimulating conversations, and the editor and three anonymous reviewers for comments that led to significant improvement of the paper.

SST Modes Derived from 20CRv2c
In this appendix, we describe the NLSA Indo-Pacific SST modes derived from monthly 20CRv2c data. As stated in section 2, we computed these modes for the period January 1851-December 2012 and the same Indo-Pacific domain as in Part I without detrending the input data. We also used the same kernel and NLSA parameters as in the analyses from Part I, aside from the embedding window length Dt, which was set to 8 yr (as opposed to the Dt 5 4 yr window used in the HadISST analysis), as we found that this value yields a cleaner representation of the TBO and the WPMM. With these parameters, NLSA applied to 20CRv2c data recovers a family of modes (presented in the same order as the CCSM4-derived modes in Part I) that includes the annual cycle and its harmonics, ENSO, ENSO-A combination modes, the TBO, a TBO-A combination mode, the WPMM, and an IPO-like mode. Thus, we recover the same classes of modes as in CCSM4, although as we will see below the TBO, TBO-A, WPMM, and IPO modes exhibit various degrees of frequency mixing. Notice that, at the same time, TBO-A modes have not FIG. B1. SST, surface wind, and z T20 anomalies for the El Niño event of 1997/98 reconstructed using the ENSO and ENSO-A modes recovered from satellite-era HadISST data via NLSA. Statistical significance of the composite patterns was assessed via a z test as in Fig. 1. been extracted successfully in other studies. Nevertheless, these patterns are substantially cleaner than their counterparts from the industrial-era HadISST analysis in Part I (see Fig. 14 in that paper). The NLSA spectrum from 20CRv2c also contains trend-like modes, as well as modulations of the trend-like modes by the annual and interannual modes; we defer study of these modes to future work. Figure A1 shows eigenfunction time series and power spectral densities for the fundamental ENSO modes (eigenfunctions ff 15 , f 16 g), the leading two combination modes between ENSO and the annual cycle (eigenfunctions ff 18 , f 19 g and ff 26 , f 27 g, respectively), a TBO mode (eigenfunction f 42 ), a TBO-A-like mode (eigenfunction f 50 ), the WPMM (eigenfunction f 49 ), and an IPO-like mode (eigenfunction f 20 ) computed from the 20CRv2c data. The ENSO and ENSO-A modes form degenerate in-quadrature pairs; only one member from each pair is plotted in Fig. A1. The fundamental ENSO modes (Figs. A1a,b) exhibit a pronounced, yet broad, spectral peak at the frequency n ENSO ' 0.25 yr 21 , while also exhibiting lower-frequency (decadal) amplitude modulations. In the modes in Figs. A1c,d and A1e,f, that frequency is shifted to 1 yr 21 2 n ENSO and 1 yr 21 1 n ENSO , respectively, in a consistent manner with combination mode theory. Moreover, the amplitude envelope of these modes correlates well with the fundamental ENSO envelope (PC coefficients equal to 0.75 in both cases, with p ' 5 3 10 26 based on the two-sided Student's t test with 26 degrees of freedom as in Fig. 13). In Part I the modes in Figs. A1c,d and A1e,f were denoted ENSO-A1 and ENSO-A2, respectively.
Meanwhile, eigenfunction f 42 (Figs. A1g,h) exhibits a discernible spectral peak at the fundamental TBO frequency, n TBO ' 0.5 yr 21 , and eigenfunction f 50 (Figs. A1i,j) exhibits spectral peaks at n TBO and the TBO-A frequency of 1 yr 21 1 n TBO . While both of these frequencies are important for capturing the TBO and its phase locking to the annual cycle (see section 4), it should be noted that these TBO modes have two discrepancies from their theoretically expected properties, namely, (i) they do not appear as doubly degenerate, in-quadrature pairs and (ii) eigenfunction f 50 mixes together the fundamental TBO and TBO-A frequencies. These discrepancies may be caused by the smaller number of samples and embedding window used in the 20CRv2c analysis compared to the CCSM4 analysis, where we recover distinct, twofold degenerate TBO and TBO-A modes. [Note that in EOF analysis, frequency mixing between the fundamental and combination frequencies occurs even for high-variance modes such as ENSO (McGregor et al. 2012;Stuecker et al. 2013;Ren et al. 2016).] It should also be noted that the TBO-A time series in Figs. A1i,j also exhibits a moderate amount of mixing with decadal frequencies.
As demonstrated in section 4, despite these issues, these modes provide an adequate representation of the TBO life cycle. As in the model data, the WPMM eigenfunction f 49 (Figs. A1k,l) carries significant power on decadal time scales, although it exhibits appreciable power spectral density (peaking at ;0.15 yr 21 ) at higher frequencies than the WPMM from CCSM4 (whose power spectral density diminishes above around 0.03 yr 21 ; see Fig. 5a in Part I).
The WPMM from 20CRv2c also exhibits a moderate amount of noise at frequencies greater than 1 yr 21 . Nevertheless, as with the TBO modes, the 20CRv2c-based WPMM yields qualitatively similar spatial composites to its CCSM4 counterpart (see section 5).
The IPO-like eigenfunction f 20 (Figs. A1m,n) displays a significant amount of mixing with interannual frequencies, although its associated SST pattern (not shown here) features the wedge-shaped eastern Pacific anomaly region characteristic of the IPO. Despite this mixing, the f 20 time series is found to correlate significantly (PC 5 0.74 at p ' 0) with the leading nontrend principal component of decadal Indo-Pacific 20CRv2c SST data, which is in turn closely related to the conventional IPO definition through EOF analysis (Power et al. 1999). Moreover, after application of a 10-yr low-pass filter, this correlation increases to PC 5 0.96 ( p ' 0). For these reasons, we interpret f 20 as an IPO mode despite its prominent interannual-frequency content.

ENSO Reconstructions from Satellite-Era HadISST Modes
In this appendix, we discuss spatiotemporal reconstructions of two historical El Niño and La Niña events based on the satellite-era HadISST ENSO and ENSO-A modes. As stated in section 2, we complement our SST reconstructions with MERRA surface wind data and C-GLORS z T20 data. Figure B1 shows monthly reconstructed snapshots, plotted every 4 months, of the period March 1997-November 1998 encompassing the strong 1997/98 El Niño event. This event is also visualized in movie 4 in the supplementary material, which spans the period March 1979-October 2009.
In Fig. B1, the evolution is broadly consistent with the 20CRv2c-derived El Niño composites in Fig. 3, exhibiting a buildup, peak, and decay of El Niño SST anomalies, along with the corresponding surface circulation and thermocline depth patterns. Notice in particular the anomalous easterlies in Figs. B1e,h,k and east-west thermocline gradient in Figs. B1f,i,l that develop as El Niño matures. Starting in late fall of 1997, and during the course of El Niño dissipation, the anomalous easterlies shift to the south and eventually diminish (Figs. B1h,k,n). Similarly to the 20CRv2c composites, the HadISST reconstructions exhibit significant qualitative differences in the Indian Ocean from the CCSM4-based El Niño composites (Fig. 1) and reconstructions (movie 1 in the supplementary material), in that they feature a positive monopole (as opposed to a dipole in CCSM4) of Indian Ocean SST anomalies in the months during and after the El Niño peak (e.g., March 1998), while a clearer dipole is present in the months leading to the peak (e.g., July and November 1997). Figure B2 shows analogous reconstructed snapshots for the 1999/2000 La Niña event, which is also visualized in movie 4 in the supplementary material. There, the behavior is again broadly consistent with the 20CRv2c-based La Niña composites in Fig. 4.