After mid-September 2021 (when Karenia blooms were first observed), surface water samples were collected daily around the Katsurakoi fishery harbor (144°26.77′E, 42°56.81′N; indicated by the star symbol in Figure 1B). Samples were collected from three fixed sites within a few hundred meters of the harbor. Buckets were used for sample collection at two of the sites, and an automated seawater intake to our shore-based laboratory was used at the other site. Water samples were filtered onto 25-mm Whatman GF/F glass fiber filters immediately after sampling to measure chlorophyll-a concentrations with a fluorometer. For enumeration of Karenia spp., water samples were fixed and preserved at 5°C with acid Lugol’s solution (final concentration, 4%) for measurements conducted prior to 28 September, and with Hepes-buffered paraformaldehyde and glutaraldehyde (Katano et al., 2009) for subsequent measurements. Karenia spp. abundances (cells mL^–1) were estimated from 46 randomly selected samples collected during 21 September–22 October, 2021.

Ocean color imagery based on Sentinel-3A/3B Ocean and Land Color Imager Level-2 Full Resolution were downloaded from Copernicus Online Data Access. The horizontal resolution was about 300 m, and the time interval between images was near-daily. The data period analyzed was from 3 October to 18 October, 2021 (Supplementary Table 1), roughly corresponding to the time period of the ship survey. For analysis, we selected imagery with a relatively small ratio of cloud coverage (< 20%) over the shelf off southeastern Hokkaido. We used normalized water-leaving reflectance (ρ_w[λ]) for a wavelength band centered at λ and chlorophyll-a concentrations that were derived from one of two algorithms: a neural network-based approach (hereafter, “NN”) for optically complex waters (Doerffer and Schiller, 2007) and a maximum band ratio semi-analytical algorithm (hereafter, “OC4ME”) for open oceans (Morel et al., 2007; Cherif et al., 2021; Moutzouris-Sidiris and Topouzelis, 2021).

Chlorophyll-a estimates obtained from NN and OC4ME had limited accuracy on the Pacific shelf (see section “Maps of Chlorophyll-a Concentration” for details), where there was considerable submesoscale variability in areas of high chlorophyll-a concentration (>20 mg m^–3). Therefore, reflectance-based maximum line height (MLH; Smith and Bernard, 2020) was used instead of NN and OC4ME to estimate chlorophyll-a concentrations at the sea surface.

Surface abundances of Karenia spp. were estimated from satellite-derived reflectance-based variables by using two methods. One method directly estimated Karenia spp. abundances through a combination of reflectance-based red band difference (RBD; Jordan et al., 2021) and in-situ Karenia spp. abundance. The other method entailed two steps: in the first step, chlorophyll-a concentrations were estimated by using a combination of MLH and in-situ chlorophyll a concentrations, as mentioned above, and in the second step, chlorophyll-a concentrations were converted into Karenia spp. abundances by using an observed relationship between them.

RBD and MLH are defined as follows:

where LH(λ) is defined as follows:

Relationships between the spectral properties of Karenia spp. and the Sentinel-3 spectral bands are illustrated in the schematic diagram in Supplementary Figure 1. Calculations of RBD and MLH were both based on the 665–753-nm spectral band, as opposed to OC4ME, which is based on the 443–560-nm spectral band. The 665–753-nm spectral band has an important advantage over the spectral band used in OC4ME in that the former is not strongly influenced by the absorption associated with colored dissolved organic matter (e.g., at around 440 nm), which is frequently found in the coastal waters of our study area (e.g., Isada et al., 2021).

Red band difference is frequently used for detecting and monitoring Karenia spp. blooms (e.g., Amin et al., 2009; Wolny et al., 2020; Jordan et al., 2021). LH(λ), which is similar to fluorescence line height (Gower et al., 1999), with a baseline formed by the reflectance between 665 and 753 nm, was applied to all spectra (Smith and Bernard, 2020). MLH was expected to be appropriate for detecting both low (<about 20 mg m^–3) chlorophyll-a concentrations and the high (>about 20 mg m^–3) chlorophyll-a concentrations associated with red tides, because the reflectance peak at 681 nm is generally related to chlorophyll-a fluorescence emissions; however, at higher biomass this peak shifts to longer wavelengths owing to the combined effects of increased phytoplankton absorption and backscattering, as well as pure water absorption (e.g., Gower et al., 1999, 2005; Gower, 2016).

Finally, after generating maps of Karenia spp. abundance, we estimated specific rates of change (hereafter, “growth rate”), assuming exponential algal growth, from two maps of satellite-derived Karenia spp. abundance. Specifically, growth rate was estimated as r = [ln (K_t2) − ln (K_t1)]/△t, where △t is the time interval between satellite images and K_t1 and K_t2 represent Karenia spp. abundances estimated in a given region at times t1 and t2, respectively (i.e., △t = t 2− t1). This estimation method was based on the work of Stumpf et al. (2008), who developed the method for satellite-derived chlorophyll-a imagery.

The oceanographic conditions in our study area during the cruise can be seen from satellite-based data captured on 9 October, 2021 (Figure 1). The northern part of the A Line was occupied by cold waters, and the southern part was occupied by warm waters (Figure 1A). In the cold-water area, the Oyashio flowed southwestward along the continental slope off southeastern Hokkaido, came in contact with the Tsugaru Warm Current south of Cape Erimo, formed a sharp thermal front, and then flowed along the edge of anticyclonic mesoscale eddies, with several meanders (Figure 1B). A part of the Oyashio then crossed the A Line again near Stn. A15.

Inshore of the Oyashio, the Coastal Oyashio flowed along the Pacific shelf off southeastern Hokkaido as a downstream extension of the Soya Warm Current in the Sea of Okhotsk. The Coastal Oyashio crossed several L lines. As explained in the subsequent section “Field Measurements,” the Coastal Oyashio transported Modified Soya Warm Current Water, which is formed by the mixing of a subtropical water mass referred to as Soya Warm Current Water that outflows from the Sea of Okhotsk with surrounding subarctic waters on the Pacific shelf.

In terms of large-scale variability along the shelf (i.e., on a scale of ∼500 km), as briefly mentioned in the section “Introduction,” the areal extent of SCCs exceeding 10 mg m^–3 around the study region (Figure 1B, blue dashed square) reached a record-breaking level in mid-to-late September 2021 (Figure 2), which was roughly comparable to the climatological mean plus six standard deviations. During the first half of our study period (i.e., 3–18 October), this record-breaking level was further maintained. The maximum areal extent of SCCs in 2021 was observed on 9 October. However, during the second half of our study, the areal extent decreased rapidly. Hence, we examined both the maintenance and decay periods of algal blooms.

Along the ship track, temperatures, salinities, and chlorophyll-a concentrations changed dramatically as the research vessel repeatedly crossed subarctic and subtropical water masses, as well as water masses created by mixture of the two (Figure 3). Roughly speaking, cold (warm) waters corresponded to low-salinity (high-salinity) waters. Cold low-salinity subarctic (warm high-salinity subtropical) waters corresponded to high (low) concentrations of chlorophyll a with high (low) spatial variability. Along the ship track around the L lines on the Pacific shelf during 10–14 October, chlorophyll-a concentrations reached particularly high values, typically >10 mg m^–3, intermittently exceeded 20 mg m^–3, and even reached up to 35 mg m^–3. High chlorophyll-a concentrations >20 mg m^–3 were distributed as patch-like local maxima along the ship track (Figure 4). The spatial scale of local maxima >20 mg m^–3 was not more than about a few kilometers and was therefore associated with submesoscale variability (Figure 4B). In addition, high chlorophyll-a concentrations roughly corresponded to high abundances of Karenia spp., as is described below.

On the Pacific shelf around the L lines, salinities were persistently lower than 33.6, and temperatures were in the range 11.0–17.0°C (Figure 3). These water masses were classified as Modified Soya Current Water (generally defined by temperatures >7.0°C and salinities of 33.0–33.7 on the Pacific shelf off southeastern Hokkaido) and Surface Coastal Oyashio Water (temperatures >2.0°C and salinities of 32.0–33.0) (Oguma et al., 2008; Kusaka et al., 2013). Hence, subarctic water masses and modified subtropical waters with particularly high chlorophyll-a concentration were present simultaneously on the Pacific shelf.

Surface chlorophyll-a concentrations were examined in relation to three environmental variables: water depth, sea surface salinity, and temperature (Figure 5). Some 83% of chlorophyll-a concentration measurements exceeding 10 mg m^–3 were obtained in places where the water depth was <1,000 m, i.e., primarily coastal shelf–slope waters (Figures 5A,B). These high chlorophyll-a concentrations were also associated with two different surface salinities (Figure 5C), with local peaks at salinities of 32.9–33.0 and 33.4–33.5. The lower and higher salinity peaks corresponded to Surface Coastal Oyashio Water and Modified Soya Warm Current Water, respectively. However, the high chlorophyll-a concentrations in the two water masses were not associated with different sea-surface temperatures (i.e., the frequency distribution of temperature showed a single mode) (Figure 5D). Chlorophyll-a concentrations for salinities >33.7 and temperatures >17.5°C were mostly lower than 3 mg m^–3, which suggests that massive algal blooms were much less common in pure subtropical waters. The above relationships were qualitatively consistent with those of SSCs as a function of simulated variables (Kuroda et al., 2021a), except for the marked in-situ discrimination between Surface Coastal Oyashio Water and Modified Soya Warm Current Water, and missing measurements of river-discharge-influenced waters with salinity <32.0 near the coast.

At CTD stations, vertical profiles of chlorophyll-a concentration tended to show a vertical maximum at the surface (Figure 6A). Also, at most of the CTD stations, the surface mixed-layer depth ranged from 1 to 54 m (Figure 6B). The modal mixed-layer depth was 20–25 m. Namely, algal blooms were not homogeneous within the surface mixed layer but were intensified at the sea surface. Mixed-layer depths were weakly correlated with the surface chlorophyll-a concentration (Figure 6C). This indicates that satellite measurements of the sea surface could appropriately capture the distribution of algal blooms during the period of the ship survey, despite the seasonal development of a surface mixed layer.

Interestingly, chlorophyll-a concentrations at depths of 0 and 10 m were correlated with each other (Figure 6D). The ratio of chlorophyll a at 0 m to that at 10 m was averaged over all profiles and estimated as 2.4 (Figure 6A); this value was used to scale Karenia spp. abundances that were sampled at 10-m depth.

Karenia spp. abundance at 10-m depth was examined in relation to the surface chlorophyll-a concentration along the ship track (Figure 4, open red circles). Measured abundances ranged from 0 to 277 cells mL^–1 during the period of the ship survey. High Karenia spp. abundance was most common in the northern part of the A Line and along sections of the L lines with high chlorophyll-a concentrations. Meanwhile, Karenia spp. were wholly absent in the southern A Line (i.e., Stns. A09–A21), which was occupied by subtropical waters, except at Stn. A15. At 10-m depth, Karenia spp. abundance was positively correlated with the chlorophyll-a concentration (Figure 7, closed red circles) after log₁₀(x + 0.1) transformation of both variables. This indicates that, to a first-order approximation, chlorophyll-a concentration can be regarded as an index of Karenia spp. abundance off southeastern Hokkaido, at least during autumn 2021.

To obtain a more robust statistical relationship between chlorophyll-a concentration and Karenia spp. abundance, we also examined the data collected from the Katsurakoi fishery harbor (Figure 7, closed blue circles). Karenia spp. abundances (chlorophyll-a concentrations) from the harbor ranged from 0 to 9,600 cells mL^–1 (1.75–123.20 mg m^–3), which was markedly higher than those obtained from the ship survey. By combining the two data sources, we obtained a robust log–log regression line spanning a wide range of values of both variables (R² = 0.75).

In this study, we generated maps of estimated Karenia spp. abundance by using two methods, the first based on RBD and the second on MLH. The MLH-derived abundances were about twice as high as those derived from RBD. First, we discuss the possible reasons for this difference below.

There are at least three weaknesses associated with the RBD-based estimation, in which RBD was directly converted into Karenia spp. abundance. First, in-situ Karenia spp. abundances measured during our ship survey ranged from 0 to 277 cells mL^–1 (Figure 11). Therefore, RBD-derived estimates of abundance >277 cells mL^–1 (before conversion to abundance at 0-m depth) were determined by extrapolating the regression line between RBD and in-situ Karenia spp. abundance (Figure 11). This extrapolation is likely associated with some degree of error. Second, the sample size (n = 68) used to estimate the regression line itself was somewhat small (Figure 11), especially in comparison with the sample size of surface chlorophyll-a concentrations along the ship track (Figures 8, 10). Third, the period of observation differed between Karenia spp. abundances determined from water samples and those determined from satellite measurements. We permitted a maximum time difference of 6.9 days for the A Line and 1 day for the L lines (Figure 11); if we had used a smaller maximum time difference, the number of data points available would have been further reduced.

Likewise, there are at least two weaknesses associated with the MLH-based estimation, in which MLH was used indirectly to estimate Karenia spp. abundance. First, although small-scale variations of chlorophyll-a concentrations on the shelf were reasonably well reproduced by MLH-derived estimates, MLH-based chlorophyll-a concentrations on the Pacific shelf were positively biased in comparison with in-situ chlorophyll-a concentrations (Figures 9B,C). This bias could be attributable to centimeter-scale vertical heterogeneities of chlorophyll a near the sea surface, which are captured differently by satellite imagery and in-situ surveys (e.g., Harvey, 1966; Mitchell and Fuhrman, 1989). Second, as with RBD, very high MLH-derived chlorophyll-a concentrations (exceeding 123.20 mg m^–3) were determined by extrapolation of a regression line (Table 1). The threshold of 123.20 mg m^–3 corresponds to an abundance of 10^4.4 (∼25,000) cells mL^–1 (Figure 7).

At present, it is difficult to determine which of the two methods more closely matched the actual Karenia spp. abundances. We therefore chose to apply both methods and evaluate the two resultant sets of Karenia spp. abundances while taking the uncertainties of each method into account. In future work, these weaknesses could be mitigated by increasing the sample size for the regressions used in each method. This could be accomplished by incorporating in-situ measurements of Karenia spp. abundance obtained by other organizations.

There were large uncertainties in RBD- and MLH-derived abundances near the southeast coast of Hokkaido, and these methods may have underestimated in-situ Karenia spp. abundances in this region. To examine this in more detail, we extended the analysis period of Sentinel-3 data and extracted RBD- and MLH-derived Karenia spp. abundances from the 300-m grid cell that was closest to the Katsurakoi fishery harbor. On 21 September, when the second-largest in-situ Karenia spp. abundance (10^3.9 cells mL^–1) was recorded, RBD- and MLH-derived abundances were estimated to be 10³ and 10^3.9 cells mL^–1, respectively. Although the MLH-derived value was reasonable, the RBD-derived value underestimated the in-situ abundance by about one order of magnitude. Interestingly, in-situ Karenia spp. abundances from the fishery harbor differed greatly among the three monitoring sites established at that location, which were all within a few hundred meters of each other. On 23 September, when Sentinel-3 imagery and in-situ Karenia spp. abundances from all three sites were available, we obtained in-situ abundances of 10^3.4, 10^2.9, and 10^2.4 cells mL^–1 from the three harbor sites and RBD- and MLH-derived abundances of 10^2.3 and 10^2.6 cells mL^–1, respectively. Thus, satellite-derived estimates failed to reproduce the highest observed in-situ abundance on that day. This indicates that the spatial distribution of in-situ Karenia spp. abundances near the coast can change drastically, even within a few hundred meters. Hence, further work is needed to evaluate the accuracy of satellite-derived abundances near the coast by examining very-nearshore in-situ Karenia spp. abundances from additional coastal sites. Simultaneously, it will be important to use ocean-color images with a finer spatial resolution. Although Sentinel-2 ocean-color imagery (which has a horizontal resolution of 10, 20, or 60 m) might be feasible (e.g., Caballero et al., 2020), some methodological changes will be required because Sentinel 2 does not capture the 681-nm spectral band.

It is noteworthy that the satellite-derived Karenia spp. abundances estimated in this study are still imperfect; they will need to be refined to accurately monitor the presence or absence of Karenia spp. or to precisely detect signs of impending Karenia spp. blooms from weak reflectance signals. A preliminary evaluation of these use cases was conducted by using Sentinel-3 images from July 2017 to October 2021 (Supplementary Figure 3). Our tentative results indicate that MLH-derived Karenia spp. abundances overestimated the true abundances in some cases, despite the likely absence of Karenia blooms in 2017–2020, particularly during massive spring diatom blooms and autumn blooms off the southeast coast of Hokkaido (e.g., Okamoto et al., 2010; Suzuki et al., 2011; Isada et al., 2019; Kuroda et al., 2019). This overestimation was likely caused by the fact that the MLH-based estimate relates Karenia spp. abundances to chlorophyll-a concentrations. A similar pattern of overestimation in spring and autumn was also identified for RBD-derived Karenia spp. abundance, although the estimated abundances during 2017–2020 never exceeded the record-breaking abundances reached in October 2021. Hence, more sophisticated algorithms are still needed to precisely discriminate harmful Karenia blooms from non-harmful algal blooms by using satellite data, e.g., Siswanto et al. (2013), El-Habashi et al. (2016), El-Habashi et al. (2019), and Martinez-Vicente et al. (2020).

To identify Karenia blooms more precisely, it is necessary to estimate the reflectance spectral properties of Karenia spp. (in particular, K. selliformis, the properties of which have not yet been reported) and the other dominant species in our study area. It should be remembered that the RBD method that we employed was originally proposed by Amin et al. (2009) for Karenia brevis and non-K. brevis species (mainly diatoms) along the west coast of Florida (Supplementary Figure 1). In fact, our cell enumeration data (not shown) revealed that diatoms existed in mixture with Karenia spp. off the Hokkaido coast and that the abundance of diatoms relative to Karenia spp. tended to increase as Karenia spp. abundance declined. This tendency can be also seen in Figure 7; deviations of Karenia spp. abundance from the log–log regression line tended to decrease as in-situ chlorophyll-a concentrations increased (i.e., as Karenia spp. abundance increased). Moreover, note that the dominant diatom species and their spectral properties probably differ between the west coast of Florida and our study area. Hence, further work is needed to identify the dominant diatom species from long-term monitoring data along the A Line and to precisely measure the spectral properties of Karenia spp. and non-Karenia spp. associated with the dominant diatoms in our study area.

We successfully generated maps of Karenia spp. abundance and described the spatiotemporal variations in abundance from 3 to 18 October, 2021 (Figure 13). However, the possible physico-biochemical dynamics responsible for spatiotemporal changes in Karenia spp. abundance remain largely unexamined. Here, we focus on the 24-h period from 12 to 13 October, when Karenia spp. abundances averaged over the shelf abruptly increased over 2-fold (Figures 13C,D). Throughout this period, as is indicated in Figures 13C,D, our study area had few clouds along the shelf, suggesting that light conditions were favorable for algal growth.

In general, any enhancement of algal abundance can be caused by biological growth, physical advection, or a combination of the two (e.g., Richardson, 1997; Stumpf et al., 2008; Thyng et al., 2013). If the abrupt increase of Karenia app. abundance from 12 to 13 October were caused solely by algal growth, the growth rate would have been 0.74 and 0.96 day^–1 for RBD- and MLH-derived abundances, respectively. These growth rates exceed the growth rates or maximum growth rates of Karenia species reported by previous studies, which range from 0.04 to 0.41 day^–1 for K. selliformis (Medhioub et al., 2009; Mardones et al., 2020), and from 0.216 to 0.727 day^–1 for K. mikimotoi (Shen et al., 2016; Zhao et al., 2017). This discrepancy indicates that physical advection might also have contributed to the abrupt abundance increase, as suggested by the interpretation of Stumpf et al. (2008).

Physical advection can occur on a variety of spatiotemporal scales. At the shelf scale (defined here by the pink polygon in Figure 13A), horizontal onshore-offshore advection probably had a negligible effect on the abrupt increase in Karenia spp. abundance during the 24-h period, because Karenia spp. abundances were higher on the shelf than on the offshore slope. At the shelf scale, therefore, we need to primarily consider sources of vertical advection such as wind-induced upwelling, which might uplift pre-existing Karenia spp. populations from the subsurface or concentrate Karenia spp. near the sea surface through a shoaling of the mixed layer/pycnocline (e.g., Pitcher et al., 1998). However, high-resolution mesoscale model data from the Japan Meteorological Agency’s weather forecast system (Saito et al., 2006) indicate that northeasterly winds with a velocity of 3.2–6.2 m s^–1 predominated at 2 m above the sea surface over the shelf throughout 12 October. This wind regime would have favored coastal downwelling, rather than upwelling. Therefore, Karenia spp. abundance increased at the sea surface in spite of gentle downwelling-favorable winds, particularly on the western shelf.

Regarding physical advection at finer scales, the streak- and patch-like submesoscale structures we observed are likely to be among the key elements for understanding any physical-biological coupling that occurred. Streak-like filaments are also frequently identified in maps of satellite-derived chlorophyll-a concentration (e.g., Malanotte-Rizzoli et al., 2014; Shulman et al., 2015; Lehahn et al., 2017). Filamentary structures are generally interpreted as Lagrangian coherent structures, which are attributed mainly to horizontal advection (Lehahn et al., 2007; Hernández-Carrasco et al., 2018, 2020). However, this interpretation is likely to be complicated in our case, because Karenia spp. growth could have been superimposed on purely Lagrangian transport processes, and strong vertical advection induced by submesoscale variations could have contributed to the abrupt shelf-scale increase of Karenia spp. abundance (e.g., Meng et al., 2020). Further research will be needed to examine how small-scale submesoscale variations combined with algal growth contributed to the abrupt shelf-scale increase of Karenia spp. abundance, including the observed western intensification of Karenia blooms. Unfortunately, although clarification of these physical–biochemical processes is urgently needed, this is beyond the scope of our study.