Four ocean reanalysis datasets, including the Ocean Reanalysis System 5 (ORAS5), Global Ocean Reanalysis and Simulations version 4 (GLORYS), Simple Ocean Data Assimilation version 3 (SODA3), and Global Ocean Data Assimilation System (GODAS), having high temporal and spatial resolutions, were selected for this study. The objective analysis product from the Institute of Atmospheric Physics (IAP) at the Chinese Academy of Sciences was also used for comparisons. All data were considered as monthly means, and the brief descriptions of the datasets are listed in Table 1.

Following Liu et al. (2023), the OHC in a grid box was calculated using the seawater potential temperature by the following equation: O H C = A ∫ z 1 z 2 ρ C p T ′ d z , (1) where A is the horizontal area of the grid box, ρ is the density of seawater, C p is the specific heat of the seawater, T ′ is the seawater potential temperature anomaly relative to the reference period of 2006–2018, and z ₁ and z ₂ are the lower and upper layer boundaries of the integration, respectively. Since the changes in the density and specific heat of the seawater are small, they were selected as constants of values ρ=1027 kgm⁻³ and C p = 4187 Jkg⁻¹K⁻¹, respectively (Jian et al., 2022). The OHC integrated within a domain is then divided by the domain area to obtain the unit of J/m².

The OHCT is defined as the first derivative of the OHC and is calculated by the central difference as follows: O H C T i = O H C i + 1 − O H C i − 1 t i + 1 − t i − 1 , (2) where i denotes a particular month, t _i-1 and t _i+1 denote the time instances one month before and after the month i, respectively.

The northward OHT at a given latitude θ can be estimated as per Hall and Bryden (1982): O H T θ = R c o s θ ∫ φ 1 φ 2 ∫ z 2 z 1 ρ C p T − T r v d z d φ , (3) where R is the radius of the Earth (in m), φ is the longitude, and v is the northward component of the current velocity. The reference temperature T r is set to the multiannual (2006–2018) domain mean temperature in accordance with Lee et al. (2004) and Verdy et al. (2023).

Similarly, the eastward OHT at a given longitude φ is estimated as O H T φ = R ∫ θ 1 θ 2 ∫ z 2 z 1 ρ C p T − T r u d z d θ , (4) where u is the eastward component of the current velocity.

Although the WPWP boundaries vary with time, a fixed region (10°S to 20°N, 110°E to 165°E) was selected (Figure 1) as the study area based on the time mean SST above 28.5°C and distance from the Australian continent to facilitate the study on changes in the OHC and OHT. The reference temperature T r is 7.1°C for ORAS5, 7.0°C for GLORYS, 7.2°C for SODA3, and 6.9°C for GODAS based on the selected domain. The multivariate ENSO index (MEI, version 2), which combines oceanic and atmospheric variables (Wolter and Timlin, 2011), was used to investigate the relationship between the WPWP and ENSO. Unless otherwise stated, the two-tailed t-test was used for the significance test of the Pearson correlation coefficient in this study at a significance level of 0.01, and the Mann–Kendall test was used to test the trend at a significance level of 0.01 (Hipel and Mcleod, 1994).

The 12-month running mean WPWP OHC and OHCT time series from the five datasets are plotted in Figure 2. The OHC is shown on the left column, and the OHCT is shown on the right column. The time series ranges from 1980 to 2020, except for the GLORYS data that spans from 1993 to 2020. The time series of the OHC and OHCT from the five datasets are very consistent over 0–300 m (Figures 2A, E); they show slight warming before 1990, following which the OHC dropped for a few years owing to the eruption of Mount Pinatubo in May 1991, and again followed by the El-Niño cycle around 1997/98, when the ocean heat was transported from the WPWP to the eastern Pacific along the thermocline (Xu et al., 2020) and accompanied by a large drop of the OHC in the WPWP. All datasets show overall increasing trends (2.22–3.15 W/m² in Table 2) during 1993–2010, followed by decreasing trends after 2010 and a large drop in the OHC around the 2010/11 ENSO event. The interannual variations were modulated by the ENSO, and the OHC and MEI have opposite phases. The correlation coefficients between the OHC and MEI are all less than or equal to −0.87 for the two selected time periods (Table 3). Below 300 m, the consistency among the five datasets is not as good as that over 0–300 m (Figures 2B, C), with the GODAS dataset showing a larger discrepancy over the others. There are large differences between the datasets in the range of 700–2000 m before 2005, and the trends from SODA3 and IAPv4 are much smaller than those from ORAS5 and GODAS, although the variability is quite consistent (high correlation coefficients in Table 2). The consistency improved after 2005, when more observations were available with the implementation of the Argo floats (Figure 2C). The GODAS has a spurious high value in January 2018 and merits further investigation. Similar to the variation in the OHC integrated from 0 to 300 m, the variability and trend of the OHC over 0–2000 m follow the MEI well and are mainly derived from the OHC in the top 0–300 m based on a comparison of the variability magnitudes. All datasets show significant warming trends over the different layers between 1993 and 2010, except for the SODA3 and IAPv4 data, which show insignificant trends in the depths between 700 and 2000 m. Most datasets show significant cooling trends in the layer between 0 and 700 m over 2010–2020, while the trends between 700 and 2000 m are significantly positive.

Similar to the OHC, the variability of the OHCT integrated from 0 to 300 m was also very consistent among the datasets (Figure 3E), as confirmed by their correlation coefficient of r = 0.99 over both 1980–2020 and 1993–2020 (Table 2). The OHCT is closely related to the MEI, and the correlation coefficients are all significant; meanwhile, the OHCT leads MEI by approximately 7 months (Table 3). It is interesting to note that the OHCT integrated over 300–700 m is also related to the ENSO cycle and leads the MEI by approximately 7–9 months depending on the selected datasets and time periods. However, the MEI leads the OHCT by 4–6 months between 700 and 2000 m. The OHCT integrated from 0 to 2000 m is similar to that integrated from 0 to 300 m, which dominates the variability of both OHC and OHCT in the top 2000 m of the WPWP.

To further investigate the causes of the OHC changes in the WPWP, the OHTs through four sections along the WPWP boundaries (Figure 3) were calculated, whose results are plotted in Figure 4. The y-axis in Figure 4 is scaled differently for each panel to show the OHT variability through each boundary. The plots with the same scaling are given in Supplementary Figure S1 for reference. The OHT was integrated from 0 to 2000 m, and the lines are 12-month running means. The eastward and northward OHTs are defined to be positive. The OHTs through the eastern, southern, western, and northern boundaries of the WPWP are denoted as OHT_E, OHT_S, OHT_W, and OHT_N, respectively. The OHT_N through the northern boundary is all positive (Figure 4A), implying that heat flows away from the WPWP. The multiannual mean OHT_N values are 0.95 (ORAS5), 1.17 (GLORYS), 1.22 (SODA3), and 1.13 PW (GODAS) (Table 4). The lines are broadly spread, and the mean standard deviations (STDs) are 0.15 PW over 1980–2020 (three datasets) and 0.58 PW over 1993–2020 (four datasets), with all of them being in phase. The correlation coefficients ( r ) between the ORAS5 and other datasets are listed in Table 5, and all the values are significant. The correlation coefficients between the OHT_N and MEI are small, but most of them are significant (Table 5).

The OHT_E through the eastern boundary shows larger amplitude variability and larger spread (Figure 4B) than those of the OHT_N, with the multidataset mean STDs being 0.58 PW (19080–2020) and 0.93 (1993–2020). The multiannual mean OHT values (ORAS5: −0.86 PW, GLORYS: −2.44 PW, SODA3: −0.86 PW, GODAS: −0.60 PW) are all negative, indicating inflow of heat to the WPWP. Therefore, the WPWP mainly gains heat from the eastern boundary and loses heat from the northern boundary. The OHT_E time series are in phase with each other and lead the MEIs by 4–7 months except for the GODAS data, where the MEI leads the OHT_E (Table 5). Compared with OHT_N and OHT_E, the multiannual mean OHT_S values are −0.12 (ORAS5), −0.12 (GLORYS), −0.22 (SODA3), and −0.005 PW (GODAS). The transport direction is southward, and the amount of heat transported is small compared to those of the north and east sections (Table 4). The STDs are 0.23 PW (1980–2020) and 0.24 PW (1993–2020). The OHT_W through the western boundary is also weak (Figure 4D), and its mean values (ORAS5: 0.05 PW, GLORYS: −0.03 PW, SODA3: −0.11 PW, GODAS: −0.02 PW) are quite small. The variability between the datasets is more or less consistent, and this is confirmed by the correlation coefficients ( r ) between ORAS5 and the other datasets, as listed in Table 5. However, the time series have large spreads compared with their interannual variabilities. The MEI leads both OHT_S and OHT_W. The net zonal OHT inflow to the WPWP is all positive, and the net meridional OHT inflow to the WPWP is all negative (Table 4). Therefore, the WPWP gains heat from the zonal flows, mainly from the flow through the eastern boundary. It is observed from Table 4 that the net OHT inflow in GLORYS is positive, while the values in the others are negative, with the value being less negative in ORAS5 than the other two. These differences between the estimates may be attributed to the increment in data assimilation, which must be taken into account to reduce uncertainties in future studies.

The net OHT inflow to the WPWP is divided by the section area to obtain the unit of W/m². The net OHT and OHCT in the WPWP from the four datasets are plotted together in Figure 5; it is seen that there is good agreement in the variability and that the correlation improves after 1995 (Table 6); for example, the correlation coefficient r = 0.88, 0.88, and 0.90 for ORAS5 over 1980–2020, 1993–2020, and 2005–2020, respectively. These results imply that although there are energy flux exchanges at the surface and at 2000 m, the heat change rate in the top 2000 m of the WPWP is dominated by horizontal heat transport. The reason for the large net OHT deviation between GODAS and others before 1995 may be attributed to the use of different surface forcings, as described in Section 2, meriting further investigation.

To track the paths of the heat flows to and away from the WPWP, the multiannual mean (2006–2018) OHTs through the latitude–depth cross sections at different longitudes are plotted in Figure 6 using ORAS5 data. There are several concentrated flow regions, which have been marked as A, B, C, and D. The positive OHT is the eastward flow, and the dashed black line is the position of the thermocline defined as the location of the 20°C isotherm. The thermocline moves toward the surface from west to east, as expected. The heat flow patterns in all panels from west to east are similar, but the concentrated flow regions move upward with elevating thermocline (Figures 6A–H). The westward south and north equatorial currents in regions A and B (Figure 6I) are due to the trade winds; the eastward countercurrent at the equator is due to the low sea-surface pressure (Hu et al., 2015) and is concentrated around the depth of the thermocline. The countercurrent is also closely related to the ENSO cycle consisting of the neutral, El-Niño, and La-Niña phases (Zhao et al., 2013; Chen et al., 2016; Webb, 2018; Wijaya and Hisaki, 2021). It has been reported that the latitude locations of these equatorial currents could change with warming climates (Zhai et al., 2014).

The time series of the OHTs through the marked regions (Figure 6I) are plotted in Figure 7, together with the MEI time series. The maximum (minimum) lead–lag correlation coefficient ( r ₁ ) between the OHT and MEI is also displayed with the leading months in the parentheses, where a negative leading month indicates that the OHT leads MEI. The selected region A includes all westward OHT (negative) grid points, while region C includes all eastward OHT (positive) grid points in a rectangular area (8°S to 3°N, 300–600 m, at 165°E); regions B (10°N to 18°N, 0–300 m, at 165°E) and D (3°S to 8°N, 0–300 m, at 165°E) include all westward OHT (negative) grid points in the corresponding rectangular areas. The OHTs show overall decreasing trends in region A following the MEI and a significant correlation coefficient of r ₁ = 0.79 with the OHT leading MEI by 4 months (Figure 7A). Both regions B and D have westward OHTs, and their magnitude variabilities are smaller than that for region A. The correlation coefficients between the OHTs, MEIs, and leading months are different in these two regions. In region C, the correlation coefficient is only 0.33, with the OHT leading MEI by 4 months. The results from the other three datasets show similar patterns (not shown here). The OHT_E through the eastern boundary of the WPWP is also plotted for comparison of the OHTs through the marked regions. The correlation coefficient r ₂ between the OHT and OHT_E is calculated and displayed; there are significant correlations for regions A, C, and D, but the correlation for B is insignificant. The corresponding scatter plots between the OHT and OHT_E are shown in Figure 8 for regions A, B, and C along with their correlation coefficients.

Similarly, the multiannual mean (2006–2018) OHTs through the longitude–depth cross sections at different latitudes are plotted in Figure 9; it is seen that the western boundary current of the Pacific Ocean dominates the heat flow and that the depth for the concentrated heat flow increases from 20°N to 35°N. The structure of the boundary current varies with latitude. There are northward and southward flows at 20°N, and the southward flow on the left is weaker than the northward flow. The flow structure from west to east becomes complicated at higher latitudes. There is a northward flow along the western boundary, followed by a southward flow and then a northward flow again from 25°N to 35°N. The northward flow is the Kuroshio current, which is mainly the diversion of the westward north equatorial current that turns eastward at approximately 35°N, and this can be seen in Figures 9E, F where the concentrated flow disappears. The time series of the OHTs through regions E and F are plotted in Figure 10 along with the MEI time series. The southward heat flow in region E shows an overall decreasing trend, and the correlation coefficient of 0.32 between the OHT and MEI is significant, with the OHT leading MEI by 14 months.