## Abstract

Changes of the Pacific decadal oscillation (PDO) under global warming are investigated by using outputs from phase 5 of the Coupled Model Intercomparison Project (CMIP5) and a theoretical midlatitude air–sea coupled model. In a warming climate, the decadal variability of the PDO is found to be significantly suppressed, with the amplitude reduced and the decadal cycle shifted toward a higher-frequency band. We used the theoretical model put forward by Goodman and Marshall (herein the GM model) to underpin the potential mechanisms. The GM model exhibits a growing coupled mode that resembles the simulated PDO. It is found that the suppression of the PDO appears to be associated with the acceleration of Rossby waves due to the enhanced oceanic stratification under global warming. This shortens the PDO period and reduces PDO amplitude by limiting the growth time of the coupled mode. The GM model also suggests that the increase of growth rate due to strengthening of oceanic stratification tends to magnify the PDO amplitude, counteracting the Rossby wave effect. This growth rate influence, however, plays a secondary role.

## 1. Introduction

As the dominant variability of the North Pacific sea surface temperature (SST) on the decadal time scale, the Pacific decadal oscillation (PDO) has received much attention since its identification in the late 1990s (Mantua et al. 1997; Newman et al. 2016). Commonly defined as the leading empirical orthogonal function (EOF) mode of SST anomalies (SSTA) in the North Pacific, the PDO displays a horseshoe-like pattern, with temperature anomalies of one sign in the central and western North Pacific encircled by anomalies of the opposite sign along the North American coastline and toward the tropics (Mantua et al. 1997). Both observational and modeling studies indicate that the PDO displays substantial decadal to multidecadal time-scale variability, with characteristic periods of approximately 15–25 and 50–70 years (Latif and Barnett 1994; Minobe 1997; Minobe 1999; Wu et al. 2003; Kwon and Deser 2007; Zhong and Liu 2009). The PDO has shifted its phase several times during the past century with far-reaching impacts on the regional marine ecosystem (Mantua et al. 1997; Mantua and Hare 2002), the climate over the surrounding continents (Cayan et al. 1998; Krishnan and Sugi 2003; Yu et al. 2015), El Niño–Southern Oscillation (ENSO; Gershunov and Barnett 1998), and the global mean temperature (Kosaka and Xie 2016).

Many studies have attempted to understand the origins of the PDO. One school of thought suggested that SST decadal variability arises from oceanic integrating effects on large-scale atmospheric stochastic forcing (e.g., Frankignoul and Hasselmann 1977; Jin 1997; Pierce et al. 2001). Other studies argued that ocean–atmosphere coupling over the midlatitudes serves as another candidate mechanism for the decadal variability of the PDO (Wu et al. 2003; Liu and Wu 2004; Fang et al. 2014) and the North Atlantic tripole (Wu and Liu 2005; Yang et al. 2012). Involving the ocean–atmosphere coupling or not, both of the above mechanisms suggest that the oceanic Rossby wave driven by the atmospheric forcing plays a crucial role in maintaining and setting the time scale of PDO decadal variability (Miller et al. 1998; Seager et al. 2001; Schneider et al. 2002). In particular, Fang et al. (2014) demonstrated that the strong decadal peak of the PDO spectrum disappears when blocking the oceanic Rossby wave.

Given the enormous impacts of the PDO on global climate and ecosystems, it is important to investigate its change in a warming climate and the underlying mechanisms. Using outputs from the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4), Furtado et al. (2011) examined the modulation of the PDO in response to global warming and found insignificant changes in either spatial pattern or temporal characteristics. However, the model data they used are too short for decadal study. Yang et al. (2012) first suggested that the decadal variability in the North Atlantic moves toward higher frequency, due to the acceleration of oceanic Rossby wave in response to the enhanced oceanic stratification in a warmer climate. This Rossby wave mechanism is further supported by the subsequent studies of PDO (Fang et al. 2014; Zhang and Delworth 2016) and the Atlantic meridional overturning circulation (AMOC; Cheng et al. 2016). Besides the time scale, the amplitude of the decadal variability is also found to be weakened in global warming (Fang et al. 2014; Zhang and Delworth 2016; Cheng et al. 2016). Tree ring proxy data also indicate less pronounced PDO amplitude after the industrial revolution (D’Arrigo et al. 2001), providing evidence for model simulations.

So far climate models appear to provide some evidence for the changes of decadal climate variability in warming climate, but the mechanisms remain elusive, given the complexities of the PDO dynamics. Fang et al. (2014) identified this change of the PDO amplitude but did not explain the underlying mechanisms. Zhang and Delworth (2016) stressed that a reduced SST meridional gradient is responsible for the weakened amplitude. Therefore, it is necessary to build a theoretical framework to understand the underlying mechanisms for such changes. In this paper, we attempt to explain the PDO modulations in a warming climate based on the simple coupled model put forward by Goodman and Marshall (1999, hereinafter GM99) and the climate models.

The paper is structured as follows. Section 2 begins with a brief description of the data and the GM model. Sections 3 and 4 present the PDO simulated in models in phase 5 of the Coupled Model Intercomparison Project (CMIP5) and its response to global warming, respectively. The mechanisms leading to PDO modulations in a warming climate are discussed in section 5. The paper concludes with a summary and discussion of the results in section 6.

## 2. Data and models

To investigate the influence of global warming, CMIP5 model outputs are analyzed (Taylor et al. 2012). We use the historical and representative concentration pathway 4.5 (RCP4.5) scenarios to represent before and after global warming, respectively. Only models with outputs from 2061 to 2300 in RCP4.5 are chosen. Twelve models are selected based on the time span and the availability of all the variables we need (see Table 1). The temporal coverage in historical simulations is from 1871 to 2005.

To evaluate the model results, we use SST data from the Hadley Centre Global Sea Ice and Sea Surface Temperature dataset (HadISST; 1° × 1° resolution; Rayner et al. 2003), the oceanic current data from the Simple Ocean Data Assimilation reanalysis product (SODA 2.2.4; 0.5° × 0.5° resolution; Carton and Giese 2008), and the atmospheric geopotential height data from the National Oceanic and Atmospheric Administration (NOAA)/Cooperative Institute for Research in Environmental Sciences (CIRES) Twentieth Century Reanalysis, version 2 (20CRv2; 2° × 2° resolution; Compo et al. 2011). Reanalysis datasets including ERSST v4 (Huang et al. 2015) and Kaplan SST v2 (Kaplan et al. 1998) are also used for comparison. All the data are during the period of 1871–2005.

We mainly focus on the cold season and decadal variability, so all variables are 8-yr low-pass-filtered wintertime [December–February (DJF)] mean data. To remove the long-term trend, the linear trend and multimodel ensemble mean (MMEM; see Xie et al. 2013) are removed from the observational data and CMIP5 outputs, respectively.

The theoretical coupled model in this paper follows that of GM99 (referred to herein as the GM model) and is schematically illustrated in Fig. 1. The atmospheric component is a two-layer quasigeostrophic model governed by quasigeostrophic potential vorticity (QGPV) equations on a midlatitude beta plane and linearized around uniform zonal winds with vertical shear. Atmospheric fluctuations are assumed to be a steady-state response to the thermal forcing, which depends on the state of both the atmosphere and the ocean (Shutts 1987). The oceanic component is a 1.5-layer reduced gravity model linearized about an at-rest basic state, consisting of a wind-driven upper layer that satisfies the quasigeostrophic dynamics and a deep motionless layer. The atmosphere and the ocean are coupled via an oceanic mixed layer with a fixed depth. A linearized SST equation is adopted to describe the mixed layer temperature evolution, where the mean entrainment, geostrophic advection, and air–sea heat exchange are three key processes. Anomalous entrainment and Ekman advection do not significantly influence the model results (not shown) and are thus neglected for algebraic simplicity. The model only considers large atmospheric zonal scales and linear perturbations to a basic oceanic state, which has been proved to be a reasonable approximation for the first baroclinic mode (Sirven and Frankignoul 2000). More details can be found in GM99.

Previously, the GM99 model has been used to study the decadal climate variability in the midlatitudes over the North Atlantic and North Pacific, and has been proved to reasonably capture some basic characteristics of the Rossby wave dynamics and midlatitude growing coupled mode (GM99; Ferreira et al. 2001; van der Avoird et al. 2002; Arzel and Huck 2003; Fang and Yang 2011). Note that the GM model adopted here is not in conflict with the important pioneering work of Frankignoul et al. (1997) that studied the response of the ocean to stochastic atmospheric forcing. If the wind curl is assumed to be a stochastic process in the GM model, then it can be reduced to the model analyzed by Frankignoul et al. (1997), which has always been viewed as the null hypothesis of midlatitude decadal-scale SST variability.

## 3. Simulated PDO in CMIP5 models

CMIP5 models reasonably simulate the observed PDO decadal variability. The PDO emerges as the first empirical orthogonal function mode (EOF1) of DJF SSTA over the North Pacific in every CMIP5 model, and therefore we use the MMEM of the EOF1 to represent the model-simulated PDO (Fig. 2). The PDO in the historical run features a typical horseshoe-like SSTA pattern, broadly similar to the observed (Fig. 2a vs Fig. 3a). For the positive phase, the PDO display negative SST anomalies over the central Pacific and positive anomalies off the coast of North America. The pattern correlation coefficient between the historical MMEM and the observations reaches 0.85, which is significant at the 95% confidence level. Discrepancies, however, exist. The spatial pattern in the MMEM has its maximum SSTA over the Kuroshio–Oyashio Extension (KOE) region whereas the observed is located in the central North Pacific. This model bias is conceived to be related to ocean dynamics rather than atmospheric processes (e.g., Yim et al. 2013).

The PDO displays decadal periodicity. To calculate the spectrum of the PDO mode, we project the unfiltered data onto the EOF1 spatial pattern to obtain the “unfiltered” principal component (UPC) for each model by following Wu et al. (2003). The MMEM of the power spectrum of the UPC (calculated by averaging the spectra of individual models) displays distinct decadal fluctuations (Fig. 2b), with enhanced variability between 20 and 25 years per cycle. Please note that the “broad” peak of MMEM result is due to the relative large spread of spectral peaks among different models, even though all of them successfully capture the decadal variability of PDO. The observed PDO also displays a spectral peak similar to the model simulation, although barely passing the 90% confidence level (Fig. 3d). This weak decadal spectral peak can also be found in other reanalysis datasets, for instance ERSST v4 and Kaplan SST v2 (not shown). The lack of a distinct decadal spectral peak perhaps is mainly due to the short data period of the observations. Long records from paleo-reconstructions appear to support the existence of decadal periodic behaviors (e.g., Biondi et al. 2001). Moreover, model results that integrated hundreds of years also indicate significant spectral peaks at the decadal time scale, including the CCSM and GFDL models (e.g., Wu et al. 2003; Zhong and Liu 2009; Zhang and Delworth 2015).

The PDO decadal variability displays coupled interaction between the ocean and the atmosphere. The corresponding MMEM atmospheric regression pattern in CMIP5 displays a monopole over the basin, demonstrating the deepening of Aleutian low (Fig. 2a). This pattern resembles the leading EOF mode of wintertime atmospheric 500-hPa geopotential height (Z500) over the North Pacific (Fig. 2c), with the correlation coefficient reaching 0.91 (significant at the 95% confidence level). The resemblance between Z500 regression and EOF1 can also be detected in observations (Figs. 3a,e). Moreover, the spectrum of Z500 shows a distinct spectral peak at a decadal time scale in CMIP5 (Fig. 2d), observations (Fig. 3f), and AMIP experiments (not shown), indicating a feature of coupled ocean–atmosphere interaction at the decadal time scale.

In short, CMIP5 reasonably captures the main characteristics of the observed PDO decadal variability, including spatial pattern, time scale, and the nature of the air–sea coupled interaction.

## 4. PDO response to global warming

The PDO spatial pattern largely remains unchanged in global warming. In the RCP4.5 scenario, the PDO shows a horseshoe-like pattern with maximum anomalies located over the KOE region, highly resembling that in the historical run with a correlation coefficient of 0.9 (exceeding the 95% confidence level; Figs. 2a and 4a). The corresponding atmospheric regression displays a strengthened Aleutian low, similar to the first EOF mode of Z500 (Figs. 4a,c).

Unlike the spatial pattern, the PDO amplitude is significantly suppressed in global warming. The maximum anomalies of the EOF mode are significantly reduced by about 25% in RCP4.5 (cf. Figs. 2a and 4a). The differences in the MMEM EOF (RCP4.5 minus the historical run) display a negative PDO pattern, featuring positive anomalies over the central Pacific and negative anomalies off the coast of North America (Fig. 5a). Also, the 8-yr low-pass-filtered SST standard deviation (STD) features a reduced variability over the entire basin, with the maximum reduction over the KOE region (Fig. 5b).

In addition, the PDO decadal variability shifts toward higher frequency in global warming. The power spectrum peaks at around 15–20 years in RCP4.5, whereas it is about 20–25 years in the historical runs (Fig. 2b vs Fig. 4b). The shift of PDO frequency is robust and can be detected in each individual model (not shown). In response to such changes, the percentage of SST variance in global warming almost doubled in the 15–20-yr band but reduced by 67% in the 20–25-yr band. Moreover, the overall energy in decadal band decreases (from 68% to 59%), consistent with the weakened variability in spatial pattern (Fig. 5b).

The frequency shift is also detected in the atmosphere. In RCP4.5, the spectrum of Z500 UPC displays a peak at 15–20 years, resembling SST variability (Figs. 2b,d). This resemblance also reflects the nature of air–sea coupling at a decadal time scale. The frequency shift of the PDO in global warming demonstrated from the collection of CMIP5 models supports previous findings from individual models with idealized greenhouse forcing (Fang et al. 2014; Zhang and Delworth 2016).

## 5. Possible mechanisms

To understand the mechanisms leading to the suppression of PDO in global warming, it is useful to build a theoretical framework. Several idealized models have been developed to understand the PDO dynamics (e.g., Frankignoul et al. 1997; Jin 1997). These models are rooted in the assumption that SST decadal variability is driven from the oceanic integrating of the atmospheric stochastic forcing. Although helpful in advancing our understanding of the PDO dynamics, these models cannot explain the role of the coupled ocean–atmosphere interaction, as demonstrated above from CMIP5 models. To incorporate the role of air–sea coupling, we use the GM model. A detailed description of the GM model is given in the appendix.

### a. Decadal SST variability in GM model

To explore the features of this model, analytical plane-wave solutions are obtained by solving the governing coupled equations [(A12)–(A15)]. Anomalies in both the ocean and atmosphere (such as SST, streamfunctions) are assumed to be in wavelike form for both space and time. These waves have the same spatial scale and frequency with only amplitude difference and phase offsets. The solution is in complex form [(A22)], with the presence of imaginary terms indicating the possibility of growth or decay of the wave.

#### 1) The coupled growing mode

The frequency and growth rate of the growing mode, represented by the real and imaginary parts of the solution [(A22)] respectively, are plotted as a function of zonal wavelength in Figs. 6a and 6b. Note that these two variables are also a function of meridional wavenumber *l*. For simplicity we set *l* as *π*/(3200 km) hereinafter, commensurate with the meridional scale of the observed PDO (Mantua et al. 1997). The frequency of the coupled mode resembles that of the free oceanic baroclinic Rossby wave (the dashed line in Fig. 6a), except within a narrow range of wavelengths (2*π*/*k* = 11 800 km) where the corresponding growth rate reaches its maximum (Fig. 6b). We call the coupled mode with the maximum growth rate the fast-growing mode.

The fast-growing mode occurs at specific zonal length, the selection of which is determined by the state of the atmosphere. The growth is most rapid when the atmosphere reaches “equilibrium state” and its anomalies locate directly over the SST forcing. At equilibration, the atmospheric response is equivalent barotropic on the scale of a free stationary Rossby wave, with highs (lows) over warm (cold) water (i.e, and , thus 2*π*/*k* < 14 462 km; see appendix for more details about the variables; Shutts 1987). In such a state, the thermal equilibration time scale is much faster than the advective–propagation time scale (i.e, , , therefore 2*π*/*k* > 10 380 km), the diabatic heating rate is vanishingly small, and the air–sea temperature difference can be neglected (*m* = 0; see more details in GM99). Please note that both the above constraint conditions are sensitive to the mean climatology of the wind speed ( and ).

The fast-growing mode captures the main characteristics of the observed PDO (Fig. 6c). The atmosphere displays an equivalent barotropic structure, with high (low) pressure anomalies above warm (cold) SSTA at both atmospheric layers [Figs. 6c(1) and 6c(2)]. The oceanic current anomalies display in phase relationship with SSTA [Fig. 6c(3)]. The frequency of the fast-growing mode is *ω*_{r} = 0.38 yr^{−1} (the period is represented as 2*π*/*ω*_{r}), resembling the observations (Figs. 3d and 6a). Some details of the observed PDO, however, are not captured due to the simplified and idealized physics of the coupled model, such as the SSTA horseshoe-like patterns.

#### 2) SST amplitude

Now we proceed to derive the analytical expression of SSTA amplitude in the model and explore the controlling mechanisms.

Given the complex form of the dispersion relation (A22), we divide the solution of the fast-growing mode into a real part and an imaginary part to represent the frequency and the growth rate , respectively:

where SST′ denotes the “initial” amplitude of SSTA defined in (A19) and determines the growth (decay) of the coupled mode.

Different physical processes work together in the development of fast-growing mode. For the entrainment, a cyclonic (low pressure) anomaly causes Ekman upwelling in the ocean, which shoals the thermocline and cools the SST [Figs. 6c(1) and 6c(2)]. Through the response of coupling, the atmosphere adjusts to barotropic low pressure system, demonstrating a positive feedback (Shutts 1987; Peng and Whitaker 1999). The heat flux forced by the atmosphere, on the other hand, is shut off (i.e., *m* = 0) due to its damping effect (GM99; Deser et al. 1996; Seager et al. 2001; Kelly and Dong 2004). The shoaling of the thermocline further propagates westward via the oceanic Rossby wave. As to the advection, the oceanic streamfunction displays a southward (northward) flow to its west (east) in response to cold SST anomalies [Fig. 6c(3)], which brings cold (warm) water from high (low) latitudes and fuels (terminates) SST cooling to the west (east). Therefore, SST and oceanic streamfunction are in phase and propagate westward. The speed of advection propagation is equivalent to Rossby wave speed (see GM99 for details).

Without continent boundaries, the fast-growing mode would keep growing and its amplitude is proportional to the developing time. In the North Pacific, this developing time is set by the cross-basin time of the Rossby wave. Therefore, the SST amplitude in this study refers to that in one propagation period (i.e., ). Although the SST amplitude keeps growing in the GM model, the oscillating pattern can be observed at a fixed point due to the westward propagation of the coupled system via oceanic Rossby wave.

The “final” amplitude of SSTA developing in one period (unless stated otherwise, the amplitude hereinafter refers to this so-called final amplitude) is thus obtained as

It is evident that the growth rate and the frequency (growing time; ) are dominant factors, but with opposite effects (Fig. 7). Specifically, the amplitude strengthens as the growth rate increases, whereas it weakens with the increase of the frequency (growing time). This could be interpreted as follows. In a given period of time, the initial perturbation with a faster growth rate would evolve into larger amplitude. Given a fixed growth rate, the initial perturbation with a longer period to grow would develop into stronger anomalies. The depth-averaged entrainment velocity and meridional SST gradient *a*, on the other hand, play minor roles. Please note that these two parameters do not represent the entrainment and advection processes.

The ratio between the growth rate and frequency is another factor that significantly modulates the SST amplitude. Mathematically, this ratio, appearing in the exponent of (4), significantly modulates the SST amplitude. Physically, given the growing time and growing speed , the ratio measures the overall growth of the initial perturbation.

#### 3) Sensitivity to model parameters

In the theoretical model, the growth rate is highly related with oceanic stratification (Fig. 8a). Without heat flux forcing from the atmosphere, the growth rate of the fast-growing mode is determined by oceanic processes including entrainment and advection. The entrainment process is dependent on the depth-averaged entrainment velocity and the temperature difference between SST and subsurface temperature *θ*_{sub}. The subsurface temperature anomalies are related to oceanic stratification (see section e of the appendix), with larger anomalies corresponding to stronger oceanic stratification. Therefore, the strengthening of oceanic stratification contributes to the increase of the growth rate.

Other parameters exert subtle influences on growth rate, for instance and *a* (Fig. 8a). Note that and *a* do not represent the oceanic entrainment and advection processes, without involving the vertical temperature difference and oceanic currents, respectively. The effect of mixed layer depth (MLD; a parameter in the GM model) is encapsulated in depth-averaged entrainment velocity .

The frequency of the fast-growing mode is mainly dependent on oceanic stratification. As presented in Chelton et al. (1998), without considering background mean flow and other forcing, the freely propagating linear Rossby wave speed can be expressed as

where denotes the *n*th-order baroclinic Rossby radius of deformation (equivalent to defined in the GM model when *n* = 1), is the buoyancy frequency, and *H* represents the depth of the local ocean. According to this formula, the Rossby wave propagates faster in more stratified water, leading to an increase of frequency.

The ratio between the growth rate and frequency is associated with the depth-averaged entrainment velocity , the meridional SST gradient *a*, and the oceanic stratification (Fig. 8b). The ratio rises with the weakening of oceanic stratification, and with the enhancement of and *a*. Among these three parameters, the oceanic stratification plays the most important role.

### b. Modulations in global warming

The weakened amplitude and faster turnover of the PDO in a warming climate are also captured by the GM model, by using parameters estimated from CMIP5 (Table 2). In the GM model, the SST amplitude significantly decreases (Fig. 5c) and the frequency increases from 0.47 to 0.66 yr^{−1} in global warming (Fig. 9a). Note that the specific zonal wavelength (2*π*/*k*) at which the fast-growing mode develops remains the same (11 800 and 11 930 km in the historical run and RCP4.5, respectively), in accordance with the unchanged spatial pattern of the PDO. Such similar wavelengths are attributable to the nearly unchanged background and after global warming (not shown).

As discussed in section 5a, the amplitude of decadal SSTA is dominated by the growth rate, the frequency, and their ratio. Thus we explore the changes of these factors under global warming and probe into the possible causes.

The MMEM growth rate of the fast-growing mode is significantly enhanced in a warming climate by about 25% (Fig. 9b). Such enhancement can also be detected in every single model, with a minimum of 18% and a maximum of 29%. The increase of growth rate is mainly attributed to changes of oceanic stratification. In global warming, the oceanic stratification is significantly intensified due to faster warming in the upper than the deep ocean (Fig. 10a). This intensification is consistent among different CMIP5 models, with the MMEM increased by 20% (Table 2). The depth-averaged entrainment velocity and SST meridional gradient *a*, however, display a slight decreasing trend (Table 2), acting to reduce the growth rate. In total, the growth rate is significantly increased, acting to intensify SST amplitude. Such a growth rate effect is in contradiction to the suppression of the PDO, suggesting the opposing effects of other factors.

In global warming, the frequency of the PDO increases in response to the enhanced oceanic stratification. The Rossby wave is accelerated with enhanced stratification [see (5)]. To calculate the speed of Rossby wave, we use data for the first 30 years (1871–1900) in historical runs and for the last 30 years (2271–2300) in RCP4.5 simulations. The specific periods of time chosen in historical and RCP4.5 simulations will not greatly affect the final results, since we mainly focus on the climatological changes of these variables. The acceleration of the Rossby wave can be found in every CMIP5 model and at all latitudes, despite with different magnitudes. The MMEM increase of Rossby wave speed is around 30% in the subtropics (Fig. 10b) where the Rossby wave crosses the basin (Fang et al. 2014). This acceleration of the Rossby wave shortens the cycle of the PDO, consistent with the shifting of spectrum toward higher frequency (Fig. 2b vs Fig. 4b). By limiting the growth time of temperature perturbations, the increase of frequency weakens the PDO amplitude.

The ratio between growth rate and frequency decreases in global warming, acting to suppress PDO amplitude. This decrease is consistent among models, with a MMEM decrease by 2.32% (Table 3). Note that any subtle variations of the ratio may raise prominent changes of SST amplitude, since the ratio appears in the exponent in (4). To explore the contribution of each parameter, we change one parameter to its RCP4.5 value at a time while fixing the others to their historical values. As a result, all of the parameters contribute to the decrease of the ratio, with the oceanic stratification playing the dominant role.

The above suggests that the strengthening of oceanic stratification is responsible for the weakened PDO amplitude in global warming. To test this suggestion, we examine the relationship between changes of stratification and PDO amplitude (Fig. 11). In CMIP5, all of the models feature suppressed amplitude and enhanced oceanic stratification. The more the oceanic stratification strengthens, the more PDO amplitude weakens.

## 6. Summary and discussion

In this paper, changes of PDO under global warming are investigated using CMIP5 model outputs as well as a theoretical model. All CMIP5 models feature a suppressed PDO in a warming climate, with a shortened period and a weakened amplitude. The MMEM period of PDO moves from 20–25 years in historical runs to 15–20 years in the RCP4.5 simulations. Although the spatial pattern of the PDO remains nearly unchanged, its amplitude is significantly suppressed by 25%.

With the aid of the GM model, the possible mechanisms of the PDO modulation are explored. In this model, the SST amplitude is determined by growth rate and frequency. Specifically, the growth rate is determined by oceanic dynamical processes including entrainment and advection. The frequency depicting the growing time of the coupled mode is associated with the cross-basin time of the oceanic Rossby wave. The oceanic stratification intensifies in global warming, which amplifies the growth rate by enhancing the entrainment process and increases frequency by accelerating the oceanic Rossby wave. Although the effects of growth rate and frequency contradict each other, their combined effect acts to weaken the PDO amplitude.

Our study presents a simplified theoretical framework for understanding the PDO modulation in a warming climate. However, factors that are not included in this simple model may also contribute. For example, the meridional displacements of the subtropical and subpolar gyre boundaries are found to influence the PDO amplitude over the KOE region (Zhang and Delworth 2016). Recent studies further highlight the important role of midlatitude mesoscale air–sea interaction in regulating the western boundary current system (Ma et al. 2016). In addition, the effects of nonlinearity and finite width of the ocean have also been discussed in other works (e.g., Ferreira et al. 2001; van der Avoird et al. 2002). Moreover, the influence of remote tropical forcing change in a warming climate also needs further investigation in the future.

Because of the simplicity of the theoretical model, some characteristics of the fast-growing mode cannot always be found in the CGCM. For the fast-growing mode, the coupled system (including SST, thermocline, and atmospheric response) increases its amplitude and propagates westward with oceanic Rossby wave. Such characteristics can be found in the study of North Atlantic tripole decadal variability, by using the Fast Ocean Atmosphere model (Fig. 6 in Yang et al. 2012). In CMIP5, however, the propagation of oceanic Rossby wave cannot be detected over Kuroshio Extension where SST anomalies reach maximum values. This might be related to the strong eastward mean flow that prominently reduces Rossby wave speed. As a result, the development of the PDO in many CMIP5 models features growing amplitude but no prominent westward propagation, such as in HadGEM2-ES (not shown), GFDL, and CCSM3 (Zhong and Liu 2009; Zhang and Delworth 2016).

The suppression of the PDO in warming climate may reduce the perturbation of the internal variability to the global warming trend. Previous studies suggest that the internal variability of global-mean surface temperature is associated with the PDO (Zhang et al. 1997; Chen and Wallace 2015; Kosaka and Xie 2016), a mode sometimes also called the interdecadal Pacific oscillation (IPO; Power et al. 1999; Henley et al. 2017). Studies have attempted to explain the recent global warming hiatus based on the PDO/IPO influences (e.g., Kosaka and Xie 2013). It can be perceived that the suppression of PDO amplitude and increase of frequency in a warming climate may reduce the perturbation of the internal variability to the greenhouse gas–induced warming, and eventually lead to indiscernible hiatus events.

## Acknowledgments

We appreciate comments from anonymous reviewers, which improved the paper substantially. Discussions with Dr. Shantong Sun were helpful. This work is supported by the National Natural Science Foundation of China (NSFC) Projects (41490643, 41490640, 41521091, U1606402, 41606008), National Key R&D Program of China (2016YFA0601803), and the Fundamental Research Funds for the Central Universities.

### APPENDIX

#### Formulation of the GM Model

##### a. Atmospheric component

The two-layer atmosphere, bounded above by a lid and below by the ocean, is governed by its quasigeostrophic potential vorticity (QGPV) equation. Atmospheric perturbations are supposed to be in a steady-state response to the diabatic heating, which is prescribed to be at the interface between two layers and is parameterized in terms of the SST anomaly and the atmospheric temperature anomaly. The governing QGPV equations are linearized at the upper and lower levels around uniform zonal winds with vertical shear and can be expressed by the following barotropic and baroclinic components (after neglecting quadratic perturbation terms):

where *U* is the zonal wind speed, is the atmospheric streamfunction, is the meridional gradient of the Coriolis parameter *f*, and the symbol "SST" is the sea surface temperature anomaly (SSTA) to which the atmospheric temperature perturbation relaxes. Here, represents the square of atmospheric baroclinic Rossby radius of deformation, with *H* being the scale height of the troposphere and the atmospheric Brunt–Väisälä frequency. Barotropic and baroclinic parts of zonal wind speed and streamfunction are expressed as the sum and difference of *U* and at levels 1 and 2, using the notations and (see also GM99). Additionally, other parameters are defined by

where is the atmospheric potential temperature perturbation with a typical value (290 K) and is a constant scaling parameter that converts atmospheric temperature to baroclinic streamfunction through the thermal wind relation.

##### b. Oceanic component

###### 1) Dynamic equations

The model ocean consists of a deep motionless layer and a moving upper layer that obeys quasigeostrophic beta-plane dynamics. The motion for the active layer is characterized by the first baroclinic Rossby mode driven by the surface wind stress. In the longwave approximation, the governing QGPV equation linearized around a basic state of rest can be written as (after neglecting quadratic perturbation terms)

where is the depth of the upper ocean, is the reference density, denotes the surface wind stress, is the upper-layer geostrophic streamfunction, and is the square of oceanic baroclinic Rossby radius of deformation with being the density difference between the two layers (i.e., oceanic stratification).

###### 2) SST evolution

Following Frankignoul (1985), the evolution of the mixed layer temperature is given by

where SST is identical to the mixed layer temperature, is the Hamilton operator, denotes the temperature anomaly entrained from the subsurface ocean, is the specific heat of seawater (4200 J kg^{−1} K^{−1}), **V**_{mix} represents the horizontal velocity in the mixed layer, is the mean entrainment velocity, and is the Heaviside step function [ = 1 if *x* > 0; = 0 if *x* < 0]. Also, denotes surface heat flux and *h*_{mix} represents the mixed layer depth that is fixed in this model.

The temperature anomaly being entrained is assumed to be generated by the adiabatic undulation of isothermal surfaces underlying the mixed layer, parameterized by

where denotes the vertical variation of isothermal surfaces and can be converted to through thermal wind relation; is the oceanic Brunt–Väisälä; frequency is the thermal expansion coefficient of seawater (2.5 × 10^{−4} K^{−1}); And is a constant parameter that links oceanic temperature with upper-layer streamfunction.

Embedding (A6) into (A5), the linearized form of SST evolution equation can be rewritten (after dropping the quadratic perturbation terms) as

Here, the mean zonal SST gradient and anomalous entrainment are neglected. Note that and are positive parameters that measure the SST meridional gradient and strength of the entrainment process, respectively. They are defined as

The mean entrainment velocity can be estimated as the mean changing rate of mixed layer depth (MLD) between winter and summer. The MLD is defined as the depth where density differs from that in the surface by 0.03kg m^{−3}, as in Weller and Plueddemann 1996 and de Boyer Montégut et al. 2004. Thus the MLD can reflect the oceanic stratification. Note that different criteria of MLD definitions do not significantly change the final results.

##### c. Air–sea coupling

In this model, heating of the atmosphere is assumed to come from the underlying SST. An initial SST anomaly heats or cools the overlying atmosphere and thus changes the atmospheric circulation. The anomalous atmospheric circulation then leads to an anomalous surface wind field that drives oceanic motions in the upper layer, which in turn alters the initial SST anomaly.

The surface heat flux *Q* in (A5) is parameterized as a function of air–sea temperature difference (SST − *θ*_{a}), expressed by

where the air–sea heat exchange coefficient can be calculated as the linear regression coefficient between the turbulent heat flux and the surface air–sea temperature difference.

The surface wind stress is assumed to be proportional to the surface wind speed, expressed as

where **k** is the vertical unit vector and represents the surface atmospheric streamfunction derived from linear extrapolation of upper levels; is the dynamical coupling constant, as defined in GM99, which is expressed in terms of the mean surface wind speed , drag coefficient , and air density , as

##### d. Dispersion relation

By substituting (A9)–(A11) into (A4) and (A7), together with (A1) and (A2), the final closed governing equations can be obtained as shown [see also (22)–(25) in GM99]:

The solutions are assumed to have the following wavelike form:

where, , , and denote the amplitude of these variables, is the frequency, and *k* and *l* are *x*- and *y*-direction wavenumbers, respectively. Inserting (A16) into (A12)–(A15), the relations of the amplitude of these quantities can be obtained as follows:

with some key parameters defined as

Here denotes the oceanic baroclinic Rossby wave frequency under the longwave approximation; determines the relative strength of atmospheric barotropic and baroclinic modes; and and represent the thermal damping time scale and advective-propagation time scale of atmospheric waves, respectively. The ratio between and controls both the phase and amplitude of the atmospheric response to the SST anomaly forcing [(A18)]. We define and for notational simplicity and is a complex. The solutions of the above coupled equations can thus be obtained by combining (A17)–(A20) together and solving the quadratic equation of , given by [see also (41) in GM99]

The solutions represent the dispersion relation of the coupled modes modulated by air–sea interactions.

For the fast-growing mode, , *m* = 0, and *k* is fixed. The frequency of (A22) is reduced to

Here, the effect of oceanic stratification, which is represented by the Brunt-Väisälä frequency , is reflected in

##### e. Dependence of subsurface temperature anomaly on oceanic stratification

Since both oceanic and atmospheric perturbation fields move westward together (i.e., with the same spatial scale and in-phase relation; see Fig. 3) at the speed of oceanic baroclinic Rossby wave , we assume the following form of variables in (A4):

Here *B*(*t*) and *A*(*t*) are positive and represent the time-dependent amplitude of oceanic geostrophic streamfunction and wind stress curl, respectively. Both of them have the same spatial pattern and move westward together at the Rossby wave speed *c*. Inserting (A25) into (A4), we can obtain

If *A* is a constant, that is, the wind forcing is always the same, then the geostrophic streamfunction increases as the oceanic stratification strengthens. Considering that subsurface temperature anomaly *θ*_{sub} is proportional to geostrophic streamfunction [see (A6)], *θ*_{sub} will be proportional to oceanic stratification, with larger subsurface temperature anomaly in response to large stratification. With the coupled positive feedback, however, the amplitude of wind will also grow (i.e., *A* becomes larger). As long as *A* is positive, *θ*_{sub} will be positively correlated with stratification.

## REFERENCES

*Earth’s Climate: The Ocean–Atmosphere Interaction*,

*Geophys. Monogr.*, Vol. 147, Amer. Geophys. Union, 347–364, https://doi.org/10.1029/147GM19.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).