Exploring controls on the timing of the phytoplankton bloom in western Baffin Bay, Canadian Arctic

Baffin Bay is an important gateway between the Arctic Ocean and the Northwestern Atlantic. Each year during the sea ice retreat, the phytoplankton spring bloom leads to biomass increasing by orders of magnitude, from its lowest values at the end of winter to a peak in summer (Perrette et al., 2011). Because phytoplankton forms the basis of the Arctic marine trophic network (Legendre and Rassoulzadegan, 1995), understanding the controls of the timing of the bloom is critical, especially given the shortness of the productive season. A mismatch between primary and secondary producers in the peak of activity and recruitment may decrease production at the higher trophic levels of the polar marine ecosystems (Søreide et al., 2010; Leu et al., 2011), and even impact the export of organic matter via the biological carbon pump (Henson et al., 2023). The biological carbon pump refers to mechanisms that export organic carbon from the ocean’s surface to its interior, where it may be sequestered (Boyd et al., 2019; Henson et al., 2023). Mismatch events have occurred in Arctic environments, but they could increase in frequency and intensity owing to climate change (Søreide et al., 2010). Some of the largest effects of climate change globally occur in the Arctic (Gutiérrez et al., 2021; Rantanen et al., 2022). Moreover, as people living along western Baffin Bay rely partly on subsistence harvest (Kenny and Chan, 2017), unexpected changes in biological production may negatively impact their access to local fishery and marine mammal stocks, and ultimately their food security.

Arctic phytoplankton bloom dynamics are split between three periods separated by two crucial days: the pre-bloom period, the day of the bloom initiation, the growth period, the day of the bloom peak and the post-bloom period (Sakshaug, 2004; Carmack and Wassmann, 2006; Wassmann and Reigstad, 2011). During the pre-bloom period, very low light prevents significant population growth even though nutrient concentrations are high. Snow accumulation allows a transmittance of less than 1% through snow-covered sea ice, but as soon as snow starts to melt transmittance increases to between 2% and 10%, thus allowing the start of the growth period (Ardyna et al., 2020a). For instance, a bloom initiation was detected by floats under 100% sea ice cover as early as February in Baffin Bay (Randelhoff et al., 2020). Data from biogeochemical-Argo floats in the Greenland Sea show that neglecting under-ice blooms would have resulted in the underestimation of the annual net community production of phytoplankton by 52% (Mayot et al., 2018). The controls of the bloom initiation appear to be linked to the optical characteristics of the ice and snow cover, and to the physiological response of phytoplankton to severe light limitation, as was also observed in the Southern Ocean (Hague and Vichi, 2021). The middle of the growth period is associated with melt ponds increasing the proportion of irradiance reaching the water column by up to 25%–31%. The growth period ends with the formation of the marginal ice zone when transmittance in the water column reaches 40% to 60%. At this point, the absence of limitation by either light or nutrients produce favourable conditions for exponential growth until the bloom peak is reached (Wassmann and Reigstad, 2011). The initial stock of phytoplankton biomass at the onset of the growth period plays a role in the timing of the bloom peak for the Southern Ocean (e.g., Sakshaug et al., 1991) and the Arctic Ocean (Christian et al., 2022). The post-bloom period follows in response to nutrient depletion and grazing pressure. This period is characterized with a negative or near equilibrium biomass accumulation rate. In the Arctic Ocean, nutrient depletion is thought to exert a stronger control than grazing (Sakshaug, 2004; Randelhoff et al., 2019).

Here we use a complex ecosystem model within an idealized one-dimensional (1-D) water column representative of the conditions in Baffin Bay to better understand the crucial processes controlling the timing of the spring bloom. In particular, we consider the role of light, nutrients and overwintering phytoplankton biomass in the phenology of the spring bloom. Light is thought to control bloom initiation (Ardyna et al., 2020b), while nitrate is the limiting nutrient in Baffin Bay where its eventual decrease is thought to cause the end of the bloom (Randelhoff et al., 2019). The phytoplankton standing stock at the end of the winter is also relevant (Sakshaug et al., 1991; Christian et al., 2022). There remain many questions on how phytoplankton survive over winter. Here we ask how important this survival is to the spring bloom dynamics relative to the other controls. To this end, we conducted sensitivity experiments to better understand and quantify the controls exerted by light, nutrient level and overwintering phytoplankton biomass at the end of winter.

The model simulations were configured to represent the location of the ice camp of the Green Edge mission in western Baffin Bay at (67.48°N, 63.79°W; Figure 1; see Oziel et al., 2019, for details) that provided rich datasets on oceanographic, biogeochemical and ecological properties of the site (Massicotte et al., 2019; Massicotte et al., 2020; see Data S3 in Benoît-Gagné et al., 2024). The camp was on seasonal landfast sea ice near Qikiqtarjuaq (Nunavut, Canada) with a water depth of 360 m (Massicotte et al., 2020). It will be referred to as the Qikiqtarjuaq ice camp from here onwards. Field observations were collected at the ice camp in 2016 from April 27 to July 22.

During the camp, 134 variables were measured including snow thickness, ice thickness, underwater photosynthetically active radiation (PAR, 400 to 700 nm), nitrate and silicic acid concentrations, chlorophyll a (Chl a) and depth of the mixing layer (in both Oziel et al., 2019, and Massicotte et al., 2020). Notations and units mentioned in the main text are described in Table S1. Additionally, there were estimates of the carbon biomass for several plankton taxonomic categories (in both Grondin, 2019, and Massicotte et al., 2020) from an Imaging FlowCytobot (IFCB; Olson and Sosik, 2007; Sosik and Olson, 2007; Moberg and Sosik, 2012; Laney and Sosik, 2014). The IFCB combines microscopy and flow cytometry to produce high-speed images of phytoplankton cells. These images can be used to identify species for cells larger than 10 µm and to identify broader taxonomic groups for cells between 3 µm and 10 μm. Further information is available in  Appendix A5.

The biogeochemical/ecosystem model follows from Dutkiewicz et al. (2021) but was modified here for the Arctic Ocean. We present a brief overview of the simulated functional groups and the size classes. Details can be found in Dutkiewicz et al. (2015; 2020) and  Appendix A1. Versions of this ecosystem model have been run and tested in 1-D configurations that were not specific for the Arctic Ocean (Hickman et al., 2010; Wu et al., 2021).

In this study, the numerical planktonic ecosystem included 26 phytoplankton types (Figure 2a) divided into four biogeochemical functional groups, each defined by a few physiological features. Picophytoplankton were the smallest (1.4 µm and 2.0 µm equivalent spherical diameter, ESD), diatoms (6.6–154 μm) required silicic acid, and mixotrophic dinoflagellates (6.6–228 μm) were capable of photosynthesis as well as grazing on other plankton. An additional group of ‘other nanophytoplankton’ were analogs of all other phytoplankton types of similar size that were neither diatoms nor mixotrophic dinoflagellates (e.g., haptophytes such as coccolithophores or Phaeocystis). The phytoplankton types differed from one another by their maximum growth rate (P^C_max; Figure 2b) and their half saturation for growth on nitrate (⁠$kNO3$⁠; Figure 2c) and were allometrically (i.e., using a relationship between cell size and the parameter) assigned within functional groups (Dutkiewicz et al., 2020; 2021). The half saturation for growth on silicic acid (⁠$kSi(OH)4$⁠), on phosphate (⁠$kPO4$⁠), on iron (⁠$kFe$⁠) and on ammonium (⁠$kNH4$⁠) for each type was calculated from the corresponding $kNO3$as described in  Appendix A1.2. The photosynthetic parameters of the model diatom analogs corresponded well to laboratory observations from Arctic diatoms (Figure S1).

Subscript j refers to a specific phytoplankton type and subscript k refers to a zooplankton type. C_j is the biomass of phytoplankton j and Z_k is the biomass of zooplankton k. The source and sinks of the biomass (in carbon) of each phytoplankton type j (S_C,j) are calculated such that

The first right-hand term represents biosynthesis. The limitation for growth by nutrients, temperature and light (γ^R, γ^T, γ^I, respectively) is between 0 (e.g., a lack of light) and 1 (e.g., unlimiting light). The calculations of γ^R, γ^T and γ^I are described in  Appendices A1.2,  A1.3 and  A1.4, respectively. The second term represents the compounded losses of biomass due to respiration, senescence, viral lysis and excretion. The relative impacts of these processes were not resolved individually. Instead, a bulk estimation of these loss processes is calculated from a constant ‘mortality’ rate at 30°C (m_p,j) of 0.1 d⁻¹. The third term represents grazing (see  Appendix A1.5 for the calculation of the grazing rate of j by k, g_jk). The fourth term represents the sinking with a constant sinking rate (⁠$ωp,j$⁠) of 0.07 m d⁻¹ for picophytoplankton, 0.36 m d⁻¹ for diatoms and 0.23 m d⁻¹ for dinoflagellates and other nanophytoplankton. The second and third terms were set to 0 when C_j dropped below a threshold of minimum biomass for phytoplankton j (C_min,j) of 10⁻² mmol C m⁻³ to simulate survival in winter. However, sinking and mixing could still dilute the biomass of phytoplankton j below this threshold. The model includes explicit parameterization of Chl a, such that each phytoplankton has a dynamic Chl:C ratio that alters due to acclimation following Geider et al. (1997; 1998). For the equations and more details, the reader is referred to Dutkiewicz et al. (2015).

We also consider 16 zooplankton numerical types differing in size, resulting in a marine ecosystem containing a complex planktonic community. The maximum grazing rate of the grazers (g_max) depends on their biogeochemical functional group (mixotrophic dinoflagellates or microzooplankton) and their size following an allometric relationship. These g_max values are parameterized from the observations by Taniguchi et al. (2014) and Jeong et al. (2010; Figure 2d). The four smallest microzooplankton are an exception, with a g_max independent of their size (following a lack of observed allometric relationship between these smallest types, as in Taniguchi et al., 2014).

We used a 1-D configuration representing the specific location of the ice camp in Baffin Bay. The configuration had 75 levels ranging in thickness from 1 m near the surface to 6 m near the bottom (which was 360 m). Temperature, salinity and vertical turbulent diffusivity (K_z) were provided as offline forcing fields. They were generated with a 1-D simulation of the LIM 3.6 sea ice model (Rousset et al., 2015) coupled to the ocean component of the general circulation model NEMO 3.6 (Madec et al., 2017). Hereafter these offline forcing fields will be referred to as NEMO-LIM3 (Data S4 in Benoît-Gagné et al., 2024). The modelled vertical diffusion, K_z was evaluated by comparing with two different metrics of the vertical mixing: the depth of the mixing layer and the depth of the equivalent mixed layer (h_BD). The term h_BD is the depth of the ‘buoyancy deficit’ as in Randelhoff et al. (2017). The depth of the mixing layer was measured at the ice camp only on June 23, 2016, and corresponded to a K_z around 10⁻⁴ m² s⁻¹. The depth at which K_z = 10⁻⁴ m² s⁻¹ was considered as the simulated mixed layer depth herein.

Preliminary experiments revealed that two major adjustments to Dutkiewicz et al. (2021) were required for the Arctic setup: setting a minimum threshold below which phytoplankton experienced no losses, especially during the harsh Arctic winter, and including light under sea ice. The sensitivity experiments that allowed us to choose the best parameters are presented in Sections 3.3 and 3.4, respectively.

The initial conditions of nitrate, silicic acid and phosphate were prescribed from in situ observations in 2015 and 2016 (Massicotte et al., 2020). Data could be averaged between mid-April and end of May in the two years (Figure 3a and Data S2 in Benoît-Gagné et al., 2024) because the observed nutrients were constant during this time (Figures 4c and S2a, b and c for 2016; not shown for 2015) and the stable sea ice conditions were likely responsible for this stability in the nutrient profiles during winter. There is no advection in a 1-D model, which prevents the supply of nutrients to the surface layer of the water column and results in their depletion. Compensation for the non-existing lateral advection and the absence of nutrient inputs by river runoff in the model was then necessary. The solution selected was the relaxation (reinitialization) of simulated nitrate, silicic acid and phosphate concentrations from January 1 to May 15 to in situ observations from mid-April to the end of May, during the ice-covered period at the ice camp. A relaxation coefficient of 1/30 d⁻¹ was used. No relaxation occurs after May 15, 20 days before the start of snow melt and biological activity in the water column, according to the in situ ice camp data (Oziel et al., 2019, their figure 10). A sensitivity experiment on the level of winter nitrate is presented in Section 3.5.

Model output from a 10-year spin-up period was used to provide initial conditions of ammonia, nitrite, total iron and dissolved organic and particulate organic matter. The same initial spin-up model results provided the total phytoplankton biomass (mmol C m⁻³). This initial spin-up biomass was then divided equally between the 26 numerical phytoplankton types to be used as initial conditions for the reference simulation. The spin-up model output also provided the total zooplankton biomass which, again, was divided equally between the 16 numerical zooplankton types for initial conditions. At initialization, Chl a is acclimated to light, temperature and nutrients following Geider et al. (1998). But Chl a and the Chl:C ratio are calculated dynamically at each time step of the model.

The model time step was 1 h. The ‘reference’ simulation was run forward for another 10 years with repeating forcing fields. Model results shown are from the last year of simulation. There was no drift in nutrient concentrations after 3 years (e.g., nitrate in Figure S3). The phytoplankton established a regular pattern after 2 years, such that we can assume a ‘quasi-steady state’ by year 10, at which time the initial conditions were no longer influencing the simulation results.

Before the sea ice break-up on July 18, the observed downwelling plane PAR just below sea ice in photon density flux, E_d,i(z = 0⁻, PAR[Q]), was converted to the scalar PAR just below sea ice in photon density flux, E_0,i(z = 0⁻, PAR[Q]), as described in  Appendix A2.1. Observations from the Qikiqtarjuaq ice camp (Matthes et al., 2019) were used to estimate the conversion factors. After the sea ice break-up, the downwelling shortwave radiation just above surface in energy units, E_s(z = 0⁺, SW), from Smith et al. (2014) was transformed into the scalar PAR just below open water in photon density flux, E_0,w(z = 0⁻, PAR[Q]), as described in  Appendix A2.2. E_0,i(z = 0⁻, PAR[Q]) and E_0,w(z = 0⁻, PAR[Q]) were used as forcing fields (Data S5 in Benoît-Gagné et al., 2024). The scalar PAR at each depth (E₀) was calculated from E_0,i(z = 0⁻, PAR[Q]) before the sea ice break-up and from E_0,w(z = 0⁻, PAR[Q]) after the sea ice break-up as described in  Appendix A2.3.

Model evaluation was performed by comparing the results of a reference simulation (EXP-0, Table 1; Data S6 in Benoît-Gagné et al., 2024) with in situ observations from the ice camp. We explored the tenth year of the reference simulation by segmenting it into: pre-bloom, bloom initiation and growth phase, and bloom peak. The day of the bloom initiation is defined as the last day of a 7-day positive accumulation period (following Boss and Behrenfeld, 2010). The bloom peak is defined as the day of maximum Chl a. We conducted a series of sensitivity experiments (Table 1) to explore the controls of the bloom timing: the magnitude of the biomass before the bloom initiation (EXP-1), treatment of light under sea ice (EXP-2) and nitrate concentration before the bloom initiation (EXP-3).

This section describes the observations measured at the ice camp and presented with dots in Figures 4 and 5. At the sea ice camp, snow melt occurred during the first half of June (Figure 4a and Oziel et al., 2019). Sea ice became thinner and melt ponds were created from the middle of June until the full sea ice break-up on July 18. The sea ice break-up necessarily led to the end of the ice camp campaign. The increase in the underwater light field in June and July (Figure 4b) corresponded to a decrease in the observed nutrient concentrations (Figure 4c) and an increase in vertically integrated phytoplankton biomass (∑C_phyto, 0–100 m; Figure 4d) and chlorophyll a (∑Chl a, 0–100 m; Figure 4e).

The early peak of accumulation rate of ∑Chl a on May 9 was likely due to the flushing of sea ice algae from the melting ice. Hence, we considered the date of bloom initiation as May 27 (Figure S4 and Table S2), a date similar to the date calculated in Oziel et al. (2019). After the snow had melted during the first half of June, ∑C_phyto and ∑Chl a increased significantly in mid-June causing a drawdown in nutrients. The highest ∑Chl a observed was on July 15, right before the sea ice cover disappeared at the ice camp. Although ∑Chl a may have reached a value higher than that of July 15 after the end of the ice camp campaign, for the purposes of this study we assume July 15 as the ‘peak’ of the bloom. From the underwater light field, the depth of a ‘reference isolume’ at 0.415 mol photons m⁻² d⁻¹ (z_0.415) was calculated (Letelier et al., 2004; Boss and Behrenfeld, 2010; Oziel et al., 2019). The significant increase in ∑Chl a in mid-June was correlated to the shoaling of this reference isolume (Figure 4e and f). This isolume followed the observed equivalent mixed layer depth (h_BD, as in Randelhoff et al., 2017; Figure 4f).

The model was evaluated by comparing the simulated variables with observations. A Taylor diagram for the observational time series summarizes the resulting statistics (Figure 6). The correlation coefficient (angular position) and normalized standard deviation (radial position) are performed on log-normalized fields. The correlation coefficient of the nutrient concentrations was greater than 0.7 for each nutrient. Observed ∑Chl a was particularly well captured by the model with a correlation coefficient of 0.89. ∑Chl a variability was also well captured with a normalized standard deviation of 1.05. The model captures the diatoms and other nanophytoplankton biomass well (correlation coefficients of 0.82 and 0.75, respectively), and slightly worse for the dinoflagellates (correlation coefficient of 0.63). We discuss the time series of the reference simulation further in terms of three phases: pre-bloom (before May 23), bloom initiation and growth, and bloom peak on July 10.

This section describes the observed and simulated quantities before the simulated bloom initiation represented by the vertical green dotted line on May 23 in Figures 4 and 5. In both model and observations, snow melt had not yet started (Figure 4a) and simulated underwater PAR was low (Figure 4b). Simulated surface nutrient concentrations were high (Figure 4c), matching observations. Simulated integrated biomass (0–100 m, ∑C_phyto) and simulated integrated Chl a (0–100 m, ∑Chl a) were low (Figure 4d and e, respectively). The simulated biomass was higher than observed, though the simulated Chl a (2 mg Chl m⁻²) was on the lower bound of the observed values of Chl a (2–23 mg Chl m⁻²).

The bloom initiation occurred on May 23 in the reference simulation, 4 days before the observations (Figure 4e). A steep increase in the underwater PAR between bloom initiation and the complete melt of the snow cover on June 15 (Figure 4b) caused a slow growth period of both observed and simulated ∑C_phyto and ∑Chl a (Figure 4d and e, respectively). This increase in biomass was not enough to decrease the nutrient concentrations significantly. In both model and observations, the depth of the reference isolume increased (Figure 4f).

The disappearance of the snow cover on June 15 (black dotted vertical line on panels a and b of Figure 4) allowed the underwater PAR to increase enough to cause a fast growth period (Figure 4a and b). This fast growth period is visible in both observed and simulated ∑C_phyto and ∑Chl a between the complete melt of the snow cover on June 15 and the simulated bloom peak on July 10 (blue dotted vertical line on panels d and e of Figure 4). Simulated and observed nitrate was depleted at the end of this fast growth period, but not silicic acid (Figure 4c). During the fast growth period, the reference isolume was deeper than observed (Figure 4f). We discuss this discrepancy in Section 3.2.3.

The simulated biomass of the diatoms increased at the same time as the observations, though the model overestimated the biomass at the start and end of the growth period (Figure 5a). Similarly, the simulated biomass of the other nanophytoplankton increased at the same time as the observations, though the biomass was too high over the full season (Figure 5e). The model captured the relative contributions (%) of diatoms and other nanophytoplankton (Figure 5b and f, respectively) during this period. In particular, the fast growth period between the disappearance of the snow cover and the bloom peak was associated with a change in the phytoplankton assemblage from other nanophytoplankton dominance to diatoms.

The depletion of simulated surface nitrate inhibited the growth of phytoplankton after July 10. The steep increase in surface PAR associated with the sea ice break-up of July 18 led to photoacclimation of the numerical phytoplankton. This decrease in the Chl:C ratio after the sea ice break-up caused a decoupling between ∑C_phyto and ∑Chl a. The ∑Chl a peak occurred on July 10, 5 days before the observations. ∑C_phyto continued to increase reaching a maximum on July 27 due to the decoupling. Following previous studies of the Qikiqtarjuaq sea ice camp (Oziel et al., 2019; Massicotte et al., 2020), maximum ∑Chl a rather than maximum ∑C_phyto was defined as the bloom peak.

The simulated maximum magnitude of ∑Chl a of 55 mg Chl m⁻² was underestimated relative to the observations of 204 mg Chl m⁻². The fact that the simulated ∑C_phyto was equal to or above the observations (Figure 4d) but that the simulated ∑Chl a was below the observations (Figure 4e) was indicative of an underestimation of the Chl:C ratio. This underestimated Chl a also led to light penetrating too deep in the simulation (Figure 4f). The trend of the simulated mixed layer depth (K_z = 10⁻⁴ m² s⁻¹) differed from the trend of the observed equivalent mixed layer depth (h_BD, as in Randelhoff et al., 2017) that was shoaling from 40 m to 10 m during this period of time. Despite underestimating the Chl a concentration, the temporal trends correlated well with the observations (Figures 4 and S2) and the subsurface chlorophyll maximum was captured (Figure S2). The biomass of the dinoflagellates was overestimated particularly after the bloom peak (Figure 5c).

We conducted a series of sensitivity experiments (Table 1) to explore controls of the bloom timing. In the model a minimum biomass threshold (C_min,j) was set below which a phytoplankton type experiences no losses such as grazing, maintenance or cell death to account for overwinter survival. In the reference simulation (EXP-0), C_min,j was set to 10⁻² mmol C m⁻³ for each phytoplankton type. This parameterization allowed the model to maintain winter Chl a like that observed at the ice camp (Figure 4e) and in other regions of the Arctic Ocean (see Section 4 for more discussion). In the first set of sensitivity experiments (EXP-1; Figures 7 and S5), this threshold was either increased (C_min,j = 10⁰ and 10⁻¹ mmol C m⁻³) or decreased (C_min,j = 10⁻³ mmol C m⁻³ and 0). A higher biomass before bloom initiation (C_min,j = 10⁰ and 10⁻¹ mmol C m⁻³) caused an earlier bloom peak on July 2 (−8 days) and July 5 (−5 days), respectively. A lower biomass before the bloom initiation (C_min,j = 10⁻³ mmol C m⁻³ and 0) delayed the bloom peak more substantially to July 16 (+6 days) and July 30 (+20 days), respectively.

In the second set of sensitivity experiments (EXP-2), we explored the importance of light in controlling bloom timing. In the reference simulation (EXP-0), the observed light just below sea ice, both for snow-covered sea ice and bare sea ice, was used as input for the model (Figures 8 and S6). In EXP-2.1, snow was considered opaque so that light just below snow-covered sea ice was set to 0 (Figure 8b). In EXP-2.2, both snow-covered sea ice and bare sea ice were considered opaque such that in this experiment there was no light below sea ice, both snow-covered sea ice and bare sea ice.

Removing light under snow-covered sea ice (EXP-2.1) delayed the bloom initiation to June 22 (+30 days) and the bloom peak to July 16 (+6 days; Figure 8c and d). Removing light under all sea ice (EXP-2.2) caused greater delays: in this case, the bloom initiation was delayed until July 25 (+63 days) and the bloom peak occurred only on August 7 (+28 days).

To explore the importance of the winter pool of nitrate on the timing of the bloom peak, a third set of sensitivity simulations considered different winter nitrate concentrations (EXP-3; Table 1). In the reference simulation EXP-0, the winter nutrients were relaxed to in situ observations at the Qikiqtarjuaq sea ice camps averaged between mid-April and end of May in 2015 and 2016 (Figure 3a). In the sensitivity simulations EXP-3 (Figures 9 and S7), the winter nitrate concentrations were instead relaxed to 1/4, 1/2, 2 and 4 times that in the reference simulation (Figure 3b).

Decreasing nitrate by a factor of 2 or 4 decreased the magnitude of the bloom peak from 55 mg Chl m⁻² to 40 mg Chl m⁻² and 27 mg Chl m⁻², respectively (Figure 9c). This decrease, in turn, caused a bloom peak slightly earlier on July 8 (−2 days) and July 6 (−4 days), respectively. The change in timing was much less than one would expect from the change in nutrients. Conversely, increasing nitrate before the bloom initiation by a factor of 2 or 4 increased marginally the magnitude of the bloom peak from 55 mg Chl m⁻² to 59 mg Chl m⁻² in both cases. Consequently, the bloom peak was barely delayed to July 11 (only +1 day).

Our study has shown that the overwintering biomass was an important control on the timing of the bloom peak (EXP-1; Figure 7). When the biomass before bloom initiation was lower than observed there was a delay in the peak Chl a. The simulated Chl a before the bloom initiation achieved with the C_min,j of the reference simulation (EXP-0) was between 10⁻² and 10⁻¹ mg Chl m⁻³, like the observed Chl a in May at the Qikiqtarjuaq ice camp (Figure S8). Randelhoff et al. (2020) measured Chl a with gliders in offshore Baffin Bay during the winters 2017–2018 and 2018–2019. Their values were also between 10⁻² and 10⁻¹ mg Chl m⁻³ (their figure 2e), supporting that a numerical winter Chl a above 10⁻² mg Chl m⁻³ was realistic. The impact of pre-bloom biomass on the timing of the phytoplankton bloom peak is not surprising and has also been highlighted for ice algal growth and bloom timing (Mortenson et al., 2017; Haddon et al., 2024).

In our model, a value range of 10⁻² to 10⁻¹ mg Chl m⁻³ for Chl a before bloom initiation was achieved by halting any further losses from metabolism, grazing or other mortality (terms 2 and 3 in Equation 1) when the biomass of one type of phytoplankton dropped below C_min,j. This parameterization was necessary to maintain an overwintering biomass, but is obviously an extreme oversimplification of complex processes controlling the survival of phytoplankton communities during winter. The numerical C_min,j, however, could compensate for excessive phytoplankton sinking and dilution in the model and could also represent the missing three-dimensional features of the sea ice such as production in leads. This simple implementation of a threshold biomass is the minimum parameterization needed in a biogeochemical model to mimic the particular conditions of the Arctic. It might also apply to the Austral Ocean given polar night there as well. A similar parameterization has been implemented by several model studies (Mortenson et al., 2018; Christian et al., 2022; Bertin et al., 2023; Haddon et al., 2024), but none of the studies provided detailed documentation of the impacts. Chl a in the oligotrophic subtropical gyres falls between 10⁻² and 10⁻¹ mg Chl m⁻³ (Kuhn et al., 2023), similar to winter observations in Baffin Bay. These Chl a values suggest that carbon biomass is low and comparable in both systems. However, these two systems differ in the seasonality of the zooplankton predation, as the pressure of zooplankton grazing is continuous in oligotrophic tropical gyres. Therefore, using a parameterization that stops grazing in oligotrophic tropical gyres would not be appropriate, and C_min,j should not be used globally. A simple workaround could be an empirical relationship between the latitude and the C_min,j threshold, while a more comprehensive solution could be the decomposition of the linear mortality into various loss terms, each one mechanistically formulated (e.g., senescence, viral lysis of phytoplanktonic cells, etc.).

However, a clearer understanding of what maintains the overwintering biomass is needed for newer and better parameterizations. By implementing a realistic overwintering biomass and demonstrating its importance, this study is a first step in its mechanistic implementation. This implementation would benefit from the results of laboratory experiments at extremely low light. For example, the limitation for growth by light (γ^I) could be modified in winter to allow more phototrophy at low levels (Hancke et al., 2018; Kvernvik et al., 2018). The mortality term (m_p,j) could be reduced in winter to represent reduced metabolism (Lacour et al., 2019; Kennedy et al., 2020; Joli et al., 2024) and sparing consumption of storage lipids (Morin et al., 2020). In winter, mixotrophy (Błachowiak-Samołyk et al., 2015; Vader et al., 2015; Marquardt et al., 2016; Kvernvik et al., 2018; Stoecker and Lavrentyev, 2018; Johnsen et al., 2020) could be enhanced or osmotrophy (Wen et al., 2002; Lavoie et al., 2018; Johnsen et al., 2020) could be activated. Winter grazing by zooplankton could be reduced by modifying the limitation of grazing by temperature (γ^T) in winter to represent the negligible losses to herbivorous grazing (Frost, 1993; Rose and Caron, 2007) or by implementing diapause on mesozooplankton (Baumgartner and Tarrant, 2017). Another strategy of survival during the polar night is resting stage formation (Johnsen et al., 2020).

Several studies have reported that Arctic diatoms can survive in the dark for months and start photosynthesis rapidly (within a few hours) when light returns (e.g., Lacour et al., 2019; Kennedy et al., 2020; Morin et al., 2020; Handy et al., 2024): a finding that has been confirmed in situ in Svalbard (Kvernvik et al., 2018). Both field experiments and laboratory experiments suggest that photosynthesis can occur at very low light from 0.17 to 1 μmol photons m⁻² s⁻¹ (Geider et al., 1986; Hancke et al., 2018; Kvernvik et al., 2018). Accumulation of biomass generally does not occur until higher light levels between 1.2 and 2.3 μmol photons m⁻² s⁻¹ are available (Geider et al., 1986; Ardyna et al., 2020b). However, biomass accumulation has been observed to occur at levels lower than previously thought, for example, in February in offshore Baffin Bay shortly after the winter solstice (Randelhoff et al., 2020). Recent field observations (Ardyna et al., 2020b; Randelhoff et al., 2020) show that even in Arctic marine ecosystems, blooms can start before stratification occurs, which is coherent with the dilution recoupling hypothesis (Behrenfeld, 2010; Behrenfeld and Boss, 2018). Importantly, the minimum light thresholds for photosynthesis and a positive accumulation rate remain uncertain. Further laboratory experiments with extremely low irradiations are needed.

The importance of under-ice PAR in triggering the initiation of the under-ice bloom is better understood (Ardyna et al., 2020a; 2020b). With the absence of light under snow-covered sea ice (EXP-2.1) and the absence of light under any sea ice (EXP-2.2), the bloom peak was delayed significantly (+6 and +28 days, respectively).

Finally, the magnitude of the bloom peak was only slightly affected by the winter pool of nitrate (EXP-3, Figure 9). Nitrate concentration before the bloom initiation had no effect on the initiation of the bloom and was only a modest control on the timing of the bloom peak. Because the accumulation of phytoplankton was exponential, a linear variation in winter nitrate and thus the magnitude of the bloom at its peak did not change the date of the bloom peak by much.

The model used in this study did not include a sympagic component. Previous modelling studies have shown some influence of ice algal production on the magnitude of the phytoplankton bloom in Arctic marine ecosystems (Hayashida et al., 2017; Mortenson et al., 2017). These studies suggested that seeding of phytoplankton production by ice algae sloughing from the sea ice in spring influenced the timing of both the bloom initiation and bloom peak by a few days. However, a review of in situ studies showed that only when extreme meteorological events trigger mass release of ice algae into the water column is seeding from sea ice an important control of the timing of under-ice blooms (Ardyna et al., 2020a). Another observation that may be interpreted as against seeding by ice algal sloughing is that the genera of polar diatoms observed in the winter water column (Vader et al., 2015; Kvernvik et al., 2018) are also those found in the spring phytoplankton bloom (Hoppe, 2022). Another impact of sea ice algae on phytoplankton production involves a shading effect. In simulations, ice algal shading can have a strong impact on the timing of the phytoplankton bloom under the ice when ice algal biomass is high (Castellani et al., 2017; Hayashida et al., 2019). In this study, observed PAR just below the (real) sea ice, and thus including the shading impact of sea ice algae, was used as a forcing field. Further research on the impact of sea ice algae is warranted but is beyond the scope of this study.

The physical variables at the location of the ice camp were close to the conditions of offshore western Baffin Bay ( Appendix A4). The conditions at the ice camp followed the dynamics described for offshore areas with phytoplankton bloom initiation occurring while there was still ice (Randelhoff et al., 2019; 2020). The processes at the ice camp and in western Baffin Bay (of which the ice camp is representative) are typical of an outflow shelf environment. Seasonal sea ice cover fosters sizeable pelagic blooms as soon as sea ice becomes translucent due to the deepening of the euphotic zone and the shoaling of the mixed layer depth (Randelhoff et al., 2019) in a manner similar to other outflow shelf environments of the Arctic Ocean, such as the Canadian Arctic Archipelago and the East Greenland shelf (Ardyna et al., 2020a). These outflow shelves are expected to react similarly to climate change (Ardyna and Arrigo, 2020). Thus, we believe that our study has broader implications than the single location studied here, although our study is specific to the polar regions.

Our study has shown that phytoplankton biomass at the end of winter is a key parameter for accurately modelling the spring phytoplankton bloom in a seasonal sea ice zone. Though a minimum threshold biomass parameterization has been used in previous studies, to our knowledge this study provides the first published sensitivity analysis. Our study also agrees with earlier results about the necessity of a reasonable representation of light transmittance through sea ice, especially for bare sea ice compared to snow-covered sea ice. Here we have shown that winter biomass and light under sea ice are comparably important controls for the timing of the bloom peak (+20 days when no winter biomass and +28 days when no light under sea ice). Research campaigns so far have generally concentrated efforts on measuring light in the water and have not focused on measuring winter biomass. Our results have shown that both observations will be important for further understanding of the spring phytoplankton bloom. This study highlights the need for field and laboratory experiments to gain a more precise understanding on the acclimation and adaptations by phytoplankton to maintain a balance between biomass growth and losses during the harsh Arctic winter. A better characterization of the underwater light field during the polar night will also be worthwhile, as darkness is never complete because of moonlight among other reasons (Cohen et al., 2020). Only with a better understanding of the mechanisms of survival of phytoplankton during winter we will be able to parameterize this aspect more realistically in models. Better models of Arctic ecosystems are urgently needed as this region has warmed four times faster than the global average since 1979, a phenomenon known as Arctic amplification (Rantanen et al., 2022). The crude parameterization of a minimum phytoplankton biomass threshold is, however, an initial step towards more accurately representing the timing of the bloom peak.

Only an overview of the model is presented here. More detail and equations can be found in Dutkiewicz et al. (2015; 2020) and at https://darwin3.readthedocs.io/en/latest/phys_pkgs/darwin.html (accessed May 30, 2024). See Tables S3, S4 and S5 for the units of the variables, the values of the constants and the coefficients for allometric scaling, respectively.

The carbon-specific photosynthesis rate of phytoplankton (P^C) is dependent on maximum growth rate (P^C_max; Figure 2b) and limitation for growth by nutrients, temperature and light (γ^R, γ^T, γ^I, respectively, between 0 and 1) such that

A limitation for growth of 0 means no growth because there are no nutrients, for example, or no light. A limitation for growth of 1 means no limitation for growth. The limitation for growth by nutrients, temperature and light is described in  Appendices A1.2,  A1.3 and  A1.4, respectively.

The half saturation for growth on nitrate (⁠$kNO3$⁠; Figure 2c) is modelled as in Dutkiewicz et al. (2020) using the model Monod formulation of growth rate (Follows et al., 2018):

The cell nutrient uptake half saturation constant on nitrate (⁠$KNO3$⁠; Figure S9a) is dependent on size but not on group as

where $aKNO3$ and $bKNO3$ are coefficients for allometric scaling and V is the cell volume. The maximum growth rate (P^C_max; Figure 2b) is dependent both on size V and on group g as

where $aPmaxC,g$ and $bPmaxC,g$ are the coefficients for allometric scaling dependent on group g. The cell minimum stoichiometric quota of nitrogen relative to carbon (⁠$QNmin$⁠; Figure S9b) is dependent on size but not on group as

where $aQNmin$ and $bQNmin$ are coefficients for allometric scaling. The cell nutrient uptake rate of nitrate relative to carbon (⁠$VNO3max$⁠; Figure S9c) is dependent on size but not on group as

where $aVNO3max$ and $bVNO3max$ are coefficients for allometric scaling. The cell elemental C:N:Si:P:Fe stoichiometry is 120:16:16:1:0.001. (Only diatoms have silicon.) The half saturation for growth on silicic acid (⁠$kSi(OH)4$⁠) is computed using this stoichiometry as

The half saturation for growth on phosphate (⁠$kPO4$⁠) is computed similarly as

The half saturation for growth on iron (⁠$kFe$⁠) is computed similarly as

The half saturation for growth on ammonia (⁠$kNH4$⁠) is computed with a factor of 1/2 because phytoplankton preferentially use ammonia:

The most limiting nutrient determines the value of the nutrient limitation for growth (γ^R, between 0 and 1) as

where $σNH4$is the coefficient for ammonia inhibition of nitrogen uptake and $RNO3$⁠, $RNO2$⁠, $RNH4$⁠, $RSi(OH)4$⁠, $RPO4$ and $RFe$ are the concentrations of nitrate, nitrite, ammonia, silicic acid, phosphate and iron, respectively.

Temperature modulation for growth (γ^T, between 0 and 1) is calculated from the ambient temperature (T) as in Dutkiewicz et al. (2015):

where τ_T is a normalization factor for temperature function, A_T is a constant and T_N is the reference temperature. The interval of ambient temperature in the forcing fields of the simulation was −1.8°C to 2.9°C. Following Equation 13, the interval of modification of growth rate by temperature was 0.27 to 0.34. The term γ^T was the same for all plankton types.

Light limitation for growth (γ^I, between 0 and 1) follows Geider et al. (1998) as described in Dutkiewicz et al. (2015):

where α^chl is a constant linear initial slope of the Chl a-specific photosynthesis versus irradiance curve in nutrient-replete conditions (Figure S1b), E₀ is the scalar PAR, θ is the Chl:C ratio, P^C_max is the maximum growth rate (Figure 2b), γ^R is the nutrient limitation for growth, γ^T is the temperature modulation for growth and M_C is the molar mass of carbon. The last two factors of Equation 14 are required to correctly convert the units of θ and α^chl. The production of new Chl a follows Geider et al. (1997) and Geider et al. (1998) as described in Dutkiewicz et al. (2015).

Grazing follows from Dutkiewicz et al. (2020) as

The subscript j is the prey, and the subscript k is the predator. The term g_jk is the grazing rate. The maximum grazing rates (g_max,k; Figure 2d) are from observations (Jeong et al., 2010; Taniguchi et al., 2014). The term σ_jk is palatability, B_j is prey biomass and G_k is the palatability-weighted total phytoplankton biomass as

The grazing half saturation rate (k_p) is a constant.

In this section, we provide a description of the calculation of the scalar PAR at each depth.  Appendix A2.1 describes the calculation of the scalar PAR just below the surface when the surface was sea ice before July 18.  Appendix A2.2 describes the calculation of the scalar PAR just below the surface when the surface was open water from July 18.  Appendix A2.3 describes the calculation of the scalar PAR at each depth from the scalar PAR just below the surface. Equations from  Appendix A2.3 were used the whole year. See Tables S6 and S7 for the units of the variables and the values of the constants, respectively.

The downwelling plane PAR just below sea ice in photon density flux, E_d,i(z = 0⁻, PAR[Q]) (Data S5 in Benoît-Gagné et al., 2024), was measured as described in Oziel et al. (2019) and Massicotte et al. (2020). The average cosine (μ_d) was used to calculate the downwelling scalar PAR just below sea ice in photon density flux, E_0d,i(z = 0⁻, PAR[Q]):

Following the observations of Matthes et al. (2019) at the Qikiqtarjuaq ice camp, μ_d was set to 0.6 under snow-covered sea ice and to 0.7 under bare sea ice. The observations of Matthes et al. (2019) also allowed the conversion of E_0d,i(z = 0⁻, PAR[Q]) into the scalar PAR just below sea ice in photon density flux, E_0,i(z = 0⁻, PAR[Q]) (Figure 8b; Data S5 in Benoît-Gagné et al., 2024), as

Downwelling shortwave radiation just above surface in energy units, E_s(z = 0⁺, SW) (Figure S10a), was from Smith et al. (2014). The processing of E_s(z = 0⁺, SW) followed that of the biogeochemical model PISCES (Aumont et al., 2015). The downwelling shortwave radiation just below open water in energy units, E_s(z = 0⁻, SW) (Figure S10b), was calculated with an albedo for open water, α_w = 0.066, as

The scalar PAR just below open water in energy units, E_0,w(z = 0⁻, PAR[W]) (Figure S10c), was calculated using a factor of 0.43 as

The visible spectrum was divided into three equal bands: blue, green and red. The scalar irradiances for each band just below open water in energy units, E_0,w,λ(z = 0⁻, PAR[W]) (Figure S10d), were considered equal as

The scalar irradiances for each band just below open water in energy units were converted into the scalar irradiances for each band just below open water in photon density flux E_0,w,λ(z = 0⁻, PAR[Q]) (Figure S10e), as

for λ = blue = 450 nm, λ = green = 550 nm and λ = red = 650 nm. E_0,w,λ(z = 0⁻, PAR[Q]) was added up to get the scalar PAR just below open water in photon density flux, E_0,w(z = 0⁻, PAR[Q]) (Figure S10f; Data S5 in Benoît-Gagné et al., 2024), as

E_0,i(z = 0⁻, PAR[Q]) was used for the forcing field for the scalar PAR just below surface in photon density flux, E₀(z = 0⁻, PAR[Q]), before sea ice break-up (before July 18). E_0,w(z = 0⁻, PAR[Q]) was used as the forcing field for E₀(z = 0⁻, PAR[Q]) from the sea ice break-up (from July 18). The water column was divided into 75 depth layers. The scalar PAR at depth z in photon density flux was computed with a Beer-Lambert law as

where K₀ was the diffuse vertical attenuation coefficient of scalar PAR. K₀ was dependent on Chl a concentration (Chl a) as

where K_w was light absorption for pure seawater and K_Chl was light absorption for Chl a. Coloured dissolved organic matter and detrital matter were not considered in the calculation of K₀.

The measured PAR was the downwelling plane PAR (E_d), while the model simulated the scalar PAR (E₀). Matthes et al. (2019) found a conversion factor of approximately 1.4 from plane to scalar PAR within the upper 20 m at the Qikiqtarjuaq ice camp. Before June 27, there were observations of Chl a only down to 40 m and extrapolation was necessary for the vertical integration to 100 m. Reaching 100 m was required after June 27 because a subsurface chlorophyll maximum was formed by that point. Extrapolation was achieved by considering that Chl a below 40 m was equal to Chl a at 40 m because, before June 27, Chl a remained unchanged with depth (Figure S2d). The depth of the mixing layer was measured during the ice camp only on June 23. The methods for the measurement of the equivalent mixed layer depth (h_BD; as in Randelhoff et al., 2017), a measure distinct from the depth of the mixing layer, are described in Oziel et al. (2019). The equivalent mixed layer depth was smoothed with a moving average of 7 days.

The relatively deep water depth of 360 m at the ice camp suggested that this site may be representative of offshore western Baffin Bay, even though it was located only a few kilometres from the coast. Data from the 2018 expedition of the Canadian Coast Guard Ship (CCGS) Amundsen was used to test this hypothesis. Six stations including one at the Qikiqtarjuaq ice camp location were sampled from July 13 to July 24, 2018 (Leg 2b; Figure S11). Salinity and temperature measurements were acquired with a conductivity-temperature-depth sensor (CTD SBD911plus, SeaBird Scientific). The temperature-salinity diagram for each station showed that the oceanographic conditions in 2018 at the ice camp location (Station 5) were like the oceanographic conditions at the offshore stations (Stations 1 to 3; Figure S12). Thus, we believe that the ice camp location could be considered representative of western Baffin Bay, allowing this modelling study to be representative of this area as well.

The materials and methods for the observation of carbon biomass at the ice camp Qikiqtarjuaq in 2016 is described in Grondin (2019) and Massicotte et al. (2020). Only a brief overview is presented here and illustrated in Figure S13. The 5 ml seawater samples were analysed with an Imaging FlowCytobot (Olson and Sosik, 2007; Sosik and Olson, 2007) manufactured by McLane®. The targeted size range was between 1 µm and 150 μm. Cells larger than 10 µm could be identified to a finer taxonomic resolution than cells between 3 µm and 10 µm due to an image resolution of approximately 3.4 pixels μm⁻¹. The Chl a in vivo fluorescence with an excitation laser at 635 nm triggered image acquisition. The resulting greyscale images were processed with a MATLAB (2013b) code (Sosik and Olson, 2007; https://github.com/hsosik/ifcb-analysis, accessed May 29, 2024). This code extracted the regions of interest and their associated features with a random forest algorithm. Each region of interest had 231 features (backscattering, Chl a, geometry, shape, symmetry, texture, etc.; see https://github.com/hsosik/ifcb-analysis/wiki/feature-file-documentation, accessed June 11, 2024). These regions and their features were then processed with the EcoTaxa application (Picheral et al., 2017) again using a random forest algorithm. The result was images annotated with one of the 35 taxonomic categories. The annotated images were converted to biovolumes (Moberg and Sosik, 2012). The biovolumes were converted to carbon (Laney and Sosik, 2014) using carbon-to-volume ratios (Menden-Deuer and Lessard, 2000). For the purposes of this study and evaluation of the model, the biomasses of the 35 taxonomic categories were binned into three functional groups: diatoms, mixotrophs and other nanophytoplankton (https://github.com/maximebenoitgagne/timing/blob/main/timing.ipynb, accessed June 11, 2024). Algaebase (Guiry and Guiry, 2022) helped in the classification of the IFCB categories into these groups (Table S8). The ‘other nanophytoplankton’ included strictly autotrophic phytoplankton. More than 99% of the mixotrophs group were dinoflagellates, and as such we refer to this group as ‘mixotrophic dinoflagellates’ in this study.

Model output and other data for the generation of the figures are openly accessible in the data directory at https://github.com/maximebenoitgagne/timing/tree/v0.1.4 (Benoît-Gagné et al., 2024). DOI: https://doi.org/10.5281/zenodo.13287110.

The physical model used here is available through http://mitgcm.org, and the generic ecosystem code used in this study is available through https://gitlab.com/jahn/gud. The specific modifications for the setup used here, the code for the new options to the ecosystem code, all parameters values, and the code to generate the figures are provided at https://github.com/maximebenoitgagne/timing/tree/v0.1.4 (Benoît-Gagné et al., 2024). DOI: https://doi.org/10.5281/zenodo.13287110. The code to generate the figures can be visualized by opening the file timing.ipynb. The code to generate the supplemental figures can be visualized by opening the file timing_supmat.ipynb.

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

• Figure S1. Photosynthetic parameters.

• Figure S2. Nutrient concentrations and chlorophyll a.

• Figure S3. Annual drift of nitrate.

• Figure S4. Accumulation rate.

• Figure S5. Sensitivity simulations by phytoplankton groups and size classes: EXP-1 prescribed minimum biomass.

• Figure S6. Sensitivity simulations by phytoplankton groups and size classes: EXP-2 with differing light under sea ice.

• Figure S7. Sensitivity simulations by phytoplankton groups and size classes: EXP-3 prescribed pre-bloom nitrate concentrations.

• Figure S8. Mean biomass and chlorophyll a (0–40m).

• Figure S9. Parameters controlling the half saturation for growth on nitrate for each numerical type.

• Figure S10. Processing of downwelling shortwave radiation just above the surface.

• Figure S11. Maps of Baffin Bay and the study site including the expedition Amundsen 2018.

• Figure S12. Observations of temperature and salinity during the expedition Amundsen 2018.

• Figure S13. Materials and methods for the observation of the carbon biomass at the ice camp Qikiqtarjuaq 2016.

• Table S1. Notations mentioned in the Main Text.

• Table S2. Features of the phytoplankton blooms at the Qikiqtarjuaq ice camps in 2015 and 2016.

• Table S3. Notations mentioned in the Appendix A1.

• Table S4. Constants and fixed parameters mentioned in the Appendix A1.

• Table S5. Coefficients for allometric scaling (unitless) mentioned in Appendix A1.

• Table S6. Notations mentioned in the Appendix A2.

• Table S7. Constants and fixed parameters mentioned in the Appendix A2.

• Table S8. Classification of the Imaging FlowCytobot (IFCB) taxonomic categories into the biogeochemical groups of the model.

We thank Joannie Ferland and Marie-Hélène Forget for organizing the Qikiqtarjuaq ice camps. We thank Virginie Galindo for snow and ice thickness, Guislain Bécu and Marcel Babin for PAR, Patrick Raimbault for nitrate and silicic acid concentrations, Joséphine Ras and Hervé Claustre for high-performance liquid chromatography analyses, Joannie Ferland, Pierre-Luc Grondin and Marcel Babin for IFCB analyses and all other scientists and technicians who contributed to fieldwork. MBG thanks Jean-Michel Campin, Oliver Jahn and Jeffery Scott for their help in mastering the code of MITgcm. MBG thanks Achim Randelhoff for useful code examples (https://github.com/poplarShift/ice-edge). MBG thanks Yannick Copin for a useful code example of a Taylor diagram (https://gist.github.com/ycopin/3342888). MBG thanks Arthur Plassart and Sarah Schembri for suggestions on the manuscript.

This research was funded by the Sentinel North program of Université Laval, funded by the Canada First Research Excellence Fund. The GreenEdge project was conducted under the scientific coordination of the CERC program and of the CNRS & Université Laval Takuvik International Research Laboratory (LRI3376). The GreenEdge project was funded by the following French and Canadian programs and agencies: ANR (Contract #111112), ArcticNet, CERC on Remote sensing of Canada's new Arctic frontier, CNES (project #131425), French Arctic Initiative, Fondation Total, CSA, LEFE and IPEV (project #1164). The Amundsen 2018 cruise was conducted using the Canadian research icebreaker CCGS Amundsen with the support of the Amundsen Science program funded by the Canada Foundation for Innovation (CFI) Major Science Initiatives Fund. The supercomputing infrastructure was provided by the Digital Research Alliance of Canada. MBG acknowledges Sentinel North for a doctoral scholarship. MBG acknowledges the Fonds de Recherche du Québec and the Sentinel North program for funding a research appointment at the Massachusetts Institute of Technology, Cambridge, MA, USA. MBG acknowledges the Fonds de Recherche du Québec for funding a research appointment at the Laboratoire des Sciences de l'Environnement Marin (CNRS), Plouzané, France. This work is a contribution to the scientific program of the interinstitutional research group Québec-Océan. CD was supported by a grant from the Weston Family Foundation awarded to FM. SD acknowledges the support of NASA Earth Science Division’s Interdisciplinary Science (IDS) program (grant number 80NM0018D0004). DD was supported by NSERC Discovery Grants (RGPIN-2013-402257, RGPIN-2019-06563). RL acknowledges the support of the SMAART program through the Collaborative Research and Training Experience program (CREATE) of the NSERC. FM was supported by NSERC Discovery Grants (RGPIN-2014-05433, RGPIN-2021-03876).

The authors have declared that no competing interests exist.

Conception and design: MBG, SD, FM.

Analysed data: MBG.

Interpreted data: MBG, SD, ID, CD, DD, LM, FM.

Acquisition and processing of temperature and salinity data during the CCGS Amundsen 2018 expedition: RL.

Processing of nutrient concentrations before May 15 from the Qikiqtarjuaq sea ice camps 2015 and 2016: GO.

NEMO-LIM3 modelling to provide offline physical forcing fields (temperature, salinity and mixing) for the MITgcm modelling: GO.

MITgcm modelling: MBG.

Manuscript preparation: MBG, RL, GO.

Critical review and revision of the manuscript: MBG, SD, ID, CD, DD, LM, FM.

Final approval of the versions to be submitted: MBG, SD, ID, CD, DD, RL, LM, GO, FM.

How to cite this article: Benoît-Gagné, M, Dutkiewicz, S, Deschepper, I, Dufresne, C, Dumont, D, Larouche, R, Mémery, L, Olivier, G, Maps, F. 2024. Exploring controls on the timing of the phytoplankton bloom in western Baffin Bay, Canadian Arctic. Elementa: Science of the Anthropocene 12(1). DOI: https://doi.org/10.1525/elementa.2024.00008

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