Unprecedented Outbreak of Harmful Algae in Paciﬁc Coastal Waters off Southeast Hokkaido, Japan, during Late Summer 2021 after Record-Breaking Marine Heatwaves

+81-154-92-1723 Abstract: Unprecedented large-scale harmful algae blooms (HABs) were reported in coastal waters off the south-eastern coast of Hokkaido, Japan, in mid-to-late September 2021, about a month after very intense and extensive marine heatwaves subsided. To understand the physical–biological processes associated with development of the HABs, we conducted analyses via a combination of realistic ocean circulation models, particle-tracking simulations, and satellite measurements. The satellite-derived chlorophyll concentrations (SCCs) and areal extent of the high SCCs associated with the HABs were the highest recorded since 1998. More speciﬁcally, the extent of SCCs exceeding 5 or 10 mg m − 3 started to slowly increase after 20 August, when the marine heatwaves subsided, intermittently exceeded the climatological daily maximum after late August, and reached record-breaking extremes in mid-to-late September. About 70% of the SCCs that exceeded 10 mg m − 3 occurred in places where water depths were <300 m, i.e., coastal shelf waters. The high SCCs were also tightly linked with low-salinity water (e.g., subarctic Oyashio and river-inﬂuenced waters). High-salinity subtropical water (e.g., Soya Warm Current water) appeared to suppress the occurrence of HABs. The expansion of the area of high SCCs seemed to be synchronized with the deepening of surface mixed layer depths in subarctic waters on the Paciﬁc shelves. That deepening began around 10 August, when the marine heatwaves weakened abruptly. However, another mechanism was needed to explain the intensiﬁcation of the SCCs in very nearshore waters off southeast Hokkaido. Particle-tracking simulations based on ocean circulation models identiﬁed three potential source areas of the HABs: the Paciﬁc Ocean east of the Kamchatka Peninsula, the Sea of Japan, and the Sea of Okhotsk east of the Sakhalin Island. Different processes of HAB development were proposed because distance, time, and probability for transport of harmful algae from the potential source areas to the study region differed greatly between the three source areas.


Introduction
Harmful algal blooms (HABs) constitute a critical worldwide problem and have devastating impacts on marine ecosystems, aquaculture, fisheries, and local economies, as well as on human health and wellbeing [1][2][3][4][5][6]. HABs are natural phenomena, but they can be stimulated by the impacts of anthropogenic activities on coastal ecosystems, primarily nutrient over-enrichment, e.g., [4,7]. In addition, local oceanographic conditions tightly linked with basin-scale climate changes have sometimes triggered HABs [8]. Examples include the extreme blooms of Pseudochattonella cf. verruculosa that occurred in the coastal waters of western Patagonia during the summer of 2016 due to the strong El Niño and the positive phase of the Southern Annular Mode [9], and massive Pseudo-nitzschia blooms 2 of 21 that occurred persistently along the entire United States West Coast during 2013-2016 in response to anomalously warm seawater referred to as "blob" [10]. HABs can be generated not only by long-term secular trends or decadal variability of oceanographic variables (e.g., sea surface temperatures) but also by weather and extreme marine conditions with short timescales (e.g., hurricanes or intra-seasonal marine heatwaves) [11][12][13].
In mid-September 2021, about a month after the subsidence of the most intense and extensive marine heatwaves over the north-western Pacific Ocean [14], an unprecedented large-scale outbreak of harmful algae occurred in Pacific coastal waters off the south-eastern coast of Hokkaido, Japan. There had been no previous reports of catastrophic disasters attributed to HABs over these coastal waters, although there were a few reports of occurrence of harmful algae in a limited area of this region during the 1970-1980s [15]. According to public disclosures from the Hokkaido Research Organization, the HABs were composed of at least four kinds of harmful algae; two species, Karenia selliformis and K. mikimotoi, had been identified by 7 October 2021, and the rest of them have not been taxonomically specified. In general, blooms of K. mikimotoi have been frequently reported from Japanese and Chinese coastal waters south of 40 • N (e.g., [6]), whereas K. selliformis [16,17] was massively identified in Japanese waters for the first time. Devastating damage to marine ecosystems by K. selliformis was reported to have occurred to the east of the Kamchatka Peninsula during September-October 2020, just a year before the HABs off southeast Hokkaido [18,19]. Likewise, along with outbreaks of harmful algae, serious damages to coastal ecosystems and fisheries were reported off southeast Hokkaido (green bold lines in Figure 1b): sea urchins experienced mass die-offs, chum salmon died in fixed nets, and fish juveniles died in rearing facilities. The relationships between the biological damage and HABs are being carefully investigated. Nevertheless, there is little understanding of the spatiotemporal development of HABs.
J. Mar. Sci. Eng. 2021, 9,1335 3 of 21 and the processes by which the HABs developed. Our methodology was based on analysis combined with state-of-the-art scientific techniques that allowed us to immediately analyse the overall features of phenomena. A short time required for analysis can be very useful, particularly at the initial stage of an emergency, because collecting, processing, and analysing in situ data generally takes a long time and does not capture all of the spatiotemporal variability of a phenomenon.

Satellite-Derived Data
Datasets of satellite-derived chlorophyll concentrations at the sea surface (hereafter There have been few reports of serious coastal pollution within the study region. The environment has been maintained in a natural, healthy condition along the Pacific coast of southeast Hokkaido. Extensive terrestrial and industrial developments have not occurred, partly because there are some terrestrial sanctuaries such as the Kushiro Marsh National Park and the Akkeshi-Kiritappu-Konbumori quasi-National Park, parts of which are registered as Ramsar Sites in Japan. Oceanographic conditions in the study area are influenced by the Oyashio, the western boundary current of the Western Subarctic Gyre in the North Pacific Ocean, which is the downstream extension of the East Kamchatka Current, which flows along the east coast of the Kamchatka Peninsula and the northern Kuril Islands (Figure 1) [20,21]. The Oyashio supplies nutrient-rich, cold, low-salinity subarctic water to the study region throughout the year [22]. Meanwhile, subtropical, oligotrophic waters are transported from the subtropics to the study region, mainly through two pathways ( Figure 1). Soya Warm Current water is supplied seasonally from the Sea of Okhotsk during summer-autumn by a downstream extension of the Soya Warm Current [23][24][25]. This water originates in the Kuroshio, which bifurcates south of Kyushu, intrudes into the Sea of Japan, flows northward as the Tsushima Warm Current, and reaches the Sea of Okhotsk. Part of this warm water outflows into the North Pacific through several straits around the Four Islands, to the east of Hokkaido. The other pathway involves subtropical Kuroshio waters that are transported to the study region in some years by clockwise mesoscale eddies referred to as Kuroshio Warm Core Rings. These eddies detach from the Kuroshio Extension, move to the north, and stagnate around the Kuril-Kamchatka Trench, off the southeast coast of Hokkaido [26][27][28]. The HABs in the late summer of 2021 occurred in the vicinity of the confluence of the subarctic and subtropical waters.
The purpose of this study was to use near-real-time satellite-derived ocean color data to describe the features of the HABs in the late summer of 2021 and to relate them to oceanic environmental conditions simulated by a realistic high-resolution (1/50 • ) ocean model. Particle-tracking simulations were conducted to infer the potential source areas and the processes by which the HABs developed. Our methodology was based on analysis combined with state-of-the-art scientific techniques that allowed us to immediately analyse the overall features of phenomena. A short time required for analysis can be very useful, particularly at the initial stage of an emergency, because collecting, processing, and analysing in situ data generally takes a long time and does not capture all of the spatiotemporal variability of a phenomenon.

Satellite-Derived Data
Datasets of satellite-derived chlorophyll concentrations at the sea surface (hereafter "SCCs") were collected from the European Union's Copernicus Marine Service. We analyzed three types of global daily datasets with a horizontal grid of about 4 km: nearreal-time L3 and L4 products for July-October in 2021, and reprocessed L4 products for July-October in 1998-2020 (the product identifiers were OCEANCOLOUR_GLO_CHL_L3_ NRT_OBSERVATIONS_009_032, OCEANCOLOUR_GLO_CHL_L4_NRT_OBSERVATIONS _009_033, and OCEANCOLOUR_GLO_CHL_L4_REP_OBSERVATIONS_009_082, respectively). These products were based on the merging of sensors SeaWiFS, MODIS, MERIS, VIIRS-SNPP, and JPSS1. The L3 product was created using single-day measurements, and the L4 products were created by a spatiotemporal interpolation for data missing due to cloud coverage. Reprocessed L4 products were employed in this study mainly to estimate climatological baselines (e.g., daily means, maximum and minimum values, and standard deviations) and to compare SCCs between 2021 and 1998-2020.

Reanalysis Data from a 1/10 • Ocean Forecast System
We used daily mean reanalysis data from an operational ocean forecast system developed by the Japan Fisheries Research and Education Agency (FRA) based on the Regional Ocean Modeling System (FRA-ROMS) [29]. In this system, 1/2 • and 1/10 • ocean models based on the ROMS (e.g., [30][31][32]) were connected by one-way nesting and coupled with three-dimensional variational objective analysis schemes in which sea level anomalies, SSTs, and profiles of in situ temperature and salinity (TS) data were assimilated at weekly time intervals. Figure 2 shows these model domains. This system successfully simulated realistic basin-scale and mesoscale variations in the northwestern Pacific Ocean [29]. In this study, the reanalysis data from the 1/10 • ocean model were used to force a realistic 1/50 • ocean model and particle-tracking simulations.

Particle-Tracking Simulation
To infer the origin of the HABs along the Pacific coast off southeast line particle-tracking simulations were carried out using daily mean ou models. Because the domain of the 1/50° model was limited to the vicinity o 2), daily mean reanalysis data from the 1/10° FRA-ROMS were combined w outputs from the 1/50° model to perform particle-tracking simulations ov (i.e., green dashed box in Figure 2). For the particle-tracking simulations, w val TRANSport (LTRANS) model [45], which is configured for ROMS o source codes were modified to suit our model's configuration and purpose Particle-tracking simulations were conducted via the following two st step, 95,900 particles were initially set at each of three depths (1, 5, and 10 sea surface in Pacific coastal waters off southeast Hokkaido (i.e., occu HABs). Initial conditions were updated at intervals of five days from 10 to 2021. Particles were passively transported at a fixed depth backward in tim uary 2021 by reversing the sign of the horizontal model velocity. A small fusion constant (5.0 m 2 s -1 ) was assumed. In the second step, 87,616 particle set at each of three depths (1, 5, and 10 m) below the sea surface east of t -resolving ocean model, respectively. For particle-tracking simulations, outputs from the 1/10 • and 1/50 • ocean models were combined within the green dashed box.

High-Resolution 1/50 • Ocean Model
Most of the configurations of the 1/50 • ocean model were identical to those used by Kuroda et al. [25,33], except for the model domain, which was extended to include the whole of the Japanese main island (light blue square in Figure 2) [34]. As shown by Kuroda et al. [25,33], the overall structures of oceanographic conditions-particularly temperature, salinity, and current velocity near the sea surface around Pacific coastal shelf waters off the south-eastern coast of Hokkaido-were reasonably reproduced.
Inputs of heat and momentum fluxes at the sea surface were based on 3-hourly mean net shortwave and downward longwave radiation data from the Japan Meteorological Agency's 55-year reanalysis dataset (JRA-55) [35]. Inputs of hourly meteorological data (wind speed and direction, relative humidity, air temperature, precipitation, and atmospheric pressure) were taken from the Grid Point Value-Mesoscale Model (GPV/MSM) [36]. Daily mean discharges of water from all freshwater systems were estimated and introduced to the 1/50 • model as described by Kuroda et al. [37]. Lateral boundary conditions were derived from the daily mean reanalysis data of the 1/10 • FRA-ROMS. Barotropic tidal currents and elevations for the eight major tidal constituents (K 1 , O 1 , P 1 , Q 1 , M 2 , S 2 , K 2 , and N 2 ) were derived from the NAO.99Jb regional tidal model [38] and added to the lateral boundary values derived from the 1/10 • FRA-ROMS. Moreover, to modify mesoscale variations in the 1/50 • model, we introduced spectral nudging [39,40], which has been adapted to many regional climate models (e.g., [41][42][43]). TS data from the 1/50 • model and the 1/10 • FRA-ROMS were low-pass-filtered by the same spatial filter to extract TS variations associated with mesoscale variability at a spatial scale greater than about 100 km. The TS differences between the two models were fed into the 1/50 • model in 3-day increments by incremental analysis updating [44]. As a result, values of TS within about 50 km from the coast (i.e., on the shelf) and submesoscale variations at spatial scales of less than about 100 km were not modified.
We analyzed outputs of daily means from the 1/50 • model to characterize oceanographic conditions and conduct particle-tracking simulations. Surface mixed layer depth was estimated as the depth at which the density in the upper layer exceeded a threshold value of 0.125 kg m −3 relative to the density at the surface.

Particle-Tracking Simulation
To infer the origin of the HABs along the Pacific coast off southeast Hokkaido, off-line particle-tracking simulations were carried out using daily mean outputs of ocean models. Because the domain of the 1/50 • model was limited to the vicinity of Japan (Figure 2), daily mean reanalysis data from the 1/10 • FRA-ROMS were combined with daily mean outputs from the 1/50 • model to perform particle-tracking simulations over a wide area (i.e., green dashed box in Figure 2). For the particle-tracking simulations, we used the Larval TRANSport (LTRANS) model [45], which is configured for ROMS output, but the source codes were modified to suit our model's configuration and purpose [25,46,47].
Particle-tracking simulations were conducted via the following two steps. In the first step, 95,900 particles were initially set at each of three depths (1, 5, and 10 m) below the sea surface in Pacific coastal waters off southeast Hokkaido (i.e., occurrence area of HABs). Initial conditions were updated at intervals of five days from 10 to 30 September 2021. Particles were passively transported at a fixed depth backward in time up to 1 January 2021 by reversing the sign of the horizontal model velocity. A small horizontal diffusion constant (5.0 m 2 s -1 ) was assumed. In the second step, 87,616 particles were initially set at each of three depths (1, 5, and 10 m) below the sea surface east of the Kamchatka Peninsula and tracked forward in time for 360 days. Initial conditions of the particles were updated at intervals of five days from 30 September to 20 October 2020, which corresponded to the period of occurrence of the HABs off the east coast of the Kamchatka Peninsula [18]. Figure 3 shows maps of daily SCCs, which were associated with the process of HAB development. In these maps, SCCs exceeding the mean plus three standard deviations at individual grid cells (i.e., extremely high SCCs) are emphasized by black polygons. At the end of August (Figure 3a), relatively high SCCs started to appear around the study region but hardly exceeded 10 mg m −3 . At the beginning of September (Figure 3b), SCCs increased abruptly south of Etorofu Island and exceeded 10 mg m −3 near the centre of patch-like local maxima. On 7 September (Figure 3c), the area of the patch-like SCCs increased, and the patches seemed to expand south-westward to the south of Shikotan Island. On 14 September (Figure 3d), the patch-like local maximum became less distinct, and the patch elongated in a clockwise direction on the continental slope, probably because of the combined effects of horizontal advection and diffusion. At the same time, high SCCs (>5 mg m −3 ) suddenly appeared along the entirety of the Pacific coast off southeast Hokkaido. On 20 September (Figure 3e), the SCCs increased dramatically and exceeded 10 mg m −3 in coastal shelf waters along the Hokkaido coast and in slope waters. On 27 September (Figure 3f), the area of high SCCs (>10 mg m −3 ) in coastal shelf waters expanded slightly in an offshore direction, and the high SCCs in slope waters decreased and were mostly lower than 10 mg m −3 . These results implied that occurrence of HABs was localized around the Four Islands and along the southeast coast of Hokkaido. The HABs were, therefore, isolated from the surrounding waters off the Kuril Islands and in the Sea of Okhotsk. 20 September (Figure 3e), the SCCs increased dramatically and exceeded 10 mg m −3 in coastal shelf waters along the Hokkaido coast and in slope waters. On 27 September (Figure 3f), the area of high SCCs (>10 mg m −3 ) in coastal shelf waters expanded slightly in an offshore direction, and the high SCCs in slope waters decreased and were mostly lower than 10 mg m −3 . These results implied that occurrence of HABs was localized around the Four Islands and along the southeast coast of Hokkaido. The HABs were, therefore, isolated from the surrounding waters off the Kuril Islands and in the Sea of Okhotsk. To examine the areal extent of high SCCs in 2021( Figure 4)., we integrated the area of a grid cell with SCCs exceeding 5 or 10 mg m −3 over the pink rectangular region (hereafter "focal region") in Figure 3 In situ chlorophyll a concentrations >5 mg m −3 have been used as typical indicators of massive spring blooms in the Oyashio area [48]. The area of high SCCs started to increase slowly after the middle of August 2021. After late August, the area intermittently became comparable to the climatological daily maximum over 23 years (1998-2020), i.e., the upper boundary of red shading in Figure 4. The area then increased dramatically in mid-to late September 2021. At that time, the area was sometimes more than twice as large as the 23-year daily maximum area. The 23-year daily maximum area was also comparable to the mean plus 3 standard deviations (blue line in Figure 4). The implication was that the area of high SCCs in mid-to-late September sometimes reached the mean plus six standard deviations. The large area of the high SCCs in mid-tolate September 2021 could, therefore, be interpreted as a very rare, extreme phenomenon.  20 September, and (f) 27 September 2021, when there were few clouds (light-grey shading) over the study region. Inside of black polygons, the daily SCC exceeded the monthly mean plus 3 standard deviations for September 1998-2020 at each grid point. The pink dotted rectangle is referred to as the "focal region" throughout this study. SCCs were based on near-real-time L3 and reprocessed L4 products.

Satellite-Derived Chlorophyll Concentration
To examine the areal extent of high SCCs in 2021( Figure 4), we integrated the area of a grid cell with SCCs exceeding 5 or 10 mg m −3 over the pink rectangular region (hereafter "focal region") in Figure 3 In situ chlorophyll a concentrations > 5 mg m −3 have been used as typical indicators of massive spring blooms in the Oyashio area [48]. The area of high SCCs started to increase slowly after the middle of August 2021. After late August, the area intermittently became comparable to the climatological daily maximum over 23 years (1998-2020), i.e., the upper boundary of red shading in Figure 4. The area then increased dramatically in mid-to late September 2021. At that time, the area was sometimes more than twice as large as the 23-year daily maximum area. The 23-year daily maximum area was also comparable to the mean plus 3 standard deviations (blue line in Figure 4). The implication was that the area of high SCCs in mid-to-late September sometimes reached the mean plus six standard deviations. The large area of the high SCCs in mid-to-late September 2021 could, therefore, be interpreted as a very rare, extreme phenomenon.
Next, to examine the duration of the high SCCs in 2021, we counted the number of days during September 2021 when the SCCs exceeded 5 or 10 mg m −3 at each grid ( Figure 5). Long durations were apparent around the Pacific coastal, shelf, and slope regions off the Four Islands and southeast of Hokkaido. The duration tended to be longer near the coast and shorter offshore, except for the elongated local maximum that extended southwestward from Etorofu Island. In coastal waters off southeast Hokkaido, the durations exceeded 20 (10) days for the threshold of 5 (10) mg m −3 . For the region illustrated in Figure 5, 74% of the areas that exceeded the threshold SCC of 5 mg m −3 and 76% of the areas that exceeded the threshold SCC of 10 mg m −3 exhibited the longest duration of exceedance of the respective thresholds during September in 1998-2021. The implication was that the persistence of the high SCCs measured around the focal region broke historical records in September 2021.
Next, to examine the duration of the high SCCs in 2021, we counted the number of days during September 2021 when the SCCs exceeded 5 or 10 mg m −3 at each grid ( Figure  5). Long durations were apparent around the Pacific coastal, shelf, and slope regions off the Four Islands and southeast of Hokkaido. The duration tended to be longer near the coast and shorter offshore, except for the elongated local maximum that extended southwestward from Etorofu Island. In coastal waters off southeast Hokkaido, the durations exceeded 20 (10) days for the threshold of 5 (10) mg m −3 . For the region illustrated in Figure  5, 74% of the areas that exceeded the threshold SCC of 5 mg m −3 and 76% of the areas that exceeded the threshold SCC of 10 mg m −3 exhibited the longest duration of exceedance of the respective thresholds during September in 1998-2021. The implication was that the persistence of the high SCCs measured around the focal region broke historical records in September 2021.    Figure 3. White circles denote climatological daily means over the 23 years (1998-2020), and the blue line indicates the climatological daily mean plus 3 standard deviations. Red shading shows the range between the daily minimum and maximum over the 23 years. SCCs were calculated from near-real-time and reprocessed L4 products.
Next, to examine the duration of the high SCCs in 2021, we counted the number of days during September 2021 when the SCCs exceeded 5 or 10 mg m −3 at each grid ( Figure  5). Long durations were apparent around the Pacific coastal, shelf, and slope regions off the Four Islands and southeast of Hokkaido. The duration tended to be longer near the coast and shorter offshore, except for the elongated local maximum that extended southwestward from Etorofu Island. In coastal waters off southeast Hokkaido, the durations exceeded 20 (10) days for the threshold of 5 (10) mg m −3 . For the region illustrated in Figure  5, 74% of the areas that exceeded the threshold SCC of 5 mg m −3 and 76% of the areas that exceeded the threshold SCC of 10 mg m −3 exhibited the longest duration of exceedance of the respective thresholds during September in 1998-2021. The implication was that the persistence of the high SCCs measured around the focal region broke historical records in September 2021.   indicates where daily SCCs were lower than the threshold value throughout the month. SCCs were calculated from the near-real-time L4 product.

Relationships to Environmental Conditions
To understand relationships between high SCCs (i.e., potential occurrence of HABs) and oceanographic conditions, we analyzed outputs from the 1/50 • ocean model. Monthly mean salinities and current velocities at the sea surface ( Figure 6a) reflected representative ocean current conditions around Hokkaido in September 2021. The Oyashio flowed from the northeast to the southwest along the continental slope. A part of the Oyashio bifurcated to the south of Shikotan Island, was deflected in a clockwise direction, and then rejoined the mainstream of the Oyashio. Meanwhile, the Soya Warm Current, which is characterized by high-salinity water (>33.5), flowed along the Okhotsk coast of northeast Hokkaido and intruded into the Kunashiri Strait. Water with relatively high salinity (~33.5) that outflowed into the Pacific shelf was then transported in a south-westward direction by the Coastal Oyashio onto the Pacific shelf, inshore of the Oyashio on the slope. The salinity along the Coastal Oyashio decreased gradually downstream from the Kunashiri Strait. The implication was that Soya Warm Current water was modified by mixing with the surrounding low-salinity waters, and so-called modified Soya Warm Current water was formed on the Pacific shelf. On the inshore side of the Coastal Oyashio, water with very low salinity (<32.5) due to river discharge was distributed along the Hokkaido coast. A sharp water-mass front along the 33.3 isohaline (white contour line in Figure 6a) was apparent around Cape Erimo between subtropical waters supplied from the Tsugaru Strait and subarctic waters around the Oyashio. The front meandered along the southern edge of the Oyashio.

Relationships to Environmental Conditions
To understand relationships between high SCCs (i.e., potential occurrence of HABs) and oceanographic conditions, we analyzed outputs from the 1/50° ocean model. Monthly mean salinities and current velocities at the sea surface ( Figure 6a) reflected representative ocean current conditions around Hokkaido in September 2021. The Oyashio flowed from the northeast to the southwest along the continental slope. A part of the Oyashio bifurcated to the south of Shikotan Island, was deflected in a clockwise direction, and then rejoined the mainstream of the Oyashio. Meanwhile, the Soya Warm Current, which is characterized by high-salinity water (>33.5), flowed along the Okhotsk coast of northeast Hokkaido and intruded into the Kunashiri Strait. Water with relatively high salinity (~33.5) that outflowed into the Pacific shelf was then transported in a south-westward direction by the Coastal Oyashio onto the Pacific shelf, inshore of the Oyashio on the slope. The salinity along the Coastal Oyashio decreased gradually downstream from the Kunashiri Strait. The implication was that Soya Warm Current water was modified by mixing with the surrounding low-salinity waters, and so-called modified Soya Warm Current water was formed on the Pacific shelf. On the inshore side of the Coastal Oyashio, water with very low salinity (<32.5) due to river discharge was distributed along the Hokkaido coast. A sharp water-mass front along the 33.3 isohaline (white contour line in Figure 6a) was apparent around Cape Erimo between subtropical waters supplied from the Tsugaru Strait and subarctic waters around the Oyashio. The front meandered along the southern edge of the Oyashio. We compared a map of monthly mean SCCs (Figure 6d) with the simulated salinities. Overall, relatively high SCCs tended to be distributed within the inshore side of the Oyashio mainstream (i.e., coastal shelf waters), except in the vicinity of the clockwise bifurcation of the Oyashio on the slope, south of the Shikotan Island. Particularly high SCCs (~10 mg m −3 ) were distributed in very nearshore waters along the southeast Hokkaido coast, where river discharges resulted in very low-salinity waters. The surface mixed layers were also very shallow in those waters (~a few meters) (Figure 6c). It is noteworthy that SCCs along the mainstream of the Soya Warm Current were quite low (<1 mg m −3 ). In addition, around the southern entrance of the Kunashiri Strait and Shikotan Island (see thick dashed ellipse in Figure 6d), where Soya Warm Current water was supplied from the Sea of Okhotsk throughout summer and autumn, SCCs were also lower (<3 mg m −3 ) than in the surrounding Pacific shelf waters. Moreover, there was an abrupt transition from low to high SCCs between the western and eastern side, respectively, of the sharp salinity front south of Cape Erimo. These results suggested that the high SCCs associated with HABs were less intense in subtropical, high-salinity waters.
To obtain a more quantitative understanding of the relationships between SCCs and environmental variables, we created four scatter plots of SCCs versus water depths, sea surface salinities, sea surface temperatures, and surface mixed layer depths. In Figure 7, we have plotted daily SCCs based on the L3 product in only Pacific waters versus the daily means of variables. Bin-averaged SCCs (filled red circles) and the frequency distributions for SCCs exceeding 10 mg m −3 (pink bar chart) overlay the scatter diagrams. Scatter plots based on monthly means (e.g., Figure 6) exhibited similar features, but there were fewer data points (not shown). We compared a map of monthly mean SCCs (Figure 6d) with the simulated salinities. Overall, relatively high SCCs tended to be distributed within the inshore side of the Oyashio mainstream (i.e., coastal shelf waters), except in the vicinity of the clockwise bifurcation of the Oyashio on the slope, south of the Shikotan Island. Particularly high SCCs (~10 mg m −3 ) were distributed in very nearshore waters along the southeast Hokkaido coast, where river discharges resulted in very low-salinity waters. The surface mixed layers were also very shallow in those waters (~a few meters) (Figure 6c). It is noteworthy that SCCs along the mainstream of the Soya Warm Current were quite low (<1 mg m −3 ). In addition, around the southern entrance of the Kunashiri Strait and Shikotan Island (see thick dashed ellipse in Figure 6d), where Soya Warm Current water was supplied from the Sea of Okhotsk throughout summer and autumn, SCCs were also lower (<3 mg m −3 ) than in the surrounding Pacific shelf waters. Moreover, there was an abrupt transition from low to high SCCs between the western and eastern side, respectively, of the sharp salinity front south of Cape Erimo. These results suggested that the high SCCs associated with HABs were less intense in subtropical, high-salinity waters.
To obtain a more quantitative understanding of the relationships between SCCs and environmental variables, we created four scatter plots of SCCs versus water depths, sea surface salinities, sea surface temperatures, and surface mixed layer depths. In Figure 7, we have plotted daily SCCs based on the L3 product in only Pacific waters versus the daily means of variables. Bin-averaged SCCs (filled red circles) and the frequency distributions for SCCs exceeding 10 mg m −3 (pink bar chart) overlay the scatter diagrams. Scatter plots based on monthly means (e.g., Figure 6) exhibited similar features, but there were fewer data points (not shown).  Water depth was a primary determinant of the occurrence of SCCs that exceeded 10 mg m −3 (Figure 7a); both bin-averaged means and frequencies of high SCCs increased toward shallow water depths. 72% of the SCCs that exceeded 10 mg m −3 occurred in water depths of 0-300 m, which corresponded to locations near the coast, shelf, and shelf break. 9% and 19% of the high SCCs were associated with water depths of 300-2000 and >2000 m, respectively. Sea surface salinity and mixed layer depth exhibited similar relationships with SCCs; their frequency distributions exhibited two local maxima. The two maxima could be interpreted to reflect the different conditions associated with high SCCs in nearshore and shelf waters. As mentioned before, the nearshore waters were characterized by very low salinities and very shallow mixed layer depths (Figure 6a,c). Those conditions corresponded to the peaks on the sides of Figure 7b,c. The high SCCs in the shelf waters corresponded to the peaks on the right sides of the same figures at a salinity of 32.8 ( Figure  7b) and mixed layer depth of 10 m (Figure 7c). In addition, high SCCs appeared less frequently in water with a high salinity (>33.3) and at locations where the mixed layer exceeded 20 m. These conditions were associated with subtropical waters (Figure 6a,c). However, part of the modified Soya Warm Current water, typically with a salinity of 33.0-33.7, contributed to the peak at high salinity ~32.8 and a mixed layer depth ~10 m. The frequency distribution of high SCCs as a function of sea surface temperature (Figure 7d) exhibited a single peak at around 16°C and differed distinctly from the bimodal salinity and mixed layer depth distributions. Sea-surface temperatures could therefore not be used to distinguish between the high SCCs in nearshore and shelf waters.
It is here emphasized that sea surface salinities simulated by the 1/50° ocean model had low-salinity bias of about 0.2 in autumn around the focal region (e.g., [25]), which was reconfirmed by comparison with in situ sea surface salinities along the A-line (e.g., [22]) that were measured on 4-10 October 2021 (not shown).
To identify the factors possibly responsible for the anomalous increase in the daily areal extent of high SCCs after late August (i.e., Figure 4), we estimated the mixed layer depths at grid cells with SCCs of 5> or >10 mg m −3 , and averaged those mixed layer depths for each day (Figure 8a-d). We refer to those averages as "high-chlorophyll mixed layer depths." The mixed layer depth is here regarded as a metric of nutrient supply from the subsurface to the vicinity of the sea surface via vertical convection/entrainment. The highchlorophyll mixed layer depths were roughly correlated with the extent of high SCCs (Figure 8b,d). High-chlorophyll mixed layer depths of approximately 10 m corresponded to the largest area, which appeared in mid-to-late September (Figure 8a,c). A mixed layer Water depth was a primary determinant of the occurrence of SCCs that exceeded 10 mg m −3 (Figure 7a); both bin-averaged means and frequencies of high SCCs increased toward shallow water depths. 72% of the SCCs that exceeded 10 mg m −3 occurred in water depths of 0-300 m, which corresponded to locations near the coast, shelf, and shelf break. 9% and 19% of the high SCCs were associated with water depths of 300-2000 and >2000 m, respectively. Sea surface salinity and mixed layer depth exhibited similar relationships with SCCs; their frequency distributions exhibited two local maxima. The two maxima could be interpreted to reflect the different conditions associated with high SCCs in nearshore and shelf waters. As mentioned before, the nearshore waters were characterized by very low salinities and very shallow mixed layer depths (Figure 6a,c). Those conditions corresponded to the peaks on the sides of Figure 7b,c. The high SCCs in the shelf waters corresponded to the peaks on the right sides of the same figures at a salinity of 32.8 ( Figure 7b) and mixed layer depth of 10 m (Figure 7c). In addition, high SCCs appeared less frequently in water with a high salinity (>33.3) and at locations where the mixed layer exceeded 20 m. These conditions were associated with subtropical waters (Figure 6a,c). However, part of the modified Soya Warm Current water, typically with a salinity of 33.0-33.7, contributed to the peak at high salinity~32.8 and a mixed layer depth 10 m. The frequency distribution of high SCCs as a function of sea surface temperature ( Figure 7d) exhibited a single peak at around 16 • C and differed distinctly from the bimodal salinity and mixed layer depth distributions. Sea-surface temperatures could therefore not be used to distinguish between the high SCCs in nearshore and shelf waters.
It is here emphasized that sea surface salinities simulated by the 1/50 • ocean model had low-salinity bias of about 0.2 in autumn around the focal region (e.g., [25]), which was reconfirmed by comparison with in situ sea surface salinities along the A-line (e.g., [22]) that were measured on 4-10 October 2021 (not shown).
To identify the factors possibly responsible for the anomalous increase in the daily areal extent of high SCCs after late August (i.e., Figure 4), we estimated the mixed layer depths at grid cells with SCCs of >5 or >10 mg m −3 , and averaged those mixed layer depths for each day (Figure 8a-d). We refer to those averages as "high-chlorophyll mixed layer depths". The mixed layer depth is here regarded as a metric of nutrient supply from the subsurface to the vicinity of the sea surface via vertical convection/entrainment. The high-chlorophyll mixed layer depths were roughly correlated with the extent of high SCCs (Figure 8b,d). High-chlorophyll mixed layer depths of approximately 10 m corresponded to the largest area, which appeared in mid-to-late September (Figure 8a,c). A mixed layer depth of 10 m was also related to the right peak of the frequency distribution of high SCCs as a function of mixed layer depth (Figure 7c). ci. Eng. 2021, 9,1335 11 of 21 depth of 10 m was also related to the right peak of the frequency distribution of high SCCs as a function of mixed layer depth (Figure 7c). The period of the most extensive and intense marine heatwaves ("MHWs period") is based on Kuroda and Setou [14] and is indicated by a double-headed arrow.
The high-chlorophyll mixed layer depths continued to be very shallow until the beginning of August (Figure 8a,c), when the marine heatwaves reached a maximum ( Figure  8e), and they deepened abruptly on 10 August, when the marine heatwaves retreated abruptly. After mid-August, when the marine heatwaves ended, the high-chlorophyll mixed layer depths tended to increase gradually to the end of September. Simulated vertical eddy diffusivity near the sea surface averaged over the focal region also continued to increase after the spike-like increase on 10 August (Figure 8e). It is noteworthy that the area of high SCCs increased slowly after mid-August, started to intermittently exceed the climatological daily maximum after late August, and reached a record-breaking extent in mid-to-late September (Figure 4). As a result, the extreme expansion of the area of high SCCs associated with the occurrence of massive HABs was tightly linked with deepening of mixed layer depths in subarctic waters beginning at the time of the abrupt end of the marine heatwaves.
However, the deepening of the mixed layer depths could not explain the abrupt increase in SCCs (i.e., intensification of HABs) in very nearshore waters along the southeast The high-chlorophyll mixed layer depths continued to be very shallow until the beginning of August (Figure 8a,c), when the marine heatwaves reached a maximum (Figure 8e), and they deepened abruptly on 10 August, when the marine heatwaves retreated abruptly. After mid-August, when the marine heatwaves ended, the high-chlorophyll mixed layer depths tended to increase gradually to the end of September. Simulated vertical eddy diffusivity near the sea surface averaged over the focal region also continued to increase after the spike-like increase on 10 August (Figure 8e). It is noteworthy that the area of high SCCs increased slowly after mid-August, started to intermittently exceed the climatological daily maximum after late August, and reached a record-breaking extent in mid-to-late September (Figure 4). As a result, the extreme expansion of the area of high SCCs associated with the occurrence of massive HABs was tightly linked with deepening of mixed layer depths in subarctic waters beginning at the time of the abrupt end of the marine heatwaves.
However, the deepening of the mixed layer depths could not explain the abrupt increase in SCCs (i.e., intensification of HABs) in very nearshore waters along the southeast Hokkaido coast (Figure 3d-f), where a very shallow mixed layer depth was related to river-influenced, very-low-salinity water that contained high SCCs (e.g., Figure 7b,c). River discharge into coastal waters were examined as a potential controlling factor (e.g., [15]), but there was little evidence of abnormal discharge conditions during the summer of 2021 (see Appendix A). As a result, the mechanism responsible for the abrupt growth of the harmful algae and a mechanism responsible for their accumulation (e.g., biological aggregation or physical convergence) needed to be identified.

Potential Sources of HABs
To identify the area that was likely the source of the HABs, we first performed retrospective particle-tracking simulations (Figure 9). Some of the particles were transported from their positions in the blooms backward in time along the main stream of the Coastal Oyashio and Soya Warm Current. More specifically, particles were tracked from their positions in the blooms in an eastward direction along the Pacific shelf (Figure 9a,b), and then from the Pacific shelf through the Kunashiri Strait to the shelf along the north-eastern coast of Hokkaido in the Okhotsk Sea (Figure 9c,d). The particles then passed through the Soya Strait ( Figure 9e) and finally reached the Sea of Japan (Figure 9f). The minimum elapsed time was 21-28 days, at least for particles to be transported retrospectively from their positions in the blooms to the Sea of Japan. Hence, one of the potential source areas was the Sea of Japan. Shimada et al. [49] have identified the presence of the dinoflagellate K. mikimotoi near the Hokkaido coast in the Tsugaru Strait during autumn of 2015 for the first time and suggested that there was horizontal transport of the harmful algae from the Sea of Japan. Our results were consistent with the potential source area that they suggested.
Hokkaido coast (Figure 3d-f), where a very shallow mixed layer depth was related to river-influenced, very-low-salinity water that contained high SCCs (e.g., Figure 7b,c). River discharge into coastal waters were examined as a potential controlling factor (e.g., [15]), but there was little evidence of abnormal discharge conditions during the summer of 2021 (see Appendix A). As a result, the mechanism responsible for the abrupt growth of the harmful algae and a mechanism responsible for their accumulation (e.g., biological aggregation or physical convergence) needed to be identified.

Potential Sources of HABs
To identify the area that was likely the source of the HABs, we first performed retrospective particle-tracking simulations (Figure 9). Some of the particles were transported from their positions in the blooms backward in time along the main stream of the Coastal Oyashio and Soya Warm Current. More specifically, particles were tracked from their positions in the blooms in an eastward direction along the Pacific shelf (Figure 9a,b), and then from the Pacific shelf through the Kunashiri Strait to the shelf along the north-eastern coast of Hokkaido in the Okhotsk Sea (Figure 9c,d). The particles then passed through the Soya Strait (Figure 9e) and finally reached the Sea of Japan (Figure 9f). The minimum elapsed time was 21-28 days, at least for particles to be transported retrospectively from their positions in the blooms to the Sea of Japan. Hence, one of the potential source areas was the Sea of Japan. Shimada et al. [49] have identified the presence of the dinoflagellate K. mikimotoi near the Hokkaido coast in the Tsugaru Strait during autumn of 2015 for the first time and suggested that there was horizontal transport of the harmful algae from the Sea of Japan. Our results were consistent with the potential source area that they suggested.
Next, we examined whether harmful algae could have been passively transported by ocean currents to the focal region from the Pacific Ocean east of the Kamchatka Peninsula, where a massive bloom of the dinoflagellate K. selliformis occurred during September-October 2020 [18], with a focus on wider area and longer tracking time ( Figure 10). Particles that were transported retrospectively to the northeast along the Oyashio and the East Kamchatka Current travelled the longest distance. It took about 4 months, at least, for the particles to be transported retrospectively from coastal waters off southeast Hokkaido to the area east of the Kamchatka Peninsula (Figure 10d-f). Next, we examined whether harmful algae could have been passively transported by ocean currents to the focal region from the Pacific Ocean east of the Kamchatka Peninsula, where a massive bloom of the dinoflagellate K. selliformis occurred during September-October 2020 [18], with a focus on wider area and longer tracking time ( Figure 10). Particles that were transported retrospectively to the northeast along the Oyashio and the East Kamchatka Current travelled the longest distance. It took about 4 months, at least, for the particles to be transported retrospectively from coastal waters off southeast Hokkaido to the area east of the Kamchatka Peninsula (Figure 10d-f).  Moreover, note that a large number of particles were transported back to the southern part of the Sea of Okhotsk ( Figure 10). In March 2021, particles were densely distributed along the east coast of the southern part of the Sakhalin Island (Figure 10f), which indicated that particles were prospectively transported southward by the East Sakhalin Current intensified in winter. Orlova et al. [50] also reported that K. mikimotoi was identified during June-October in 1994 and 1998-2000 near the east coast of the Sakhalin coast between 46°30′ N and 53° N. Therefore, the vicinity of this area could be also regarded as potential source areas.
The numbers of particles passing through the following four transects were counted at daily intervals ( Figure 11a): (1) a meridional transect in the Tsugaru Strait, (2) a meridional transect in the Soya Strait, (3) a zonal transect from 143°30′ E to 146° E east of the Sakhalin Island, and (4) a zonal transect east of the Kamchatka Peninsula (see Figure 10a). For particles passing through the transect in the Soya Strait, the number of particles started to increase after the 20-day tracking period, reached a maximum during the 30-60-day tracking period (i.e., July-August), and then decreased abruptly. Meanwhile, the number of particles passing through the transect in the Tsugaru Strait was nearly zero throughout the particle-tracking period. The number of particles passing through the transect east of the Sakhalin Island increased remarkably from 135-day tracking period (i.e., mid-April). This result indicated that particles prospectively left the east coast of the Sakhalin Island in winter, when particles were actively transported southward together with sea ice-related very cold waters with a nearly freezing temperature by the East Sakhalin Current. Moreover, the number of particles passing through the transect east of the Kamchatka Peninsula was zero until roughly the 100-day tracking period (i.e., mid-May). It then increased and fluctuated. Hence, if the ocean east of the Kamchatka Peninsula had been the Moreover, note that a large number of particles were transported back to the southern part of the Sea of Okhotsk ( Figure 10). In March 2021, particles were densely distributed along the east coast of the southern part of the Sakhalin Island (Figure 10f), which indicated that particles were prospectively transported southward by the East Sakhalin Current intensified in winter. Orlova et al. [50] also reported that K. mikimotoi was identified during June-October in 1994 and 1998-2000 near the east coast of the Sakhalin coast between 46 • 30 N and 53 • N. Therefore, the vicinity of this area could be also regarded as potential source areas.
The numbers of particles passing through the following four transects were counted at daily intervals ( Figure 11a): (1) a meridional transect in the Tsugaru Strait, (2) a meridional transect in the Soya Strait, (3) a zonal transect from 143 • 30 E to 146 • E east of the Sakhalin Island, and (4) a zonal transect east of the Kamchatka Peninsula (see Figure 10a). For particles passing through the transect in the Soya Strait, the number of particles started to increase after the 20-day tracking period, reached a maximum during the 30-60-day tracking period (i.e., July-August), and then decreased abruptly. Meanwhile, the number of particles passing through the transect in the Tsugaru Strait was nearly zero throughout the particle-tracking period. The number of particles passing through the transect east of the Sakhalin Island increased remarkably from 135-day tracking period (i.e., mid-April). This result indicated that particles prospectively left the east coast of the Sakhalin Island in winter, when particles were actively transported southward together with sea ice-related very cold waters with a nearly freezing temperature by the East Sakhalin Current. Moreover, the number of particles passing through the transect east of the Kamchatka Peninsula was zero until roughly the 100-day tracking period (i.e., mid-May). It then increased and fluctuated. Hence, if the ocean east of the Kamchatka Peninsula had been the source of the HABs off southeast Hokkaido, a transport period of more than four months would have been needed.
source of the HABs off southeast Hokkaido, a transport period of more than four months would have been needed. To investigate linkages between HABs that developed east of the Kamchatka Peninsula in late summer 2020 to those that developed in the focal region in late summer 2021, we analyzed results from particle-tracking simulations forward in time that released a large number of particles from the vicinity of the Pacific Ocean east of the Kamchatka Peninsula during 30 September to 20 October 2020 ( Figure 12). Most particles started to be transported to the southwest at the beginning of the particle-tracking period ( Figure  12a). A large percentage of the particles were then greatly dispersed horizontally and transported outside of the particle-tracking domain through the eastern boundary until the tracking time reached 160 days (approximately mid-March 2021) (Figure 12b,c). The rest of the particles were partially transported south-westward along the East Kamchatka Current and the Oyashio, and some of them reached the focal region of the study (Figure  12d-f). In July-September 2021, only a small percentage of the particles (~0.3% relative to the number of initially released particles) remained in the focal region ( Figure 13). These results, therefore, implied that growth of harmful algae during the period of passive transport would have been essential if HABs in the Pacific Ocean to the east of the Kamchatka Peninsula in the late summer of 2020 were the cause of the blooms in the focal region in the late summer of 2021, but that possibility was not included in our particletracking simulation. To investigate linkages between HABs that developed east of the Kamchatka Peninsula in late summer 2020 to those that developed in the focal region in late summer 2021, we analyzed results from particle-tracking simulations forward in time that released a large number of particles from the vicinity of the Pacific Ocean east of the Kamchatka Peninsula during 30 September to 20 October 2020 ( Figure 12). Most particles started to be transported to the southwest at the beginning of the particle-tracking period (Figure 12a). A large percentage of the particles were then greatly dispersed horizontally and transported outside of the particle-tracking domain through the eastern boundary until the tracking time reached 160 days (approximately mid-March 2021) (Figure 12b,c). The rest of the particles were partially transported south-westward along the East Kamchatka Current and the Oyashio, and some of them reached the focal region of the study (Figure 12d-f). In July-September 2021, only a small percentage of the particles (~0.3% relative to the number of initially released particles) remained in the focal region ( Figure 13). These results, therefore, implied that growth of harmful algae during the period of passive transport would have been essential if HABs in the Pacific Ocean to the east of the Kamchatka Peninsula in the late summer of 2020 were the cause of the blooms in the focal region in the late summer of 2021, but that possibility was not included in our particle-tracking simulation.

Possible Linkages from Algal Transport to Outbreak
Particle-tracking simulations identified three potential source areas: the Pacific Ocean east of the Kamchatka Peninsula, the Sea of Japan, and the Sea of Okhotsk east of the Sakhalin Island. It is noteworthy that the high SCCs apparent around the focal region in late summer 2021 were isolated from the surrounding waters of the Pacific Ocean and Sea of Okhotsk (e.g., Figure 3). High SCCs were also not observed along the transport pathways from the potential sources to the focal region. These facts allowed us to hypothesize that harmful algae or their seeds were transported from the potential source areas to the focal region without generating blooms, and the algae finally generated massive blooms within the focal region when they encountered optimal environmental conditions. This hypothesis is probably reasonable for the algae transported from the Sea of Japan, because our results indicated that high SCCs were less frequently observed in subtropical, high-salinity waters (e.g., Figure 7b). Along the transport pathway from the Sea of Japan to the focal region, harmful algae or seeds could be transported together with oligotrophic subtropical waters by the Tsushima Warm Current and Soya Warm Current. Figure 13. The number of particles that were distributed within the focal region of the study in prospective particle-tracking simulations. The left axis is the ratio of the number of particles to the total number of released particles (i.e., 1,314,240).

Possible Linkages from Algal Transport to Outbreak
Particle-tracking simulations identified three potential source areas: the Pacific Ocean east of the Kamchatka Peninsula, the Sea of Japan, and the Sea of Okhotsk east of the Sakhalin Island. It is noteworthy that the high SCCs apparent around the focal region in late summer 2021 were isolated from the surrounding waters of the Pacific Ocean and Sea of Okhotsk (e.g., Figure 3). High SCCs were also not observed along the transport pathways from the potential sources to the focal region. These facts allowed us to hypothesize that harmful algae or their seeds were transported from the potential source areas to the focal region without generating blooms, and the algae finally generated massive blooms within the focal region when they encountered optimal environmental conditions. This hypothesis is probably reasonable for the algae transported from the Sea of Japan, because our results indicated that high SCCs were less frequently observed in subtropical, high-salinity waters (e.g., Figure 7b). Along the transport pathway from the Sea of Japan to the focal region, harmful algae or seeds could be transported together with oligotrophic subtropical waters by the Tsushima Warm Current and Soya Warm Current.
The cases for the Sea of Okhotsk east of the Sakhalin Island and the Pacific Ocean east of the Kamchatka Peninsula are more complex, principally because the transport time of the algae would have been longer (>4 months), and the algae would have been transported together with subarctic waters by the East Sakhalin Current, East Kamchatka Current, and Oyashio. The long transport time forced us to consider biological growth of the harmful algae, as mentioned in Section 3.3. If these two pathways were feasible, the subarctic waters along the transport pathways might have suppressed the growth of the algae. In general, subarctic waters exhibit nutrient-rich conditions, associated with enhancement of algal growth. Therefore, light-deficient associated with a persistently deep mixed layer induced by strong tidal mixing [51,52], and/or reduction in growth attributed to very cold temperatures could have inhibited outbreaks of harmful algae during the transport period.
We want to stress that ocean circulation near the sea surface in the focal region is characterized by high retention of particles. As shown in Figure 11b, a large number of particles remained within the focal region, even after a long period of particle tracking. More specifically, about 60%, 40%, and 20% of the particles initially released remained within the focal region after retrospective tracking periods of 60, 120, and 180 days, respectively. In addition, the number of particles distributed within the focal region during the prospective particle-tracking simulation became nearly constant after June of 2021 ( Figure 12). The implication is that the focal region very efficiently retained particles. In other words, once particles were transported into the focal region, it was difficult for them to exit. The high particle retention by the focal region might have allowed the harmful algae to remain in the focal region until environmental conditions became favourable for massive blooms.
HABs that included both K. mikimotoi and K. selliformis have been identified off the Hokkaido coast. According to Shimada et al. [49], a potential source of the K. mikimotoi distributed around Hokkaido is the Sea of Japan. Orlova et al. [50] also identified K. mikimotoi near the east coast of the Sakhalin Island. These results are consistent with two potential source areas identified by our particle-tracking simulations. Moreover, our particle-tracking simulations revealed that the Pacific Ocean to the east of the Kamchatka Peninsula could be the other potential source. This possibility is consistent with the outbreaks of K. selliformis in September-October 2020 [18,19]. These observations lead to the hypothesis that different algal species may have been separately transported from the three potential source areas and contributed simultaneously to the occurrence of the HABs off the southeast coast of Hokkaido in late summer of 2021.

Relationship of Algal Outbreaks to Marine Heatwaves
It is very instructive to remember that the massive HABs in the Pacific Ocean near the east coast of the Kamchatka Peninsula in the late summer of 2020 followed extreme marine heatwaves that occurred in both June and July of 2020, when the maximum anomalies of sea surface temperature reached about 6 • C [18]. The magnitude and timing of the second marine heatwave and the time lag between the second marine heatwave and the HABs were very consistent with the sequence of events in this study. In the case of the Pacific Ocean to the east of the Kamchatka Peninsula, a large influx of nutrients came from river discharge. In our focal region, entrainment of subsurface waters by deepening of the mixed layer over the shelf was an important source of nutrients. The 10 m depth of the mixed layer in the low-salinity, subarctic waters likely provided optimal conditions for blooming of harmful algae during September 2021 (Figure 7b,c), i.e., there was a sufficient supply of both nutrients and light. The importance of uplifted nutrients to the occurrence of HABs has likewise been stressed by McCabe et al. [53], who attributed large HABs along the North American west coast associated with marine heatwaves to seasonal coastal upwelling.
Numerous studies have also reported that marine heatwaves have triggered massive HABs in many places globally [18,[53][54][55][56][57][58][59][60][61]. However, the process and mechanism, namely, how marine heatwaves trigger HABs, has not been accurately understood. Additional studies for clarification are thus essential to mitigate the impacts of HABs on human life. It is difficult, however, to repeatedly monitor HABs triggered by marine heatwaves in a specific ecosystem, because marine heatwaves are climate extremes (i.e., rare events) in the ocean systems. Interregional studies are the best way to compare HABs triggered by marine heatwaves among different marine ecosystems and to understand the relevant mechanisms and predict the occurrence of HABs.

Conclusions
How can we most efficiently develop an understanding of the impact of extreme events on marine ecosystems and relate those impacts to environmental conditions? One of the most robust methodologies is analysis by combining near-real-time ocean monitoring with realistic ocean modelling. Unprecedented outbreaks of harmful algae were reported in the late summer of 2021 in coastal waters off southeast Hokkaido, one month after the most intense and extensive marine heatwaves in the north-western Pacific Ocean had subsided. Just after the marine heatwaves subsided (~20 August), the areal extent of the high SCCs began to slowly increase. The area of high SCCs occasionally exceeded the climatological maxima after late August and finally reached the most extreme record-breaking level (~mean plus 6 standard deviations) in mid-to-late September. The expansion of the area was accompanied by a deepening of surface mixed layer depths that began after the terminal phase of the marine heat waves. That deepening of the mixed layer was associated with an influx of nutrients from the subsurface via vertical convection/entrainment. About 70% of the high SCCs that exceeded 10 mg m −3 were observed in coastal shelf waters with water depths < 300 m. Characteristics of the high SCCs clearly differed between nearshore waters and shelf waters. The nearshore waters were strongly influenced by river discharge, i.e., the salinity was very low, and the mixed layer was very shallow (~a few meters). The waters over the shelf with high SCCs corresponded to subarctic waters with mixed layer depths of about 10 m. In other words, the high SCCs rarely occurred in subtropical waters.
Retrospective particle-tracking simulations revealed three potential source areas of the HABs that occurred in the coastal waters of the Pacific Ocean off southeast Hokkaido: the Pacific Ocean to the east of the Kamchatka Peninsula, the Sea of Japan, and the Sea of Okhotsk east of the Sakhalin Island. The transport time from the potential source areas differed greatly: about 1-2 months for the second potential source area and >4 months for the others. Prospective particle-tracking simulations also indicated that when particles were released from coastal waters off the east coast of the Kamchatka Peninsula in October 2020 on the assumption that a massive HAB had occurred there,~0.3% of the particles initially released reached the focal region in July-September of 2021 and stagnated there because the focal area seems to very efficiently retain particles. The very small percentage of~0.3% implied that algal growth during the period of passive transport over the long distance would have been essential to cause the subsequent massive HABs, although that possibility was completely neglected in our simulations. In any case, it is needed for future work to verify the validity of each potential source area by examining in situ distribution of Karenia species in 2020-2021 around the three potential transport pathways.
Finally, we want to note two weaknesses of this study. The first weakness was that the study focused merely on the development stage of the HABs. The HABs seemed to be declining as of 15 November 2021. The processes associated with this decline need to be examined in subsequent work. The second weakness of this study was that it included some biases. According to Kuroda et al. [25], for instance, sea surface salinity estimated with a 1/50 • ocean model in coastal shelf waters off the southeast coast of Hokkaido in autumn 2016 exhibited a low-salinity bias by about 0.2. This bias suggested that the biological-physical relationships revealed by this study should be based on relative values of environmental conditions (e.g., Figure 7). These model biases need to be corrected by comparisons with in situ measurements in future work to precisely understand biological phenomena in relation to absolute values of environmental variables. In addition, the SCCs are associated with interpretative biases. Specifically, high SCCs might not be entirely attributable to the occurrence of HABs. It is possible that non-harmful algae might increase in abundance, particularly in offshore waters. Our study should, therefore, be confirmed and could be improved by accumulating scientific data based on in situ measurements. Nevertheless, we emphasize that analyses combining satellite monitoring and realistic ocean modelling are very useful for immediately overviewing the entirety of a phenomenon, even if it is an unprecedented, extreme event.  Data Availability Statement: Gridded data of chlorophyll concentrations from ocean color sensors were downloaded from a webserver of the European Union's Copernicus Marine Service (accessed on 12 October 2021). Daily reanalysis data from FRA-ROMS are available at http://fm.dc.affrc.go. jp/fra-roms after user registration (accessed on 12 October 2021). The source code of the 1/50 • model based on Rutgers University's ROMS is available at https://www.myroms.org (accessed on 17 February 2010). All external forcings for the 1/50 • model (NAO.99Jb, GPV/MSM, JRA55, and river discharge data) are also publicly available. The original LTRANS code was downloaded from https://northweb.hpl.umces.edu/LTRANS.htm (accessed on 4 September 2012).