Assessing the Role of High‐Frequency Winds and Sea Ice Loss on Arctic Phytoplankton Blooms in an Ice‐Ocean‐Biogeochemical Model

The long‐term trend of increasing phytoplankton net primary production (NPP) in the Arctic correlates with increasing light penetration due to sea ice loss. However, recent studies suggest that 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 ice‐ocean‐biogeochemical model configured for the North Atlantic and the Arctic. In the model, the presence of high‐frequency winds stimulates nutrient upwelling by producing an earlier and longer autumn‐winter mixing period with deeper mixing layer. The early onset of autumn mixing results in nutrients being brought‐up to near‐surface waters before the light becomes the dominant limiting factor, which leads to the autumn bloom. The enhanced mixing results in higher nutrient concentrations in spring and thus a large spring bloom. The model also shows significant iron limitation in the Labrador Sea, which is intensified by high‐frequency winds. The effect of sea ice loss on NPP was found to be regionally 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 limitation and that the isolated effect of sea ice loss on light plays a comparatively minor role.


Introduction
Accelerated loss of Arctic sea ice in summer (Comiso et al., 2008(Comiso et al., , 2017Serreze et al., 2007) can increase the Arctic Ocean's primary production by enhancing the amount of light available at the ocean surface Pabi et al., 2008). This sea ice loss, however, is also accompanied by increasing freshwater input from sea ice melt. While the former decreases light limitation on phytoplankton productivity, the later increases the nutrient limitation. The input of buoyant waters stratifies the water column, which inhibits the upward transport of nutrients, limiting phytoplankton growth in the well-lit surface ocean waters (Harrison 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 phytoplankton 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.

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

Numerical Model
We simulate the physics using a coupled ocean sea ice system based on the Nucleus for European Modelling 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), including both thermodynamic and dynamic components (Fichefet & Maqueda, 1997;Hunke & Dukowicz, 1997). 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 dependency. 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.

Model Configuration: ANHA4
We use a regional configuration of NEMO for the Arctic and the North Atlantic (ANHA4) with a nominal 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 boundaries, 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 vertical 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 m 2 /s for 10.1029/2018JG004869 the vertical eddy diffusivity (K z ), which is dependant on the Richardson number (Ri) and on sea ice concentration. That is, K z 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.

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. Biological 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., 2014aGarcia et al., , 2014, and for dissolved iron and organic phosphate fields we used climatologies 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 latter simulation is a global configuration at 1 • nominal resolution that ran for 100 years. The climatologies were built from the last 20 years.

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 Simulations (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.

Atmospheric Forcing Data
From January 1958 to December 2001 we forced the model with the reanalysis product from the Coordinated 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.

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 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 significantly 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 atmospheric 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 concentration 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 comparison 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.

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 simulations, 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.

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, V z ; that is, Ri ≡ N 2 ∕V 2 z . 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): The vertical shear is indicative of the capacity of the water to become turbulent and is computed as 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.

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, 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, K z is the vertical diffusivity coefficient, A R 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, K z , 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 (K z = 10 −4 m 2 /s); otherwise, K z = 10 −5 m 2 /s. In addition, the background K z is reduced as a linear function of sea ice concentration toward a minimum (K z = 10 −6 m 2 /s) when the grid cell has 100% ice cover. This definition of K z follows NEMO's turbulent kinetic energy mixing scheme for subgrid-scale vertical processes (Madec & the NEMO team, 2008). By using K z = 10 −5 m 2 /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 K z ∼ 10 −6 m 2 /s (Waterhouse et al., 2014).

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 radiation (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 temperature 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 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 mixing 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.

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 biogenic 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 spatially 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 Canadian 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.

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 (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 radiation 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 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 limitation 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 estimation 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 relatively 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.

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 phytoplankton when PAR is still high (i.e., PAR ≥ 1 W/m 2 ). 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/m 2 and see light limitation term (triangles) in Figure 4; PAR = −5 W/m 2 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/m 2 , 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/m 2 ). 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 induce an earlier onset turbulence that allows for the increase in chl-a concentration in autumn, when light conditions can still favor phytoplankton growth.

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 summer (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 magnitude 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 turbulence 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 obduction 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  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 freshwater 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 vertical 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 advection 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 obduction 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.

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 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 Canadian 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 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 interannual changes that are not directly related to increasing radiative forcing. The consideration of long-term sea ice loss is left for future research.

Discussion
Our analysis suggests that high-frequency winds play a major role in present day (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(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/m 3 , while the same storm in spring only increased the bloom by a maximum of 1 mg/m 3 . 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 stratification. 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 nutrient 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. Overall, 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 production (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.
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)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(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 biogeochemical 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 processes (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.

Conclusion
Our study, comparing a CONTROL simulation with high-frequency winds against a CALM simulation without 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)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(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. 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 improvement 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.

Appendix A: Model Evaluation
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 time period of the sensitivity experiment (January 2002 andDecember 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(Comiso, updated 2000(Comiso, updated 2015. We define the ice-free period as the number of days with <25% sea ice concentration. 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., 1997Comiso et al., , 2017. The model (red in Figure A1) captures the seasonality 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 concentration 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 conditions, 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 ( 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/m 3 and the autumn bloom by 0.7 mg/m 3 . We suggest the discrepancies may arise from a combination of a contamination 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 underestimation 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/m 3 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/m 3 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/m 3 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. 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  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/m 3 . 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 (NO 3 ) and phosphate (PO 4 ) 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 (NO 3 ) and phosphate (PO 4 ). Following Galbraith et al. (2010), we compare the simulated PO 4 to an observed DIM by averaging the observed PO 4 and NO 3 , weighted by its P:N Redfield ratio; that is, DIM = 1∕2(PO 4 + NO 3 ∕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/m 3 , respectively. In the Barents Sea, the model also underestimates by 0.16 mmol/m 3  (Mignot et al., 2018) Labrador Sea spring 600 800 (Mignot et al., 2018) Labrador Sea and Baffin Bay autumn 300 298 (Harrison & Cota, 1991) Baffin Bay autumn 230 227 (Harrison et al., 1982) Hudson Bay autumn 250 200 (Lapoussière et al., 2013) Beaufort Sea/Darnley Bay autumn 30 20-310 (Mundy et al., 2007) relative to WOA13. DIM underestimations in both regions are likely related to the high influence of runoff and/or glacial melt water containing DIM (Bhatia et al., 2013;Holmes et al., 2012;Markussen et al., 2016). Meanwhile, these land DIM sources are absent in our simulation.
In the seasonally ice-covered region Hudson Bay ( Figure A3c), data are limited to August. DIM differences between the model and observation on this month are of ∼0.2 mmol/m 3 , with the model underestimating the observation. The model underestimation may also result from the lack of DIM in river runoff. Meanwhile, in Baffin Bay ( Figure A3d), the model simulates DIM concentrations within the range of the in situ measurements but underestimates by ∼0.4 mmol/m 3 relative to WOA13. The difference between WOA13 and in situ data suggests that in situ observations are probably preferable for model evaluation in this region. In September, in situ data have a peak in nutrient concentration that is not present in our simulation ( Figure A3e), but the lack of an increasing response in the observed chl-a concentration from August to September ( Figure A2d) could suggest that the peak in DIM is specific to the sampling year.
In the perennially ice-covered Beaufort Gyre ( Figure A2e), the few observations from July to September show a good agreement between simulated and in situ DIM. The similarity between the observations and model suggests the overall effect of the denitrification that characterizes this region (Devol et al., 1997) is partially captured with BLING's generic macronutrient. In the Canadian Arctic Archipelago ( Figure A2f) the only two months of data available suggests a large underestimation of simulated DIM concentration of ∼0.5 mmol/m 3 . However, the simulation captures the observed decline between August and September.

A.4. Evaluation of Simulated NPP
Observed NPP in the Labrador Sea and Baffin Bay is 298 mg C·m −2 ·day −1 (Harrison & Cota, 1991), which is similar to the simulated estimates derived by averaging the same two regions (300 mg ·m −2 ·day −1 ; Table A1). In Baffin Bay alone, autumn estimates of primary production (227 mg C·m −2 ·day −1 ) in Harrison et al. (1982) are also very similar to our modeled estimate (230 mg C·m −2 ·day −1 ). Likewise, estimates of primary production derived from bio-optical floats near the Labrador Sea in spring (∼800 mg C·m −2 ·day −1 ) and autumn (∼200 mg C·m −2 ·day −1 ; values from Figure 4 of Mignot et al., 2018) are similar to simulated values in spring (∼600 mg C·m −2 ·day −1 ) and autumn (∼250 mg C·m −2 ·day −1 ) in the Labrador Sea (Table A1). In Hudson Bay observed estimates of early fall mean primary production are ∼200 mg C·m −2 ·day −1 (derived from 12 stations in Figure 4a in Lapoussière et al, 2013). This estimate is only 50 mg C·m −2 ·day −1 lower than the simulated regional autumn mean (Table A1). In Darnley Bay, southwest of the Beaufort Gyre, field measurements of phytoplankton productivity under the sea ice, integrated over the upper 50 m, ranged from 20 to 310 mg C·m −2 ·day −1 (Mundy et al., 2007). In our CONTROL simulation, mean values in the Beaufort Gyre are at the low end of this observed range, 30 mg C·m −2 ·day −1 (Table A1). The generally lower simulated NPP could be a result of our climatology and regional average, while observations are derived from single to few measurements in time and space. Centers Global Deterministic Prediction System (CMC-GDPS) atmospheric forcing data (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) are freely available at the Environment and Climate Change Canada website (https://weather.gc.ca/grib/grib2_glb_ 25km_e.html). The Coordinate Ocean-ice Reference Experiment (CORE) interannual atmospheric forcing data  are freely available online (https://data1.gfdl. noaa.gov/nomads/forms/core/ COREv2/CIAF_v2.html). The CORE data sets are collaboratively supported by the National Center for Atmospheric Research (NCAR) and the Geophysical Fluid Dynamics Laboratory (GFDL) under the umbrella of the Climate Variability and Predictability (CLIVAR) Working Group on Ocean Model Development (WGOMD). We acknowledge the effort the Beaufort Sea Exploration Program based at the Woods Hole Oceanographic Institution in collaboration with researchers from Fisheries and Oceans Canada at the Institute of Ocean Sciences for providing the in situ bottle data within the Beaufort Sea (http://www.whoi. edu/beaufortgyre). We thank Dr. Kumiko Azetzu-Scott and her team at the Bedford Institute of Ocean Sciences (BIO), Fisheries and Oceans Canada for providing the nitrate and phosphate bottle data within the 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 analyzed macronutrient (nitrate and phosphate) and oxygen climatology for model evaluation and initial conditions are extracted from the World Ocean Atlas (WOA13) version 2 data. These data are freely available through the National Centers for