Extensive Marine Heatwaves at the Sea Surface in the Northwestern Paciﬁc Ocean in Summer 2021

: In July–August 2021, intense marine heatwaves (MHWs) occurred at the sea surface over extensive areas of the northwestern Paciﬁc Ocean, including the entire Sea of Japan and part of the Sea of Okhotsk. In extent and intensity, these MHWs were the largest since 1982, when satellite measurements of global sea surface temperatures started. The MHWs in summer 2021 were observed at the sea surface and occurred concomitantly with a stable shallow oceanic surface boundary layer. The distribution of the MHWs was strongly related to heat ﬂuxes at the sea surface, indicating that the MHWs were generated mainly by atmospheric forcing. The MHWs started to develop after around 10 July, concurrent with an extreme northward shift of the atmospheric westerly jet. The MHWs developed rapidly under an atmospheric high-pressure system near the sea surface, associated with a northwestward expansion of the North Paciﬁc Subtropical High. The MHWs exhibited peaks around 30 July to 1 August. Subsequently, following the southward displacement of the westerly jet, the MHWs weakened and then shrank abruptly, synchronously with rapid deepening of the oceanic surface boundary layer. By 18 August, the MHWs had disappeared.


Introduction
As generally defined, an extreme event occurs when the values of a weather or climate variable either exceed a threshold near the upper end of an observed range, or are lower than a threshold near the lower end [1]. Many extreme weather and climate events result from natural climate variability, but some extreme events have occurred as a result of anthropogenic climate change and ongoing global warming (e.g., [2][3][4]). Marine heatwaves (MHWs) are an example of an extreme climate in an oceanic system [5][6][7]. In the last century, from 1925 to 2016, the global average frequency and duration of MHWs increased by 34% and 17% respectively, resulting in a 54% increase in annual marine heatwave days [8]. According to Oliver [9], an increase in mean sea surface temperature (SST), rather than its variability, has been the dominant driver of the increasing frequency of MHW days over approximately two-thirds of the world's ocean, and it has been the dominant driver of changes in MHW intensity over approximately one-third of the ocean. Moreover, it has been projected that these historical trends in MHW properties will continue over the coming decades under global warming caused by anthropogenic climate change [10]. Furthermore, MHWs have devastatingly affected marine ecosystems in the past, and they are expected to do so in the future, through abrupt changes in biological habitat, mortality, reproduction, community structure, and so on [10][11][12][13][14][15].
Most published studies of MHWs in the North Pacific Ocean have investigated conditions in northeastern waters. A MHW known as "the blob," which occurred persistently in 2014-2016 in the northeastern Pacific, is especially famous (e.g., [16]). In contrast, a search

Materials and Methods
We analyzed gridded "Merged satellite and in situ data Global Daily Sea Surface Temperature" (MGDSST) datasets, comprising data from 1 January 1982 to 31 August 2021. These datasets, which were compiled by the Japan Meteorological Agency, have a spatial resolution of 0.25 • (latitude) × 0.25 • (longitude) and a daily temporal resolution [21].
We also used the "Roemmich-Gilson Argo Climatology," a global monthly dataset of seawater temperature in the subsurface created by the Scripps Institute of Oceanography [22] with a horizontal resolution of 1 • (latitude) × 1 • (longitude). We included data from January 2004 to July 2021 in our analysis.
We also used global atmospheric products at multiple vertical levels from the Japanese 55-year Reanalysis dataset (JRA-55), the second global atmospheric reanalysis project of the Japan Meteorological Agency [23]. These data have a horizontal resolution of 1.25 • (latitude) × 1.25 • (longitude). In our analysis, we used 3-hourly or 6-hourly outputs for the period from 1 January 1990 to 31 August 2021.
To infer the oceanic surface boundary layer conditions during the period of MHWs in 2021, numerical experiments were conducted with realistic ocean circulation models that are part of the operational ocean forecast system of the Japan Fisheries Research and Education Agency, "FRA-ROMS" [24]. The ocean circulation models, which are based on the Regional Ocean Modeling System (ROMS) [25][26][27], consist of a 1/2 • resolution parent model for the North Pacific and a 1/10 • resolution child model for the northwestern Pacific, connected by one-way nesting to simulate dominant basin-scale and mesoscale variations around Japan [28]. Daily mean reanalysis data for 1 July 2021, assimilated by using a three-dimensional variational scheme, were utilized for the initial conditions of 2-month hindcast and sensitivity experiments conducted without further data assimilation. To sequentially estimate momentum and heat fluxes at the sea surface during each model run, 31 sets of atmospheric elements near the sea surface (i.e., wind vectors at 10 m, air temperature and specific humidity at 2 m above the sea surface, air pressure, precipitation, and net shortwave and downward longwave radiation at the sea surface) were derived from the 3-hourly outputs of the JRA-55 system for July-August of 1990-2019 and 2021.
Climatological daily means were estimated for each element (e.g., SST, air temperature, sea level pressure), except subsurface ocean temperature, by averaging their data for 1990-2019. This 30-year averaging period, which was selected by referring to Hobday et al. [5], does not include the 1980s, during which a typical cold SST regime prevailed around Japan. Therefore, because the climatological averaging period excluded data from the 1980s, MHWs in recent years were less frequently detected than would have been the case had the selected climatological period included the 1980s. Further, because of limited data availability, climatological monthly means of the Argo-based subsurface Remote Sens. 2021, 13, 3989 3 of 13 ocean temperature were estimated by averaging data over just 15 years (2004-2018) instead of 30 years. Throughout this study, an anomaly was defined as a difference from the climatological mean.
The MHWs were detected from daily SST data and categorized according to intensity following Hobday et al. [6]. Intensity categories were defined based on the SST difference, Tdiff, between the 30-year climatological average (Tavg) and the 90th percentile value for each grid point. If the SST on a given day was ≥(Tavg + N × Tdiff ) and <(Tavg + (N + 1) × Tdiff ), it was categorized into intensity category N. In addition, the MHW index was defined as (SST -Tavg)/Tdiff. Note that in this paper, the duration of MHWs is referred to as their age.

Surface Ocean Conditions
3.1.1. Sea Surface Temperatures Figure 1a,b show intensity and age respectively, of MHWs on 30 July 2021, when the MHWs were at their maximum strength and extent, as explained later. The MHWs in the northwestern Pacific occurred within a zonal band, between 38 • N and 50 • N and extending from 120 • E to 170 • W, including the marginal seas, that is, the northern part of the Yellow Sea, the whole of the Sea of Japan, and the southern half of the Sea of Okhotsk. The MHW intensity on 30 July was classified into categories 1-4. Some areas of the MHWs corresponded to large SST anomalies > 6 • C (Figure 1d). The MHWs exhibited maximum intensity (i.e., category 4) on 30 July in the northern part of the Sea of Japan and tended to weaken gradually toward the east. The age of the MHWs in the northwestern Pacific on 30 July was within 1-30 days: it was longest (~30 days) in the northern part of the Sea of Japan and tended to be shorter further east. It is worth noting that intense MHWs appeared to start on the northwestern shelf of the Sea of Okhotsk around 10 July 2021 and subsequently expand from there (Figure 1c,d). However, the SST data analyzed in this study did not strictly discriminate between SSTs in the presence or absence of sea ice; thus, if sea ice remained on the northwestern shelf at the beginning of summer, the estimated MHWs would include large uncertainties. Therefore, whether the above feature is a real phenomenon is questionable.
We focused on MHWs east of Japan in the northwestern Pacific and excluded from our analysis the Sea of Japan, where MHWs that were more intensive than those in 2021 occurred locally after 1982 (e.g., January 2020). Thus, we estimated daily time series of SST anomalies and the MHW index averaged over the rectangular area within 143 It is worth noting that intense MHWs appeared to start on the northwestern shelf of the Sea of Okhotsk around 10 July 2021 and subsequently expand from there (Figure 1c,d). However, the SST data analyzed in this study did not strictly discriminate between SSTs in the presence or absence of sea ice; thus, if sea ice remained on the northwestern shelf at the beginning of summer, the estimated MHWs would include large uncertainties. Therefore, whether the above feature is a real phenomenon is questionable.
We focused on MHWs east of Japan in the northwestern Pacific and excluded from our analysis the Sea of Japan, where MHWs that were more intensive than those in 2021 occurred locally after 1982 (e.g., January 2020). Thus, we estimated daily time series of SST anomalies and the MHW index averaged over the rectangular area within 143°E-180° and 40°-50°N ( Figure 2). Both time series show historical maxima in July-August 2021 ( Figure  2a). The spatial mean of the daily MHW index exceeded 1.0 on a total of 207 days during 1982-2021, and 23 of these days were in July-August 2021. The spatial mean also exceeded 1.5 on 15 days during 1982-2021, of which 14 were in July-August 2021. These results indicate that the MHWs in summer 2021 were the most intense and extensive in the northwestern Pacific during the last four decades.
The transition to MHWs in summer 2021 can be interpreted as follows (Figure 2b). The MHWs started to develop on around 10 July 2021 (Figure 1c), when the spatial mean SST anomaly was close to zero ( Figure 2b). The MHWs reached maximum intensity and extent around 30 July to 1 August (Figure 1d), and then immediately started to weaken. By 18 August (Figure 1e), the spatial mean SST anomaly was again near zero ( Figure 2b) and the MHWs had almost entirely disappeared.  Figure  1a). As indices for evaluating the variability of the daily SST anomaly, the +3 standard deviation line (dashed blue line) is shown in (a), and the daily SST anomaly normalized by the standard deviation at each grid point and spatially averaged over the focal area (green line, corresponding to the right green axis) is shown in (b).

Heat Fluxes at the Sea Surface
Net heat flux anomalies at the sea surface during 10-30 July 2021 were mostly positive within the focal area, that is, the central area of MHWs ( Figure 3a); thus, the sea surface was heated by the atmosphere, relative to the climatology. The distribution of net heat flux anomalies (Figure 3a) also correlated with that of SST anomaly differences between 10 and 30 July (Figure 3b) (r = 0.69), suggesting that the MHWs were attributable primarily to increased net heat fluxes at the sea surface over the northwestern Pacific.
The net heat fluxes were decomposed into four components: net shortwave radiation, net longwave radiation, latent heat, and sensible heat ( Figure 4). Within our focal area, net shortwave radiation anomalies exhibited three positive maxima, centered around (145°E, 40°N), (170°E, 40°N), and (160°E, 45°N) ( Figure 4b). These three positive maxima corresponded spatially to three negative minima of the total cloud cover and low cloud cover anomalies (see Section 3.3.1). Further, the spatial patterns of positive and negative net  Figure 1a). As indices for evaluating the variability of the daily SST anomaly, the +3 standard deviation line (dashed blue line) is shown in (a), and the daily SST anomaly normalized by the standard deviation at each grid point and spatially averaged over the focal area (green line, corresponding to the right green axis) is shown in (b).

Heat Fluxes at the Sea Surface
Net heat flux anomalies at the sea surface during 10-30 July 2021 were mostly positive within the focal area, that is, the central area of MHWs ( Figure 3a); thus, the sea surface was heated by the atmosphere, relative to the climatology. The distribution of net heat flux anomalies ( Figure 3a) also correlated with that of SST anomaly differences between 10 and 30 July (Figure 3b) (r = 0.69), suggesting that the MHWs were attributable primarily to increased net heat fluxes at the sea surface over the northwestern Pacific.
The net heat fluxes were decomposed into four components: net shortwave radiation, net longwave radiation, latent heat, and sensible heat ( focal area (Figure 4d,h) corresponded to positive values of the respective components (Figure 4c,g). The positive values thus indicate that the sea surface was actually heated. Although MHWs with high SSTs were expected to cause large upward longwave radiation (i.e., negative heat fluxes), the net longwave radiation was nevertheless positive in some areas within the MHW area, for example, around the northwestern corner of the focal area (Figure 4c). This result suggests that strong downward longwave radiation from clouds exceeded the large upward longwave radiation from the sea surface. In fact, as described in Section 3.3.1, dense clouds tended to cover areas with positive net longwave radiation (Figure 4c). In any case, it should be emphasized that there were inter-regional differences in the origin of positive net heat flux anomalies over the area of MHWs; that is, the SST heating process differed inter-regionally within the area of MHWs.

Subsurface Ocean Conditions
In July 2021, temperatures in the subsurface at 20 m depth, just below the oceanic surface boundary layer, exhibited positive anomalies, mostly over the whole of the focal area (Figure 5c), although relatively large positive anomalies (>1 °C) were limited to its southern half. These relatively large positive anomalies were also distributed continu- Note that positive anomalies of net longwave radiation and sensible heat around the focal area (Figure 4d,h) corresponded to positive values of the respective components (Figure 4c,g). The positive values thus indicate that the sea surface was actually heated. Although MHWs with high SSTs were expected to cause large upward longwave radiation (i.e., negative heat fluxes), the net longwave radiation was nevertheless positive in some areas within the MHW area, for example, around the northwestern corner of the focal area ( Figure 4c). This result suggests that strong downward longwave radiation from clouds exceeded the large upward longwave radiation from the sea surface. In fact, as described in Section 3.3.1, dense clouds tended to cover areas with positive net longwave radiation (Figure 4c). In any case, it should be emphasized that there were inter-regional differences in the origin of positive net heat flux anomalies over the area of MHWs; that is, the SST heating process differed inter-regionally within the area of MHWs.

Oceanic Surface Boundary Layer
In this study, the depth of the oceanic surface boundary layer, defined as the shallowest depth where the bulk Richardson number is equal to its critical value, was derived from model simulations instead of the more conventional mixed layer depth (i.e., the depth at which the water density in the upper layer exceeds a threshold value, typically

Oceanic Surface Boundary Layer
In this study, the depth of the oceanic surface boundary layer, defined as the shallowest depth where the bulk Richardson number is equal to its critical value, was derived from model simulations instead of the more conventional mixed layer depth (i.e., the depth at which the water density in the upper layer exceeds a threshold value, typically 0.125 kg m −3 , relative to the 10 m density), primarily because the mixed layer was very thin during summer 2021 and could not be accurately evaluated by the conventional mixed layer depth.
Hindcast and sensitivity experiments conducted with realistic ocean circulation models revealed that thickness anomalies of the oceanic surface boundary layer during 10-30 July 2021 were mainly negative within the focal area, except for an area of positive anomalies that extended southwest-northeast, from (165 • E, 40 • N) to (180 • , 46 • N) (Figure 6a). The thickness of the oceanic surface boundary layer remained small in the focal area during the period from 13 July to 5 August 2016 (Figure 6b), during most of which MHWs were rapidly developing (Figure 2b). From 7 August, however, the thickness increased abruptly, at the same time that the MHWs were weakening rapidly (Figure 2b). These results suggest that positive net heat fluxes ( Figure 3a) and relatively weak wind stress (as explained in Section 3.3.1) near the sea surface stably maintained strong stratification, suppressed vertical convection, and accelerated the SST increase, that is, the rapid development of MHWs.

Conditions Near the Sea Surface
Sea level pressure anomalies during 10-30 July 2021 were widely positive in the northwestern Pacific, except in the northern Sea of Okhotsk and the western Bering Sea (Figure 7b). These positive anomalies were attributable to a northwestward expansion of the North Pacific Subtropical High (Figure 7a). The distribution of sea level pressure anomalies was also similar to MHW distribution (Figure 1a), although the distribution of air temperature anomalies near the sea surface (Figure 7c) more closely resembled the

Atmospheric conditions 3.3.1. Conditions Near the Sea Surface
Sea level pressure anomalies during 10-30 July 2021 were widely positive in the northwestern Pacific, except in the northern Sea of Okhotsk and the western Bering Sea (Figure 7b). These positive anomalies were attributable to a northwestward expansion of the North Pacific Subtropical High (Figure 7a). The distribution of sea level pressure anomalies was also similar to MHW distribution (Figure 1a), although the distribution of air temperature anomalies near the sea surface (Figure 7c) more closely resembled the MHW distribution (Figure 1a). Local maxima of the air temperature anomalies exceeded 3 • C.
As mentioned in Section 3.1.2, the distributions of positive and negative anomalies of total and low cloud cover (Figure 7e,f) were almost opposite to those of positive and negative net shortwave radiation anomalies (Figure 4b). In addition, in the focal area, the positive anomalies of total and low cloud cover (Figure 7e

Variation of the Westerly Jet at 200 hPa
Here, we focus on the westerly jet variation at 200 hPa as a potential driver or trigger of the atmospheric variations and MHWs in summer 2021. Maps of eastward wind velocity at 200 hPa averaged over 10-30 July (Figure 8a) revealed that the strong eastward wind (>20 m s −1 ) associated with the westerly jet [29] migrated northward, largely to the north of Japan, and skirted the focal area. The position of the velocity maximum reached 55 • N at around 140 • E. The velocity maximum between 120 • E and 180 • during 10-30 July 2021 was clearly located north of the climatological mean position (Figure 8b).
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 14 July, the meridional position of the jet in 143°E-180° stayed north of 50°N until the beginning of August, although it was intermittently disturbed (e.g., by typhoon events in late July): 50°N corresponded to the northern boundary of the MHWs in the northwestern Pacific (Figure 1a). After around 1-8 August, the axial position of the jet was displaced to the south. Simultaneously, the MHWs in the northwestern Pacific started to weaken and shrink.
The large northward shift of the westerly jet identified in early to mid-July exceeded the climatological mean plus 3 standard deviations ( Figure 9b); therefore, it was interpreted as a rare event. This result suggests that the abnormally large northward shift of the westerly jet triggered warm atmospheric conditions near the sea surface and thus caused the strongest, most extensive MHWs yet recorded over the northwestern Pacific.   July, the meridional position of the jet in 143°E-180° stayed north of 50°N until the beginning of August, although it was intermittently disturbed (e.g., by typhoon events in late July): 50°N corresponded to the northern boundary of the MHWs in the northwestern Pacific (Figure 1a). After around 1-8 August, the axial position of the jet was displaced to the south. Simultaneously, the MHWs in the northwestern Pacific started to weaken and shrink.
The large northward shift of the westerly jet identified in early to mid-July exceeded the climatological mean plus 3 standard deviations ( Figure 9b); therefore, it was interpreted as a rare event. This result suggests that the abnormally large northward shift of the westerly jet triggered warm atmospheric conditions near the sea surface and thus caused the strongest, most extensive MHWs yet recorded over the northwestern Pacific.   The northward shift of the westerly jet began suddenly at the beginning of July: large northward shifts of the westerly jet to the north of 50 • N started at 70 • E on 5 July, propagated eastward, and reached 120 • -140 • E around Japan on 11-15 July (Figure 9a). After 15 July, the meridional position of the jet in 143 • E-180 • stayed north of 50 • N until the beginning of August, although it was intermittently disturbed (e.g., by typhoon events in late July): 50 • N corresponded to the northern boundary of the MHWs in the northwestern Pacific (Figure 1a). After around 1-8 August, the axial position of the jet was displaced to the south. Simultaneously, the MHWs in the northwestern Pacific started to weaken and shrink.
The large northward shift of the westerly jet identified in early to mid-July exceeded the climatological mean plus 3 standard deviations ( Figure 9b); therefore, it was interpreted as a rare event. This result suggests that the abnormally large northward shift of the westerly jet triggered warm atmospheric conditions near the sea surface and thus caused the strongest, most extensive MHWs yet recorded over the northwestern Pacific.

Discussion
In this article, we reported on the most extensive and intense MHWs in the historical record of the last 40 years, which occurred in the northwestern Pacific during July-August 2021. However, we mainly limited our explanation to certain features: SST, the oceanic surface boundary layer, heat flux at the sea surface, and atmospheric conditions. We only briefly mentioned the redistribution of heat from the sea surface to the ocean subsurface by ocean advection and diffusion. In addition, although differences in SST anomalies were positively correlated with net heat flux anomalies (r = 0.69, Figure 3), anomalies of SST and net heat flux were not strongly correlated with each other. In considering these problems, the uncertainties of the SST and heat flux data should first of all be meticulously validated by comparing different datasets with each other and with in situ data (e.g., surface buoy and vessel data). Moreover, it should be noted that the core area of the MHWs was at 43 • -46 • N (not shown), whereas the net heat flux anomalies did not show corresponding maxima at 43 • -46 • N. Additionally, as described in Section 3.1.2, SST heating processes differed inter-regionally within the area of MHWs. Hence, we anticipate that a feedback process from the ocean to the atmosphere contributed to localized enhancement of MHWs. Air-sea processes that might be involved in the enhancement of MHWs at mid-latitudes include SST-cloud feedback (e.g., [30][31][32]), surface heat flux feedback dominated by turbulent flux (e.g., [33]), or an air-sea process related to the imbalance between incoming and outgoing fluxes through the sea surface (e.g., [34]). A useful first step to clarify such processes would be to conduct numerical experiments with a coupled air-sea model.
In this article, we also focused mainly on the development of MHWs in the northwestern Pacific, but the MHWs also decayed rapidly after the westerly jet was displaced southward at the beginning of August 2021. Therefore, we provide here a brief description of atmospheric conditions near the sea surface during 1-18 August 2021. In this period, the North Pacific Subtropical High retreated southward, and the Okhotsk High strengthened (Figure 10a). The sea level pressure valley between the two high-pressure systems was located in the southeastern half of the focal area (143 • E-180 • , 40 • -50 • N) (Figure 10a,b). Air temperature greatly decreased over the northwestern Pacific after 1 August (Figure 10c), and the negative net heat flux anomalies at the sea surface in the focal area (Figure 10d) contributed to sea surface cooling. Simultaneously, as shown by hindcast experiments, the oceanic surface boundary layer deepened (Figure 6b), and this deepening implies that entrainment of subsurface water into the surface boundary layer accelerated cooling at the sea surface. As a consequence, the extensive and intense MHWs had disappeared by 18 August.

Conclusions
In July-August 2021, the largest and most intense MHWs of the satellite era occurred at the sea surface over extensive areas of the northwestern Pacific Ocean. These were observed at the sea surface and were accompanied by a stable shallow oceanic surface boundary layer. Their spatial relationship to heat fluxes at the sea surface indicates that the MHWs were generated mainly by atmospheric forcing associated with northwestward expansion of the North Pacific Subtropical High and northward displacement of the westerly jet. Additional air-sea studies are needed to understand inter-regional differences in SST heating processes within the area of MHWs.
Moreover, the MHWs in summer 2021 were limited to the vicinity of the sea surface (<20 m depth), where primary production is especially high (e.g., [35][36][37]). In fact, some studies have reported impacts of MHWs on lower trophic levels of marine ecosystems (e.g., [38,39]). Commercially important fisheries' resources include Japanese sardine Sardinops melanostictus and Pacific saury Cololabis saira, which utilize the northwestern Pacific waters around the focal area as a nursery ground in summer (e.g., [40,41]). Hence, for sustainable management of ecosystems in the northwestern Pacific Ocean, physical and biological studies are essential, not only to evaluate influences of the 2021 MHWs but also to project the influences of future MHWs on marine ecosystems.  Agency of Japan. This work was also partially supported by the Fisheries Agency (promotion project to precisely estimate fish stock size around Japan (no grant number)), but the contents of this study do not necessarily reflect the views of the Fisheries Agency.
Data Availability Statement: All the data analyzed in this study are publicly available. SST data were collected from a data server of NEAR-GOOS (https://ds.data.jma.go.jp) (accessed on 2 September 2021). JRA-55 data were downloaded from the University Corporation for Atmospheric Research data server (https://rda.ucar.edu) (accessed on 4 September 2021). Gridded Argo-based seawater temperatures in the subsurface were obtained from Scripps Institute of Oceanography

Conclusions
In July-August 2021, the largest and most intense MHWs of the satellite era occurred at the sea surface over extensive areas of the northwestern Pacific Ocean. These were observed at the sea surface and were accompanied by a stable shallow oceanic surface boundary layer. Their spatial relationship to heat fluxes at the sea surface indicates that the MHWs were generated mainly by atmospheric forcing associated with northwestward expansion of the North Pacific Subtropical High and northward displacement of the westerly jet. Additional air-sea studies are needed to understand inter-regional differences in SST heating processes within the area of MHWs.
Moreover, the MHWs in summer 2021 were limited to the vicinity of the sea surface (<20 m depth), where primary production is especially high (e.g., [35][36][37]). In fact, some studies have reported impacts of MHWs on lower trophic levels of marine ecosystems (e.g., [38,39]). Commercially important fisheries' resources include Japanese sardine Sardinops melanostictus and Pacific saury Cololabis saira, which utilize the northwestern Pacific waters around the focal area as a nursery ground in summer (e.g., [40,41]). Hence, for sustainable management of ecosystems in the northwestern Pacific Ocean, physical and biological studies are essential, not only to evaluate influences of the 2021 MHWs but also to project the influences of future MHWs on marine ecosystems. Funding: This research was supported by KAKENHI, Grant Number 19K06216, from the Ministry of Education, Culture, Sports, Science and Technology, and by the Environment Research and Technology Development Fund (4-2102) of the Environmental Restoration and Conservation Agency of Japan. This work was also partially supported by the Fisheries Agency (promotion project to precisely estimate fish stock size around Japan (no grant number)), but the contents of this study do not necessarily reflect the views of the Fisheries Agency.
Data Availability Statement: All the data analyzed in this study are publicly available. SST data were collected from a data server of NEAR-GOOS (https://ds.data.jma.go.jp) (accessed on 2 September 2021). JRA-55 data were downloaded from the University Corporation for Atmospheric Research data server (https://rda.ucar.edu) (accessed on 4 September 2021). Gridded Argo-based seawater temperatures in the subsurface were obtained from Scripps Institute of Oceanography (http://sioargo.ucsd.edu/RG_Climatology.html) (accessed on 17 August 2021). Daily reanalysis data from FRA-ROMS are available at http://fm.dc.affrc.go.jp/fra-roms after user registration (accessed on 15 August 2021). The source code of the ocean circulation models based on Rutgers University's ROMS is available at https://www.myroms.org (accessed on 15 August 2021).