Phytoplankton Bloom Changes under Extreme Geophysical Conditions in the Northern Bering Sea and the Southern Chukchi Sea

: The northern Bering Sea and the southern Chukchi Sea are undergoing rapid regional biophysical changes in connection with the recent extreme climate change in the Arctic. The ice concentration in 2018 was the lowest since observations began in the 1970s, due to the unusually warm southerly wind in winter, which continued in 2019. We analyzed the characteristics of spring phytoplankton biomass distribution under the extreme environmental conditions in 2018 and 2019. Our results show that higher phytoplankton biomass during late spring compared to the 18-year average was observed in the Bering Sea in both years. Their spatial distribution is closely related to the open water extent following winter-onset sea ice retreat in association with dramatic atmospheric conditions. However, despite a similar level of shortwave heat ﬂux, the 2019 springtime biomass in the Chukchi Sea was lower than that in 2018, and was delayed to summer. We conﬁrmed that this difference in bloom timing in the Chukchi Sea was due to changes in seawater properties, determined by a combination of northward oceanic heat ﬂux modulation by the disturbance from more extensive sea ice in winter and higher surface net shortwave heat ﬂux than usual.


Introduction
Over the past few decades, the Arctic Ocean has experienced significant warming and a dramatic decline in sea ice extent. Changes associated with diminishing sea ice include long-term thinning trends, a lengthening of the summer melt season, and a shift from primarily perennial multiyear ice to seasonal first-year ice [1]. In particular, the Pacific Arctic region, encompassing the Bering Sea and the Chukchi Sea, has experienced the most dramatic sea-ice changes in decades [2].
The northern Bering Sea and the southern Chukchi Sea, investigated in this study, are regions connecting the Pacific and Arctic Oceans through the Bering Strait and have been recognized as regions with high productivity in the world [3][4][5]. Biophysical states in these areas are influenced by the flow of heat, freshwater, and nutrients from the Pacific water [6]. Cold and saline Anadyr Water flowing along the western side of the northern shelf continually provides abundant nutrients in the northern Bering Sea and the Chukchi Sea [6], promoting high growth rates of phytoplankton during the spring season [7,8]. The spring phytoplankton bloom during the sea ice melting season is vital to annual high primary production in a given region [9]. However, due to shifts in regional atmospheric and hydrographic forcings, the areas have experienced a significant biological change [5,10]. The recent and rapid decline in sea ice across the region, associated with geophysical forcings such as temperature and wind, may significantly impact seasonal phytoplankton Table 1 shows the datasets used in this analysis. The fundamental variable to estimate chlorophyll-a concentration (Chl-a, unit: mg m −3 ), daily gridded remote sensing reflectance (Rrs(λ), unit: sr −1 ) at five wavelengths (λ = 412, 443, 488, 547, and 667 nm) at 4 km spatial resolution was obtained from Moderate Resolution Imaging Spectroradiometer (MODIS), managed by NASA's Ocean Biology Processing Group (OBPG Maryland USA) (https://oceancolor.gsfc.nasa.gov/, accessed on 7 July 2021). Although NASA OBPG also provides the daily Chl-a data, we used Rrs data to apply a regional algorithm for the appropriate calculation of Chl-a in the northern Bering and southern Chukchi Seas.
In addition, satellite-based datasets used for cloud-filling in this analysis include photosynthetically active radiation (PAR, unit: Einstein m −1 d −1 ) from MODIS. Daily sea ice concentration (SIC, unit: %) data were obtained from the National Snow and Ice Data Center (NSIDC, Colorado Boulder, USA) (https://nsidc.org/, accessed on 13 July 2021), and SIC data were derived from the NASA Team algorithm, which has a resolution of 25 km on a polar stereographic grid [14]. For daily sea surface temperature (SST, unit: • C), we used the Optimal Interpolation Sea Surface Temperature (OISST) [15] provided by the National Oceanic and Atmospheric Administration (NOAA, Washington, USA) (https://www.ncdc.noaa.gov/oisst, accessed on 13 July 2021). The OISST is recorded daily at a 25 km resolution. Furthermore, for the atmospheric parameters such as 10-m winds (zonal (U10) and meridional (V10) components, units: m s −1 ), 2-m atmospheric temperature (T2M, unit: • C), and net surface shortwave flux (unit: W m −2 ), the European Centre for Medium-Range Weather Forecasts Reanalysis v5 (ERA5) hourly datasets with a 31 km spatial resolution were used (https://www.ecmwf.int/, accessed on 15-July-2021). In addition, ocean floor depth data were provided by the general bathymetric chart of the ocean (GEBCO), and have a spatial resolution of approximately 30 arc-seconds (https://www.gebco.net, accessed on 23 March 2021). All data were remapped into the resolution on the MODIS data using nearest-neighbor interpolation. In addition, for the mixed layer depth estimation, the model data has been taken from the Copernicus Marine Environment Monitoring Service (CMEMS, Europe) global weekly coupled ocean-atmospheric data assimilation, using 3D-var, which utilized GLO-CPL for 3 d ocean forecasts at 0.25 degree, and the Met Office Unified Model coupled to an hourly NEMO v 3.4 [16]. Ocean observations include satellite SST data, in-situ SSTs from moored buoys, drifting buoys and ships, sea level anomaly observations from CMEMS, subsurface temperature, and salinity profiles (including Argo, moored buoys, marine mammals, and gliders), and satellite-based sea ice concentration.

Chl-a Estimation
The abundance of phytoplankton can convert inorganic carbon dioxide into organic carbon through photosynthesis in the upper layers of the ocean, and plays a considerable role in the global carbon cycle through its function as a biological pump. Therefore, understanding phytoplankton is of paramount importance in climate change research [17]. Chl-a is a proxy to estimate phytoplankton biomass. Consequently, it is widely used in temporal and spatial phytoplankton dynamic analysis as it can be frequently measured from ocean color at a large scale. Therefore, the accuracy of Chl-a must be guaranteed. In this paper, a three-step process was performed for a more accurate Chl-a estimation. First, to resolve the uncertainty issue of satellite-derived blue wavelength products (e.g., Rrs(412), Rrs(443)), the blue-band estimation algorithm suggested by [18] was used. This is a method to solve the problem that reliable Chl-a values cannot be produced because the blue band used for Chl-a calculation tends to have large uncertainties in areas such as coastal and inland waters and due to sensor degradation. In brief, blue bands are estimated using the values of the three wavelengths of 48×, 55×, and 67× nm and 1763 observation-based spectra references. Second, the quality assurance system proposed by [19] was applied for selecting Rrs(λ) data that is superior in quality. They established Rrs spectra references for 23 optical water types through observation and developed a method to calculate the final score through similarity between the remote sensing Rrs spectrum and the reference. We applied this system to five bands of MODIS data, and only data with a total score of 0.8 or higher were used. Lastly, the regional Chl-a algorithm developed in the Arctic Ocean by [20] was used in this study. We considered this study area as the "inflow" area they classified and applied this algorithm as follows.

Filling Gaps on Chl-a Data
We used a machine learning technique to fill the cloud-induced gaps in ocean color, which allows for a more accurate interpretation of phytoplankton biomass in time and space. At high latitudes, ocean color products such as Chl-a contain massive missing data due to various causes, such as the presence of clouds and sea ice [21]. If the spatial and temporal mean is constructed with missing data, the accuracy may be significantly affected [22]. The machine learning model used in this study is the random forest (RF) [23]. The overall framework is similar to previous works [21,24]. However, there were several different preprocessing steps for Chl-a images: (i) a 3 × 3 median filter was applied to Chl-a data, (ii) Chl-a images with only 0.1% pixels were excluded, and (iii) the pixels that were not observed 15 days before and after were also removed for continuity. There were a total of nine input variables used for model training: SST, T2M, PAR, U10, V10, water depth, longitude, latitude, and days of the year. Unlike the previous versions, we excluded the climatology of Chl-a from input variables because the model performance for this study region deteriorated when including it.
As a result, 17,192,761 pixels were used for model training, and 5,730,921 pixels were used for validation. Furthermore, the critical parameters determining ultimate model performance, number of trees, and number of features were set to 60 and 3, consistent with those in [21]. As such, the R-squared (R 2 ), root mean square error (RMSE), and mean absolute error (MAE) for the training set were 0.99, 0.02 mg m -3 , and 0.01, respectively ( Figure 1). Those for the validation set were 0.99, 0.05 mg m -3 , and 0.03, respectively. Finally, Chl-a images reconstructed during winter (November to January), when few standard Chl-a data are used for training, were excluded from the analysis.

Filling Gaps on Chl-a Data
We used a machine learning technique to fill the cloud-induced gaps in ocean color, which allows for a more accurate interpretation of phytoplankton biomass in time and space. At high latitudes, ocean color products such as Chl-a contain massive missing data due to various causes, such as the presence of clouds and sea ice [21]. If the spatial and temporal mean is constructed with missing data, the accuracy may be significantly affected [22]. The machine learning model used in this study is the random forest (RF) [23]. The overall framework is similar to previous works [21,24]. However, there were several different preprocessing steps for Chl-a images: (i) a 3 × 3 median filter was applied to Chla data, (ii) Chl-a images with only 0.1% pixels were excluded, and (iii) the pixels that were not observed 15 days before and after were also removed for continuity. There were a total of nine input variables used for model training: SST, T2M, PAR, U10, V10, water depth, longitude, latitude, and days of the year. Unlike the previous versions, we excluded the climatology of Chl-a from input variables because the model performance for this study region deteriorated when including it.
As a result, 17,192,761 pixels were used for model training, and 5,730,921 pixels were used for validation. Furthermore, the critical parameters determining ultimate model performance, number of trees, and number of features were set to 60 and 3, consistent with those in [21]. As such, the R-squared (R 2 ), root mean square error (RMSE), and mean absolute error (MAE) for the training set were 0.99, 0.02 mg m -3 , and 0.01, respectively (Figure 1). Those for the validation set were 0.99, 0.05 mg m -3 , and 0.03, respectively. Finally, Chl-a images reconstructed during winter (November to January), when few standard Chl-a data are used for training, were excluded from the analysis.

Monthly Variations in Satellite-Derived Chl-a
The satellite-derived monthly climatological (from January 2003 to December 2020) Chl-a images in the northern Bering Sea and the southern Chukchi Sea are shown in Figure 2a. In January, sea ice and clouds cover most of the area, making it entirely impossible to observe Chl-a. Simultaneously with the seasonal retreat of the sea ice edge, Chl-a values are provided in confined regions. The spring bloom in the northern Bering Sea begins in March and peaks in May, particularly in the Gulf of Anadyr and Bering Strait. At that time, the southern Chukchi Sea also exhibited high biomass. These blooms start to terminate rapidly from June, and Chl-a increases weakly around the Bering Strait in the late summer season (August). In October, due to recurrent sea ice advance, the Chl-a data only cover the northern Bering Sea, and from November, the area is again completely covered with clouds and sea ice. The time series associated with such change in Chl-a is presented in Figure 2b. The number of Chl-a data available through the reconstruction process in this study has increased by six times; in particular, 6.67 times in the Bering Sea and 5.28 times in the Chukchi Sea. Due to this increased Chl-a data, the time series is generally more stable than the fluctuant standard satellite Chl-a data.

Monthly Variations in Satellite-Derived Chl-a
The satellite-derived monthly climatological (from January 2003 to December 2020) Chl-a images in the northern Bering Sea and the southern Chukchi Sea are shown in Figure 2a. In January, sea ice and clouds cover most of the area, making it entirely impossible to observe Chl-a. Simultaneously with the seasonal retreat of the sea ice edge, Chl-a values are provided in confined regions. The spring bloom in the northern Bering Sea begins in March and peaks in May, particularly in the Gulf of Anadyr and Bering Strait. At that time, the southern Chukchi Sea also exhibited high biomass. These blooms start to terminate rapidly from June, and Chl-a increases weakly around the Bering Strait in the late summer season (August). In October, due to recurrent sea ice advance, the Chl-a data only cover the northern Bering Sea, and from November, the area is again completely covered with clouds and sea ice. The time series associated with such change in Chl-a is presented in Figure 2b. The number of Chl-a data available through the reconstruction process in this study has increased by six times; in particular, 6.67 times in the Bering Sea and 5.28 times in the Chukchi Sea. Due to this increased Chl-a data, the time series is generally more stable than the fluctuant standard satellite Chl-a data.  Figure 3a shows the annual variation in Chl-a and SIC. The Chl-a variation is mainly contributed to by the spring bloom in April-June. It has a distinct negative relationship  Figure 3a shows the annual variation in Chl-a and SIC. The Chl-a variation is mainly contributed to by the spring bloom in April-June. It has a distinct negative relationship with SIC (R = −0.62, p < 0.01 from student t-test at a significance level of 95%). This high correlation was not calculated with unreconstructed standard Chl-a data (R = −0.36, p = 0.13). The annual Chl-a shows a decreasing trend during the early decade since the start of MODIS observation. The lowest biomass recorded was in 2013, with 0.96 mg m −3 across the entire period. After 2013, the biomass gradually increases, reaching the greatest Chl-a with 1.  Figure 3b shows the SIC phenology in 2017-2019, with the climatological phenology during 2003-2020. SIC in 2018 is abnormally and continuously lower than normal from winter (January) to the beginning of summer (July). In addition, the seasonal ice retreat in 2019 also starts two months earlier, unlike the ordinary retreat beginning in April. In 2017, although the winter ice advance was delayed, it quickly recovered to a normal level in April.

Annual Variation in Chl-a Linked with Sea Ice Change
with SIC (R = −0.62, p < 0.01 from student t-test at a significance level of 95%). This high correlation was not calculated with unreconstructed standard Chl-a data (R = −0.36, p = 0.13). The annual Chl-a shows a decreasing trend during the early decade since the start of MODIS observation. The lowest biomass recorded was in 2013, with 0.96 mg m −3 across the entire period. After 2013, the biomass gradually increases, reaching the greatest Chl-a with 1.22 mg m −3 in 2017, followed by 1.21 mg m −3 in 2018 and 1.16 mg m −3 in 2019, and then in 2020 dramatically falling back to 2003 levels. In association with the annual Chl-a variation, sea ice reaches its maximum in 2012 (47.9%), limited to this study period, and then declines sharply, reaching its lowest point in 2019 (26.3%). Even so, 2017 and 2018 show significant ice loss compared to normal years as much as that of 2019 (30.0% in 2017 and 26.6% in 2018). Figure 3b shows the SIC phenology in 2017-2019, with the climatological phenology during 2003-2020. SIC in 2018 is abnormally and continuously lower than normal from winter (January) to the beginning of summer (July). In addition, the seasonal ice retreat in 2019 also starts two months earlier, unlike the ordinary retreat beginning in April. In 2017, although the winter ice advance was delayed, it quickly recovered to a normal level in April.
In connection with annual ice changes, Chl-a had a weak negative trend until almost 2012, thereafter gradually increasing and experiencing a significant drop again in 2020 along with the 2020 turnover of sea ice (Figure 3c). This pattern appears throughout the northern Bering Sea and eastern Chukchi Sea (Figure 3d), and is especially evident adjacent to the Gulf of Anadyr (R > 0.9, near 175° W and 63° N). In the southern Chukchi Sea, such a pattern is evident on the eastern side (R > 0.7).  In connection with annual ice changes, Chl-a had a weak negative trend until almost 2012, thereafter gradually increasing and experiencing a significant drop again in 2020 along with the 2020 turnover of sea ice (Figure 3c). This pattern appears throughout the northern Bering Sea and eastern Chukchi Sea (Figure 3d), and is especially evident adjacent to the Gulf of Anadyr (R > 0.9, near 175 • W and 63 • N). In the southern Chukchi Sea, such a pattern is evident on the eastern side (R > 0.7).

Response of Chl-a to Episodic Climatic Patterns: 2018-2019
In 2018 and 2019, the ice retreat occurred earlier than ever (Figure 3b), SIC remained below average levels until summer, and the associated abundance of phytoplankton was high ( Figure 3a). In order to examine the causes of the extreme changes in SIC and Chl-a during these two years, the spatial distributions of the anomalies of SIC, T2M, and wind fields from winter to spring bloom season (January to June) in both years were investigated (Figure 4), including climatological monthly patterns. Typical prevailing wind patterns over the study region during January-February are dominated by northeasterly, with mostly very low atmospheric temperatures (Figure 4a). Therefore, the area is mainly covered by sea ice, with the sea ice edge located at the southwestern end of the domain. During MA, the predominant northeasterly winds and low T2M were maintained, and the sea ice edge slightly advanced ( Figure 4b). Thereafter, seasonal changes, such as elevated temperatures and weakened northerly winds, removed most of the sea ice covering the Bering Sea ( Figure 4c). In January-February 2018 (Figure 4d), due to a strong southwesterly wind anomaly and warm T2M anomaly in the Bering Sea, a significant negative SIC anomaly (<−60%) was formed, consistent with the spatial distribution of a high air temperature anomaly (near    Figure 5 shows the distribution of Chl-a during the early spring (March to April), late spring (May to June), and summer (July to August) seasons, including the positions of the marginal sea ice edges. In March and April 2018, the ice edge receded far from the climatological ice edge position, and most of the northern Bering Sea was exposed to the atmosphere ( Figure 5a). Nevertheless, phytoplankton biomass was at a level similar to that of normal years. A weak positive Chl-a anomaly patch existed only in part of the eastern Bering shelf. Instead, due to the substantial displacement of the ice edge, a strong negative anomaly was observed for both years in the areas associated with the bloom that had previously appeared in the marginal sea ice zone (near climatological ice edge). The spatial distributions of the Chl-a anomaly and open water in March-April 2019 are not substantially different from those in 2018 ( Figure 5b). However, the eastern Bering shelf exhibited the opposite Chl-a anomaly pattern (Figure 5c). The rapid ice retreat and high positive Chl-a anomaly were observed in May-June as a result of weakened wind forcing, warmer temperature (Figure 4c), and a seasonal increase in solar radiance [25]. In 2018, a strong  Figure 5 shows the distribution of Chl-a during the early spring (March to April), late spring (May to June), and summer (July to August) seasons, including the positions of the marginal sea ice edges. In March and April 2018, the ice edge receded far from the climatological ice edge position, and most of the northern Bering Sea was exposed to the atmosphere ( Figure 5a). Nevertheless, phytoplankton biomass was at a level similar to that of normal years. A weak positive Chl-a anomaly patch existed only in part of the eastern Bering shelf. Instead, due to the substantial displacement of the ice edge, a strong negative anomaly was observed for both years in the areas associated with the bloom that had previously appeared in the marginal sea ice zone (near climatological ice edge). The spatial distributions of the Chl-a anomaly and open water in March-April 2019 are not substantially different from those in 2018 (Figure 5b). However, the eastern Bering shelf exhibited the opposite Chl-a anomaly pattern (Figure 5c). The rapid ice retreat and high positive Chl-a anomaly were observed in May-June as a result of weakened wind forcing, warmer temperature (Figure 4c), and a seasonal increase in solar radiance [25]. In 2018, a strong positive Chl-a anomaly was primarily and widely observed across most of the Bering Sea (Figure 5d correlate strongly with open water extent. The phytoplankton bloom during summertime (July-August) 2018 was similar to normal, except for a weak negative Chl-a anomaly in the Bering Strait (Figure 5g). Instead, a strong positive Chl-a anomaly was observed in the Bering Strait in 2019 (Figure 5h), and the Chukchi Sea also exhibited higher Chl-a than average. This delayed summertime phytoplankton bloom in the Chukchi Sea in 2019 was markedly intense (~0.14 mg m −3 ) throughout the entire period ( Figure 6).   [1,11]. However, despite the similar open water formation in the Chukchi Sea, we found a significant difference in springtime phytoplankton bloom between the two years, with the Chukchi bloom delayed to summer in 2019. Therefore, we theorized that other factors predominantly drive the phytoplankton bloom in the Chukchi Sea.
One plausible factor regulating bloom formation in the eastern Chukchi Sea could be seawater temperature [28].   [1,11]. However, despite the similar open water formation in the Chukchi Sea, we found a significant difference in springtime phytoplankton bloom between the two years, with the Chukchi bloom delayed to summer in 2019. Therefore, we theorized that other factors predominantly drive the phytoplankton bloom in the Chukchi Sea.
One plausible factor regulating bloom formation in the eastern Chukchi Sea could be seawater temperature [28]. Figure 7 shows the distribution of SST anomalies from March to June of both years, including SIC and net shortwave heat flux anomalies. The spatial distributions of ice melting and positive SST anomaly are closely related to the positive (downward) heat flux anomaly in 2018. The Chukchi Sea exhibited a positive SST anomaly from April, especially along the Alaskan coasts (Figure 7b). It then intensified significantly, associated with a positive shortwave heat flux anomaly from May to June (Figure 7c,d). In July, there was no significant difference from normal years (not shown). In 2019, despite the magnitude and distribution of shortwave heat flux similar to those of 2018, the high SST anomaly patch was limited to the Bering Sea until April, with a slightly negative SST anomaly to the north (Figure 7e,f). It was not until May that a narrow positive SST anomaly formed along the Alaskan coast, barely reaching the Chukchi Sea. Even at this point, these narrow, north-south positive SST anomaly patches were surrounded by the negative anomaly (Figure 7g). However, in June, a rather explosively strong positive SST anomaly appeared over most of the Chukchi Sea and the Bering Sea (Figure 7h).

Factor Impeding SST Increase in the Chukchi Sea in 2019
We assumed that the distributional patterns of the SST anomaly in the Chukchi Sea ultimately determined biomass and tried to find the cause of the different SST anomaly distribution between both years. We illustrated the basin-scale distribution of the monthly SST anomaly in Figure 8 with sea ice motion and ice edge location. For ice motion data, Polar Pathfinder Daily 25 km EASE-Grid Sea Ice Motion Vectors, Version 4 was obtained from the NSIDC [29]. Unlike 2018, which recorded a minimum ice extent in winter, 2019 recorded a winter SIC that exceeded the average for 18 years from 2003 to 2020 (Figure 3b). Due to the quantitative level of winter sea ice, there were low SST anomaly distributions in the entire northern Bering Sea and along the western coasts from January to February (Figure 8a,b). The sea ice during this period appeared to be retreating in January, but was moving southward again in February, hindering high water temperatures from advancing northward. In March, SST anomalies began to increase slowly outside the southern part of the study area (i.e., the southeastern part of the Bering continental shelf) compared to normal years. In addition, although it seemed to push the sea ice away with a slight northward flow, there was still more sea ice in the northern Bering Sea than in other years, resulting in a lower SST anomaly. As warming intensified in the southeastern part of the region from April (Figure 8d), it gradually expanded to the north and entered the Chukchi Sea along the Alaska coastal regions in May (Figure 8e), consistent with the distribution of the shortwave heat flux (Figure 7g). At the same time, whether it was the result of sea ice itself or cold water from melting, SST anomalies were pushed into the Chukchi Sea due to the northward movement of the anomalous warm pool and likely to reduce the ambient temperature.

Factor Impeding SST Increase in the Chukchi Sea in 2019
We assumed that the distributional patterns of the SST anomaly in the Chukchi Sea ultimately determined biomass and tried to find the cause of the different SST anomaly distribution between both years. We illustrated the basin-scale distribution of the monthly SST anomaly in Figure 8 with sea ice motion and ice edge location. For ice motion data, Polar Pathfinder Daily 25 km EASE-Grid Sea Ice Motion Vectors, Version 4 was obtained from the NSIDC [29]. Unlike 2018, which recorded a minimum ice extent in winter, 2019 recorded a winter SIC that exceeded the average for 18 years from 2003 to 2020 ( Figure  3b). Due to the quantitative level of winter sea ice, there were low SST anomaly distributions in the entire northern Bering Sea and along the western coasts from January to February (Figure 8a,b). The sea ice during this period appeared to be retreating in January, but was moving southward again in February, hindering high water temperatures from advancing northward. In March, SST anomalies began to increase slowly outside the southern part of the study area (i.e., the southeastern part of the Bering continental shelf) compared to normal years. In addition, although it seemed to push the sea ice away with a slight northward flow, there was still more sea ice in the northern Bering Sea than in other years, resulting in a lower SST anomaly. As warming intensified in the southeastern part of the region from April (Figure 8d), it gradually expanded to the north and entered the Chukchi Sea along the Alaska coastal regions in May (Figure 8e), consistent with the distribution of the shortwave heat flux (Figure 7g). At the same time, whether it was the The propagation patterns on Hovmöller diagrams of zonally averaged SST anomalies in 2018 and 2019 are shown in Figure 9. In 2018, a positive SST anomaly in the Bering Sea was propagated to the Chukchi Sea for a month from April, and the heat in the Bering Sea transferred again in May. In June, since most of the area was exposed to the open ocean (Figures 4 and 7), a strong positive SST anomaly occurred overall due to the combination of an intense positive shortwave heat flux (Figure 7) and the seasonal temperature increase. On the other hand, in March 2019, there was a low SST anomaly due to the presence of sea ice south of the Bering Strait, which had a strong temperature gradient with the south. The relatively warm SST anomaly south of the Bering Sea began to move northward at a slower pace from April and gradually accelerated until June. At the same time, the negative SST anomaly was progressively pushed to the north. After June, a strong positive SST anomaly occurred throughout the study area, similar to 2018, though it was more robust than in 2018.
The spatial difference in solar heating in the study area is not significant (refer to Figure 7). The seasonal solar heating due to shortwave heat flux is most vigorous in mid-June and does not differ significantly year to year compared to the northward oceanic heat flux [30]. In fact, in June of both years, the shortwave heat flux ranged from 151.6 to 283 W m −2 , with an average of 213.1 W m −2 , and the average in 2019 was 213.9 W m −2 (158.9-262.3 W m −2 ), which is highly similar between the two years ( Figure 9c).  To examine the relationship between the extent of winter sea ice in the northern Bering Sea and monthly SST in the Chukchi Sea, we calculated the correlation between the year-to-year SIC in the Bering Sea from January to March and the SST in the Chukchi Sea during the months after that ( Table 2). Our results show that winter SIC in the Bering Sea significantly affects the SST in the Chukchi Sea in June. From January to March, the Bering The spatial difference in solar heating in the study area is not significant (refer to Figure 7). The seasonal solar heating due to shortwave heat flux is most vigorous in mid-June and does not differ significantly year to year compared to the northward oceanic heat flux [30]. In fact, in June of both years, the shortwave heat flux ranged from 151.6 to 283 W m −2 , with an average of 213.1 W m −2 , and the average in 2019 was 213.9 W m −2 (158.9-262.3 W m −2 ), which is highly similar between the two years ( Figure 9c).
To examine the relationship between the extent of winter sea ice in the northern Bering Sea and monthly SST in the Chukchi Sea, we calculated the correlation between the year-to-year SIC in the Bering Sea from January to March and the SST in the Chukchi Sea during the months after that ( Table 2). Our results show that winter SIC in the Bering Sea significantly affects the SST in the Chukchi Sea in June. From January to March, the Bering Sea SIC has a very strong negative relationship with the Chukchi Sea SST in June, −0.73, −0.86, and −0.74, respectively. SIC is generally linked with Chukchi Sea SST from April to July, although it is most strongly associated with Chukchi Sea SST in June. Another distinct feature is that wintertime Bering Sea SIC has a strong correlation with SST in autumn, in particular, the Bering Sea SIC in February (R = −0.85, p < 0.01) and March (R = −0.83, p < 0.01).

Chl-a Distributional Features during Extreme Atmospheric Events
During the MODIS era (2003-2020), Chl-a in the northern Bering Sea and southern Chukchi Sea (especially eastern parts) has undergone significant changes, alongside changes in sea ice (Figure 3). Chl-a in this region has risen since 2010, with a high abundance between 2017-2019 when significant sea-ice loss occurred [31,32]. In particular, the ice extent in 2018 was the lowest since observations in the 1970s [9]. Consequently, over three years, such extreme ice-cover changes have impacted oceanographic and biological conditions in the northern Bering Sea and the Chukchi Sea, including impacts on lower and upper trophic levels [33,34]. Annual Chl-a during this period has been higher than ever, which has been associated with annual SIC reduction (R = −0.62, p < 0.01). This relationship implies that increased open water extent has relaxed light restrictions, promoting phytoplankton growth, and increasing open ocean primary production [11].
We identified the atmospheric characteristics of 2018 and 2019, with unusual sea ice loss from winter to spring (Figure 4). We excluded 2017 from our analysis because, although the sea ice advance started late due to warm southerly winds in winter, it rapidly recovered to the average level of ice extent. On the other hand, in 2018 and 2019, the sea ice extent continued to be far below the average from the winter season until the season when the sea ice completely melted. Usually, by March, when sea ice reaches a maximum in the Bering Sea [1,35], the region is typically dominated by northerly winds (Figure 4a,b) [9]. However, during those two years, the seasonal advance of sea ice was hampered by the exceptionally warm air temperature and the strengthening of southerly winds (Figure 4d,g). These unprecedented warm southerly winds lasted until early spring (March and April) and not only resulted in the acceleration of the Bering Sea ice retreat (Figure 4e,h) but also affected the open water extent of the Chukchi Sea in summer (Figure 4f,i). In conjunction with such anomalous atmospheric forcing for both years, the Aleutian Low in winter can be considered a critical factor in the study region [36]. The position and intensity of this atmospheric pressure system have a significant correlation with regional circulation and, in turn, sea ice distribution in the Bering Sea [37]. In particular, when this low moves west, it pushes warm Pacific air toward the eastern Bering Sea [37,38]. In fact, in 2018, the core of Aleutian Low shifted about 20 degrees westward from its climatological longitudinal position (180 • ) [9].
Although the atmospheric properties of the two years under extreme atmospheric conditions are not significantly different, the spatiotemporal distributions of Chl-a are markedly different ( Figure 5). Blooms, which generally occur in the marginal ice zone in the Bering Sea in March and April, did not appear due to the unusual sea ice retreat in both years, resulting in a strongly negative Chl-a anomaly across the region (Figure 5a,b). In May-June 2018, an abnormal increase in phytoplankton biomass was observed over most of the Bering Sea, with exceptionally high biomass in the western part of the Bering Sea shelf, on the typical route of the Anadyr Current (centered on the Gulf of Anadyr) (Figure 5d). The Anadyr Current flows clockwise along the Gulf of Anadyr, supplying highly nutritious Anadyr Water to the western Bering Sea, and, in turn, passes the Anadyr Strait, the Chirikov Basin, the Bering Strait, and into the Chukchi Sea [38,39]. The magnitude of the spring bloom that generally occurs in April-May (sometimes even in June [4]), along with high nutrient supply and a seasonal increase in solar radiation in the Adryna Bay, is about 10-13 mg m −3 and is much larger than that of the eastern Bering shelf area [40]. In addition, in 2018, a strongly positive Chl-a anomaly was observed in the vicinity of St. Lawrence Island, probably being influenced by episodic Anadyr Water [38] or related to a polynya formed to the south [4]. However, in 2019, despite a strong negative Chl-a anomaly which formed along the Russian coast, a strong positive Chl-a anomaly patch existed in the southeast of Anadyr Bay (Figure 5e). The spatial characteristics of the anomalous late spring phytoplankton blooms in the Bering Sea appear to be closely related to the distribution of sea ice modulated by atmospheric forcings. In particular, the distribution of Chl-a anomalies in May-June (Figure 5d,e) is closely related to the spatial distribution of open water formed by anomalous ice retreat from winter to spring (Figure 4e,h). If sea ice retreats earlier, ice-edge bloom can initiate due to solar radiation limitations [41]. However, Since May-June is a period of sufficiently high seasonal shortwave radiation, the positive Chl-a anomaly in the Bering Sea was caused by an open ocean bloom. Therefore, we argue that the Bering Sea for two years showed an increased phytoplankton biomass, especially in areas with longer growing seasons, in association with the early timing of ice retreat [11,42].
During May-June 2018, a positive Chl-a anomaly in the Chukchi Sea was observed along the Alaskan coasts, however the biomass in 2019 was at the level of a typical year. If it is assumed that this area is greatly affected by the longer growing season and the larger open water extent by early ice retreat like the Bering Sea, the characteristics of atmospheric parameters and ice edge position between the two years should be significantly different. However, our results did not show such a fact. Although ice distribution in the Chukchi Sea may have been involved in the biomass to some extent, it is likely that other factors also had a significant impact. In the late spring of 2019, a positive Chl-a anomaly was not be observed in the Chukchi Sea. However, a relatively high phytoplankton abundance was observed in the eastern part of the Chukchi Sea in summer (July and August) 2019, including the Bering Strait. We will discuss this in the following section.

Impact of Winter Bering Sea Ice on Seawater Temperature in the Chukchi Sea
One possible cause of the delayed summertime bloom in 2019 is seawater temperature. Ocean temperature can be related to the growth of phytoplankton in the upper ocean directly by physiological factors or indirectly by stratification facilitating light availability. The primary sources of influence on the temperature of the Chukchi Sea are heat flux flowing from the Pacific Ocean or solar radiation. In March, the Bering Sea SST in both years rose due to relatively lower ice extent and higher net shortwave heat flux than those in average years (Figure 7a,e). However, in April, a robust downward shortwave heat flux formed along the Alaskan coast due to the warm Alaskan Coastal Current [43]. Although the degree and distribution of shortwave heat flux were similar in both years, the temperature in 2019 was only at the average level (Figure 7b,f). Then, due to continued solar heating, the SST in the Chukchi Sea increased from May to June in 2018 (Figure 7c,d).
On the other hand, in 2019, SST increased slightly in May ( Figure 7g) and more explosively in June (Figure 7h). This late increase in SST might be attributed to the degree of northward heat transport from the Pacific water and is likely to be closely related to the delayed and weakened bloom in the Chukchi Sea. The northward heat transport of Pacific water through the Bering Strait has more energy than the solar heating that occurs in the Chukchi Sea. Furthermore, the interannual variation in solar heating is not considerable compared to that of the northward heat flux [30]. Thus, it appears that by April 2019, the northward heat flux was hampered by relatively cold water or ice forming along the pathway of the Anadyr Coastal Current (refer to Figures 8 and 9). In May, heat transport to the Chukchi Sea began slowly, and in June, the heat was rapidly amplified when combined with solar heating in the Chukchi Sea [39,43]. In 2019, compared to 2018, it seems that a large amount of sea ice formed in the Bering Sea in winter (refer to Figure 3b), preventing the northward movement of the unusual warm pool that occurred on the eastern Bering continental shelf. In fact, SIC in the Bering Sea from January to March was closely and negatively related to SST in the Chukchi Sea throughout spring and summer (especially in June) ( Table 2); this relationship has been maintained to some extent since then and again has a close relationship with the sea surface temperature of the Chukchi Sea in October when the sea ice retreated the most.

Role of Seawater Temperature Regulating Phytoplankton Biomass
Seawater temperature could impact phytoplankton biomass both directly and indirectly. Direct pathways include physiological effects, with higher metabolic rates appearing at higher temperatures to promote the growth of phytoplankton [44]. It has long been recognized that phytoplankton grows faster in warm water within a temperature range detrimental to macromolecules [44][45][46]. The previous physical-biological simulation result showed that sea ice is vital for the seasonal timing of the phytoplankton bloom in the Bering Sea shelf and that the bloom is controlled by the temperature of the surrounding seawater in the Chukchi Sea [28], consistent with our results. The results also suggest that water temperature can directly promote the growth of phytoplankton in the Chukchi Sea.
With the direct water temperature effect, seawater temperature can modulate the light availability of phytoplankton through thermal-based ocean mixed layers. Figure 10 shows the seasonal changes in the mixed layer depth from May, when the ice cover starts to be removed, to August for both years in the Chukchi Sea (160-170 • W, 66-72 • N), where the positive Chl-a anomaly occurred. Both years show marked seasonal variation, with a minimum (~10 m) in late June [47]. The 2018 MLD started getting shallower through mid-May, associated with increased SST (Figure 7). Up until this point, the MLD in 2018 was 1-2 m shallower than the MLD in 2019. Since then, MLD in 2019 has rapidly shallowed until about 10 m in mid-June, which appears to be associated with a significant heat supply from the south (Figures 8 and 9). Due to a month later and stronger stratification compared to 2018, the bloom in the Chukchi Sea in 2019 seems to have been delayed into the summer season compared to 2018. Although there is likely to be some relationship between the bloom timing in the Chukchi sea and the change in MLD, the year-round difference in MLD of several meters may not be significant. However, given that the nutrient profiles observed from July to September in the Chukchi Sea showed the depletion to about 20 m [48], and that the average depth is about 50 m, this degree of difference can be considered significant. considered significant.
As a result, the ocean temperature in the Chukchi Sea is simultaneously affected by northward oceanic heat transport according to the level of winter sea ice (especially in February) and solar heating by shortwave heat flux, and directly or indirectly affects the timing of phytoplankton bloom.  As a result, the ocean temperature in the Chukchi Sea is simultaneously affected by northward oceanic heat transport according to the level of winter sea ice (especially in February) and solar heating by shortwave heat flux, and directly or indirectly affects the timing of phytoplankton bloom.

Conclusions
We compared the characteristics of the distribution of Chl-a in the northern Bering Sea and the southern Chukchi Sea in 2018 and 2019, which experienced extreme atmospheric and oceanographic changes and recorded high annual phytoplankton biomass. As a result, the northern Bering Sea exhibited relatively high biomass during late spring in both years, associated with the exceptionally early and broad early spring open water extents caused by the warm southerly winds, blowing continuously from winter to early spring. On the other hand, in the Chukchi Sea, we found that the seawater temperature controlled by the northward oceanic heat transport and shortwave heat flux is a significant factor for phytoplankton abundance rather than the open water extent. According to the modulation of such a factor, the delayed summer bloom in the Chukchi Sea occurred in 2019, unlike the late spring bloom in 2018. However, these changes are not limited to just these two years. Arctic temperatures in the Pacific have warmed more than twice the global average in recent decades, associated with a robust Arctic amplification. Anomalous atmospheric events have become more frequent and more substantial in recent years [31], and model results report that this phenomenon will increase over the next few decades [31,49]. Furthermore, the extended ice-free season in the Bering Sea and the Chukchi Sea caused by climate change might be attributed to more robust flows and, accordingly, more transport of heat and nutrients from the Pacific waters into the Arctic. This could eventually lead to high production, carbon sequestration, and a negative feedback on global warming in the regions [10]. Funding: This study was financially supported by the the modulation of such a factor, the delayed summer bloom in the Chukchi Sea occurred in 2019, unlike the late spring bloom in 2018. However, these changes are not limited to just these two years. Arctic temperatures in the Pacific have warmed more than twice the global average in recent decades, associated with a robust Arctic amplification. Anomalous atmospheric events have become more frequent and more substantial in recent years [31], and model results report that this phenomenon will increase over the next few decades [31,49]. Furthermore, the extended ice-free season in the Bering Sea and the Chukchi Sea caused by climate change might be attributed to more robust flows and, accordingly, more transport of heat and nutrients from the Pacific waters into the Arctic. This could eventually lead to high production, carbon sequestration, and a negative feedback on global warming in the regions [10]. the modulation of such a factor, the delayed summer bloom in the Chukchi Sea occurred in 2019, unlike the late spring bloom in 2018. However, these changes are not limited to just these two years. Arctic temperatures in the Pacific have warmed more than twice the global average in recent decades, associated with a robust Arctic amplification. Anomalous atmospheric events have become more frequent and more substantial in recent years [31], and model results report that this phenomenon will increase over the next few decades [31,49]. Furthermore, the extended ice-free season in the Bering Sea and the Chukchi Sea caused by climate change might be attributed to more robust flows and, accordingly, more transport of heat and nutrients from the Pacific waters into the Arctic. This could eventually lead to high production, carbon sequestration, and a negative feedback on global warming in the regions [10].