Disentangling Global Warming, Multidecadal Variability, and El Niño in Pacific Temperatures

1 Introduction

Internal variability of Pacific sea surface temperatures (SSTs) and ocean heat content has been implicated in the decadal modulation of global mean surface temperature, in particular during the global warming “hiatus” between 1998 and 2013 (England et al., 2014; Kosaka & Xie, 2013; Meehl et al., 2011; Trenberth & Fasullo, 2013). However, long-term predictions of Pacific SSTs are limited by our understanding of the interactions between El Niño–Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO) and the role both play in the decadal modulation of global warming. ENSO and the PDO, as traditionally defined, have similar spatial patterns despite differences in time scale (Mantua et al., 1997; Zhang et al., 1997), complicating diagnosis of their mechanisms and interactions.

The PDO is the dominant mode of SST variability in the North Pacific, traditionally defined as the first principal component of Pacific SST north of 20°N (Mantua et al., 1997). Its positive phase exhibits warm anomalies that extend along the west coast of North America and into the tropical Pacific, concurrent with cold anomalies extending from the coast of Japan within the Kuroshio-Oyashio extension. The PDO exhibits variability on interannual to multidecadal time scales and is characterized by rapid phase shifts (Mantua et al., 1997; Minobe, 1997). It has been suggested that PDO modulates the atmospheric response to El Niño (Gershunov & Barnett, 1998; McCabe & Dettinger, 1999; Wang et al., 2008), influencing the variability of precipitation in East Asia and western North America. However, because PDO as traditionally defined is strongly correlated with ENSO, some studies have suggested that it is primarily a high-latitude response to El Niño (Newman et al., 2003) or an average over ENSO precursors and postcursors (Vimont, 2005). Climate models show a more limited connection between PDO and tropical Pacific SSTs (Newman et al., 2016), instead emphasizing midlatitude gyre dynamics (Latif & Barnett, 1994; Schneider et al., 2002; Zhang & Delworth, 2015). The emerging picture is that PDO is not one process but many, with each process having distinct mechanisms and time scales (Newman et al., 2003, 2016; Schneider et al., 2002). To study the links between PDO, ENSO, and global warming, it is useful to define PDO in a way that does not mix it with other modes of climate variability.

More generally, studies of climate variability and change rely on reducing the dimensionality of spatially and temporally varying fields. The dimension reduction often takes one of two approaches: reducing the temporal complexity by linear trend analysis or reducing the spatial complexity by analyzing the time evolution of large-scale averages, such as a global average. The former presupposes a response that varies linearly in time; the latter loses valuable information about the spatial structure of changes. Principal component analysis (PCA) takes advantage of the spatial structure of covariation in climate data to find a basis of Empirical Orthogonal Functions (EOFs) that are ordered by the fraction of total variance they capture. The associated principal components give the evolution of each EOF in time, such that PCA includes information about both spatial and temporal variations. However, because PCA maximizes the variance captured by the first EOF, it can group together multiple processes and suggest spurious connections that are not rooted in shared physical mechanisms (Deser, 2000), as is apparent in the mixture of global warming and ENSO found in the first EOF of monthly Pacific SST anomalies (Figure S1 in the supporting information). Pairwise rotation of principal components is effective at eliminating mode mixing in some cases and has been used to define uncorrelated indices of PDO, ENSO, and global warming (Chen & Wallace, 2016; Chen et al., 2017). In this work, we present a method that transforms any number of principal components at once, based on minimization of their low-frequency covariance, in order to eliminate mode mixing and isolate modes of low-frequency variability.

We identify modes of low-frequency variability by solving for linear combinations of EOFs that maximize the ratio of low-frequency to total variance, using a method adapted from pattern recognition (Ripley, 1996; Schneider & Held, 2001). We define low-frequency variance as the variance remaining after the pointwise application of a Lanczos lowpass filter with a 10 year lowpass cutoff. By considering high-frequency as well as low-frequency variance, this method takes advantage of the spatial structure of covariance in high-frequency variability to optimally filter it out. Our method—low-frequency component analysis (LFCA)—identifies low-frequency patterns (LFPs) and corresponding uncorrelated low-frequency components (LFCs) that are ordered by their ratio of low-frequency to total variance, effectively sorting modes of variability based on their dominant time scale (see Methods, Section 5). In this way, LFCA isolates modes of low-frequency variability without obscuring fast variability that may contribute to their evolution. We apply this method to study monthly Pacific SST anomalies between 45°S and 70°N over the period 1900–2016 using data from the National Oceanic and Atmospheric Administration Extended Reconstructed Sea Surface Temperature data set version 3b (Smith et al., 2008).

2 Disentangling Global Warming, PDO, and ENSO

We first demonstrate LFCA on the simple case where only three EOFs are included, before illustrating the utility of retaining further EOFs in the next section. The first three EOFs of Pacific SST show a mixture of long-term warming, PDO, and ENSO (Table S1 and Figure S1). LFCA finds a linear combination of these three EOFs that eliminates their covariance on decadal time scales, helping to correct for this mode mixing. The resulting three LFPs/LFCs show a clean separation of SST variability into components associated with long-term global warming, PDO, and ENSO (Figure 1). Since only three EOFs are included, the third LFP/LFC shows the variability that remains after the low-frequency variability is filtered out, illustrating that this method is filtering out interannual ENSO variability.

The first LFP shows a nearly uniform warming pattern, corresponding to a 0.67°C per century warming of Pacific SSTs throughout the record (Figure 1a). The warming of tropical Pacific SSTs is notably uniform, without a clear El-Niño-like or La-Niña-like warming trend. There is detailed temporal structure to the global temperature rise, with periods of reduced global warming between the mid-1940s and late 1980s and between 1998 and 2013. This analysis gives a clearer picture of the long-term increase in Pacific SSTs through the twentieth century than could be obtained with linear trend analysis or regression against the basin-mean temperature change. In particular, the ratio r of low frequency (greater than 10 years) to total variance is 8% greater than for a Pacific basin mean (Table S1), demonstrating that LFCA is a more effective temporal filter than a large-scale average.

The second LFP represents the PDO (Figure 1b); its LFC is strongly correlated (94%) with the traditional definition of PDO as the first principal component of North Pacific SST (Mantua et al., 1997), even though it is obtained here from an analysis of SST anomalies over the full Pacific basin (45°S to 70°N). The associated SST variability is strongest in the Northern Hemisphere midlatitudes and is weaker in the tropical Pacific and the Southern Hemisphere. The third LFP represents ENSO, the dominant mode of SST variability left over after global warming and PDO variability are removed (Figure 1c). It is highly correlated with traditional indices of ENSO such as Niño3.4 (84% correlation), Niño3 (89% correlation), and the cold tongue index (Deser & Wallace, 1987) (86% correlation). The corresponding SST pattern is confined to the eastern equatorial Pacific to an even greater extent than SST regressions on traditional ENSO indices. The PDO and ENSO indices obtained from LFCA are uncorrelated by construction, in contrast to the 50% correlation between the traditional PDO index (Mantua et al., 1997) and Niño3.4. Despite this construction, our indices are largely consistent with the canonical definitions in terms of spatial pattern and time evolution, providing a new framework for understanding the relationship between PDO and ENSO. We refer to these LFCA-based PDO and ENSO indices as PDO* and ENSO* throughout the rest of the text.

The power spectra of PDO* and ENSO* are marginally more separated from each other than the power spectra of the traditional PDO and Niño3.4 (Figure 2a). PDO* has 10% more power at periods greater than 10 years than the traditional PDO. ENSO* has more power at periods of 3–4 years than Niño3.4, and less at periods between 8 and 30 years, suggesting that some of the multidecadal variability of ENSO has been filtered out and is actually associated with PDO. While most of these power spectra differences are not statistically significant, due to the short length of the record, these results suggest that the PDO* and ENSO* indices are as good as or better than the traditional indices (i.e., they are more representative of the time scales typically associated with PDO and ENSO). Time-lagged correlations between PDO* and ENSO* differ substantially from time-lagged correlations between the traditional PDO index and Niño3.4, with a weaker correlation when ENSO leads and a stronger correlation when PDO leads by more than 1 year (Figure 2b). By removing the correlation between PDO and ENSO and finding that the resulting PDO* and ENSO* indices are broadly consistent with canonical definitions, LFCA suggests that the mechanistic connection between PDO and ENSO may be weaker than previously thought. Moreover, ENSO* exhibits a much weaker relationship than Niño3.4 with the atmospheric circulation over the North Pacific, as characterized by the strength of the Aleutian low in the Twentieth Century Reanalysis (Figure S2) (Compo et al., 2011). This factor of 3 reduction in the ENSO teleconnection to the Aleutian low may help bring observational estimates in line with model-based estimates (Sterl et al., 2007).

3 Isolating Multidecadal Variability

By increasing the number of EOFs retained in the analysis, we include variability at smaller spatial scales that contains significant low-frequency power, further isolating patterns of multidecadal variability and filtering out variability that acts at shorter time scales. The key difference between our method and that of Chen et al. (2017) is that it allows for the case where small-scale processes captured by higher-order EOFs (such as the dynamics of ocean frontal zones) are important for low-frequency variability.

The pattern of long-term warming (LFP 1) is nearly unchanged across a wide range of EOF truncations, suggesting that we have converged on a good approximation of the true greenhouse gas forced temperature response. As more EOFs are retained, high-frequency contributions to PDO* are filtered out of the most PDO-like mode (determined based on its temporal correlation with PDO*), leading to an increase in its ratio of low-frequency to total variance (Figure 3a) and a decrease in the fraction of the total variance associated with this mode (Figure 3b). Note, however, that because the spatial patterns (LFPs) are not orthogonal, the fraction of variance in the modes need not add to the total explained variance. The PDO-like mode maintains a high temporal correlation with PDO* across a range of LFCA parameters (Figure 3c) and remains nearly independent of ENSO* as long as the lowpass cutoff is greater than 6 years (Figure 3d).

With 30 EOFs retained, the most PDO-like mode (LFC 4, based on its temporal correlation with PDO*: 68%) is contained within the midlatitude North Pacific (Figure 4d). It has more than four times the persistence of the traditional PDO (assessed at an autocorrelation of 0.2, Figure 2c). It shows 0.4°C warming per standard deviation along the west coast of North America, where PDO was first identified (Mantua et al., 1997). This is comparable to the 0.5°C warming associated with the traditional PDO index in the same region (see side-by-side comparison in Figure S3). Its spatial pattern differs from the traditional PDO definition primarily in the Kuroshio extension and the equatorial Pacific. Notably, the PDO-like mode (LFC 4) captures the major PDO phase shifts in 1924/1925, 1947/1948, 1976/1977, and 1998/1999, and 2013/2014 (Figure 4d). This mode resembles the response of the North Pacific subpolar and subtropical gyres to wind stress forcing, which has been discussed in the context of the 1976/77 PDO shift (Seager et al., 2001) and the PDO in coupled climate models (Schneider et al., 2002; Zhang & Delworth, 2015). Its time-lagged correlation with ENSO* is less than 0.1 at all lags except at a lead time of 2.5 years (Figure 2b), implying that this mode is mostly independent of ENSO, with possibly some (weak) ability to influence ENSO at 2.5 year lead times.

The other leading LFCs in the 30 EOF analysis describe global warming, coherent multidecadal variations of SST in the Kuroshio extension and eastern equatorial Pacific, and decadal variability linked to the central equatorial Pacific (Figure 4). The long-term warming mode (LFC 1, Figure 4a) is almost identical to that obtained in Figure 1 (temporal correlation of 97%), suggesting that our analysis robustly identifies the influence of global warming on the Pacific. The multidecadal mode (LFC 2, Figure 4b) shows a pronounced La-Niña-like warming trend since 1977, as has been observed (Hansen et al., 2006; Karnauskas et al., 2009; Kohyama et al., 2017). It may in part be the result of forcing from anthropogenic aerosols (Takahashi & Watanabe, 2016). Interestingly, this mode explains much of the multidecadal variability of ENSO; large El Niño events (greater than 2σ anomaly in ENSO* for three consecutive months) did not occur during the period from 1920 to 1970, when LFC 2 was in a strong positive phase. It does not contribute substantially to PDO* variability (correlation of −7%). The central equatorial Pacific mode (LFC 3, Figure 4c) shows variability on 10–20 year time scales. It contributes to the variability of PDO* (correlation of 33%) but less than the North Pacific PDO-like mode (LFC 4; correlation of 68%). ENSO is filtered out of the leading LFCs and is split across a few of the remaining (higher frequency) LFCs. In this analysis, variability characterized by the traditional PDO is split across multiple LFCs; however, the PDO itself is thought to be the superposition of several different physical mechanisms (Newman et al., 2016), and there is no physical motivation for it to be characterized by a single statistical mode.

Across a range of LFCA parameters, the most PDO-like mode illustrates the properties that are essential to the PDO: warm SST anomalies along the west coast of North America, cold anomalies extending from Japan eastward within the Kuroshio-Oyashio extension, variability on 10–40 year time scales, and marked phase transitions (Figures 4d and S4). The most tenuous feature of the PDO-like mode is its connection to tropical Pacific SSTs, which is strongest when fewer EOFs are included in the analysis or when a shorter lowpass cutoff time scale is used (6 years or less). Where there is a strong connection between the PDO and equatorial SSTs, such as in the traditional PDO definition, this connection comes not from decadal variability but from averaging in the effects of a few large El Niño years (Figure S5). Models that show a weaker connection between the PDO and the tropics (Newman et al., 2016) are not necessarily wrong, as the apparent connection to the tropics in observations comes mostly from a few large events. This does not preclude a causal relationship between equatorial SSTs and the PDO, but it does suggest that the warm phase of the PDO can occur irregardless of the state of the eastern equatorial Pacific.

4 Discussion and Conclusions

We have demonstrated a means of disentangling the influences of global warming, multidecadal internal variability, and ENSO on Pacific SSTs, using a statistical methodology (LFCA) that separates modes of climate variability and change based on time scale. Our work quantifies the long-term warming of the Pacific that is independent of the dominant modes of internal variability, without relying on climate models, strengthening existing evidence for the role of global warming in shaping twentieth century temperatures. Furthermore, we characterize the modes of internal variability that act on decadal time scales, recovering an index of PDO variability from an analysis of SST anomalies in the full Pacific basin and showing that its evolution is largely independent of ENSO.

The traditional definition of the PDO (Mantua et al., 1997) aimed to capture coherent decadal SST variability along the west coast of North America, with impacts on salmon production. PDO-like modes based on LFCA capture this SST variability, as shown by their similar spatial pattern in the North Pacific and their coherence with the traditional PDO index at decadal time scales (Figure S6), yet they have a higher ratio of low-frequency (decadal) to total variance, implying increased persistence (Figure 2c). The PDO* and ENSO* indices (Figure 1) provide useful definitions of PDO and ENSO for impact studies, because they provide uncorrelated indices of PDO, ENSO, and long-term Pacific warming. PDO* and ENSO* are similar to canonical definitions of PDO and ENSO despite the fact that they are uncorrelated by construction. This suggests that PDO and ENSO may not be as mechanistically linked as has been thought previously. By retaining more EOFs in the LFCA, we isolate the decadal component of PDO variability and show that it is confined to midlatitudes. This suggests a path forward for long lead-time predictions of Pacific temperature: future work should focus on SST variability in midlatitudes and the impact of sharp SST fronts, such as the Kuroshio-Oyashio extension, on atmospheric circulations.

The two main examples presented in this study represent opposite end-members of the LFCA. For a small number of EOFs (as in Figure 1), LFCA rearranges the variability contained in the included EOFs based on the ratio of low-frequency to total variance and is useful for fixing problems with mode mixing. When more EOFs are included (as in Figure 4), LFCA converges on a few modes that describe the low-frequency variability in a data set and are relatively insensitive to further parameter changes (Figure S7). Our LFCA method is more straightforward and more powerful than existing approaches to filter climate variability based on time scale and provides a framework for separating modes of climate variability in other contexts.

5 Methods

5.1 Low-Frequency Component Analysis

Our method is based on linear discriminant analysis, a statistical method widely used in pattern recognition and machine learning to find linear combinations of properties that best separate groups of data by maximizing the intergroup variance relative to the total variance (Ripley, 1996). Schneider and Held (2001) suggest an application of linear discriminant analysis to spatiotemporal data that isolates modes of low-frequency variability by identifying linear combinations of the first N EOFs that maximize the ratio of low-frequency to total variance. While they focused on isolating externally forced climate change from internal variability, their method also provides a framework for characterizing internal variability on different time scales, since the resulting modes tend to be ordered by time scale. For example, their method has been used to characterize the competing influences of the 11 year solar cycle, the Quasi-Biennial Oscillation, and ENSO on the polar stratosphere (Camp & Tung, 2007a, 2007b) and the influence of the solar cycle on surface temperatures (Tung & Camp, 2008), but it has not been widely used elsewhere. We simplify Schneider and Held's (2001) methodology by using a lowpass filter to define the low-frequency variance, in place of the decadal grouping of data that they used. We call this method low-frequency component analysis (LFCA). In LFCA, the intergroup variance to be maximized is generalized to the variance remaining after application of a lowpass filter. This is related to a broader class of spectral discriminant or optimal filtering analyses, based on the maximization of a particular type of variance representing a “signal” compared to “noise” that exists within internal variability or among realizations (Allen & Smith, 1997; Chang et al., 2000; Schneider & Griffies, 1999; Ting et al., 2009; Venzke et al., 1999).

LFCA starts with a conventional PCA to identify a small number of spatial degrees of freedom that explain much of the variance in the data. For an n × p spatiotemporal data matrix X with zero time mean (e.g., n months of SST anomalies at p grid points), we compute the EOFs, which are the eigenvectors a_k of the sample covariance matrix

(1)

such that . Note that for our analysis of monthly SST anomalies we have subtracted the climatological seasonal cycle. The data matrix X is weighted by the square root of grid cell area such that the covariance is area weighted. The EOFs a_k are normalized such that ||a_k||=1, which means the eigenvalue gives the variance associated with the kth EOF. The total variance of the data set is . The projection of the kth EOF onto the data matrix X defines the kth principal component, PC , where the factor ensures that the principal component has unit variance.

LFCA is based on the projection of filtered data onto a truncated basis of EOFs. Here we apply a linear Lanczos lowpass filter L(T) with cutoff frequency T⁻¹ and periodic boundary conditions to the data, obtaining a lowpass filtered data matrix . The projection of the lowpass filtered data onto the kth EOF, , is equivalent to lowpass filtering the kth principal component because the filter is linear. We look for linear combinations

(2)

of the first N EOFS such that the ratio r_k of low-frequency to total variance is maximized when the data are projected onto them:

(3)

The coefficient vectors e_k are normalized such that ||e_k||=1. The normalization factors in 2 serve to give the linear combinations u_k unit variance, such that the covariance in the denominator of 3 is equal to 1. Using 2-3, and the definition of a principal component, we find that the coefficient vectors e_k are eigenvectors of the covariance (cov) matrix of the first N lowpass filtered principal components,

(4)

The matrix R has N eigenvectors, Re_k=r_ke_k, which we sort by eigenvalue r_k, the fraction of the total variance in the kth mode that occurs at low frequencies (i.e., makes it through the lowpass filter). By projecting the unfiltered data onto the linear combination vectors u_k, we find the low-frequency components (LFCs):

(5)

The LFCs are dominated by low-frequency variability but can still exhibit rapid transitions. The regression of the unfiltered data onto the kth LFC,

(6)

defines the spatial pattern of the kth mode of low-frequency variability, which we call the kth low-frequency pattern (LFP).

This methodology is a transformation of the first N EOFs such that the resulting modes are ranked by the fraction of their variance that occurs at low frequencies (based on lowpass cutoff time scale T). The transformation relaxes the constraint on orthogonality of the LFPs v_k but requires the LFCs to be uncorrelated, whereas orthogonal EOF rotations, such as Varimax Rotation (Kaiser, 1958; Richman, 1986), maintain the orthogonality of the EOFs but relax the orthogonality of the principal components. LFCA can be generalized to maximize the variance in any range of frequencies. In the case of lowpass filtering used here, there are only two parameters: the lowpass cutoff T and the EOF truncation number N. We use a 10 year lowpass cutoff to discriminate between modes of decadal variability and modes of interannual variability (such as ENSO). The sensitivity of our results to lowpass cutoff T and EOF truncation number N is explored in Figures 3 and S4. In general, our method also exhibits some sensitivity to the choice of filter boundary conditions. Here we use periodic boundary conditions for simplicity, but we can recover similar results by filtering only departures from the linear trend and using reflected boundary conditions. This latter method may be useful in data sets with strong linear trends.

The three LFCs resulting from the transformation of the first three EOFs of Pacific SST describe SST variability associated with global warming, PDO, and ENSO. All choices of lowpass cutoff greater than 5 years give indices that are greater than 99% correlated with the indices shown in Figure 1 (based on a 10 year cutoff). The ordering can be switched by replacing the lowpass filter with a highpass or bandpass filter. The high-frequency components are indistinguishable from the low-frequency components (temporal correlations greater than 99.9% for the Lanczos filter used in this study), though the order is reversed. This suggests a broader class of spectral discriminant analysis, where other types of filters are used in place of the lowpass filter used in LFCA. This could be useful for characterizing higher-frequency modes of climate variability.

LFCA is closely related to optimally persistent pattern (OPP) analysis, another time series analysis method based on linear discriminant analysis (DelSole, 2001, 2006; DelSole et al., 2011), which identifies modes of variability with maximal decorrelation time. This is accomplished by using linear combination coefficients e_k that are eigenvectors of the matrix of principal component decorrelation times

(7)

The main difference between OPP analysis and LFCA is the presence of a specified filter time scale T in LFCA, compared to the integration over all time scales in OPP. We have repeated our analysis with the OPP method and found that it primarily optimizes interannual persistence, resulting in modes of variability with lower ratios of low frequency (decadal) to total variance than LFCA, especially when more EOFs are retained. This suggests that LFCA is a better tool for filtering out ENSO and studying decadal variability, while OPP is useful in other contexts. Both LFCA and OPP provide substantially more effective spatiotemporal filtering than the common practice of applying PCA to filtered data.

Erratum

The supporting information has been revised since this article was originally published. The colorbar in Figure S3 and Figure S5 has been corrected. This corrected supporting information may be considered the authoritative version of record.