Emergence of the North Pacific heat storage pattern delayed by decadal wind-driven redistribution

Domain: earth_science
Published: 2025-01-14
Editorial outcome: minor_revision

Abstract
Storage of anthropogenic heat in the oceans is spatially inhomogeneous, impacting regional climates and human societies. Climate models project enhanced heat storage in the mid-latitude North Pacific (MNP) and much weaker storage in the tropical Pacific. However, the observed heat storage during the past half-century shows a more complex pattern, with limited warming in the MNP and enhanced warming in the northwest tropical Pacific. Here, based on observational datasets, ocean model experiments, and climate models, we show that the emergence of human-induced heat storage is likely postponed in the North Pacific by natural variability until the late-21st century. Specifically, phase shifts of the Pacific Decadal Oscillation have vitally contributed to trends in the North Pacific winds during recent decades. Changes in surface winds drove meridional heat redistribution via Rossby wave dynamics, leading to regional warming and cooling structures and a more complex historical heat storage than models predict. Despite this, enhanced anthropogenic warming has already been emerging in marginal seas along the North Pacific basin rim, for which we shall prepare for the pressing consequences such as increasing marine heatwaves.

Introduction
The storage of excess heat caused by anthropogenic greenhouse warming in oceans1 is geographically inhomogeneous2–4. Climate models5,6 predict that the Southern Ocean3,7–9, the North Atlantic3,8, and the mid-latitude North Pacific (MNP) will warm up more quickly than other oceans (Fig. 1a). These ocean warming features accompany accelerated ocean currents8,10, poleward shift of storm tracks and westerlies11,12, increased extreme atmospheric rivers13, and altered marine biodiversity patterns14. While the enhanced heat storage in the Southern Ocean and North Atlantic has been witnessed in nature, that in the MNP has hardly emerged (Fig. 1b). Indeed, despite significant warming trends in marginal seas along the basin rim, there were weak warming or cooling trends in the MNP since the 1950s.

The tropical Pacific is among regions with the weakest heat storage in CMIP6 multi-model mean (MMM) (Fig. 1a). However, the Northwest Tropical Pacific (NWTP) has exhibited the strongest warming rate of the Pacific Ocean during the past decades (Fig. 1b). Increased ocean heat content (OHC) in the NWTP has led to rapid sea-level rise15,16, increasing marine heatwave and coral bleaching events17–19, and altered tropical cyclone behaviors20–22. It also exerted impacts on the downstream marginal seas through western boundary currents23–25 and the Indian Ocean through the Indonesian Throughflow26,27. The Pacific is projected to be the leading reservoir of anthropogenic heat among all oceans by the latter half of the 21st century3. Since prediction of climate change relies heavily on the ensemble of climate models, correctly interpreting the discrepancies between MMM and observation represents a vital scientific issue.

Although natural variability exists in individual models, the CMIP6 MMM mainly represents externally forced change28. Therefore, the contrasts between observations and CMIP6 MMM indicate either systematic model biases in simulating externally forced heat storage or substantial impacts from natural variability on observed heat storage. Although anthropogenic fingerprints have emerged in many aspects of the ocean3,29,30, natural variability remains influential in the observed OHC changes31,32. For instance, the persistent negative phase of the Pacific Decadal Oscillation (PDO) since 1998 has led to heat pile-up in the western tropical Pacific via the strengthened Pacific Walker Circulation33–35. Successful attribution of model-observation discrepancies in the North Pacific heat storage to natural variability may enhance our confidence in the projected future change by climate models.

Here, we set out to quantify the contributions of various processes to the observed North Pacific heat storage pattern, addressing in particular whether natural variability can largely explain the model-observation discrepancy. Our analysis is based on (1) five observational datasets36–40 to characterize the historical heat storage, (2) experiments of a forced ocean model to isolate the effects of surface heat fluxes and wind-driven ocean dynamics, and (3) 20 CMIP6 models to estimate anthropogenic fingerprints and when they may emerge (Methods). We demonstrate that surface wind changes arising from phase shifts of the PDO have driven basin-scale heat redistribution through Rossby waves and modulations in western boundary currents. This effect complicates the observed heat storage pattern by creating regional warming/cooling structures that conceal anthropogenic fingerprints. According to model projections, the human-induced heat storage pattern will likely hide until the late-21st century. These results provide useful implications for climate prediction in the North Pacific and marginal seas.

Results
Heat storage in the North Pacific
We begin with more details of the simulated and observed heat storage patterns in the North Pacific during the 1958–2021 period. In the MMM of 20 CMIP6 models, consisting of the historical and Shared Socioeconomic Pathway (SSP) 5–8.5 simulations before and after 2014, respectively, the entire MNP is characterized by enhanced 0–2000 m heat storage, in stark contrast to the warming minimum in the western tropical Pacific (Fig. 1a). The observation-based heat storage shows more regional structures (Fig. 1b and Supplementary Fig. 1); for example, there are alternating warming and cooling regions in the North Pacific, particularly in the western basin. The NWTP (e.g., 125°-180°E, 8°–18°N) stands out with a warming maximum, whereas a “warming hole” occurs in the MNP with insignificant trends and is surrounded by significant warming trends in marginal seas, such as the Bering Sea and the Gulf of Alaska. There is another cooling region near the western boundary of the 18°–30°N band, sandwiched by warming areas of the NWTP and the Kuroshio extension (KE) between 30°–40°N. Cooling trends are also seen beyond the North Pacific, such as the subpolar North Atlantic, southwestern subtropical Pacific, and southwestern subtropical Indian Ocean. CMIP6 MMM also shows weakened warming trends in these regions (Fig. 1a). Albeit with weaker intensity and smaller spatial range, the subpolar North Atlantic warming hole41 is clearly discernible in CMIP6 MMM. Therefore, models and observations are broadly consonant in all major ocean basins except for the North Pacific. The model-observation discrepancies in the North Pacific heat storage are worthy of in-depth investigation. These discrepancies are also seen in the 0–700 m OHC change (Supplementary Fig. 2), which is better covered by historical data in the 20th century than the deep ocean below 700 m.

We further examine the temporal evolution of OHC in key regions (Fig. 1c, d). CMIP6 models show an average heat storage rate (i.e., trend of regional average OHC anomalies) of 2.20±2.08×107 J m-2 yr-1 (equivalent to 0.71±0.67 W m-2; ± indicates the one standard deviation range of 20 models) in the MNP of 170°E-155°W, 40°-52°N, close in magnitude to that of the Southern Ocean (0.66±0.29 W m-2 in MMM for 55°-33°S). In comparison, a much weaker rate of 0.07±0.08 W m-2 is obtained from observation-based datasets for the MNP (Fig. 1c; ± indicates the one standard deviation range of 5 datasets). Meanwhile, the NWTP shows enhanced warming of 0.58±0.29 W m-2 in observations, one order stronger than the simulated rate of 0.06±0.33 W m-2 in CMIP6 models (Fig. 1d). Interestingly, CMIP6 agrees with observation in the OHC change of the entire North Pacific (120°E-80°W, 0°-60°N), which are 0.32±0.08 W m-2 and 0.32±0.18 W m-2, respectively, implying a heat redistribution over the North Pacific in the observed realization relative to the simulated pattern.

The observed OHC anomalies show prominent interannual and decadal variabilities in the MNP and NWTP (Fig. 1c, d). Meanwhile, we notice a large inter-model spread relative to the MMM change in the historical simulation of CMIP6 before 2014, which also indicates significant influence from natural variability. The effect of decadal natural variability is particularly notable from the 1990s through the mid-2010s. During this period, the observed OHC of the NWTP was substantially elevated and exceeded the one standard deviation range of CMIP6 models, while that of the MNP did not rise significantly as models expected. Given that anomalies during this period may represent an extreme case of natural climate variability, 20 models are likely insufficient to capture them. Also, we cannot ignore the possibility that CMIP6 models are incapable of reproducing such decadal variability. Due to the significant influence of natural variability, the contrasts between observations and CMIP6 MMM are sensitive to the time length but become stable when the time length exceeds 50 years (Supplementary Fig. 3). In the following, we examine how decadal natural variability affects the observed North Pacific heat storage.

Heat redistribution driven by wind changes
To understand how the heat storage pattern is formed, we performed experiments using the Hybrid Coordinate Ocean Model (HYCOM) with a coarse horizontal resolution of 0.5° (Methods) and prescribed reanalysis of atmospheric fields42 as the surface forcing. Despite regional simulation errors, the control run (CTRL) of HYCOM captured broad-scale features of the observed heat storage—the enhanced warming of the NWTP and the insignificant trends in the MNP (Fig. 2a). CTRL also well reproduced the prominent interannual and decadal variabilities (Fig. 2d, e). During the simulation period of 1958–2019, the correlation coefficients between CTRL and IAP are 0.44 and 0.81 for OHC anomalies in the MNP and NWTP (Fig. 2d, e) derived from CTRL, HTFL, and WND of HYCOM 0.5° simulation, referenced to those derived from IAP. Anomalies are obtained relative to a 1960–1980 baseline. HTFL+WND denotes the sum of the OHC anomalies of HTFL and WND.

With the aid of HYCOM experiments, we can separate the effects of heat redistribution induced by wind-driven ocean circulation changes and heat uptake through surface heat fluxes. The heat-flux run (HTFL) retains changes in surface heat fluxes and keeps wind stress and precipitation invariant (Methods), representing the effect of heat uptake. HTFL produces substantially stronger warming in the MNP than in the tropical Pacific (Fig. 2b) – a pattern resembling CMIP6 MMM (Fig. 1a). This highlights the dominance of heat uptake in shaping the model-projected heat storage. Alternatively, the wind stress run (WND), retaining changes only in surface wind stress and keeping other forcing fields unchanged (Methods), produces basin-wide cooling trends in the MNP and enhanced warming of the NWTP (Fig. 2c). These wind-driven changes greatly modify the pattern shaped by heat uptake and vitally contribute to the total storage in CTRL. Checking the temporal evolution clearly suggests that while heat fluxes drive quasi-monotonic warming trends in both regions (stronger in MNP), interannual and decadal fluctuations arise mainly from winds (stronger in NWTP) (Fig. 2d, e). Critically, winds have induced an OHC decrease in the MNP and an abrupt OHC increase in the NWTP since the late 1990s, which greatly altered the overall trends of 1958–2019. The linear separation of the wind and heat flux forcing effects inevitably lead to model errors in HTFL and WND. To assess the potential influence of these artificial errors, we compare the sum of HTFL and WND with CTRL (gray and red lines in Fig. 2d, e). The sum of HTFL and WND broadly agrees with CTRL in both OHC changes of the MNP and NWTP (Fig. 2d, e) and the heat storage pattern (figure not shown). Despite considerable model errors arising from the linear separation, HTFL and WND still provide useful insights into the underlying processes.

Then, we explore changes in surface winds. Reanalysis datasets42–46 suggest basin-scale trends of anti-cyclonic winds over the North Pacific since 1958 (Fig. 3a and Supplementary Fig. 4). Correspondingly, there are easterly winds and negative Ekman pumping velocity ω E (indicating downwelling; Methods) in the NWTP, which causes convergence of the upper-layer warm water47 and enhances the heat storage there. Meanwhile, we also see westerly winds and positive ω E (upwelling) in the eastern tropical Pacific, which might also affect the NWTP through Rossby waves. However, this effect fails to dampen the NWTP warming induced by local wind changes, probably owing to the dissipation of Rossby waves during their transition across the Pacific basin48. The anti-cyclone also involves strengthening westerlies north of 40°N, which slackens the heat storage in the subpolar North Pacific through Ekman upwelling and enhances heat pile-up in the 30°-40°N band through Ekman downwelling.

In CMIP6 MMM, trends of anti-cyclonic winds are confined north of 30°N, and westerly trends occupy the entire tropical Pacific (Fig. 3b), amplifying the model-projected heat storage pattern driven by heat uptake. Some key features of the observed wind trend are missed, including the easterlies in the NWTP and westerlies in the subpolar North Pacific. These differences in surface winds, critically accounting for the model-observation discrepancies in heat storage, likely arise from natural variability, given that the CMIP6 MMM is assumed to represent externally forced changes. The PDO, as the leading mode of decadal natural variability in the North Pacific, was in its positive phase during 1977-1997 and then shifted to its negative phase of 1998–2014 (Supplementary Fig. 5). A regression onto the negative PDO index (-PDO) (Fig. 3c) shows easterlies in the NWTP and westerlies in the northeastern subpolar Pacific – key features seen in observation but missed in CMIP6 MMM, suggesting that the PDO can largely explain the data-model difference in wind trends (Fig. 3a compared to Fig. 3b). The negative PDO also induces positive ω E anomalies in the western and central parts of the 18°-30°N band – a feature of obvious data-model discrepancy that explains the absence of the observed cooling trend near the western boundary in CMIP6 MMM (Fig. 1a, b).

The regression of OHC onto the negative PDO show some features resembling the observed heat storage, particularly the warming of the NWTP and the 30°-40°N band and the cooling of the MNP and the 18°-30°N band (Fig. 3d). The NWTP heat content shows correlation coefficients of -0.42 and -0.66 with the 8-year low-passed PDO index and the unfiltered December-January-February Oceanic Niño Index (ONI), respectively, both significant at the 95% confidence level (Supplementary Fig. 5a). This suggests a strong modulation of natural variability on the NWTP on interannual and decadal timescales, with heat pile-up in the NWTP under negative PDO and La Niña conditions.

In the mid-latitudes, the regression cannot explain the observed heat storage. For example, cooling trends occur mainly in the marginal seas in the regression (Fig. 3d) rather than in the MNP as in the trend pattern (Fig. 1b); the meridional dipole-like structure observed east of Japan, linked to the strengthened geostrophic transport of the KE jet (Fig. 1b), is replaced by prevailing warming in Fig. 3d. These discrepancies can be reconciled by considering delayed response of the ocean to wind forcing49. Lagged regressions show that the strengthened KE east of Japan is established ~4 years after the negative peak of PDO, and the MNP cooling takes ~9 years (Fig. 3e, f). The wind-driven OHC anomaly of the MNP, as measured by the WND run of HYCOM, also shows a lagged negative correlation with the PDO index (Supplementary Fig. 5b). The PDO shifted from the negative to the positive phase during the 1970s, and correspondingly there was a warming trend in the MNP throughout the 1980s and 1990s; subsequently, the opposite transition from the late 1990s through ~2010 led to a cooling trend of the MNP persisting through the late 2010s. This lagged response of the MNP region to PDO is also seen in observation (Supplementary Fig. 6). The lagged changes to the negative PDO phase are essential in shaping the OHC changes in the MNP and KE regions. Dynamical processes underlying this delayed response demand in-depth understanding, which is the focus of the following subsection.

Regional ocean dynamics
The above analysis points to the vital role of PDO phase shifts in shaping the historical heat storage pattern through wind-driven redistribution. However, correctly simulating the PDO-induced variability, which is critical for predicting regional OHC changes in the upcoming decades and their climatic and environmental impacts, remains a challenging task for state-of-the-art models50–52. Even with prescribed reanalysis winds, our 0.5° simulation of HYCOM fails to fully reproduce the observed heat storage in the MNP (Fig. 2a); for instance, the simulated cooling in the MNP is stronger in amplitude, broader in spatial extent, and shifted to lower latitudes in CTRL compared to that in observation. These discrepancies probably arise from complex regional ocean dynamics that are not properly represented by coarse-resolution models-such as the “bimodal” variability of the KE49,53.

In the mid-latitudes, a negative phase of the PDO drives downwelling Rossby waves through negative ω E anomalies between 150°-140°W (Fig. 3c). Along the latitudinal band of the KE (e.g., 32°-34°N), these Rossby waves propagate across the Pacific basin to the KE region east of Japan by ~4 years, as manifested in satellite-based sea surface height (SSH) anomalies (Fig. 4). These downwelling waves lead to a “stable” state of the KE system characterized by a strengthened KE jet and weakened mesoscale eddy variability49,53,54, as indicated by the enhanced anticyclonic “southern recirculation” locating south of the jet55,56 (Fig. 4a). The strengthened geostrophic KE jet manifests as a meridional dipole in OHC east of Japan (Fig. 3e). In addition to the PDO, the KE is also modulated by the upstream Kuroshio path south of Japan. The PDO shifted to a positive phase in ~2014, and correspondingly, an unstable dynamical state of the KE was established in early 2017, as indicated by the weakened southern recirculation (Fig. 4a). However, the Kuroshio large meander occurred in August 2017-with the Kuroshio following an offshore meandering path in the Shikoku basin south of Japan57–compelled the KE to switch back to the stable state in late 2017 and remains stable since then53,58. In contrast, earlier phase shifts of the PDO, such as those in 1992 and 2003, also led to an unstable state of the KE after ~3 years, which persisted for many years (Fig. 4a and Supplementary Fig. 7). Albeit interrupted in 2017, there has been a persistently stable KE since ~2009. This led to an overall trend of the KE toward its stable state during the past decades, given the more unstable KE during the 1980s and 1990s due to the positive PDO phase (Supplementary Fig. 7). This trend manifests as a dipole-like structure in heat storage east of Japan (Fig. 1b). By affecting the overlying storm tracks and large-scale winds, the stable KE also contributed to the cooling in the MNP and the warming in the eastern basin53,59.

The response of KE to Rossby waves is successfully reproduced by a 0.1° simulation of HYCOM (Methods; Supplementary Fig. 7) that can better resolve the Kuroshio path and mesoscale eddies. As a result, this 0.1° simulation can well represent the heat redistribution driven by PDO winds, including the meridional dipole east of Japan with a lag time of ~4 years (Fig. 5a) and the cooling in the MNP with a lag time of ~9 years (Fig. 5b) to the negative peak of the PDO. This ~9-year lag time reflects the slow propagation of upwelling Rossby waves between 40°–52°N (Supplementary Figs. 8 and 9) generated by positive ω E anomalies near the west coasts of North America (Fig. 3c). In addition, the delayed impact of the KE system also contributes to the ~9 year lag time (Supplementary Fig. 9). There are also detailed discrepancies between the 0.1° simulation and observation in the KE region (Fig. 5a, b compared to Fig. 3e, f), which may result from the lack of feedback between mesoscale eddies and the atmosphere in our standard-alone HYCOM simulation. The eddy-atmosphere feedback has been demonstrated to efficiently damp mesoscale eddies and amend the simulated strength and pathway of the KE60.

Most CMIP6 models are of horizontal resolutions of ~50 or ~100 km in their ocean components, close to our 0.5° HYCOM simulation. We examine how the PDO-induced heat redistribution is represented in these models. Regressions using individual models (Methods) suggest that CMIP6 can capture some large-scale features of the PDO-induced changes in surface winds and OHC (Fig. 5c, d and Supplementary Fig. 10). The cooling of the MNP is seen at a lag of ~9 years, but the spatial structure does not correlate with the observed regression (Fig. 5d compared to Fig. 3f). Owing to coarse resolutions, CMIP6 models cannot realistically represent the complex dynamics governing the lagged response to the PDO. The PDO’s signatures in the NWTP are weaker in models than in observations, linked to uncertain tropical Pacific wind changes (Supplementary Fig. 10). This indicates a weaker coupling between the PDO and the tropical Pacific climate in models61. By prescribing the observed (reanalysis) surface winds, the 0.5° HYCOM simulations can well reproduce the response of NWTP changes associated with the PDO winds (Supplementary Fig. 11). Furthermore, the improved fidelity of the 0.1° HYCOM (Fig. 5a, b) relative to the 0.5° HYCOM and CMIP6 models highlights the necessity of sufficient model resolutions to account for the complex ocean dynamics involved in heat redistribution.

To further investigate the PDO’s impact, we select 5 model realizations with the simulated PDO showing the highest positive correlations with the observed PDO (+5 models) and 5 model realizations showing the highest negative correlation with the observation (-5 models). As such, contrasting the +5 and -5 models mimics the impacts of the observed PDO phase shifts on the OHC trends since 1958 (Supplementary Fig. 12). The composite of +5 minus -5 models reassures that the PDO phase shifts in observation can dampen the heat storage in the MNP and enhance the warming of the NWTP (Supplementary Fig. 12c) through wind-driven heat redistribution (Supplementary Fig. 12d). Therefore, despite unresolved regional structures due to coarse model resolutions, discrepancies between CMIP6 models and observations in the North Pacific heat storage pattern, particularly the contrasts in the NWTP and the MNP, are primarily explained by natural variability represented by the PDO.

Emergence of anthropogenic heat storage
The above analysis indicates that the emergence of anthropogenic heat storage has been delayed by natural variability in some key regions of the North Pacific and will eventually take place in the future as the emission of greenhouse gases continues. Therefore, it is instructive to estimate the time of emergence (ToE) of anthropogenic heat storage based on CMIP6 simulations (Fig. 6), although we fully acknowledge that such estimation is subject to potential influence from inaccurate ocean dynamics as discussed earlier. Here, we primarily utilize the piControl simulation, historical simulation, and projections under three emission scenarios (SSP1-2.6, SSP2-4.5, and SSP5-8.5) of 20 CMIP6 models to determine the “signal-to-noise” ratio and compute the ToE for the 0–2000 m OHC change (Methods29,62).

Under the “middle-of-the-road” scenario of SSP2-4.5, a scenario viewed as where social and economic trends do not shift markedly from historical patterns, the anthropogenic heat storage (Fig. 6a) has already emerged (e.g., ToE earlier than the 2020s) in the Southern Oceans between 55°–33°S and the tropical and northern Atlantic. In the North Pacific, the ToE is as early as the 2010s in marginal seas along the basin rim, including the western coasts of North America, the Bering Sea, the Okhotsk Sea (most of which is shallower than 2000 m and not estimated for OHC and ToE), and the Japan Sea (Fig. 6a). By contrast, anthropogenic signals have not emerged by present (the 2020s) in the basin interior and even until the 2050s in the MNP, KE, and NWTP where we see obvious model-observation discrepancies. The ToEs under the SSP1-2.6 and SSP5-8.5 scenarios are similar in spatial distribution to those under SSP2-4.5 but postponed and advanced in average time (Supplementary Fig. 13a, b), respectively. In addition, the signal-to-noise ratios calculated from the observational datasets are generally in line with the ToEs based on CMIP6 simulations in spatial pattern (Supplementary Fig. 13), suggesting the reasonableness of the ToE estimation.

In marginal seas of the North Pacific, such as the Bering Sea and the eastern boundary coasts, anthropogenic heat storage has already been emerging. The multi-model median ToE is during the 2010s and 2020s in these regions regardless of the emission scenario (Fig. 6b, c). By contrast, the MNP, KE, and NWTP regions show delayed ToEs and relatively high sensitivities to the emission scenario (Fig. 6d–f). Under the SSP1-2.6 and SSP2-4.5 scenarios, anthropogenic heat storage is projected to emerge in the latter half of this century in these regions, such as the 2080s in the KE region. Under the SSP5-8.5 scenario, the ToE is as early as the 2040s in the MNP, the 2050s in the KE, and the 2060s in the NWTP. To summarize, despite uncertainties arising from inter-model spread, the North Pacific heat storage is projected to remain under a significant impact of natural variability until the late 21st century, except for marginal seas along the basin rim. Therefore, we expect that the heat storage pattern of the North Pacific will alter but probably still differ from that projected by climate models in the upcoming several decades.

Summary and Implications
Based on observational datasets, ocean model experiments, and CMIP6 model simulations, we demonstrate a strong modification effect by the PDO on the observed heat storage pattern in the North Pacific since the mid-20th century. Specifically, phase shifts of the PDO in recent decades have altered the trends in surface winds over the North Pacific. Surface wind changes induced meridional heat redistribution through Rossby waves and variability of the KE system, effectively erasing the warming in the MNP and fueling the heat pile-up in the NWTP. These effects led to a more complex heat storage pattern in observations than in models by creating regional warming/cooling structures that mask the human-induced fingerprints.

Recognition of the strong influence of the PDO on heat storage provides critical insights into near-term climate prediction when anthropogenic signals have not emerged. The persistent negative phase of PDO since 1998 terminated in the mid-2010s (Fig. 4c). If the PDO switches to a positive phase in the upcoming decade, the MNP “warming hole” shall vanish soon, and the basin-scale heat storage will evolve toward the pattern projected by models. However, the recent triple La Niña events in 2020–202363,64 indicate a possibility of a lengthened negative phase, which acts to maintain the existing heat storage pattern. The broad impacts of the PDO on surface atmospheric and ocean conditions have been well established, such as storm tracks, atmospheric rivers, and Pacific salmon production65–67. Here, we add that the PDO is also critical for the transient response of the North Pacific to anthropogenic forcing. Through modulating wind-driven ocean circulations, the PDO fundamentally regulates the distributions of anthropogenic heat, CO2, and the Fukushima nuclear effluents in the North Pacific and their spread into the Indian and Arctic Oceans via the Indonesian through-flow26,27 and Bering Strait68, respectively. These implications highlight the urgent need for an accurate initialized prediction of the PDO and its far-reaching impacts.

Climate models tend to underestimate the PDO variability and its impacts over the historical period51. As a result, none of the 20 CMIP6 models has produced a historical heat storage pattern resembling the observation (not shown). Increasing model resolution seems promising. The representation of natural variability in the North Pacific is considerably improved in ~10km simulations compared to ~100km simulations69,70; our analysis based on the 0.1° HYCOM simulation also supports this notion. In addition, increasing ensemble members can better account for the extreme cases of natural variability, such as the 1998–2014 negative phase71, and improve the ToE estimate.

The human-induced heat storage pattern in the Pacific Ocean is projected to eventually emerge by the late 21st century, which is much later than in the Atlantic and Southern Oceans. Apart from the strong impacts of natural variability, the postponed ToEs in the Pacific are linked to the lack of a deep-reaching meridional overturning circulation, which also holds true for the Indian Ocean. In the North Atlantic and Southern Oceans, overturning circulations can efficiently communicate anthropogenic heat to the deep ocean, enhancing the basin-average heat uptake72 and the signal-to-noise ratio in heat storage. Nevertheless, human-induced warming has already been emerging in marginal seas along the North Pacific basin rim, such as the Bering Sea, Okhotsk Sea, Japan Sea, Gulf of Alaska, and western coasts of North America. Increasing frequency and severity of marine heatwaves have been extensively reported recently73–75, with adverse influences on marine ecosystems and regional socioeconomics76,77. Anthropogenic fingerprints in these marginal seas could be enhanced by influence from lands where the air temperature warms quicker than in oceans78 or coastal processes that are insufficiently resolved by climate models. This calls for extensive investigations to reveal the mechanisms underlying the warming of marginal seas and multi-disciplinary implications, which shall aid the prediction and decision-making.

Methods
Datasets
To estimate the 0–2000 m ocean heat content (OHC) trend since the mid-20th century, ocean temperature data since 1958 from four observational analyses and one reanalysis product are utilized. The four observational analyses are the Institute of Atmospheric Physics ocean analysis (IAP) for 1958–2021 provided by the Chinese Academy of Sciences36, the World Ocean Atlas (WOA) for 1958–2018 provided by the National Oceanic and Atmospheric Administration’s (NOAA’s) National Centers for Environmental Information of the United States38; the Ishii analysis for 1958–202137; and the version 4.2.2 of the Met Office Hadley Centre “EN” series of datasets (EN4) for 1958–2021 provided by the Met Office of the United Kingdom39. The one ocean reanalysis product is the Ocean Reanalysis System 4 (ORA-S4) for 1958–201740 from the European Centre for Medium-Range Weather Forecasts (ECMWF). All five datasets provide 1°×1° ocean temperature fields. WOA provides pentad-mean data in annual intervals, and monthly fields of IAP, Ishii, EN4, and ORA-S4 are averaged into annual-mean data. Surface winds of four atmospheric reanalysis products are used, including the 1.25°×1.25° Japanese 55year Reanalysis (JRA55) product for 1958–202142, the 1°×1° ORA-S4 product as a combination of the 40year ECMWF Re-Analysis (ERA-40) of 1958-198843 and the ECMWF Interim Re-Analysis (ERA-Interim) of 1989-201744, the 2.5°x2.5° National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR) reanalysis for 1958–202145, and the 0.25°x0.25° fifth generation ECMWF reanalysis (ERA5) for 1958–202146. To examine the changes in the Kuroshio Extension (KE) jet, the 0.25°x0.25° sea surface height data distributed by the Copernicus Marine and Environment Monitoring Service (CMEMS) for 1993–2021 is used.

The coupled model simulations
Monthly outputs of the pre-industrial control (piControl) (the last 100 years), historical (from 1958–2014), Shared Socioeconomic Pathways 126 (SSP1-2.6), SSP2-4.5, and SSP5-8.5 (from 2015–2100) simulations of twenty models of the Coupled Model Intercomparison Project Phase 6 (CMIP6) are analyzed. These models include the following: ACCESS-CM2, ACCESS-ESM1-5, BCC-CSM2-MR, CAMS-CSM1-0, CESM2-WACCM, CNRM-CM6-1, CanESM5, EC-Earth3, EC-Earth3-Veg, FGOALS-g3, FIO-ESM-2-0, GFDL-ESM4, GISS-E2-1-G, INM-CM4-8, INM-CM5-0, MPI-ESM1-2-LR, MRI-ESM2-0, NESM3, NorESM2-LM, NorESM2-MM. All CMIP6 models integrate for at least 500 years under the pre-industrial CO2 level (piControl) and then from 1850–2014 to form the historical simulation. In the historical simulation, all models are forced by the observed external forcing (such as CO2, aerosols, and ozone)5. Owing to the nonlinearity of the climate system, the phase of internal variability (such as El Niño-Southern Oscillation and PDO) is different across models and different from observation. Therefore, anomalies arising from internal variability cancel each other when we average these 20 models to obtain the multi-model mean (MMM). As such, CMIP6 MMM primarily represents the externally (mainly anthropogenic) forced change, and the inter-model spread represents internal variability and differences in model physics. Selection of these 20 models is dictated by data availability for the three emission scenarios when our work begins. Monthly fields of CMIP6 simulations are interpolated onto the same set of 1°×1° grids and averaged into annual-mean data to match observational datasets.

HYCOM experiments
We use the Hybrid Coordinate Ocean Model (HYCOM) version 2.3.0179 to simulate and understand the ocean heat storage pattern. The HYCOM 0.5° simulation adopts a quasi-global ocean domain (78°S-75°N, 30°-389.5°E) with resolutions of 0.5°×0.5°×50 levels. The model is spun up from a state of rest for 600 years under the daily JRA55 atmospheric forcing of 1958, and the atmospheric surface forcing fields include the wind stress, net shortwave and longwave radiations (SWR and LWR), wind speed, precipitation, and air temperature and humidity. The control run (CTRL) is forced with realistic daily atmospheric forcings of JRA55 reanalysis from 1958 to 2019. CTRL contains complete physical processes of forming ocean heat storage pattern and is compared with observational data to evaluate the model performance. Two sensitivity runs are conducted to explore the effects of different drivers on the ocean heat storage pattern. Only wind stress adopts realistic daily fields in wind stress run (WND), while all the other forcings, such as wind speed and radiations, are fixed to the 1958 fields as in spin-up. In the WND, OHC trends primarily arise from heat redistribution due to wind stress-driven changes in ocean circulation processes. The heat flux run (HTFL) only adopts realistic SWR and LWR, 2m air temperature, specific humidity, and wind speed for 1958–2019 but fixes wind stress and precipitation to the 1958 fields. As such, changes in these fields of HTFL affect OHC trends mainly by altering the heat fluxes (SWR, LWR, sensible, and latent heat fluxes), representing the role of heat-driven processes or ocean heat uptake. This linear separation of the wind and heat flux forcing effects inevitably leads to nonlinear interactions or model errors, but HTFL and WND still provide useful insights into the underlying physical processes. To examine the impacts of decadally varying KE jet and its associated mesoscale eddies, we also conducted HYCOM 0.1° simulation using the same domain as the HYCOM 0.5° simulation but with resolutions of 0.1°×0.1°×50 levels. After a 20-year spin-up, the HYCOM 0.1° simulation is integrated forwards from 1979 to 2021 under hourly ERA5 atmospheric forcing. The comparison of the HYCOM 0.1° simulation and the HYCOM 0.5° CTRL allows us to test the possible impact of mesoscale oceanic processes on heat storage pattern in the North Pacific.

Definitions
The OHC is calculated by integrating ocean temperature T within the upper 2000 m as

OHC = \int_{0}^{2000} \rho c_p T dz (1)

where c_p = 4096 J (kg°C)-1 is thermal capacity, and \rho = 1025 kg m-3 is seawater density. The Ekman pumping velocity is calculated as \omega_E = curl(\tau/f) \rho-1, where \tau is the wind stress vector, and f is the Coriolis parameter. The least square method is used to estimate the linear trend. Unless otherwise stated, the statistical significance of the trend is defined as >95% confidence level based on a modified Mann-Kendall test. The estimation of the confidence level of the correlation is based on a two-tailed Student’s t-test using the effective degrees of freedom, estimated as N* = N(1-r1r2)/(1+r1r2), where N is the number of years, and r1 and r2 are the lag-1 autocorrelation of the two variables. The PDO index is calculated as the leading principle component from the empirical orthogonal function (EOF) analysis of North Pacific (poleward of 20°N) SST anomalies66, using the Met Office’s Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST)81 for observation and historical and SSP5-8.5 simulations of the 20 CMIP6 models. The long-term trend is removed at each grid point before conducting the EOF analysis. The regressions of OHC and wind anomalies on the negative PDO index (-PDO) are first calculated from each of the 20 CMIP6 models, and then the MMM regression maps are derived from the ensemble mean of their coefficients. The Oceanic Niño Index (ONI) is the 3-month running-mean detrended temperature anomalies in the Niño-3.4 region (5°N–5°S, 120°–170°W) based on version 5 of the Extended Reconstructed Sea Surface Temperature82, and the monthly ONI is averaged to the annual mean.

Time of Emergence (ToE)
A “signal-to-noise” method29,62 is adopted to compute the time of emergence (ToE) of the OHC trend. We use the piControl simulation, historical simulation, and projections under three emission scenarios (SSP1-2.6, SSP2-4.5, and SSP5-8.5) of 20 CMIP6 models to determine the “signal-to-noise” ratio (SNR). The noise is defined as the interannual standard deviation of the quadratically detrended83 OHC variability of piControl simulation in the last 100 years, representing the natural noise envelope. The anthropogenic signals for 1958–2100 are computed as the OHC differences between the historical plus three emission scenarios simulations and the time-mean piControl value. The ToE is defined as the time when the signal last exceeds and remains above a threshold of two times of noise (i.e., SNR > 2.0). The threshold for SNR is chosen as 2.0 because it represents a 95% confidence level of signal emergence84. For the regional ToEs, the noise and signal of selected regions are first averaged and then the time with SNR > 2.0 is searched. We also estimated the observed SNR to examine the possible spatial pattern of the observed anthropogenic OHC, as using observational data to calculate the ToE is infeasible, which relies on pre-industrial simulations from climate models. The signal is defined as the observed OHC trend (change) for 1958–2021, with the noise being one standard deviation of the detrended OHC.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this study, a mismatch between multi-decadal ocean heat storage patterns from observations and climate model simulations is found in the North Pacific and further attributed to the phase shifts of the Pacific Decadal Oscillation (PDO). Meanwhile, based on climate model simulations, the earlier emergence of ocean warming is identified in marginal seas along the North Pacific basin rim. Generally speaking, this is a comprehensive analysis examining contrast between observations and climate model simulations. However, there are some key concerns that the authors need to address before further consideration.\n\nThe authors should pay more attention to data quality issue and relevant uncertainties. (1) It might not be suitable to focus on multidecadal ocean heat content trends integrated over the upper 2000m, as ocean temperature measurements under 700m were very limited before the Argo project. Would 0-700m OHC reflect the same features? (2) Regarding the wind, previous studies have demonstrated that while there is large consistency among various wind products in terms of interannual and decadal variability, the long-term wind stress changes differ substantially. (3) the time of emergence analysis is purely based on climate model simulations. Would you also see similar spatial contrast from observational record?\n\nIt is well acknowledged that the PDO plays a dominant role in the North Pacific interannual and decadal variability. I believe this is the motivation for the authors trying to link ocean heat storage patterns to the PDO. However, the provided evidence is not solid enough to justify this statement. (1) The MNP (or the subpolar gyre) region is actually less affected by the PDO. The authors tried to explain using the delayed effect of Rossby wave propagation. However, the propagation feature is not clear at the latitude of MNP (Fig. S4; from Fig. 3e to 3f, also statistically insignificant). Figure 5d seems more like a reversal pattern from Lag0. (2) The PDO impacts on wind stress and Ekman pumping are mostly at lower latitudes (Figure 3c) than seen in the observed wind trend patterns (Figure 3a). (3) Would you also see a close relation between PDO index and MNP OHC from observational datasets (as Figure S2)? (4) Extended Fig.7: why the models are selected based on PDO index since 1998?\n\nIt seems that the HTFL run could reproduce, and thus the heat flux should dominate, the CMIP6 MMM patterns. However, this is not the way summarized in left figures of Extended Fig. 9. Moreover, it is not shown what’s the effect of wind in CMIP6 MMM – does it play little role or enhance the heat flux driven patterns?", "reviewer_2": "Comments on Duan et al. “Emergence of the North Pacific heat storage pattern delayed by decadal wind-driven redistribution”\nThe authors tray to decompose the temporal development of the ocean heat storage, which is modified by spatial processes, by general anthropogenic forcing, and by decadal climate variability. The claim that the asymmetric internal variability over the recent decades explain a distinct signal in the ocean heat content.\nEven if the approach is trivial in principle, it remains difficult to derive the contributions from the individual processes. This makes the topic interesting and has the potential to expand our understanding of the future spatial and temporal evolution of ocean heat content.\nThe authors use a variety of data and change the focus of the analysis frequently and very quickly. It is therefore not easy to follow the train of thought of the manuscript. Even if the individual analyses are justified, I would like to see a more stringent presentation of the results in order to facilitate readability and better understanding.\nHoever, my main concerns are based on the too superficial consideration of the quality of the model results. Unfortunately, the authors do not show to what extent the model data reflect a decadal variability of the PDO. It is therefore not clear which variability patterns are caused by the model deficiencies and which are caused by the decadal variability itself. This severely limits the confidence in the results of the study.\nI recommend a revision of the structure for a more stringent presentation and a clear distinction between the effects from model errors and those from the delay in warming patterns due to internal variability.\nComments:\nL 46-47: There may be more reasons for a difference between model and observations, e.g. the inability of simulating the natural internal variability (mainly in the ocean), or the initialization of the ocean state concerning the presentation of multi decadal and centennial variability in the models.\nL 52-54: This thesis is bold, there are many regions known where climate model results show strong biases in the mean state (North Atlantic Ocean circulation, intensity of ENSO, etc.).\nL 80: Why only 20 models, according to which criteria is the selection made?\nL93-94: Unfortunately, this is not comprehensible to the reader from the shown figure.\nL 99-110 (Fig 1c+d):\n- The figure caption names the content OHC. However, these are anomalies. I expect an explanation of the underlying value here.\n- The text deals with “average heat storage rate” and refers to the figure which present total ocean heat content anomalies. Please help the reader to relate these two values for a better understanding.\nL 111-119 (Fig 1c+d): Unfortunately, this paragraph does not help for understanding.\n- The scale of fig. 1c+d is different, which makes it difficult to compare the internal variability.\n- - The trend pattern of the MNP region shows different signs in the observation, which poses an additional problem for interpretation. It therefore appears to make more sense to only analyze mean values from homogeneous regions.\n- Keeping in might the temporal evolution of the observations in fig 1c (MNP), it might be that the difference between 1a und 1b in this region are random. The observed anomalies show very similar values at the beginning and at the end of the timeseries, with an strong negative trend until 1990 and strong positive afterwards. The MMM simulates a positive trend after 1980. It would have to be clarified what influence the length of the time series has on the derived result\n- Fig.1d NWTP, the observed anomalies are covered by strong internal variability, which was not captured by the models. Additionally, to this the variability increases over time. Here I would expect an explanation, why the models fail, and what process in the observations drives the gain in variability.\nL 126-150, fig. 2: Please indicate that between of figure 1a+b and figure 2a+b the scale differs by factor two. The model experiments presented here are interesting. However, the results obviously show clear problems in the model interaction between local energy balance and wind-driven ocean circulation. The breakdown of the processes in Figures 2b and 2c shows that the local equilibrium between the processes in the model deviates considerably from the observations, which is due to the large-scale uniform and common trends in MNP. This gives cause for concern that model errors dominate here and that the statement on the physical processes is strongly limited.\nL:187-189: Calculating a correlation between a low-passed filtered index with an annual index is really daring. Please provide a confidence level the method of calculating this level.\nL 196-205: Unfortunately, I could not follow the authors' reasoning based on the materials provided, please rephrase the paragraph and sharpen the train of thought.\nL 208-211: The authors would have to answer this statement or question before conducting the study. This would indicate whether the differences between model simulations and observations are due to the influence of decadal variability or to the deficiency of the model formulations.\nL 233-235: Please compare this event with earlier PDO phase shifts e.g. 1992, 2003\nL 265-269: The value of the study with regard to the influence of the PDO after LAG=9 years is clearly shown in the comparison of Figures 3f, 5b and 5d. If one assumes that the observations reach the highest confidence level, then the HYCOM simulation shows certain deficits. The results of the MMM do not correlate with the observations, especially in the high latitudes. Thus, a statement about the relationship between PDO and trend due to external forcing is highly problematic.", "reviewer_3": "Review of “Emergence of the North Pacific heat storage pattern 1 delayed by decadal wind-driven redistribution by Duan et al.” submitted to Nature Geosciences.\nGeneral comments:\nThis paper uses models and observations to try to understand interactions with the global warming trend and decadal variability in the North Pacific. The increase in ocean heat content projected with climate models for the last 60 years or so has not been observed. The authors use model simulations to try to separate the effect of wind forcing and thermal forcing with the purpose of interpreting the discrepancies. To begin they argue that the discrepancies can be associated with model inadequacies or associated with natural (decadal) variability processes. They decide to trust the models and use them to explain how natural variability is responsible for the differences. Natural variability through wave propagation is in conclusion argued to delay the excess warming in the North Pacific depicted in the climate models. I found this interpretation plausible, but I do not believe the authors have been able to demonstrate that this is the case mainly because of the concerns (that I list below). I would have expected a more critical interpretation of their results. Otherwise the paper reads well and the presentation of figures and organization is good.\nMain concerns:\n1) The linear approach is questionable\nThe approach used to separate the wind and thermal forcing assumes that these effects (dynamic and thermodynamic) can be decoupled in a linear fashion. Although modelers sometimes use this approach, in reality both dynamic and thermodynamic forcing may not linearly be separable. For instance, assuming wind stress variable but wind speed constant over time to isolate the effect of wind forcing is not realistic because all three air-sea fluxes (momentum flux, latent heat flux, and sensible heat flux) depend on wind speed. Thus, the conclusions that may be drawn from these modeling experiments may be questionable. I think that a more critical interpretation of the methodology should be included in the paper.\n2) Shortcoming in the climate models\nThe discrepancies between the modeled and observed warming may simply be related to inadequacies in the climate models in capturing coupled modes of natural variability. This is similar to the case of ENSO in the tropical Pacific where despite the efforts in obtaining oceanic and atmospheric observations both in the ocean and atmosphere over the years, the models have struggled to predict the main features of the phenomenon such us time of onset, intensity and duration.\nIn short, I think this paper does not convincingly explain why the climate models are not able to reproduce the observed heat storage patterns. Thus, the inferences drawn in the paper regarding long term heat storage are mostly speculative in my opinion. Perhaps climate models are not up to the task to address this problem at this time.\nSpecific Comments:\n- Lines 99-110 from the opening sentence seems to pertain to Figs. 1c and 1d but they do not really address these figures. They seem to just refer to Figs. 1a and 1b.\n- Line 447, the integral should be defined for z pointing upwards and from -2000 to 0. This is the convention in geophysical fluid dynamics."}, "author_response": "Point-by-Point Response to Reviewers’ Comments (our response in black)\nReviewer #1\nWe are grateful to Reviewer #1 for his/her effort in reading our paper and making the\ncomments. These comments/suggestions are very helpful for improving our work. We\nhave made necessary modifications in our paper to address them. Below we present our\nresponse and describe how the paper is revised.\n1. The authors should pay more attention to data quality issue and relevant uncertainties.\n(1) It might not be suitable to focus on multidecadal ocean heat content trends integrated\nover the upper 2000m, as ocean temperature measurements under 700m were very\nlimited before the Argo project. Would 0-700m OHC reflect the same features? (2)\nRegarding the wind, previous studies have demonstrated that while there is large\nconsistency among various wind products in terms of interannual and decadal\nvariability, the long-term wind stress changes differ substantially. (3) the time of\nemergence analysis is purely based on climate model simulations. Would you also see\nsimilar spatial contrast from observational record?\nThank you for pointing this out! We agree with the reviewer that data quality issues\nand relevant uncertainties should be better discussed. This has been addressed in the\nrevised paper.\n(1) In the revised version, we have added a supplementary figure of the multi-decadal\nOHC trends integrated over the upper 700 m (Supplementary Fig. 1), which are close\nto those of 0-2000 m (Fig. 1), despite detailed discrepancies. Similar to the 0-2000 m\nOHC, the 0-700 m OHC also shows a substantial increase in the MNP (170°E-155°W,\n40°-52°N), in stark contrast to the warming minimum in the NWTP (125°-180°E, 8°-\n18°N) (Supplementary Fig. 1a and Fig. 1a). In addition, the observed heat storage also\nshows a more complex structure with limited warming in the MNP and enhanced\nwarming in the NWTP (Supplementary Fig. 1b and Fig. 1b).\nSupplementary Fig. 1 | Ocean heat storage pattern and the evolution in key regions. a, b, Patterns\nof 0-700 m ocean heat content (OHC) trend (W m-2) of 1958-2021 based on the multi-model mean\n(MMM) of 20 CMIP6 historical and SSP5-8.5 simulations (a) and the Institute of Atmospheric\nPhysics ocean analysis (IAP) (b). Stippling indicates significant trends at the 95% confidence level\nbased on a Mann-Kendall test. The black boxes remark the regions with enhanced heat uptake: the\nmid-latitude North Pacific (MNP; 170°E-155°W, 40°-52°N) and the northwest tropical Pacific\n(NWTP; 125°-180°E, 8°-18°N). c, d, Regional OHC anomalies (GJ m-2, 1 GJ = 109 J) of MNP (c) and\nNWTP (d) derived from observations and CMIP6 simulations, relative to a 1960-1980 baseline.\nThe observations include IAP, WOA, Ishii, EN4, and ORA-S4. Shadings denote the one standard\ndeviation range.\n(2) Indeed, we need to examine the consistency in long-term wind stress change across\nvarious wind products. We have added a supplementary figure showing surface wind\ntrends of four atmospheric reanalysis products (JRA55, ERA5, NCEP-NCAR, and\nORA-S4) (Supplementary Fig. 2). They are consistent in showing critical features of\ntheir ensemble mean wind trend pattern (Fig. 3a), such as basin-scale trends of anti-\ncyclonic winds over the North Pacific since 1958, easterly winds and negative Ekman\npumping velocity in the NWTP, and strengthening westerlies north of 40°N.\nSupplementary Fig. 2 | Surface wind trends in different datasets. a, b, c, d, Trends of wind stress\n(arrows; Pa) and Ekman pumping velocity ω (color shading; 10-7 m s-1 yr-1) derived from JRA55\nE\n(a), ORA-S4 (b), ERA5 (c), and NCEP-NCAR (d). The ORA-S4 product over 1958-2017 combines\nthe ERA-40 for 1958-1988 and the ERA-Interim for 1989-2017. The ERA5 data of 1958 is derived\nfrom the ERA5 preliminary version.\n(3) We also agree that it is necessary to check the spatial pattern of the ToE in\nobservation. Although it is impossible to calculate the ToE using observational data as\nin climate models (pre-industrial simulation is necessary), we calculate the signal-to-\nnoise ratio (SNR) that mimics the ToE. The results based on the IAP and observation\nensemble (ensemble mean of IAP, WOA, Ishii, EN4, and ORA-S4) are shown in\nExtended Data Fig. 9c and 9d. The signal is defined as the observed OHC trend from\n1958 to 2021, with the noise being one standard deviation of the detrended OHC. The\nobserved SNR are generally in line with the ToE based on CMIP6 simulations in spatial\npattern (Extended Data Fig. 9). For instance, there are lower SNRs (corresponding to\nlater ToEs) in the MNP, KE, and NWTP regions than the Atlantic and Southern Oceans\nand regions along the North Pacific basin rim.\nExtended Data Fig. 9 | Time of emergence (ToE) for anthropogenic heat storage. Multi-model\nmedian ToE of the anthropogenic OHC trend under the SSP1-2.6 (a) and SSP5-8.5 (b) scenarios.\nThe ToE is the year when the climate change trend exceeds the range of the twofold natural\ninterannual standard deviation (Methods). Grey shading denotes the signal that has not emerged by\n2099. The boxes from west to east represent the NWTP, Kuroshio Extension (141°-165°E, 33°-40°N),\nthe Bering Sea (160°E-170°W, 52°-60°N), MNP (170°E-155°W, 40°-52°N), and Eastern boundary\n(132°-117°W, 30°-55°N). c, d, Same as a, b, but for signal-to-noise ratio (SNR) derived from IAP\nand the mean of five observations. The observations include IAP, WOA, Ishii, EN4, and ORA-S4.\nThe signal is defined as the observed OHC trend from 1958 to 2021, with the noise being one\nstandard deviation of the detrended OHC.\nWe have also revised the text to describe these new results (see the revised paper for\ndetails).\n2. It is well acknowledged that the PDO plays a dominant role in the North Pacific\ninterannual and decadal variability. I believe this is the motivation for the authors trying\nto link ocean heat storage patterns to the PDO. However, the provided evidence is not\nsolid enough to justify this statement. (1) The MNP (or the subpolar gyre) region is\nactually less affected by the PDO. The authors tried to explain using the delayed effect\nof Rossby wave propagation. However, the propagation feature is not clear at the\nlatitude of MNP (Fig. S4; from Fig. 3e to 3f, also statistically insignificant). Figure 5d\nseems more like a reversal pattern from Lag0. (2) The PDO impacts on wind stress and\nEkman pumping are mostly at lower latitudes (Figure 3c) than seen in the observed\nwind trend patterns (Figure 3a). (3) Would you also see a close relation between PDO\nindex and MNP OHC from observational datasets (as Figure S2)? (4) Extended Fig.7:\nwhy the models are selected based on PDO index since 1998?\nThese are important issues! We agree with the reviewer that we need to clarify the\nrelationship between North Pacific heat storage and PDO better. It is well established\nthat the PDO is characterized by decadal modulations in the Aleutian low and surface\nwinds over the North Pacific (e.g., Latif and Barnett 1994; Newman et al. 2016). This\nmotivates us to attribute the observed heat storage pattern to the PDO through wind-\ndriven heat redistribution.\n(1) Indeed, in the previous version, the propagation at the latitudes of the MNP is not\nclear enough. This is because we used a wide latitude range of 40°-55°N to obtain the\ntime-longitude plot. The phase speed of Rossby waves varies widely with latitude. In\nthe revised MS, we have narrowed the latitudinal range to 40°-52°N (Extended Data\nFig. 5). Following the suggestion from another reviewer, we have narrowed down the\nMNP region from (155°E-150°W, 40°-55°N) to (170°E-155°W, 40°-52°N). Now, we\ncan see a more consistent westward propagation feature, such as the dashed arrow.\nExtended Data Fig. 5 | Simulated impacts of the PDO on MNP. a, The normalized OHC anomaly\nof the MNP region (170°E-155°W, 40°-52°N). The linear trend of the OHC anomaly is removed. b,\nTime-longitude plot of SSH anomaly along the zonal band of 40°-52°N. To highlight interannual-\nto-interdecadal variations, the SSH anomalies are detrended and smoothed by LOWESS (locally\nweighted scatterplot smoothing) with a span width of 2 years. a, b are based on monthly HYCOM\n0.1° simulation of 1979-2021. The black arrow in (b) indicates the propagation of Rossby waves. c,\nNormalized monthly PDO index based on HadISST data.\nIn addition, we have added a supplementary figure showing the spatial evolution of the\nOHC derived from HYCOM 0.1° simulation lag -PDO from 0 to 9 years. As shown in\nSupplementary Fig. 4, negative anomalies near the eastern boundary north of 25°N\npropagate to the MNP and subsequently influence the KE system. Therefore, the\nrelatively unclear westward propagation may also be related to the influence of the KE\nsystem.\nThe westward propagation is further confirmed by the maximum lagged correlation\nbetween the zonal average SLA or OHC anomaly at 40°-52°N and the PDO index\nincreases from east to west (black circle-line in Figure R1). Significant simultaneous\npositive correlations exist mainly at the eastern boundary, indicating a quick response\nto local wind forcing. Figure 5a is not a reversal pattern from Lag = 0, as the fact that\nlag = 0 of the CMIP6 MMM differs from the observed lag = 0 pattern (Figure R2). For\nexample, in the NWTP, lag = 0 and lag = 9 show consistent negative anomalies.\nSupplementary Fig. 4 | Evolution of simulated impacts of the PDO on OHC in the North Pacific.\nLagged regression map of OHC on -PDO derived from HYCOM 0.1° simulation of 1979-2021, with\nthe former lagging the latter by 0 (a), 1 (b), 2 (c), 3 (d), 4 (e), 5 (f), 6 (g), 7 (h), 8 (i), and 9 years\n(j). Stippling indicates the fields exceeding the 95% confidence level.\nFigure R1. a, Maximum lagged correlation coefficient and lag year between the annual PDO index\nderived from HadISST and the meridional average (40°-52°N) SLA at different longitudes derived\nfrom the HYCOM 0.1° simulation from 1979-2021, with the latter lagging the former. The dashed\nline indicates the 95% confidence level. b, Same as a, but between annual PDO and OHC.\nFigure R2. a, Regression map of OHC anomalies derived from the MMM of 20 CMIP6 models on\n-PDO for 1958-2021. Stippling indicates the fields exceeding the 95% confidence level based on an\nF-test. The black boxes denote the NWTP and the MNP. b, c Same as a, but with the OHC anomalies\nlagging -PDO by 4 (b) and 9 years (c).\n(2) The observed surface wind trend contains both externally forced changes and\nnatural variability (PDO impacts). Discrepancies between Figures 3a (trend) and 3c\n(PDO regression) may largely reflect externally forced changes. Therefore, the PDO\nregression mainly explains the data-model difference between Figures 3a and 3b, as\nshown in Figure R3. Indeed, many observed features missed by CMIP6 MMM do show\nup in the PDO regression, such as easterlies in the NWTP and westerlies in the\nnortheastern subpolar Pacific, which has been stated in our MS.\nRegarding the enhanced changes in the tropics, comparing the zonal averages of Figures\n3a and 3c reveals that both have enhanced changes in Ekman pumping velocity in the\ntropics (Figure R4).\nFigure R3. Difference between the trends of wind stress (arrows; Pa yr-1) and Ekman pumping\nvelocity ω (color shading; 10-7 m s-1 yr-1) for 1958-2021 of the observation ensemble (mean of\nE\nJRA55, ERA5, NCEP-NCAR, and ORA-S4) and the CMIP6 MMM.\nFigure R4. a, Zonal-mean (120°E-90°W) trends of Ekman pumping velocity ω (10-8 m s-1 yr-1) for\nE\n1958-2021 based on observation ensemble (mean of JRA55, ERA5, NCEP-NCAR, and ORA-S4). b,\nSame as a, but for regressions of observed ω (10-6 m s-1) on the negative PDO index (-PDO) based\nE\non HadISST for 1958-2021. The shading denotes the one standard deviation range.\n(3) The MNP OHC derived from observational datasets also shows a lagged\nrelationship with the PDO, as that of HYCOM WND, as shown in Supplementary Fig.\n3. For example, the PDO shifted from the negative to the positive phase in the 1970s,\nand correspondingly there was a warming trend in the MNP throughout the 1980s and\n1990s (dashed line); subsequently, the opposite transition from the late 1990s to ~2010\nled to a cooling trend in the MNP that persisted until the late 2010s (dashed line).\nSupplementary Fig. 3 | Relationship between the evolution of OHC in the North Pacific and\nmodes of natural variability. a, Observed Oceanic Niño Index (ONI) and PDO indices and regional\nOHC in the NWTP derived from the ensemble mean of five observations and for 1958-2019. The\nobservations include IAP, WOA, Ishii, EN4, and ORA-S4. b Same as a, but for regional OHC of the\nMNP and the PDO index. The linear trends of regional OHC anomalies of the NWTP and MNP are\nremoved. All timeseries are normalized by their standard deviations. The PDO index is smoothed\nby LOWESS (locally weighted scatterplot smoothing) with a span width of 8 years. The monthly\nONI is downloaded from the Climate Prediction Center of the National Weather Service website and\naveraged to the annual mean. The blue dashed line represents the linear trend.\n(4) The year 1998 was chosen because the negative PDO phase started in 1998. It has\nhad a decisive influence on the global surface warming slowdown (Meehl et al. 2013;\nEngland et al. 2014; Nieves et al. 2015; Liu et al. 2016). This persistent negative phase\nis representative of the impact of the PDO phase shift on the long-term trends in surface\nwinds and OHC in the North Pacific. Therefore, models that agree with observation in\nthis negative phase are expected to better reproduce the observed trends in surface\nwinds and OHC in the North Pacific than other models.\nWe have added the descriptions to the text to reflect the above issues.\nLines 199-206: A regression onto the negative PDO index (-PDO) (Fig. 3c) shows\neasterlies in the NWTP and westerlies in the northeastern subpolar Pacific – key\nfeatures seen in observation but missed in CMIP6 MMM, suggesting that the PDO can\nlargely explain the data-model difference in wind trends (Fig. 3a compared to Fig. 3b).\nThe negative PDO also induces positive w anomalies in the western and central parts\nE\nof the 18°-30°N band – a feature of obvious data-model discrepancy that explains the\nabsence of the observed cooling trend near the western boundary in CMIP6 MMM (Fig.\n1a, b).\nLines 220-228: Lagged regressions show that the strengthened KE east of Japan is\nestablished ~4 years after the negative peak of PDO, and the MNP cooling takes ~9\nyears (Fig. 3e, f). The wind-driven OHC anomaly of the MNP, as measured by the WND\nrun of HYCOM, also shows a lagged negative correlation with the PDO index\n(Extended Data Fig. 3b). The PDO shifted from the negative to the positive phase\nduring the 1970s, and correspondingly there was a warming trend in the MNP\nthroughout the 1980s and 1990s; subsequently, the opposite transition from the late\n1990s through ~2010 led to a cooling trend of the MNP persisting through the late\n2010s. This lagged response of the MNP region to PDO is also seen in observation\n(Supplementary Fig. 3).\nLines 272-281: The response of KE to Rossby waves is successfully reproduced by a\n0.1° simulation of HYCOM (Methods; Extended Data Fig. 4) that can better resolve\nthe Kuroshio path and mesoscale eddies. As a result, this 0.1° simulation can well\nrepresent the heat redistribution driven by PDO winds, including the meridional dipole\neast of Japan with a lag time of ~4 years (Fig. 5a) and the cooling in the MNP with a\nlag time of ~9 years (Fig. 5b) to the negative peak of PDO. This ~9-year lag time\nreflects the slow propagation of upwelling Rossby waves between 40°-52°N (Extended\nData Fig. 5 and Supplementary Fig. 4) generated by positive w anomalies near the\nE\nwest coasts of North America (Fig. 3c). In addition, the delayed impact of the KE system\nalso contributes to the ~9-year lag time (Supplementary Fig. 4).\nLines 808-812: The year 1998 was chosen to reflect the 1997-1998 phase shift of PDO.\nThe persistent negative PDO phase of 1998-2014 had vital impacts on the overall\ntrends of OHC during the past several decades. Therefore, models that agree with\nobservation in this negative phase are expected to better reproduce the observed trends\nin surface winds and OHC in the North Pacific than other models.\n3. It seems that the HTFL run could reproduce, and thus the heat flux should dominate,\nthe CMIP6 MMM patterns. However, this is not the way summarized in left figures of\nExtended Fig. 9. Moreover, it is not shown what’s the effect of wind in CMIP6 MMM\n– does it play little role or enhance the heat flux driven patterns?\nThank you for pointing this out! We agree that the heat flux is the dominant factor in\nthe CMIP6 MMM patterns. We are sorry that the Extended Data Fig. 9 was not accurate.\nWe have removed this figure due to the number limit of the Extended Data figures.\nRegarding the role of wind in the CMIP6 MMM, as shown in Fig. 3b, trends of anti-\ncyclonic winds north of 30°N and westerly trends in the entire tropical Pacific can lead\nto OHC increase in the mid-latitudes and OHC decrease in the tropical Pacific. This\nwind-driven heat storage pattern aligns with the pattern driven by heat flux. Therefore,\nthe wind in the CMIP6 MMM can amplify the heat flux-driven patterns of heat storage.\nWe have added related descriptions.\nLines 153-154: This highlights the dominance of heat uptake in shaping the model-\nprojected heat storage.\nLines 191-193: In CMIP6 MMM, trends of anti-cyclonic winds are confined north of\n30°N, and westerly trends occupy the entire tropical Pacific (Fig. 3b), amplifying the\nmodel-projected heat storage pattern driven by heat uptake.\nReferences\nEngland MH, Mcgregor S, Spence P, et al (2014) Recent intensification of wind-driven circulation\nin the Pacific and the ongoing warming hiatus. Nat Clim Chang 4:222–227.\nhttps://doi.org/10.1038/nclimate2106\nLatif M, Barnett TP (1994) Causes of decadal climate variability over the North Pacific and North\nAmerica. Science (80- ) 266:634–637. https://doi.org/10.1126/science.266.5185.634\nLiu W, Xie SP, Lu J (2016) Tracking ocean heat uptake during the surface warming hiatus. Nat\nCommun 7:1–9. https://doi.org/10.1038/ncomms10926\nMeehl GA, Hu A, Arblaster JM, et al (2013) Externally Forced and Internally Generated Decadal\nClimate Variability Associated with the Interdecadal Pacific Oscillation. J Clim 26:7298–7310.\nhttps://doi.org/10.1175/JCLI-D-12-00548.1\nNewman M, Alexander MA, Ault TR, et al (2016) The Pacific Decadal Oscillation, Revisited. J\nClim 29:4399–4427. https://doi.org/10.1175/JCLI-D-15-0508.1\nNieves V, Willis JK, Patzert WC (2015) Recent hiatus caused by decadal shift in Indo-Pacific\nheating. Science (80- ) 349:532–535. https://doi.org/10.1126/science.aaa4521\nReviewer #2\nWe are very grateful to Reviewer #2 for his/her insightful comments. Apparently, this\nreviewer is of great expertise in this research topic, and these comments are super\nhelpful for improving our work. The reviewer also took a good effort in carefully\nreading our paper and provide detailed/specific suggestions for both scientific and\nwriting aspects. Below we describe our response and describe how the paper is revised.\n1. L 46-47: There may be more reasons for a difference between model and\nobservations, e.g. the inability of simulating the natural internal variability (mainly in\nthe ocean), or the initialization of the ocean state concerning the presentation of multi\ndecadal and centennial variability in the models.\nThank you for pointing this out! We agree that there are multiple sources of differences\nbetween climate models and observations. However, our focus here is to explain the\ncontrasts between the CMIP6 multi-model mean (MMM) and observations. In fact, in\nCMIP6 MMM, the effect of internal variability has been largely averaged out. The\nmotivation of the CMIP project is to simulate the response of the climate to external\nforcing (anthropogenic CO2, aerosols, ozone, etc). Attribution of the data-model\ncontrasts is of importance because the prediction of future climate change under\nanthropogenic forcing relies largely on the MMM (such as the IPCC reports).\nAll CMIP6 models integrate for at least 500 years under the pre-industrial CO2 level\n(piControl) and then begin the historical simulation from 1850 to 2014. In the historical\nsimulation, all models are forced by the observed external forcing (CO2, aerosols,\nozone, etc) (Eyring et al. 2016), but each model has its own internal variability owing\nto the nonlinearity of the climate system. The phase of internal variability (such as\nENSO and PDO) is different across these models and different from observation. For\nexample, in most models, there was no super El Niño occurred in 1997-1998. This does\nnot indicate that these models are incompetent in simulating ENSO. Even a “perfect”\nmodel cannot predict the phase (timing) of internal variability because of the\nnonlinearity of climate. Similarly, most models do not produce the negative phase of\nPDO during 1998-2014, which does not mean models fail to capture the PDO.\nGiven the different phases of internal variability across models, anomalies arising from\ninternal variability cancel each other when we average all models to obtain the MMM.\nAs such, CMIP6 MMM primarily represents the externally (mainly anthropogenic)\nforced change (Table R1). In this sense, even if the simulation of internal variability is\nseverely flawed in models, this is not assumed to contribute to the contrast between\nobservation and CMIP6 MMM in long-term change of OHC. It is the absence rather\nthan the bias of the simulated internal variability that may lead to discrepancies between\nobservation and CMIP6 MMM.\nGiven the above, the contrasts between observations and CMIP6 MMM can only arise\nfrom two origins: systematic model biases in simulating the externally forced heat\nstorage or the absence of internal variability in the CMIP6 MMM. If the former is true,\nthe predicted heat storage under anthropogenic climate change in the 21st century by\nCMIP6 models contains substantial errors. If the latter is true, we can still somehow\ntrust models in predicting externally forced heat storage in the future.\nThis reasoning is widely adopted among existing studies addressing the data-model\ndiscrepancies in historical change (e.g., Zhou et al. 2024).\nWe have comprehensively explained these issues in the revised paper.\nLines 44-46: Since prediction of climate change relies heavily on the ensemble of\nclimate models, correctly interpreting the discrepancies between MMM and\nobservation represents a vital scientific issue.\nLines 47-57: Although natural variability exists in individual models, the CMIP6\nMMM mainly represents externally forced change28. Therefore, the contrasts between\nobservations and CMIP6 MMM indicate either systematic model biases in simulating\nexternally forced heat storage or substantial impacts from natural variability on\nobserved heat storage. Although anthropogenic fingerprints have emerged in many\naspects of the ocean3,29,30, natural variability remains influential in the observed OHC\nchanges31,32. For instance, the persistent negative phase of the Pacific Decadal\nOscillation (PDO) since 1998 has led to heat pile-up in the western tropical Pacific via\nthe strengthened Pacific Walker Circulation33–35. Successful attribution of model-\nobservation discrepancies in the North Pacific heat storage to natural variability may\nenhance our confidence in the projected future change by climate models.\nLines 193-197: Some key features of the observed wind trend are missed, including the\neasterlies in the NWTP and westerlies in the subpolar North Pacific. These differences\nin surface winds, critically accounting for the model-observation discrepancies in heat\nstorage, likely arise from natural variability, given that the CMIP6 MMM is assumed\nto represent externally forced changes.\nLines 446-455: All CMIP6 models integrate for at least 500 years under the pre-\nindustrial CO level (piControl) and then from 1850 to 2014 to form the historical\n2\nsimulation. In the historical simulation, all models are forced by the observed external\nforcing (such as CO , aerosols, and ozone)5. Owing to the nonlinearity of the climate\n2\nsystem, the phase of internal variability (such as ENSO and PDO) is different across\nmodels and different from observation. Therefore, anomalies arising from internal\nvariability cancel each other when we average these 20 models to obtain the multi-\nmodel mean (MMM). As such, CMIP6 MMM primarily represents the externally\n(mainly anthropogenic) forced change, and the inter-model spread represents internal\nvariability and differences in model physics.\nTable R1. Major component of heat storage in observations and CMIP6 models\nDataset Major component of OHC trend or heat storage\nObservation Externally forced changes + Observational natural variability\nCMIP6 MMM Externally forced changes\nIndividual model of CMIP6 Externally forced changes + Simulated natural variability\n2. L 52-54: This thesis is bold, there are many regions known where climate model\nresults show strong biases in the mean state (North Atlantic Ocean circulation, intensity\nof ENSO, etc.).\nThank you for pointing this out! We agree that there are many regions known where\nclimate model results show strong biases in the mean state, such as the North Atlantic\nOcean circulation and intensity of ENSO (e.g., Wang et al. 2014; Yeh et al. 2018;\nXu et al. 2023). Therefore, even if model-observation discrepancies of heat storage pattern\nin the North Pacific can be successfully attributed to natural variability, the predictions\nof climate models are not considered reliable. We have revised the statement.\nLines 55-57: Successful attribution of model-observation discrepancies in the North\nPacific heat storage to natural variability may enhance our confidence in the projected\nfuture change by climate models.\n3. L 80: Why only 20 models, according to which criteria is the selection made?\nThank you for pointing this out! This is mainly dictated by data availability. These 20\nmodels provide all the variables we need for the analysis in all three emission scenarios\n(SSP5-8.5, SSP2-4.5, and SSP1-2.6) at the time our work began. For other models,\nsome variables of some scenarios were unavailable.\nLines 455-456: Selection of these 20 models is dictated by data availability for the\nthree emission scenarios when our work begins.\n4. L 93-94: Unfortunately, this is not comprehensible to the reader from the shown\nfigure.\nThank you for pointing this out! Sorry for that. We have revised the description.\nLines 96-97: CMIP6 MMM also shows weakened warming trends in these regions (Fig.\n1a).\n5. L 99-110 (Fig 1c+d):\n- The figure caption names the content OHC. However, these are anomalies. I expect\nan explanation of the underlying value here.\n- The text deals with “average heat storage rate” and refers to the figure which present\ntotal ocean heat content anomalies. Please help the reader to relate these two values for\na better understanding.\nThank you for pointing this out! We have changed the caption names of Figures 1c and\n1d to “OHC anomalies” and added the description. “Average heat storage rate” refers\nto the trend of regional average OHC anomalies. We apologize for not clearly stating\nthat the units of these values have been converted from J m-2 yr-1 to W m-2. We have\nmodified the text in the revised paper.\nLines 64-66: c, d, Regional OHC anomalies (GJ m-2, 1 GJ = 109 J) of MNP (c) and\nNWTP (d) derived from observations and CMIP6 simulations, relative to a 1960-1980\nbaseline.\nLines 104-116: We further examine the temporal evolution of OHC in key regions (Fig.\n1c, d). CMIP6 models show an average heat storage rate (i.e., trend of regional average\nOHC anomalies) of 2.20±2.08×107 J m-2 yr-1 (equivalent to 0.71±0.67 W m-2; ±\nindicates the one standard deviation range of 20 models) in the MNP of 170°E-155°W,\n40°-52°N, close in magnitude to that of the Southern Ocean (0.66±0.29 W m-2 in MMM\nfor 55°-33°S). In comparison, a much weaker rate of 0.07±0.08 W m-2 is obtained from\nobservation-based datasets for the MNP (Fig. 1c; ± indicates the one standard\ndeviation range of 5 datasets). Meanwhile, the NWTP shows enhanced warming of\n0.58±0.29 W m-2 in observations, one order stronger than the simulated rate of\n0.06±0.33 W m-2 in CMIP6 models (Fig. 1d). Interestingly, CMIP6 agrees with\nobservation in the OHC change of the entire North Pacific (120°E-80°W, 0°-60°N),\nwhich are 0.32±0.08 W m-2 and 0.32±0.18 W m-2, respectively, implying a heat\nredistribution over the North Pacific in the observed realization relative to the\nsimulated pattern.\nFig. 1 | Ocean heat storage pattern and the evolution in key regions. a, b, Patterns of 0-2000 m\nocean heat content (OHC) trend (W m-2) of 1958-2021 based on the multi-model mean (MMM) of\nanalysis (IAP)36 (b). Stippling indicates significant trends at the 95% confidence level based on a\nMann-Kendall test. The black boxes remark the regions with enhanced heat uptake: the\nmid-latitude North Pacific (MNP; 170°E-155°W, 40°-52°N) and the northwest tropical Pacific\n(NWTP; 125°-180°E, 8°-18°N). c, d, Regional OHC anomalies (GJ m-2, 1 GJ = 109 J) of MNP (c) and\nNWTP (d) derived from observations and CMIP6 simulations, relative to a 1960-1980 baseline. The\nobservations include IAP, WOA, Ishii, EN4, and ORA-S4. Shadings denote the one standard\ndeviation range.\n6. L 111-119 (Fig 1c+d): Unfortunately, this paragraph does not help for understanding.\n- The scale of fig. 1c+d is different, which makes it difficult to compare the internal\nvariability.\n- - The trend pattern of the MNP region shows different signs in the observation, which\nposes an additional problem for interpretation. It therefore appears to make more sense\nto only analyze mean values from homogeneous regions.\n- Keeping in might the temporal evolution of the observations in fig 1c (MNP), it might\nbe that the difference between 1a und 1b in this region are random. The observed\nanomalies show very similar values at the beginning and at the end of the timeseries,\nwith an strong negative trend until 1990 and strong positive afterwards. The MMM\nsimulates a positive trend after 1980. It would have to be clarified what influence the\nlength of the time series has on the derived result\n- Fig.1d NWTP, the observed anomalies are covered by strong internal variability,\nwhich was not captured by the models. Additionally, to this the variability increases\nover time. Here I would expect an explanation, why the models fail, and what process\nin the observations drives the gain in variability.\nThank you! These are important issues and have been addressed in the revised paper!\n(1) In the revised version, we adopt the same scale in Figures 1c and 1d for better\ncomparison of the internal variability.\n(2) Thank you for your nice suggestion! We have narrowed down the MNP region from\n(155°E-150°W, 40°-55°N) to (170°E-155°W, 40°-52°N), where most changes are\nnegative (Figure 1b and Extended Data Fig. 1).\n(3) To clarify the influences of the length of the timeseries on the model-observation\ndiscrepancies, we compare the heat storage patterns at different start and end years from\nobservations and CMIP6 MMM. As shown in Extended Data Fig. 2, although sensitive\nto the length of timeseries, model-observation discrepancies in the North Pacific are\nbasically stable when the time length exceeds 50 years (Extended Data Fig. 2a-c, k,\nl).\n(4) According to Figure 1d, the observed OHC anomalies in the NWTP show\nsignificant internal variability, and the effect of decadal variability is particularly\nnotable from the 1990s through the mid-2010s. During this period, the observed OHC\nof the NWTP was substantially elevated and exceeded the +1 standard deviation range\nof CMIP6 models. The reasons are twofold. First, given that the Interdecadal Pacific\nOscillation (IPO)/PDO negative phase of this period is likely an extreme case of natural\nvariability (Han et al. 2014; Bordbar et al. 2017; Chung et al. 2019; Geng et al. 2023),\n20 ensemble members are insufficient to capture it. Second, the CMIP6 models have\nsystematic shortcomings in simulating tropical Pacific climate, such as the mean state\nbiases in sea surface temperature (e.g., Zhang et al. 2023; Erickson and Patricola 2023).\nWe have revised the descriptions and added discussions.\nLines 86-92: The observation-based heat storage shows more regional structures (Fig.\n1b and Extended Data Fig. 1); for example, there are alternating warming and cooling\nregions in the North Pacific, particularly in the western basin. The NWTP (e.g., 125°-\n180°E, 8°-18°N) stands out with a warming maximum, whereas a “warming hole”\noccurs in the MNP with insignificant trends and is surrounded by significant warming\ntrends in marginal seas, such as the Bering Sea and the Gulf of Alaska.\nLines 117-130: The observed OHC anomalies show prominent interannual and\ndecadal variabilities in the MNP and NWTP (Fig. 1c, d). Meanwhile, we notice a large\ninter-model spread relative to the MMM change in the historical simulation of CMIP6\nbefore 2014, which also indicates significant influence from natural variability. The\neffect of decadal natural variability is particularly notable from the 1990s through the\nmid-2010s. During this period, the observed OHC of the NWTP was substantially\nelevated and exceeded the one standard deviation range of CMIP6 models, while that\nof the MNP did not rise significantly as models expected. Given that anomalies during\nthis period may represent an extreme case of natural climate variability, 20 models are\nlikely insufficient to capture them. Also, we cannot ignore the possibility that CMIP6\nmodels are incapable of reproducing such decadal variability. Due to the significant\ninfluence of natural variability, the contrasts between observations and CMIP6 MMM\nare sensitive to the time length but become stable when the time length exceeds 50 years\n(Extended data Fig. 2). In the following, we examine how decadal natural variability\naffects the observed North Pacific heat storage.\nLines 299-301: The PDO’s signatures in the NWTP are weaker in models than in\nobservations, linked to uncertain tropical Pacific wind changes (Extended Data Fig.\n6). This indicates a weaker coupling between the PDO and the tropical Pacific climate\nin models60.\nExtended Data Fig. 1 | Heat storage pattern in different datasets. Patterns of 0-2000 m ocean heat\ncontent (OHC) trend (W m-2) based on WOA (a), ORA-S4 (b), Ishii (c), and EN4 (d). Stippling\nindicates significant trends at the 95% confidence level. The black boxes remark the regions with\nenhanced heat uptake: the mid-latitude North Pacific (MNP; 170°E-155°W, 40°-52°N) and the\nnorthwest tropical Pacific (NWTP; 125°-180°E, 8°-18°N).\nExtended Data Fig. 2 | Sensitivity of OHC trend to the start and end year. a-f, Difference in OHC\ntrends starting from 1960 (a), 1965 (b), 1970 (c), 1975 (d), 1980 (e), and 1985 (f), derived from IAP\nminus CMIP6 MMM. g-l, Same as a-f, but for OHC trends ending in 1990 (g), 1995 (h), 2000 (i),\n2005 (j), 2010 (k), and 2015 (l).\n7. L 126-150: fig. 2: Please indicate that between of figure 1a+b and figure 2a+b the\nscale differs by factor two. The model experiments presented here are interesting.\nHowever, the results obviously show clear problems in the model interaction between\nlocal energy balance and wind-driven ocean circulation. The breakdown of the\nprocesses in Figures 2b and 2c shows that the local equilibrium between the processes\nin the model deviates considerably from the observations, which is due to the large-\nscale uniform and common trends in MNP. This gives cause for concern that model\nerrors dominate here and that the statement on the physical processes is strongly limited.\nThank you for pointing this out! We have added descriptions to clarify the difference\nin scale between Figures 1a,b and 2a,b,c.\nWe agree that it is necessary to discuss the nonlinear interaction between the heat\nbalance (controlled by heat fluxes, represented by HTFL) and ocean circulation change\n(driven by winds, represented by WND) in HYCOM experiments. A linear separation\nof heat flux change and wind-driven circulation change will inevitably lead to model\nerrors. We need to examine whether these errors dominate the heat storage patterns in\nHTFL and WND. To pursue this, we compare the sum of HTFL and WND\n(HTFL+WND) with that of CTRL. We find that the HTFL+WND broadly agrees with\nCTRL in OHC changes in the MNP and NWTP (grey and red lines in Fig. 2d, e) and\nspatial pattern of heat storage (Fig. R4). In addition, the OHC climatology in HTFL\nand WND, are consistent with that in CTRL (Fig. R5). This indicates that model errors\ndo not fundamentally affect the regional energy balance in HTFL and WND. Therefore,\nalthough the model errors arising from the nonlinear interaction are considerable, the\nidealized experiments of HTFL and WND still provide useful insights into the\nunderlying physical processes.\nWe have added the descriptions and supplementary materials.\nLines 163-169: The linear separation of the wind and heat flux forcing effects\ninevitably leads to model errors in HTFL and WND. To assess the potential influence\nof these artificial errors, we compare the sum of HTFL and WND (grey and red lines in\nFig. 2d, e) with CTRL. The sum of HTFL and WND broadly agrees with CTRL in both\nOHC changes of the MNP and NWTP (Fig. 2d, e) and the heat storage pattern (figure\nnot shown). Despite considerable model errors arising from the linear separation,\nHTFL and WND still provide useful insights into the underlying processes.\nFig. 2 | Effect of wind forcing on heat storage. a-c, Linear OHC trends of 1958-2019 derived from\nCTRL (a), HTFL (b), and WND (c) runs of the HYCOM 0.5° simulation. Stippling indicates\nsignificant trends at the 95% confidence level. The range of values for a-c is twice that of Fig. 1a,\nb. d, e, Regional OHC anomalies of the MNP (d) and NWTP (e) derived from CTRL, HTFL, and\nWND of HYCOM 0.5° simulation, referenced to those derived from IAP. Anomalies are obtained\nrelative to a 1960-1980 baseline. HTFL+WND denotes the sum of the OHC anomalies of HTFL and\nWND.\nFig. R4. Linear OHC trends of 1958-2019 derived from the sum of HTFL and WND runs of the\nHYCOM 0.5° simulation. Stippling indicates significant trends at the 95% confidence level.\nFig. R5. Climatological mean OHC of 1958-2019 derived from CTRL (a), HTFL (b), and WND (c)\nruns of the HYCOM 0.5° simulation.\n8. L 187-189: Calculating a correlation between a low-passed filtered index with an\nannual index is really daring. Please provide a confidence level the method of\ncalculating this level.\nThank you for your nice suggestion! We have added the confidence level and the related\ncalculation method. In the revised paper, the estimation of confidence level is based on\na two-tailed Student's t-test with the effective degree estimated as 𝑛\n(1−𝑟1𝑟2)\n, where n is\n(1+𝑟1𝑟2)\nthe data length in years, 𝑟 and 𝑟 are the autocorrelations of two variables (Bretherton\n1 2\net al. 1999). As such, the correlation coefficients of the NWTP heat content with the 8-\nyear low-pass PDO index and the unfiltered December-January-February Ocean Nino\nIndex (ONI) were -0.42 and -0.66, respectively, and were significant at 95% confidence.\nWe have revised related descriptions.\nLines 209-212: The NWTP heat content shows correlation coefficients of -0.42 and -\n0.66 with the 8-year low-passed PDO index and the unfiltered December-January-\nFebruary Oceanic Niño Index (ONI), respectively, both significant at the 95%\nconfidence level (Extended Data Fig. 3a).\nLines 494-497: The estimation of the confidence level of the correlation is based on a\ntwo-tailed Student's t-test using the effective degree of freedom81 estimated as 𝑛\n(1−𝑟1𝑟2)\n,\n(1+𝑟1𝑟2)\nwhere n is the data length in years, 𝑟 and 𝑟 are the autocorrelations of two variables.\n1 2\n9. L 196-205: Unfortunately, I could not follow the authors' reasoning based on the\nmaterials provided, please rephrase the paragraph and sharpen the train of thought.\nThank you for pointing this out! This paragraph explains why the observed OHC trend\n(Fig. 1b) and Lag = 0 regression onto the PDO (Fig. 3d) are different in the KE and\nMNP regions. Lagged changes to the negative PDO phase are critical in determining\nthe OHC trends of these two regions. We have rewritten the paragraph in the revised\npaper.\nLines 219-231: These discrepancies can be reconciled by considering delayed response\nof the ocean to wind forcing48. Lagged regressions show that the strengthened KE east\nof Japan is established ~4 years after the negative peak of PDO, and the MNP cooling\ntakes ~9 years (Fig. 3e, f). The wind-driven OHC anomaly of the MNP, as measured by\nthe WND run of HYCOM, also shows a lagged negative correlation with the PDO index\n(Extended Data Fig. 3b). The PDO shifted from the negative to the positive phase\nduring the 1970s, and correspondingly there was a warming trend in the MNP\nthroughout the 1980s and 1990s; subsequently, the opposite transition from the late\n1990s through ~2010 led to a cooling trend of the MNP persisting through the late\n2010s. This lagged response of the MNP region to PDO is also seen in observation\n(Supplementary Fig. 3). The lagged changes to the negative PDO phase are essential\nin shaping the OHC changes in the MNP and KE regions. Dynamical processes\nunderlying this delayed response demand in-depth understanding, which is the focus of\nthe following subsection.\n10. L 208-211: The authors would have to answer this statement or question before\nconducting the study. This would indicate whether the differences between model\nsimulations and observations are due to the influence of decadal variability or to the\ndeficiency of the model formulations.\nThank you for pointing this out! This is indeed an important question to answer. Climate\nmodels are known to have long-standing biases in simulated PDO and its impacts on\nthe ocean. As the reviewer pointed out, this affects the initialized prediction by models\nfor the coming decades. This has been clearly stated in the revised paper. In addition,\nwe have also detailed the biases in simulated PDO-induced winds and the OHC\nresponse to PDO of CMIP6 models in the subsection titled Regional Ocean Dynamics.\nLines 234-237: However, correctly simulating the PDO-induced variability, which is\ncritical for predicting regional OHC changes in the upcoming decades and their\nclimatic and environmental impacts, remains a challenging task for state-of-the-art\nmodels49–51.\nLines 292-301: We examine how the PDO-induced heat redistribution is represented\nin these models. Regressions using individual models (Methods) suggest that CMIP6\ncan capture some large-scale features of the PDO-induced changes in surface winds\nand OHC (Fig. 5c, d and Extended Data Fig. 6). The cooling of the MNP is seen at a\nlag of ~9 years, but the spatial structure does not correlate with the observed regression\n(Fig. 5d compared to Fig. 3f). Owing to coarse resolutions, CMIP6 models cannot\nrealistically represent the complex dynamics governing the lagged response to the PDO.\nThe PDO’s signatures in the NWTP are weaker in models than in observations, linked\nto uncertain tropical Pacific wind changes (Extended Data Fig. 6). This indicates a\nweaker coupling between the PDO and the tropical Pacific climate in models60.\nFor individual CMIP6 models, the phase and simulation skills of the PDO significantly\naffect the heat storage over the period 1958-2021. As such, none of the 20 CMIP6\nmodels has produced a historical heat storage pattern resembling the observation (not\nshown), which has been stated in our MS. However, as explained in our response to\nComment 1, the CMIP6 MMM contains primarily externally forced changes, and the\nPDO-induced changes have been largely removed by averaging all models into the\nMMM. It is the lack rather than the bias in the simulated PDO impacts that may lead to\ndiscrepancies in long-term OHC trends between observations and CMIP6 MMM.\nTherefore, biases in the simulated PDO mainly affect the near-term prediction when\nanthropogenic signals have not emerged.\n11. L 233-235: Please compare this event with earlier PDO phase shifts e.g. 1992, 2003\nThank you for your good suggestion! Earlier PDO phase shifts in ~1992 and 2003 from\nnegative to positive phase also led to the unstable state of KE after ~3 years (Fig. 4a\nand Extended Data Fig. 4). However, this unstable state persists for many years, unlike\nthe event of ~2014. In late 2017, the unstable KE returned to a stable state, dictated by\nan unexpected occurrence of Kuroshio large meander in August (Qiu et al. 2020). We\nhave revised the description.\nLines 258-265: The PDO shifted to a positive phase in ~2014, and correspondingly, an\nunstable dynamical state of the KE was established in early 2017, as indicated by the\nweakened southern recirculation (Fig. 4a). However, the Kuroshio large meander\noccurred in August 2017 - with the Kuroshio following an offshore meandering path in\nthe Shikoku basin south of Japan56 – compelled the KE to switch back to the stable state\nin late 2017 and remain stable since then52,57. In contrast, earlier phase shifts of the\nPDO, such as those in 1992 and 2003, also led to an unstable state of the KE after ~3\nyears, which persisted for many years (Fig. 4a and Extended Data Fig. 4).\nExtended Data Fig. 4 | Simulated impacts of the PDO on the Kuroshio Extension (KE). a, The\nnormalized OHC anomaly of the southern recirculation region (140°-165°E, 31°-36°N). b, Time-\nlongitude plot of SSH anomaly along the zonal band of 32°-34°N. a, b are based on monthly\nHYCOM 0.1° simulation of 1979-2021. The black arrows in (b) indicate the propagation of Rossby\nwaves. c, The normalized monthly PDO index based on HadISST data.\n12. L 265-269: The value of the study with regard to the influence of the PDO after\nLAG=9 years is clearly shown in the comparison of Figures 3f, 5b and 5d. If one\nassumes that the observations reach the highest confidence level, then the HYCOM\nsimulation shows certain deficits. The results of the MMM do not correlate with the\nobservations, especially in the high latitudes. Thus, a statement about the relationship\nbetween PDO and trend due to external forcing is highly problematic.\nThanks for pointing this out! We agree that HYCOM shows deficits, while CMIP6\nmodels cannot reproduce the lagged response to -PDO. The statement has been revised.\nLines 293-299: Regressions using individual models (Methods) suggest that CMIP6\ncan capture some large-scale features of the PDO-induced changes in surface winds\nand OHC (Fig. 5c, d and Extended Data Fig. 6). The cooling of the MNP is seen at a\nlag of ~9 years, but the spatial structure does not correlate with the observed regression\n(Fig. 5d compared to Fig. 3f). Owing to coarse resolutions, CMIP6 models cannot\nrealistically represent the complex dynamics governing the lagged response to the PDO.\nReferences\nBordbar MH, Martin T, Latif M, Park W (2017) Role of internal variability in recent decadal to\nmultidecadal tropical Pacific climate changes. Geophys Res Lett 44:4246–4255.\nhttps://doi.org/10.1002/2016GL072355\nBretherton CS, Widmann M, Dymnikov VP, et al (1999) The Effective Number of Spatial Degrees\nof Freedom of a Time-Varying Field. J Clim 12:1990–2009. https://doi.org/10.1175/1520-\n0442(1999)012<1990:TENOSD>2.0.CO;2\nChung ES, Timmermann A, Soden BJ, et al (2019) Reconciling opposing Walker circulation trends\nin observations and model projections. Nat Clim Chang 9:405–412.\nhttps://doi.org/10.1038/s41558-019-0446-4\nErickson NE, Patricola CM (2023) Future Projections of the El Niño—Southern Oscillation and\nTropical Pacific Mean State in CMIP6. J Geophys Res Atmos 128:\nhttps://doi.org/10.1029/2022JD037563\nEyring V, Bony S, Meehl GA, et al (2016) Overview of the Coupled Model Intercomparison Project\nPhase 6 (CMIP6) experimental design and organization. Geosci Model Dev 9:1937–1958.\nhttps://doi.org/10.5194/gmd-9-1937-2016\nGeng T, Jia F, Cai W, et al (2023) Increased occurrences of consecutive La Niña events under global\nwarming. Nature 619:774–781. https://doi.org/10.1038/s41586-023-06236-9\nHan W, Meehl GA, Hu A, et al (2014) Intensification of decadal and multi-decadal sea level\nvariability in the western tropical Pacific during recent decades. Clim Dyn 43:1357–1379.\nhttps://doi.org/10.1007/s00382-013-1951-1\nQiu B, Chen S, Schneider N, et al (2020) On the Reset of the Wind-Forced Decadal Kuroshio\nExtension Variability in Late 2017. J Clim 33:10813–10828. https://doi.org/10.1175/JCLI-D-\n20-0237.1\nWang C, Zhang L, Lee S-K, et al (2014) A global perspective on CMIP5 climate model biases. Nat\nClim Chang 4:201–205. https://doi.org/10.1038/nclimate2118\nXu Y, Li T, Xu M, et al (2023) Evaluation of the Pacific Decadal Oscillation from 1901 to 2014 in\nCMIP6 models. Clim Res 90:1–15. https://doi.org/10.3354/cr01711\nYeh SW, Cai W, Min SK, et al (2018) ENSO Atmospheric Teleconnections and Their Response to\nGreenhouse Gas Forcing. Rev Geophys 56:185–206. https://doi.org/10.1002/2017RG000568\nZhang Q, Liu B, Li S, Zhou T (2023) Understanding Models’ Global Sea Surface Temperature Bias\nin Mean State: From CMIP5 to CMIP6. Geophys Res Lett 50:1–11.\nZhou W, Leung LR, Lu J (2024) Steady threefold Arctic amplification of externally forced warming\nmasked by natural variability. Nat Geosci. https://doi.org/10.1038/s41561-024-01441-1\nReviewer #3\nWe are very grateful to Reviewer #3 for his/her effort in reading our paper and making\nthe comments. Apparently, this reviewer is of great expertise in this research topic, and\nthese comments are super helpful for improving our work. We have made necessary\nmodifications in our paper to address them. Below we present our response and describe\nhow the paper is revised.\nMain concerns:\n1. The linear approach is questionable\nThe approach used to separate the wind and thermal forcing assumes that these effects\n(dynamic and thermodynamic) can be decoupled in a linear fashion. Although modelers\nsometimes use this approach, in reality both dynamic and thermodynamic forcing may\nnot linearly be separable. For instance, assuming wind stress variable but wind speed\nconstant over time to isolate the effect of wind forcing is not realistic because all three\nair-sea fluxes (momentum flux, latent heat flux, and sensible heat flux) depend on wind\nspeed. Thus, the conclusions that may be drawn from these modeling experiments may\nbe questionable. I think that a more critical interpretation of the methodology should be\nincluded in the paper.\nThank you for pointing this out! Indeed, it is necessary to discuss the nonlinear\ninteraction between dynamic (ocean circulation change driven by wind stresses\nrepresented by WND) and thermodynamic (heat balance controlled by heat fluxes\nrepresented by HTFL) effects in HYCOM experiments. The linear separation of wind\nand heat flux forcing will inevitably lead to nonlinear effects or model errors. We need\nto examine whether these nonlinear effects or model errors dominate the heat storage\npatterns in HTFL and WND. To pursue this, we compare the sum of HTFL and WND\n(HTFL+WND) with that of CTRL. We find HTFL+WND broadly agrees with CTRL\nin OHC changes in the MNP and NWTP (grey and red lines in Fig. 2d, e) and the spatial\npattern of heat storage (Fig. R4). Furthermore, the OHC climatology in HTFL and\nWND are consistent with that in CTRL (Fig. R5). This indicates that nonlinear effects\nand model errors do not fundamentally affect the underlying physical processes of OHC\nin HTFL and WND. Therefore, although the nonlinear interaction and model errors are\nconsiderable, the idealized experiments of HTFL and WND still provide useful insights\ninto the underlying physical processes. We have added the related descriptions and\nmodified the interpretation of the linear separation methodology.\nLines 163-169: The linear separation of the wind and heat flux forcing effects\ninevitably leads to model errors in HTFL and WND. To assess the potential influence\nof these artificial errors, we compare the sum of HTFL and WND (grey and red lines in\nFig. 2d, e) with CTRL. The sum of HTFL and WND broadly agrees with CTRL in both\nOHC changes of the MNP and NWTP (Fig. 2d, e) and the heat storage pattern (figure\nnot shown). Despite considerable model errors arising from the linear separation,\nHTFL and WND still provide useful insights into the underlying processes.\nLines 468-480: Two sensitivity runs are conducted to explore the effects of different\ndrivers on the ocean heat storage pattern. Only wind stress adopts realistic daily fields\nin wind stress run (WND), while all the other forcings, such as wind speed and\nradiations, are fixed to the 1958 fields as in spin-up. In the WND, OHC trends primarily\narise from heat redistribution due to wind stress-driven changes in ocean circulation\nprocesses. The heat flux run (HTFL) only adopts realistic SWR and LWR radiations, 2-\nm air temperature, specific humidity, and wind speed from 1958 to 2019 but fixes wind\nstress and precipitation to the 1958 fields. As such, changes in these fields of HTFL\naffect OHC trends mainly by altering the heat fluxes (SWR, LWR, sensible, and latent\nheat fluxes), representing the role of heat-driven processes or ocean heat uptake. This\nlinear separation of the wind and heat flux forcing effects inevitably leads to nonlinear\ninteractions or model errors, but HTFL and WND still provide useful insights into the\nunderlying physical processes.\nFig. 2 | Effect of wind forcing on heat storage. a-c, Linear OHC trends of 1958-2019 derived from\nCTRL (a), HTFL (b), and WND (c) runs of the HYCOM 0.5° simulation. Stippling indicates\nsignificant trends at the 95% confidence level. The range of values for a-c is twice that of Fig. 1a,\nb. d, e, Regional OHC anomalies of the MNP (d) and NWTP (e) derived from CTRL, HTFL, and\nWND of HYCOM 0.5° simulation, referenced to those derived from IAP. Anomalies are obtained\nrelative to a 1960-1980 baseline. HTFL+WND denotes the sum of the OHC anomalies of HTFL and\nWND.\nFig. R4. Linear OHC trends of 1958-2019 derived from the sum of HTFL and WND runs of the\nHYCOM 0.5° simulation. Stippling indicates significant trends at the 95% confidence level.\nFig. R5. Climatological mean OHC of 1958-2019 derived from CTRL (a), HTFL (b), and WND (c)\nruns of the HYCOM 0.5° simulation.\n2. Shortcoming in the climate models\nThe discrepancies between the modeled and observed warming may simply be related\nto inadequacies in the climate models in capturing coupled modes of natural variability.\nThis is similar to the case of ENSO in the tropical Pacific where despite the efforts in\nobtaining oceanic and atmospheric observations both in the ocean and atmosphere over\nthe years, the models have struggled to predict the main features of the phenomenon\nsuch us time of onset, intensity and duration.\nIn short, I think this paper does not convincingly explain why the climate models are\nnot able to reproduce the observed heat storage patterns. Thus, the inferences drawn in\nthe paper regarding long term heat storage are mostly speculative in my opinion.\nPerhaps climate models are not up to the task to address this problem at this time.\nThank you for pointing this out! We agree that climate models have shortcomings in\ncapturing coupled modes of natural variability, such as ENSO and PDO. CMIP6 models\nare known to have long-standing biases in simulating the PDO and its impacts on the\nOHC, and the subsection entitled Regional Ocean Dynamics has been specifically\ndevoted to elucidating these shortcomings in the revised MS. Indeed, the surface wind\nchanges induced by the PDO are generally inaccurate (Extended Data Fig. 6), and\nalmost none of the 20 CMIP6 models realistically reproducing the spatial structure of\nthe observed regression (Extended Data Fig. 6 compared to Fig. 3c). Specifically, the\nPDO’s signatures in the NWTP in CMIP6 models exhibit a general inconsistency with\nobservations, and the central position of the anti-cyclonic wind anomalies in the mid-\nlatitudes indicates a significant misalignment with observations. This is likely related\nto the inability of climate models to successfully resolve key atmospheric processes of\ngenerating the PDO, such as underestimating the atmospheric noise forcing (e.g.,\nNewman et al. 2016; Cassou et al. 2018). Furthermore, there are evident biases in\nCMIP6 models regarding the OHC response to PDO-induced winds. For instance, the\nspatial structure of the CMIP6 MMM at a lag of 9 years is significantly inconsistent\nwith the observed regression, although cooling anomalies are also present in the MNP\n(Fig. 5d compared to Fig. 3f). Owing to coarse resolutions, CMIP6 models cannot\nrealistically represent the complex dynamics governing the OHC response to the PDO.\nThis is also supported by our coarse resolution (0.5°) simulation of HYCOM, which\nfails to fully reproduce the observed heat storage in the North Pacific (Fig. 2a), even\nwith prescribed reanalysis winds. For individual CMIP6 models, the phase and\nsimulation skills of the PDO significantly affect the heat storage over the period 1958-\n2021. As such, none of the 20 CMIP6 models has produced a historical heat storage\npattern resembling the observation (not shown), which has been stated in our MS.\nExtended Data Fig. 6 | Impacts of the PDO on simulated winds. Regressions of wind stress (Pa)\nand ω (10-6 m s-1) on the negative PDO index (-PDO) derived from CMIP6 MMM and 20 CMIP6\nE\nmodels for 1958-2021.\nHowever, our focus here is to explain the contrasts between the CMIP6 multi-model\nmean (MMM) and observations. In fact, in CMIP6 MMM, the effects of internal\nvariability have been averaged out. The motivation of the CMIP project is to simulate\nthe response of the climate to external forcing (anthropogenic CO2, aerosols, ozone,\netc). All CMIP6 models integrate for at least 500 years under the pre-industrial CO2\nlevel (piControl) and then begin the historical simulation from 1850 to 2014. In the\nhistorical simulation, all models are forced by the observed external forcing fields (CO2,\naerosols, ozone, etc) (Eyring et al. 2016), but each model has its own internal variability\nowing to the nonlinearity of the climate system. The phase of internal variability (such\nas ENSO and PDO) is different across these models and different from observation. For\nexample, in most models, there was no super El Niño occurred in 1997-1998. Even a\n“perfect” model cannot predict the phase (timing) of internal variability because of the\nnonlinearity of climate. Similarly, most models do not produce the negative phase of\nPDO during 1998-2014.\nGiven the different phases of internal variability across models, anomalies arising from\ninternal variability cancel each other when we average all models to obtain the MMM.\nTherefore, CMIP6 MMM primarily represents the externally (mainly anthropogenic)\nforced change (Table R1). In this sense, even if the simulation of internal variability is\nseverely flawed in models, this is not assumed to contribute to the contrast between\nobservation and CMIP6 MMM in long-term change of OHC. It is the absence rather\nthan the bias of the simulated internal variability that may lead to discrepancies between\nobservation and CMIP6 MMM. By contrast, the simulation of the PDO significantly\ninfluences the initialized decadal climate prediction for the upcoming decades when\nanthropogenic signals have not emerged, which has been stated in our MS.\nWe have comprehensively explained these issues in the revised paper.\nLines 47-50: Although natural variability exists in individual models, the CMIP6\nMMM mainly represents externally forced change28. Therefore, the contrasts between\nobservations and CMIP6 MMM indicate either systematic model biases in simulating\nexternally forced heat storage or substantial impacts from natural variability on\nobserved heat storage.\nLines 193-197: Some key features of the observed wind trend are missed, including the\neasterlies in the NWTP and westerlies in the subpolar North Pacific. These differences\nin surface winds, critically accounting for the model-observation discrepancies in heat\nstorage, likely arise from natural variability, given that the CMIP6 MMM is assumed\nto represent externally forced changes.\nLines 234-242: However, correctly simulating the PDO-induced variability, which is\ncritical for predicting regional OHC changes in the upcoming decades and their\nclimatic and environmental impacts, remains a challenging task for state-of-the-art\nmodels49-51. Even with prescribed reanalysis winds, our 0.5° simulation of HYCOM fails\nto fully reproduce the observed heat storage in the MNP (Fig. 2a); for instance, the\nsimulated cooling in the MNP is stronger in amplitude, broader in spatial extent, and\nshifted to lower latitudes in CTRL compared to that in observation. These discrepancies\nprobably arise from complex regional ocean dynamics that are not properly represented\nby coarse-resolution models - such as the “bimodal” variability of the KE48,52.\nLines 292-301: We examine how the PDO-induced heat redistribution is represented\nin these models. Regressions using individual models (Methods) suggest that CMIP6\ncan capture some large-scale features of the PDO-induced changes in surface winds\nand OHC (Fig. 5c, d and Extended Data Fig. 6). The cooling of the MNP is seen at a\nlag of ~9 years, but the spatial structure does not correlate with the observed regression\n(Fig. 5d compared to Fig. 3f). Owing to coarse resolutions, CMIP6 models cannot\nrealistically represent the complex dynamics governing the lagged response to the PDO.\nThe PDO’s signatures in the NWTP are weaker in models than in observations, linked\nto uncertain tropical Pacific wind changes (Extended Data Fig. 6). This indicates a\nweaker coupling between the PDO and the tropical Pacific climate in models60.\nLines 392-394: Climate models tend to underestimate the PDO variability and its\nimpacts over the historical period50. As a result, none of the 20 CMIP6 models has\nproduced a historical heat storage pattern resembling the observation (not shown).\nLines 446-455: All CMIP6 models integrate for at least 500 years under the pre-\nindustrial CO2 level (piControl) and then from 1850 to 2014 to form the historical\nsimulation. In the historical simulation, all models are forced by the observed external\nforcing (such as CO2, aerosols, and ozone)5. Owing to the nonlinearity of the climate\nsystem, the phase of internal variability (such as ENSO and PDO) is different across\nmodels and different from observation. Therefore, anomalies arising from internal\nvariability cancel each other when we average these 20 models to obtain the multi-\nmodel mean (MMM). As such, CMIP6 MMM primarily represents the externally\n(mainly anthropogenic) forced change, and the inter-model spread represents internal\nvariability and differences in model physics.\nTable R1. Major component of heat storage in observations and CMIP6 models\nDataset Major component of OHC trend or heat storage\nObservation Externally forced changes + Observational natural variability\nCMIP6 MMM Externally forced changes\nIndividual model of CMIP6 Externally forced changes + Simulated natural variability\nSpecific Comments:\n- Lines 99-110 from the opening sentence seems to pertain to Figs. 1c and 1d but they\ndo not really address these figures. They seem to just refer to Figs. 1a and 1b.\nThank you for pointing this out! It was somewhat misleading in our previous MS\nbecause we converted the unit from J m-2 yr-1 to W m-2. We have modified related\ndescriptions.\nLines 104-116: We further examine the temporal evolution of OHC in key regions (Fig.\n1c, d). CMIP6 models show an average heat storage rate (i.e., trend of regional average\nOHC anomalies) of 2.20±2.08×107 J m-2 yr-1 (equivalent to 0.71±0.67 W m-2; ±\nindicates the one standard deviation range of 20 models) in the MNP of 170°E-155°W,\n40°-52°N, close in magnitude to that of the Southern Ocean (0.66±0.29 W m-2 in MMM\nfor 55°-33°S). In comparison, a much weaker rate of 0.07±0.08 W m-2 is obtained from\nobservation-based datasets for the MNP (Fig. 1c; ± indicates the one standard\ndeviation range of 5 datasets). Meanwhile, the NWTP shows enhanced warming of\n0.58±0.29 W m-2 in observations, one order stronger than the simulated rate of\n0.06±0.33 W m-2 in CMIP6 models (Fig. 1d). Interestingly, CMIP6 agrees with\nobservation in the OHC change of the entire North Pacific (120°E-80°W, 0°-60°N),\nwhich are 0.32±0.08 W m-2 and 0.32±0.18 W m-2, respectively, implying a heat\nredistribution over the North Pacific in the observed realization relative to the\nsimulated pattern.\n- Line 447, the integral should be defined for z pointing upwards and from -2000 to 0.\nThis is the convention in geophysical fluid dynamics.\nThank you！We have corrected it.\nReferences\nCassou C, Kushnir Y, Hawkins E, et al (2018) Decadal Climate Variability and Predictability:\nChallenges and Opportunities. Bull Am Meteorol Soc 99:479–490.\nhttps://doi.org/10.1175/BAMS-D-16-0286.1\nEyring V, Bony S, Meehl GA, et al (2016) Overview of the Coupled Model Intercomparison Project\nPhase 6 (CMIP6) experimental design and organization. Geosci Model Dev 9:1937–1958.\nhttps://doi.org/10.5194/gmd-9-1937-2016\nNewman M, Alexander MA, Ault TR, et al (2016) The Pacific Decadal Oscillation, Revisited. J\nClim 29:4399–4427. https://doi.org/10.1175/JCLI-D-15-0508.1"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I appreciate that the authors have spent great effort to address my previous comments. I'm now happy to see the publication of this work.", "reviewer_2": "After reviewing all the reviewers' comments and the authors' responses, I am completely satisfied with the scientific content and quality of the manuscript.\nThe answers to the questions in particular show that although there is an interesting discussion between all those involved, this is only marginally reflected in substantial changes to the manuscript.\nI therefore suggest that the authors be given a certain amount of freedom to determine the focus of the argument independently.\nIn conclusion, I consider the manuscript to be mature and am in favor of publishing the current version.", "reviewer_3": "General comments:\nI am satisfied by how the authors have addressed by comments and have no other comments at this time. I am pleased by their openness to accept and deal with my criticism and am thankful for that. I think the reviewed paper does offer a insightful view at this difficult problem and therefore it makes a positive contribution."}, "author_response": ""}]

Decision letter