The Hudson Bay System (HBS), the world’s largest inland sea, has experienced disproportionate atmospheric warming and sea-ice decline relative to the whole Arctic Ocean during the last few decades. The establishment of almost continuous positive atmospheric air temperature anomalies since the late 1990s impacted its primary productivity and, consequently, the marine ecosystem. Here, four decades of archived satellite ocean color were analyzed together with sea-ice and climatic conditions to better understand the response of the HBS to climate forcing concerning phytoplankton dynamics. Using satellite-derived chlorophyll-a concentration [Chla], we examined the spatiotemporal variability of phytoplankton concentration with a focus on its phenology throughout the marginal ice zone. In recent years, phytoplankton phenology was dominated by two peaks of [Chla] during the ice-free period. The first peak occurs during the spring-to-summer transition and the second one happens in the fall, contrasting with the single bloom observed earlier (1978–1983). The ice-edge bloom, that is, the peak in [Chla] immediately found after the sea-ice retreat, showed substantial spatial and interannual variability. During the spring-to-summer transition, early sea-ice retreat resulted in ice-edge bloom intensification. In the northwest polynya, a marine wildlife hot spot, the correlation between climate indices, that is, the North Atlantic Oscillation and Arctic Oscillation (NAO/AO), and [Chla] indicated that the bloom responds to large-scale atmospheric circulation patterns in the North Hemisphere. The intensification of westerly winds caused by the strong polar vortex during positive NAO/AO phases favors the formation of the polynya, where ice production and export, brine rejection, and nutrient replenishment are more efficient. As a result, the winter climate preconditions the upper layer of the HBS for the subsequent development of ice-edge blooms. In the context of a decline in the NAO/AO strength related to Arctic warming, primary productivity is likely to decrease in the HBS and the northwest polynya in particular.

## 1. Introduction

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 km3y–1 and –287 km3y–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 × 106 km2). For example, Stroeve and Notz (2018) showed that Hudson Bay sea-ice cover in September has decreased at the rate of –1,046 km2 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.

## 2. Materials and methods

### 2.1. Satellite data

#### 2.1.1. Ocean color

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 (bbp) 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).

#### 2.1.2. SIC

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.

### 2.2. Climate index and reanalysis

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.

### 2.3. Data analysis

#### 2.3.1. Phytoplankton phenology during the ice-free season

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, Chlclim(x, y, t), for each 4.60-km pixel (x, y), as:

$ChlN(x,y,t)=Chlclim(x,y,t)−Chl‾(x,y)σChl(x,y)$
1

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).

#### 2.3.2. Phytoplankton phenology at the marginal ice zone

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, tR, 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 (tR).

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.

Figure 1.

Schematic of phytoplankton phenologies throughout marginal ice zones. The five possible phenologies are consistently low biomass indicating undetectable (old) under-ice bloom or an oligotrophic set up (I: blue line or dotted), probable (recent) under-ice bloom (II: red line), consistently high biomass indicating efficient nutrient replenishment (III: green line), bloom triggered in ice-free waters after the sea-ice retreat (VI: orange line), and bloom triggered under-ice that develop a peak just after the sea-ice retreat (V: cyan line). DOI: https://doi.org/10.1525/elementa.039.f1

Figure 1.

Schematic of phytoplankton phenologies throughout marginal ice zones. The five possible phenologies are consistently low biomass indicating undetectable (old) under-ice bloom or an oligotrophic set up (I: blue line or dotted), probable (recent) under-ice bloom (II: red line), consistently high biomass indicating efficient nutrient replenishment (III: green line), bloom triggered in ice-free waters after the sea-ice retreat (VI: orange line), and bloom triggered under-ice that develop a peak just after the sea-ice retreat (V: cyan line). DOI: https://doi.org/10.1525/elementa.039.f1

Table 1.

Remote sensing criteria defining the types of phytoplankton phenology throughout marginal ice zones. DOI: https://doi.org/10.1525/elementa.039.t1

Phenology TypeaCriteriab
I—Oligotrophic or old under-ice bloom $[Chla(tR)]
$[Chla(tR)]<[Chla(t1,2,3...)max]
II—Probable (recent) under-ice bloom $Ol<[Chla(tR)]
$[Chla(t1,2,3...)max]<[Chla(tR)]$
III—Mesotrophic/nutrients replenished $[Chla(tR)]≥B$
$[Chla(t1,2,3...)max]≥B$
$[Chla(t1,2,3...)min]≥B$
VI—Bloom triggered in ice-free waters $[Chla(tR)]<[Chla(t1,2,3...)max]$
$[Chla(tR)]
$[Chla(t1,2,3...)max]≥B$
V—Bloom triggered under-ice $[Chla(tR)]<[Chla(t1,2,3...)min]$
$[Chla(tR)]≥B$
$[Chla(t1,2,3...)min]
Phenology TypeaCriteriab
I—Oligotrophic or old under-ice bloom $[Chla(tR)]
$[Chla(tR)]<[Chla(t1,2,3...)max]
II—Probable (recent) under-ice bloom $Ol<[Chla(tR)]
$[Chla(t1,2,3...)max]<[Chla(tR)]$
III—Mesotrophic/nutrients replenished $[Chla(tR)]≥B$
$[Chla(t1,2,3...)max]≥B$
$[Chla(t1,2,3...)min]≥B$
VI—Bloom triggered in ice-free waters $[Chla(tR)]<[Chla(t1,2,3...)max]$
$[Chla(tR)]
$[Chla(t1,2,3...)max]≥B$
V—Bloom triggered under-ice $[Chla(tR)]<[Chla(t1,2,3...)min]$
$[Chla(tR)]≥B$
$[Chla(t1,2,3...)min]

aDefined based on the time evolution of Chla after the sea-ice retreat (with Chla in mg m–3).

bB is the bloom threshold (B), here set to 0.5 mg m–3; the oligotrophic threshold (Ol) was set to 0.2 mg m–3; tR is the day of sea-ice retreat defined as the day at which sea-ice concentration (SIC) is below 10% for at least 24 days; $t1,2,3...$ refers to weeks after the sea-ice retreat.

#### 2.3.3 Correlation analysis between MIZ Chla and environmental forcing

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 (tR), 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.

## 3. Results

### 3.1. Air temperature and phytoplankton phenology: CZCS versus MODIS era

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 (Tair) 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, Tair 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 Tair during most of the year.

Figure 2.

Surface air temperature anomalies and ocean color satellite coverage. Monthly anomalies of air temperature (Tair) spanning the period from 1948 to 2018 (negative anomalies in blue and positive in red bars). The temporal coverage of CZCS (1978–1983), MODIS (2002–2014), and Globcolour (1998–2018) Chla observations are marked. During the modern era of ocean color radiometry, positive values for Tair were often observed during the spring-to-summer transition. DOI: https://doi.org/10.1525/elementa.039.f2

Figure 2.

Surface air temperature anomalies and ocean color satellite coverage. Monthly anomalies of air temperature (Tair) spanning the period from 1948 to 2018 (negative anomalies in blue and positive in red bars). The temporal coverage of CZCS (1978–1983), MODIS (2002–2014), and Globcolour (1998–2018) Chla observations are marked. During the modern era of ocean color radiometry, positive values for Tair were often observed during the spring-to-summer transition. DOI: https://doi.org/10.1525/elementa.039.f2

The distinct Tair 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.

Figure 3.

Phytoplankton phenology changes between CZCS (1978–1983) and MODIS eras (2002–2014). Cluster-derived maps of phytoplankton phenological regimes for the (A) CZCS era (1978–1983) and the (B) MODIS era (2002–2014). (C) The plot shows the three types of phenology obtained using the K-means cluster analysis. (D) Open water season duration (number of days with sea-ice concentration lower than 10%) for each phenological domain in the CZCS era (blue boxplots) and MODIS era (yellow boxplots). The boxplots show the median (red lines), 25th and 75th percentiles at the hinges, and the whiskers extend to show ±1 standard deviation. DOI: https://doi.org/10.1525/elementa.039.f3

Figure 3.

Phytoplankton phenology changes between CZCS (1978–1983) and MODIS eras (2002–2014). Cluster-derived maps of phytoplankton phenological regimes for the (A) CZCS era (1978–1983) and the (B) MODIS era (2002–2014). (C) The plot shows the three types of phenology obtained using the K-means cluster analysis. (D) Open water season duration (number of days with sea-ice concentration lower than 10%) for each phenological domain in the CZCS era (blue boxplots) and MODIS era (yellow boxplots). The boxplots show the median (red lines), 25th and 75th percentiles at the hinges, and the whiskers extend to show ±1 standard deviation. DOI: https://doi.org/10.1525/elementa.039.f3

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.

Table 2.

Changes in the number of days (Δt) of open water season between the CZCS (1978–1983) and MODIS (2002–2014) eras, according to phenological cluster. DOI: https://doi.org/10.1525/elementa.039.t2

ClusteraMODIS SBMODIS FBMODIS DB
CZCS SB 39 24 34
CZCS FB 27 15 21
CZCS DB 24 18 23
ClusteraMODIS SBMODIS FBMODIS DB
CZCS SB 39 24 34
CZCS FB 27 15 21
CZCS DB 24 18 23

aSB indicates spring bloom; FB, fall bloom; and DB, double blooms.

### 3.2. Marginal ice zone blooms between 1998 and 2018

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.

Figure 4.

Climatologies of Chla, percentage of open water, and photosynthetic active radiation (PAR). (A) Chla climatology from Globcolour (GSM; merged product) and open water fraction (blue area) as a function of day of year (DOY). (B) $PAR(0+)$ climatology (Bélanger et al., 2013; Laliberté et al., 2016) (black lines) their respective standard deviation (shaded) observed between 2003 and 2010 in HBS. Sea-ice retreat was typically followed by a pelagic bloom during the spring-to-summer transition with a broad peak in early June (DOY 152) that was preceded by a smaller peak in February (DOY 120) and followed by a peak in the fall (DOY 243). DOI: https://doi.org/10.1525/elementa.039.f4

Figure 4.

Climatologies of Chla, percentage of open water, and photosynthetic active radiation (PAR). (A) Chla climatology from Globcolour (GSM; merged product) and open water fraction (blue area) as a function of day of year (DOY). (B) $PAR(0+)$ climatology (Bélanger et al., 2013; Laliberté et al., 2016) (black lines) their respective standard deviation (shaded) observed between 2003 and 2010 in HBS. Sea-ice retreat was typically followed by a pelagic bloom during the spring-to-summer transition with a broad peak in early June (DOY 152) that was preceded by a smaller peak in February (DOY 120) and followed by a peak in the fall (DOY 243). DOI: https://doi.org/10.1525/elementa.039.f4

To illustrate the recurrence of ice-edge blooms in the HBS, we plotted the frequency distribution of Chla (N approximately 106 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).

Figure 5.

Frequency distribution of Chla after the sea-ice retreat. Frequency distribution of Chla as a function of the number of days after the sea-ice retreat. The marginal ice zone (red box) is defined by the first 24 days after the sea-ice retreat (Perrette et al., 2011). The subsequent days define the open water season. Relations are extracted from annual time series of Chla for all pixels and all years between 1998 and 2018 in the HBS. DOI: https://doi.org/10.1525/elementa.039.f5

Figure 5.

Frequency distribution of Chla after the sea-ice retreat. Frequency distribution of Chla as a function of the number of days after the sea-ice retreat. The marginal ice zone (red box) is defined by the first 24 days after the sea-ice retreat (Perrette et al., 2011). The subsequent days define the open water season. Relations are extracted from annual time series of Chla for all pixels and all years between 1998 and 2018 in the HBS. DOI: https://doi.org/10.1525/elementa.039.f5

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 (tR) for the same period (1998–2018) based on annual maps shown in Figure S2. The central HB had the latest tR, ranging from DOY 180 to 210 (month of July), relatively low interannual variability (STD < 10). The highest (lowest) interannual variability of tR occurred in the coastal (offshore) areas of SW and NE Hudson Bay. Interestingly, the tR 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.

Figure 6.

Climatic maps of maximum Chla in the ice-edge zone and of sea-ice retreat timing. Average (A) and standard deviation (B) of the maximum Chla observed during the 24-day period after the sea-ice retreat. The statistics were calculated using the annual maps between 1998 and 2018 (n = 21 for each grid cell) shown in Figure S1. Similar statistics (C and D) were made for the DOY of sea-ice retreat (tR) shown in Figure S2. DOI: https://doi.org/10.1525/elementa.039.f6

Figure 6.

Climatic maps of maximum Chla in the ice-edge zone and of sea-ice retreat timing. Average (A) and standard deviation (B) of the maximum Chla observed during the 24-day period after the sea-ice retreat. The statistics were calculated using the annual maps between 1998 and 2018 (n = 21 for each grid cell) shown in Figure S1. Similar statistics (C and D) were made for the DOY of sea-ice retreat (tR) shown in Figure S2. DOI: https://doi.org/10.1525/elementa.039.f6

Table 3 presents the statistics and frequency of occurrence of the surface extension (in km2 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 km2 (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.

Table 3.

Annual extension of distinct trophic states in the ice-edge zone estimated using Chla thresholds and sea-ice retreat between 1998 and 2018 in the HBS. DOI: https://doi.org/10.1525/elementa.039.t3

Ice-edge Zone Trophic StateExtension in 104 km2 (%a) and Extreme Yearsb
CategoriesChla (mg m–3)MedianMax., yearMin., yearTrendc,dRMSEe
Oligotrophy <0.2 3.7 (4.6%) 18.2 (22.7%), 1999 0.5 (0.6%), 2008 –0.4 (–0.4%) 4.3
Moderate oligotrophy 0.2–0.5 28.8 (35.9%) 36.4 (45.4%), 2000 13.6 (17.0%), 2015 –0.2 (–0.3%) 4.9
Bloom 0.5–2.0 33.1 (41.3%) 45.0 (56.2%), 2008 22.3 (27.9%), 1999 +0.4 (+0.5%) 5.6
Major bloom >2.0 4.0 (5.0%) 7.2 (9.0%), 2015 1.8 (2.3%), 2010 0.0 (0.0%) 1.5
Sea-ice retreatf 183 194, 2004 174, 2006 +0.1 5.7
Ice-edge Zone Trophic StateExtension in 104 km2 (%a) and Extreme Yearsb
CategoriesChla (mg m–3)MedianMax., yearMin., yearTrendc,dRMSEe
Oligotrophy <0.2 3.7 (4.6%) 18.2 (22.7%), 1999 0.5 (0.6%), 2008 –0.4 (–0.4%) 4.3
Moderate oligotrophy 0.2–0.5 28.8 (35.9%) 36.4 (45.4%), 2000 13.6 (17.0%), 2015 –0.2 (–0.3%) 4.9
Bloom 0.5–2.0 33.1 (41.3%) 45.0 (56.2%), 2008 22.3 (27.9%), 1999 +0.4 (+0.5%) 5.6
Major bloom >2.0 4.0 (5.0%) 7.2 (9.0%), 2015 1.8 (2.3%), 2010 0.0 (0.0%) 1.5
Sea-ice retreatf 183 194, 2004 174, 2006 +0.1 5.7

a% of total area of Hudson Bay and Hudson Strait deeper than 50 m: 80 × 104 km2 (Amante and Eakins, 2009);

bExtreme events are reported according to year of maximum and minimum coverage;

cIn an interval of confidence of 95% (P < 0.05) for the extension of trophic states;

dIn an interval of confidence of 99% (P < 0.01) for the day of sea-ice retreat;

eRoot mean square error;

fFirst day of the year when sea-ice concentration is below 10%.

Figure 7.

Maximum Chla in ice-edge zone and ice-retreat timing during extreme years and the BaySys expedition. Annual map of maximum Chla following the sea-ice retreat (left panels) and the timing (DOY) of sea-ice retreat (right panels) for selected years (same as Figures S1 and S2). Top panels are for 2008 and 2000, which showed the largest and smallest ice-edge bloom extension, respectively, as reported in Table 2. Bottom panels are for 2015, the year showing the largest ice-edge bloom in the NW polynya, and 2018, the year of the BaySys expedition. DOI: https://doi.org/10.1525/elementa.039.f7

Figure 7.

Maximum Chla in ice-edge zone and ice-retreat timing during extreme years and the BaySys expedition. Annual map of maximum Chla following the sea-ice retreat (left panels) and the timing (DOY) of sea-ice retreat (right panels) for selected years (same as Figures S1 and S2). Top panels are for 2008 and 2000, which showed the largest and smallest ice-edge bloom extension, respectively, as reported in Table 2. Bottom panels are for 2015, the year showing the largest ice-edge bloom in the NW polynya, and 2018, the year of the BaySys expedition. DOI: https://doi.org/10.1525/elementa.039.f7

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.

Figure 8.

Frequency of occurrence of four trophic stages in the ice-edge zone. Maps of frequency occurrence of four categories of chlorophyll-a concentration in the ice-edge zone based on 21 years (1998–2018) of ocean color observations defining MIZ trophic conditions: (A) oligotrophy, (B) moderate oligotrophy, (C) bloom, and (D) major bloom. DOI: https://doi.org/10.1525/elementa.039.f8

Figure 8.

Frequency of occurrence of four trophic stages in the ice-edge zone. Maps of frequency occurrence of four categories of chlorophyll-a concentration in the ice-edge zone based on 21 years (1998–2018) of ocean color observations defining MIZ trophic conditions: (A) oligotrophy, (B) moderate oligotrophy, (C) bloom, and (D) major bloom. DOI: https://doi.org/10.1525/elementa.039.f8

### 3.3. Phytoplankton phenology at the MIZ

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.

Figure 9.

Frequency of occurrence of five phytoplankton phenological types throughout the marginal ice zone. Maps of frequency of occurrence for the five types of MIZ phenology between 1998 and 2018, as depicted in Figure 1: (A) Oligotrophy/old under-ice blooms, (B) probable (recent) under-ice bloom, (C) mesotrophic/nutrient replenished, (D) bloom triggered in ice-free waters, and (E) bloom triggered under-ice. DOI: https://doi.org/10.1525/elementa.039.f9

Figure 9.

Frequency of occurrence of five phytoplankton phenological types throughout the marginal ice zone. Maps of frequency of occurrence for the five types of MIZ phenology between 1998 and 2018, as depicted in Figure 1: (A) Oligotrophy/old under-ice blooms, (B) probable (recent) under-ice bloom, (C) mesotrophic/nutrient replenished, (D) bloom triggered in ice-free waters, and (E) bloom triggered under-ice. DOI: https://doi.org/10.1525/elementa.039.f9

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).

### 3.4. The magnitude of the ice-edge bloom and timing of sea-ice breakup

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 (tR). The frequency occurrence of tR 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.

Figure 10.

Influence of sea-ice retreat timing on Chla along the ice edge in Hudson Bay subsystems. Boxplots of maximum chlorophyll-a concentration ($[Chla]max$) in the marginal sea-ice zone (MIZ) in relation to day of the year of sea-ice retreat (tR) in each subregion of the Hudson Bay System. The central red line mark is the $[Chla]max$ median, the edges of the blue boxes are the 25th and 75th percentiles, the whiskers extend to the most extreme data points between ±0.26 σ considering a normal distribution and the percentage of valid pixels for each boxplot in relation to each subregion is marked above the respective upper whisker. The significant differences between these boxes were ensured by the ANOVA statistical test at a confidence interval of 99.9%. Chla above the threshold of 0.5 mg m–3 (green line) defines the marginal ice bloom occurrence (Perrette et al., 2011). DOI: https://doi.org/10.1525/elementa.039.f10

Figure 10.

Influence of sea-ice retreat timing on Chla along the ice edge in Hudson Bay subsystems. Boxplots of maximum chlorophyll-a concentration ($[Chla]max$) in the marginal sea-ice zone (MIZ) in relation to day of the year of sea-ice retreat (tR) in each subregion of the Hudson Bay System. The central red line mark is the $[Chla]max$ median, the edges of the blue boxes are the 25th and 75th percentiles, the whiskers extend to the most extreme data points between ±0.26 σ considering a normal distribution and the percentage of valid pixels for each boxplot in relation to each subregion is marked above the respective upper whisker. The significant differences between these boxes were ensured by the ANOVA statistical test at a confidence interval of 99.9%. Chla above the threshold of 0.5 mg m–3 (green line) defines the marginal ice bloom occurrence (Perrette et al., 2011). DOI: https://doi.org/10.1525/elementa.039.f10

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 tR 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 tR in most of the central HB.

Figure 11.

Relationships between the timing of sea-ice retreat and the maximum Chla along the ice edge. Linear regressions between tR and the maximum Chla in the ice-edge zone (i.e., $[Chla]IEZ(x,y)=slope(x,y)tR(x,y)+offset(x,y)$) were performed for each grid cell (x, y) using the annual maps shown in Figures S1 and S2 (n = 21 for each grid cell). The slope of the regressions is presented in A and its associated root mean square error (RMSE) in B. DOI: https://doi.org/10.1525/elementa.039.f11

Figure 11.

Relationships between the timing of sea-ice retreat and the maximum Chla along the ice edge. Linear regressions between tR and the maximum Chla in the ice-edge zone (i.e., $[Chla]IEZ(x,y)=slope(x,y)tR(x,y)+offset(x,y)$) were performed for each grid cell (x, y) using the annual maps shown in Figures S1 and S2 (n = 21 for each grid cell). The slope of the regressions is presented in A and its associated root mean square error (RMSE) in B. DOI: https://doi.org/10.1525/elementa.039.f11

### 3.5. The magnitude of the ice-edge bloom and atmospheric forcing

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.

Figure 12.

Relationships between anomalies of winter air temperature and maximum Chla in the marginal ice zone. Linear regressions between the winter air temperature anomaly ($Tair′$) and the maximum Chla in the ice-edge zone (i.e., $[Chla]IEZ(x,y)=slope(x,y)tR(x,y)+offset(x,y)$) were performed for each grid cell (x, y) using the annual map shown in Figures S1 (n = 21 for each grid cell). The slope of the regressions is presented in A and its associated root mean square error (RMSE) in B. Note that a single value of $Tair′$ (Figure S4) was used for the whole HBS due to the lack of spatial variability. DOI: https://doi.org/10.1525/elementa.039.f12

Figure 12.

Relationships between anomalies of winter air temperature and maximum Chla in the marginal ice zone. Linear regressions between the winter air temperature anomaly ($Tair′$) and the maximum Chla in the ice-edge zone (i.e., $[Chla]IEZ(x,y)=slope(x,y)tR(x,y)+offset(x,y)$) were performed for each grid cell (x, y) using the annual map shown in Figures S1 (n = 21 for each grid cell). The slope of the regressions is presented in A and its associated root mean square error (RMSE) in B. Note that a single value of $Tair′$ (Figure S4) was used for the whole HBS due to the lack of spatial variability. DOI: https://doi.org/10.1525/elementa.039.f12

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).

Figure 13.

Teleconnection between climatic indices (NAO/AO) and phytoplankton ice-edge blooms. Map of correlation coefficients (P > 95%) between Chla maxima in the ice-edge zone and AO (A) and NAO (B) between annual 1998 and 2018 generated using the Chla GSM algorithm from Globcolour and climatic indexes obtained from CPC/NOAA. Pearson’s correlation with t-test smaller than 95% confidence interval was removed. A 95% confidence interval was applied to the time series of climatic indexes AO and NAO and normalized mean values of annual Chla maxima in the ice-edge zone for the whole HBS (C) and in the Northwest Hudson Bay polynya (D). DOI: https://doi.org/10.1525/elementa.039.f13

Figure 13.

Teleconnection between climatic indices (NAO/AO) and phytoplankton ice-edge blooms. Map of correlation coefficients (P > 95%) between Chla maxima in the ice-edge zone and AO (A) and NAO (B) between annual 1998 and 2018 generated using the Chla GSM algorithm from Globcolour and climatic indexes obtained from CPC/NOAA. Pearson’s correlation with t-test smaller than 95% confidence interval was removed. A 95% confidence interval was applied to the time series of climatic indexes AO and NAO and normalized mean values of annual Chla maxima in the ice-edge zone for the whole HBS (C) and in the Northwest Hudson Bay polynya (D). DOI: https://doi.org/10.1525/elementa.039.f13

## 4. Discussion

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.

### 4.1. Satellite ocean color of the HBS

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/Kd), 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.

### 4.2. Winter air temperatures and ice retreat impact on ice-edge blooms

What role do winter air temperatures ($Tair′$) and sea-ice retreat (tR) 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, tR 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 tR 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.

### 4.3. Ice-edge phytoplankon blooms response to global and local forcing

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).

Figure 14.

Conceptual model of teleconnections of NAO/AO to phytoplankton dynamics in the NW HB Polynya. DOI: https://doi.org/10.1525/elementa.039.f14

Figure 14.

Conceptual model of teleconnections of NAO/AO to phytoplankton dynamics in the NW HB Polynya. DOI: https://doi.org/10.1525/elementa.039.f14

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.

### 4.4. Potential role of under-ice blooms

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.

## 5. Conclusions and perspectives

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.

## Supplemental files

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 (tR). 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.

## Data Accessibility Statement

The day of the year of Sea-ice retreat (tR) 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).

## Acknowledgments

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.

## Funding information

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.).

## Competing interests

The authors have no competing interests to declare.

## Author contributions

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.

## References

Amante
,
C
,
Eakins
,
BW
.
2009
. Etopo1 1 Arc-Minute global relief model: Procedures, data sources and analysis.
Boulder, CO
:
National Geophysical Data Center, NOAA
.
Ambaum
,
MHP
,
Hoskins
,
BJ
,
Stephenson
,
DB
.
2001
.
Arctic oscillation or North Atlantic oscillation?
J Clim
14
(
16
):
3495
3507
. DOI: http://dx.doi.org/10.1175/1520-0442(2001)014<3495:AOONAO>2.0.CO;2.
Andrews
,
J
,
Babb
,
D
,
Barber
,
DG
.
2017
.
Climate change and sea ice: Shipping in Hudson Bay, James Bay, Hudson Strait, and Foxe Basin
.
Elem Sci Anth
6
(
19
). DOI: http://dx.doi.org/10.1525/elementa.281.
Antoine
,
D
,
Morel
,
A
,
Gordon
,
HR
,
Banzon
,
VF
,
Evans
,
RH
.
2005
.
Bridging ocean color observations of the 1980s and 2000s in search of long-term trends
.
J Geophys Res Oceans
110
(
C6
). DOI: http://dx.doi.org/10.1029/2004JC002620.
Ardyna
,
M
,
Babin
,
M
,
Gosselin
,
M
,
Devred
,
E
,
Bélanger
,
S
,
Matsuoka
,
A
,
Tremblay
,
J-É
.
2013
.
Parameterization of vertical chlorophyll-a in the Arctic Ocean: Impact of the subsurface chlorophyll maximum on regional, seasonal, and annual primary production estimates
.
Biogeosci
10
:
4383
4404
. DOI: http://dx.doi.org/10.5194/bg-10-4383-2013.
Ardyna
,
M
,
Babin
,
M
,
Gosselin
,
M
,
Devred
,
E
,
Rainville
,
L
,
Tremblay
,
J-É
.
2014
.
Recent Arctic Ocean sea ice loss triggers novel fall phytoplankton blooms
.
Geophys Res Lett
41
(
17
):
6207
6212
. DOI: http://dx.doi.org/10.1002/2014GL061047.
Ardyna
,
M
,
Claustre
,
H
,
Sallée
,
JB
,
D’Ovidio
,
F
,
Gentili
,
B
,
van Dijken
,
G
,
D’Ortenzio
,
F
,
Arrigo
,
KR
.
2017
.
Delineating environmental control of phytoplankton biomass and phenology in the Southern Ocean
.
Geophys Res Lett
44
(
10
):
5016
5024
. DOI: http://dx.doi.org/10.1002/2016GL072428.
Ardyna
,
M
,
Gosselin
,
M
,
Michel
,
C
,
Poulin
,
M
,
Tremblay
,
.
2011
.
Environmental forcing of phytoplankton community structure and function in the Canadian High arctic: Contrasting oligotrophic and eutrophic regions
.
Mar Ecol Prog Ser
442
:
37
57
. DOI: http://dx.doi.org/10.3354/meps09378.
Arrigo
,
KR
,
Perovich
,
DK
,
Pickart
,
RS
,
Brown
,
ZW
,
van Dijken
,
GL
,
Bates
,
NR
,
Arrigo
,
KR
.
2014
.
Phytoplankton blooms beneath the sea ice in the Chukchi sea
.
Deep Sea Res
105
:
1
16
. DOI: http://dx.doi.org/10.1016/j.dsr2.2014.03.018.
Arrigo
,
KR
,
Perovich
,
DK
,
Pickart
,
RS
,
Brown
,
ZW
,
van Dijken
,
GL
,
Lowry
,
KE
,
Mills
,
MM
,
Palmer
,
MA
,
Balch
,
WM
,
Bahr
,
F
,
Bates
,
NR
,
Benitez-Nelson
,
C
,
Bowler
,
B
,
Brownlee
,
E
,
Ehn
,
JK
,
Frey
,
KE
,
Garley
,
R
,
Laney
,
SR
,
Lubelczyk
,
L
,
Mathis
,
J
,
Matsuoka
,
A
,
Mitchell
,
BG
,
Moore
,
GWK
,
Ortega-Retuerta
,
E
,
Pal
,
S
,
Polashenski
,
CM
,
Reynolds
,
RA
,
Schieber
,
B
,
Sosik
,
HM
,
Stephens
,
M
,
Swift
,
JH
.
2012
.
Massive phytoplankton blooms under Arctic Sea Ice
.
Science
336
(
6087
). DOI: http://dx.doi.org/10.1126/science.1215065.
Arrigo
,
KR
,
van Dijken
,
GL
.
2015
.
Continued increases in Arctic Ocean primary production
.
Prog Oceanogr
136
:
60
70
. DOI: http://dx.doi.org/10.1016/j.pocean.2015.05.002.
Assmy
,
P
,
Fernández-Méndez
,
M
,
Duarte
,
P
,
Meyer
,
A
,
Randelhoff
,
A
,
Mundy
,
CJ
,
Olsen
,
LM
,
Kauko
,
HM
,
Bailey
,
A
,
Chierici
,
M
,
Cohen
,
L
,
Doulgeris
,
AP
,
Ehn
,
JK
,
Fransson
,
A
,
Gerland
,
S
,
Hop
,
H
,
Hudson
,
SR
,
Hughes
,
N
,
Itkin
,
P
,
Johnsen
,
G
,
King
,
JA
,
Koch
,
BP
,
Koenig
,
Z
,
Kwasniewski
,
S
,
Laney
,
SR
,
Nicolaus
,
M
,
Pavlov
,
AK
,
Polashenski
,
CM
,
Provost
,
C
,
Rösel
,
A
,
Sandbu
,
M
,
Spreen
,
G
,
Smedsrud
,
LH
,
Sundfjord
,
A
,
,
T
,
Tatarek
,
A
,
Wiktor
,
J
,
Wagner
,
PM
,
Wold
,
A
,
Steen
,
H
,
Granskog
,
MA
.
2017
.
Leads in Arctic pack ice enable early phytoplankton blooms below snow-covered sea ice
.
Sci Rep
7
(
January
):
1
9
. DOI: http://dx.doi.org/10.1038/srep40850.
Babin
,
M
,
Bélanger
,
S
,
Ellingsen
,
I
,
Forest
,
A
,
Fouest
,
VL
,
Lacoura
,
T
,
Ardynaa
,
M
,
,
D
.
2015
.
Estimation of primary production in the Arctic Ocean using ocean colour remote sensing and coupled physical-biological models: Strengths, limitations and how they compare
.
Prog Oceanogr
139
. DOI: http://dx.doi.org/10.1016/j.pocean.2015.08.008.
Baldwin
,
MP
,
Dunkerton
,
TJ
.
1999
.
Propagation of the Arctic Oscillation from the stratosphere to the troposphere
.
J Geophys Res
104
(
D24
):
30937
30946
. DOI: http://dx.doi.org/10.1029/1999JD900445.
Barber
,
DG
,
Hop
,
H
,
Mundy
,
CJ
,
Else
,
B
,
Dmitrenko
,
IA
,
Tremblay
,
JE
,
Ehn
,
JK
,
Assmy
,
P
,
Daase
,
M
,
Candlish
,
LM
,
Rysgaard
,
S
.
2015
.
Selected physical, biological and biogeochemical implications of a rapidly changing Arctic Marginal Ice Zone
.
Prog Oceanogr
139
:
122
150
. DOI: http://dx.doi.org/10.1016/j.pocean.2015.09.003.
Barthélemy
,
A
,
Fichefet
,
T
,
Goosse
,
H
,
,
G
.
2015
.
Modeling the interplay between sea ice formation and the oceanic mixed layer: Limitations of simple brine rejection parameterizations
.
Ocean Model
86
:
141
152
. DOI: http://dx.doi.org/10.1016/j.ocemod.2014.12.009.
Beckers
,
JM
,
Rixen
,
M
.
2003
.
EOF calculations and data filling from incomplete oceanographic datasets
.
J Atmos Ocean Technol
20
(
12
):
1839
1856
Bélanger
,
S
,
Babin
,
M
,
Tremblay
,
.
2013
.
Increasing cloudiness in Arctic damps the increase in phytoplankton primary production due to sea ice receding
.
Biogeosci
10
(
6
):
4087
4101
. DOI: http://dx.doi.org/10.5194/bg-10-4087-2013.
Ben Mustapha
,
S
,
Bélanger
,
S
,
Larouche
,
P
.
2012
.
Evaluation of ocean color algorithms in the southeastern Beaufort Sea, Canadian Arctic: New parameterization using SeaWiFS, MODIS, and MERIS spectral bands
.
Can J Remote Sens
38
(
5
):
535
556
. DOI: http://dx.doi.org/10.5589/m12-045.
Brand
,
U
,
Came
,
RE
,
Affek
,
H
,
Azmy
,
K
,
Mooi
,
R
,
Layton
,
K
.
2014
.
Climate-forced change in Hudson Bay seawater composition and temperature, Arctic Canada
.
Chem Geol
388
:
78
86
. DOI: http://dx.doi.org/10.1016/j.chemgeo.2014.08.028.
Bricaud
,
A
,
Morel
,
A
,
Prieur
,
L
.
1981
.
Absorption by dissolved organic matter of the sea (yellow substance) in the UV and visible domains
.
Limnol Oceanogr
26
(
1
):
43
53
. DOI: http://dx.doi.org/10.4319/lo.1981.26.1.0043.
Brody
,
SR
,
Lozier
,
MS
,
Dunne
,
JP
.
2013
.
A comparison of methods to determine phytoplankton bloom initiation
.
J Geophys Res Oceans
118
(
5
):
2345
2357
. DOI: http://dx.doi.org/10.1002/jgrc.20167.
Burt
,
WJ
,
Thomas
,
H
,
Miller
,
LA
,
Granskog
,
MA
,
Papakyriakou
,
TN
,
Pengelly
,
L
.
2016
.
Inorganic carbon cycling and biogeochemical processes in an Arctic inland sea (Hudson Bay)
.
Biogeosci
13
:
4659
4671
. DOI: http://dx.doi.org/10.5194/bg-13-4659-2016.
Campbell
,
EC
,
Wilson
,
EA
,
Moore
,
GWK
,
Riser
,
SC
,
Brayton
,
CE
,
Mazloff
,
MR
,
Talley
,
LD
.
2019
.
Antarctic offshore polynyas linked to Southern Hemisphere climate anomalies
.
Nature
570
(
7761
):
319
325
. DOI: http://dx.doi.org/10.1038/s41586-019-1294-0.
Carmack
,
E
,
Barber
,
D
,
Christensen
,
J
,
Macdonald
,
R
,
Rudels
,
B
,
Sakshaug
,
E
.
2006
.
Climate variability and physical forcing of the food webs and the carbon budget on panarctic shelves
.
Prog Oceanogr
71
(
2–4
):
145
181
. DOI: http://dx.doi.org/10.1016/j.pocean.2006.10.005.
Chaudhuri
,
AH
,
Ponte
,
RM
,
Nguyen
,
AT
.
2014
.
A comparison of atmospheric reanalysis products for the Arctic Ocean and implications for uncertainties in air-sea fluxes
.
J Clim
27
(
14
):
5411
5421
. DOI: http://dx.doi.org/10.1175/JCLI-D-13-00424.1.
Comiso
,
J
,
Cavalieri
,
DJ
,
Parkinson
,
C
,
Gloersen
,
P
.
1997
.
Passive microwave algorithms for sea ice concentration: A comparison of two techniques
.
Remote Sens Environ
60
(
96
):
357
384
. DOI: http://dx.doi.org/10.1016/S0034-4257(96)00220-9.
Comiso
,
JC
.
2000
.
Bootstrap sea ice concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS, Version 2
. DOI: http://dx.doi.org/10.5067/J6JQLS9EJ5HU.
Cota
,
GF
,
Wang
,
J
,
Comiso
,
JC
.
2004
.
Transformation of global satellite chlorophyll retrievals with a regionally tuned algorithm
.
Remote Sens Environ
90
(
3
):
373
377
. DOI: http://dx.doi.org/10.1016/j.rse.2004.01.005.
D’Ortenzio
,
F
,
Antoine
,
D
,
Martinez
,
E
,
Ribera D’Alcal
,
M
.
2012
.
Phenological changes of oceanic phytoplankton in the 1980s and 2000s as revealed by remotely sensed oceancolor observations
.
Global Biogeochem Cycles
26
(
4
):
1
16
. DOI: http://dx.doi.org/10.1029/2011GB004269.
D’Ortenzio
,
F
,
D’Alcalà
,
MR
.
2009
.
On the trophic regimes of the Mediterranean Sea: A satellite analysis
.
Biogeosci Discuss
5
(
4
):
2959
2983
. DOI: http://dx.doi.org/10.5194/bgd-5-2959-2008.
Déry
,
SJ
,
Stieglitz
,
M
,
McKenna
,
EC
,
Wood
,
EF
.
2005
.
Characteristics and trends of river discharge into Hudson, James, and Ungava Bays, 1964–2000
.
J Clim
18
(
14
):
2540
2557
. DOI: http://dx.doi.org/10.1175/JCLI3440.1.
Déry
,
SJ
,
Wood
,
EF
.
2004
.
Teleconnection between the Arctic Oscillation and Hudson Bay river discharge
.
Geophys Res Lett
31
(
18
). DOI: http://dx.doi.org/10.1029/2004GL020729.
Ding
,
Q
,
Wallace
,
JM
,
Battisti
,
DS
,
Steig
,
EJ
,
Gallant
,
AJE
,
Kim
,
H-J
,
Geng
,
L
.
2014
.
Tropical forcing of the recent rapid Arctic warming in the northeastern Canada and Greenland
.
Nature
509
:
209
210
. DOI: http://dx.doi.org/10.1038/nature13260.
Ferguson
,
SH
,
Loseto
,
LL
,
Mallory
,
ML
.
2010
.
A little less Arctic: Top predators in the world’s largest northern inland sea, Hudson Bay
.
Springer
:
the Netherlands
. DOI: http://dx.doi.org/10.1007/978-90-481-9121-5.
Ferland
,
J
,
Gosselin
,
M
,
Starr
,
M
.
2011
.
Environmental control of summer primary production in the Hudson Bay system: The role of stratification
.
J Mar Syst
88
(
3
):
385
400
. DOI: http://dx.doi.org/10.1016/j.jmarsys.2011.03.015.
Francis
,
JA
,
Vavrus
,
SJ
.
2015
.
Evidence for a wavier jet stream in response to rapid Arctic warming
.
Environ Res Lett
10
(
1
). DOI: http://dx.doi.org/10.1088/1748-9326/10/1/014005.
Galbraith
,
PS
,
Larouche
,
P
.
2011
.
Reprint of “Sea-surface temperature in Hudson Bay and Hudson Strait in relation to air temperature and ice cover breakup, 1985–2009
.
J Mar Syst
88
(
3
):
463
475
. DOI: http://dx.doi.org/10.1016/j.jmarsys.2011.06.006.
Garver
,
SA
,
Siegel
,
DA
.
1997
.
Inherent optical property inversion of ocean color spectra and its biogeochemical interpretation 1. Time series from the Sargasso Sea
.
J Geophys Res
102
(
C8
). DOI: http://dx.doi.org/10.1029/96JC03243.
Granskog
,
MA
,
Kuzyk
,
ZZA
,
Azetsu-Scott
,
K
,
Macdonald
,
RW
.
2011
.
Distributions of runoff, sea-ice melt and brine using δ 18O and salinity data—A new view on freshwater cycling in Hudson Bay
.
J Mar Syst
88
(
3
):
362
374
. DOI: http://dx.doi.org/10.1016/j.jmarsys.2011.03.011.
Granskog
,
MA
,
Macdonald
,
RW
,
Kuzyk
,
ZZA
,
Senneville
,
S
,
Mundy
,
CJ
,
Barber
,
DG
,
Stern
,
GA
,
Saucier
,
F
.
2009
.
Coastal conduit in southwestern Hudson Bay (Canada) in summer: Rapid transit of freshwater and significant loss of colored dissolved organic matter
.
J Geophys Res Oceans
114
(
C8
). DOI: http://dx.doi.org/10.1029/2009JC005270.
Granskog
,
MA
,
Macdonald
,
RW
,
Mundy
,
CJ
,
Barber
,
DG
.
2007
.
Distribution, characteristics and potential impacts of chromophoric dissolved organic matter (CDOM) in Hudson Strait and Hudson Bay, Canada
.
Cont Shelf Res
27
(
15
):
2032
2050
. DOI: http://dx.doi.org/10.1016/j.csr.2007.05.001.
Guéguen
,
C
,
Granskog
,
MA
,
McCullough
,
G
,
Barber
,
DG
.
2011
.
Characterisation of colored dissolved organic matter in Hudson Bay and Hudson Strait using parallel factor analysis
.
J Mar Syst
88
(
3
):
423
433
. DOI: http://dx.doi.org/10.1016/j.jmarsys.2010.12.001.
Heikkila
,
M
,
Pospelova
,
V
,
Forest
,
A
,
Stern
,
GA
,
Fortier
,
L
,
Macdonald
,
RW
.
2016
.
Dinoflagellate cyst production over an annual cycle in seasonally ice-covered Hudson Bay
.
Mar Micropaleontol
125
:
1
24
. DOI: http://dx.doi.org/10.1016/j.marmicro.2016.02.005.
Hennig
,
C
.
2007
.
Cluster-wise assessment of cluster stability
.
Comput Stat Data Anal
52
(
1
):
258
271
. DOI: http://dx.doi.org/10.1016/j.csda.2006.11.025.
Hirawake
,
T
,
Shinmyo
,
K
,
Fujiwara
,
A
,
Saitoh
,
SI
.
2012
.
Satellite remote sensing of primary productivity in the Bering and Chukchi Seas using an absorption-based approach
.
ICES J Mar Sci
69
(
7
):
1194
1204
. DOI: http://dx.doi.org/10.1093/icesjms/fss111.
Hochheim
,
KP
,
Barber
,
DG
.
2010
.
Atmospheric forcing of sea ice in Hudson Bay during the fall period, 1980-2005
.
J Geophys Res Oceans
115
(
5
):
1
20
. DOI: http://dx.doi.org/10.1029/2009JC005334.
Hochheim
,
KP
,
Barber
,
DG
.
2014
.
An update on the ice climatology of the Hudson Bay system
.
Arct Antarct Alp Res
46
(
1
):
66
83
. DOI: http://dx.doi.org/10.1657/1938-4246-46.1.66.
Hochheim
,
KP
,
Lukovich
,
JV
,
Barber
,
DG
.
2011
.
Atmospheric forcing of sea ice in Hudson Bay during the spring period, 1980–2005
.
J Mar Syst
88
(
3
):
476
487
. DOI: http://dx.doi.org/10.1016/j.jmarsys.2011.05.003.
Hoover
,
C
,
Pitcher
,
T
,
Christensen
,
V
.
2013
.
Effects of hunting, fishing and climate change on the Hudson Bay marine ecosystem: I. Re-creating past changes 1970–2009
.
Ecol Model
264
:
130
142
. DOI: http://dx.doi.org/10.1016/j.ecolmodel.2013.02.005.
Horvat
,
C
,
Jones
,
DR
,
Iams
,
S
,
Schroeder
,
D
,
Flocco
,
D
,
Feltham
,
D
.
2017
.
The frequency and extent of sub-ice phytoplankton blooms in the Arctic Ocean
.
3
(
3
Hurrell
,
JW
,
Kushnir
,
Y
,
Ottersen
,
G
,
Visbeck
,
M
.
2003
.
An overview of the North Atlantic Oscillation
.
American Geophysical Union (AGU)
;
1
35
.
(The North Atlantic Oscillation: Climatic significance and environmental impact)
. DOI: http://dx.doi.org/10.1029/134GM01.
Hurrell
,
JW
,
Kushnir
,
Y
,
Visbeck
,
M
.
2001
.
The North Atlantic oscillation
.
Science
291
(
5504
):
603
605
. DOI: http://dx.doi.org/10.1126/science.1058761.
IOCCG
.
2015
.
Ocean colour remote sensing in polar seas
. In:
Babin
M
,
Arrigo
,
K
,
Bélanger
,
S
,
Forget
,
MH
, editors.
Dartmouth, NS, Canada (Reports of the International Ocean Colour Coordinating Group)
.
:
International Ocean Colour Coordinating Group
.
Available
: http://dx.doi.org/10.25607/OBP-107.
Jönsson
,
P
,
Eklundh
,
L
.
2004
.
TIMESAT–a program for analyzing time-series of satellite sensor data
.
Comput Geosci
30
(
8
):
833
845
. DOI: http://dx.doi.org/10.1016/j.cageo.2004.05.006.
Kahru
,
M
,
Brotas
,
V
,
Mazzano-Sarabia
,
M
,
Mitchell
,
BG
.
2011
.
Are phytoplankton blooms occurring earlier in the Arctic?
Glob Chang Biol
17
(
4
):
1733
1739
. DOI: http://dx.doi.org/10.1111/j.1365-2486.2010.02312.x.
Kahru
,
M
,
Lee
,
Z
,
Mitchell
,
BG
,
Nevison
,
CD
.
2016
.
Effects of sea ice cover on satellite-detected primary production in the Arctic Ocean
.
Biol Lett
12
(
11
). DOI: http://dx.doi.org/10.1098/rsbl.2016.0223.
Keller
,
WB
,
Paterson
,
AM
,
Rühland
,
KM
,
Blais
,
JM
.
2014
.
Introduction—Environmental change in the Hudson and James Bay region
.
Arct Antarct Alp Res
46
(
1
):
2
5
. DOI: http://dx.doi.org/10.1657/1938-4246-46.1.2.
Kowal
,
S
,
Gough
,
WA
,
Butler
,
K
.
2017
.
Temporal evolution of Hudson Bay sea ice (1971–2011)
.
Theor Appl Climatol
127
:
753
760
. DOI: http://dx.doi.org/10.1007/s00704-015-1666-9.
Kwok
,
R
.
2018
.
Arctic sea ice thickness, volume, and multiyear ice coverage: Losses and coupled variability (1958–2018)
.
Environ Res Lett
13
(
10
). DOI: http://dx.doi.org/10.1088/1748-9326/aae3ec.
Laliberté
,
J
,
Bélanger
,
S
,
Frouin
,
R
.
2016
.
Evaluation of satellite-based algorithms to estimate photosynthetically available radiation (PAR) reaching the ocean surface at high northern latitudes
.
Remote Sens Environ
184
:
199
211
. DOI: http://dx.doi.org/10.1016/j.rse.2016.06.014.
Landy
,
JC
,
Ehn
,
JK
,
Babb
,
DG
,
Thériault
,
N
,
Barber
,
DG
.
2017
.
Sea ice thickness in the Eastern Canadian Arctic: Hudson Bay Complex & Baffin Bay
.
Remote Sens Environ
200
(
August
):
281
294
. DOI: http://dx.doi.org/10.1016/j.rse.2017.08.019.
Lee
,
YJ
,
Matrai
,
PA
,
Friedrichs
,
MAM
,
Saba
,
VS
,
Antoine
,
D
,
Ardyna
,
M
,
Asanuma
,
I
,
Babin
,
M
,
Bélanger
,
S
,
Benoît-Gagné
,
M
,
Devred
,
E
,
Fernández-Méndez
,
M
,
Gentili
,
B
,
Hirawake
,
T
,
Kang
,
S-H
,
Kameda
,
T
,
Katlein
,
C
,
Lee
,
SH
,
Lee
,
Z
,
Mélin
,
F
,
Scardi
,
M
,
Smyth
,
TJ
,
Tang
,
S
,
Turpie
,
KR
,
Waters
,
KJ
,
Westberry
,
TK
.
2015
.
An assessment of phytoplankton primary productivity in the Arctic Ocean from satellite ocean color/in situ chlorophyll-a based models
.
J Geophys Res Oceans
120
(
9
):
6508
6541
. DOI: http://dx.doi.org/10.1002/2015JC011018.
Leu
,
E
,
Mundy
,
CJ
,
Assmy
,
P
,
Campbell
,
K
,
Gabrielsen
,
TM
,
Gosselin
,
M
,
Juul-Pedersen
,
T
,
,
R
.
2015
.
Arctic spring awakening—Steering principles behind the phenology of vernal ice algal blooms
.
Prog Oceanogr
139
:
151
170
. DOI: http://dx.doi.org/10.1016/j.pocean.2015.07.012.
Lewis
,
KM
,
Arrigo
,
KR
.
2020
.
Ocean color algorithms for estimating chlorophyll a, CDOM absorption, and particle backscattering in the Artic Ocean
.
J Geophy Res Oceans
125
(
6
):
1
23
. DOI: https://doi.org/10.1029/2019JC015706.
Leung
,
A
,
Gough
,
W
.
2016
.
Air mass distribution and the heterogeneity of the climate change signal in the Hudson Bay/Foxe Basin region, Canada
.
Theor Appl Climatol
125
:
583
592
. DOI: http://dx.doi.org/10.1007/s00704-015-1523-x.
Lowry
,
KE
,
van Dijken
,
GL
,
Arrigo
,
KR
.
2014
.
Evidence of under-ice phytoplankton blooms in the Chukchi Sea from 1998 to 2012
.
Deep Sea Res
105
:
105
117
. DOI: http://dx.doi.org/10.1016/j.dsr2.2014.03.013.
Marchese
,
C
,
Albouy
,
C
,
Tremblay
,
,
Dumont
,
D
,
D’Ortenzio
,
F
,
Vissault
,
S
,
Bélanger
,
S
.
2017
.
Changes in phytoplankton bloom phenology over the North Water (NOW) polynya: A response to changing environmental conditions
.
Polar Biol
40
:
1721
1737
. DOI: http://dx.doi.org/10.1007/s00300-017-2095-2.
Maritorena
,
S
,
D’Andon
,
OHF
,
Mangin
,
A
,
Siegel
,
DA
.
2010
.
Merged satellite ocean color data products using a bio-optical model: Characteristics, benefits and issues
.
Remote Sens Environ
114
(
8
):
1791
1804
. DOI: http://dx.doi.org/10.1016/j.rse.2010.04.002.
Maritorena
,
S
,
Siegel
,
DA
,
Peterson
,
AR
.
2002
.
Optimization of a semianalytical ocean color model for global-scale applications
.
Appl Opt
41
(
15
):
2705
2714
. DOI: http://dx.doi.org/10.1364/AO.41.002705.
Markus
,
T
,
Stroeve
,
JC
,
Miller
,
J
.
2009
.
Recent changes in Arctic sea ice melt onset, freezeup, and melt season length
.
J Geophys Res Oceans
114
(
12
):
1
14
. DOI: http://dx.doi.org/10.1029/2009JC005436.
Maynard
,
NG
,
Clark
,
DK
.
1987
.
Satellite color observations of spring blooming in Bering Sea shelf waters during the ice edge retreat in 1980
.
J Geophys Res
92
(
C7
). DOI: http://dx.doi.org/10.1029/JC092iC07p07127.
McClain
,
CR
.
2009
.
A decade of Satellite Ocean Color Observations
.
Ann Rev Mar Sci
1
:
19
42
. DOI: http://dx.doi.org/10.1146/annurev.marine.010908.163650.
Meleshko
,
VP
,
Johannessen
,
OM
,
Baidin
,
AV
,
Pavlova
,
T
,
Govorkova
,
VA
.
2016
.
Arctic amplification: Does it impact the polar jet stream?
Tellus A: Dynamic Meteorology and Oceanography
0870
. DOI: http://dx.doi.org/10.3402/tellusa.v68.32330.
Michel
,
C
,
Legendre
,
L
,
Demers
,
S
,
Therriault
,
JC
.
1988
.
Photoadaptation of sea-ice microalgae in springtime: Photosynthesis and carboxylating enzymes
.
Mar Ecol Prog Ser
50
:
177
185
. DOI: http://dx.doi.org/10.3354/meps050177.
Michel
,
C
,
Legendre
,
L
,
Therriault
,
JC
,
Demers
,
S
,
Vandevelde
,
T
.
1993
.
Springtime coupling between ice algal and phytoplankton assemblages in southeastern Hudson Bay, Canadian Arctic
.
Polar Biol
13
(
7
):
441
449
. DOI: http://dx.doi.org/10.1007/BF00233135.
Mitchell
,
BG
,
Brody
,
EA
,
Yeh
,
EN
,
Mcclain
,
C
,
Comiso
,
J
,
Maynard
,
NG
.
1991
.
Meridional zonation of the Barents Sea ecosystem inferred from satellite remote sensing and in situ bio-optical observations
.
Polar Res
10
(
1
):
147
162
. DOI: http://dx.doi.org/10.3402/polar.v10i1.6734.
Monti
,
D
,
Legendre
,
L
,
Therriault
,
JC
,
Demers
,
S
.
1996
.
Horizontal distribution of sea-ice microalgae: Environmental control and spatial processes (southeastern Hudson Bay, Canada)
.
Mar Ecol Prog Ser
133
(
1–3
):
229
240
. DOI: http://dx.doi.org/10.3354/meps133229.
Moritz
,
RE
,
Bitz
,
CM
,
Steig
,
EJ
.
2002
.
Dynamics of recent climate change in the Arctic
.
Science
297
(
5586
):
1497
1502
. DOI: http://dx.doi.org/10.1126/science.1076522.
Mundy
,
CJ
,
Gosselin
,
M
,
Ehn
,
J
,
Gratton
,
Y
,
Rossnagel
,
A
,
Barber
,
DG
,
Martin
,
J
,
Tremblay
,
J-É
,
Palmer
,
M
,
Arrigo
,
KR
,
Darnis
,
G
,
Fortier
,
L
,
Else
,
B
,
Papakyriakou
,
T
.
2009
.
Contribution of under-ice primary production to an ice-edge upwelling phytoplankton bloom in the Canadian Beaufort Sea
.
Geophys Res Lett
36
(
17
). DOI: http://dx.doi.org/10.1029/2009GL038837.
Mundy
,
CJ
,
Gosselin
,
M
,
Gratton
,
Y
,
Brown
,
K
,
Galindo
,
V
,
Campbell
,
K
,
Levasseur
,
M
,
Barber
,
D
,
Papakyriakou
,
T
,
Bélanger
,
S
.
2014
.
Role of environmental factors on phytoplankton bloom initiation under landfast sea ice in Resolute Passage, Canada
.
Mar Ecol Prog Ser
497
:
39
49
. DOI: http://dx.doi.org/10.3354/meps10587.
Mundy
,
CJ
,
Gosselin
,
M
,
Starr
,
M
,
Michelc
,
C
.
2010
.
Riverine export and the effects of circulation on dissolved organic carbon in the Hudson Bay system, Canada
.
Limnol Oceanogr
55
(
1
):
315
323
. DOI: http://dx.doi.org/10.4319/lo.2010.55.1.0315.
Nguyen
,
AT
,
Menemenlis
,
D
,
Kwok
,
R
.
2009
.
Improved modeling of the arctic halocline with a subgrid-scale brine rejection parameterization
.
J Geophys Res Oceans
114
(
11
):
1
12
. DOI: http://dx.doi.org/10.1029/2008JC005121.
O’Reilly
,
JE
,
Maritorena
,
S
,
Mitchell
,
BG
,
Siegel
,
DA
,
Carder
,
KL
,
Garver
,
SA
,
Kahru
,
M
,
McClain
,
C
.
1998
.
Ocean color chlorophyll algorithms for SeaWiFS
.
J Geophys Res
103
(
C11
):
24937
24953
. DOI: http://dx.doi.org/10.1029/98JC02160.
O’Reilly
,
JE
,
Maritorena
,
S
,
O’Brien
,
M
,
Siegel
,
D
,
Toole
,
D
,
Menzies
,
D
,
Smith
,
RC
,
Mueller
,
JL
,
Mitchell
,
BG
,
Kahru
,
M
,
Chavez
,
FP
,
Strutton
,
P
,
Cota
,
GF
,
Hooker
,
SB
,
McClain
,
CR
,
Carder
,
KL
,
Müller-Karger
,
F
,
Harding
,
L
,
Magnuson
,
A
,
Phinney
,
D
,
Moore
,
GF
,
Aiken
,
J
,
Arrigo
,
KR
,
Letelier
,
R
,
Culver
,
M
.
2000
.
SeaWIFS postlaunch calibration and validation analyses, Part 3
, in
Hooker
SB
,
Firestone
ER
eds.,
SeaWiFS postlaunch technical report series
., vol.
11
.
Greenbelt, Maryland
:
NASA
. https://oceancolor.gsfc.nasa.gov/docs/technical/seawifs_reports/postlaunch/post_vol11_abs/.
Accessed 01 June 2020
.
Overland
,
JE
,
Spillane
,
MC
,
Percival
,
DB
,
Wang
,
M
,
Mofjeld
,
HO
.
2004
.
Seasonal and regional variation of pan-arctic surface air temperature over the instrumental record
.
J Clim
17
:
3263
3282
. DOI: http://dx.doi.org/10.1175/1520-0442(2004)017<3263:SARVOP>2.0.CO;2.
Palmer
,
MA
,
Saenz
,
BT
,
Arrigo
,
KR
.
2014
.
Impacts of sea ice retreat, thinning, and melt-pond proliferation on the summer phytoplankton bloom in the Chukchi Sea, Arctic Ocean
.
Deep Sea Res
105
. DOI: http://dx.doi.org/10.1016/j.dsr2.2014.03.016.
Perrette
,
M
,
Yool
,
A
,
Quartly
,
GD
,
Popova
,
EE
.
2011
.
Near-ubiquity of ice-edge blooms in the Arctic
.
Biogeosci
8
(
2
):
515
524
. DOI: http://dx.doi.org/10.5194/bg-8-515-2011.
Piecuch
,
CG
,
Ponte
,
RM
.
2015
.
A wind-driven nonseasonal barotropic fluctuation of the Canadian inland seas
.
Ocean Sci
11
:
175
185
. DOI: http://dx.doi.org/10.5194/os-11-175-2015.
Pithan
,
F
,
Mauritsen
,
T
.
2014
.
Arctic amplification dominated by temperature feedbacks in contemporary climate models
.
Nat Geosci
7
(
3
):
181
184
. DOI: http://dx.doi.org/10.1038/ngeo2071.
Platt
,
T
,
White
,
GN
,
Zhai
,
L
,
Sathyendranath
,
S
,
Roy
,
S
.
2009
.
The phenology of phytoplankton blooms: Ecosystem indicators from remote sensing
.
Ecol Model
220
(
21
):
3057
3069
. DOI: http://dx.doi.org/10.1016/j.ecolmodel.2008.11.022.
Prinsenberg
,
S
.
1988
.
Ice-cover and ice-ridge contributions to the freshwater contents of Hudson Bay and Foxe Basin
.
Arctic
41
(
1
). DOI: http://dx.doi.org/10.14430/arctic1686.
Qian
,
M
,
Jones
,
C
,
Laprise
,
R
,
Caya
,
D
.
2008
.
The influences of NAO and the Hudson Bay sea-ice on the climate of eastern Canada
.
Clim Dyn
31
:
169
182
. DOI: http://dx.doi.org/10.1007/s00382-007-0343-9.
Renaut
,
S
,
Devred
,
E
,
Babin
,
M
.
2018
.
Northward expansion and intensification of phytoplankton growth during the early ice-free season in Arctic
.
Geophys Res Lett
45
(
19
). DOI: http://dx.doi.org/10.1029/2018GL078995.
Rogers
,
J
,
McHugh
,
M
.
2002
.
On the separability of the North Atlantic oscillation and Arctic oscillation
.
Clim Dyn
19
:
599
608
. DOI: http://dx.doi.org/10.1007/s00382-002-0247-7.
Runge
,
JA
,
Therriault
,
JC
,
Legendre
,
L
,
Ingram
,
RG
,
Demers
,
S
.
1991
.
Coupling between ice microalgal productivity and the pelagic, metazoan food web in southeastern Hudson Bay: A synthesis of results
.
Polar Res
10
(
2
):
325
338
. DOI: http://dx.doi.org/10.3402/polar.v10i2.6750.
Saucier
,
FJ
,
Senneville
,
S
,
Prinsenberg
,
S
,
Roy
,
F
,
Smith
,
G
,
Gachon
,
P
,
Caya
,
D
,
Laprise
,
R
.
2004
.
Modelling the sea ice-ocean seasonal cycle in Hudson Bay, Foxe Basin and Hudson Strait, Canada
.
Clim Dyn
23
(
3–4
):
303326
. DOI: http://dx.doi.org/10.1007/s00382-004-0445-6.
Screen
,
JA
,
Simmonds
,
I
.
2010
.
The central role of diminishing sea ice in recent Arctic temperature amplification
.
Nature
464
(
7293
):
1334
1337
. DOI: http://dx.doi.org/10.1038/nature09051.
Serreze
,
MC
,
Barrett
,
AP
,
Stroeve
,
JC
,
Kindig
,
DN
,
Holland
,
MM
.
2009
.
The emergence of surface-based Arctic amplification
.
The Cryosphere
3
(
1
):
11
19
. DOI: http://dx.doi.org/10.5194/tc-3-11-2009.
Sibert
,
V
,
Zakardjian
,
B
,
Gosselin
,
M
,
Starr
,
M
,
Senneville
,
S
,
LeClainche
,
Y
.
2011
.
3D bio-physical model of the sympagic and planktonic productions in the Hudson Bay System
.
J Mar Syst
88
(
3
):
401
422
. DOI: http://dx.doi.org/DOI10.1016/j.jmarsys.2011.03.014.
Sibert
,
V
,
Zakardjian
,
B
,
Saucier
,
F
,
Gosselin
,
M
,
Starr
,
M
,
Senneville
,
S
.
2010
.
Spatial and temporal variability of ice algal production in a 3D ice-ocean model of the Hudson Bay, Hudson Strait and Foxe Basin system
.
Polar Res
29
(
3
):
353
378
. DOI: http://dx.doi.org/10.1111/j.1751-8369.2010.00184.x.
Stewart
,
D
,
Lockhart
,
W
.
2005
.
An overview of the Hudson Bay marine ecosystem
.
(Canadian technical report of fisheries and aquatic sciences no. 2586)
.
Winnipeg, Manitoba
.
Available at
Government of Canada, Department of Fisheries and Oceans
website
Accessed 19 November 2020
.
St-Laurent
,
P
,
Straneo
,
F
,
Dumais
,
JF
,
Barber
,
D
.
2011
.
What is the fate of the river waters of Hudson Bay?
J Mar Syst
88
(
3
):
352
361
. DOI: http://dx.doi.org/10.1016/j.jmarsys.2011.02.004.
Stroeve
,
J
,
Notz
,
D
.
2018
.
Changing state of Arctic sea ice across all seasons
.
Environ Res Lett
13
(
10
Stroeve
,
JC
,
Markus
,
T
,
Boisvert
,
L
,
Miller
,
J
,
Barrett
,
A
.
2014
.
Changes in Arctic melt season and implications for sea ice loss
.
Geophys Res Lett
41
:
1216
1225
. DOI: http://dx.doi.org/10.1002/2013GL058951.
Stroeve
,
JC
,
Maslanik
,
J
,
Serreze
,
MC
,
Rigor
,
I
,
Meier
,
W
,
Fowler
,
C
.
2011
.
Sea ice response to an extreme negative phase of the Arctic Oscillation during winter 2009/2010
.
Geophys Res Lett
38
(
2
). DOI: http://dx.doi.org/10.1029/2010GL045662.
Sutherland
,
DA
,
Straneo
,
F
,
Lentz
,
SJ
,
Saint-Laurent
,
P
.
2011
.
Observations of fresh, anticyclonic eddies in the Hudson Strait outflow
.
J Mar Syst
88
(
3
):
375
384
. DOI: http://dx.doi.org/10.1016/j.jmarsys.2010.12.004.
Thompson
,
DW
,
Wallace
,
JM
.
1998
.
The Arctic oscillation signature in the wintertime geopotential height and temperature fields
.
Geophys Res Lett
25
(
9
):
1297
1300
. DOI: http://dx.doi.org/10.1029/98GL00950.
Tremblay
,
,
Lee
,
J
,
Gosselin
,
M
,
Belanger
,
S
.
2019
.
Nutrient dynamic and marine biological productivity in the greater Hudson Bay marine region
, in
Kuzyk
SZZ
,
Candlish
L
eds.,
An integrated regional impact study (IRIS) Arcticnet
.
University of Manitoba and ArcticNet
:
Chap. II Ecossys
:
225
244
.
Wallance
,
JM
.
2007
.
North Atlantic Oscillatiodannular mode: Two paradigms-one phenomenon
.
Q J R Meteorol Soc
126
:
791
805
. DOI: http://dx.doi.org/10.1002/qj.49712656402.
Wang
,
D
,
Wang
,
C
,
Yang
,
X
,
Lu
,
J
.
2005
.
Winter Northern Hemisphere surface air temperature variability associated with the Arctic Oscillation and North Atlantic Oscillation
.
Geophys Res Lett
32
(
16
):
1
4
. DOI: http://dx.doi.org/10.1029/2005GL022952.
Wang
,
J
,
Mysak
,
LA
,
Ingram
,
RG
.
1994
.
A three-dimensional numerical simulation of Hudson Bay summer ocean circulation: Topographic gyres, separations, and coastal jets
.
J Phys Oceanogr
24
(
12
):
2496
2514
. DOI: http://dx.doi.org/10.1175/1520-0485(1994)024<2496: ATDNSO>2.0.CO;2.
Webb
,
DJ
.
2014
.
On the tides and resonances of Hudson Bay and Hudson Strait
.
Ocean Sci
10
(
3
):
411
426
. DOI: http://dx.doi.org/10.5194/os-10-411-2014.
Winsor
,
P
,
Björk
,
G
.
2000
.
Polynya activity in the Arctic Ocean from 1958 to 1997
.
J Geophys Res Oceans
105
(
C4
):
8789
8803
. DOI: http://dx.doi.org/10.1029/1999jc900305.
Xi
,
H
,
Larouche
,
P
,
Michel
,
C
,
Tang
,
S
.
2015
.
Beam attenuation, scattering and backscattering of marine particles in relation to particle size distribution and composition in Hudson Bay (Canada)
.
J Geophys Res Oceans
120
(
5
):
3286
3300
. DOI: http://dx.doi.org/10.1002/2014JC010668.
Xi
,
H
,
Larouche
,
P
,
Tang
,
S
,
Michel
,
C
.
2013
.
Seasonal variability of light absorption properties and water optical constituents in Hudson Bay, Canada
.
J Geophys Res Oceans
118
(
6
):
3087
3102
. DOI: http://dx.doi.org/10.1002/jgrc.20237.
Xi
,
H
,
Larouche
,
P
,
Tang
,
S
,
Michel
,
C
.
2014
.
Characterization and variability of particle size distributions in Hudson Bay, Canada
.
J Geophys Res Oceans
119
(
6
):
3392
3406
. DOI: http://dx.doi.org/10.1002/2013JC009542.

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

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (CC-BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. See http://creativecommons.org/licenses/by/4.0/.