Decadal changes in the basin-wide heat budget of the mid-latitude North Pacific Ocean

Research has indicated that the meridional heat transport (MHT) in the North Pacific Ocean (NP) across 24°N increased in the 1980s and 1990s, resulting in different heat distributions: the ocean heat content (OHC) increased in the 1980s, while the net surface heat release was strengthened in the 1990s; however, the reasons for these differences remained unclear. The authors revisited the investigation of the heat budget in the NP mainly using hydrographic observations to understand why the heat distribution was different between these 2 decades and extend the analysis to the 2010s. The OHC in the upper 700 m north of 24°N and east of 137°E exhibited sharp increases around 1990 as well as in the 2010s, while it was nearly stable in between. The northward retreat of the subarctic gyre boundary coincided with the spin-up of the subtropical gyre in the late 1980s, thereby allowing warm anomaly from the subtropics to propagate northeastward. Meanwhile, the concurrent weakening of the wintertime westerlies resulted in the suppression of the surface heat loss in the western NP. In contrast, the southward shift of the subarctic front suppressed the OHC rise, despite the MHT increase in the late 1990s. In the 2010s, unprecedented warming occurred in the eastern NP. The MHT estimation based on hydrographic observations indicates that the net surface release must have been suppressed since the MHT did not increase; however, the latest atmospheric reanalysis datasets failed to reproduce this.


1. Introduction
The Earth's climate has changed drastically since the midtwentieth century. Reports have indicated that the global mean temperatures of both the atmosphere and ocean have exhibited significant warming trends in recent decades (e.g. Levitus et al. 2012). A recent study indicated using a combination of repeat hydrography and Argo data that the global ocean is warming at all pressure levels, with local warming maxima being evident at the surface, and at 1000 and 4200 m depth (Desbruyères et al. 2017). According to their results, warming in the upper 700 m accounted for 44% of the global full-depth ocean heat content (OHC) trend. In the upper layer, warming trends can be seen in most of the globe, while cooling trends are evident in the subtropics in the Pacific Ocean and north of 30°N in the North Atlantic Ocean (Desbruyères et al. 2017;Johnson and Lyman 2020). The redistribution of OHC is crucial for spatiotemporal changes in the upper-layer temperature. OHC changes are caused by two different processes: sea surface heating and the convergence of horizontal heat transport. Both need to be considered to understand ocean temperature variations.
The meridional heat transports (MHTs) of the ocean and atmosphere are among the most important factors determining the Earth's climate, and several studies have elucidated their climatologies (e.g. Hsiung et al. 1989;Trenberth and Caron 2001;Ganachaud and Wunsh 2003).

3
Capturing temporal meridional transport variations across zonal sections in wide basins based on in situ observations is challenging: nevertheless, some studies have demonstrated such variations. For example, in the North Atlantic Ocean, an extensive mooring array has been deployed along 26.5°N to monitor the Atlantic meridional overturning circulation, and a time series of transport across this latitude has been obtained since 2004. The temporal mean MHT across 26.5°N in 8.5 years was 1.25 PW with a standard deviation of 0.36 PW, but it fell to about 0.5 PW in the winters of 2009/10 and 2010/11 (McCarthy et al. 2015. Cunningham et al. (2013) indicated that the OHC in the upper 2000 m between 26.5°N and 41.0°N in the Atlantic Ocean largely decreased during 2010-2011, and that both the surface cooling and MHT divergence were responsible for the OHC change.
For the North Pacific Ocean (NP), Roemmich et al. (2001) evaluated the MHT using a high-resolution expendable bathythermograph (XBT)/expendable conductivity-temperaturedepth (XCTD) transect around 22°N. The mean MHT across this transect was 0.83 ± 0.12 PW during 1993-1999. They indicated that the interannual MHT range was more than 0.40 PW. Kawai et al. (2008) (hereafter, K08) analyzed observational OHC data, surface heat fluxes of various atmospheric reanalysis datasets, and numerical model data, and concluded that the 5-year running mean of the MHT across 24°N grew by 0.1-0.2 PW in the 1980s and 1990s. Both the Kuroshio and Ekman transports were responsible for the increase in the 1980s, and the Kuroshio predominated regarding the change in the 1990s. While the MHT increase in the 1980s resulted in warming in the mid-latitudes of the NP, the heat release from the ocean to the atmosphere was enhanced in the 1990s. An increase in the OHC was not observed north of 24°N in the NP in the 1990s. The reason for the differences between the 1980s and 1990s, however, is yet to be explained.
Sea surface temperature (SST) in the NP has a dominant basin-wide variation on a decadal scale, known as the Pacific Decadal Oscillation (PDO) (e.g. Newman et al. 2016). SST becomes lower (higher) in the central NP and higher (lower) in the eastern NP along the North America continent in the positive (negative) phase of the PDO. OHC changes also resembled the PDO pattern, especially in the central NP, but not in the Kuroshio recirculation (KR) region (Figs. 2 and 3 of K08). And, the effects of the global warming are eliminated in the extraction of the PDO pattern. This study examines the OHC variation independently of the PDO index. Newman et al. (2016) indicated that the PDO was the result of a combination of different physical processes on different time scales; changes in surface heat fluxes and Ekman transport, ocean memory, and decadal changes in the Kuroshio-Oyashio system. Although we have a different perspective and do not pay our attention to the periodicity, our interpretations also consist of these three factors.
In this paper, the authors revisit the investigation of the heat balance among the OHC, MHT, and surface heat flux in the NP mainly based on hydrographic observations. This study aims to (1) elucidate why the heat distribution in the NP was different between the 1980s and 1990s, and (2) extend the analysis of the heat budget to the 2010s and assess the surface heat flux of several datasets. Section 2 describes the data and methods used in this study. The spatiotemporal variations of the OHC are described in Section 3. Estimated MHT and net surface heat flux are presented in Section 4 and 5. Section 6 discusses the differences in heat budget between the 1980s and 1990s. Section 7 summarizes our findings.

Gridded historical temperature and salinity dataset
We used a historical gridded dataset produced by the Japan Meteorological Agency (JMA), version 7.3.1, (Ishii et al. 2017) to calculate the OHC. The dataset had a horizontal resolution of 1° × 1° and 28 levels from the surface to 3000 m depth, and was produced from observations collected by the World Ocean Data 2013, the Global Temperature-Salinity Profile Program, and the Argo project. COBE-SST2, an objectively analyzed SST dataset by Hirahara et al. (2014), provided sea surface observations. The advanced depth bias correction by Ishii and Kimoto (2009) was applied to the XBT/mechanical bathythermograph (MBT) data, although no depth bias correction was used in version 6.2 by Ishii et al. (2006). The analyzed values on 1° × 1° grids were deviations from background fields, optimally estimated using a three-dimensional variational minimization technique. The background was the climatology of version 6.7 by Ishii and Kimoto (2009) plus a prescribed trend, which was defined as the leading empirical orthogonal function computed from 5-year 5° means.

137°E repeat hydrography
The JMA has been conducting hydrographic surveys along 137°E south of 34°N since 1967-once a year until 1971, and basically twice a year (i.e., in winter and summer) since 1972. The history and details of these surveys have been summarized in a review paper by Oka et al. (2018). Their meridional resolution is 1° south of 30°N (south of 32°N until 2002) with finer resolutions north of that. CTDs have been used since 1990, but Nansen bottles and reversing thermometers were used until 1989, when the vertical resolution was 10 m near the surface and approximately 250 m below 1000 m depth. The JMA has produced and released optimally interpolated data with a vertical resolution of 1 dbar and a unified meridional resolution of 20' at 31-34°N, 30' at 30-31°N, and 1° south of 30°N (Kawakami et al. 2022). Profiles have been vertically interpolated with the Akima spline (Akima 1970). The vertical coverage of the interpolated data was limited to 1250 m depth in and before 1991; hence, throughout this study geostrophic current is calculated relative to 1250 dbar. The Gibbs-SeaWater Oceanographic Toolbox is used for calculating the geostrophic currents (McDougall and Barker 2011).

Surface heat and momentum fluxes
Although the air-sea heat flux integrated over a wide basin significantly differs among datasets, different datasets may exhibit consistent temporal changes (K08). To detect consistent decadal changes over the NP, we used three kinds of latest atmospheric reanalysis dataset: the Japanese 55-Year Reanalysis (JRA-55) provided by the JMA (Kobayashi et al. 2015), the Climate Forecast System Reanalysis/Climate Forecast System version 2 (CFSR/CFSv2) provided by the National Centers for Environmental Prediction (NCEP) (Saha et al. 2010(Saha et al. , 2014, and the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis version 5 (ERA5) (Hersbach et al. 2020). The horizontal resolutions of JRA-55, CFSR, and ERA5 were 1.25°, 0.31°, and 0.25°, respectively. We mainly show the spatial distributions of ERA5, which have the finest resolution; however, we note that there are few differences in the spatial patterns of the fluxes among the different datasets.
In addition to the atmospheric reanalysis products mentioned above, we used the latest satellite-based surface heat flux product, i.e., the third-generation Japanese Ocean Flux Data Sets with Use of Remote-Sensing Observations (J-OFURO3) (Tomita et al. 2019). The main products of J-OFURO3 are surface turbulent heat fluxes using multi-satellite retrievals of surface winds, surface specific humidity, and sea surface temperature. J-OFURO3 also includes surface radiative flux data from the International Satellite Cloud Climatology Project (ISCCP; Rossow and Schiffer 1991) and the Clouds and the Earth's Radiant Energy System project (CERES; Loeb et al. 2018). The first version of J-OFURO3 involved a problem with long-term variations due to a systematic difference between ISCCP and CERES, but version 1.1 minimized the problem and extended the period. The spatial resolution was 0.25°; here we used monthly mean products for 1988-2017.

Calculation of meridional heat transport
Because the northern boundary of the NP is almost closed, we assumed that the MHT across a basin section is nearly equal to the sum of the OHC change rate and the surface net heat flux north of the section, and can be written in the following form: where ρ, C p , v, θ, and Q net are the density of seawater, specific heat of seawater at constant pressure, meridional velocity, potential temperature, and net surface heat flux, while α and γ represent the western and eastern ends of the section, respectively. Shallow circulation is responsible for most of the heat transport in the NP (Roemmich et al. 2001), and the observation coverage at depths below 700 m was quite poor before the Argo era (Ishii et al. 2017). Hence, we ignored variations below 700 m depth. Equation (1) approximates the convergence of the internal energy transport into the target region as long as the mass is conserved: we refer to this as "heat transport" by convention. We used 137°E observations to evaluate the Kuroshio transport, and the target region of the NP was partitioned by the 24°N and 137°E lines. Equation (1) can be rewritten for the target region: where u is zonal velocity. The locations of stations α, β, and γ were, respectively, the following: 34°N, 137°E; 24°N, 137°E; and 24°N, 107°W.
is the 137°E section, which includes the Kuroshio. We neglected heat transport through any straits on the border of the region.
The MHT consists of the Kuroshio transport V K , Ekmanlayer transport V E and interior transport V I (1) where Θ K , Θ I , and Θ E are the transport-weighted mean temperatures in the Kuroshio region, the Ekman layer, and the interior, respectively, τ x and τ y are the zonal and meridional wind stress components along the section. f is the Coriolis parameter, and θ e is the mean temperature of the Ekman layer, which we assume that approximates SST. For the Kuroshio region , we calculate the zonal transport instead of the meridional one; however, we expediently refer to the sum of the meridional transport across 24°N and the zonal transport across 137°E as "meridional" transport. We note that the axis of the Kuroshio is located north of 30°N around 137°E in most cases, and the Kuroshio region defined in this formulation is too wide, including a part of the North Equatorial Current (see Section 4). The entire Japan Sea and East China Sea were excluded from the target region, but the Sea of Okhotsk was included ( Fig. 1). Previous studies indicated that the transports through the Bering and Tsushima Straits are approximately 1.0 (Woodgate 2018) and 2.1 Sv (Kida et al. 2021), respectively (1 Sv = 10 6 m 3 s −1 ). The net transport across the section of + needs to be equal to the inflow through the Tsushima Strait minus the outflow through the Bering Strait, V R , assuming that river runoff and transport through the Strait of Tartary are negligible: If V R is included in V I , Eq. 3 can be transformed into the following form: The first and second terms on the right-hand side represent the net heat transports by the Kuroshio and Ekman currents, respectively.
V K and V I are obtained by integrating the geostrophic velocity above 1250 dbar (Sects. 2.1.1 and 2.1.2), V E is calculated from the wind stress of ERA5, and V R is assumed to be constant and 1.1 Sv. The sum of these four components does not become exactly zero (see Sect. 4); this imbalance could be due to the uncertainty of the historical gridded data or the limited number of samples of the 137°E repeat hydrography. To close the mass budget, we uniformly add a correction to the geostrophic velocity in the interior each year so that the imbalance almost vanishes (Roemmich et al. 2001). V R is regarded as a part of the imbalance in Eq. 4.

Overall temporal change
The OHC variations are remarkable east of the Philippines, with local minimum values being evident around 24°N (Fig. 1). In the mid-latitudes, the temporal OHC variation is most noticeable east of Japan around the beginning of the Kuroshio Extension (KE), and large along the Kuroshio, KE, subarctic boundary (SAB), and KE northern branch (KENB), which is also known as the Kuroshio bifurcation front (Kida et al. 2015). The area with a large standard deviation extends from the KE toward the Gulf of Alaska; this large variability in the subarctic frontal zone corresponds to variations in the density-compensated temperature anomaly (i.e., spiciness component), as indicated by Taguchi et al. (2017).
The global OHC in the upper 700 m has been increased nearly monotonously after 1970 (e.g. Ishii and Kimoto 2009), and the same was true for the whole NP (Fig. 2a). The OHC anomaly north of 24°N and east of 137°E, however, exhibits pronounced differences (Fig. 2b-c); including areas west of 137°E made little difference. The anomaly decreased until the mid-1980s, followed by an abrupt rise around 1990, and then remained nearly constant until the beginning of the 2010s. A sharp increase occurred again in the mid-2010s. Compared with the dataset of Ishii et al. (2006) (Fig. 1 of  K08), the OHC anomaly in the 1970s was smaller in the updated dataset, probably due to the correction of the XBT/ MBT data (Ishii and Kimoto 2009), but it did not change significantly from the 1980s to the 2000s. Hanawa (2002, 2003) indicated that a regime shift occurred Fig. 1 Standard deviation of annual-mean OHC above 700 m depth during 1955-2020. Dashed lines are the borders of the target region of this study. The target region is partitioned with the 24°N and 137°E lines. While the Japan Sea is excluded, the Sea of Okhotsk is included. See the text for the abbreviations of region names during 1988/89, corresponding to the OHC jump in Fig. 2b. The mid-latitude SST change during the 1988/89 regime shift was not associated with tropical variations, unlike the previous regime shifts (Yasunaka and Hanawa 2003). This suggests that the reason for the OHC rise around 1990 was within the middle and high latitudes, independent of tropical phenomena such as the El Niño-Southern Oscillation and Indian Ocean Dipole. Conversely, the effect from the tropics in 2014/15 indicated by Di Lorenzo and Mantua (2016) and Joh and Di Lorenzo (2019) resembled the regime-shift SST patterns before 1988/89 that were associated with tropical variations (Yasunaka and Hanawa 2003).
Hereafter, we focus on 5-year running means to smooth short interannual variations. The OHC jump from the late 1980s to the early 1990s was mainly attributed to warming in the central NP around 40°N, and the KR region, as already indicated by K08 (Fig. 3a). The OHC largely decreased around 1990 south of 24°N, especially in the western tropics. The decrease in the tropics and increase in the mid-latitudes canceled each other, resulting in the suppression of the OHC over the entire NP around 1990 (Fig. 2a). Conversely, in the 2010s, noticeable warming was evident in the eastern NP along North America, the southwestern Bering Sea, and the western NP, especially southeast of Hokkaido, although the central NP exhibited cooling (Fig. 3b). Both OHC jumps involved warming in the KR region in common. While the subsurface layer between about 100-m and 400-m depth, where subtropical mode water (STMW) is distributed, was especially warmed around 1990 (Fig. 3c), warming was mainly evident near the surface and the layer around 200-m depth was not warmed in the 2010s (Fig. 3d). Their spatial warming patterns over the entire NP between these periods were quite different, and it does not appear that an identical phenomenon repeated periodically. Another important point is that the 24°N line is a node that separates the OHC variations in the mid-latitudes from those in the tropics.

OHC jump around 1990
The temporal evolution of OHC from the late 1980s to the early 1990s is shown in Fig. 4. The OHC started to increase in the coastal area southeast of Japan and in the KE. The warm anomaly propagated to the KENB around 38°N, 160°E, extending east-northeast toward the Gulf of Alaska. Meanwhile, the KR region gradually warmed. The average salinity in the upper 700 m increased in these regions, especially around the subarctic boundary, during the same period (Fig. 5), suggesting advection of subtropical warm, saline waters.
For the period from the late 1980s to the early 1990s, the warm anomaly extended to more than 600 m depth at 30-35°N west of 180° (Fig. 3c), indicating that the KE axis shifted northward. Previous studies have indicated that the  Fig. 1). c Change rate of the anomaly in panel b. Thin and bold lines represent annual and 5-year means, respectively KE becomes stable and its front moves northward (e.g. Seo et al. 2014;Qiu et al. 2014;Joh and Di Lorenzo 2019), and that the surface temperature and velocity tend to be higher around the KENB when the KE front shifts northward (Sugimoto 2014). The positive OHC anomaly spread along the KE and KENB regions during this period, as a result of the stabilization of the KE. Furthermore, Yasuda and Kitamura (2003) used a numerical model to show that a dynamical adjustment to the wind stress fields in the central NP caused an increase in the Kuroshio and KE transport in the late 1980s with a lag of three years (their Fig. 8), leading to an increase in the STMW core temperature through the advection of heat by the KE. The negative wind stress curl strengthened in the central NP from 1983 to 1989, three years before the period shown in Fig.3a, 4 and 5, which is consistent with Yasuda and Kitamura's indication (Fig. 6). The increase in Kuroshio transport during this period was confirmed in both the numerical model (K08) and routine 137°E in situ observations (see Sect. 4). Overall, the spin-up of the subtropical gyre was responsible for the OHC changes in the KR, KE, and KENB regions. Figure 6 is also close to the 2-month lagged regression coefficient between the KE index and Ekman pumping velocity anomaly shown by Qiu et al. (2014), suggesting that the preceding KE unstable situation in the mid-1980s triggered the gyre spin-up.
The warming extending from the KENB to the Gulf of Alaska is accounted for by the northward shift of the subarctic front. The OHC within box B (Fig. 3a) corresponds well to the latitude of the maximum horizontal gradient of temperature at 100 m depth with a correlation coefficient of 0.771 (p < 0.01), and to the North Pacific Index (NPI; Trenberth and Hurrell 1994) with a correlation coefficient of 0.599 (p < 0.01) on an annual mean basis (Figs. 7 and 8a).This is consistent with the indication in a previous study that SST in the central NP was correlated with the NPI in the cold season (Newman et al. 2016). The simultaneous correlation coefficient between the annual mean latitude of the maximum gradient and the NPI was 0.385  (Kuroda et al. 2015). The abrupt OHC increase in the target region was attributed to nearly synchronized OHC increases in regions A and B, which were observed only during this period. In the other periods, the OHC changes nearly canceled each other out (Fig. 8a).

OHC jump in the 2010s
While the OHCs of boxes A and B exhibited no or slight increases in the 2010s, the rest of the NP warmed unprecedentedly (Fig. 8a). The KE was unstable in the late 2000s  1986 (i.e., 1984-1988) to 1992 (i.e., 1990-1994) Fig. 6 Difference in 5-year means of surface wind stress curl of ERA5 between 1983ERA5 between (i.e., 1981ERA5 between -1985ERA5 between ) and 1989ERA5 between (i.e., 1987ERA5 between -1991. Bold black line denotes the zero line of the climatological mean. Unit is 1 × 10 −7 N m −3 and stabilized in the 2010s (e.g. Joh and Di Lorenzo 2019; Qiu et al. 2020). The northward shift of the KE axis caused a high-temperature anomaly east of Japan west of 160°E (box F in Figs. 3b and 8b), but did not expand east of 160°E, obstructed by the southward shift of the subarctic front (Fig. 7). Besides in the vicinity of the KE axis, Miyama et al. (2021) reported that the summertime SST of the Oyashio region, 40-43°N and 143-147°E, abruptly jumped in 2010 and remained high until 2016, and the warm anomaly extended down to 200 m depth. The temperature became anomalously high even in winter from 2016 onward. Salinity also increased with temperature, and they concluded that the retreat of the Oyashio first intrusion and the increase in anticyclonic eddies from the KE caused a warm, saline anomaly.
Meanwhile, warming is evident over wide areas in the eastern NP and southwestern Bering Sea (boxes D and E in Figs. 3b and 8b). There is a clear decadal cycle in the OHC in box D, but the rise in the 2010s is unparalleled. The marine heatwave event called "Blob" started in the Gulf of Alaska in the winter of 2013/2014 (Bond et al. 2015), and the largescale warm anomaly expanded to the waters off the coast of California and the eastern Bering Sea by the autumn of 2014 (Di Lorenzo and Mantua 2016); the event ended in 2016. The near-surface warm anomaly in the Gulf of Alaska in the beginning of the event was caused by both reduced heat loss from the sea surface and horizontal heat advection (Bond et al. 2015), and propagated to depths greater than 300 m across the base of the mixed layer (Cummins and Masson 2018). The OHC change in Fig. 3b resembles the SST anomaly pattern of the Pacific Meridional Mode (PMM), except for the western part of the NP. Joh and Di Lorenzo (2019) indicated that the wind stress anomaly, which formed downstream of the KE in response to the KE stabilization, excited the PMM through weakening the trade wind, thereby modifying the Blob event and making it persistent (Di Lorenzo and Mantua 2016). The PMM can affect the OHC over a wide area, but we note that it does not always cause a drastic change in the OHC over the entire NP basin. The southward shift of the subarctic gyre boundary played a role in suppressing the OHC increase, but was overwhelmed by warming.

Meridional heat transport
The time series of the Kuroshio and Ekman transports, and the imbalance (V K + V E + V R + uncorrected V I ) are shown in Fig. 9. The transport across 137°E between 27°N and 34°N   Fig. 7 Latitude of the maximum horizontal gradient of temperature at 100 m depth (5-year running mean; black line). The latitude is averaged between 155°E and 150°W. Blue line indicates 5-year running means of the NPI. The NPI is defined as the sea level pressure anomaly in the region of 30-65°N, 160°E-140°W from November in the previous year to March normalized by standard deviation, and calculated by the JMA using JRA-55 (https:// www. data. jma. go. jp/ gmd/ kaiyou/ data/ db/ clima te/ pdo/ npiwin. txt) Fig. 8 Five-year running mean of the OHC in the boxes shown in Fig. 3. a A, B, and the rest of the target region. b C, D, E, and F exhibits variations on a time scale of 5-10 years with an average of 40.3 Sv, which is larger than the net Kuroshio transport in the Shikoku Basin south of Japan evaluated in previous studies, i.e., approximately 33.0 Sv (Qiu and Joyce 1992;Sugimoto et al. 2010;Nagano et al. 2013) (Fig. 9a). The transport between 24°N and 34°N was smaller, especially around 1990, owing to the westward flow in the south. We expediently refer to the 137°E transport up to 24°N as the Kuroshio transport, although it contains part of the North Equatorial Current, which should be included in the interior, V I . The Ekman transport, consisting of the northward flow across 24°N and the eastward flow across 137°E, increased in the 1980s and the late 1990s, as K08 indicated (Fig. 9b). The annual-mean imbalance ranged between −9.1 and 15.1 Sv, with an average of 2.1 Sv (Fig. 9c). The outflow across 24°N in the interior was forcibly strengthened uniformly as a correction to compensate for the excessive inflow into the target region.
The mean MHT throughout the period was 0.25 PW (Fig. 10a), and was smaller compared with previous studies (e.g. Roemmich et al. 2001;K08). MHT decreased in the late 1980s, which was earlier compared to the findings of K08. One of the reasons for this discrepancy was be the exclusion of the regions west of 137°E. MHT exhibits an increase of more than 0.1 PW in the 1980s and 1990s based on 5-year running means. Both the net Kuroshio and Ekman transports, defined as the first and second terms in Eq. 4, contributed to these increases (Fig. 10b).

Overall temporal change and assessment of heat flux
The net surface heat flux can be estimated by subtracting the OHC change rate (Fig. 2c) from the MHT (Eq. 1) and used as a benchmark. The net heat flux dropped by approximately 0.2 PW in the late 1980s and early 2010s, and gradually increased from around 1990 to 2010 with a slight decrease in the early 2000s ( Fig. 10a; negative values denote an increase in downward flux). Figure 11 depicts the four kinds Fig. 9 a Transport across 137°E. Blue and black lines indicate the transport integrated from 34°N to 27°N and to 24°N, respectively. b Ekman transport across the 24°N and 137°E lines. c Imbalance of the transport (V K + V E + V R + uncorrected V I ). Thin line in panel c shows annual means of 5-year-running-mean net surface heat flux integrated over the target region. These time series exhibit common temporal tendencies in the periods that we focus on: decreases in the late 1980s and early 2000s, and increase in the 1990s. The net flux of the three kinds of reanalysis also increased in the 2010s. J-OFURO3, however, did not indicate a decrease in the 2000s or an increase in the 2010s. In the late 2000s, both CFSR/CFSv2 and ERA5 indicated an increase in heat release from the ocean, while JRA-55 exhibited the opposite tendency. CFSR/CFSv2 exhibited decadal variations with the largest amplitudes among the datasets. Comparisons with the net heat flux derived from the heat budget indicated the following: (1) all the reanalyses were able to reproduce the decrease in the late 1980s, but they underestimated it; (2) ERA5 was the closest to the benchmark in the 1990s and 2000s; and (3) all the reanalyses failed to reproduce the decline in the 2010s, and only J-OFURO3 indicated a slight decrease. ERA5 data are used hereafter to examine the spatial variations.

Decrease in the late 1980s
The weakening of the net heat release occurred in the western NP (Fig. 12a). Differences in net heat flux mainly reflected that in turbulent heat flux, i.e., sensible plus latent heat fluxes (Fig. 12b). SST rose in the western and central NP during this period (Fig. 12e), but the increase in surface air temperature surpassed that, resulting in a decrease in the air-sea temperature difference (Fig. 12f); the surface wind speed weakened simultaneously (Fig. 12g). Decreases in wind speed and air-sea temperature differences were obvious in winter and corresponded to the weakening of the cold air outbreak from the continent. The change in shortwave radiation also contributed to the net heat flux decrease: warming by insolation was strengthened in the western subtropics (Fig. 12c). The contribution of longwave radiation was small (Fig. 12d).

Increase in the 1990s
Unlike in the 1980s, the SST rise in the western NP magnified the air-sea temperature difference (Fig. 13e, f) and increased the turbulent heat flux (Fig. 13b), mainly accounting for the change in net heat flux (Fig. 13a). The spatial pattern of the difference in downward shortwave Fig. 10 a Total heat transport across the 24°N and 137°E lines calculated by Eq. 2 (black), and net surface heat flux estimated by subtracting the OHC change rate from the total heat transport (gray). b Fiveyear running means of the net heat transport of the Kuroshio (bold black) and Ekman (gray dashed) transport defined as the first and second terms in Eq. 4;they are shown as deviations from the averages of the whole period   ERA5 from 1983ERA5 from (i.e., 1981ERA5 from -1985ERA5 from ) to 1989ERA5 from (i.e., 1987ERA5 from -1991. a Surface net heat flux, b turbulent (sensible plus latent) heat flux, c net shortwave radiation flux, d net longwave radiation flux, e SST, f SST-minus-air temperature difference, and g surface wind speed. Arrows in panel g indicate the difference of the wind vector radiation during this period was opposite to that in the 1980s (Fig. 13c), and the variation in longwave radiation was small (Fig. 13d).

Variation in the 2010s
The spatial pattern of the net heat flux reflected that of the turbulent heat flux (Fig. 14a, b), and the net longwave radiation slightly contributed to intensifying the heat release from the ocean (Fig. 14d). While the turbulent heat flux variation was dominant in the western and central NP, the heat release by longwave radiation was larger in the eastern NP.
Although the SST rise in the eastern NP was much higher than that in the western NP (Fig. 14e), the changes in air-sea temperature difference and wind speed were small or negative due to the southerly wind anomaly accompanied by the PMM (Di Lorenzo and Mantua 2016; Joh and Di Lorenzo 2019) (Fig. 14f, g), indicating that the atmosphere drove the ocean warming. Warming by shortwave radiation was noticeable off the North American coasts, which was almost offset by the change in net longwave radiation (Fig. 14c, d).
Although the OHC largely increased in the 2010s, the MHT into the target region decreased because the Kuroshio transport weakened. This means that ocean warming was ascribed to the weakening of surface heat loss. The decline in the surface heat flux estimated from the hydrographic observations (gray line in Fig. 11) may be too large, but at least the basin-wide increase shown by the three kinds of reanalysis dataset is unreasonable. The moderate decrease revealed by J-OFURO3 is reasonable.
While the turbulent heat flux increased from the end of the 2000s both in the reanalyses and J-OFURO3 (Fig. 15a), net radiation (i.e., longwave plus shortwave radiation) decreased in all datasets (Fig. 15b). The spatial distributions of J-OFURO3 in net and turbulent heat flux increments are similar to those of ERA5 (Fig. 16a, b), but J-OFURO3 exhibits greater decreases in shortwave and longwave radiation over the basin (Fig. 16c, d). The larger downward radiation decrease in J-OFURO3 might be more realistic. It is almost impossible to validate surface heat fluxes on a basin scale and to determine which flux component is responsible for the overestimation of the net surface heat loss. Although the OHC increase east of northern Japan was mainly caused by the advection of Kuroshio water (Miyama et al. 2021), the anomalous warming in the eastern NP resulted from large-scale atmospheric teleconnections (Di Lorenzo and Mantua 2016;Liang et al. 2017;Joh and Di Lorenzo 2019). Improved surface heat flux estimation is critical for understanding and reproducing the OHC change in the 2010s.

Discussions on the difference in heat distribution between the 1980s and 1990s
OHC leaped in the late 1980s, following the MHT increase. In contrast, the MHT increase was associated with heat release into the atmosphere in the late 1990s (Fig. 10a). What has yielded this contrast? The decrease in net heat flux was synchronous with the increase in the OHC change rate in the 1980s. Both the horizontal transport amplified by the spin-up of the subtropical gyre and the suppression of surface heat loss contributed to the anomalous OHC jump, although the OHC change rate lagged behind the net Kuroshio heat transport. One of the crucial factors is the strength of the westerlies and the Aleutian Low, which is represented by the NPI. The NPI increased from the mid-1980s to the early 1990s, indicating weakening of the wintertime Aleutian Low. A positive sea-level pressure anomaly appeared in the central NP during this period, resulting in the weakening of the cold air outbreak around Japan (Fig. 12g). The effects of the warmer air temperature and weaker surface wind dominated the SST increase in determining the turbulent heat flux, and the downward solar radiation also increased (Fig. 12b, c). They played an essential role in keeping the upper ocean warm. The trade wind around 24°N was enhanced and resulted in an increase in northward Ekman transport (Figs. 9b and 12g), which is consistent with K08. Meanwhile, the northward retreat of the subarctic gyre boundary in response to the weakening of the Aleutian Low (Fig. 7) paved the way for the OHC anomalies in the KE and KENB regions to propagate northeastward and spread over the NP (Figs.4 and  8a). The response of surface heat fluxes would be earlier than the warming of the KR and KE regions because the adjustment of the subtropical gyre to the wind stress curl requires a few years (e.g. Yasuda and Kitamura 2003;Qiu et al. 2014;Joh and Di Lorenzo 2019).
In contrast, the subarctic front moved southward in the late 1990s and prevented the warm anomaly from advecting north. The OHC anomalies in the subarctic and subtropics offset each other (Fig. 8a). The westerlies in the western NP became stronger in response to the NPI decrease ( Figs. 7 and 13g). The anomaly transported by the Kuroshio was contained within the subtropics and around Japan, and was released into the atmosphere via turbulent diffusion. Although the MHT increased in both the late 1980s and 1990s, the meridional shift in the subarctic gyre boundary was out of phase with the meridional transport in the subtropics, leading to a difference in heat distribution between these two decades.

Summary and conclusion
This study revisited the heat budget in the NP north of 24°N mainly using hydrographic observations to clarify why the heat distribution was different between the 1980s and 1990s, and extended the heat budget analysis to the 2010s. The target region was the upper 700 m, partitioned by the 24°N and 137°E lines. The OHC in the target region exhibited jumps around 1990 and in the 2010s, and it was nearly stable in between, unlike the nearly monotonous warming over the whole NP or globally (Fig. 2). The OHC jump around 1990 corresponded to the regime shift occurred in 1989/90, which was not associated with tropical variations unlike the other regime shifts (Yasunaka and Hanawa 2003). In contrast, the OHC surge in the 2010s was affected by the PMM from the tropics (Fig. 14;Di Lorenzo and Mantua 2016;Joh and Di Lorenzo 2019).
From the late 1980s to the early 1990s, the OHC increased in the coastal area southeast of Japan and the KE, and subsequently in the KENB; eventually, the anomaly propagated toward the Gulf of Alaska along the boundary of the subarctic gyre (Fig. 4). Meanwhile, the KR was also warmed (Fig. 3a,  c). The northward retreat of the subarctic gyre boundary in response to the weakening of the Aleutian Low (Fig. 7) paved the way for the OHC anomaly, which was caused by the spinup of the subtropical gyre (Yasuda and Kitamura 2003), to propagate northeastward and spread across the NP (Fig. 8a). On the contrary, the southward shift of the subarctic front suppressed the OHC rise, despite the MHT rise in the late 1990s (Figs. 7 and 8a). The MHT balanced the heat flux released to the atmosphere through surface turbulence instead of being stored in the ocean (Fig. 10a).
In the 2010s, unprecedented warming occurred in the eastern NP, southwestern Bering Sea, and east of northern Japan (Figs. 8a and 14). The PMM was mainly responsible for large-scale warming in the eastern part, and the advection of the Kuroshio water and retreat of the Oyashio resulted in warming in the western part (Miyama et al. 2021). The MHT estimation based on the hydrographic observations indicated that net surface heating must have been strengthened, but none of the latest atmospheric reanalysis data agreed with it (Fig. 11).
A study using a large ensemble of anthropogenic-forced model simulations indicated that the variability of the PMM would increase in the future (Liguori and Di Lorenzo 2018). Qiu et al. (2020) demonstrated that decadal-scale, windforced KE variability was disrupted in 2017 by the large Kuroshio meander. The MHT and heat distribution in the NP may change unexpectedly and irregularly due to anthropogenic forcing, the Kuroshio meander, or other factors, and it is indispensable to monitor the MHT and OHC, as well as the development of surface heat flux accuracy.