Abstract

A slowdown in the rate of surface warming in the early 2000s led to renewed interest in the redistribution of ocean heat content (OHC) and its relationship with internal climate variability. We use the Community Earth System Model version 1 to study the relationship between OHC and the interdecadal Pacific oscillation (IPO), a major mode of decadal sea surface temperature variability in the Pacific Ocean. By comparing the relative contributions of surface heat flux and ocean dynamics to changes in OHC for different phases of the IPO, we try to identify the underlying physical processes involved. Our results suggest that during IPO phase transitions, changes of 0–300-m OHC across the northern extratropical Pacific are positively contributed by both surface heat flux and oceanic heat transport. By contrast, oceanic heat transport appears to drive the OHC changes in equatorial Pacific whereas surface heat flux acts as a damping term. During a positive IPO phase, weakened wind-driven circulation acts to increase the OHC in the equatorial Pacific while the enhanced evaporation acts to damp OHC anomalies. In the Kuroshio–Oyashio Extension region, a dipole anomaly of zonal heat advection amplifies an OHC dipole anomaly that moves eastward, while strong turbulent heat fluxes act to dampen this OHC anomaly. In the northern subtropical Pacific, both the wind-driven evaporation change and the change of zonal heat advection along Kuroshio Extension contribute to the OHC change during phase transition. For the northern subpolar Pacific, both surface heat flux and enhanced meridional advection contribute to the positive OHC anomalies during the positive IPO phase.

1. Introduction

In the early 2000s the observed rate of increase in global mean surface temperature (GMST) slowed relative to the rapid warming period of 1970–98 (e.g., Medhaug et al. 2017; Yan et al. 2016; Xie and Kosaka 2017; Fyfe et al. 2016). This slowdown period is bracketed between two major El Niño events that occurred in 1997/98 and again in 2015/16. However, slight differences in the choice of start and end years lead to significant differences in the interpretation of the warming trend in this period. Some studies suggest that GMST trends during this period show no pause in warming and do not deviate significantly from the long-term warming trend (Medhaug et al. 2017; Lewandowsky et al. 2016), while other studies support a slowdown in the relative increase of GMST (Fyfe et al. 2016). After this period was initially referred to as the “surface warming hiatus,” the word slowdown was considered to be a better term to describe the period, which is now generally referred to as the “surface warming slowdown.”

In the early 2010s, phase 5 of the Coupled Model Intercomparison Project (CMIP5) showed that most coupled climate models, including model ensemble means, failed to capture the slowdown, showing instead that simulated warming in the early 2000s did not differ significantly from previous decades. Further analysis of the CMIP5 models by Meehl et al. (2014) suggests that only 10% of the 262 twentieth-century CMIP5 simulations were able to realistically reproduce the slowdown in the early twenty-first century. Meehl et al. (2014) attribute the skill of this subset of models to improvements in the simulated interdecadal Pacific oscillation (IPO) and conclude that to close the gap between observations and model simulations of trends in GMST, we need a better understanding of observational errors (Karl et al. 2015; Cowtan and Way 2014), external forcing (Santer et al. 2014), and internal variability (e.g., Meehl et al. 2013a; Dai et al. 2015). By using simulations with the correct phase of internal variability in addition with correcting external forcing and errors in observational data, Medhaug et al. (2017) significantly reduce the gap between simulated and observed GMST trends. In the following paragraphs, we narrow our focus on the discussion of internal variability.

Variability in the IPO over the eastern equatorial Pacific plays an important role in modulating GMST on decadal time scales (e.g., Dai et al. 2015; Kosaka and Xie 2016; Meehl et al. 2016; Dong and McPhaden 2017). The IPO is usually defined as the first mode of the empirical orthogonal function (EOF) of SST over the entire Pacific, after the removal of the externally forced trend and application of a low-pass (13-yr cutoff) filter on the SST data. Analysis shows that the IPO can have an influence on both global and regional climate, modulating precipitation over Australia and droughts over the southwestern United States (Power et al. 1999; Meehl and Hu 2006). Model-based studies (Meehl et al. 2011, 2013a) suggest that decades with a slowdown in warming are typically correlated with a pattern in Pacific SST trends similar to that of the negative phase of IPO (IPO−), with La Niña–like SST anomalies (such as cooling over the eastern equatorial Pacific). Based on observations, England et al. (2014) also show that the slowdown is related to the IPO− and strengthened trade winds. Kosaka and Xie (2013) show that climate models are indeed able to reproduce the slowdown in warming if forced by observed SST over the eastern equatorial Pacific.

The recent slowdown in surface warming was not only an atmospheric phenomenon, but also indicated a change in ocean circulation and a redistribution of ocean heat content (OHC). During the slowdown period, the radiation imbalance at the top of the atmosphere was nearly the same as previous decades (~1 W m−2; Hansen et al. 2011; Trenberth et al. 2014), indicating a constant energy flux into the Earth system. Due to the large heat capacity of the ocean relative to the atmosphere and the land, the majority of this excess heat was stored in the ocean (Yan et al. 2016). Multiple processes can modulate how this excess heat is redistributed within the ocean. Meehl et al. (2011) suggest that in the case of the recent warming slowdown, more excess heat has been stored in the deeper ocean than in the upper ocean. However, the relative roles of different ocean basins in the ocean heat uptake are under intense debate (Meehl et al. 2011; Chen and Tung 2014; Lee et al. 2015; Nieves et al. 2015; Liu et al. 2016), primarily due to large discrepancies among observations, reanalysis data, and model simulations. By comparing multiple observational datasets, Wang et al. (2018) suggest that all of the four major ocean basins play a role in the uptake of ocean heat during the recent warming slowdown.

Wang et al. (2018) show a high correlation (0.5–0.7) between the IPO index and the principal components of the first and second EOF modes of OHC in the upper 1500 m, indicating that the IPO plays an important role in global ocean heat content variability. By correlating changes in the distribution of global OHC in a preindustrial control simulation with changes in the phase of the IPO, Hu et al. (2018) show that global mean OHC redistributes vertically between the surface (defined as 0–100 m) and subsurface (defined as 100–300 m) layers during different phases of the IPO. During a positive phase of the IPO (IPO+), OHC increases at the surface and decreases in the subsurface, with a reciprocal response during a negative phase of the IPO (IPO−), although the OHC change at 0–100 m slightly dominates over the opposite change at 100–300 m. Hu et al. (2018) asserts that this redistribution of ocean heat is mainly contributed by the Pacific Ocean and is related to changes in the strength of subtropical cells in different phases of the IPO.

This paper attempts to further explore the mechanisms behind IPO-related changes in OHC in the Pacific basin. Specifically, we want to assess the relative roles of surface heat flux and ocean heat transport in the change of ocean heat. We find that in tropical region the changes in 0–300-m OHC are mainly driven by ocean circulation change, while surface heat flux just acts to damp the OHC changes. In the extratropics, both the surface heat flux and ocean dynamics positively contribute to the OHC changes at different locations. Important processes involved in OHC evolution are further identified in different regions. The rest of this paper is organized as follows: section 2 describes the model and data. Section 3 describes evolution of OHC changes and the heat budget during different IPO phases and transition periods. Section 4 describes the mechanisms behind the OHC and heat budget changes. In section 5 we summarize and discuss our findings.

2. Model and data description

The model used for our analysis is the Community Earth System Model version 1 (CESM1) developed at the National Center for Atmospheric Research (NCAR) in collaboration with scientists at laboratories of the U.S. Department of Energy (DOE) and universities (Hurrell et al. 2013). The CESM1 model uses the Community Atmosphere Model version 5 (CAM5; Neale et al. 2010), the Parallel Ocean Program version 2 (POP2; Danabasoglu et al. 2012), the Community Land Model version 4 (CLM4; Lawrence et al. 2011), and the Community Sea Ice Code version 4 (CICE4; Hunke and Lipscomb 2008). The twentieth-century climate simulated by CESM1 agrees reasonably well with observations (Meehl et al. 2013b).

The model data we use for this analysis are from the 2200-yr preindustrial (PI) control simulation of the CESM1 Large Ensemble (LENS) Project (Kay et al. 2015). The PI control is conducted with a horizontal resolution of 1° for all model components. We use the PI control for this analysis to avoid the potential interaction between internally generated variability and the climate system’s response to time-evolving external forcing. To avoid the model spinup period, we focus only on the last 1000 years of the simulation, and the model data are linearly detrended to remove any remaining model drift. All model variables analyzed in this paper are smoothed with a 13-yr low-pass filter to remove high-frequency fluctuation. As shown in Hu et al. (2018), this control simulation captures spatial and temporal characteristics of the observed IPO reasonably well in the Pacific, but the simulated SST spatial pattern does not agree well with observations in other basins.

To investigate the processes affecting regional OHC during IPO phase transition period, we build the IPO time evolution composites, such that the peaks with a normalized IPO index greater than 1 standard deviation were chosen (Fig. 1a). Then for each of these peaks, the time series extends forward and backward for 12 years to build a 25-yr cycle of the IPO. This 25-yr cycle is based on Hu et al. (2018), who indicate that the modeled IPO has a significant frequency of about 25 years. To minimize the possible overlaps between consecutive two periods, only 18 of the 30 potential IPO periods are used (gray lines in Fig. 1b and black bars at the upper edge of Fig. 1a indicating these 18 periods) to build the 25-yr-long composite IPO cycle (hereafter cpIPO; the letter “p” here indicates that this IPO cycle is constructed based on IPO positive peaks at year 0; black line in Fig. 1b). Using these same 18 IPO cycles, we reconstruct the OHC in the Pacific Ocean basin for each cycle, and then create a composite mean of the OHC to match the cpIPO composite. As shown in Fig. 1b, this 25-yr period roughly represents the time evolution of an IPO cycle.

To test the symmetry of the composite IPO cycle, a composite IPO cycle centered around negative IPO peaks (with IPO index less than negative one standard deviation) is also constructed and named the cnIPO composite (see Fig. S1 in the online supplemental material). In this composite, fifteen 25-yr-long IPO cycles are used.

3. How OHC and heat budgets vary

a. Time evolution of the composite OHC distribution

The time evolution of the composite OHC (Fig. 2 for the upper 0–300 m) in the cpIPO cycle shows clear phase transitions, first from an IPO negative phase (cpIPO−) to an IPO positive phase (cpIPO+), and then from a cpIPO+ to a cpIPO−. The magnitude of the OHC anomaly in the prior and subsequent cpIPO− phases (Figs. 2a,h) is weaker than in the IPO+ phase (Fig. 2e). The magnitude of the OHC anomaly in the later cpIPO− phase (Fig. 2h) is stronger than in the earlier cpIPO− phase (Fig. 2a). As shown in Fig. 1b, the normalized composite IPO time series shows a negative peak (with a value of −0.55) preceding the composite positive peak by 12 years. The next composite negative peak is 10 years later after the cpIPO+ peak (with a value of −1.05). This difference between these two IPO negative phases is due to the fact that only a limited number of IPO peaks were used and the cycle of each IPO period is not exactly 25 years. Hereafter, we refer to the middle year of the 25-yr time series composite as the cpIPO+peak, and we identify the prior and subsequent composite negative peaks as cpIPO−prior and cpIPO−post peaks, respectively.

For the upper 300-m ocean, the pattern of OHC anomaly at the cpIPO+peak (Fig. 2e) resembles the SST anomaly during the IPO+ phase: positive anomalies in the tropics, the subpolar region, and near the west coast of North and South Americas, and negative anomalies at the north and south subtropical regions. OHC anomalies near the cpIPO−prior and cpIPO−post negative peaks (Figs. 2a,h) are opposite to OHC anomalies at the cpIPO+peak. During the transition from the cpIPO+peak to the cpIPO−post, positive OHC anomalies near the Kuroshio Extension (Fig. 2e) grow in amplitude and extend downstream of the Kuroshio (Figs. 2f–h).

For simplicity, we did not separate the 0–300-m layer into surface 0–100-m and subsurface 100–300-m layers here, because the patterns of OHC evolution in these two layers are quite similar. One exception is at western equatorial Pacific where at the cpIPO+peak there are positive OHC anomalies near the surface but negative OHC anomalies below at the subsurface (see Figs. S2 and S3).

b. Heat budget during IPO phase transition

To study which processes contribute to the OHC change from one IPO phase to the opposite phase, we consider the composite evolution of heat budget terms around the cpIPO+peak described in section 3a. To build the phase transition climatology and to avoid any overlapping effect, only the data between cpIPO−prior and cpIPO+peak (11 years for each transition; i.e., from year −11 to year −1 prior to each of the 18 IPO+ peak years) are used to represent the transition from a negative to a positive phase (hereafter, Transition+), and the data between cpIPO−post and cpIPO+peak (nine years for each transition; i.e., from year 1 to year 9 after each of the 18 IPO+ peak years) is used to represent the transition from positive to negative phase (hereafter, Transition−). In this section, by ignoring the asymmetry between these two periods, their difference (Transition+ minus Transition−) is used to maximize the signal and to discuss the common features between them. In most regions, the pattern of heat budget terms in the two transition periods are just opposite to each other (not shown), but in certain regions, such as near the Kuroshio–Oyashio Extension, the asymmetry cannot be ignored and the role of different heat budget terms in transition periods must be considered with their time series (see section 4b).

According to the budget equation of OHC within volume V (OHC = cpρVT),

 
dOHCdtVcpρTtdV=VcpρdV[QsuT+D(T)],
(1)

where cp is ocean heat capacity, ρ is potential density, T is potential temperature, Qs is the heating rate due to surface heat flux, u is the Eulerian-mean velocity, and D(T) represents subgrid mixing. Subgrid mixing primarily includes vertical mixing (using K-profile parameterization, including the effects such as local static instability, shear instability, internal wave breaking, tides, etc.; see Large et al. 1994) and mesoscale eddy-induced mixing along isopycnal surface (using Gent–McWilliams parameterization; see Gent and McWilliams 1990). Full details of the subgrid mixing parameterization are documented in Smith et al. (2010). Since the relative change in density is much smaller than the relative change in temperature (Δρ/ρ ~ 0.1/1000, while ΔT/T ~ 0.5/300), we ignore the term related to ∂ρ/∂t, and briefly divide the total tendency of OHC into three terms:

 
dOHCdt=dOHCdt|shf+dOHCdt|adv+dOHCdt|mix,
(2)

where dOHC/dt|shf=VcpρdVQs represents the contribution from surface heat flux; dOHC/dt|adv=VcpρdV(uT)VcpρdV(uT)=VdV(HA) represents the contribution from convergence of heat advection (HA = cpρTu, representing heat advection); and dOHC/dt|mix=VcpρD(T)dV represents the contribution from subgrid mixing.

For a column of upper 300-m ocean with unit horizontal area, we can rewrite the terms in Eq. (2) as follows:

 
dOHCdt|shf=SHF=FSW+FLW+FLH+FSH,
(3)
 
dOHCdt|adv=(HA)dz,
(4)
 
dOHCdt|mix=cpρ(T)dz,
(5)

where SHF is total surface heat flux which is the summation of the shortwave (FSW) and longwave (FLW) radiation, and the sensible (FSH) and latent (FLH) heat fluxes. Equation (3) represents the ocean heat content changes due to changes of the total surface heat flux, which is the summation of the net shortwave and longwave radiation at the surface, and the sensible and latent heat flux due to the air–sea temperature difference and the saturation of the moisture in air. Equation (4) represents the ocean heat content changes induced by the ocean advection process (both horizontal and vertical), such as the divergence or convergence of heat into/out of a region/grid cell. Equation (5) is the heat content changes due to oceanic subgrid mixing processes.

Figure 3 show the spatial pattern of the terms in Eqs. (2)(5), as a difference between Transition+ and Transition− for the upper 300-m layer (for reference, the climatological mean of these terms is shown in Fig. S4). The sum of the contribution from surface heat flux, convergence of heat advection, and subgrid mixing is very close to the total OHC tendency derived directly from OHC time series (see the black contours in Fig. 3f whichs show the residual OHC tendency), although the residual is not tiny near the equator.

Now we examine the relative roles of surface heat flux (Figs. 3a–c), heat convergence due to advection (Fig. 3d), and subgrid mixing (Fig. 3e) in the total OHC tendency during IPO phase transitions (Fig. 3f). In total surface heat flux, the radiative flux (Fig. 3b) is dominated by shortwave flux, while longwave flux is usually opposite and weaker in amplitude than the shortwave flux. The surface turbulent flux (Fig. 3c) is dominated by latent heat flux in low latitudes, except in subpolar regions where the sensible heat flux change is comparable in amplitude to latent heat flux change. In the tropics, heat convergence in the surface layer due to advection (Fig. 3d) is the primary contributor to the positive total OHC tendency (Fig. 3f), especially over the eastern equatorial region. Although surface heat flux (Fig. 3a) contributes to the positive OHC tendency in the tropical central Pacific, surface heat fluxes tend to cool the surface layer through enhanced evaporation over the eastern equatorial Pacific (Figs. 3a,c). For the subtropical Pacific between 15° and 45°N, the negative total OHC tendency at central Pacific near 30°N resembles the changes in surface heat flux (Fig. 3a), primarily turbulent flux (Fig. 3c). The heat convergence due to advection contributes to the positive OHC tendency near the west coast of North America. Surface heat fluxes, primarily in the form of turbulent flux (Figs. 3a,c), contribute significantly to the total OHC tendency in the subpolar North Pacific.

Figure 4 shows the relative contribution of surface heat flux and ocean dynamics (as a sum of convergence of heat advection and subgrid mixing) to the absolute amplitude of total OHC tendency at different regions. The convergence of heat advection is the primary contributor to the total OHC tendency in the tropical Pacific, especially eastern equatorial Pacific, and along the west coast of North America. In the subtropics the surface heat flux and ocean dynamics both contribute to the total OHC tendency, while in the subpolar region the surface heat flux dominates over ocean dynamics. At the regions where surface heat flux dominates over ocean dynamics, turbulent heat dominates over radiative flux at most of those regions.

One caveat of the current composition method of representing phase transition is that the composite negative peaks (cpIPO−prior and cpIPO−post negative peaks) are more of a mathematical result rather than physically real negative peaks. Because the IPO time series is not purely sinusoidal, the real negative peaks prior/after the local positive peak in cpIPO composite members do not overlap with each other (Fig. 1b). As a result, the absolute value of composite mean IPO index at cpIPO−prior and cpIPO−post is much weaker than the composite mean IPO index at cpIPO+peak. To avoid this caveat, we modified our composite as below to test the robustness of the above results. We choose the real negative peaks prior/after the local positive peak, and we use the years between the positive peak and the prior (subsequent) negative peak to represent transition period from negative (positive) peak to positive (negative) peak, as shown in red (blue) periods in Fig. 1a. We first calculate the averaged heat budget terms within each negative-to-positive (positive-to-negative) transition period in each composite member, and then average them across the 18 composite members to represent our modified Transition+ (Transition−). Using modified Transition+ and Transition− (hereafter, mTransition+ and mTransition−), we get very similar results to those in Fig. 3 (see Fig. S5), suggesting these results are robust to the definition of transition period.

The results in this section and the previous section do not change too much when the composite evolution is constructed with a composite IPO negative peak at the center (see the online supplemental material). Figure S1 shows the constructed time series. Figure S6 shows the corresponding heat budget terms as mTransition+ minus mTransition− (the mTransition+ and mTransition− periods are marked as red and blue in Fig. S1). These heat budget terms are consistent (with opposite sign) with the results when the composite evolution is constructed with a composite IPO positive peak at the center (comparing Fig. 3 and Fig. S6).

c. Heat budgets during the IPO peaks

In addition to the phase transition period, it is important to also know whether the same physical processes are at work during the climatological IPO positive or negative phases, so that we can understand the evolution of heat budgets during the full IPO period. To build the mean state of the positive or negative IPO states, we use the average between 18 positive peaks used in constructing cpIPO (hereafter, cnIPO+peak) and the 15 negative peaks used in constructing cnIPO (hereafter, cnIPO−peak) to represent the heat budget at IPO peaks.

We plot the individual components of the heat budget as differences between the IPO positive peaks and negative peaks (cpIPO+peak minus cnIPO−peak) in Fig. 5. Total surface heat flux (Fig. 5a) has a negative anomaly near the equator and a strong dipole anomaly near the Kuroshio–Oyashio Extension, acting to damp the OHC anomaly in these regions (Fig. 2e). In these regions, the heat convergence due to advection is opposite in sign to the surface heat flux anomaly, acting to enhance the OHC anomaly (Fig. 5d). In the subpolar regions, both surface heat flux and heat convergence due to advection show generally positive anomalies, although near the Bering Sea the surface heat flux anomaly is negative due to the loss of sea ice. Subgrid mixing positive anomalies over the subtropical region (Fig. 5e) acting to damp the negative OHC anomalies in this region.

Among all the components of total surface heat flux, both radiative flux (Fig. 5b) and turbulent flux (Fig. 5c) contribute to the pattern of total surface heat flux change at low latitudes, while for high latitudes the change of total surface heat flux is dominated by turbulent flux. The change of turbulent flux below 50°N is mainly due to change in latent heat flux, but the positive anomaly of turbulent flux in north subpolar region is mainly contributed by changes in sensible heat flux (not shown).

Comparing the distribution of heat budget terms at IPO peaks (Fig. 5) to those during phase transition (Fig. 3), we find that the spatial patterns are very different. This difference is not surprising because the transition periods we defined and the IPO peaks do not overlap. However, we can see some similarity, such as for total surface heat flux at the north subpolar region, or for the convergence of heat advection near the west coast of North America. This similarity may suggest that for some regions the heat budget pattern in the transition period is dominated by the pattern at one or two years near the peak, which resembles the pattern at the peak. More details about the time evolution of heat budget terms are shown in section 4b.

4. Mechanisms behind the OHC changes and the heat budget

a. Decomposing individual heat budget terms

In this section, we explore the physical mechanisms behind the change of heat budget terms we see in the previous sections. We will connect the change of surface heat flux and heat advection to changes of cloudiness, surface air condition, winds, and ocean circulations. To make these connections, we primarily focus on the mechanism at the mean states of IPO peaks in this section in which the signals are much clearer. For the IPO phase transitions, lead–lag analysis is used in order to figure out the causes for the phase transition. However, this portion of the analysis is more challenging than the analysis of the IPO mean states due to the complicated interactions among many processes.

The radiative flux changes are mainly due to the changes of surface shortwave radiation reaching the surface. The shortwave radiation changes are mostly due to changes in cloud amounts, except at very high latitudes where changes in sea ice fraction impacts surface albedo, which further influences the amount of shortwave radiation absorbed by the ocean. Near the east side of the subtropical North Pacific, the positive shortwave flux is due to reduced low cloud amount which is related to a reduced lower-troposphere stability (figure not shown), consistent with Clement et al. (2009).

To understand the change in turbulent flux (which is mainly dominated by latent heat flux) at IPO positive/negative peaks, and also during phase transitions, we decompose the change in latent heat flux (LHF) into the separate contributions from changes in moisture and in wind speed: δLHF=δ(αU10Qdif)α(U10¯δQdif+Qdif¯δU10), where δ represents the deviation from the climatological mean, and U10 is the 10-m wind speed. The term Qdif represents the deviation from near surface saturation: Qdif = Qbottom − 0.98Qsaturation(TS), where Qbottom is the specific humidity at the lowest model layer, while 0.98Qsaturation is the saturation specific humidity inferred from surface temperature TS (0.98 is due to the effect of salinity), and α is a positive coefficient which relies on wind speed, although here we assume the change of α is small. LHF is positive for downward flux, and Qdif¯ is generally negative over open ocean. Figure 6 shows the change in U10¯δQdif, Qdif¯δU10, and their sum for IPO peaks and during phase transitions. The sum of these two terms (Figs. 6c,f) captures the pattern of turbulent flux change very well (comparing to Fig. 3c and Fig. 5c), with a pattern correlation coefficient (uncertered) of 0.92 to the latent heat flux change for IPO peaks and 0.96 for transition period. For IPO positive peak (negative peak), U10¯δQdif is negative (positive) where surface ocean becomes warmer (cooler) (Fig. 6a). Changes in relative humidity show a similar pattern (with small differences in the midlatitude central Pacific, figure not shown), which simply suggests that when SST increases, air in the near surface becomes less saturated. The term Qdif¯δU10 shows a quite similar (but opposite in sign because Qdif¯ is negative over most regions) pattern to δU10 (Fig. 7b). At the tropics, the negative anomaly of Qdif¯δU10 is consistent with the weakened trade winds (Fig. 7b). Changes in latent heat flux are mainly explained by the change of moisture in the extratropics (cf. Fig. 5c with Figs. 6a,b), while changes in wind speed also matter in the tropics. The relative contributions from changes in moisture and wind speed in the tropics and extratropics also applies to changes in latent heat flux during the IPO phase transitions (comparing Fig. 3c with Figs. 6c,d). The Qdif¯δU10 term and the near-surface wind change in the transition period show qualitative similarity to the pattern at IPO peaks in the Northern Hemisphere (comparing Figs. 7c to 7b and Figs. 6e to 6b), with an anomalous cyclonic circulation (and an anomalous low pressure) at the central Pacific south of the Aleutian low (Figs. 7a,c).

To understand the effect of oceanic advection on heat convergence, we further examine changes in ocean heat advection (Figs. 8 and 9). Because the upper branch of the subtropical cells (STCs) in the ocean is very shallow, we decide to separate the upper 300 m later to two layers, 0–100 and 100–300 m, in order to identify the contribution of the STCs to horizontal (zonal and meridional) heat advection. Changes in heat advection are mainly due to the change in velocity (not shown), which has a very similar pattern to the heat advection pattern in Figs. 8 and 9.

At IPO positive peak, zonal heat advection in the tropics is weakened in both surface and subsurface layers (Figs. 8a and 9a, compared to climatology, Figs. 8b and 9b). This is consistent with weakened trade winds (Fig. 7b). Meridional and vertical heat advection show weakened STCs. The climatological mean state of upper-100-m meridional heat advection exhibits poleward heat transport along each side of equator (Fig. 8d), while meridional heat advection in the 100–300-m layer indicates heat transport toward equator (Fig. 9d). Also at IPO positive peak, anomalies of meridional heat transport in the surface layer (Fig. 8c) in the central and eastern tropics indicate anomalous heat advection toward the equator, while at the subsurface the meridional heat advection in the central and western tropical Pacific indicates anomalous poleward heat transport (Fig. 9c). Upwelling along the equator (Fig. 8e) is also weakened as the result of the weakened trade winds. All of these results are consistent with weaker STCs at IPO positive peak than at IPO negative peak.

Zonal heat advection in both surface and subsurface layers at the IPO positive peak shows positive anomalies in the midlatitude Pacific between 30° and 40°N, and negative anomalies north of 40°N (Figs. 8a and 9a), contributing to the dipole anomalies of the convergence of heat advection near the Kuroshio–Oyashio Extension. Meridional heat advection in the surface and subsurface both show positive anomalies over large portions of the subpolar region north of 45°N and negative anomalies over the western and central Pacific between 20° and 40°N (Figs. 8c and 9c). These changes in meridional heat advection are coupled to cyclonic winds, which show westward wind anomalies north of 50°N and eastward winds south of 40°N (Fig. 7b). The cyclonic wind anomalies also lead to anomalous upwelling between 40° and 50°N (Figs. 8e and 9e).

b. Time series of different heat budget terms

In this section, we combine the results in sections 3b, 3c, and 4a and give a quantitative description of the evolution of different heat budget terms in several different regions identified by boxes in Fig. 2e. The reason for choosing these boxes is that the OHC anomalies in these boxes are the main contributors to the zonally integrated OHC pattern in Pacific basin at different latitudes. For each region, we use the composite method described in section 2 to get the time evolution of the heat budget terms in Eqs. (3)(5) for the 25-yr IPO cycle. Then these heat budget terms are integrated vertically and horizontally within each selected region to obtain the time series of these heat budget terms centered around the composite cpIPO+peak. For Eq. (4), we further decompose ∇ ⋅ (−HA) into the horizontal advection term [(−∂HAx/∂x) − (∂HAy/∂y)] and the vertical advection term (−∂HAz/∂z) to isolate the role of these two terms in the heat budget within each region. In each box, the climatological mean of different heat budget terms is shown in Fig. S7.

In the eastern equatorial Pacific (the purple rectangle in Fig. 2e; 150°–100°W, 15°S–15°N), at the cpIPO+peak, the total OHC tendency becomes zero and the accumulated OHC reaches its positive maximum (Fig. 10a). The anomalies of the oceanic heat convergence reach its maximum magnitude 2–3 years prior to the cpIPO+peak, contributing to the OHC increase. The total surface heat flux anomaly becomes the most negative around the cpIPO+peak, acting to damp the OHC positive anomaly. The damping effect of total surface heat flux is mainly due to the intensified evaporation (as a response of warmer SST and more negative Qdif; Fig. 6a) leading to an elevated latent heat flux loss (Fig. 11a). The changes of the oceanic advection terms are consistent with weakened STCs, which leads to a weakened poleward divergence flow and a reduced upwelling in this region. As a result, there are positive anomalies in horizontal heat advection and a negative anomaly in vertical heat advection (Fig. 10a).

In the subtropical North Pacific (the red rectangle in Fig. 2e; 170°E–160°W, 25°–35°N), the evolution of total OHC tendency in the surface layer is similar to the evolution of surface heat flux, with a negative anomaly during the transition from the cpIPO−prior to the cpIPO+peak phase and a positive anomaly during the transition from a cpIPO+peak to a cpIPO−post phase (Fig. 10b). Latent heat flux is the main contributor to this evolution of the surface heat flux (Fig. 11b). The evolution pattern of the latent heat flux results from the competing effect between near-surface saturation and wind (U10¯δQdif and Qdif¯δU10). By separately examining the time series of U10¯δQdif and Qdif¯δU10 averaged in this region (not shown), the negative anomaly in the latent heat flux from the cpIPO−prior to the cpIPO+peak phase results from Qdif¯δU10. The wind speed change is related to an anomalous northwesterly wind (Fig. 7b) that starts to grow around five years before the cpIPO+peak. On the contrary, U10¯δQdif acts to damp the surface SST anomaly, with positive value while the SST anomaly is negative (Fig. 6a). Besides surface heat flux, horizontal advection (mainly from zonal advection) also positively contributes to the total OHC tendency (Fig. 10b). The enhanced eastward heat transport at around 30°N near the Kuroshio (Figs. 8a and 9a) lasts for 5 years after the cpIPO+peak, which positively contributes to the eastward development of the positive OHC anomaly along the Kuroshio in Figs. 2e–g.

In the Kuroshio–Oyashio Extension region (the orange rectangle in Fig. 2e; 146°–170°E, 39°–44°N), the total OHC tendency is negative before the cpIPO+peak and becomes positive two years after the cpIPO+peak (Fig. 10c). Around the cpIPO+peak, the convergence of horizontal advection tends to amplify the negative OHC anomaly, consistent with the change of zonal heat advection (Fig. 10a). Comparing with Fig. 4, the positive contribution of surface heat flux to the total OHC change during phase transition (the difference between Transition+ and Transition−) shown in Fig. 4 primarily represents the Transition− period. The surface heat flux reaches a positive maximum two years after the cpIPO+peak and acts to damp the negative OHC anomaly. The change in surface heat flux is mainly contributed by latent and sensible heat fluxes (Fig. 11c).

In the subpolar North Pacific (blue rectangle in Fig. 2e; 160°E–130°W, north of 50°N), the total OHC tendency in the surface is positive before the cpIPO+peak, and becomes negative roughly 1 year later (Fig. 10d). Both the surface heat flux and the heat convergence due to horizontal advection contribute to the total OHC tendency. Among all the components of surface heat flux, shortwave flux and sensible heat flux positively contribute to the total surface heat flux in this region during the transition from the cpIPO−prior to the cpIPO+peak phase (Fig. 11d), while latent heat flux and longwave flux act to damp the positive OHC anomaly around cpIPO+peak. The positive horizontal heat advection convergence is related to anomalous poleward meridional heat advection around the cpIPO+peak (Figs. 8c and 9c).

5. Conclusions and discussion

We investigate the contribution of the IPO to the time evolution of OHC and the related heat budget in the Pacific Ocean using a CESM1 preindustrial control simulation. Based on the model-simulated IPO frequency, we construct a 25-yr nonoverlapping ensemble composite for both OHC distribution and heat budget analysis around the composite cpIPO+peak in order to examine the underlying physical processes affecting the OHC changes during IPO phase transitions and at IPO positive/negative peaks. Specifically, we assess the relative role of surface heat flux and oceanic heat transport in the temporal evolution of OHC distribution in the Pacific Ocean.

Around the cpIPO+peak year, the pattern of OHC anomalies in the upper 300-m layer shows a similar structure to SST anomalies in a typical IPO¯+ mean state, with positive OHC anomalies over the equatorial Pacific and subpolar North Pacific, and negative anomalies over the subtropical North/South Pacific. In the Kuroshio–Oyashio Extension region, a dipole pattern of OHC anomalies (negative in the north and positive in the south) may represent a shift in the position of the Kuroshio under different phases of the IPO. During the IPO phase transition from positive to negative (from cpIPO+peak to cpIPO−post), this dipole pattern enhances and move eastward. However, during the IPO phase transition from negative to positive (from cpIPO−prior to cpIPO+peak), the enhancement and eastward movement of this dipole pattern is not as clear.

We specifically analyze the relative contribution of surface heat flux and ocean dynamics to the change in OHC during IPO phase transitions by temporally averaging the heat budget terms over the IPO phase transition period. Our analysis shows that surface heat flux can positively contribute to the OHC changes over the west and central subtropical Pacific and subpolar region, suggesting that the OHC changes at the extratropics can be partly driven by surface heat flux change. Ocean dynamics also partly contributes to OHC change over the extratropics, especially off the west coast of North America, and dominates over the eastern equatorial Pacific region. Latent heat flux is the dominant contributor to total surface heat flux, indicating a coupling of air–sea interactions through SST, wind, and evaporation feedback processes (Xie 1999).

We also examined the temporal evolution of these heat budget terms in specific regions, attempting to identify underlying physical processes. Near the equator, weakened subtropical cells act to amplify the positive OHC anomaly in the central and eastern Pacific by reducing the divergence of heat from the equatorial tropics during the IPO+ phase, while enhanced evaporation tends to dampen the positive OHC anomaly. For the Kuroshio–Oyashio Extension region, the dipole change in zonal ocean heat advection acts to amplify the dipole OHC anomaly during the IPO phase transition from positive to negative while surface turbulent flux acts to damp this OHC anomaly. In the subtropical North Pacific, enhanced evaporation due to anomalous northwesterly wind during the transition from negative to positive IPO leads to elevated latent heat loss and a negative OHC anomaly during the IPO+ phase. Changes in zonal heat advection along the Kuroshio Extension also contribute to the OHC change during phase transition. In the subpolar North Pacific, both surface heat flux and anomalous meridional heat advection during the IPO+ phase contribute to the positive OHC anomaly in the surface layer.

One caveat of this work is that the analysis is based on one model simulation. The mechanisms mentioned here may vary for other models and may also differ from those in observations. Due to the lack of systematic observational data, it is nearly impossible to do a similar type of analysis using observed data. For the future work, we would like to investigate whether the same mechanism is at work for other coupled models, such as the models participating in phase 5 or phase 6 of the Coupled Model Intercomparison Project (CMIP5 or CMIP6).

Acknowledgments

We thank three anonymous reviewers for their most helpful suggestions. Z. Hu and Y. Hu are supported by the National Natural Science Foundation of China (NSFC), under Grants 41530423 and 41888101. Portions of this study were supported by the Regional and Global Model Analysis (RGMA) component of the Earth and Environmental System Modeling Program of the U.S. Department of Energy’s Office of Biological and Environmental Research (BER) Cooperative Agreement DE-FC02-97ER62402, and the National Science Foundation. The National Center for Atmospheric Research is sponsored by National Science Foundation of the United States of America.

REFERENCES

REFERENCES
Chen
,
X.
, and
K.-K.
Tung
,
2014
:
Varying planetary heat sink led to global-warming slowdown and acceleration
.
Science
,
345
,
897
903
, https://doi.org/10.1126/science.1254937.
Clement
,
A. C.
,
R.
Burgman
, and
J. R.
Norris
,
2009
:
Observational and model evidence for positive low-level cloud feedback
.
Science
,
325
,
460
464
, https://doi.org/10.1126/science.1171255.
Cowtan
,
K.
, and
R. G.
Way
,
2014
:
Coverage bias in the HadCRUT4 temperature series and its impact on recent temperature trends
.
Quart. J. Roy. Meteor. Soc.
,
140
,
1935
1944
, https://doi.org/10.1002/qj.2297.
Dai
,
A.
,
J. C.
Fyfe
,
S.-P.
Xie
, and
X.
Dai
,
2015
:
Decadal modulation of global surface temperature by internal climate variability
.
Nat. Climate Change
,
5
,
555
559
, https://doi.org/10.1038/nclimate2605.
Danabasoglu
,
G.
,
S. C.
Bates
,
B. P.
Briegleb
,
S. R.
Jayne
,
M.
Jochum
,
W. G.
Large
,
S.
Peacock
, and
S. G.
Yeager
,
2012
:
The CCSM4 ocean component
.
J. Climate
,
25
,
1361
1389
, https://doi.org/10.1175/JCLI-D-11-00091.1.
Dong
,
L.
, and
M. J.
McPhaden
,
2017
:
The role of external forcing and internal variability in regulating global mean surface temperatures on decadal timescales
.
Environ. Res. Lett.
,
12
,
034011
, https://doi.org/10.1088/1748-9326/aa5dd8.
England
,
M. H.
, and et al
,
2014
:
Recent intensification of wind-driven circulation in the Pacific and the ongoing warming hiatus
.
Nat. Climate Change
,
4
,
222
227
, https://doi.org/10.1038/nclimate2106.
Fyfe
,
J. C.
, and et al
,
2016
:
Making sense of the early-2000s warming slowdown
.
Nat. Climate Change
,
6
,
224
228
, https://doi.org/10.1038/nclimate2938.
Gent
,
P. R.
, and
J. C.
McWilliams
,
1990
:
Isopycnal mixing in ocean circulation models
.
J. Phys. Oceanogr.
,
20
,
150
155
, https://doi.org/10.1175/1520-0485(1990)020<0150:IMIOCM>2.0.CO;2.
Hansen
,
J.
,
M.
Sato
,
P.
Kharecha
, and
K.
Schuckmann
,
2011
:
Earth’s energy imbalance and implications
.
Atmos. Chem. Phys.
,
11
,
13 421
13 449
, https://doi.org/10.5194/acp-11-13421-2011.
Hu
,
Z.
,
A.
Hu
, and
Y.
Hu
,
2018
:
Contributions of interdecadal Pacific oscillation and Atlantic multidecadal oscillation to global ocean heat content distribution
.
J. Climate
,
31
,
1227
1244
, https://doi.org/10.1175/JCLI-D-17-0204.1.
Hunke
,
E.
, and
W.
Lipscomb
,
2008
: CICE: The Los Alamos sea ice model user’s manual version 4.1. Los Alamos National Laboratory Tech. Rep. LA-CC-06-012, 76 pp.
Hurrell
,
J. W.
, and et al
,
2013
:
The Community Earth System Model: A framework for collaborative research
.
Bull. Amer. Meteor. Soc.
,
94
,
1339
1360
, https://doi.org/10.1175/BAMS-D-12-00121.1.
Karl
,
T. R.
, and et al
,
2015
:
Possible artifacts of data biases in the recent global surface warming hiatus
.
Science
,
348
,
1469
1472
, https://doi.org/10.1126/SCIENCE.AAA5632.
Kay
,
J.
, and et al
,
2015
:
The Community Earth System Model (CESM) large ensemble project: A community resource for studying climate change in the presence of internal climate variability
.
Bull. Amer. Meteor. Soc.
,
96
,
1333
1349
, https://doi.org/10.1175/BAMS-D-13-00255.1.
Kosaka
,
Y.
, and
S.-P.
Xie
,
2013
:
Recent global-warming hiatus tied to equatorial Pacific surface cooling
.
Nature
,
501
,
403
407
, https://doi.org/10.1038/nature12534.
Kosaka
,
Y.
, and
S.-P.
Xie
,
2016
:
The tropical Pacific as a key pacemaker of the variable rates of global warming
.
Nat. Geosci.
,
9
,
669
673
, https://doi.org/10.1038/ngeo2770.
Large
,
W. G.
,
J. C.
McWilliams
, and
S. C.
Doney
,
1994
:
Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization
.
Rev. Geophys.
,
32
,
363
403
, https://doi.org/10.1029/94RG01872.
Lawrence
,
D. M.
, and et al
,
2011
:
Parameterization improvements and functional and structural advances in version 4 of the Community Land Model
.
J. Adv. Model. Earth Syst.
,
3
,
M03001
, https://doi.org/10.1029/2011MS00045.
Lee
,
S.-K.
,
W.
Park
,
M. O.
Baringer
,
A. L.
Gordon
,
B.
Huber
, and
Y.
Liu
,
2015
:
Pacific origin of the abrupt increase in Indian Ocean heat content during the warming hiatus
.
Nat. Geosci.
,
8
,
445
449
, https://doi.org/10.1038/ngeo2438.
Lewandowsky
,
S.
,
J. S.
Risbey
, and
N.
Oreskes
,
2016
:
The “pause” in global warming: Turning a routine fluctuation into a problem for science
.
Bull. Amer. Meteor. Soc.
,
97
,
723
733
, https://doi.org/10.1175/BAMS-D-14-00106.1.
Liu
,
W.
,
S.-P.
Xie
, and
J.
Lu
,
2016
:
Tracking ocean heat uptake during the surface warming hiatus
.
Nat. Commun.
,
7
,
10926
, https://doi.org/10.1038/ncomms10926.
Medhaug
,
I.
,
M. B.
Stolpe
,
E. M.
Fischer
, and
R.
Knutti
,
2017
:
Reconciling controversies about the ‘global warming hiatus.’
Nature
,
545
,
41
47
, https://doi.org/10.1038/nature22315.
Meehl
,
G. A.
, and
A.
Hu
,
2006
:
Megadroughts in the Indian monsoon region and southwest North America and a mechanism for associated multidecadal Pacific sea surface temperature anomalies
.
J. Climate
,
19
,
1605
1623
, https://doi.org/10.1175/JCLI3675.1.
Meehl
,
G. A.
,
J. M.
Arblaster
,
J. T.
Fasullo
,
A.
Hu
, and
K. E.
Trenberth
,
2011
:
Model-based evidence of deep-ocean heat uptake during surface-temperature hiatus periods
.
Nat. Climate Change
,
1
,
360
364
, https://doi.org/10.1038/nclimate1229.
Meehl
,
G. A.
,
A.
Hu
,
J. M.
Arblaster
,
J.
Fasullo
, and
K. E.
Trenberth
,
2013a
:
Externally forced and internally generated decadal climate variability associated with the interdecadal Pacific oscillation
.
J. Climate
,
26
,
7298
7310
, https://doi.org/10.1175/JCLI-D-12-00548.1.
Meehl
,
G. A.
, and et al
,
2013b
:
Climate change projections in CESM1(CAM5) compared to CCSM4
.
J. Climate
,
26
,
6287
6308
, https://doi.org/10.1175/JCLI-D-12-00572.1.
Meehl
,
G. A.
,
H.
Teng
, and
J. M.
Arblaster
,
2014
:
Climate model simulations of the observed early-2000s hiatus of global warming
.
Nat. Climate Change
,
4
,
898
902
, https://doi.org/10.1038/nclimate2357.
Meehl
,
G. A.
,
A.
Hu
,
B. D.
Santer
, and
S.-P.
Xie
,
2016
:
Contribution of the Interdecadal Pacific Oscillation to twentieth-century global surface temperature trends
.
Nat. Climate Change
,
6
,
1005
1008
, https://doi.org/10.1038/nclimate3107.
Neale
,
R. B.
, and et al
,
2010
: Description of the NCAR Community Atmosphere Model (CAM 5.0). NCAR Tech. Note NCAR/TN-486+STR, 268 pp., www.cesm.ucar.edu/models/cesm1.1/cam/docs/description/cam5_desc.pdf.
Nieves
,
V.
,
J. K.
Willis
, and
W. C.
Patzert
,
2015
:
Recent hiatus caused by decadal shift in Indo-Pacific heating
.
Science
,
349
,
532
535
, https://doi.org/10.1126/science.aaa4521.
Power
,
S.
,
T.
Casey
,
C.
Folland
,
A.
Colman
, and
V.
Mehta
,
1999
:
Inter-decadal modulation of the impact of ENSO on Australia
.
Climate Dyn.
,
15
,
319
324
, https://doi.org/10.1007/s003820050284.
Santer
,
B. D.
, and et al
,
2014
:
Volcanic contribution to decadal changes in tropospheric temperature
.
Nat. Geosci.
,
7
,
185
189
, https://doi.org/10.1038/ngeo2098.
Smith
,
R. D.
, and et al
,
2010
: The Parallel Ocean Program (POP) reference manual. Los Alamos National Laboratory Tech. Rep. LAUR-10-01853, 140 pp.
Trenberth
,
K. E.
,
J. T.
Fasullo
, and
M. A.
Balmaseda
,
2014
:
Earth’s energy imbalance
.
J. Climate
,
27
,
3129
3144
, https://doi.org/10.1175/JCLI-D-13-00294.1.
Wang
,
G.
,
L.
Cheng
,
J.
Abraham
, and
C.
Li
,
2018
:
Consensuses and discrepancies of basin-scale ocean heat content changes in different ocean analyses
.
Climate Dyn.
,
50
,
2471
2487
, https://doi.org/10.1007/s00382-017-3751-5.
Xie
,
S.-P.
,
1999
:
A dynamic ocean–atmosphere model of the tropical Atlantic decadal variability
.
J. Climate
,
12
,
64
70
, https://doi.org/10.1175/1520-0442-12.1.64.
Xie
,
S.-P.
, and
Y.
Kosaka
,
2017
:
What caused the global surface warming hiatus of 1998–2013?
Curr. Climate Change Rep.
,
3
,
128
140
, https://doi.org/10.1007/s40641-017-0063-0.
Yan
,
X.-H.
,
T.
Boyer
,
K.
Trenberth
,
T. R.
Karl
,
S.-P.
Xie
,
V.
Nieves
,
K.-K.
Tung
, and
D.
Roemmich
,
2016
:
The global warming hiatus: Slowdown or redistribution?
Earth’s Future
,
4
,
472
482
, https://doi.org/10.1002/2016EF000417.
For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).