Atmospheric forcing and photo-acclimation of phytoplankton fall blooms in Hudson Bay

In polar seas, ongoing climatic changes have subjected phytoplankton dynamics to a new configuration of physical forcings that control incoming light, nutrient availability, and water column structure, with impacts on all marine trophic levels (Leu et al., 2011; Wassmann, 2011; Ardyna et al., 2014). Nowadays, sea-ice retreat is occurring earlier and recovering later (Markus et al., 2009; Stroeve et al., 2011; Stroeve et al., 2014), lengthening the open water season (Arrigo and van Dijken, 2015). Consequently, annual phytoplankton production rises because more light is available in the upper ocean layer (Bélanger et al., 2013; Arrigo and van Dijken, 2015). In turn, phytoplankton phenology, i.e., the seasonal succession of phytoplankton processes, impacts the carbon pump, pelagic ecosystems, the benthic realm, and valuable fish stocks (Platt et al., 2003; Friedland et al., 2008; Barber et al., 2015). The synchronization between phytoplankton bloom timing and organisms in early life and sensitive developmental stages (eggs or larval) is essential for species survival and recruitment (match/mismatch hypothesis: Cushing, 1990). Satellite observations have suggested that fall bloom occurrence has risen in high northern latitudes during the last decades (Ardyna et al., 2014). In low productive ecosystems like Hudson Bay (Tremblay et al., 2019), the world’s most extensive inland sea, an increase of phytoplankton abundance in the fall season can supply extra energy stock for polar species, increasing their overwintering potential.

Field and satellite observations in the middle- and high-northern latitudes have confirmed the importance of fall blooms. For instance, Rho and Whitledge (2007) showed that sediment traps recorded a high flux of organic carbon during the summer–fall transition in the Bering Sea. Sigler et al. (2014) reported that large crustacean zooplankton taxa respond to interannual variability in the time interval between spring and fall blooms in the Eastern Bering Sea. In general, short-term physical events, like storms, are assumed to affect the dynamics of fall blooms, as observed over the Nova Scotia Shelf (Greenan et al., 2004) and the sub-polar Japan/East Sea (Kim et al., 2007). Strength of wind-fetch and sea surface freezing delays have intensified phytoplankton blooms triggered by the transport of nutrients into the euphotic zone, as reported by Rho and Whitledge (2007) in the Bering Sea shelf, Fujiwara et al. (2018) in the Chukchi Sea, and Ardyna et al. (2014) and Castro de la Guardia et al. (2019) in the Arctic Ocean. Although recent studies have sought to understand the dynamics of phytoplankton fall blooms, information is still lacking about processes that trigger and maintain the effective productivity of fall blooms in north polar seas.

The Hudson Bay System (HBS) drainage basin alone equates to 18% of the total freshwater input of the Arctic Ocean (Déry et al., 2005; Granskog et al., 2011), while melting sea ice supplies twice this freshwater volume during a short melting period throughout the spring-to-summer transition (Mysak et al., 1996). The stratification resulting from this large amount of freshwater produces persistent oligotrophy in the summer (Harvey et al., 1997; Ferland et al., 2011; Lapoussiere et al., 2013; Xi et al., 2013; Tremblay et al., 2019). Ferland et al. (2011) reported low primary production rates in late summer (August and early September) dominated by small-sized phytoplankton adapted to oligotrophic and stratified conditions. In contrast, Xi et al. (2013) measured higher chlorophyll-a concentration ([chla]) and phytoplankton absorption coefficient in the fall (end of September to early October) compared to summer (July), suggesting an increase in primary productivity later during the open water season in the HBS. Recently, a systematic analysis of satellite-derived [chla] suggested that phytoplankton biomass systematically increases during the summer-to-fall transition (Barbedo et al., 2020), but the mechanisms explaining this phenomenon were not examined in detail.

Atmospheric forcing (convection- and wind-driven turbulence) trigger vertical mixing during the fall season, promoting the seeding of the upper-water layer with phytoplankton communities that developed at the subsurface chlorophyll maximum (SCM) during the summer and pumping nutrients in the euphotic zone. SCM phytoplankton is adapted to low light levels (photo-acclimation) with large cells having high intracellular chla content (Behrenfeld et al., 2002; Westberry et al., 2008; Cullen, 2015; Halsey and Jones, 2015; Marañón, 2015). Phytoplankton communities with high intracellular chla content in an environment of low photosynthetically active radiation (PAR) and excessive turbulence can be saturated quickly in terms of net carbon fixation (low maximum rate of photosynthesis, P_max; Huot et al., 2013) and can result in a satellite overestimation of primary production (Sathyendranath et al., 2020). To gain insights about physiological state of the fall bloom, we investigated how satellite-based approaches can help to address this question.

Phytoplankton community size structure and photo-physiological aspects (i.e., photo-acclimation and photoadaptation) are key aspects for a better understanding of pelagic ecosystems in terms of productivity, food web complexity, and the biological pump (Antoine et al., 2011; Brewin et al., 2014; Sathyendranath et al., 2020). Behrenfeld et al. (2005) used a satellite-derived particle backscattering coefficient (b_bp) as a proxy of phytoplankton carbon content (C_phy) allowing an estimation of the ratio of phytoplankton carbon to chlorophyll-a (C_phy:[chla]), a proxy of phytoplankton photoacclimation (Behrenfeld et al., 2002). Devred et al. (2011) derived phytoplankton size classes using chlorophyll-specific spectral absorption coefficient in the Canadian east coast and northwest Atlantic. Brewin et al. (2012) examined the chlorophyll-specific backscattering coefficient in relation with small (< 20 μm) and large (> 20 μm) phytoplankton size classes. In contrast, in Arctic waters where phytoplankton is not always the dominant optical component, spectral variations in remote sensing reflectance arise mainly from differences in the bio-optical environment in which specific communities are found (Reynolds and Stramski, 2019). To understand the phytoplankton photoacclimation state and its seasonal evolution, satellite-derived C_phy:[chla] was estimated with a method that considered the non-algal contribution to b_bp (Bellacicco et al., 2019).

Our study focuses on the physical mechanisms responsible for the fall bloom phenomenon as well as potential photo-physiological adaptations of phytoplankton to late summer/fall conditions. Specifically, the objectives were to i) understand the role of atmospheric forcings and water column structure on fall bloom dynamics, and ii) assess the physiological state of phytoplankton communities blooming in the fall season using satellite-derived optical proxies.

To investigate the role of atmospheric and oceanographic processes in triggering and maintaining phytoplankton fall blooms, we employed in synergy the outputs of an ice-ocean dynamic model, climatic reanalysis, the satellite ocean color constellation, and in situ bio-optic and radiometry. We discuss the seasonality of atmospheric and oceanographic interplays, water column evolution, and their effects on phytoplankton dynamics in Hudson Bay. In Southern Hudson Bay, inter-annual and seasonal variability of freshwater input results in remarkable shifts on stratification drivers, i.e., salinity and temperature vertical gradients. To better understand the counterpart of atmospheric forcing on the water column and phytoplankton dynamics, we focused on South Hudson Bay where more pronounced ice-edge and fall blooms were observed ([chla] > 2 mg m⁻³), and where riverine input can contribute to the upper-layer stratification during the melting season.

To enlighten phytoplankton dynamics in the fall season, we compared two decades (1998–2018) of ocean color satellite observations with bio-optical conditions collected in situ during a recent Canadian expedition throughout marginal ice zones during the spring-to-summer transition in Hudson Bay (2018), when phytoplankton is expected to bloom. Therefore, spring in situ observations can support the satellite ocean color observations to fill in the lack of information in fall. We analyzed in detail both satellite-derived [chla] and b_bp to explore the potential role of photo-acclimation on phytoplankton and its productivity during the fall season.

Vertical profiles of potential temperature and salinity were modeled using version 3.4 of the Nucleus for European Models of the Ocean (NEMO) general circulation model. NEMO is a coupled ice-ocean dynamic model. The experiments used the Arctic and Northern Hemisphere Atlantic configuration run at 1/4° resolution (ANHA4) and 5-day averaged output. The vertical resolution was 50 unevenly spaced levels, with the majority of the resolution in the upper column to capture the mixed layer and thermocline processes. The model inputs for atmospheric fields are taken from the Canadian Meteorological Centre’s Global Deterministic Prediction System Re-Forecasts (CGRF; Smith et al., 2014), and the coupled ocean-ice Coordinated Ocean Research Experiments (CORE) bulk formula was used to compute fluxes between the atmosphere and the ocean and to assimilate wind stress into the momentum equations (Large and Yeager, 2004; 2009). River runoff within HBS was provided by calibrated output from the HYdrological Predictions for the Environment (HYPE; Lindström et al., 2010; Andersson et al., 2015; Gelfan et al., 2017). This runoff dataset has not been integrated with observations but included the improvements to the model, such as a better representation of La Grande Rivière discharge and Nelson River regulated. Outside HBS, runoff was compiled by Dai and Trenberth (2002), with the addition of Greenland melt by Bamber et al. (2012). Further details on the model setup and evaluation can be found in Ridenour et al. (2019).

Heat fluxes (in W m⁻²), i.e., the difference in short- and long-wave radiation, latent and sensible heat between the atmosphere and the ocean, were computed using the bulk formulas developed by Large and Yeager (2004; 2009). The calculations considered near-surface atmospheric state (wind, potential temperature, specific humidity, and density), and the ocean state (ocean surface current and temperature). Briefly, sea-surface currents and sea-surface temperature (SST) were calculated from NEMO-ANHA4 model simulations. Atmospheric potential temperature, humidity, wind (10 m), precipitation, short and long-wave radiation were obtained from CGRF, which is assumed suitable to parameterize physical processes of polar seas and has a finer spatial-temporal scale (i.e., 33 km daily) than other current available global reanalysis (Smith et al., 2014).

Values for [chla] (in mg m⁻³), the particle backscattering coefficient at 443 nm (b_bp(443), in m⁻¹), and the diffuse attenuation coefficient at 490 nm (k_d(490), in m⁻¹) were downloaded from the Globcolour Project (http://www.globcolour.info/). Globcolour Level-3 [chla] (European Space Agency, 2020) and b_bp(443) were estimated using the semi-analytical algorithm of Garver-Siegel-Maritorenna (GSM: Garver and Siegel, 1997; Maritorena et al., 2002). We performed an in situ evaluation of the GSM algorithm to ensure the robustness of our choice of for [chla] (Figure S1). The value for b_bp(555) was calculated from b_bp(443) using a power law function with a spectral dependency of -1.0337 (Maritorena et al., 2002). The value for k_d(490) was estimated using the quasi-analytical algorithm of Lee et al. (2005). The Globcolour Level-3 (i.e., binned and mapped) merged products have temporal and spatial resolutions of 1 day and 4.63 km, respectively, and cover the period of 1998 to 2018. The merged products improve the spatial-temporal coverage, diminishing gaps due to cloud cover and sea ice (Maritorena et al., 2010). The binning methodology combines radiometry from ocean color if the sensors are available, which include the Sea Wide Field-of-view Sensor (SeaWiFS: 1998–2010), Moderate Resolution Imaging Spectroradiometer (MODIS: 2002–2018) aboard the Aqua satellite, Medium-Resolution Imaging Spectrometer (MERIS: 2002–2011), Visible Infrared Imaging Radiometer Suite (VIIRS: 2012–2018), and Ocean and Land Colour Instrument (OLCI: 2016–present). Coastal areas with depths shallower than 50 m were removed from analysis to avoid turbidity and bottom effect. Satellite data were restricted to complex water non-scatter dominated and case water-1 (Lee and Hu, 2006) to avoid the influence of non-algal particle backscattering on estimation of C_phy.

We evaluated the effect of atmospheric forcing on the water column structure and phytoplankton dynamics in spatial bin areas of 100 × 100 ocean color pixels (approximately 4 km) centered on 84°W 57°N in southern Hudson Bay. Time series of heat flux, wind speed, light level, and [chla] were combined with water column profiles. The first optical depth (z_OC: in m), the water column layer inside the ocean color satellites range, and the euphotic zone depth (z_eu: in m), defined as the depth at which PAR reduces to 1% of the incoming surface value, were estimated using k_d(490), i.e., $zeu=−4.5/kd(490)$ and $zOC=−1/kd(490)$⁠, respectively. For that, NEMO simulations were used to produce a high resolution, basin-scale hindcast of ocean vertical structure of potential temperature, salinity, stratification, and mixed layer depth. The mixed layer thickness was calculated based on sea water density (ρ), i.e., the depth of a threshold value of density (⁠$Δρ=0.01 kg m−3$⁠) from a near-surface value (de Boyer Montégut et al., 2004). Fluid stratification was quantified using the buoyancy frequency (N²), known as the Brunt-Vaissala frequency, i.e., the frequency required to break up the stratification (Equation 1):

where g is the gravity, ρ is the density, p is the pressure, and z is depth. The heat-flux inversion, i.e., day of the year when atmosphere began to cool the ocean, was detected by the first day in which the heat-flux sign was negative between June and November.

To indicate the ice-edge zone, the sea-ice retreat, t_R, was defined as the day at which sea-ice concentration (SIC) is below 10% for at least 24 days (Barbedo et al., 2020). SIC was obtained from the National Snow and Ice Data Center (NSIDC), which is based on daily multichannel passive microwave radiometry sensors clustered using the Bootstrap algorithm at 25-km resolution (Comiso et al., 1997; Comiso, 2000).

In situ bio-optical measurements were performed aboard the Canadian Icebreaker NCGC Amundsen during the BaySys expedition in Hudson Bay in June 2018 (Ahmed et al., 2021). The sampling was distributed over recently opened waters in western Hudson Bay in the marginal ice zone (MIZ), i.e., the area along the edge of the ice pack that is affected by open ocean processes (Barber et al., 2015), in leads surrounded by landfast ice loading/releasing sediment in the south, and along a transect in the Nelson River plume (Figure 1). The field measurements included in situ particle backscattering, pigment analysis and in-water radiometry, which are used to support our interpretation of the satellite-derived phytoplankton photoacclimation state described in the next section.

The total spectral backscattering coefficient, $bb(λ)$⁠, is the amount of light scattered by water molecules or suspended particles into the hemisphere from which light has originated. Total $bb(λ)$ can be further divided in terms of additive components of pure water (b_bw) and particles (b_bp). The latter is a good proxy of particulate organic carbon (POC) in the open ocean (Stramski et al., 1999). Total $bb(λ)$ was calculated from the volume scattering function (β: in sr⁻¹) measured at a fixed scattering angle (⁠$θb∘$⁠) in the backward direction relative to the laser beam using an ECO BB9 (Sea-Bird Scientific). The instrument was calibrated by the manufacturer before the field campaign. The BB9 operated at nine wavelengths (412, 440, 488, 510, 532, 595, 660, 675, and 715 nm) with $β(θb∘=117∘)$⁠.

The measurements were performed at 14 stations in Hudson Bay. The data were processed using an open source R package (DOI: https://github.com/belasi01/Riops; see also Bélanger et al., 2017) following the manufacturer’s protocol (Twardowski et al., 2007). Briefly, the spectral backscattering coefficient was obtained from the measured $β(θb∘,λ)$ through integration over the backward hemisphere based on a phase function according to the Mie theory (wavelength-selective scattering) using the relation (Equation 2):

where χ_p, a spectrally constant value of 1.1 (Boss and Pegau, 2001), represents the relationship between $β(θb∘)$ and b_bp. The spectral contribution of pure seawater to scattering (β_sw) at $θb∘$ was computed using temperature and salinity at the same time and depth with a CTD sensor (Zhang et al., 2009). The dependence on absorption, a, caused by attenuation along with the sensor and light source path was corrected based on simultaneous measurements of total absorption by an absorption and attenuation spectral (ac-s: Sea-Bird Scientific) sensor according to recommendations of Doxaran et al. (2016). Profiles of $bbp(λ,z)$ were binned at 0.1-m depth intervals after applying a local polynomial regression smoothing function (Cleveland et al., 1992); b_bp(555) was calculated at each depth by fitting a power law function to $bbp(λ,z)$⁠. Finally, b_bp at 555 nm was selected for further analysis because pigment absorption exhibits weak contribution in this spectral region, as reported by Reynolds et al. (2016).

Chlorophyll-a concentration from discrete samples were measured by the high-performance liquid chromatography (HPLC) technique (⁠$[chla]HPLC$⁠: in mg m⁻³). Seawater discrete samples of 1–3 L were collected using Niskin bottles and filtered onto Whatman GF/F filters (25 mm in diameter) under a gentle vacuum (3.98 Hg). The filters were placed immediately in 2-mL cryovials and flash-frozen in liquid nitrogen, then stored at −80°C until analysis. After the cruise, filters were extracted in 100% methanol at −20°C, disrupted by sonication, and clarified later by vacuum filtration through Whatman GF/F filters. Extraction time lasted 2 hours, and HPLC analysis was carried out the same day (Ras et al., 2008; Tran et al., 2013; Robinson et al., 2018). More details on the analytical methods can be found in Matthes et al. (2021) for the BaySys, and Ras et al. (2008) for the GreenEdge data sets. In this study, $[chla]HPLC$ corresponds to the sum of chla and chlorophyllide-a.

Downwelling in-water irradiance profiles, $Ed(z,λ)$⁠, were measured by Compact-Optical Profiling System (C-OPS; Hooker et al., 2013) built by Bio-Spherical Instruments Inc. (San Diego, California, USA). C-OPS has 19 bands distributed between the spectral range of 320 to 865 nm measuring the planar downwelling irradiance coupled with auxiliary sensors for water temperature, pressure, pitch and roll. A similar reference sensor in the ship deck measured above-water surface downwelling irradiance $Ed(0+,λ)$⁠. An above-water reference sensor was coupled to the Bioshade to measure the ratio of direct sun to diffuse sky irradiation (⁠$Edsun/Edsky$⁠; Morrow et al., 2010), which is used for instrument self-shadow correction.

We profiled at least three consecutive downcasts between the sea surface and a variable water depth corresponding to the approximate 1% light level. In-water radiometric measurements were normalized to account for the illumination variations during each cast. Underwater casts trespassing a vertical tilt of 5° were removed from the analysis. In-water radiometry was processed using the R-package COPS available in DOI: https://github.com/belasi01/Cops. A detailed discussion about the C-OPS instrument and the processing protocols can be found elsewhere (Mueller et al., 2003; Morrow et al., 2010; Antoine et al., 2013; Hooker et al., 2013; Bélanger et al., 2017).

The daily average of photosynthetic active radiation (⁠$PAR¯$⁠) was computed after applying the factor $cos(θs(t))/cos(θoncast)$⁠, where θ_s is the solar zenith range in the daytime (t) and θ_oncast is the solar zenith on cast of the light profile (Reda and Andreas, 2004). PAR (in μmol photons m⁻²s⁻¹) was calculated from trapezoidal integration of $Ed(λ,z)$ in the visible spectrum from 400 to 700 nm (Equation 3):

where ℏ is the plank constant (6.623×10⁻³⁴ J s⁻¹), c is the constant for the light speed (299,792,458 m s⁻¹), and N is Avogadro’s number (6.022×10²³). $Ed(λ,z)$ was measured in μW m⁻² nm⁻¹. The conversion from mol to μmol used a factor of 1×10⁶.

The spectral coefficient for diffuse downwelling irradiation attenuation, $kd(λ)$⁠, at a geometric depth, defined as (Equation 4):

was determined from in-water $Ed(z,λ)$ profiles following Mueller et al. (2003). The value for $kd(λ)$ was computed as the local slope of $lnEd(z)$ using a linear regression fit, and assuming a constant k_d within the interval $(zm−Δz)≤z≤(zm+Δz)$ where z_m is the center depth (Equation 5):

As proposed by Lee et al. (2005), $kd(z)$ in the euphotic zone may be replaced by the integrated value of k_d from surface to a depth where E_d is reduced to approximately 1% or 10% of its surface value, i.e., k_d value in the layer between $Ed(0−)$ and 1 or 10% of $Ed(0−)$⁠. Here, $kd(Ed10%)$ was denominated as k_d because this layer contributes most to the water column photosynthesis and signals measured by remote sensors.

Nelson River and landfast ice introduce considerable quantities of minerogenic particles into Hudson Bay. To reduce the contribution to non-chlorophyllous particles that efficiently backscatter light, our analysis was restricted to the MIZ light environment. Water masses were classified using a cluster analysis applied to $kd(λ)$⁠, which is closely related to inherent optical properties (IOPs) and [chla] (Solonenko and Mobley, 2015).

As mentioned in the Introduction, ocean optical properties can be used as proxies to assess phytoplankton photoacclimation or size structure proxies (Huot et al., 2008; Devred et al., 2011; Brewin et al., 2012; Marañón, 2015; Fujiwara et al., 2018; Reynolds and Stramski, 2019). Here, we examined the potential of satellite-derived phytoplankton photo-acclimation proxies to better understand the state of the phytoplankton community characterizing the fall blooms. Our assessment is based on [chla] and b_bp, which are important proxies for phytoplankton dynamics and pelagic ecosystems. While [chla] is a direct proxy for phytoplankton, b_bp may be influenced by a wide range of non-chlorophyllous particles such as bacteria, micrograzers, heterotrophic nano-flagellates, ciliates, and viruses; organic particles of detrital origin such as fecal pellets and cell debris (non-living organic detritus derived from the breakdown of micro-organisms); mineral particles of both biogenic (e.g., calcite liths and shells) and terrestrial origins (e.g., clays and sand); bubbles; and plastics (Vaillancourt et al., 2004; Bellacicco et al., 2019). Total b_bp therefore may be partitioned into two additive components:

where $bbp[chla]$ is the backscattering component that co-varies with phytoplankton biomass and $bbpk$⁠, the backscattering background that is due to non-algal particles (NAP) that do not co-vary with phytoplankton biomass.

To support our interpretations of b_bp as a phytoplankton proxy, we first needed to estimate $bbpk$⁠. Behrenfeld et al. (2005) for example, estimated $bbpk$ as the intercept of the least-squares regression analysis of the [chla] versus b_bp relationship. The backscattering background was considered as constant over the global ocean and fixed to a value of 0.00035 m^-1. More recently, however, Bellacicco et al. (2016) demonstrated that $bbpk$was not constant in either space or time scales in the Mediterranean Sea. Similar conclusions were drawn from an extensive analysis of global satellite ocean color data and in situ bio-ARGO data by Bellacicco et al. (2018; 2019), suggesting that $bbpk$ variability must be considered in order to assess phytoplankton-related backscattering (or carbon content assuming a scaling factor, e.g., Behrenfeld et al., 2005). Therefore, we computed $bbpk$ at 555 nm on a pixel-to-pixel basis using all daily satellite observations available for the 1998 to 2018 period for each month of open waters (May to October).

Values of $bbpk$ were calculated as the background value of the relationship between b_bp and [chla] by fitting three methods: i) the linear model proposed by Behrenfeld et al. (2005); ii) the model presented by Bellacicco et al. (2019), which was first proposed by Brewin et al. (2012) to account non-linearity driven by phytoplankton size fractions; and iii) an exponential model that we propose with an offset to account for $bbpk$ (Equation 7):

where a and b are fitting parameters. The three fitting procedures were first tested on our in situ observations of b_bp and $[chla]HPLC$ described above, as well as another data set from the Canadian Arctic (i.e., from Baffin Bay; results not shown).

Phytoplankton physiological adjustments to environmental conditions (e.g., light, nutrient status, taxonomy, stratification, and grazing) result in changes in the ratio of phytoplankton-carbon biomass (C_phy: in mg C m⁻³) to [chla] (Geider, 1987; Behrenfeld et al., 2005; Westberry et al., 2008; Sathyendranath et al., 2020). Therefore, assuming that phytoplankton dominates the covarying part of b_bp, as proposed by Behrenfeld et al. (2005), we estimated C_phy by multiplying $bbp[chla]$ by a fixed scaling factor (SF: 13,000 mg C m⁻²) (Equation 8):

Satellite ocean color monthly climatologies (1998–2018) of [chla] revealed the seasonal and spatial variability of phytoplankton in the HBS (Figure 2). In May, high [chla] (> 1 mg m⁻³) in open waters was a manifestation of ice-edge bloom. In June, concentrations remained high in most of the HBS also due to late ice-edge blooms, except in the central HB where oligotrophy was already established. The lowest [chla] were observed in August, except in coastal zones influenced by riverine inputs or where strong mixing processes (e.g., tidal mixing) operated year-long (i.e., in Hudson Strait). The trophic conditions revealed by satellite-derived [chla] are based on Barbedo et al. (2020) and Perrette et al. (2011) who defined the ice-edge bloom with [chla] greater than a threshold of 0.5 mg m⁻³ just after the ice retreat. The term oligotrophy is adopted for low [chla], as observed by Ferland et al. (2011) in summer in the HBS. Concentrations clearly began to increase again in September almost everywhere and continued to rise in October, including in the central HBS where oligotrophic conditions relaxed ([chla] increased from < 0.2 mg m⁻³ to > 0.5 mg m⁻³). We further examine in the next sections the physical processes potentially involved in explaining the late summer/fall increase in [chla] in the southern part of the oligotrophic area.

Figure 3 shows time-series of [chla], b_bp(555), wind speed, and heat flux between 2002 and 2009 in the south Hudson Bay (i.e., red box on Figure 2 centered at location 84°W and 57°N). High [chla] during the spring-to-summer transition (June–July) was related to ice-edge blooms and, as a consequence, the stratification onset following the ice breakup. A reduction in [chla] was abrupt after the ice-edge bloom peak, following a decline in the ocean heating rate until south Hudson Bay reached a surface oligotrophic state (i.e., lowest [chla]). Chlorophyll-a concentration trends (⁠$∂ [chla]∂t−1$⁠) indicated a systematic increase in [chla] between summer, with maximum oligotrophic state, and fall co-occurring with a seasonal process of atmospheric warming ocean decay (i.e., heat flux losing intensity). These trends were obtained by linear fit in the interval of confidence of 99% in 2003 (0.08 mg m⁻³ month⁻¹), 2005 (0.73 mg m⁻³ month⁻¹), 2006 (0.11 mg m⁻³ month⁻¹), 2007 (0.21 mg m⁻³ month⁻¹), 2008 (0.15 mg m⁻³ month⁻¹), and 2009 (0.25 mg m⁻³ month⁻¹). In 2002 and 2004, in contrast, [chla] time series remained flat at about 0.65 ± 0.26 and 0.48 ± 0.16 mg m⁻³, respectively. In addition to that, a sustained [chla] incremental increase during the summer-to-fall transition was observed in 2003, 2005 and 2007, just after the heat-flux inversion (positive to negative), i.e., when the ocean surface begin to cool (heat escaping the ocean).

The time evolution of b_bp was distinct from [chla], as b_bp presented a slight decrease during the spring-to-summer transition and remained relatively constant in summer. In fall, increases of b_bp had a small amplitude compared to [chla] (e.g., October 12 to 25, 2005 and 2007) or the b_bp time-series remained almost flat (e.g., 2002, 2003, 2006, 2008, and 2009).

Wind speed showed high variability, but its signature did not show a clear influence on [chla] (Figure 3). This mismatch between wind and [chla] were observed when wind components (zonal and meridional) were analyzed (Figure S2). The mismatch is particularly evident when high wind speed events (> 8 m s⁻¹) occurred and heat fluxes were positive. Such events were evident in late summer (August–September) and associated with a relatively low [chla] (e.g., 2002, 2003, 2004 and 2005) or followed by an ephemeral incremental increase in [chla]. Although strong wind events in summer, during positive heat flux, can result in an increase in [chla], as observed in 2003 (around August 13), they can also result in a decrease in [chla] (e.g., August 28 and September 27, 2003). For most years analyzed, peaks of wind speed during summer (between August 15 and September 15) were associated with [chla] lower than their respective phenological peaks. For example, when high wind speed (exceeding 6 m s⁻¹) was maintained between August 20 and September 12, 2007, [chla] remained low (< 0.8 mg m⁻³) relative to the spring and fall blooms (2.0 mg m⁻³).

In summary, wind speed versus [chla] links in late summer or fall were hard to detect because of the high wind-speed variability, i.e., ranging from 2 m s⁻¹ to above 10 m s⁻¹ during a short time period, and because a large increment in [chla] often happened in the presence of moderate wind. However, strong wind events had more impact on [chla] when the ocean was cooling (e.g., 2005, 2007). Other factors such as heat fluxes and freshwater inputs have to be considered.

We extracted correlation coefficients between [chla] and wind speed (WS) or heat flux (HF) using a moving average filter with a time window ranging from 1 to 30 days applied to daily time series between 1998 and 2008. For each window range, WS and HF generated a pair of coefficients of correlation with [chla]. We segmented time series between pre-fall and fall seasons. The first day that the atmosphere-cooling ocean began (i.e., the day of heat-flux-inversion) defined the fall seasons. As a result, Figure 4 shows a correlation analysis of [chla] with atmospheric forcing (WS and HF) on different time scales in southern Hudson Bay. In the spring-to-summer transition period, before the heat-flux inversion, wind speed and heat flux had a weak correlation outside the confidence interval of 95% (p < 0.05). After the heat-flux inversion, [chla] and wind speed reached a coefficient of correlation of 0.42 (p < 0.05), while [chla] and heat flux reached a high negative coefficient of correlation of -0.71 (p < 0.05) for the moving average filter of a 7-day window. The coefficient of correlation between [chla] and wind speed only had the confidence interval of 95% for a filter length of 7 days or longer, which indicated the low impact of high-frequency wind events on the phytoplankton dynamic.

Figure 5A shows well-marked spring and fall blooms reaching 2 mg m⁻³ separated by an oligotrophic period of [chla] below 0.5 mg m⁻³ in 2005. The annual lowest sea surface salinity (< 29.5) and highest stratification (N² ∼ 1.5 × 10³ s⁻¹) were reached on June 15 (Figure 5B), as a consequence of freshwater lenses produced by sea-ice melting and riverine input. After that date, there was a transition in the control of vertical stratification (i.e., β to α transition: Carmack and Wassmann, 2006; Carmack et al., 2006) triggered on June 23, associated with the freshwater lens dissipation, mixed layer expansion, and pycnocline ventilation. The water column stratification, which was controlled mainly by the salinity gradient becomes controlled by the temperature gradient between July and August (Figure 5C and D). Pycnocline ventilation was fed continually with freshwater lens dissipation (Figure 5B and C). On August 10, the mixed layer depth (dashed line) definitively exceeded the first optical depth (z_OC, grey line), which was stably positioned around 10 m. The beginning of mixed layer expansion corresponded to the oligotrophic period of 2005 with lowest [chla] observed (< 0.2 mg m⁻³). The [chla] decreased systematically following the ice-edge bloom dissipation from June 15 to August 18.

The phytoplankton fall bloom was observed as a systematic increase in [chla], coinciding with cooling throughout the water column and the heat-flux inversion observed on September 11 (Figure 5A). The vertical gradients of temperature and salinity remained weak in the upper ocean layer (top 30 m) when an atmosphere-cooling ocean began. The heat-flux inversion also was coincident with the intercept between pycnocline and z_eu (Figure 5A and B). The abrupt increase in [chla] coincide with a decreasing trend in wind speed starting on October 5 (Figure 5A). Although [chla] observed during the fall blooms was similar to the spring-to-summer season in ice-edge booms, the dispersion of daily compared to 7-day averages revealed that the variability of [chla] in fall was clearly higher than in the ice-edge blooms (Figure 5A). Finally, an intensification of fall bloom was marked by another abrupt increase in [chla] (> 2 mg m⁻³) after October 13 (Figure 5A), which was triggered by an intersection between the mixed layer and euphotic depth (Figure 5B–D), and ceased/collapsed by ocean freezing and sea-ice recovery.

To gain further insights about the phytoplankton growing in the late summer and fall seasons, we first report the relationship between particle backscattering and chlorophyll-a measured in the field in summer. Figure 6A shows the relationship between in situ $[chla]HPLC$ and b_bp(555) obtained from field measurements in Hudson Bay in 2018. Stations were grouped in distinct light environments based on $kd(λ)$ with high, moderate, and low light attenuation in marginal ice zones. Values for $[chla]HPLC$ ranged from 0.1 to 15 mg m⁻³, b_bp(555) from 0.0004 to 0.02 m⁻¹, and $bbp*(555)$ from 0.0007 to 0.0506 m² (mg chla)⁻¹. These observations were used to test different fitting algorithms to estimate $bbpk$⁠, before applying it to satellite ocean color data. As mentioned in Section 2.5, we compared the linear fit of Behrenfeld et al. (2005), the nonlinear least-squares algorithm of Bellacicco et al. (2019) and an exponential model (Equation 7). The later model performed better with r² of 0.760 (rms.e = 0.001, n. = 37), compared to r² of 0.707 and 0.30 for the Bellacicco et al. (2019) and Behrenfeld et al. (2005) models, respectively. The estimated $bbpk$ using the exponential model yield values of 0.0023 m⁻¹, which is slightly higher than the value obtained using Bellacicco et al. (2019), and similar to the values obtained with the linear fit. Similar results and conclusions were obtained in Baffin Bay (results not shown).

The carbon-to-chlorophyll ratio estimated after the application of Equations 8 and 7, C_phy:$[chla]HPLC$⁠, yielded a mean value of 17.31 and maximum value of 120. Figure 6B shows chlorophyll-specific backscattering at 555 nm (⁠$bbp*(555)$ in m² (mg chla)^-1) as a function of $[chla]HPLC$ in distinct light levels (colored symbols). As reported for high latitudes (Reynolds et al., 2001; Stramska et al., 2003; Wang et al., 2005; Zhuang et al., 2020), we observed an inverse relationship between $bbp*(555)$ and $[chla]HPLC$⁠.

Besides [chla] (Figure 2), ocean color satellite monthly climatologies revealed the seasonal and spatial variability of particle backscattering in the HBS, which is assumed to respond both phytoplankton and non-algal particle concentrations. The monthly climatology of b_bp(555) (Figure 7) shows the highest values (> 0.002 m⁻¹) near the coast, simultaneously with the seasonal peak of river runoff and ice melting between May and June. The b_bp was maximum in May (mean ± standard deviation: 0.0025 ± 0.0005 m⁻¹, n = 12312) and June (0.0022 ± 0.0006 m⁻¹, n = 60957), decreased in July (0.0019 ± 0.0006 m⁻¹, n = 69238), and remained constant in August (0.0016 ± 0.0005, n = 70813), September (0.0016 ± 0.0005 m⁻¹, n = 70841) and October (0.0016 ± 0.0005 m⁻¹, n = 62450).

Similar to Bellacicco et al. (2020), $bbpk(555)$ was estimated by applying Equation 7 to all satellite daily observations available for a given month between 1998 and 2018, and for a given 4.5-km resolution pixel. The spatially resolved $bbpk(555)$ monthly maps are shown in Figure 8 and statistics for the whole HB in boxplot in Figure S9. The seasonal minimum coincided with fall blooms in September (0.0004 ± 0.0004 m⁻¹, n = 70839) and October (0.0002 ± 0.0003 m⁻¹, n = 62442), and ice-edge blooms in May (0.0004 ± 0.0005 m⁻¹, n = 12304). Seasonal low contribution of NAP to total backscattering (⁠$bbpk$⁠:b_bp) also occurred in May (14.6%) and October (13.8%). The $bbpk(555)$ reached the seasonal peak after the ice-edge bloom season in July (0.0006 ± 0.0007 m⁻¹, n = 69238), simultaneous with its peak of NAP contribution to total backscattering (⁠$bbpk$⁠:b_bp of approximately 36.2%).

Monthly phytoplankton carbon, C_phy, was estimated by subtracting $bbpk$ (Figure 8) from b_bp (Figure 7) and multiplying by a scaling factor of 13,000 mg C m⁻² (Equation 8) (Behrenfeld et al., 2005; Bellacicco et al., 2020). The seasonal patterns in C_phy (Figure 9) slightly differ from [chla] (Figure 2) and b_bp (Figure 7). In general, C_phy was higher than 27.9 mg C m⁻³ in May and 21.1 mg C m⁻³ in June during the ice-edge bloom, then decreased to reach a minimum in August (< 12.4 mg C m⁻³). During the summer months, C_phy remained high in the southern HBS and in coastal waters of Hudson Strait. An overall increase in C_phy was observed in September (15.7 ± 5.9 mg C m⁻³, n = 70841, and during fall blooms in October (18.3 ± 5.9 mg C m⁻³, n = 62450; Figure S9).

The seasonal variability of [chla] and C_phy resulted in considerable seasonal variability in the $Cphy:[chla]$ ratio (Figure 10). Here, we present the seasonal evolution of this parameter, which is known to reflect phytoplankton photoacclimation, nutrient availability, and temperature (Geider, 1987; Geider et al., 1997; Sathyendranath et al., 2009; Jakobsen and Markager, 2016). In May, the phytoplankton found in the marginal ice zone presented high concentration of carbon (C_phy: 27.2 ± 7.9 mg C m⁻³, n = 12312) and chlorophyll-a ([chla]: 1.13 ± 0.44 mg m⁻³, n = 12312). The [chla] higher than 2 mg m⁻³ dominated the open water regions near the coasts, with maximum C_phy (> 20 mg C m⁻³) observed in polynyas, i.e., in the Northwest Hudson Bay and south of Belcher Islands. C_phy:[chla] ratios were particularly low along the ice edge in the northwest polynya. C_phy:[chla] increased systematically offshore from May to June. High C_phy:[chla] (> 40) coincided with maximum annual solar flux (summer solstice), and remained high in oligotrophic regions in the central HBS until August. Relatively low C_phy:[chla] values were found in eastern HB and in coastal waters of the Hudson Strait where high [chla] was also evident (Figure 2). From August to October, C_phy:[chla] tended to decrease everywhere in the HBS, except in some coastal zones (e.g., Ungava Bay). The seasonal minimum was reached in October (22.0 ± 7.2, n = 62450) with values lower than those observed during the spring bloom in May (26.5 ± 9.7, n = 12312). At this time of the year, the phytoplankton contribution to backscattering (b_phy:b_bp) was largest (approximately 86%).

How is the phytoplankton fall bloom triggered and maintained in Hudson Bay? The synergy of the ocean color satellite data, sea-ice dynamic model output, climatic reanalysis, and in situ bio-optic profiles underscores the influence of atmospheric forcing and water-column evolution on phytoplankton dynamics.

As observed in Figure 5, analogous to other years presented in Figures S3–S8, the seasonal evolution of the water-column structure follows these stages: 1) freshwater lens dissipation during summer; 2) vertical displacement between mixed layer and pycnocline; 3) mixed layer expansion and pycnocline ventilation; 4) a transition between a β (haline) to α (thermal) stratification; and 5) euphotic depth reaching the pycnocline and, subsequently, the depth of the mixed layer.

Climatic maps of [chla], C_phy, and C_phy:[chla] (Figures 2, 9 and 10), revealed that the ice-edge blooms in the spring-to-summer transitions and the rise in phytoplankton in the summer-to-fall transitions were productive features interrupted by a relatively severe oligotrophy. This oligotrophic period was defined by a high seasonal C_phy:[chla] ratio in regions of [chla] below 0.5 mg m⁻³, which occupied more than 50% of the bay. Ice-edge and fall blooms are very distinct from each other in terms of water column structure, stratification, and heat flux. In spring, the freshwater lens produced the highest annual stratification near the sea surface, while mixed layer expansion and pycnocline ventilation affected the whole euphotic layer in fall.

As presented in Figure 5, the annual maximum of stratification coincided with the seasonal peak of riverine input and ice melting in June and July (Prinsenberg, 1984). According to Déry et al. (2005), riverine inflow can exceed 4 km³ day^-1 in midsummer. In the upper layer, freshwater lenses (salinity < 30) produced strong stratification, which held both the mixed layer and pycnocline close to the surface layer (approximately 10 m; Stewart and Lockhart, 2005; Granskog et al., 2011). After the seasonal peaks of ice-melting and continental drainage (Prinsenberg, 1984; Déry et al., 2005), the dissipation of freshwater lenses triggered a continuous process of pycnocline ventilation and mixed-layer expansion until sea surface freezing at the end of autumn (Figure 5). Meanwhile, after the peak of heat flux (atmospheric warming ocean) and sea-ice retreat, the vertical stratification, previously controlled by salinity gradients, was then switched by temperature gradient, i.e., it switched from β to α stratification (Carmack and Wassmann, 2006; Carmack et al., 2006). As discussed by Carmack and Wassmann (2006) and Carmack et al. (2006), phytoplankton production and vertical distribution of chlorophyll in thermal-stratified (α stratification) polar seas is more sensitive to mixing and convection than salinity-stratified (β stratification) shelves dominated by freshwater.

The pre-fall-bloom stage occurred when the euphotic zone became shallower until it found the maximum stratification layer, which was getting deeper response to pycnocline ventilation, and when [chla] had risen from 0.4 to 0.9 mg m⁻³ between September 11 and 19, 2005, the beginning of the heat-flux inversion period. Then, a massive increase in [chla] was triggered when mixed layer expansion reached z_eu. When this stage was reached, the euphotic layer and whole phytoplankton productivity layer became confined to the mixed layer, characterized by weak stratification and homogeneous profiles of salinity and temperature. Similar to observations reported by Fujiwara et al. (2018) in the Chukchi shelf, wind-driven and convective turbulence intensified the erosion of water column structure in south Hudson Bay. Subsequently, the fall bloom development, defined by a significant increase in [chla], followed homogenization of the whole productive layer (Figures 5 and S3–S8).

In south Hudson Bay, NEMO modeling revealed high interannual variability of freshwater input, with years when surface salinity was low (2006, 2007, and 2008; Figures S6–S8) and years when it was high (2002 and 2004; Figures S3 and S5). As a result, the response of phytoplankton to atmospheric forcing will depend on freshwater inputs, which vary from year to year. Riverine input can increase the stocks of nutrients and remineralization in the whole bay. However, freshwater counteracts brine rejection blocking the upward flux of nutrients from the deep nutrient pool in winter (Eastwood et al., 2020). Phytoplankton can rapidly assimilate nutrients from rivers, which confines high production only to coastal areas and river plumes in summer (Kuzyk et al., 2010; Jacquemot et al., 2021). Indeed, the interannual variability of complex interplays of nutrient inventories, an atmospheric-heating ocean, and wind-driven turbulence throughout summer are pre-conditions for the phytoplankton dynamics in the fall season. For example, years of predominantly low upper-layer stratification (e.g., 2006) can allow nutrient replenishment in the euphotic layer, sustaining an increase of chlorophyll throughout the summer but weakening an eventual fall bloom triggered by convection due to prior depletion of the subsurface nutrient pool. On the other hand, strong stratification counteracts vertical transport of nutrients from below the euphotic zone and reinforces oligotrophy. The SCM responds directly to water column structure in Hudson Bay as discussed by Sibert et al. (2010; 2011), Ferland et al. (2011), Estrada et al. (2012), and Lapoussiere et al. (2013). In summer, water column structure can restrict phytoplankton production in the sub-surface chlorophyll-a maximum causing a nutricline deepening at levels below the range of convection-driven turbulence, which will dampen the fall bloom process (e.g., 2002 and 2004).

Turbulent mixing is a key to understanding the bloom onset. As reported by Falkowski and Oliver (2007), the climate-driven processes that influence turbulent mixing in the ocean seem to have strongly influenced the diversity and relative abundance of the major eukaryotic phytoplankton taxa in the geological (Phanaerozoic Eon) and contemporary ocean (Anthropocene). Sources of turbulence in the upper layer occur over a wide range of spatial and temporal scales, and result in non-homogeneous vertical distribution of turbulence in the mixed layer. As detailed by Franks (2015), turbulence in the ocean is a result of energy dissipation through the following processes: wind, waves, Langmuir circulation, and convectively driven turbulence. Both wind-driven and convection-driven (heat-fluxes) processes were examined for fall bloom onset in HB (Figures 4 and S2). The convection impacts uniformly and strongly the mixed layer defining the seasonal pycnocline. Otherwise, wind-driven turbulence is likely dominant in the upper layer, decays with depth, and dissipates in a short time scale (hours or less). Wind-driven turbulence can impact negatively phytoplankton production in the upper layer, especially in the first optical depth, (Behrenfeld, 2010; Taylor and Ferrari, 2011; Fischer et al., 2014; Franks, 2015). In polar seas, an increase in annual phytoplankton production is assumed, in part, to be caused by storm events that drive turbulence in the upper ocean layer (Ardyna et al., 2013; Fujiwara et al., 2018). However, a direct link between storms and fall blooms is not a general consensus, as pointed by Sigler et al. (2014). In Hudson Bay, sea-surface wind events only have a direct impact on [chla] after heat-flux inversion, i.e., when the atmosphere began to cool the upper-ocean layer (Figure 3). However, the increase in [chla] after the HF inversion is not systematic and probably depends on the nature of the stratification (α versus β stratification), i.e., the presence of freshwater in the upper ocean.

As discussed by Kim et al. (2007) for the Japan/East Sea and by Olita et al. (2014) in the Mediterranean Sea, the phytoplankton response takes some time to assimilate a new supply of nutrients introduced into the upper ocean layer. However, lagged correlations between wind speed and [chla] were not statistically significant in this study. Otherwise, convection-driven turbulence has a homogeneous effect on the whole mixed layer that determines its seasonal depth. As reported by Ferland et al. (2011) in summer (August and September), the depths of the nutricline and sub-surface chlorophyll-a maximum are very close to the euphotic zone depth in Hudson Bay. Therefore, convection-driven turbulence results in the mixed layer expansion and pycnocline ventilation, potentially releasing nutrients and phytoplankton communities from the SCM into the upper euphotic layer. The freshwater input and sea-ice melting impact the resilience of the water column structure against wind-driven turbulence because freshwater shields the upper-layer structure, and the absence of convection dampens the bloom development, as observed in situ during the summer-to-autumn transition on the northern Chukchi Sea (Nishino et al., 2015). Although nutrient resupply by Ekman pumping is an important factor in increasing phytoplankton production, excessive turbulence and mixing can have a negative effect. As demonstrated by Franks (2015), wind-driven turbulence strongly affects the upper layer but loses efficiency with depth. In addition, the restructuring of stratification in response to the wind-driven turbulence is a much faster process than the restructuring in response to convection (Franks, 2015).

The negligible impact on [chla] in response to persistent strong winds and high frequency wind events (Figure 3), coastal upwelling (Enriquez and Friehe, 1995) forced by winds parallel to coastal line orientation (Figure S2) and in the ice-edge (Dumont et al., 2010), or in the absence of convection prior to heat flux inversion (Figure 4), supports the idea that fall blooms are triggered when convective-driven turbulence begins to erode stratification. Specifically, wind-driven mixing alone can be more effective in limiting phytoplankton production in the first optical depth than in promoting it by nutrient transport via Ekman pumping. As presented in the section on the ice-ocean model and heat-flux calculations (Methods), the wind has a direct effect on heat flux. Although wind speed is key to understanding how atmospheric forcing impacts phytoplankton dynamics, our results highlight the importance of atmospheric cooling of the surface ocean on stratification. We found that the water column stratification resisted wind-driven turbulence events until the atmosphere effectively began to cool the ocean, producing a systematic increase in [chla].

As mentioned in Methods, phytoplankton physiological adjustments to environmental conditions result in changes in the ratio of phytoplankton-carbon biomass to [chla]. We estimated C_phy from satellite-derived particle backscattering after considering the non-algal backscattering background. We first discuss this backscattering background, followed by the use of C_phy:[chla] as a proxy to investigate the seasonal evolution of phytoplankton photoacclimation, and how this ratio informs the productivity of fall blooms.

All matter in the ocean scatters light, whether water molecules, sea-salt ions, or particles (Zhang et al., 2020). Natural aquatic assemblages have a variable pool of particle size, shape, and composition (organic/inorganic), which limits our ability to correctly interpret ocean color satellite observation of b_bp in terms of phytoplankton photo-physiological aspects (Organelli et al., 2018; Koestner et al., 2020). Uncertainties on the use of b_bp to infer phytoplankton community composition derive from the intricate effect of non-algal particles composed of heterotrophic organisms, organic detritus (e.g., fecal pellets, and cell debris), and inorganic sediments (e.g., clays, sand, calcite, and liths). For example, inorganic particles generally dominated the backscatter of assemblages because of the high refractive index of minerogenic particles, compared with organic-dominated assemblages (Effler et al., 2013; Reynolds et al., 2016; Koestner et al., 2020). Optically complex waters exhibit diverse particle types where phytoplankton and NAP represent varying proportions of the suspended particle pool, subject to potentially large and independent variations (Reynolds and Stramski, 2019). In order to gain more insights about phytoplankton-related backscattering, we assessed from satellite ocean color data the NAP contribution on a pixel basis using an approach similar to that proposed by Bellacicco et al. (2018; 2020). The satellite-derived $bbpk$ is the contribution of NAP to b_bp that does not covary with [chla]. Based on field observations, $bbpk$ was calculated using an exponential equation with offset (Figure 6), which performed better than the non-linear model proposed by Bellacicco et al. (2019) or the linear model of Behrenfeld et al. (2005). More importantly, we showed that $bbpk$ varies in space and time (Figures 8, S9, and S10) with the highest values obtained after the ice-edge bloom (in June or July, excluding the nearshore coastal waters) reaching a contribution to b_bp of about 35%. This result may reflect the production of non-algal particles such as micrograzers, bacteria, viruses, cell debris and fecal pellets during post-bloom conditions. This contribution of NAP, however, is relatively low compared to the global ocean, but similar to that reported in the productive waters of the North Atlantic (Bellacicco et al., 2018). At the ice edge in May, for example, we observed low $bbpk$⁠:b_bp (Figure S9) and high [chla], which could have involved two processes: i) fast-sinking aggregates of ice algae and sympagic communities that quickly remove NAP from the pelagic to the benthic realm (Lannuzel et al., 2020; Trudnowska et al., 2021); and ii) early ice retreat that results in a faint manifestation of under-ice and ice-edge blooms (Barbedo et al., 2020).

Provided that a satellite climatology successfully described the seasonality of $bbpk$ and b_bp in HBS, C_phy:[chla] appeared as a robust proxy to investigate the seasonal evolution of phytoplankton community size structure and photo-physiology in the summer-to-fall transition. However, we recognized that the conversion of b_bp to phytoplankton carbon may be controversial, as isolating the backscattering signal from phytoplankton cells alone (or phytoplankton carbon from POC) from bulk measurements in the natural environment is not possible (Sathyendranath et al., 2009). Therefore, C_phy includes the contribution of all particles that co-varied with [chla]. Absolute values of C_phy:[chla] reported here (Figure 10) in spring and fall seasons, or in coastal waters, are in the lower range compared to the literature in the open ocean (Sathyendranath et al., 2009; Bellacicco et al., 2016), but similar to the range reported for coastal waters of the North Baltic Sea (same latitude as the HB; ratios of 10–60; Jakobsen and Markager, 2016). A phytoplankton population composed of smaller cells backscatters more light than one composed of larger cells, for an equal amount of chlorophyll present (Vaillancourt et al., 2004; Reynolds et al., 2016; Koestner et al., 2020). Low values of C_phy:[chla] probably indicate the presence of larger cells with low backscattering efficiency and high chlorophyll-a content that are typically found in low light environments and/or nutrient-replenished (mesotrophic or eutrophic) environments (Devred et al., 2006; Barbieux et al., 2018).

As shown in Figure 10, the highest seasonal values of C_phy:[chla] were associated with oligotrophy in the summer. Therefore, satellite ocean color observations highlighted the photo-acclimation strategy of phytoplankton communities to deal with the high summer irradiance and nutrient depletion in the upper layer of the HBS. Satellite observations and vertical profiles acquired in the spring-to-summer transition help to improve our understanding about phytoplankton dynamics. In situ empirical relationships between b_bp(555), [chla], and PAR in distinct light-attenuation environments of the MIZ revealed that phytoplankton communities in the upper layer at high PAR levels had a relatively low intracellular pigment content and small size structure and, consequently, high $bbp*(555)$ and C_phy:[chla] (Figure 6). As reported by Westberry et al. (2008), for example, upper layer phytoplankton communities generally have low pigment content and small size structure cells to better deal with relatively high PAR, poor nutrient supply, and high wind-driven turbulence. On the other hand, the phytoplankton community around the SCM is photo-acclimatized to a balance between low PAR, relative high nutrients, high stratification and low turbulence, which results in relatively large cells and high intracellular chlorophyll-a content, i.e., relatively low $bbp*(555)$ and C_phy:[chla] (Brewin et al., 2012; Barbieux et al., 2018).

During the fall season, a decrease of C_phy:[chla] and negligible stratification in the productive layer may indicate an advection/seeding of SCM communities to the first optical depth. Otherwise, climatology-depicted photoacclimation as a seasonal process was synchronized with solar elevation and nutrients in Hudson Bay. For example, low C_phy:[chla] in an ice-edge bloom is likely due to high nutrient availability in a shallow mixed layer under the maximum seasonal light incoming in May and June. In summer, high C_phy:[chla] indicated the presence of small phytoplankton cells with relatively low intracellular chlorophyll pigments in oligotrophic domains under high PAR (Figures 2 and 9). Although very distinct water column structures are involved in ice-edge and fall blooms, in general similar and relatively low C_phy:[chla] characterized spring and fall, intercalated with a smooth transition through maximum C_phy:[chla] in domains of severe oligotrophy during summer. Phytoplankton photoacclimation, size structure, and community were similar in ice-edge and fall blooms, (Figures 2, 9, and 10). In addition, a seasonal increment of C_phy together with [chla] (Figures 2 and 9) and a decrease of $bbpk$ (Figures 8) and its relative contribution to b_bp (Figures S10 and S9) support the hypothesis of a productive fall bloom.

The present study achieved a synergy between hydrodynamic modeling, climate reanalysis, ocean-color imagery, and in situ observations to better understand fundamental aspects of phytoplankton dynamics, with special attention dedicated to the fall bloom phenomenon. Atmospheric cooling of the ocean and wind-driven turbulence can facilitate the transport of nutrients by mixed-layer expansion and pycnocline ventilation after the heat-flux inversion, when the atmospheric cooling of the ocean began. However, the interplay of phytoplankton and atmospheric forcing is not straightforward in a region where both haline and thermal processes affect water-column stratification.

Satellite-based optical proxies revealed that phytoplankton communities with a low carbon-to-chlorophyll ratio observed at the surface in the fall may be advected from the subsurface chlorophyll-a maximum to the depth visible by spaceborne sensors. However, the significant increase in phytoplankton-related carbon concentration (while non-algal backscattering remained constant) indicated that these fall blooms are potentially productive. To quantify the relative importance of fall blooms in Arctic and sub-Arctic seas, the carbon-to-chlorophyll ratio can improve the parameterization of phytoplankton primary production from satellite or coupled biogeochemical models (International Ocean Colour Coordinating Group, 2015; Sathyendranath et al., 2020).

The persistence of fall blooms can enhance the carbon pump and mitigate Arctic warming by removing greenhouse gases (carbon dioxide). However, the biogeochemical capability of the Arctic and sub-Arctic seas to assimilate both a high content of organic matter discharged by rivers and that produced in situ by phytoplankton blooms in a rapidly changing Arctic scenario remains an open topic.

The day of in situ backscattering coefficient and chlorophyll-a concentration in this research are made available through the BaySys data repository DOI: https://dev.uni-manitoba.links.com.au/data/project/baysys.

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

Text S1. Evaluation of chlorophyll-a concentration ocean color satellite algorithms.

Figure S1. In situ evaluation of satellite-derived chlorophyll-a concentration algorithms.

Text S2. Wind direction effect on chlorophyll-a concentration.

Text S3. Water column structure evolution in the South Hudson Bay.

Figure S2. South Hudson Bay: phytoplankton dynamic and wind components.

Figure S3. Atmosphere, ocean and phytoplankton, 2002.

Figure S4. Synergy of the sea-ice model, wind speed reanalysis, and ocean color satellites in south Hudson Bay, 2003.

Figure S5. Synergy of the sea-ice model, wind speed reanalysis, and ocean color satellites in South Hudson Bay, 2004.

Figure S6. Synergy of the sea-ice model, wind speed reanalysis, and ocean color satellites in South Hudson Bay, 2006.

Figure S7. Synergy of the sea-ice model, wind speed reanalysis, and ocean color satellites in South Hudson Bay, 2007.

Figure S8. Synergy of the sea-ice model, wind speed reanalysis, and ocean color satellites in South Hudson Bay, 2008.

Text S4. Seasonal phytoplankton photo-acclimation.

Figure S9. Seasonal influence of bio-optical properties in the Hudson Bay.

Figure S10. Seasonal influence of non-algal particles on light backscattering coefficient.

Supplementary references

This work is a contribution to the Natural Sciences and Engineering Council of Canada (NSERC) Collaborative Research and Development (CRD) project titled BaySys led by Dr David Barber. 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.

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

The authors have declared that no competing interests exist. Jean-Éric Tremblay is an associate editor at Elementa. He was not involved in the review process of this article.

Conception and design: LB, SB.

Analysed, interpreted data, manuscript preparation: LB, SB, JVL, PGM.

Acquisition and processing of marine bio-optical data: LB, SB.

Sea-ice dynamic modeling: PGM.

Critical review of the manuscript: LB, SB, JVL, PGM, JET.

Final approval of the versions to be submitted: LB, SB, JVL, PGM, JET.

How to cite this article: Barbedo, L, Bélanger, S, Lukovich, JV, Myers, PG, Tremblay, J-E. 2022. Atmospheric forcing and photo-acclimation of phytoplankton fall blooms in Hudson Bay. Elementa: Science of the Anthropocene 10(1). DOI: https://doi.org/10.1525/elementa.2021.00067

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

Associate Editor: Kevin R. Arrigo, Department of Earth System Science, Stanford University, Stanford, CA, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: The Hudson Bay System Study (BaySys)