Climate control of sea-ice edge phytoplankton blooms in the Hudson Bay system

The Arctic Ocean and its surrounding seas are facing the most pronounced climatic changes on the Earth. Several regional positive feedback processes amplify the warming of the Arctic (Overland et al., 2004; Serreze et al., 2009), such as the sea-ice albedo (Screen and Simmonds, 2010), decline of sea ice in all seasons (Stroeve and Notz, 2018), and the lapse rate and the thermal radiative balance (Pithan and Mauritsen, 2014). Those changes also increase the frequency of extreme climatic events in the northern latitudes (Pithan and Mauritsen, 2014). As a consequence, the timing of Arctic-wide sea-ice melt season has advanced at a rate of 5 days per decade since 1979 (Stroeve et al., 2011, 2014). Arctic sea-ice volume has declined at a rate of –513 km³y^–1 and –287 km³y^–1 during winter and fall season, respectively, between 2002 and 2018 (Kwok, 2018). Climatic changes have been even more accentuated in the Hudson Bay System (HBS), a complex ecosystem embracing Hudson Bay, James Bay, Foxe Basin, and Hudson Strait, which forms the world’s most extensive inland sea (1.24 × 10⁶ km²). For example, Stroeve and Notz (2018) showed that Hudson Bay sea-ice cover in September has decreased at the rate of –1,046 km² y^–1 between 1979 and 2018, a loss of –93.6% relative to the average of 45.2% for the whole Arctic. As reported by Hochheim et al. (2011), sea-ice concentration (SIC) losses range from –15.1% to –20.4% per decade since 1980 in the northwest and southwest sectors of Hudson Bay (HB), respectively. As a consequence, the duration of the open water season increased by 12 days per decade between 1980 and 2005, which is almost twice the rate observed in the Arctic Ocean (i.e., 6.4 days per decade; (Markus et al., 2009). Finally, the mean annual sea surface temperature trend is about 3.7 °C since the postindustrial era, which is much larger than the trend observed over the Arctic Ocean (approximately 2 °C; (Brand et al., 2014).

The ongoing changes in the sea-ice cover extent and thickness have impacted the primary producers of the arctic and subarctic marine ecosystems (Kahru et al., 2016). Changes in the cryosphere have resulted in an imbalance in primary production between ice algae, under-ice, and open water phytoplankton in the Arctic Ocean and adjacent polar seas (Arrigo et al., 2012; Leu et al., 2015). This imbalance has been driven by decreasing sea-ice coverage and time shifts in seasonal sea-ice dynamics (i.e., ice and snow thickness, melt onset, pond onset, etc.; Arrigo et al. 2014; Horvat et al. 2017). Recent field observations (Arrigo et al., 2012) and model simulations (Horvat et al., 2017) suggest that the warming trends at high latitudes may have increased the occurrence of under-ice blooms as a consequence of the rise in under-ice solar radiation during the spring season due to thinner first-year ice, early melting, a significant fraction of leads, and melt ponds at the ice surface (Mundy et al., 2009, 2014; Palmer et al., 2014; Leu et al., 2015; Assmy et al., 2017; Horvat et al., 2017).

High phytoplankton biomass and primary production rates in the marginal ice zone (MIZ), as defined hereafter as the area along the edge of the ice pack that is affected by open ocean processes, have been reported in the literature over several decades (Barber et al., 2015, and references therein). Physical, chemical, and biological processes occurring at the MIZ trigger important ecological successions in the marine ecosystem and impact the coupling between sympagic, pelagic, and benthic realms (Leu et al., 2015). Satellite ocean color observations have proven a handy tool to map chlorophyll-a concentration [Chla], a proxy for sea surface pelagic phytoplankton biomass and abundance, in open waters found in the edge zone of the MIZ (see Barber et al., 2015) in both poles since the launch of the Coastal Zone Color Scanner (CZCS) (Maynard and Clark, 1986, 1987; Mitchell et al., 1991). Recently, Perrette et al. (2011) developed a method to assess the phytoplankton bloom along the edge of the ice pack, which is referred to as ice-edge bloom hereinafter, at the pan-Arctic scale using the Sea Wide Field-of-view Sensor (SeaWiFS). Based on the year 2007, they concluded that the ice-edge blooms were ubiquitous in the northern latitudes (>66.6 °N). Also using SeaWiFS, Lowry et al. (2014) reported that ice-edge blooms on the Chukchi Sea shelf depend on the timing of sea-ice retreat and further speculated that massive under-ice blooms were widespread in nutrient-rich waters of Pacific origin. Renaut et al. (2018) applied the method of Perrette et al. (2011) to 11 years of chlorophyll-a observation of the Moderate-Resolution Imaging Spectrometer Sensor (MODIS) and reported an intensification and a northward expansion of the ice-edge bloom in many subarctic seas. Altogether, these results suggest that nonlinear processes regulate ice-edge blooms. Among them, sea-ice production and snow cover thickness affect convective mixing and nutrient replenishment during winter, light availability at the onset of the growing season, and the amount of meltwater available for the subsequent emergence of meltwater stratification.

Are ice-edge blooms a recurrent feature throughout the HBS? To our knowledge, ice-edge blooms in the HBS have not been examined explicitly, using satellite ocean color or with field observations. The HBS was excluded from recent ocean color analyses conducted in the Arctic (Perrette et al., 2011; Renaut et al., 2018). Except in river plumes and some upwelling spots in Hudson Strait and near the Foxe Peninsula and the Belcher Islands, the HBS has been assumed to be a low-productivity oligotrophic system (Ferland et al., 2011; Tremblay et al., 2019). A few studies investigated under-ice phytoplankton and pelagic community interactions during the spring season in a few coastal locations (i.e., near the Belcher Islands and Great Whale River in the southeastern coastal HB; Michel et al., 1988, 1993; Runge et al., 1991; Monti et al., 1996). These land-fast ice camp-based studies reported many critical ecological processes occurring during the transition from the under-ice to a pelagic-dominated system: photo-adaptation of ice-algal communities (Michel et al., 1988); ice-algal communities seeding the pelagic ecosystem (Michel et al., 1993); the coupling between under-ice grazers, interfacial, and pelagic communities (Runge et al., 1991); and the impact of river plume dynamics on algal community composition (Monti et al., 1996). None of these studies found evidence of under-ice phytoplankton blooms during the melting season at the Hudson Bay basin scale. Model-based investigations by Sibert et al. (2010, 2011), however, reported that under-ice production timing is mainly controlled by surface melting (snow and ice) processes that determine the light levels, even under conditions of sufficient nutrient availability. A bay-wide assessment of [Chla] is needed to obtain a large-scale perspective of the phytoplankton productivity of the HBS.

In summary, there are many indications that the MIZ is a biologically productive feature in the HBS. Still, no systematic observations of ice-edge phytoplankton blooms have been reported in this subarctic region. The main objective of this study was to assess the temporal and spatial variability of the ice-edge bloom based on a systematic analysis of available satellite time series. Therefore, we analyzed two decades of satellite remote sensing data combining ocean color observation of [Chla] and passive microwave imagery for SIC. First, we compared the phytoplankton phenology assessed from archived observations of the CZCS (1978–1983) to the modern phenology obtained from MODIS-Aqua (2002–2014). Next, we characterized phytoplankton dynamics based on temporal evolution of satellite [Chla] and timing of sea-ice retreat to obtain climatological conditions, trophic categories, and phytoplankton phenological types in marginal ice zones. We documented the interannual variability of the surface [Chla] at the MIZ and examined how it is impacted by the timing of sea-ice retreat and winter air temperature. Finally, we examined the influence of large-scale climate variability patterns, particularly the Arctic Oscillation (AO) and North Atlantic Oscillation (NAO), on oceanographic processes and ice-edge bloom intensity in the HBS.

Multi-sensor merged [Chla] Level-3 (i.e., binned and mapped) 8-day composites from the Globcolour Project (http://www.globcolour.info/) were used as a proxy for phytoplankton biomass. Globcolour products have a spatial resolution of 4.63 km and cover the 1998–2018 period. The merged product was selected to improve the spatial-temporal coverage diminishing gaps due to cloud cover and sea-ice coverage (Maritorena et al., 2010). The binning methodology combines the normalized water-leaving radiances from different ocean color sensors whenever they are available, which includes SeaWiFS (1998–2010), MODIS-Aqua (2002–2018), Medium-Resolution Imaging Spectrometer (MERIS: 2002–2011), and Visible Infrared Imaging Radiometer Suite (VIIRS: 2012–2018). [Chla] was estimated from normalized water-leaving radiances merged using the Garver-Siegel-Maritorena (GSM) semi-analytical model (Garver and Siegel, 1997; Maritorena et al., 2002). GSM also yields particle backscattering (b_bp) and colored detrital matter (CDM) coefficients at 443 nm.

Colored dissolved organic matter (CDOM) and non-algal particles can deteriorate the performance of ocean color algorithms because, as phytoplankton pigments, they absorb the blue part of the visible spectrum (Bricaud et al., 1981). Consequently, in Arctic and sub-Arctic seas, [Chla] overestimates (underestimates) in the lower (higher) range of [Chla] have often been reported (Cota et al., 2004; Hirawake et al., 2012; International Ocean Colour Coordinating Group [IOCCG], 2015). CDOM is known to be a dominant optical component in most of the HBS, making this region optically complex (Granskog et al., 2007; Mundy et al., 2010; Guéguen et al., 2011; Xi et al., 2013; Burt et al., 2016; Heikkila et al., 2016). As discussed by Ben Mustapha et al. (2012), Bélanger et al. (2013), and Lewis and Arrigo (2020), the GSM algorithm can minimize the impact of other optical constituents on Chla retrievals because the more mechanistic GSM algorithm is able to take into account seasonal and regional variability of bio-optical properties compared to empirical algorithms for Polar shelves. Therefore, the GSM was selected because it can better represent the optically complex waters of the HBS (Xi et al., 2013, 2014, 2015), and data from depths shallower than 50 m were excluded from the analysis to avoid turbid or CDM-rich waters and river plumes.

We performed a cluster analysis (see below) on 8-day composites (multiyear averages) of Chla to detect potential changes in phytoplankton phenology between the 1980s and the 2000s. We used the CZCS for the period of 1978–1983 and MODIS for 2002–2014. For consistency in the data processing for these two ocean color missions, Chla was calculated using the case-1 waters standard empirical Ocean Color algorithm CZCS OC3 and MODIS OC3v5 (O’Reilly et al., 1998, 2000). The use of band ratio algorithms allows a better interconnection between distinct sensors (Antoine et al., 2005; McClain, 2009). We tested the sensitivity of the cluster analysis to the choice of Chla algorithms and available ocean color products. Specifically, we used the Globcolour merged chlorophyll-a products obtained using both GSM (as above) and empirical algorithms for the same period as MODIS OC3v5 alone (2002–2014). We found some spatial discrepancy in the cluster distribution between products, but the main conclusion drawn from the analysis remained unchanged (i.e., increased occurrence of double bloom, see Results section).

SIC was obtained from the National Snow and Ice Data Center. It is based on daily passive microwave radiometry processed using the Bootstrap algorithm (Comiso, 2000) at 25 km resolution. The Bootstrap technique clusters the multichannel passive microwave sensors: Scanning Multichannel Microwave Radiometer on the Nimbus-7 satellite, Special Sensor Microwave/Imager and Special Sensor Microwave Imager/Sounder from the Defense Meteorological Satellite Program’s satellites, and the Advanced Microwave Scanning Radiometer (Comiso et al., 1997). SIC was interpolated onto the same Chla grid using the nearest neighborhood scheme implemented in Matlab.

The Arctic Oscillation is a climate pattern characterized by winds circulating counterclockwise around the Arctic at around 55 °N latitude. The AO index is calculated from the first component of the empirical orthogonal function of monthly anomaly variations in sea-level atmospheric pressure north of 40 °N, which explains 22% of total variance (Thompson and Wallace, 1998). The North Atlantic Oscillation index is extracted from the difference between monthly anomaly variations in sea-level atmospheric pressure over Greenland (low pressure) and The Azores Island (high pressure) (Hurrell et al., 2001, 2003; Qian et al., 2008). NAO/AO data were obtained from the Climatic Prediction Center/the National Oceanic and Atmospheric Administration (http://www.cpc.ncep.noaa.gov). Here, we used the average index for the winter months of January, February, and March.

We obtained the monthly air temperatures from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) Reanalysis Project. Anomalies of air temperature were calculated as the difference between each monthly average and its corresponding monthly climatology and normalized by the standard deviation, both of which were calculated using data from the 1948 to 2018 period.

We investigated the potential changes in phytoplankton phenology between the 1980s and 2000s using a K-means cluster analysis applied on the annual time series of Chla. The approach used here is similar to that adopted by D’Ortenzio and D’Alcalà (2009) in the Mediterranean Sea, D’Ortenzio et al. (2012) at the global scale, Ardyna et al. (2017) in the Southern Ocean, and Marchese et al. (2017) at regional Arctic scales.

Briefly, we analyzed the seasonal variability (or shape) of the annual Chla time series (i.e., a multiyear average of 8-day Chla available for a given period, hereafter referred to as climatology) to determine the number of phytoplankton blooms during the annual cycle, as well as their timing. Therefore, we focused on the relative change in Chla over the year rather than the absolute values of Chla. In addition, we standardized the individual climatologies, Chl_clim(x, y, t), for each 4.60-km pixel (x, y), as:

where t is the time referring to the 8-day composite; $ChlN(x,y,t)$ is the normalized time series; and $Chl‾(x,y)$ and $σChl(x,y)$ are the annual mean and standard deviation values of climatological time series for the pixel (x, y). $Chlclim(x,y,t)$ was calculated with a temporal resolution of 8 days for CZCS and MODIS-Aqua. Data applied to phenology analysis were also filtered by a spatial-temporal median filter (3 × 3 × 3), and the gaps were filled using the data interpolating empirical orthogonal functions scheme implemented in R (DINEOF; Beckers and Rixen, 2003). The ideal number of clusters was determined by Davies-Bouldin criterion, and cluster stability was accessed by Jaccard, Bootstrap, and Noised coefficients (Hennig, 2007; D’Ortenzio and D’Alcalà, 2009; D’Ortenzio et al., 2012).

To assess the impacts of sea-ice retreat timing on MIZ phytoplankton blooms (also refers to phytoplankton spring blooms or ice-edge blooms), we analyzed both Chla and SIC variability in parallel. The method is similar to that of Perrette et al. (2011), which was also adopted by Lowry et al. (2014) and Renaut et al. (2018). The sea-ice retreat, t_R, is defined as the day at which SIC is below 10% for at least 24 days. This time interval is longer than the 20 days applied by Perrette et al. (2011) and Renaut et al. (2018) and the 14 days by Lowry et al. (2014) because we used 8-day composites instead of daily maps. However, to avoid sub-pixel contamination in ice-infested regions near the ice edge (Bélanger et al., 2013), we opted to be more conservative by applying a 10% threshold on SIC, as did Perrette et al. (2011) and Renaut et al. (2018) instead of 50% as applied by Lowry et al. (2014). The maximum Chla observed in the MIZ was extracted for each pixel for each year, yielding one map of MIZ Chla per year.

We recognized that sea-ice cover hides an important component of phytoplankton phenology that is out of the reach of ocean color satellites (Arrigo et al., 2012). To address this problem, we developed a trophic predictor of under-ice and pelagic phytoplankton phenology throughout the MIZ by analyzing a seasonal time series of satellite-derived SIC and Chla. This trophic predictor is defined by the temporal evolution of Chla since the sea-ice retreat (t_R).

Figure 1 illustrates five typical Chla phenologies expected during the spring-to-summer transition at the MIZ, assuming that the temporal evolution of Chla during phytoplankton blooms has a Gaussian-like curve (Jönsson and Eklundh, 2004; Platt et al., 2009). Each of these phenologies is defined in terms of three metrics: (1) the very first Chla observation after the sea-ice retreat (⁠$[Chla(tR)]$⁠), (2) the maximum Chla reached within five weeks after the sea-ice retreat (⁠$[Chla(t1,2,3,...)max]$⁠), and (3) the minimum reach within five weeks after the sea-ice retreat (⁠$[Chla(t1,2,3,...)min]$⁠). Table 1 presents the criteria that define each type of MIZ phenology. The bloom threshold, B, was set at 0.5 mg m^–3 as in Perrette et al. (2011). The oligotrophy threshold, Ol, was set at 0.2 mg m^–3. We applied these criteria to each pixel and for each year, yielding one thematic map per year.

To investigate the impact of climate and environmental forcing on ice-edge blooms, we calculated the Pearson correlation coefficient (r) of the linear relationships between the MIZ Chla and (1) the annual values of the AO, (2) the sea-ice retreat day (t_R), and (3) the winter (December–January–February, DJF) air temperature anomaly. The analysis was made on each grid cell of the HBS and between 1998 and 2018. We standardized the maximum Chla values in the MIZ (Perrette et al., 2011) using climatology mean and standard deviation as above. Only significant slopes, obtained from a least-squares fit at an interval of confidence of 95% (P > 0.05) and for pixels having more than 10 valid years between 1998 and 2018, were mapped.

Warmer air temperatures in the HBS began around 1998 (Figure 2), which coincides with the beginning of the modern era of ocean color radiometry (i.e., the launch of SeaWiFS in August 1997). This recent period contrasted with the CZCS era (1978–1983), which was characterized by the predominance of a negative air temperature anomaly (T_air) relative to the average air temperature calculated over the 1948–2018 period. One extreme cold event was observed at the beginning of the CZCS era in 1978. During the MODIS observational era, T_air was mainly positive, but colder winters were often seen, as in 2002, 2003, 2004, 2011, 2014, 2015, 2016, and 2018. In contrast, most summers showed positive anomalies, except 2004, 2007, 2011, and 2018. The latter stands out in the modern era with a noticeable negative T_air during most of the year.

The distinct T_air between the CZCS and MODIS eras coincided with remarkable changes in phytoplankton phenological patterns in the HBS (Figure 3). Our analysis resulted in three statistically distinct bioregions characterized by the shape and the number of observable Chla peaks, hereafter referred to as phytoplankton “blooms,” over an annual cycle. The spring–summer bloom (SB) type (red pixels and curve in Figure 3) exhibited only one Chla peak, which reached its maximum a few weeks after the summer solstice (June 21). The double bloom (DB) type (green pixels) was characterized by a first Chla peak occurring in late May/early June, and a second peak at the very end of the open water season, when ocean color observations were still available (i.e., mid-October). The last type, the fall bloom (FB) (blue pixels), was initiated in mid-summer, peaked in mid-September, and declined until the sea-ice recovered (Figure 3C). FB contrasted with the second peak of the DB phenology, in which Chla increased continually until the sea-ice recovery.

Three types of bloom were observed during the CZCS era, but with no distinct spatial patterns (Figure 3A). Of 65,873 valid observations over that period, the SB, DB, and FB occupied 28%, 38%, and 24%, respectively. During the MODIS era, for which 70,972 valid pixels were found, SB remained at about 28%, but DB increased to 54%, and FB decreased to 18% of the total surface area. This result suggests that a remotely sensed spring or summer bloom is a common feature in the HBS. Whether this increase in Chla at the beginning of the open water season is an actual ice-edge bloom or not is examined further in the next sections.

Single Chla peak phenology (SB and FB) might be seen as a variation of DB phenology. The start date and the peak amplitude of single SB occurred shortly after the spring bloom of DB phenology. Similarly, the single FB peak occurred before the fall bloom of the DB phenology. The fact that the DB phenology has become a predominant feature in the HBS implies that an earlier spring–summer bloom is now occurring. We examined whether or not the duration of the open-water-season (OW), which was defined as the number of days with SIC below 10% during an annual cycle, could explain the fact that more DB are now observed. The OW duration increased from 15 to 39 days between CZCS and MODIS eras, reaching about 125 days (4 months) in modern times (Figure 3D and Table 2). The regions where SB in the CZCS era was replaced by DB in the MODIS era had an increase in OW of 34 days (Table 2). The OW duration in areas where FB was replaced by DB or remained as DB had increased by 21 and 23 days, respectively. However, these increases in OW were not more significant than for pixels showing SB in both CZCS and MODIS eras (39 days). Therefore, no clear link was identified between the OW duration and the change (or the lack of change) in the phytoplankton phenology patterns.

Figure 4A shows the climatology of the 8-day bins of Chla for the whole HBS from the Globcolour data set (1998–2018), together with the average daily percentage of OW in the HBS (blue bars). Interestingly, we observed two distinct peaks in Chla: one in spring around day-of-year (DOY) 75 (mid-March) and a second around DOY 150 (beginning of June). The first peak, reaching values around 0.7 mg m^–3, is observed when sea ice covers the bay almost entirely (90% of the area). During that time, we noted early phytoplankton blooms in very restricted locations known as polynyas, such as the northwestern part of Hudson Bay (Ferguson et al., 2010). The second increase in Chla begins as soon as the OW area starts to increase and can be considered as an ice-edge bloom, which is a typical feature along the marginal ice zone during the spring-to-summer transition (Mitchell et al., 1991). During this transition period, the Chla reached a relatively high level (1.0 mg m^–3) with lower variability, as depicted by the lower standard deviation (gray-shaded area on Figure 4A) compared to the rest of the year. This bloom is probably productive, as it is synchronized almost perfectly with the peak of incident photosynthetic available radiation (⁠$PAR(0+)$⁠, where 0+ indicates PAR just above the sea surface) (Figure 4B). The end of the OW season is characterized by a continuous increment of Chla, referred to as the fall bloom, which ended when sea ice recovered in late October (DOY 300). The high standard deviation during this season indicates extremes of low and high Chla, pointing toward high spatial variability.

To illustrate the recurrence of ice-edge blooms in the HBS, we plotted the frequency distribution of Chla (N approximately 10⁶ pixels) as a function of number of days after the sea-ice breakup between 1998 and 2018 (Figure 5). Higher Chla values are generally restricted to the first 8-day period after the sea-ice retreat, with Chla values ranging from 1 to 3 mg m^–3. In contrast, observations of high chlorophyll-a later during the open water season were relatively low. These results suggest that ice-edge blooms, which are characterized by an increase in phytoplankton abundance propagating along the ice edge, could be observed frequently by ocean color satellites in the HBS. However, spatial analysis is required to confirm whether ice-edge blooms are present everywhere in the HBS or whether they are restricted to some parts of the Bay, where ocean conditions are more favorable, as suggested by the phenology analysis based on the climatology (Figure 3).

Panels A and B of Figure 6 show maps of the average and the standard deviation of the maximum Chla, respectively, observed in the ice-edge zone of the MIZ between 1998 and 2018. The statistics were calculated using the 21 annual maps depicted in Figure S1. The limits of the subregions are those used by Landy et al. (2017) who reported satellite-based sea-ice thickness and snow depth in the HBS. In general, we noted higher Chla in the surrounding of the Bay, with the highest values found in the Hudson Strait (HS) and Southeastern (SE) HB. These two regions also show a higher standard deviation (STD) in the maximum Chla at the MIZ, indicating relatively high interannual variability. We also observed a high STD in the northwestern (NW) polynya. Southwestern (SW) and Northeastern (NW) HB have moderate to high Chla (around 1 mg m^–3) in the MIZ and relatively low interannual variability. The central HB shows both low Chla (<0.5 mg m^–3) and STD interannual variability, suggesting that ice-edge blooms are not a dominant feature offshore. Figure 6C and D illustrates the average and standard deviation of sea-ice retreat date (t_R) for the same period (1998–2018) based on annual maps shown in Figure S2. The central HB had the latest t_R, ranging from DOY 180 to 210 (month of July), relatively low interannual variability (STD < 10). The highest (lowest) interannual variability of t_R occurred in the coastal (offshore) areas of SW and NE Hudson Bay. Interestingly, the t_R and Chla were somewhat very similar, suggesting that the magnitude of the ice-edge bloom is linked to the timing of sea-ice retreat, as explored further below.

Table 3 presents the statistics and frequency of occurrence of the surface extension (in km² and % of the total surface area) for four subcategories of Chla, which is a good proxy for the trophic state of the MIZ. On average, an ice-edge bloom (0.5 to 2.0 mg m^–3) or major ice-edge bloom (here major bloom is defined as Chla > 2 mg m^–3) occupied 46.3% of the Bay, but blooms can range from 27.9% (1999) and 56.2% (2008) of the total area (see also Figure S1). High interannual variability was also observed in the area of major blooms ranging from 2.3% (2010) to 9.0% (2015) (Figure 7). During the study period (1998–2018), the surface area of the ice-edge blooms increased by 10.5% (i.e., 0.5% per year times 21 years). In 2015, a major ice-edge bloom occupied as much as 7,200 km² (or 9.0% of the HBS surface) but was found in the polynyas (Figure 7). In 2018, the timing of the sea-ice breakup showed large spatial variability with an early breakup in the NW and a late breakup in the NE (Figure 7). The ice-edge bloom occupied a large portion of the NW, but with modest Chla compared to 2015. This result suggests that the magnitude of the ice-edge bloom is explained not only by the timing of the sea-ice retreat.

Figure 7 shows annual maps of the maximum Chla observed during the MIZ time window for selected years (2008, 2011, 2015, and 2018), showing interesting features in terms of ice-edge blooms. Several key observations can be made from these maps. First, the Chla in the central HB remained relatively low (<0.5) after the breakup for most of the years, except in 2005, 2008, 2012 and 2016, which means that no ice-edge blooms occurred in the central HB for most of the years analyzed. Moderate to high values of Chla were observed during the MIZ period along the surrounding of the bay with obvious variability in both northwestern and northeastern parts of the bay, which sometimes showed opposite patterns (e.g., compare 2015 and 2016; Figure S1). Higher Chla in the MIZ was observed systematically in the Hudson Strait where mixing and nutrient replenishment are more intense relative to the rest of HB (Ferland et al., 2011). However, the largest ice-edge bloom was observed in the NW HB in 2015 with Chla higher than 2.0 mg m^–3.

To gain insight into the variability in the MIZ Chla (Figure 7), we considered the timing of the sea-ice retreat. Figure S2 depicts the DOY of the sea-ice breakup or, in other words, the beginning of the MIZ period in the HBS for each year between 1998 and 2018. The breakup occurred within a wide range of 4 months, that is, from DOY 120 to 240. An early breakup can be found either in the NW or NE HB, where polynyas are common. Sea ice often accumulated in the southern part of HB delaying the breakup to the end of the summer, as observed in 2000, 2004, 2008, 2009, and 2015. A visual comparison of the annual maps presented in Figures 7, S1, and S2 suggests that, at least in some parts of Hudson Bay, early sea-ice breakup resulted in more intense ice-edge blooms observable from space. This relationship, however, is far from applicable to the whole HBS, which we examine further in the next section.

Figure 8 shows maps of the frequency of occurrence of the four trophic states in the ice-edge zone defined above (Table 3). These maps provide an additional point of view about the spatial and interannual variability of ice-edge blooms in the HBS between 1998 and 2018. Oligotrophy occurrence (<0.2 mg m^–3) was more frequent offshore, covering a vast region in the central part of HB in 30% of the years analyzed. On the other hand, the occurrence of oligotrophy conditions across the bay was less than 5% (Figure 8A). As observed in Figure 8B, moderate oligotrophy (0.2 ≤ Chla < 0.5 mg m^–3) conditions were more frequent in the central HB with a frequency higher than 60% covering a large proportion of the bay (almost half of it). Ice-edge bloom (0.5 < Chla ≤ 2.0 mg m^–3) occurrence in the surrounding of the bay and HS generally surpassed 65%. A presence higher than 80% was reached in the domains of the NW polynya, proximal areas of the most significant river plumes, the eastern bay, and HS (Figure 8C). Recurrently productive hot spots (>50% of occurrence; Figure 8D) were found in relatively small regions, south of Belcher Island, some places in the northern HS, and river-influenced coastal areas. We also observed occurrence above 10% in the NW HB polynya and south of Southampton Island.

We performed a time-series analysis of SIC and Chla at the pixel level to assess the potential that under-ice blooms escaped satellite range. We considered four scenarios for the temporal evolution of Chla at the MIZ, as shown in Figure 1. Applying the criteria of Table 1 resulted in 21 thematic maps presented in Figure S3. To synthesize the results, we calculated the frequency of occurrence of each type based on 21 years of observations. Figure 9 shows the frequency of occurrence of the five phenology types between 1998 and 2018.

Consistently low biomass observed in the MIZ after the ice retreat (Figure 9A) could be due to an under-ice bloom that escaped the satellite range entirely or to an oligotrophic state at the sea surface. Similarly, the probable (recent) under-ice bloom (Figure 9B) may also indicate that oligotrophic domains remain in the ice-edge zone. However, these domains involve moderate concentrations between 0.2 and 0.5 mg m^–3 after the sea-ice retreat followed by a decline in Chla. Oligotrophy/old under-ice bloom (Type I) and probable (recent) under-ice blooms (Type II) can be considered similar features, though with distinct magnitudes and bloom-timing. This scenario (Types I and II) dominated the whole central HB, extending north of Southampton Island with a frequency of occurrence ranging from 50% to 100%. In contrast, the situation characterized by Chla consistently above 0.5 mg m^–3 (Figure 9C), suggesting mesotrophic conditions where surface waters are efficiently replenished in nutrients, was frequent only in coastal waters or at the ends of the Hudson Strait. The last scenario, that is, when the bloom is triggered in ice-free waters, was infrequent, except at the margin between the coastal domains and the central HB, from Churchill to NE HB (Figure 9D). Finally, the typical ice-edge bloom scenario, that is, when Chla peaks just after the ice breakup and then vanishes to low concentration, was relatively frequent (30%–70% of the time) all around HB and in the central part of HS (Figure 9E).

As shown above, a high spatiotemporal variability characterized the dynamic of the MIZ, in terms of both Chla and sea-ice breakup in HBS (Figures 7, S1, and S2). Here, we examined the potential link between these two variables by aggregating data in the seven subregions presented above (Figure 10): northwestern (NW), southwestern (SW), central and offshore (Central), narrows at the Hudson Bay entrance (Narrows), Hudson Strait (HS), northeastern (NE), and southeastern (SE). The boxplots of Figure 10 show the Chla in the MIZ as a function of the date of the sea-ice retreat (t_R). The frequency occurrence of t_R is also shown on top of each boxplot. These distributions allow an understanding of the dynamic of the ice-edge bloom in each subregion of the HBS.

The western part of the bay (NW and SW; Figure 10) showed a similar relationship between the timing of the ice breakup and the magnitude of the ice-edge bloom. The sooner the sea ice breaks up, the higher the Chla in the ice edge. The NW subregion, known to host wind-driven polynyas (Landy et al., 2017), shows more frequent sea-ice retreat in June (86.4% of the time), though it was also often observed in May (21.7% of the time). The median Chla values were >1 mg m^–3 when the breakup was in May, and approximately 0.7 mg m^–3 if it occurred during the first half of June. We observed oligotrophic conditions whenever the ice melt happened after June 15, with Chla well below the threshold of 0.5 mg m^–3. In the SW in May, however, ice-edge blooms only occurred 4.2% of the time. Sea-ice breakup happened mainly between mid-June and mid-July (60.3% of the time) and showed ice-edge Chla very close to the 0.5 mg m^–3 threshold. The influence of the Churchill and the Nelson rivers may explain the larger Chla in summer compared to the NW subregion where river inputs were minimum in the North. The increase in Chla in the ice edge after August 1 only represented 0.3% and 2.1% of the time in the NW and SW, respectively, and thus was negligible.

The Narrows of HB entrance and the central HB (Figure 10) showed a similar relationship between the timing of the breakup and the magnitude of the ice-edge bloom. First, the sea-ice retreat happened mainly in late June (approximately 50% of the time) or early July (approximately 30% of the time) but never before mid-May. Second, when the sea-ice retreat between mid-May and mid-June (approximately 11% of the time), moderate phytoplankton blooms characterized the ice edge with a median Chla around 1.0 mg m^–3. Third, when the sea-ice retreat occurred after mid-June, it was followed by Chla approximately 0.5 mg m^–3. Finally, slightly higher Chla was found in the ice edge when the sea ice retreated in late July, but this scenario represents only 7% or 8% of cases.

In the eastern part of the bay (NE and SE; Figures 10), Chla remained well above the bloom threshold regardless of the date of the sea-ice retreat. We found higher Chla in the SE, compared to NE, just south of the Belcher Islands (Figure 6A). Similarly, the Hudson Strait should be considered as a highly productive system because the median Chla in the ice-edge zone remained close to 1.5 mg m^–3 when sea ice retreated between mid-May and the end of July. There, early ice retreat (before mid-May) resulted in lower Chla compared to the rest of the season.

In summary, the timing of sea-ice retreat seems to play a crucial role in setting the magnitude of the ice-edge bloom observed from space in most of the HBS (NW, SW, Central, and Narrows), except in the eastern part of the bay (SE, NE, and HS). Significant ice-edge blooms were more restricted to May in the NW and SW subregions where the sea-ice retreat was earlier compared to the rest of the HBS (Figure 6C).

The response of ice-edge blooms to the timing of sea-ice retreat showed some distinct spatial patterns. Figure 11A shows the slope $(∂ChlaIEZ∂tR)$ of the linear regression between the maximum Chla in the ice-edge zone and t_R on a pixel basis. The negative (positive) slope indicates that the sooner (later) the ice retreat, the higher the Chla in the ice-edge zone. Interestingly, we found strong negative slopes in the Northwest HB, the Narrow at the HB entrance, and in the central HS. Negative slopes were also found in about half of pixels in the NE subregion. Some coastal spots in the south presented positive slopes, which are associated with river discharge. The increment of Chla with delays in the sea-ice retreat was also observed southeast of Belcher Islands, where tidal mixing is particularly intense (Webb, 2014). However, these domains were also associated with relatively high root mean square errors (RMSE above 1.00 mg m^–3), suggesting large uncertainty in the relationship. Finally, we noted a neutral response of ice-edge blooms to t_R in most of the central HB.

We computed a winter temperature index using the normalized air temperature anomalies relative to the 1948–2018 mean (⁠$Tair′$⁠). A similar approach was used by Hochheim and Barber (2014) to investigate the impact of surface air temperatures on sea-ice breakup in the HBS. Here, the response of ice-edge bloom magnitude is discussed in terms of the positive, negative, and neutral slope of the relationship between the winter air temperature anomaly and the maximum Chla in the ice-edge zone (⁠$∂ChlaIEZ/∂Tair′$⁠; Figure 12A). The positive (negative) slope indicates that the warmer (colder) the winter, the higher (lower) the Chla in the ice-edge zone in the spring–summer season. Based on this analysis, we found a positive relationship between air temperature and ice-edge bloom in the Hudson Bay entrance and NE (red pixels in Figure 12A), with NW and SE subregions appearing more sensitive than the rest of the bay (blue pixels). This positive relationship means that higher magnitude ice-edge blooms are expected after a cooler winter. Some spots, for example, near large river plumes (Churchill, Nelson), had a negative relation between ice-edge blooms and $Tair′$⁠, but were also associated with the highest RMSE. This result may be related to river discharge variability (timing or volume). In contrast to the HB, the Hudson Strait presented a complex spatial distribution concerning the relation between ice-edge Chla and $Tair′$⁠, with strong negative slopes observed in most of the area.

We subsequently examined the impact of climatic teleconnections, such as the North Atlantic Oscillation and Arctic Oscillation, on the phytoplankton dynamic in MIZ, by correlating the maximum Chla in the ice-edge zone with the NAO/AO index (Thompson and Wallace, 1998). As illustrated in Figure 13A and B, the main pattern observed is a positive correlation between NAO/AO and ice-edge zone Chla in the NW and south of Belcher Islands polynyas, where sea ice usually retreats earlier (Landy et al., 2017). As reported in the Figure 13C (D), the time series of normalized ice-edge Chla and NAO (AO) index correlated significantly, with a Pearson correlation coefficient of 0.27 (0.29) in a confidence interval of 95% for the whole HBS and 0.46 (0.41) for the NW polynya. The correlation maps between ice-edge blooms, AO (Figure 13A), and NAO (Figure 13B) showed a similar pattern: a vast region of high correlation controlled by the NW polynya dynamic, negative correlation in the HB and NE bay, and high correlation south of Belcher Island. The maximum value of both Chla, AO, and NAO indices during the 1998–2018 period occurred in 2015 when we observed a massive phytoplankton bloom in the whole HBS and the NW polynya (Figures 7 and 13C, D). In contrast, the lowest NAO/AO indices occurred in 2010, preceding a negative anomaly in ice-edge Chla in the HBS by one year but simultaneous with a local minimum in the NW polynya (Figure 13C vs. 13D). The primary mismatches with NAO/AO indices occurred in 2005 and 2017 (Figure 13C) when considering the whole HBS and 2000 and 2017 in the NW polynya (Figure 13D).

We first discuss the inherent limitations of ocean color data for describing the phytoplankton dynamic in an ice-infested region under the influence of significant riverine inputs. Then, we explore the variability of Chla in the ice-edge zone in terms of potential physical forcings that control nutrient availability and incoming light in the upper layer and consequently the magnitude of the bloom. We further examine the impact of climatic changes on the phytoplankton dynamic throughout the MIZ. To achieve this objective, we explore the relationships between ice-edge bloom intensity and climate indices (NAO/AO), winter air temperature, and sea-ice retreat. Climate indices are used to describe the planetary teleconnections between the northern polar vortex. Finally, we focus on the HB NW polynya, an essential hot spot for marine life, as it shows high sensitivity to large-scale climate forcing.

Satellite remote sensing provides high revisit frequency, global coverage, and continued operational monitoring capability for access-limited regions such as HBS (Babin et al., 2015; Lee et al., 2015). Despite inherent limitations (IOCCG, 2015), that is, restricted vertical range to the first optical depth (1/K_d), lack of data under clouds or sea ice, presence of terrigenous optical components and sea-ice contamination of ocean color products, satellite data have proven efficient to describe the phytoplankton dynamic at high latitudes, including in the MIZ (Maynard and Clark, 1987; Perrette et al., 2011; Lowry et al., 2014; Renaut et al., 2018).

Notwithstanding, standard empirical Chla algorithms for global ocean data, and even Arctic-adapted ones, can be biased (IOCCG, 2015), calling into question the choice of a fixed threshold to detect phytoplankton blooms. In the Arctic and sub-Arctic seas, empirical algorithms tend to overestimate (underestimate) Chla in the lower (higher) range (Cota et al., 2004; Hirawake et al., 2012; IOCCG, 2015).

In particular, CDOM and non-algal particles can diminish the performance in these algorithms because, as phytoplankton pigments, they absorb the blue part of the visible spectrum (Bricaud et al., 1981). CDOM is known to be a dominant optical component in most of the HBS, making this region optically complex (Granskog et al., 2007; Mundy et al., 2010; Guéguen et al., 2011; Xi et al., 2013; Burt et al., 2016; Heikkila et al., 2016). The GSM algorithm can better deal with the presence of CDOM compared to empirical algorithms and was selected to minimize the impact of other optical constituents on Chla retrievals (Ben Mustapha et al., 2012). The nearshore region where river inputs contribute more to the water reflectance variability than phytoplankton was avoided by masking shallow water pixels (i.e., depth < 50 m). Granskog et al. (2009) found that river runoff and CDOM are constrained within the coastal domain (100–150 km from the shore), where they are transported along the coast counterclockwise around the HBS. Xi et al. (2013) reported that phytoplankton was the main optically significant component driving the seasonal variability of light absorption in the HBS, while CDOM remained relatively stable over the annual cycle. These results provide confidence in the present Chla time-series analyses. However, we cannot exclude an influence of riverine-derived CDOM on Chla, in particular along the eastern coast of HB, where most river runoff is transported before its exit through the Hudson Strait (Saucier et al., 2004). A high CDOM background in the east of HB (Xi et al., 2013) explains the relatively high Chla observed in the ice edge regardless of the sea-ice breakup in the eastern part of HB (NE and SE in Figure 10). An evaluation of the GSM algorithm based on in situ measurements in the HBS beyond the scope of this study but would be necessary to evaluate the potential bias in Chla and consequently the choice of 0.5 mg m^–3 as a threshold for detecting pelagic blooms.

Here, we adopted a fixed Chla-based threshold of 0.5 mg m^–3 to detect the ice-edge bloom (Perrette et al., 2011) instead of dynamic methods like, for example, rate of change method, cumulative biomass-based thresholds (Brody et al., 2013), or model fitting (Ardyna et al., 2014; Marchese et al., 2017). Because the ice-edge bloom can reach its peak just after the sea-ice retreat (Figure 1 and type V in Figure 9), the pelagic cycle often began with the highest annual values, without a well-defined Gaussian-like seasonal evolution in Chla. Therefore, we could not apply the Gaussian fitting method to the HBS.

What role do winter air temperatures (⁠$Tair′$⁠) and sea-ice retreat (t_R) play in the magnitude and variability of ice-edge blooms? Can these relationships be affected by climate change? We used 21 years of SIC and Chla to address these questions.

As illustrated in Figures 11 and 12, t_R and $Tair′$ can have a very strong but local impact on the balance between under-ice and pelagic production throughout the MIZ (Palmer et al., 2014). This impact can be explained by the interactions of various oceanographic processes in the MIZ, such as tidal resonance (Webb, 2014), water mass exchange throughout the Hudson Strait (Sutherland et al., 2011), polynya dynamics (Landy et al., 2017), and freshwater inputs (St-Laurent et al., 2011). Notably, the relatively strong negative correlation between Chla and t_R in the northern part of the HBS illustrated in Figure 11 is very similar to the Arctic water intrusion pattern in the system depicted by Wang et al. (1994), Saucier et al. (2004), and Piecuch and Ponte (2015).

In winter, primary production in the HBS is limited by photosynthetic available radiation, as confirmed by the absence of dinoflagellate cysts (a proxy for under-ice production) reported by Heikkila et al. (2016) in two coastal sites in both west and east HB. On the other hand, physical processes controlling the cryosphere, that is, brine production and vertical mixing, are essential to bringing deep nutrient-rich waters to the upper ocean. These processes can be critical in the HBS because it is a relatively shallow sea (Stewart and Lockhart, 2005; Granskog et al., 2011). Also, compared to other polar domains of the Arctic, the sea ice is thinner, which may result in less brine rejection and relatively weak vertical mixing and nutrient replenishment during winter.

Winter air temperature is a good proxy for sea-ice production and the overall winter severity. Hochheim and Barber (2014) reported that a positive anomaly of only 1°C decreased the sea-ice coverage by 14% and advanced the breakup by one week. However, this generalization can be ambiguous concerning its effect on phytoplankton productivity because winter air temperature can impact many different key processes. Those include the nutrient replenishment due to sea-ice or brine production (Landy et al., 2017) or winter convection in the open waters (Stewart and Lockhart, 2005; Granskog et al., 2011; St-Laurent et al., 2011), the PAR attenuation by sea ice and snow cover (Horvat et al., 2017), the timing and length of the open water season (Ardyna et al., 2014; Hochheim and Barber, 2014), and the riverine input (Déry and Wood, 2004; Déry et al., 2005). The latter controls the amount of CDOM and nutrients (Granskog et al., 2007), which have a contrasting impact on phytoplankton production.

The complexity of air, sea, and ice interactions at high latitude makes the interpretation of the relationship between ice-edge Chla and winter air temperature difficult (Figure 12A). Indeed, various feedbacks are in play between polynya dynamics, air temperature, vertical mixing, brine production, and phytoplankton blooms (Chaudhuri et al., 2014). For example, Saucier et al. (2004) reported that the NW polynya dynamic causes significant sensible heat loss at a rate of 100 W m^–2 in winter. The presence of a polynya during the cold season tends to intensify the sea-to-air heat flux, which potentially warms the atmosphere locally (Saucier et al., 2004). Similarly, the Weddell Sea polynya, Antarctica, can effectively heat the atmosphere at the regional scale (Campbell et al., 2019). However, the polynya may only affect local air temperatures, which may be difficult to assess, given the coarse resolution of the NCEP/NCAR reanalysis data (Winsor and Björk, 2000).

The strong negative relationship between ice-edge Chla and winter air temperature in the HB polynyas (i.e., NW and Belcher Island) reflects the intensification of polynya dynamics, which leads to efficient nutrient replenishment in the euphotic zone (Figure 12). Atmospheric cooling of the ocean also influences the amount of ice and brine production, which further increases vertical mixing and weakens the under-ice stratification (Stewart and Lockhart, 2005; Granskog et al., 2011). Consequently, the nutrient stock available at bloom onset may be greater (Nguyen et al., 2009; Barthélemy et al., 2015). Indeed, thick sea ice stocks a large quantity of freshwater, which provides low-density water and stabilizes the upper ocean upon sea-ice melting (Galbraith and Larouche, 2011).

In summary, an inheritance effect of winter sea-ice production could indirectly impact the level of phytoplankton primary production that will occur in the spring–summer season. Therefore, we conclude that cooler air temperatures intensify the wind-driven polynyas of the HB due to steady atmospheric cooling of the ocean, which supports ice (and brine) production, vertical mixing, and deep nutrient replenishment during winter. As a result, as soon as the stratification is reestablished and light becomes available in spring, phytoplankton can reach higher biomass.

A warmer climate has already been established in the HBS due to warmer air and less frequent cold and dry polar air masses (Leung and Gough, 2016). Based on air temperature analysis between 1948 and 2018 (Figure 2), positive anomalies have been predominant in most seasons since 1998. Impacts on the cryosphere dynamics have been observed: thinner sea-ice cover, early sea-ice retreat, late ice recover, and longer open water season (Hochheim and Barber, 2014; Kowal et al., 2017). These changes in the sea-ice dynamics may have driven significant changes in the marine ecosystems, with potential cascading impacts on all trophic levels in the HBS (Hoover et al., 2013; Keller et al., 2014; Andrews et al., 2017).

Mechanisms controlling the spring–summer and fall blooms are distinct. Our results suggest that pelagic production may have increased in the HBS since the open water season has increased in duration and double blooms are more frequent (Figure 2). Sea-ice dynamics and incoming PAR in the upper ocean layer control the spring–summer bloom. Kahru et al. (2011) reported a pronounced trend toward earlier phytoplankton blooms in about 11% of the North Polar seas area between 1997 and 2009, including the northern part of the HBS. The latter can be related to earlier ice breakup trends (Kowal et al., 2017). In contrast, the lengthening of the open water season and the increase in nutrient pumping by wind-driven turbulence (Ardyna et al., 2011) and the inversion of heat fluxes drive the fall blooms. Notwithstanding, the increase in annual phytoplankton primary production is correlated with the length of the open water season (Arrigo and van Dijken, 2015).

The HBS is located in a transition between subpolar and polar domains. However, its complex topography and its semi-confined configuration by the continental shelf result in relative isolation with regard to global ice exchange and ocean circulation (Sutherland et al., 2011). Its primary connection to the global ocean takes place through Hudson Strait. Therefore, HBS trends and dynamics are especially impacted by local atmospheric processes (Hochheim and Barber, 2010, 2014; Hochheim et al., 2011) and by climatic teleconnections (Déry and Wood, 2004). Therefore, climatic forcing and global scale teleconnections (e.g., NAO/AO) can impact both sea-ice retreat and recovery in HBS (Hochheim and Barber, 2010; Hochheim et al., 2011).

Climate indices express global or hemispherical teleconnections (Carmack et al., 2006). Part of the atmospheric and oceanic variability of the HBS has been correlated with both AO (Déry and Wood, 2004; Hochheim and Barber, 2010; Hochheim et al., 2011) and NAO indices (Qian et al., 2008). NAO is recognized as a significant descriptor of the winter-to-winter variability over northeastern America (Wallance, 2007; Qian et al., 2008). It accounts for the direct North–South dipole patterns between the North Atlantic and the Arctic oceans. Meanwhile, AO (annular paradigm) also includes the Pacific sector processes, such as the Aleutian High (Wang et al., 2005; Wallance, 2007). Nevertheless, there is an overlap between NAO and AO patterns in the Atlantic sector. Therefore, NAO and AO are highly correlated and related to similar phenomena (Ambaum et al., 2001; Wallance, 2007), especially during winter (Rogers and McHugh, 2002). Both NAO and AO indices are associated with the intensity of the north polar vortex (Wallance, 2007), with a positive phase reflecting the strengthening of the winter-time polar vortex and westerly intensification between 50 °N and 70 °N (Moritz et al., 2002). For example, Hurrell et al. (2001) reported that positive NAO strongly affects the Atlantic Ocean, causing substantial changes in surface wind patterns with stronger westerly winds encircling the North Pole. In contrast, Arctic warming and a wavier polar jet resulting in an intensification of southward advection of polar air masses in North America were observed during negative NAO/AO phases (Baldwin and Dunkerton, 1999; Ding et al., 2014; Francis and Vavrus, 2015; Leung and Gough, 2016; Meleshko et al., 2016). We found that the ice-edge zone Chla and NAO/AO indices are relatively well correlated, especially in the NW polynya (Figure 13). The extreme event of the highest (lowest) NAO/AO in 2015 (2010) corresponded to a strong (weak) ice-edge bloom in NW HB polynya (Figure 13). These results suggest that global atmospheric circulation patterns, depicted by the NAO/AO, strongly influence phytoplankton blooms in the NW HB polynya.

Figure 14 illustrates schematically how the strength of westerlies affects the NW HB polynya processes controlling ice-edge blooms during winter and spring-to-summer transition. As discussed by Saucier et al. (2004) and Landy et al. (2017), the NW polynya is maintained by strong westerly winds opening up areas of water along the northwestern coast by coastal divergence and enhancing ice production. In winter, the NW polynya acts as an “ice factory,” where ice growth is favored by thermodynamic processes (ocean cooling) and then the ice is exported by winds (Landy et al., 2017). Brine rejection due to sea-ice production plays a fundamental role in vertical stratification and deep convection (Prinsenberg, 1988) and results in the nutrient replenishment of the euphotic zone through vertical mixing. Meanwhile, we found that early sea-ice retreat results in more intense ice-edge blooms (NW in Figures 10 and 11). This finding suggests that nutrient stocks set up in winter and the timing of polynya expansion in the spring-to-summer transition are the main drivers of the balance between pelagic and under-ice production throughout the MIZ. Notably, NAO/AO events can control the feedback between westerlies, the NW HB polynya dynamic, and the ice-edge blooms because initially brine rejection forced by ice production and export sets up the preconditions of nutrients stocks in the euphotic zone. After that, early polynya expansion due to intensification of westerlies results in ice-edge bloom intensification during the spring-to-summer transition (Figure 10-NW).

Déry and Wood (2004) highlighted that AO could explain 70%–90% of the variance in HBS river discharge anomalies. According to Déry and Wood, the origin of dominant air masses over the HBS controls this pattern. During the positive (negative) AO index, north easterly (north westerly) winds that advect relatively warm (cool), moist (dry) air from the Labrador Sea to the HBS increase the river discharge (Déry and Wood, 2004; Hochheim and Barber, 2014). We found weak correlations between AO and ice-edge zone Chla in the southern HB, where riverine input into the HBS is more important (St-Laurent et al., 2011). We also found no significant correlations between the maximum Chla in the ice-edge zone and riverine discharge, even at local scales (results not shown), suggesting a marginal role of rivers in the primary productivity of the HBS.

Recent observations have suggested that under-ice phytoplankton blooms can draw down the surface nutrient pool before the open water season starts (Mundy et al., 2009; Arrigo et al., 2014). Then, the subsurface chlorophyll-a maximum sinks near the nitracline before or just after the ice retreat, allowing the phytoplankton to achieve a balance between PAR and nutrients, as illustrated in Figure 14. As a result, the ocean color satellites detect a low Chla in the ice-edge zone (Horvat et al., 2017; blue and red curves in Figure 1). The frequency occurrence map depicting this situation (Figure 9A and B) suggests that a vast portion of the central-western part of the HBS is experiencing under-ice blooms. However, this scenario may also reflect the oligotrophic nature of the system.

Only in situ observations or model simulations can help to resolve the above situation (under-ice bloom vs. oligotrophy). Most field observations before the BaySys project were based on ice camps, which focused more on ice algae or open waters later in summer. Michel et al. (1993) and Michel et al. (1988) reported a maximum of vertically integrated ice-algal Chla of 23.6 mg m^–2 in April, while Monti et al. (1996) reported algal concentrations at the ice/water interface as high as 100.6 mg m^–2 between April and May. High ice-algal concentrations were followed by an increase in water column Chla reaching moderate concentrations (i.e., 2 and 4 mg m^–3) in May during the ice melt (Runge et al., 1991; Michel et al., 1993), suggesting an initiation of the bloom under ice (type V in Figure 1). Model simulations also suggest that nutrient draw-down in the surface waters begins in May, but the maximum diatom biomass is reached after the ice retreat in June or July (Sibert et al., 2011). Recently, however, Tremblay et al. (2019) concluded that, except in river plumes and some upwelling spots (e.g., Foxe Peninsula, Belcher Islands and Hudson Strait), the HBS should be considered as an oligotrophic system. More observations in the MIZ during the spring-to-summer transition are still needed.

We showed that phytoplankton phenology in the HBS is subject to substantial spatiotemporal variability and closely linked to large-scale atmospheric forcings. In recent years, the phenology has been characterized by two peaks in Chla, one after the sea-ice breakup in the marginal ice zone (May–June) and one in the fall. Here, we focused our effort on the ice-edge bloom. In the western part of the HB, the magnitude of the ice-edge bloom depends on the timing of ice breakup, with more intense blooms occurring when ice retreats early. The northwestern polynya stands out as a significant feature in the HBS in terms of primary production, and its variability is highly coupled to the northern hemisphere climate variability (NAO and AO indices). At the bay scale, we found no direct evidence that river discharge, which can supply the surface waters with nutrients and impact the vertical stratification, influenced the MIZ.

The relative contributions of the under-ice, MIZ, and fall blooms to total annual primary production remain to be established. Regionally tuned satellite-derived primary production models (Ardyna et al., 2013; Bélanger et al., 2013; Lee et al., 2015) could be used at least to assess the contribution of the MIZ and fall blooms. The combination of in situ observations, satellite monitoring (e.g., ocean color in synergy with sea-ice thermodynamic stages and/or albedo), and 3D ocean models coupled to biological models will be most promising to predict how all components of the HBS marine ecosystem will be impacted by climate change, including primary productivity.

The supplemental files for this article can be found as follows:

Figure S1. Ice-edge blooms between 1998 and 2018 in the Hudson Bay System. Ice-edge blooms detected by the maximum Chla during the 24 days after sea-ice retreat (first day of continuous SIC < 10%) (Perrette et al., 2011) between 1998 and 2018. Days of the year and their respective calendar date: 120 (Apr 30), 135 (May 15), 150 (May 30), 165 (Jun 14), 180 (Jun 29), 195 (Jul 14), 210 (Jul 29), 225 (Aug 13), and 240 (Aug 28).

Figure S2. Sea-ice retreat between 1998 and 2018 in the Hudson Bay System. Sea-ice retreat is detected when SIC (Comiso, 2000) reaches the threshold of 10% (Perrette et al., 2011).

Figure S3. Maps of four phytoplankton phenological categories throughout marginal ice zone between 1998 and 2018. Annual classification maps of MIZ phenology based on evolution of Chla following sea-ice retreat (t_R). The 4 categories are oligotrophic state or old under-ice bloom (blue); probable (recent) under-ice bloom (red); mesotrophic system where efficient nutrient replenishment is in place (green), bloom triggered in ice-free waters (orange); and bloom triggered under ice (cyan).

Figure S4. Winter air temperature anomalies in the Hudson Bay System. Winter anomalies of surface air temperature (⁠$Tair′$⁠) calculated using air temperatures from National Centers for Environmental Prediction/Atmospheric Research (NCEP/NCAR) Reanalysis Project as the difference between each winter (January, February, and March) average and its corresponding climatology and normalized by the standard deviation, using data from the 1948 to 2018 period.

The day of the year of Sea-ice retreat (t_R) and maximum chlorophyll-a concentration in the sea-ice edge zone produced in this research are made available through the BaySys data repository (open access).

This work is a contribution to the Natural Sciences and Engineering Council of Canada (NSERC) Collaborative Research and Development (CRD) project titled BaySys. Funding for this research was graciously provided by Manitoba Hydro, NSERC, Amundsen Science, and the Canada Research Chairs program. In addition, this research contributes to the ArcticNet Networks of Centers of Excellence and the Arctic Science Partnership (ASP). We thank the NSIDC, NASA-OBPG, and ESA-GlobColour for providing satellite data freely.

The authors would like to express their gratitude to Prof. Christopher Horvat, Brown University (Rhode Island/USA), because his review and directions were very important to improve this study.

Additional information: http://umanitoba.ca/faculties/environment/departments/ceos/research/BaySys.html.

LBDF is supported by the BaySys grant as well as UQAR and Québec-Ocean grants. This study was also supported by individual grants from ArcticNet and NSERC (355774-2009 and RGPIN-2014-03680 to S.B.).

The authors have no competing interests to declare.

Contributed to conception and design: LBF and SB.

Contributed to analysis and interpretation of data: all authors.

Drafted and/or revised the article: all authors.

How to cite this article: Barbedo, L, Bélanger, S, Tremblay, JE. 2020. Climate control of sea-ice edge phytoplankton blooms in the Hudson Bay System. Elem Sci Anth, 8: 1. DOI: https://doi.org/10.1525/elementa.039

Domain Editor-in-Chief: Jody W. Deming, University of Washington, WA, USA

Associate Editor: Kevin Arrigo, Stanford University, CA, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: BaySys