RESEARCH ARTICLE Assessing the Role of High-Frequency Winds and Sea Ice 10.1029/2018JG004869 Loss on Arctic Phytoplankton Blooms in an Key Points: Ice-Ocean-Biogeochemical Model • High-frequency winds stimulate surface mixing and nutrient upwelling enhancing primary L. Castro de la Guardia1 , Y. Garcia-Quintana1 , M. Claret2,3 , X. Hu1,4 , production and biogenic carbon E. D. Galbraith2,5,6 , and P. G. Myers1 export • The model shows significant iron 1 Department of Earth and Atmospheric Sciences, University of Alberta, Edmonton, Alberta, Canada, 2 Department of limitation in the Labrador Sea, which is intensified by high-frequency Earth and Planetary Sciences, McGill University, Quebec, Montréal, Canada, 3 Joint Institute for the Study of the winds Atmosphere and the Ocean, University of Washington, Seattle, WA, USA, 4 Bedford Institute of Oceanography, • Sea ice loss contributes to increasing Department of Fisheries and Oceans, Dartmouth, Nova Scotia, Canada, 5 Institució Catalana de Recerca i Estudis primary production only in Avançats, Barcelona, Spain, 6 Institut de Ciència i Tecnologia Ambientals, Universitat Autònoma de Barcelona, Hudson Bay and in the presence of high-frequency winds Cerdanyola del Vallès, Spain Correspondence to: Abstract The long-term trend of increasing phytoplankton net primary production (NPP) in the Arctic L. Castro de la Guardia, correlates with increasing light penetration due to sea ice loss. However, recent studies suggest that castrode@ualberta.ca enhanced stormy wind mixing may also play a significant role enhancing NPP. Here, we isolate the role of sea ice and stormy winds (hereafter high-frequency winds) using an eddy-permitting Citation: Castro de la Guardia, L., Garcia‐ ice-ocean-biogeochemical model configured for the North Atlantic and the Arctic. In the model, the Quintana, Y., Claret, M., Hu, X., presence of high-frequency winds stimulates nutrient upwelling by producing an earlier and longer Galbraith, E. D., & Myers, P. G. (2019). autumn-winter mixing period with deeper mixing layer. The early onset of autumn mixing results in Assessing the role of high‐frequency winds and sea ice loss on Arctic nutrients being brought-up to near-surface waters before the light becomes the dominant limiting factor, phytoplankton blooms in an ice‐ocean‐ which leads to the autumn bloom. The enhanced mixing results in higher nutrient concentrations in biogeochemical model. Journal of spring and thus a large spring bloom. The model also shows significant iron limitation in the Labrador Sea, Geophysical Research: Biogeosciences, 124, 2728–2750. https://doi.org/ which is intensified by high-frequency winds. The effect of sea ice loss on NPP was found to be regionally 10.1029/2018JG004869 dependent on the presence of high-frequency winds. This numerical study suggests high-frequency winds play significant role increasing NPP in the Arctic and sub-Arctic by alleviating phytoplankton nutrient Received 12 OCT 2018 limitation and that the isolated effect of sea ice loss on light plays a comparatively minor role. Accepted 26 JUL 2019 Accepted article online 23 AUG 2019 Plain Language Summary The amount of phytoplankton (i.e., small floating algae) in the Published online 5 SEP 2019 Arctic Ocean has been increasing in response to sea ice losses, which increases light penetration and Author Contributions reduces light limitation on growth. However, the parallel increase in storm activity can also enhance Conceptualization: L. phytoplankton growth by reducing nutrient limitation through increasing surface ocean mixing. Here, we Castro de la Guardia use a numerical ocean—sea ice model with marine biogeochemistry to separate the contribution of stormy Funding Acquisition: P. G. Myers Methodology: L. winds and sea ice loss on the phytoplankton growth. In the model, stormy winds induce more ocean Castro de la Guardia, Y. mixing. This results in a seasonally larger amount of nutrient in surface waters with the capacity of Garcia-Quintana, M. Claret, X. Hu, E. sustaining a larger phytoplankton bloom in spring and autumn if light conditions are optimal for growth. D. Galbraith, P. G. Myers Validation: L. Castro de la Guardia The model shows a significant iron limitation in the Labrador Sea, which is intensified by high-frequency Writing - Original Draft: L. winds. The contribution of sea ice loss to phytoplankton growth was only significant in the presence of Castro de la Guardia stormy winds. This numerical study suggests high-frequency winds play significant role increasing Formal Analysis: L. Castro de la Guardia, Y. phytoplankton growth in the Arctic and sub-Arctic by alleviating phytoplankton nutrient limitation and Garcia-Quintana, M. Claret that the isolated effect of sea ice loss on light plays a comparatively minor role. Investigation: Y. Garcia-Quintana Resources: P. G. Myers Supervision: E. D. Galbraith, P. G. Myers 1. Introduction Visualization: L. Castro de la Guardia Accelerated loss of Arctic sea ice in summer (Comiso et al., 2008, 2017; Serreze et al., 2007) can increase Writing - review & editing: Y. the Arctic Ocean's primary production by enhancing the amount of light available at the ocean surface Garcia-Quintana, M. Claret, X. Hu, E. (Arrigo et al., 2008; Pabi et al., 2008). This sea ice loss, however, is also accompanied by increasing freshwater D. Galbraith, P. G. Myers input from sea ice melt. While the former decreases light limitation on phytoplankton productivity, the later ©2019. American Geophysical Union. increases the nutrient limitation. The input of buoyant waters stratifies the water column, which inhibits the All Rights Reserved. upward transport of nutrients, limiting phytoplankton growth in the well-lit surface ocean waters (Harrison CASTRO DE LA GUARDIA ET AL. 2728 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 & Cota, 1991; McLaughlin & Carmack, 2010). This stratification is amplified in summer due to surface warming. Hence, a year with high sea ice loss may need to be accompanied by a source of energy to erode stratification and resupply nutrients to upper ocean layers, in order to result in high productivity. Winds are one of the primary sources of energy for upper ocean mixing (Huang & Qiao, 2009). Arctic sea ice losses, along with increased heat and moisture fluxes, have enhanced storm genesis within the Arctic during summer (Day et al., 2017; Deser et al., 2000; Simmonds & Rudeva, 2012). Therefore, the increase in phytoplankton productivity in recent years may not only be the result of increased light at the ocean surface but also be the result of increasing nutrient supply by more intense wind-driven mixing. Yang et al. (2004) using buoy data collected in winter and summer within the Arctic Ocean showed that intense Arctic storms can lead to abnormally deep mixing events, which can temporally break the Arctic halocline stratification and lead to an exchange of heat and salt with the mixed layer. Stormy wind events (i.e., wind speed ≥10 m/s) over the Arctic waters are indeed increasing (Yang et al., 2004) and are correlated with higher autumn chlorophyll-a (chl-a) concentration (Ardyna et al., 2014). Using satellite observations, Ardyna et al. (2014) proposed that the wind-driven mixing breaks the summer strat- ification and brings-up nutrients into the mixed layer, which in turn increases the chl-a concentration in autumn. Yet they were unable to confirm the exact mechanisms by which the wind-induced mixing impacts chl-a concentration because satellite observations are limited to the surface ocean. In year-round ice-free waters, large wind speeds (≥6 m/s) can indeed cause substantial macronutrient entrainment, leading to an increase in surface phytoplankton growth (Lin et al., 2003; Pan et al., 2017; Rumyantseva et al., 2015; Wu et al., 2007; Yin et al., 1996). In contrast, in seasonally ice-covered waters, a correlation between stormy wind events and chl-a is not necessarily conclusive of higher phytoplankton productivity. In particular, Carranza and Gille (2015) pointed out that in the Southern Ocean, increases in chl-a concentration after a wind-driven mixing event may be the result of chl-a from the subsurface chl-a maximum being brought-up into the mixed layer, and not necessarily from a local enhancement on phy- toplankton growth as a result of a nutrient-stress relief. Hence, the role of stormy wind events on the net annual productivity of seasonally ice-covered regions of the Arctic Ocean still remains ambiguous. The present study uses two numerical experiments described in section 2 and evaluated in Appendix A to quantify the role of stormy wind events and sea ice losses on the phenology and magnitude of phytoplankton blooms and the net annual primary production and biogenic carbon export over six basins in the Arctic and sub-Arctic. We suggest that high-frequency winds contribute to the development of the autumn bloom through stimulating an early onset of autumn mixing that brings-up nutrients and releases phytoplankton from the summer nutrient limitation (section 3.1–3.2). We also observe that in our simulations the loss of sea ice without the enhanced wind-driven mixing by high-frequency winds does not result in higher phytoplankton productivity (section 3.3). Finally, we offer a discussion of these results in section 4 and state our conclusions in section 5. 2. Methods Arctic and subarctic regions can be characterized according to the seasonality of the sea ice cover as ice-free, seasonally ice-covered, or perennially ice-covered regions. Here, we define ice-free regions as having less than 25% sea ice concentrations year-round; seasonally ice-covered regions as having sea ice concentrations less than 25% for a minimum of 3 months, but more than 40% the rest of the year; and perennially ice-covered regions as having sea ice concentration above 40% year-round. This study has a focus on the Canadian side of the Arctic, and based on these definitions, we chose two representative regions for each ice regime. For ice-free regions, we chose the Labrador Sea and Barents Sea; for seasonally ice-covered regions Baffin Bay and Hudson Bay; and for perennially ice-covered regions, the Beaufort Gyre and the Canadian Arctic Archipelago (Figure 1). 2.1. Numerical Model We simulate the physics using a coupled ocean sea ice system based on the Nucleus for European Mod- elling of the Ocean (NEMO, https://www.nemo-ocean.eu) Version 3.4 (Madec & the NEMO team, 2008). The sea ice module used is the Louvain-la-neauve Ice Model Version 2 (Vancoppenolle et al., 2009), includ- ing both thermodynamic and dynamic components (Fichefet & Maqueda, 1997; Hunke & Dukowicz, 1997). CASTRO DE LA GUARDIA ET AL. 2729 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 1. Selected Arctic regions. The selection is based on simulated sea ice concentration climatology 2002–2015. (a) A map identifying the six regions selected. (b) Simulated sea ice fraction at minimum sea ice extent (i.e., September). (c) Simulated sea ice fraction at maximum sea ice extent (i.e., March). Regions can be characterized according to the seasonality of the sea ice cover as ice-free (Labrador Sea and Barents Sea), seasonally ice-covered (Hudson Bay and Baffin Bay), and perennially ice-covered (Beaufort Gyre and Canadian Arctic Archipelago). The detailed configuration on sea ice and ocean components are documented in Hu et al. (2018) and Garcia-Quintana et al. (2019). The marine biogeochemical model is the Biogeochemistry with Light Iron and Nutrient limitation and Gases (BLING) Version 0 (Galbraith et al., 2010), which is a new addition to the NEMO modeling framework. BLING version 0 is a reduced complexity biogeochemical model with four prognostic tracers: inorganic phosphate, dissolved organic phosphate, oxygen, and iron. It diagnoses chl-a, phytoplankton production, and particle export considering light, macronutrient, and iron limitations as well as a temperature depen- dency. BLING version 0 does not explicitly model the nitrogen cycle and therefore does not explicitly represent the loss of nitrate by denitrification. Although denitrification has been identified as an important nitrogen sink process in the Beaufort Sea and the Arctic shelves (Devol et al., 1997; McTigue et al., 2016), the similarities between modeled and observed macronutrient concentration (see Appendix A.3) suggest that the lack of explicit denitrification within the domain does not exert a large overall influence on the model results. The benefit of using BLING version 0 over other more complex biogeochemical models is that it reduces computational costs, which represents an advantage when integrating high-resolution simulations over several decades. 2.1.1. Model Configuration: ANHA4 We use a regional configuration of NEMO for the Arctic and the North Atlantic (ANHA4) with a nomi- nal horizontal resolution of 0.25◦ . The three-dimensional primitive equations are discretized on a tripolar grid and 50 vertical levels using the z coordinate. This configuration has two open boundaries: one near the Bering Strait in the Pacific Ocean and the other one at 20◦ S across the Atlantic Ocean. At the open bound- aries, there is an exchange of heat, momentum, and tracers, while within the interior the model evolves without restoring. The ocean module represents the dynamics and thermodynamics of the ocean using the spherical Earth approximation, the thin-shell approximation, and the hydrostatic approximation. We parameterize the ver- tical subgrid-scale mixing with a 1.5-th order turbulent kinetic energy closure scheme that consists of a prognostic equation for turbulent kinetic energy (TKE) with a length scale closure assumption to diagnose the length scales for turbulent mixing and dissipation. This scheme uses a background value of 10−5 m2 /s for CASTRO DE LA GUARDIA ET AL. 2730 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 the vertical eddy diffusivity (Kz ), which is dependant on the Richardson number (Ri) and on sea ice concen- tration. That is, Kz increases in a region where Ri becomes critical and it is reduced under sea ice cover (i.e., as linear function of ice concentration to a minimum of 1/10 of the regular value under a 100% ice cover). In ANHA4, we do not use the Gent-McWilliams parametrization because the configuration is eddy permitting. 2.1.2. Initial Conditions We initialized the model on January 1958 using climatological ocean temperature, salinity from the Polar Science Center Hydrographic Climatology 3.0 (Steele et al., 2001), with the ocean at rest and a constant sea ice cover of 3-m thickness in all grid cells where the ocean temperature was close to the freezing point. Bio- logical fields were initialized using a combination of observed climatologies and model output. For dissolved oxygen and inorganic phosphate fields we used climatologies from the World Ocean Atlas 2013 Version 2 (WOA13; Garcia et al., 2014a, 2014), and for dissolved iron and organic phosphate fields we used climatolo- gies from model output of the Geophysical Fluid Dynamics Laboratory (GFDL) Earth System Model Version 2 with z-based vertical coordinate (ESM2M) coupled with BLING Version 0 (Galbraith et al., 2015). The lat- ter simulation is a global configuration at 1◦ nominal resolution that ran for 100 years. The climatologies were built from the last 20 years. 2.1.3. Boundary Conditions For temperature, salinity, and horizontal velocities we used model output from two global simulations of NEMO at 0.25◦ resolution. Between January 1958 and December 2001 we used the output from a simulation run by GEOMAR, Helmholtz Center for Ocean Research Kiel (ORCA025-K3415; http://www.geomar.de/). Between January 2002 and December 2015 we used the output from the Global Ocean Reanalysis and Sim- ulations (GLORYS2v3; Masina et al., 2015). Boundary conditions for biological tracers were derived from monthly climatologies built from GFDL-ESM2M. All boundary conditions were interpolated on the fly to the model time step of 18 min. A source of iron at the surface ocean was added following the relation between dust deposition and iron concentrations described in Galbraith et al. (2010). The climatological monthly dust deposition input at the surface of the ANHA4 domain was derived from the Global Ozone Chemistry Aerosol Radiation and Transport model (Ginoux et al., 2001). To represent freshwater sources from the land, we used the monthly interannual river runoff of Dai et al. (2009) and the Greenland meltwater runoff of Bamber et al. (2012). These sources were remapped onto ANHA4 grid to give a more realistic freshwater input to the ocean, and they did not contain nutrients or dissolved organic products. Runoff was treated as zero salinity freshwater with a temperature equal to the local sea surface temperature. 2.1.4. Atmospheric Forcing Data From January 1958 to December 2001 we forced the model with the reanalysis product from the Coordi- nated Ocean-ice Reference Experiment version 2-Interannual Atmospheric Forcing (COREv2-IAF; Large & Yeager, 2009). This data set consists of interannually varying monthly precipitation, daily shortwave and longwave radiation, and 6-hourly air temperature, specific humidity, and zonal and meridional wind velocities. All 6-hourly data are taken at 10 m above ground. From January 2002 to December 2015, we switch to the high spatial and temporal resolution forecast-reanalysis data from the Canadian Meteorological Centers Global Deterministic Prediction System (CMC-GDPS; Smith et al., 2014). This data set includes hourly shortwave and longwave radiation fluxes, precipitation, wind velocities at 10 m above ground, and air temperature and specific humidity at 2 m above ground. 2.2. Setup of High-Frequency Winds Sensitivity: CONTROL and CALM We ran a spinup simulation with COREv2-IAF for 43 years (starting in January 1958) in order to minimize the sensitivity of our results to the initial conditions. Starting in January 2002 until December 2015, two twin simulations branched out from this spinup, differing only in the atmospheric forcing. One simulation (CONTROL) was forced with CMC-GDPS atmospheric forcing. A second simulation (CALM) was forced with a low-pass time-filtered version of the CMC-GDPS atmospheric forcing. The low-pass time-filtered was applied to the CMC-GDPS fields of wind and air temperature. The filtered forcing was created by applying a running mean with a 10-day window length to the 2-m air temperature and the zonal and meridional components of the 10-m winds following the methodology in Holdsworth and Myers (2015) and Garcia-Quintana et al. (2019). This window length captures the mean lifetime of CASTRO DE LA GUARDIA ET AL. 2731 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 2. Example of the wind and temperature forcing in the original Canadian Meteorological Centers-Global Deterministic Prediction System (red) and the low-pass-filtered Canadian Meteorological Centers-Global Deterministic Prediction System (green) in the Labrador Sea. Time series of (a) daily mean wind speed, (b) maximum wind speed, and (c) daily maximum, minimum, and mean temperature in 2012. (d and e) Autumn climatology (October–December 2002–2015) of regional wind speed pattern. Original fields (red) are used to force the CONTROL simulation while filtered fields (green) are used to force the CALM simulation. This example focuses on the Labrador Sea as defined in Figure 1, but the same pattern was found in all six basins (not shown). Arctic storms, which is commonly 3 days, ranging from 1 to 13 days (Simmonds & Rudeva, 2012). This low-pass filter has been proven effective at removing the high-frequency atmospheric phenomena in the North Atlantic (Garcia-Quintana et al., 2019; Holdsworth & Myers, 2015). As a result, the CALM simulation does not include high-frequency atmospheric phenomena such as stormy winds and cold-air outbreaks. The daily mean and maximum wind speed of the filtered forcing (green in Figures 2a and 2b) are signif- icantly reduced in comparison to the original CMC-GDPS fields (red in Figures 2a and 2b) but the wind speed pattern of direction is preserved (Figures 2d and 2e). Noticeably, most of the winds with wind speed ≥10 m/s, which are commonly associated with transitory storms (Wang et al., 2017), are smoothed out by the low-pass filter. The temperature peaks associated with high-frequency atmospheric phenomena are also reduced, but the daily mean temperature is not significantly affected by the filter (Figure 2c) in agreement with Holdsworth and Myers (2015). Hence, with the purpose of highlighting the difference in the atmo- spheric forcing between the CONTROL and CALM simulations, we hereafter adopt “high-frequency winds” to refer the “high-frequency atmospheric phenomena,” which are absent in the CALM simulation. We note that filtering the atmospheric temperature and wind components did not affect the sea ice concen- tration nor the ice-free period in our study regions, except for the Labrador Sea. In the Labrador Sea, the absence of high-frequency winds increased the ice concentration along the ice edge, but it did not affect the length of the ice-free season. Furthermore, most of this difference occurred during winter, which does not interfere with phytoplankton growth as phytoplankton is light-limited during this period. A seasonal com- parison of the light limitation in the CALM and CONTROL simulation (more details in section 3.1), revealed that light limitation remains similar in both simulations for most of the year, with a stronger light limitation in spring and autumn in the CONTROL simulation driven by a deeper mixing layer depth. CASTRO DE LA GUARDIA ET AL. 2732 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 3. Selected years with anomalously more (filled squares) and less sea ice (filled circles) in CONTROL (red) and CALM (green) simulations. In perennially ice-covered regions, (a) Beaufort Gyre and (b) the Canadian Arctic Archipelago, sea ice loss is measured using the anomaly of the open water area. In seasonally ice-free regions, (c) Baffin Bay and (d) Hudson Bay, sea ice loss is measured using the anomaly of ice-free period. Years with anomalously more (less) sea ice are one standard deviation (dashed lines) above (below) the anomaly mean. Anomalies are relative to the 2002–2015 climatology. For the geographical locations of these regions please refer to Figure 1a. BB = Baffin Bay; BG = Beaufort Gyre; CAA = Canadian Arctic Archipelago; HB = Hudson Bay. 2.3. Setup of Sea Ice Loss Sensitivity: Climatologies for Less and More Sea Ice Sensitivity to sea ice loss was addressed in perennially and seasonally ice-covered regions by comparing years with less sea ice to years with more sea ice. In perennially ice-covered regions (Figures 3a and 3b), sea ice loss was measured as a function of sea ice cover by comparing years with anomalously large open water area (less sea ice) to years with anomalously small open water area (more sea ice). In seasonally ice-covered regions (Figures 3c and 3d), sea ice loss was measured as a function of the length of the ice-free period by comparing years with an anomalously long ice-free season (less sea ice) to years with an anomalously short ice-free season (more sea ice). Climatologies with “more” and “less” sea ice were constructed as follows. First, sea ice loss anomalies were computed by subtracting the climatology (2002–2015) of sea ice cover or ice-free season length. Years with anomalously more (less) sea ice were selected to be one standard deviation below (above) the anomaly mean (dashed line in Figure 3). A climatology for more (less) sea ice was finally obtained by averaging at least 3 years with anomalously high (low) sea ice conditions. Exceptions were made when less than 3 years meet this criterion (i.e., we selected years which anomaly was closest to the standard deviation). In both simu- lations, the climatology with less sea ice has a longer ice-free period (in Hudson Bay: 18% longer in CALM and 16.2% longer in CONTROL; in Baffin Bay: 16.6% longer in CALM, and 11.4% longer in CONTROL) and a larger open water area (in the Canadian Arctic Archipelago: 27.9% larger in CALM, and 16.6% larger in CONTROL; and in the Beaufort Sea: 11.1% larger in CALM, and 11.4% larger in CONTROL) relative to the climatology with more sea ice. 2.4. Model Output Analyses 2.4.1. Estimation of Turbulence The turbulent regime is evaluated in the model using the Richardson number (Ri), defined as the ratio between the squared buoyancy frequency, N, and the squared vertical shear, Vz ; that is, Ri ≡ N 2 ∕Vz2 . The buoyancy frequency is indicative of the stratification of the water column and is computed in consistency with the formulation used in NEMOv3.4, based on Jackett and Mcdougall (1995): N2 ≡ g (𝛼𝜕 z T − 𝛽𝜕 z S). CASTRO DE LA GUARDIA ET AL. 2733 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Where g = 9.8 m/s2 is the gravity acceleration, 𝛼 = 2 × 10−4 ◦ C−1 is the thermal expansion coefficient, T is the potential temperature in degrees Celsius, 𝛽 = 7.7 × 10−4 is the saline expansion coefficient, and S is salinity (unitless). The operator 𝜕 z is the partial derivative with respect to the vertical. The vertical shear is indicative of the capacity of the water to become turbulent and is computed as Vz ≡ |𝜕z u⃗ h |, where u⃗ h is the horizontal velocity vector. We use a conventional definition to define a turbulent regime, that is, when the Ri number reaches a value below a critical threshold of 0.25, or in other words, when the vertical shear is twice the buoyancy frequency. 2.4.2. Estimation of Obduction We define obduction of macronutrients (O, in units of moles of phosphate per square meter per second) as the area-weighted average flux of dissolved macronutrients into the mixing layer due to physical transport. We derived O following Karleskind et al.'s (2011) and Levy et al.'s (2013) equation for the subduction of biogeochemical tracers, namely, [ ( ) ( ) ] 1 𝜕b 𝜕C O= −Cb − Cb u⃗ b · 𝛥⃗h b − Cb wb + Kz dA, (1) AR ∫A 𝜕t 𝜕z b where subscript b denotes values at the base of the mixed layer, which depth is b, C is any given nutrient concentration, u⃗ is the model horizontal velocity, w is the model vertical velocity, Kz is the vertical diffusiv- ity coefficient, AR is the regional area, and A is the area of each model grid cell. We note that the vertical dimension z points upward, such that positive values of O (obduction) indicate a supply of nutrients into the mixing layer, and negative values of O (subduction) indicate that nutrients leave the mixing layer. Each term on the right-hand side of equation (1) describes a different physical process that contributes to the exchange of macronutrients through the base of the mixing layer. The first term describes the flux of macronutrients due to the seasonal fluctuation of the mixing layer depth (i.e., entrainment). The second term describes the lateral transfer across the sloping base of the mixing layer (i.e., horizontal advection). The third term describes the vertical transfer across the base of the mixing layer by the vertical velocity (i.e., vertical advection). The fourth term describes the vertical diffusive flux as a result of a macronutrient concentration gradient at the base of the mixing layer. We computed the obduction value at each grid cell and then computed the area-weighted averages for each region. The vertical diffusion coefficient, Kz , varies in space and time as a function of Ri and sea ice concentration. When the water is turbulent (Ri ≤ 0.25), there is enhanced diffusion (Kz = 10−4 m2 /s); otherwise, Kz = 10−5 m2 /s. In addition, the background Kz is reduced as a linear function of sea ice concentration toward a minimum (Kz = 10−6 m2 /s) when the grid cell has 100% ice cover. This definition of Kz follows NEMO's turbulent kinetic energy mixing scheme for subgrid-scale vertical processes (Madec & the NEMO team, 2008). By using Kz = 10−5 m2 /s as the background diffusivity we may induce higher mixing in strongly stratified regions of the ocean where vertical diffusivity is several orders of magnitude lower Kz ∼ 10−6 m2 /s (Waterhouse et al., 2014). 3. Results 3.1. Seasonal Phytoplankton Growth Limitation Is Sensitive to the Presence of High-Frequency Winds Phytoplankton growth in the BLING model is limited by light, temperature, and nutrients. Seasonal changes that affect ambient nutrient concentration, sea water temperature, and photosynthetically available radia- tion (PAR) modulate which variable limits the growth during the course of a year. In Figure 4 we normalized the growth limitation using the seasonal maximum, such all three limiting factors range from 0 to 1. Three lines are depicted for each simulation, nutrient limitation (circles), light limitation (triangles), and tem- perature limitation (squares). The marker color represent the two simulation: CONTROL (red) and CALM (green). The line closest to the yellow shade identifies the factor that limits phytoplankton growth in each simulation. Figure 4 also depicts the sea ice (blue shade) and light (white-gray shade) regime of each region. In our six northern regions, phytoplankton growth is strongly limited by light from late autumn to early spring (e.g., line with triangle closest to yellow shade). The exact length of the period of light limitation varies with the geographic location and the regional ice conditions. For example, in the Beaufort Gyre at 72◦ N and CASTRO DE LA GUARDIA ET AL. 2734 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 4. Seasonal phytoplankton growth limitation by light (filled triangles), temperature (filled squares), iron (filled circles), and phosphate (unfilled circles) in CONTROL (red) and CALM (green) simulations. Sea ice fraction in CONTROL (blue area) and in CALM (dashed blue line) are plotted above each region. Light regime (gray to white area) is the solar radiation at ocean/sea ice surface in the CONTROL simulation. Five-daily climatology (2002–2015) for the six regions selected in Figure 1a. BB = Baffin Bay; BG = Beaufort Gyre; BS = Barents Sea; CAA = Canadian Arctic Archipelago; DIM = dissolved inorganic macronutrient; HB = Hudson Bay; IC = ice concentration; LS = Labrador Sea. perennial sea ice conditions, light limits growth 10 months of the year; meanwhile, in the Labrador Sea at 61◦ N and ice-free conditions light limits growth only 8 months of the year. Summer growth limitation is influence by temperature (only in Hudson Bay and early in summer) and nutrients (all six regions). In spring and autumn in ice-free and seasonally ice-covered regions, the phytoplankton growth in the CON- TROL simulation is more light limited than in the CALM simulation (line with red triangles is above line with green triangles in Figure 4). This is the result of BLING averaging the PAR over the depth of the mix- ing layer. Therefore, PAR in surface water is more attenuated in the CONTROL simulation, which has the deeper mixing layer (not shown). In ice-free regions, nutrient limitation is strongest in the CALM relative to the CONTROL simulation throughout the year, while in seasonally ice-covered and perennially ice-covered regions nutrients limit growth only in autumn. The lower nutrient limitation in the CONTROL simulation is the result of higher nutrient replenishment to surface waters due to higher surface turbulence. Additionally, while the nutrient limitation in summer is a result of macronutrient limitation in most basins (open circles in Figure 4), in the Labrador Sea the nutrient limitation results of iron limitation (filled circles in Figure 4). Furthermore, this iron limitation in the Labrador Sea is more pronounced in the CONTROL simulations than in the CALM simulation (more red filled circles than green filled circles in Figure 4b). In our model, the major sources of iron are dust deposition and diffusion from the seafloor, while the major sink of iron is scavenging by sinking particles (Galbraith et al., 2010). The Labrador Sea is a deep basin such that the seafloor iron supply is small compared to shallower basins, and dust input is low. In addition, deep mixing in the Labrador Sea produces a large input of macronutrients to the surface ocean, supporting large particulate export, which results in rapid removal of iron from the water column via scavenging. Thus, the weak iron source, combined with large scavenging rate, leads to low iron-to-macronutrient ratio in waters of the Labrador Sea. This iron limitation is further strengthened in the presence of high-frequency winds because by enhancing the macronutrient supply, high-frequency winds increase the particulate organic export, which leads to higher iron scavenging rates. In summary, iron limitation occurs in surface waters of the Labrador Sea due to the combination of a low input of iron to the basin with rapid macronutrient resupply and consequent export-driven scavenging. CASTRO DE LA GUARDIA ET AL. 2735 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 5. Averaged (2002–2015) NPP integrated from the surface to 50 m deep and biogenic carbon export at 50 m. Maps in (a) and (b) show the percentage difference between the CONTROL and the CALM simulations per model gridpoint (positive meaning higher in CONTROL). Bar plots in (c) and (d) show the regional average NPP and biogenic carbon export in CONTROL (red) and CALM (green) simulations. The percentage value over each bar highlights the contribution of the high-frequency winds relative to the CALM. Regions are identified in Figure 1a. Error bars in (c) and (d) show the standard deviation. Asterisks denote that CONTROL and CALM estimates are significantly difference within the 95% confidence interval (*) or within the 99% confidence interval (**) using a Student's test. BB = Baffin Bay; BG = Beaufort Gyre; BS = Barents Sea; CAA = Canadian Arctic Archipelago; HB = Hudson Bay; LS = Labrador Sea; NPP = net primary production. 3.2. High-Frequency Winds Increase the Mean Annual Phytoplankton Productivity and Biogenic Carbon Export High-frequency winds increase the annual mean phytoplankton net primary production (NPP) and bio- genic carbon export throughout the Arctic and sub-Arctic (positive values in Figures 5a and 5b), except in some parts of the southern Beaufort Gyre and the central Canadian Arctic Archipelago (negative values in Figures 5a and 5b). The strongest positive response to high-frequency winds is in the eastern and central Arctic (>50% increase in CONTROL relative to CALM; Figures 5a and 5b). In ice-free and seasonally ice-covered regions, the high-frequency winds account for 30–46% of the mean annual primary production and carbon export. The highest contribution of high-frequency winds is in the seasonally ice-covered region Hudson Bay (HB in Figures 5c and 5d), where the NPP has increased by 46% and the biogenic carbon export by 41% (red bar) relative to the CALM simulation (green bar). In perennially ice-covered regions, the Canadian Arctic Archipelago and Beaufort Gyre, the response is spa- tially dependent, with parts of these regions having lower primary production and biogenic carbon export in the CONTROL simulation (negative values (blue) in Figures 5a and 5b). Basin-average values in the Cana- dian Arctic Archipelago show that the high-frequency winds increase the annual mean primary production by 17% and biogenic carbon export by 11%. In contrast, basin-average values in the Beaufort Gyre show no net impact of the high-frequency winds on the mean annual net primary productivity and carbon export. 3.2.1. High-Frequency Winds Enhance the Phytoplankton Bloom in Spring and Autumn The larger net annual primary production and carbon export in the presence of high-frequency winds is primarily supported by a moderate enhancement of the phytoplankton productivity in spring (Figures 6a–6d), accompanied by a significant enhancement in autumn (Figures 6e–6h). The positive impact of high-frequency winds is subject to large spatial variations including the negative influence of high-frequency winds in some areas of the Canadian Arctic Archipelago and the Beaufort Gyre (Figures 6a, 6b, 6e, and 6f). In ice-free and seasonally ice-covered regions, high-frequency winds increase the NPP and CASTRO DE LA GUARDIA ET AL. 2736 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 6. Spring (March–June) and autumn (September–December) climatology (2002–2015) of NPP integrated from the surface to 50 m depth and biogenic carbon export at 50 m. Maps show the differences between the CONTROL and the CALM simulations per model grid-point in percentage (positive is higher in CONTROL) in (a and b) spring and (e and f) autumn. Bar plots show the regional (c and d) spring and (g and h) autumn climatology of NPP and biogenic carbon export in the CONTROL (red) and CALM (green). Regions are identified in Figure 1a. Percentage value over each set of bars highlight the change in the CONTROL relative to CALM simulation. Error bars in (c) and (d) show the standard deviation. Asterisks denote that CONTROL and CALM estimates are significantly difference within the 95% confidence interval (*) or within the 99% confidence interval (**) using a Student's test. BB = Baffin Bay; BG = Beaufort Gyre; BS = Barents Sea; CAA = Canadian Arctic Archipelago; HB = Hudson Bay; LS = Labrador Sea; NPP = net primary production. carbon export by 17-47% in spring (Figures 6c and 6d) and by 42–59% in autumn (Figures 6g and 6h). In perennially ice-covered regions, the impact of high-frequency winds on NPP and carbon export is much smaller for both seasons, and not statistically significant in the Beaufort Gyre (BG in Figures 6c, 6d, 6g, and 6h). A comparison of the seasonality of chl-a concentration in the CONTROL and CALM simulations allows us to further investigate how high-frequency winds influence the primary production in each region. In Figure 7 we show the evolution of 5-day climatologies of chl-a concentration in CONTROL (red area) and CALM (green area) along with environmental characteristics, namely, sea ice concentration (blue area and dashed line), shortwave radiation and PAR (gray-white background), upper layer turbulence (Ri solid red/green line), and standardized macronutrient and iron limitation (red/green circles ranging from minimum 0 to maximum of 1). The gray-white background above the sea ice concentration curve is the shortwave radia- tion from the CMC-GDPS forcing, and below the sea ice concentration curve is the PAR of the CONTROL simulation averaged over the top 50 m of the water column or the mixing layer depth, whichever is deeper. In ice-free and seasonally ice-covered regions (Figures 7a–7d), the presence of high-frequency winds drive a larger spring bloom and a second bloom in autumn in the CONTROL simulation (red). This contrasts CASTRO DE LA GUARDIA ET AL. 2737 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 7. Seasonal dynamics of simulated chlorophyll-a (chl-a) concentration and regional environment in the CONTROL (red) and the CALM (green) and in each region identified in Figure 1a: (a and b) ice-free regions, (c and d) seasonally ice-covered, and (e and f) perennially ice covered. Five-day climatologies of chl-a, surface water turbulence defined by the Richardson Number (Ri; water is turbulent (turb.) when Ri ≤ 0.25, and strongly stratified (strat.) when Ri is large), nutrient-limited growth rate standardized to the seasonal maximum and differentiating between macronutrient (DIM, open circles) and dissolved iron (DFed, filled circles) limitation (range: min = 0 to max = 1), and sea ice concentration (IC). The light regime above the sea ice is the shortwave radiation (showing a range from 0–150 W/m2 ), and below is the photosynthetically available radiation from CONTROL simulation as reference (showing range: 0–10 W/m2 ) since light limitation and photosynthetically available radiation are similar in both simulations. The Ri turbulent threshold of 0.25 is indicated with a dotted-black line. All proprieties, but IC, are averaged over the top 50 m. BB = Baffin Bay; BG = Beaufort Gyre; BS = Barents Sea; CAA = Canadian Arctic Archipelago; DIM = dissolved inorganic macronutrient; HB = Hudson Bay; LS = Labrador Sea. with the CALM simulation (green), in which phytoplankton has primarily a spring bloom. The onset of autumn turbulence (first day when Ri ≤ 0.25 after summer maximum) is particularly important in the development of the autumn bloom because it releases the phytoplankton from summer nutrient limitation (circles in Figures 7a–7d descend from its maximum limitation value). In the CONTROL simulation, upper waters become turbulent earlier in autumn (September) compared to the CALM simulation (late October). Since light is not limiting growth in early autumn, phytoplankton blooms as soon as high-frequency winds start eroding the summer stratification of the water column in the CONTROL simulation. In contrast, in the CALM simulation this erosion occurs later in time, when light (PAR) is already limiting phytoplankton growth (line with triangles in Figures 4a–4d). We note that the seasonality of temperature and light limita- tion is similar between the two simulations (Figure 4) and that the modeled sea ice concentration compares well with observations in the summer-autumn period (Figures A1a–A1d), making our sensitivity simulation robust and increasing our confidence on the link between high-frequency winds and autumn bloom. In perennially ice-covered regions (Figures 7e and 7f), phytoplankton is light-limited throughout most of the year (October–June; Figures 4e and 4f). This light limitation is likely exacerbated by a regional over esti- mation of sea ice concentration in the model during summer and autumn (Figures A1e and A1f). Still, in the CONTROL simulation (red), the high-frequency winds generate more turbulent mixing (i.e., lower Ri in CONTROL) resulting in relatively low nutrient limitation from June to October, relative to the CALM simulation. However, the high surface water stratification and low surface shear causes Ri values to be rela- tively large throughout the year and rarely below the turbulent threshold (Ri ≤ 0.25), specially in the CALM simulation. Overall, the high-frequency winds act to reduce the seasonal nutrient limitation in all the six regions, which results in the increase in NPP and carbon export in CONTROL relative to the CALM simulation. In particular, macronutrients (i.e., dissolved inorganic macronutrient, DIM) limits phytoplankton growth throughout summer in all regions except the Labrador Sea, where iron is important (circles, filled vs. unfilled in Figure 7). This summer iron limitation in Labrador Sea is exacerbated by high-frequency winds (more red than green filled circles in Figure 7b). Iron limitation in the Labrador Sea has received little attention, likely because of difficulties in sampling methods and a tendency of the research community to believe that the North Atlantic is not iron limited. We suggest that this could be a worthy topic of future research. CASTRO DE LA GUARDIA ET AL. 2738 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 8. Controls on the phytoplankton bloom in spring (March–June) and autumn (September–December) in CONTROL (filled circles) and CALM (filled squares) color coded for the six regions identified in Figure 1a. (a and b) Basin-averaged spring chl-a concentration in a parameter space defined by the turbulent period (number of weeks in the year with Richardson number Ri ≤ 0.25) and the spring (a) DIM concentration or (b) dissolved iron concentration from January–February. (c and d) Basin-averaged autumn chl-a concentration in a parameter space defined by the depth averaged PAR on the first day of turbulence (tRi ; first day of Ri ≤ 0.25 after summer maximum) and the autumn (c) DIM concentration or (d) DFed concentration. Negative PAR (PAR = −5 W/m2 ) indicates that the region did not become turbulent in autumn. The marker type indicates the experiment, the marker size is proportional to the chl-a concentration, and the marker color defines the region (see legend on the left). DIM = dissolved inorganic macronutrient; PAR = photosynthetically available radiation. 3.2.2. Turbulence and PAR as Drivers of Bloom Dynamics The spring bloom starts in response to the increasing light availability (i.e., declining light limitation on growth Figure 4), but the magnitude of the spring bloom is more sensitive to the length of the autumn-winter mixing period. The turbulent period is generally 10 weeks longer in the CONTROL simulation (filled circles in Figures 8a and 8b) compared to the CALM simulation (filled squares in Figures 8a and 8b). The longer turbulent period in CONTROL results in larger DIM and iron concentrations in the upper ocean waters (y axis in Figures 8a and 8b), which fuels a larger spring bloom relative to CALM (i.e., larger circles size relative to squares size in Figures 8a and 8b). The autumn bloom is instead conditioned by the timing of autumn turbulent mixing (i.e., the day of Ri ≤ 0.25 after summer maximum in Figure 7). The seasonally early onset of the autumn turbulence under high-frequency winds (filled circle in Figures 8c and 8d) results in nutrients becoming available to phyto- plankton when PAR is still high (i.e., PAR ≥ 1 W/m2 ). In contrast, in the absence of high-frequency winds (filled squares in Figures 8c and 8d), the seasonally late onset of autumn turbulence results in nutrients becoming available to phytoplankton when PAR is already limiting growth (PAR ≤ 1 W/m2 and see light limitation term (triangles) in Figure 4; PAR = −5 W/m2 indicates that the region did not experienced turbulence in autumn). In perennially ice-covered regions there is no autumn bloom because phytoplankton is light-limited in autumn in both simulations (i.e., green and brown symbols in Figure 8b are associated with PAR ≤ 1 W/m2 , and see also filled triangles in Figure 4). In the ice-free Labrador Sea, which is in the southern limits of the subarctic, PAR is available in late autumn (i.e., purple symbols in Figures 8c and 8d are associated with PAR ≥ 5 W/m2 ). This allows for a small autumn phytoplankton bloom even when turbulence starts in late October in the CALM simulation (Figures 7c and 7d). In general, we find that high-frequency winds CASTRO DE LA GUARDIA ET AL. 2739 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 9. Five-day climatologies (2002–2015) of net obduction (positive) and subduction (negative) of DIM at the base of the mixing layer (shaded area) for the CONTROL (red) and CALM (green) simulation. (a and b) Perennially ice-covered, (c and d) seasonally ice-covered, and (e and f) ice-free regions as defined in Figure 1a. All fluxes are expressed in phosphate units. Vertical line indicate the day the onset of turbulence in autumn with the absence of a line indicating that the region did not become turbulent in autumn. Note the different y axis range among basins. BB = Baffin Bay; BG = Beaufort Gyre; BS = Barents Sea; CAA = Canadian Arctic Archipelago; DIM = dissolved inorganic macronutrient; HB = Hudson Bay; LS = Labrador Sea. induce an earlier onset turbulence that allows for the increase in chl-a concentration in autumn, when light conditions can still favor phytoplankton growth. 3.2.3. Obduction of DIM Is Modulated by High-Frequency Forcing In previous sections we have illustrated that high-frequency winds impact productivity by influencing the timing when upper waters become turbulent. In this section, we unveil the physical mechanisms driving the obduction of DIM (see section 2). Obduction of DIM (i.e., nutrient enter the mixing layer) occurs from autumn to spring (positive values in Figure 9), while DIM subduction (i.e., nutrient exit from the mixing layer) occurs in spring and early sum- mer (negative values in Figure 9). Exceptions to this seasonal variability are observed in the seasonally ice-covered Hudson Bay, where obduction occurs mainly in summer, and perennial ice-covered regions, where obduction starts in summer and continues throughout winter. In general, the magnitude of the DIM obduction is 0.7 to 2.3 times larger when high-frequency winds are present (positive red shade) than when they are absent (positive green shade in Figure 9). Similarly, the mag- nitude of the DIM subduction is 1.2 to 3.0 times larger when high-frequency winds are present (negative red shade) than when they are absent (negative green shade). In particular, in ice-free and seasonally ice-covered regions (Figures 9a–9d), the presence of high-frequency winds (red shade) leads to an earlier onset of turbu- lence in autumn, which leads to an earlier onset of DIM obduction compared to when high-frequency winds are absent (green shade). In perennially ice-covered regions (Figures 9e and 9f), the magnitude of obduc- tion is significantly reduced compared to ice-free and seasonally ice-covered regions, but DIM obduction is less sensitive to high-frequency winds. The latter is in agreement with waters remaining nearly stratified all year-round (solid lines in Figures 7e and 7f). The processes dominating the flux of nutrients into and out of the mixing layer are different in each region but similar between simulations (Figure 10). In all regions, but the Beaufort Gyre which does not experience turbulence in autumn, the beginning of autumn turbulence (vertical black line in Figure 10) marks the beginning of DIM obduction by means of a deepening of the mixing layer (i.e., positive 𝛥MLD in Figure 10). In ice-free regions (Figures 10a–10d), DIM obduction and subduction is mostly driven by deepening and shallowing of the mixing layer depth, respectively. There is only a small DIM obduction by vertical advection (blue shade) from winter to summer. In the Barents Sea (Figures 10c and 10d), horizontal advection (black shade) drives a small additional DIM subduction in the winter-summer period. In seasonally ice-covered regions (Figures 10e–10h), DIM obduction over the autumn-winter period is driven by the deepening of the mixing layer, although this flux is largely reduced compared to ice-free CASTRO DE LA GUARDIA ET AL. 2740 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 10. Five-day regional climatologies (2002–2015) of the obduction of DIM at the base of the mixed layer by diffusion (gray shade), horizontal advection (H. Advection, black shade), vertical advection (W. Advection, blue shade) and entrainment by the seasonal fluctuation of the mixed layer depth (𝛥MLD, yellow shade). Obduction components are calculated as in the obduction equation (1) described in section 2.4.2. All fluxes are expressed in phosphate units. Positive flux denotes DIM obduction (flux enters the mixing layer), and negative flux denotes DIM subduction (flux leaves the mixing layer). (a–d) Ice-free regions, (e–h) seasonally ice-covered regions, and (i–l) perennially ice-covered regions as define in Figure 1a. BB = Baffin Bay; BG = Beaufort Gyre; BS = Barents Sea; CAA = Canadian Arctic Archipelago; DIM = dissolved inorganic macronutrient; HB = Hudson Bay; LS = Labrador Sea. CASTRO DE LA GUARDIA ET AL. 2741 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure 11. Impact of sea ice loss on bloom dynamics in Hudson Bay in the presence (a) and in the absence (b) of high-frequency winds. Comparing a climatology of years with high sea ice loss (red) with years with low sea ice loss (green). Five-day climatologies of chl-a, surface water turbulence defined by the Richardson Number (Ri; water is turbulent when Ri ≤ 0.25, and strongly stratified when Ri is large), nutrient-limited growth rate standardized to the seasonal maximum and differentiating between macronutrients (DIM, open circles) and dissolved iron (Fe, filled circles) limitation, and sea ice concentration (IC). The light regime above the sea surface is the shortwave radiation (showing a range from 0–150 W/m2 ), and within the water column, it is the photosynthetically available radiation (showing range: 0–10 W/m2 ). The Ri turbulent threshold of 0.25 is indicated with a dotted-black line. Photosynthetically available radiation, chl-a, Ri, and DIM Lim. are averaged over the top 50 m. regions. This reduction is compensated by a summer DIM obduction driven by diffusion (gray shade in Figures 10e–10h), especially in Hudson Bay. The large summer diffusion may be associated to freshwater input from sea ice melt with low DIM content. This effect would be amplified in our simulation since fresh- water fluxes have zero DIM concentration. Additionally, the fluxes may be larger in Hudson Bay relative to Baffin Bay, because Hudson Bay is a shallow basin with a maximum depth of ∼250 m, while Baffin Bay has a maximum depth greater than 2,000 m. The shallow depth in Hudson Bay can lead to a stronger DIM ver- tical gradient (and hence stronger DIM flux) because sinking particulate organic matter accumulates and remineralizes at the seafloor, which is relatively close to the mixing layer depth (40 m in winter and 15 m in summer). Finally, DIM subduction in the seasonally ice-covered regions considered is primarily driven by shallowing of the mixing layer in spring and with a small contribution by horizontal advection throughout the year. In perennially ice-covered regions (Figures 10i–10l), the shallowing of the mixing layer depth drives a strong DIM subduction in the spring-summer period and the deepening of the mixing layer drives a weak DIM obduction in winter. In the Canadian Arctic Archipelago (Figures 10i and 10j), horizontal and vertical advec- tion drives an additional DIM subduction during the winter-spring period. While in the Beaufort Gyre (Figures 10k and 10l), a summer-autumn diffusion drives DIM obduction which is stronger in the absence of high-frequency winds. Overall, the obduction analysis shows that high-frequency winds are most important driving a DIM obduc- tion in the autumn-winter period. The physical mechanism by which high-frequency winds enhance this DIM obduction is by contributing to deepening of the mixing layer depth. 3.3. Sea Ice Loss in the Absence of High-Frequency Winds Does Not Increase NPP Nor Biogenic Carbon Export Thus far we have shown that high-frequency winds play a major role in nutrient resupply, but the question remaining is whether the loss of sea ice can increase the phytoplankton production and biogenic carbon export by means of only increasing surface water PAR independently of high-frequency winds (Arrigo et al., 2008; Pabi et al., 2008). To test this hypothesis in our seasonally and perennially ice-covered regions, we built climatologies for high and low sea ice loss conditions (details in section 2.3). In the presence of high-frequency winds (showing only Hudson Bay in Figure 11a), less sea ice in Hudson Bay increases the phytoplankton productivity by ∼11%, but not in Baffin Bay, the Beaufort Sea, or the Cana- dian Arctic Archipelago (not shown). In the absence of high-frequency winds, less sea ice does not increase phytoplankton productivity nor carbon export in any of the four regions (Figure 11b). Although, iron is the CASTRO DE LA GUARDIA ET AL. 2742 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 limiting nutrient in late spring, the phytoplankton growth in Hudson Bay during this period is limited by temperature and not by nutrients (i.e., filled squares closest to yellow shading in Figure 4c). The response of phytoplankton to interannual mean sea ice losses (as presented here) may differ from long-term trends in sea ice losses due to global warming. The processes leading to the difference in the length of the ice-free period between our high and low sea ice concentration climatologies could result from inter- annual changes that are not directly related to increasing radiative forcing. The consideration of long-term sea ice loss is left for future research. 4. Discussion Our analysis suggests that high-frequency winds play a major role in present day (2002–2015) phytoplankton productivity and biogenic carbon export throughout the Arctic and sub-Arctic. High-frequency winds are shown to be crucial for generating the autumn bloom in seasonally ice-covered and ice-free regions, in particular, in Hudson Bay and the Labrador Sea. High-frequency winds are responsible for over 50% of the autumn primary production and carbon export and between 20 and 40% of the spring primary production and carbon export (Figure 6). The higher sensitivity of the phytoplankton to high-frequency winds in autumn was also found by Wu et al. (2007) in the Labrador Sea, where a simulated storm in autumn was shown to enhance the autumn bloom by up to 5 mg/m3 , while the same storm in spring only increased the bloom by a maximum of 1 mg/m3 . Wu et al. (2007) suggest that the magnitude of the phytoplankton response to the passage of a storm is directly related with the degree of water column stratification. While this may apply to ice-free regions like the Labrador Sea, we find that in the perennially ice-covered regions like the Beaufort Gyre the high-frequency winds has no effect on primary production even though the region is characterized by strong summer strat- ification. The latter is in agreement with a field survey in 2008, which reported that strong summer winds did not enhance phytoplankton primary production at an off-shelf station in the Beaufort Gyre because the wind-induced mixing energy was insufficient to break the strong summer stratification (Tremblay et al., 2011). Meanwhile, others found that in 1996, an intense stormy events with winds speed ≥16 m/s broke the summer stratification in the Beaufort Gyre and induced heat and salt exchanges into the mixed layer (Yang et al., 2004). Overall, this suggest that the impact of high-frequency winds depends on the local stratification and the intensity of the wind-induced mixing. The high-frequency winds trigger a larger autumn bloom by stimulating DIM obduction via deepening of the mixing layer in early September. This DIM obduction released phytoplankton from the summer nutri- ent limitation, and because light is still available in September, it favored the generation of the autumn bloom. In spring, the high-frequency winds support a larger bloom by stimulating a stronger and longer mixing period leading to higher nutrient accumulation in surface water between autumn and winter. Over- all, the high-frequency winds increased nutrient flux to surface waters by enhancing turbulence (Figure 9) with an annual maximum in the ice-free Labrador Sea of 2.3 times the value in the simulation without high-frequency winds. This value is similar to an observed estimate from the ice-free Pacific waters, where a transient storm induced a significant surface turbulence that increased the daily mean nutrient flux to surface waters by 3 times the background value (Rumyantseva et al., 2015). The simulation revealed that late spring-early autumn productivity in the Labrador Sea is limited by iron and not macronutrients as in the other basins (Figure 7) and that this iron limitation was exacerbated by high-frequency winds. In our simulation the iron limitation occurs due to the combination of a low input of iron to the surface waters in the Labrador Sea with high export-driven scavenging due to rapid macronutrient resupply. Iron limitation in the western North Atlantic has received little attention, but field measurements in summer 2010 captured surface iron concentrations near zero in Labrador Sea (Rijkenberg et al., 2014). Likewise, field measurements in the high-latitude eastern North Atlantic and the central Arctic indicate that iron indeed limits phytoplankton productivity in some regions (Birchill et al., 2019; Nielsdóttir et al., 2009; Rijkenberg et al., 2018). These results suggest that further exploring iron limitation in the Labrador Sea could be a worthy target of future research. We find that Hudson Bay is the only region in which less sea ice increases the annual mean primary pro- duction (by ∼11%), and only when high-frequency winds are present (Figure 11). This suggests that in the near future, when longer ice-free periods are predicted in Hudson Bay (Castro de la Guardia et al., 2013), a parallel increase in high-frequency winds may be necessary to maintain a high regional primary production. CASTRO DE LA GUARDIA ET AL. 2743 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Meanwhile in the Canadian Arctic Archipelago, Beaufort Gyre, and Baffin Bay, neither a decreasing sea ice cover (by ∼11–27%) nor an increasing ice-free period (by ∼14–16%) has a significantly effect on regional phytoplankton productivity or carbon export. These results contrast Arrigo et al. (2008) who estimated that a recession in sea ice cover of 20–25% increased Arctic Ocean NPP by ∼7%. However, our study focuses on a relatively short and recent period (2002–2015) where sea ice losses are small compared to those experienced since 1979, and it is possible that increased light could be more important with long-term sea ice losses. Also, our sea ice model overestimates the summer sea ice concentration in the central Arctic resulting in strong light limitation even in the climatology with less sea ice and, thus, reducing the effectiveness of the sea ice sensitivity experiment. We caution that the quantitative details of the response presented here may differ from other biogeochemi- cal models. For instance, BLING version 0 has a relatively shallow remineralization profile that can lead to a more pronounced vertical gradient of DIM in the upper water column. BLING also uses an implicit biomass tracer, which can lead to a very rapid phytoplankton response to shallowing mixed layer depths. These two features were shown by Galbraith et al. (2015) to produce overly intense high-latitude blooms. In addition, the absence of a zooplankton class results in a phytoplankton bloom that is controlled by bottom-up pro- cesses (i.e., light, temperatures, and nutrient conditions) as in the critical depth hypothesis (Sverdrup, 1953). Meanwhile, in some parts of the subarctic Atlantic (Behrenfeld, 2010) including the southern Labrador Sea (Marchese et al., 2019), the balance between grazing and optimal physical conditions control the initiation of the phytoplankton bloom in spring (i.e., dilution hypothesis of Behrenfeld, 2010). The absence of an explicit grazer class in BLING could thus result in our model having too large of a phytoplankton bloom that starts too earlier in the season. Nonetheless, BLING's mechanistic and qualitative response are robust. 5. Conclusion Our study, comparing a CONTROL simulation with high-frequency winds against a CALM simulation with- out high-frequency winds, supports Ardyna et al.'s (2014) attribution of more frequent autumn blooms to stormy wind events in the Arctic. By using a numerical approach, we are able to uncover the mechanisms behind the positive relationship between autumn blooms and stormy wind events. We also demonstrate that over the time period of our numerical experiments (2002–2015), sea ice loss alone is not associated with increased regional primary production. The principal mechanism by which high-frequency winds influence phytoplankton productivity relates to the stimulation of an early onset of autumn DIM obduction by deepening of the mixing layer, and a larger winter DIM obduction due to a longer mixing period. The earlier timing of the DIM obduction is essential for the generation of the autumn bloom, while larger winter DIM obduction increases the magnitude of the spring bloom. Our numerical experiments suggest that increasing phytoplankton primary production in the Arctic Ocean may not only be driven by sea ice loss but also by processes influencing surface ocean mixing and nutrient upwelling. In particular, we find that high-frequency winds play significant role increasing primary production by alleviating phytoplankton nutrient limitation. Appendix A: Model Evaluation Galbraith et al. (2015) performed an in-depth multimodel evaluation of BLING Version 0 coupled to the global GFDL-ESM2M ocean model. They showed that BLING estimates the large-scale biogeochemical cycling relatively well at the global scale despite its relatively small number of tracers. However, BLING tended to overestimate particulate export and productivity at high latitudes. The authors suggested that these differences were associated with the deep mixed layer depth at the time of the year when light became available. In BLING, the biomass is a diagnostic variable, which means it can adjust quickly to an improve- ment in growth conditions. Thus, when light conditions increase even if the mixed layer is still deep, the phytoplankton has access to the high nutrient concentrations of the deep ocean, growing into a large bloom, and leading to higher than observed export rates. Taken into consideration the caveats highlighted in the global simulation of BLING-GFDL-ESM2M, we focused on a regional evaluation of the first BLING-NEMO simulation on the ANHA4 configuration. We evaluated the performance of the CONTROL simulation over the six representative regions. We compare observed fields of sea ice concentration, chl-a, and DIM to model climatologies constructed using the CASTRO DE LA GUARDIA ET AL. 2744 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure A1. Five-day climatologies of sea ice concentration (IC) in each basin as per the CONTROL simulation (red shade) and the satellite-derived PMW data (solid black line). Dashed black line marks the 25% threshold used to define ice-free conditions. The bias was calculated as the difference between the model mean and the observed mean. For the geographical locations of the regions please refer to Figure 1a of the main manuscript. BB = Baffin Bay; BG = Beaufort Gyre; BS = Barents Sea; CAA = Canadian Arctic Archipelago; HB = Hudson Bay; LS = Labrador Sea; PMW = Passive Microwave. time period of the sensitivity experiment (January 2002 and December 2015). In addition, we qualitatively compare simulated NPP with available estimates in published literature. A.1. Evaluation of Simulated Sea Ice To evaluate the simulated sea ice concentration and the length of the ice-free period, we used Passive Microwave (PMW) data from 2002–2015 derived using the Bootstrap Algorithm Version 2 (Comiso, 2000 updated 2000 updated 2015). We define the ice-free period as the number of days with <25% sea ice con- centration. This threshold was chosen because PMW sea ice concentration data has an error of 15% to 20% throughout the Arctic, especially near the ice margin and in seasonally ice-covered regions (Cavalieri et al., 1991; Comiso & Kwok, 1996; Comiso et al., 1997, 2017). The model (red in Figure A1) captures the sea- sonality of the sea ice concentration presented in the PMW data (black solid line) in each region. The best agreement is in winter and spring (January to May), and lowest agreement is in late summer and autumn (August to November). In ice-free regions, Barents Sea and the Labrador Sea (Figures A1a and A1b), the simulated sea ice concen- tration is ≤25% year-round, while the PMW data have a winter sea ice concentrations up to 30%. As a result, the simulated ice-free period is ∼12 days longer than that of PMW data. These differences in sea ice condi- tions, however, have a minor effect on annual productivity because, at this high latitude, light is unavailable during winter. In the summer-autumn period, the simulated sea ice concentration is within 2% of PMW estimates. In seasonally ice-covered regions, Hudson Bay and Baffin Bay (Figures A1c and A1d), the simulated ice-free period is in good agreement with the PMW data (differences ≤5 days). The model slightly underestimates the winter-spring sea ice concentration in Baffin Bay (15 ± 1.1%) and in Hudson Bay (4.1 ± 2.0%). In perennially ice-covered regions (Figures A1e and A1f), the model simulates the sea ice concentration with relatively good accuracy during spring and winter (absolute differences of ≤1%). However, the model overestimates the sea ice concentration during summer and autumn by 31 ± 11% in the Beaufort Gyre and 12 ± 1.6% in the Canadian Arctic Archipelago. Given that sea ice caps the ocean surface, these higher sea ice concentration likely results in less light available to phytoplankton growth in the simulation. A.2. Evaluation of Simulated Chlorophyll-a (chl-a) We evaluate the simulated seasonality of the chl-a concentration from 2003 to 2015 using chl-a data derived from ocean color measurements from the Moderate Resolution Imaging Spectroradiometer sensors on board of the Aqua satellite (MODIS-A; Ocean Biology Processing Group, 2003). We also obtained a bottled-derived in situ chl-a concentration data set in the Beaufort Gyre from data collected as part of the Beaufort Sea Exploration Program. Simulation and in situ data are averaged over the top 50 m of the water column. The seasonality of chl-a concentration is captured by the model in all basins (Figure A2). In ice-free (Figures A2a and A2b) and seasonally ice-covered (Figures A2c and A2d) regions, there are two peaks in chl-a, one in spring and one in autumn (i.e., two phytoplankton blooms). In perennially ice-covered regions CASTRO DE LA GUARDIA ET AL. 2745 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure A2. Monthly climatology of simulated chl-a concentration (red area) evaluated against satellite data from MODIS-A (solid line with black triangles) and in situ measurements when available (filled black square ±1 standard deviation). Chl-a from each in situ data and simulation have been averaged over the top 50 m of the water column in each region identified in Figure 1a of the main manuscript. BB = Baffin Bay; BG = Beaufort Gyre; BS = Barents Sea; CAA = Canadian Arctic Archipelago; HB = Hudson Bay; LS = Labrador Sea; MODIS = Moderate Resolution Imaging Spectroradiometer. (Figures A2e and A2f), there is a flat bloom that has approximately the same chl-a magnitude throughout summer. In the ice-free Barents Sea (Figure A2a) the model underestimates the spring bloom by 1.7 mg/m3 and the autumn bloom by 0.7 mg/m3 . We suggest the discrepancies may arise from a combination of a contamina- tion by colored dissolved organic matter in satellite chl-a data and the lack of nutrients in riverine input by our model. The Barents Sea receives the outflow of two large rivers, the Northern Dvina and the Pechorae. The high color index found in these rivers year-round (Brekhovskikh et al., 2009) indicates the presence of large amounts of colored dissolved organic matter, which is known to contaminate chl-a satellite retrievals causing spuriously high values of satellite chl-a concentrations (Camara Lins et al., 2017; Goyens et al., 2013; Griffin et al., 2018). Meanwhile, river runoff also contains large amounts of nutrients (Holmes et al., 2012), which can sustain large phytoplankton growth in the coastal zone as seen in the satellite chl-a retrievals. However, our simulation does not contain nutrients in river waters, which could be resulting in an underes- timation of the regional chl-a concentration due to an unrealistic nutrient limitation. Indeed, our simulation underestimates the spring nutrient concentration in the Barents Sea (Figure A3a). In the Labrador Sea (Figure A2b), the simulated chl-a concentration has the observed seasonality with the magnitude of the spring and autumn bloom being similar to that of the satellite data. The simulated spring bloom, however, starts earlier in the season causing a small overestimation of 0.2 mg/m3 from February to April. In seasonally ice-covered regions, Hudson Bay and Baffin Bay (Figures A2c and A2d), the simulated spring bloom starts in March, rather than in April, but peaks in May like in the satellite data. The model captures the magnitude of the spring bloom with only a slight overestimation of 0.5 mg/m3 in Baffin Bay. In both regions, the chl-a concentration reaches a minimum in August but rises again in September, alluding to the presence of an autumn bloom. In Hudson Bay, the simulation estimates an autumn bloom ∼0.3 mg/m3 larger than in the satellite data. The overestimation of chl-a concentration in spring in Baffin Bay and in autumn in Hudson Bay could be associated with the underestimation of sea ice concentration (Figures A1c and A1d). The simulated longer ice-free season allows phytoplankton to be less light limited at the beginning and end of the growing season. However, satellite data limitations in spring and autumn due to ice cover and cloudy conditions may also influence the low observed satellite chl-a estimate. CASTRO DE LA GUARDIA ET AL. 2746 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Figure A3. Monthly climatology of simulated DIMs concentration (red area) compared to WOA13 (black triangles) and in situ bottle data (black squares ±1 standard deviation). On the right axis, gray triangles show the percentage of WOA13 pixels containing data each month. DIM concentration, from each data sources, has been averaged over the top 50 m of the water column. Regions are identified in Figure 1a of the main manuscript. BB = Baffin Bay; BG = Beaufort Gyre; BS = Barents Sea; CAA = Canadian Arctic Archipelago; DIM = dissolved inorganic macronutrient; HB = Hudson Bay; LS = Labrador Sea. In perennially ice-covered regions, the Beaufort Gyre and the Canadian Arctic Archipelago (Figures A2e and A2f), the simulated monthly chl-a concentration captures the observations relatively well, with biases ≤ −0.1 mg/m3 . The slight underestimation of chl-a in the model could be caused by severe light limitation due to the overestimation of the regional sea ice concentration in the model relative to observations (Figures A1e and A1f). A.3. Evaluation of Simulated DIM We obtained dissolved inorganic nitrate (NO3 ) and phosphate (PO4 ) concentrations from gridded WOA13 climatology. WOA13 data are limited in regions with scarce number of observations (Garcia et al., 2014a). Therefore, we chose to present WOA13 estimates only in regions and months when there is at least 5% area coverage by observations. Additionally, we built a climatology from available in situ bottled data of dissolved inorganic nitrate and phosphate within Baffin Bay and the Labrador Sea from 1958 to 2009 provided by the Bedford Institute of Oceanography, Department of Fisheries and Oceans Canada, and within the Beaufort Gyre region from 2003 to 2015 provided by the Beaufort Sea Exploration Program. The biogeochemical model BLING Version 0 uses a generic DIM in units of phosphate to represent the mean concentrations of the two major macronutrients in the ocean: nitrate (NO3 ) and phosphate (PO4 ). Following Galbraith et al. (2010), we compare the simulated PO4 to an observed DIM by averaging the observed PO4 and NO3 , weighted by its P:N Redfield ratio; that is, DIM = 1∕2(PO4 + NO3 ∕16). In all basins, the simulated DIM seasonality has a sinusoidal shape with a maximum between April and May and a minimum between August and September (Figure A3). DIM data are significantly scarce in the Beaufort Sea, the Canadian Arctic Archipelago, and Hudson Bay, with only three, two and one months of data, respectively, to evaluate the model in these regions. In ice-free regions, Labrador Sea and Barents Sea (Figures A3a and A3b), the observations and simulated DIM follow similar seasonality. In the Labrador Sea, where both WOA13 and in situ measurements are available, the model underestimates the DIM concentration compared to both data sets with mean biases of 0.21 and 0.11 mmol/m3 , respectively. In the Barents Sea, the model also underestimates by 0.16 mmol/m3 CASTRO DE LA GUARDIA ET AL. 2747 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Table A1 Mean or Range of Phytoplankton Primary Production (mg C·m−2 ·day−1 ) in Studied Regions Where Acknowledgments We thank Westgrid and Compute Observed Estimates Were Available Canada for computational resources, Region Season CONTROL Observed (reference) where all model experiments were performed and are archived (http: Labrador Sea autumn 250 200 (Mignot et al., 2018) www.computecanada.ca). We thank Labrador Sea spring 600 800 (Mignot et al., 2018) Dr. Nathan Grivault for helping with the setup of the simulations. The Labrador Sea and Baffin Bay autumn 300 298 (Harrison & Cota, 1991) authors will also like to thank the two Baffin Bay autumn 230 227 (Harrison et al., 1982) anonymous reviewers and the Editors for helping improve this manuscript. Hudson Bay autumn 250 200 (Lapoussière et al., 2013) This work was supported by Natural Beaufort Sea/Darnley Bay autumn 30 20–310 (Mundy et al., 2007) Science and Engineering Research Council (NSERC) of Canada through International Create grants relative to WOA13. DIM underestimations in both regions are likely related to the high influence of runoff (ArcTrain-432295), Climate Change and/or glacial melt water containing DIM (Bhatia et al., 2013; Holmes et al., 2012; Markussen et al., 2016). and Atmospheric Research Grant Meanwhile, these land DIM sources are absent in our simulation. (VITALS-RGPCC 433898), a Discovery Grant awarded to Paul G. Myers In the seasonally ice-covered region Hudson Bay (Figure A3c), data are limited to August. DIM differences (rgpin227438-09), and the European Research Council (ERC) funding to between the model and observation on this month are of ∼0.2 mmol/m3 , with the model underestimating Eric D. Galbraith under the European the observation. The model underestimation may also result from the lack of DIM in river runoff. Mean- Unions Horizon 2020 research and while, in Baffin Bay (Figure A3d), the model simulates DIM concentrations within the range of the in situ innovation program (Grant Agreement 682602). Canadian Meteorological measurements but underestimates by ∼0.4 mmol/m3 relative to WOA13. The difference between WOA13 Centers Global Deterministic and in situ data suggests that in situ observations are probably preferable for model evaluation in this Prediction System (CMC-GDPS) region. In September, in situ data have a peak in nutrient concentration that is not present in our simula- atmospheric forcing data (2002–2015) are freely available at the Environment tion (Figure A3e), but the lack of an increasing response in the observed chl-a concentration from August and Climate Change Canada website to September (Figure A2d) could suggest that the peak in DIM is specific to the sampling year. (https://weather.gc.ca/grib/grib2_glb_ 25km_e.html). The Coordinate In the perennially ice-covered Beaufort Gyre (Figure A2e), the few observations from July to September Ocean-ice Reference Experiment show a good agreement between simulated and in situ DIM. The similarity between the observations and (CORE) interannual atmospheric forcing data (1958–2001) are freely model suggests the overall effect of the denitrification that characterizes this region (Devol et al., 1997) is available online (https://data1.gfdl. partially captured with BLING's generic macronutrient. In the Canadian Arctic Archipelago (Figure A2f) noaa.gov/nomads/forms/core/ the only two months of data available suggests a large underestimation of simulated DIM concentration of COREv2/CIAF_v2.html). The CORE data sets are collaboratively supported ∼0.5 mmol/m3 . However, the simulation captures the observed decline between August and September. by the National Center for Atmospheric A.4. Evaluation of Simulated NPP Research (NCAR) and the Geophysical Fluid Dynamics Laboratory (GFDL) Observed NPP in the Labrador Sea and Baffin Bay is 298 mg C·m−2 ·day−1 (Harrison & Cota, 1991), which is under the umbrella of the Climate similar to the simulated estimates derived by averaging the same two regions (300 mg ·m−2 ·day−1 ; Table A1). Variability and Predictability In Baffin Bay alone, autumn estimates of primary production (227 mg C·m−2 ·day−1 ) in Harrison et al. (1982) (CLIVAR) Working Group on Ocean Model Development (WGOMD). We are also very similar to our modeled estimate (230 mg C·m−2 ·day−1 ). Likewise, estimates of primary pro- acknowledge the effort the Beaufort duction derived from bio-optical floats near the Labrador Sea in spring (∼800 mg C·m−2 ·day−1 ) and autumn Sea Exploration Program based at the (∼200 mg C·m−2 ·day−1 ; values from Figure 4 of Mignot et al., 2018) are similar to simulated values in spring Woods Hole Oceanographic Institution in collaboration with researchers from (∼600 mg C·m−2 ·day−1 ) and autumn (∼250 mg C·m−2 ·day−1 ) in the Labrador Sea (Table A1). In Hudson Fisheries and Oceans Canada at the Bay observed estimates of early fall mean primary production are ∼200 mg C·m−2 ·day−1 (derived from 12 Institute of Ocean Sciences for stations in Figure 4a in Lapoussière et al, 2013). This estimate is only 50 mg C·m−2 ·day−1 lower than the providing the in situ bottle data within the Beaufort Sea (http://www.whoi. simulated regional autumn mean (Table A1). In Darnley Bay, southwest of the Beaufort Gyre, field measure- edu/beaufortgyre). We thank Dr. ments of phytoplankton productivity under the sea ice, integrated over the upper 50 m, ranged from 20 to Kumiko Azetzu-Scott and her team at 310 mg C·m−2 ·day−1 (Mundy et al., 2007). In our CONTROL simulation, mean values in the Beaufort Gyre the Bedford Institute of Ocean Sciences (BIO), Fisheries and Oceans are at the low end of this observed range, 30 mg C·m−2 ·day−1 (Table A1). The generally lower simulated NPP Canada for providing the nitrate and could be a result of our climatology and regional average, while observations are derived from single to few phosphate bottle data within the measurements in time and space. Labrador Sea and Baffin Bay (these data are available by request; http:// www.bio.gc.ca/science/data-donnees/ base/run-courir-en.php). Objectively References analyzed macronutrient (nitrate and phosphate) and oxygen climatology for Ardyna, M., Babin, M., Gosselin, M., Devred, E., Rainville, L., & Tremblay, J. E. (2014). Recent Arctic Ocean sea ice loss triggers novel fall model evaluation and initial phytoplankton blooms. Geophysical Research Letters, 41, 6207–6212. https://doi.org/10.1002/2014GL061047 conditions are extracted from the Arrigo, K. R., van Dijken, G., & Pabi, S. (2008). Impact of a shrinking Arctic ice cover on marine primary production. Geophysical Research World Ocean Atlas (WOA13) version 2 Letters, 35, L19603. https://doi.org/10.1029/2008GL035028 data. These data are freely available Bamber, J., den Broeke, M., Ettema, J., Lenaerts, J., & Rignot, E. (2012). Recent large increases in freshwater fluxes from Greenland into through the National Centers for the North Atlantic. Geophysical Research Letters, 39, L19501. https://doi.org/10.1029/2012GL052552 CASTRO DE LA GUARDIA ET AL. 2748 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Environmental Information (NCEI Behrenfeld, M. J. (2010). Abandoning Sverdrup's critical depth hypothesis on phytoplankton blooms. Ecology, 91(4), 977–989. https://doi. https://www.nodc.noaa.gov/OC5/ org/10.1890/09-1207.1 Bhatia, M. P., Kujawinski, E. B., Das, SB., Breier, CF., Henderson, P. B., & Charette, M. A. (2013). Greenland meltwater as a significant and woa13/woa13data.html), which is potentially bioavailable source of iron to the ocean. Nature Geoscience, 6, 274–278. https://doi.org/10.1038/ngeo1746 funded by the National Oceanic and Birchill, A. J., Hartner, N. T., Kunde, K., Siemering, B., Daniels, C., González-Santana, D., & Lohan, MC. (2019). The eastern extent of sea- Atmospheric Administration (NOAA). sonal iron limitation in the high latitude North Atlantic Ocean. Scientific reports, 9(1), 1435. https://doi.org/10.1038/s41598-018-37436-3 Initial ocean temperature and salinity Brekhovskikh, V. F., Volkova, Z. V., & Lomova, D. V. (2009). Assessment of water quality in the rivers of northern European Russia by using are derived from the Polar Science oxygen demand indices. Russian Meteorology and Hydrology, 34(5), 321–330. https://doi.org/10.3103/S1068373909050082 Center Hydrographic Climatology Camara Lins, R., Martinez, J. M., Motta Marques, D. d., Cirilo, J. A., & Fragoso, C. R. (2017). Assessment of chlorophyll-a remote sensing (PHC), which can be accessed for free algorithms in a productive tropical estuarine-lagoon system. Remote Sensing, 9, 516. https://doi.org/10.3390/rs9060516 at the Research Data Archive (https:// Carranza, M. M., & Gille, S. T. (2015). Southern Ocean wind-driven entrainment enhances satellite chlorophyll-a through the summer. rda.ucar.edu/datasets/ds285.2/; Journal of Geophysical Research: Oceans, 120, 304–323. https://doi.org/10.1002/2014JC010203 https://doi.org/10.5065/V0EZ-AH27). Castro de la Guardia, L., Derocher, A. E., Myers, P. G., Terwisscha van Scheltinga, A. D., & Lunn, N. J. (2013). Future sea ice conditions in This Research Data Archive is Western Hudson Bay and consequences for polar bears in the 21st century. Global Change Biology, 19(9), 2675–2687. https://doi.org/10. managed by the Data Engineering and 1111/gcb.12272 Curation Section of the Computational Cavalieri, D. J., Crawford, J. P., Drinkwater, M. R., Eppler, D. T., Farmer, L. D., Jentz, R. R., & Wackerman, C. C. (1991). Aircraft active and Information Systems Laboratory and passive microwave validation of sea ice concentration from the Defense Meteorological Satellite Program special sensor microwave at the National Center for Atmospheric imager. Journal of Geophysical Research, 96(C12), 21,989–22,008. https://doi.org/10.1029/91JC02335 Research (NCAR), which is sponsored Comiso, J. C. (2000 updated 2015). Bootstrap sea ice concentration from Nimbus-7 SMMR and DMSP SSM/I, Version 2 [1979–2016]. Boulder, by the National Science Foundation CO: NASA DAAC at National Snow and Ice Data Centre. (NSF). Passive Microwave (PMW) sea Comiso, J. C., Cavalieri, D. J., Parkinson, C. L., & Gloersen, P. (1997). Passive microwave algorithms for sea ice concentration: A comparison ice concentration data (2002–2015) of two techniques. Remote Sensing of Environment, 60(3), 357–384. https://doi.org/10.1016/S0034-4257(96)00220-9 derived using the Bootstrap Algorithm Comiso, J. C., & Kwok, R. (1996). Surface and radiative characteristics of the summer Arctic sea ice cover from multisensor satellite Version 2 can be downloaded free from observations. Journal of Geophysical Research, 101(C12), 28,397–28,416. https://doi.org/10.1029/96JC02816 the National Snow and Ice Data Center Comiso, J. C., Meier, W. N., & Gersten, R. (2017). Variability and trends in the Arctic sea ice cover: Results from different techniques. (NSIDC; https://nsidc.org/data/ Journal of Geophysical Research: Oceans, 122, 6883–6900. https://doi.org/10.1002/2017JC012768 nsidc-0079/versions/2). Chlorophyll Comiso, J. C., Parkinson, C. L., Gersten, R., & Stock, L. (2008). Accelerated decline in the Arctic sea ice cover. Geophysical Research Letters, (chl-a) data derived from the Moderate 35, L01703. https://doi.org/10.1029/2007GL031972 Resolution Imaging Spectroradiometer Dai, A., Qian, T., Trenberth, K. T., & Milliman, J. D. (2009). Changes in continental freshwater discharge from 1948–2004. Journal of sensors on board of the Aqua satellite Climate, 22, 2773–2791. (MODIS-A) can be freely access from Day, J. J., Holland, M. M., & Hodges, K. I. (2017). Seasonal differences in the response of Arctic cyclones to climate change in CESM1. the NASA's OceanColor Web, which is Climate Dynamics, 50, 3885–3903. https://doi.org/10.1007/s00382-017-3767-x supported by the Ocean Biology Deser, C., Walsh, J. E., & Timlin, M. S. (2000). Arctic sea ice variability in the context of recent atmospheric circulation trends. Journal of Processing Group (OBPG) at NASA's Climate, 13(3), 617–633. https://doi.org/10.1175/1520-0442(2000)0130617:ASIVIT2.0.CO;2 Goddard Space Flight Center (https:// Devol, A. H., Codispoti, L. A., & Christensen, J. P. (1997). Summer and winter denitrification rates in western Arctic shelf sediments. oceandata.sci.gsfc.nasa.gov/ Continental Shelf Research, 17(9), 1029–1050. https://doi.org/10.1016/S0278-4343(97)00003-4 MODIS-Aqua/Mapped/Monthly/9km/ Fichefet, T., & Maqueda, M. A. M. (1997). Sensitivity of a global sea ice model to the treatment of ice thermodynamics and dynamics. chlor_a/). Simulations are saved using Journal of Geophysical Research, 102(C6), 12,609–12,646. https://doi.org/10.1029/97JC00480 specific coding system (http://knossos. Galbraith, E. D., Dunne, J. P., Gnanadesikan, A., Slater, R. D., Sarmiento, J. L., Dufour, C. O., & Marvasti, S. S. (2015). Complex functionality eas.ualberta.ca/anha/anhatable.php): with minimal computation: Promise and pitfalls of reduced-tracer ocean biogeochemistry models. Journal of Advances in Modeling Earth spinup simulation (ANHA4-ENG3), Systems, 7, 2012–2028. https://doi.org/10.1002/2015MS000463 CONTROL simulation Galbraith, E. D., Gnanadesikan, A., Dunne, J. P., & Hiscock, M. R. (2010). Regional impacts of iron-light colimitation in a global (ANHA4-ENG4), and CALM biogeochemical model. Biogeosciences, 7(3), 1043–1064. https://doi.org/10.5194/bg-7-1043-2010 simulation (ANHA4-ELC004). These Garcia, H. E., Locarnini, R. A., Boyer, T. P., Antonov, J. I., Baranova, O., Zweng, M., & Johnson, D. (2014). World Ocean Atlas 2013: Dissolved data are available by request (Paul G. oxygen, apparent oxgyen utilization, and oxygen saturation Edited by S. Levitus, & A. Mishonov Technical, vol. 3, pp. 27: NOAA Atlas Myers, pgmyers@ualberta.ca; http:// NESDIS. knossos.eas.ualberta.ca/anha/ Garcia, H. E., Locarnini, R. A., Boyer, T. P., Antonov, J. I., Baranova, O., Zweng, M., & Johnson, D. (2014a). World Ocean Atlas 2013: Dissolved anhatable.php). inorganic nutrients (phosphate, nitrate, silicate) Edited by S. Levitus, & A. Mishonov Technical (Vol. 4, pp. 25): NOAA Atlas NESDIS. Garcia-Quintana, Y., Courtois, P., Hu, X., Pennelly, C., Kieke, D., & Myers, P. G. (2019). Sensitivity of Labrador Sea Water formation to changes in model resolution, atmospheric forcing, and freshwater input. Journal of Geophysical Research: Oceans, 124, 2126–2152. https://doi.org/10.1029/2018JC014459 Ginoux, P., Chin, M., Tegen, I., Prospero, J. M., Holben, B., Dubovik, O., & Lin, S. J. (2001). Sources and distributions of dust aerosols simulated with the GOCART model. Journal of Geophysical Research, 106(D17), 20,255–20,273. https://doi.org/10.1029/2000JD000053 Goyens, C., Jamet, C., & Schroeder, T. (2013). Evaluation of four atmospheric correction algorithms for MODIS-Aqua images over contrasted coastal waters. Remote Sensing of Environment, 131, 63–75. https://doi.org/10.1016/j.rse.2012.12.006 Griffin, C., McClelland, J., Frey, K., Fiske, G., & Holmes, R. (2018). Quantifying CDOM and DOC in major Arctic rivers during ice-free conditions using Landsat TM and ETM+ data. Remote Sensing of Environment, 209, 395–409. https://doi.org/10.1016/j.rse.2018.02.060 Harrison, W. G., & Cota, G. F. (1991). Primary production in polar waters: Relation to nutrient availability. Polar Research, 10(1), 87–104. https://doi.org/10.3402/polar.v10i1.6730 Harrison, W. G., Platt, T., & Irwin, B. (1982). Primary production and nutrient assimilation by natural phytoplankton populations of the eastern Canadian Arctic. Canadian Journal of Fisheries and Aquatic Sciences, 39(2), 335–345. https://doi.org/10.1139/f82-046 Holdsworth, A. M., & Myers, P. G. (2015). The influence of high-frequency atmospheric forcing on the circulation and deep convection of the Labrador Sea. Journal of Climate, 28(12), 4980–4996. https://doi.org/10.1175/JCLI-D-14-00564.1 Holmes, R. M., McClelland, J. W., Peterson, B. J., Tank, S. E., Bulygina, E., Eglinton, T. I., & Zimov, S. A. (2012). Seasonal and annual fluxes of nutrients and organic matter from large rivers to the Arctic Ocean and surrounding seas. Estuaries and Coasts, 35(2), 369–382. https:// doi.org/10.1007/s12237-011-9386-6 Hu, X., Sun, J., Chan, T. O., & Myers, P. G. (2018). Thermodynamic and dynamic ice thickness contributions in the Canadian Arctic Archipelago in NEMO-LIM2 numerical simulations. The Cryosphere, 12(4), 1233–1247. https://doi.org/10.5194/tc-12-1233-2018 Huang, C., & Qiao, F. (2009). The relationship between sea surface temperature anomaly and wind energy input in the Pacific Ocean. Progress in Natural Science, 19(10), 1409–1412. https://doi.org/10.1016/j.pnsc.2009.03.004 Hunke, E. C., & Dukowicz, J. K. (1997). An elastic-viscous-plastic model for sea ice dynamics. Journal of Physical Oceanography, 27(9), 1849–1867. https://doi.org/10.1175/1520-0485(1997)0271849:AEVPMF2.0.CO;2 Jackett, D. R., & Mcdougall, T. J. (1995). Minimal adjustment of hydrographic profiles to achieve static stability. Journal of Atmospheric and Oceanic Technology, 12(2), 381–389. https://doi.org/10.1175/1520-0426(1995)0120381:MAOHPT2.0.CO;2 CASTRO DE LA GUARDIA ET AL. 2749 Journal of Geophysical Research: Biogeosciences 10.1029/2018JG004869 Karleskind, P., Lévy, M., & Memery, L. (2011). Subduction of carbon, nitrogen, and oxygen in the northeast Atlantic. Journal of Geophysical Research, 116, C02025. https://doi.org/10.1029/2010JC006446 Lapoussière, A., Michel, C., Gosselin, M., Poulin, M., Martin, J., & Tremblay, J. E. (2013). Primary production and sinking export during fall in the Hudson Bay system, Canada. Continental Shelf Research, 52, 62–72. https://doi.org/10.1016/j.csr.2012.10.013 Large, W. G., & Yeager, S. G. (2009). The global climatology of an interannually varying air–sea flux data set. Climate Dynamics, 33(2), 341–364. https://doi.org/10.1007/s00382-008-0441-3 Levy, M., Bopp, L., Karleskind, P., Resplandy, L., Ethe, C., & Pinsard, F. (2013). Physical pathways for carbon transfers between the surface mixed layer and the ocean interior. Global Biogeochemical Cycles, 27, 1001–1012. https://doi.org/10.1002/gbc.20092 Lin, I., Liu, W. T., Wu, C. C., Wong, G. T. F., Hu, C., Chen, Z., & Liu, K. K. (2003). New evidence for enhanced ocean primary production triggered by tropical cyclone. Geophysical Research Letters, 30(13), 1718. https://doi.org/10.1029/2003GL017141 Madec, G., & the NEMO team (2008). NEMO ocean engine. Institut Pierre-Simon Laplace (IPSL), France, No 27, ISSN No 1288-1619: Note du Pôle de modélisation. Marchese, C., Castro de la Guardia, L., Myers, P. G., & Bélanger, S. (2019). Regional differences and inter-annual variability in the timing of surface phytoplankton blooms in the Labrador Sea. Ecological Indicators, 96, 81–90. https://doi.org/10.1016/j.ecolind.2018.08.053 Markussen, T. N., Elberling, B., Winter, C., & Andersen, T. J. (2016). Flocculated meltwater particles control Arctic land-sea fluxes of labile iron. Scientific Reports, 6, 24033. https://doi.org/10.1038/srep24033 Masina, S., Storto, A., Ferry, N., Valdivieso, M., Haines, K., Balmaseda, M., & Parent, L. (2015). An ensemble of eddy-permitting global ocean reanalyses from the MyOcean project. Climate Dynamics, 49, 1–29. https://doi.org/10.1007/s00382 McLaughlin, F. A., & Carmack, E. C. (2010). Deepening of the nutricline and chlorophyll maximum in the Canada Basin interior, 2003–2009. Geophysical Research Letters, 37, L24602. https://doi.org/10.1029/2010GL045459 McTigue, N. D., Gardner, W. S., Dunton, K. H., & Hardison, A. K. (2016). Biotic and abiotic controls on co-occurring nitrogen cycling processes in shallow Arctic shelf sediments. Nature Communications, 7, 13145. https://doi.org/10.1038/ncomms13145 Mignot, A., Ferrari, R., & Claustre, H. (2018). Floats with bio-optical sensors reveal what processes trigger the North Atlantic bloom. Nature Communications, 9(1), 190. https://doi.org/10.1038/s41467-017-02143-6 Mundy, C. J., Ehn, J. K., Barber, D. G., & Michel, C. (2007). Influence of snow cover and algae on the spectral dependence of transmitted irradiance through Arctic landfast first-year sea ice. Journal of Geophysical Research, 112, C03007. https://doi.org/10.1029/2006JC003683 Nielsdóttir, M. C., Moore, C. M., Sanders, R., Hinz, D. J., & Achterberg, E. P. (2009). Iron limitation of the postbloom phytoplankton communities in the Iceland Basin. Global Biogeochemical Cycles, 23, GB3001. https://doi.org/10.1029/2008GB003410 Ocean Biology Processing Group (2003). MODIS Aqua Level 3 global daily mapped 9 km chlorophyll a.Ver. 6. DAAC CA USA. Dataset accessed [2017-02-20]. Pabi, S., van Dijken, G. L., & Arrigo, K. R. (2008). Primary production in the Arctic Ocean, 1998–2006. Journal of Geophysical Research, 113, C08005. https://doi.org/10.1029/2007JC004578 Pan, S., Shi, J., Gao, H., Guo, X., Yao, X., & Gong, X. (2017). Contributions of physical and biogeochemical processes to phytoplankton biomass enhancement in the surface and subsurface layers during the passage of Typhoon Damrey. Journal of Geophysical Research: Biogeosciences, 122, 212–229. https://doi.org/10.1002/2016JG003331 Rijkenberg, M. J. A., Middag, R., Laan, P., Gerringa, L. J. A., van Aken, H. M., Schoemann, V., & de Baar, H. J. W. (2014). The distribution of dissolved iron in the West Atlantic Ocean. PLOS ONE, 9(6), 1–14. https://doi.org/10.1371/journal.pone.0101323 Rijkenberg, M. J. A., Slagter, H. A., Rutgers van der Loeff, M., van Ooijen, J., & Gerringa, L. J. A. (2018). Dissolved Fe in the deep and upper Arctic Ocean with a focus on Fe limitation in the Nansen Basin. Frontiers in Marine Science, 5, 88. https://doi.org/10.3389/fmars.2018. 00088 Rumyantseva, A., Lucas, N., Rippeth, T., Martin, A., Painter, SC., Boyd, TJ., & Henson, S. (2015). Ocean nutrient pathways associated with the passage of a storm. Global Biogeochemical Cycles, 29, 1179–1189. https://doi.org/10.1002/2015GB005097 Serreze, M. C., Holland, M. M., & Stroeve, J. (2007). Perspectives on the Arctic's shrinking sea-ice cover. Science, 315(5818), 1533–1536. https://doi.org/10.1126/science.1139426 Simmonds, I., & Rudeva, I. (2012). The great Arctic cyclone of August 2012. Geophysical Research Letters, 39, L23709. https://doi.org/10. 1029/2012GL054259 Smith, G. C., Roy, F., Mann, P., Dupont, F., Brasnett, B., Lemieux, J. F., & Bélair, S. (2014). A new atmospheric dataset for forcing ice–ocean models: Evaluation of reforecasts using the Canadian global deterministic prediction system. Quarterly Journal of the Royal Meteorological Society, 140(680), 881–894. https://doi.org/10.1002/qj.2194 Steele, M., Morley, R., & Ermold, W. (2001). PHC: A global ocean hydrography with a high-quality Arctic Ocean. Journal of Climate, 14(9), 2079–2087. Sverdrup, H. U. (1953). On conditions for the vernal blooming of phytoplankton. ICES Journal of Marine Science, 18(3), 287–295. https:// doi.org/10.1093/icesjms/18.3.287 Tremblay, J. E., Bélanger, S., Barber, D. G., Asplin, M., Martin, J., Darnis, G., & Gosselin, M. (2011). Climate forcing multiplies biological productivity in the coastal Arctic Ocean. Geophysical Research Letters, 38, L18604. https://doi.org/10.1029/2011GL048825 Vancoppenolle, M., Fichefet, T., Goosse, H., Bouillon, S., Madec, G., & Maqueda, M. A. M. (2009). Simulating the mass balance and salinity of Arctic and Antarctic sea ice. 1. Model description and validation. Ocean Modelling, 27(1–2), 33–53. https://doi.org/10.1016/j.ocemod. 2008.10.005 Wang, J., Kim, H. M., & Chang, E. K. M. (2017). Changes in Northern Hemisphere winter storm tracks under the background of Arctic amplification. Journal of Climate, 30, 3705–3724. https://doi.org/10.1175/JCLI-D-16-0650.1 Waterhouse, A. F., MacKinnon, J. A., Nash, J. D., Alford, M. H., Kunze, E., Simmons, H. L., & Lee, C. M. (2014). Global patterns of diapycnal mixing from measurements of the turbulent dissipation Rate. Journal of Physical Oceanography, 44(7), 1854–1872. https://doi.org/10. 1175/JPO-D-13-0104.1 Wu, Y., Platt, T., Tang, C. C. L., & Sathyendranath, S. (2007). Short-term changes in chlorophyll distribution in response to a moving storm: A modelling study. Marine Ecology Progress Series, 335, 57–68. https://doi.org/doi:10.3354/meps335057 Yang, J., Comiso, J., Walsh, D., Krishfield, R., & Honjo, S. (2004). Storm-driven mixing and potential impact on the Arctic Ocean. Journal of Geophysical Research, 109, C04008. https://doi.org/10.1029/2001JC001248 Yin, K., Harrison, P. J., Goldblatt, R. H., & Beamish, R. J. (1996). Spring bloom in the central Strait of Georgia: Interactions of river discharge, winds and grazing. Marine Ecology Progress Series, 138(1/3), 255–263. https://doi.org/doi:103354/meps138255 CASTRO DE LA GUARDIA ET AL. 2750