A- A+
Alt. Display

# A model-based analysis of physical and biological controls on ice algal and pelagic primary production in Resolute Passage

## Abstract

A coupled 1-D sea ice-ocean physical-biogeochemical model was developed to investigate the processes governing ice algal and phytoplankton blooms in the seasonally ice-covered Arctic Ocean. The 1-D column is representative of one grid cell in 3-D model applications and provides a tool for parameterization development. The model was applied to Resolute Passage in the Canadian Arctic Archipelago and assessed with observations from a field campaign during spring of 2010. The factors considered to limit the growth of simulated ice algae and phytoplankton were light, nutrients, and in the case of ice algae, ice melt. In addition to the standard simulation, several model experiments were conducted to determine the sensitivity of the simulated ice algal bloom to parameterizations of light, mortality, and pre-bloom biomass. Model results indicated that: (1) ice algae limit subsequent pelagic productivity in the upper 10 m by depleting nutrients to limiting levels; (2) light availability and pre-bloom biomass determine the onset timing of the ice algal bloom; (3) the maximum biomass is relatively insensitive to the pre-bloom biomass, but is limited by nutrient availability; (4) a combination of linear and quadratic parameterizations of mortality rate is required to adequately simulate the evolution of the ice algal bloom; and (5) a sinking rate for large detritus greater than a threshold of ∼10 m d–1 effectively strips the surface waters of the limiting nutrient (silicate) after the ice algal bloom, supporting the development of a deep chlorophyll maximum.

##### Knowledge Domain: Ocean Science
How to Cite: Mortenson, E., Hayashida, H., Steiner, N., Monahan, A., Blais, M., Gale, M.A., Galindo, V., Gosselin, M., Hu, X., Lavoie, D. and Mundy, C.J., 2017. A model-based analysis of physical and biological controls on ice algal and pelagic primary production in Resolute Passage. Elem Sci Anth, 5, p.39. DOI: http://doi.org/10.1525/elementa.229
Published on 24 Jul 2017
Accepted on 07 Jun 2017            Submitted on 09 Feb 2017
Domain Editor-in-Chief: Jody W. Deming; University of Washington, USA
Associate Editor: Jean-Éric Tremblay; Université Laval, Canada

## 1 Introduction

Satellite records indicate that the minimum annual sea ice extent in the Arctic has been decreasing by more than 10% per decade over the last half century (Vaughan et al., 2013), which results in a longer and more widespread open-water season (Barber et al., 2015). In addition to the loss of sea ice, there has been a general shift in ice type, from thicker multiyear ice to younger and thinner first-year ice (Lindsay and Schweiger, 2015). These trends in ice type, cover, and timing have significant consequences on marine and sea-ice ecosystems and air-sea exchange, as well as broader implications for the global climate. To reproduce recent changes and project future changes of sea ice related primary production in models we need to be able to understand the driving processes and develop adequate model parameterisations. 1-D models are excellent tools to develop such parameterisations and test the system sensitivity to parameter variations.

In the Arctic, ice algae live in a relatively sheltered environment concentrated within the bottom few centimeters of the sea ice (Smith et al., 1990; Galindo et al., 2014; Brown et al., 2015–1). Ice algal blooms occur at high polar latitudes where snow and ice-cover substantially reduce incident light to the bottom of the ice column. This environment, and the timing of ice algal blooms, suggest that they are shade-acclimated to low-light conditions (Kirst and Wiencke, 1995). The algae within the ice can reach very high biomass (exceeding 1000 mg Chl a m–3) that is up to two orders of magnitude greater than the underlying phytoplankton biomass (Galindo et al., 2014; Leu et al., 2015). Previous observational studies indicate that primary production by ice algae can make a substantial contribution to the total (sea ice and pelagic) primary production at various locations in the Arctic Ocean (Legendre et al., 1992; Gosselin et al., 1997). Ice algae are dependent on the ice as a habitat and also affect the ice through light absorption and its subsequent conversion to heat, and through production of extracellular polymeric substances (Riedel et al., 2006; Krembs et al., 2011). In addition, the termination of the ice algal bloom translates to nutrient release to, and possible seeding of, the phytoplankton bloom (Galindo et al., 2014) in the surface ocean.

One challenge for model studies of Arctic sea ice is that observations from the field are sparse due to the remote location and harsh environment. As a result, many parameters required to simulate biogeochemical processes in ice-covered regions are poorly constrained. In this modeling study, we have been able to take advantage of observations of ice algal blooms and environmental variables from several recent field campaigns at one location in order to better understand the processes constraining the simulation. To address the impact of remaining uncertainties, the modelled ice algal growth can be tested against variations in relevant parameters, with ranges based on measured or inferred uncertainty. Sensitivity analyses are a common way to assess the impact of specific processes or parameters on the whole system and evaluate the variables to which the system is most sensitive. Testing the model’s sensitivity over a certain parameter range, based on observations, allows for an estimate of the importance of a given process, compared to others, and identification of parameters that need to receive focused observational attention to reduce the overall uncertainty of the system (Steiner et al., 2016). Several 1-D sea ice algal models have been developed in order to reproduce observations at particular locations (Lavoie et al., 2005; Pogson et al., 2011). Some include focused sensitivity studies, e.g., Arrigo and Sullivan (1994), show that adjustments lowering the ice algal nutrient supply (via a nutrient transport coefficient) can cause the ice algal system to become nutrient-limited, and identify a high sensitivity to the ice algal growth rate. Jin et al. (2006) identified a strong correlation between net primary production of ice algae and the initial nutrient concentration in the water column. Steiner et al. (2016) highlighted several components and parameters that lack either full understanding or observational constraints. Based on these previous studies, the following parameters were selected for testing in this study: the amount of algae in the ice during the winter (pre-bloom biomass), photosynthetic efficiency of the ice algae in low light conditions, the strength of nutrient flushing during the ice algal bloom period, and the magnitude and form of specific mortality of the ice algae. While model studies suggest that ice algal seeding of an ice-associated pelagic bloom mainly affects the timing rather than the magnitude of the pelagic bloom (Jin et al., 2007; Tedesco et al., 2012) the link between ice algal and pelagic production remains an area of uncertainty and that we also address here.

Another challenge for both 1-D and 3-D modelling of sea ice ecosystems is the treatment of (subgrid-scale) heterogeneous snow cover and how this heterogeneity affects the light penetration to the bottom of the ice (where Arctic ice algae are most prominent). In order to represent a grid cell average over multiple square kilometers, this heterogeneity needs to be taken into account in the model. This challenge has been the focus of Abraham et al. (2015). They compare light penetration through a Rayleigh-distributed snow cover to a uniformly distributed snow cover, identifying substantial improvement to the grid-cell mean light transmission compared to observations. Light transmission to the bottom of the sea ice has been identified as a major problem in simulating ice algal growth particularly during the period of snow decline (Arrigo Sullivan 1994; Lavoie et al., 2005; Pogson et al., 2011). In the present study, we test the impact of the newly-developed parameterization for light transmission through sea ice (Abraham et al., 2015) on ice algal growth.

With the broader objective of establishing a set of parameterizations that can be transferred into a 3-D regional Arctic model (coupling sea-ice and the ocean along with associated ecosystems), this study uses a 1-D coupled sea ice-ocean physical-biogeochemical model to analyze the physical and biological controls on simulated ice algae and phytoplankton blooms. The analysis contains three distinct components: 1) Investigation of the impacts of subgrid-scale non-uniform snow depth distributions on the growth of ice algae by applying a new parameterization for light transmission through sea ice (Abraham et al., 2015); 2) assessment of the influences of ice algae on the simulated phytoplankton bloom by coupling and decoupling the sympagic and pelagic ecosystems; and 3) evaluating the sensitivity of the simulated ice algal bloom to a set of selected parameters and parameterizations following recommendations by Steiner et al. (2016). The test location for our model study is set in Resolute Passage in the Canadian Arctic Archipelago, based on the availability of a comparatively rich observational dataset at this location.

## 2 Methods

### 2.1 Model description

#### 2.1.1 Physical model

The sea ice component of the coupled sea ice-ocean physical model is the 1-D thermodynamic model of Flato and Brown (1996) with most recent updates from Abraham et al. (2015). These updates include new parameterizations for the light fields and heat fluxes through sea ice by accounting for a subgrid-scale snow depth distribution, melt ponds, and temperature-dependent extinction and transmissivity coefficients (see Appendix A1 for a synopsis of these updates). These new parameterizations improved the evolution of the simulated light fields under the landfast ice in Resolute Passage during the melt period of 2002 (Abraham et al., 2015). In the present study, some of the optical parameters of the sea ice model were modified to improve the fit of the simulated results to observations. A set of retuned optical parameters is provided with references for justification in Table 1. Although seasonal changes to the properties of snowfall have not been included in the present study, the snowfall rate has been varied with time based on specified precipitation data, in contrast to a prescribed constant rate as in earlier studies (Flato and Brown, 1996; Abraham et al., 2015).

Table 1

Extinction and transmissivity coefficients, as well as surface albedos used in this study. DOI: https://doi.org/10.1525/elementa.229.t1

Symbol Quantity Value Reference

κs,f Extinction coefficient for freezing snow 14 m–1 Grenfell and Maykut (1977)
κs,m Extinction coefficient for melting snow 7.5 m–1 Grenfell and Maykut (1977)
κi,m Extinction coefficient for melting sea ice 0.8 m–1 Light et al. (2008)
κm Extinction coefficient for melt ponds 0.5 m–1 Abraham et al. (2015)
κia Extinction coefficient for ice algae 0.017 (mmol N m–3)–1 m–1 McDonald et al. (2015)
κpd Extinction coefficient for phytoplankton and detritus 0.03 (mmol N)–3)–1 m–1 Lavoie et al. (2009)
i0,s,m Transmissivity coefficient for melting snow 0.15 Vancoppenolle et al. (2010)
i0,i,f Transmissivity coefficient for freezing sea ice 0.5 Lavoie et al. (2005)
i0,i,m Transmissivity coefficient for melting sea ice 0.5 Lavoie et al. (2005)
i0,m Transmissivity coefficient for melt ponds 0.5 Abraham et al. (2015)
αs,f Surface albedo of freezing snow 0.8 Vancoppenolle et al. (2010)
αs,m Surface albedo of melting snow 0.7 Lavoie et al. (2005)
αi,f Surface albedo of freezing sea ice 0.6 Within the range between Vancoppenolle et al. (2010) and Perovich et al. (2002)
αi,m Surface albedo of melting sea ice 0.5 Vancoppenolle et al. (2010)
αm Surface albedo of melt ponds 0.3 Light et al. (2008)

The physical processes in the water column are simulated by the General Ocean Turbulence Model (GOTM) of Burchard et al. (2006). GOTM provides the physical quantities required for computation of biogeochemical variables in the water column, such as horizontal velocity fields, turbulent transports, photosynthetically active radiation (PAR), and temperature. Details of model parameterizations for these quantities are provided in the GOTM website (http://www.gotm.net).

#### 2.1.2 Biogeochemical model

A biogeochemical model representing the lower-trophic level of sea ice and pelagic ecosystems in the Arctic was developed within the Framework for Aquatic Biogeochemical Models (Bruggeman and Bolding 2014) to facilitate the coupling with the physical model described above. The schematic diagram of the biogeochemical model is shown in Figure 1. The sea ice component of the biogeochemical model simulates the temporal evolution of four state variables (ice algae, nitrate, ammonium, and silicate) in the sea ice skeletal layer. The ice algae module is based on Lavoie et al. (2005). It was updated in this study by incorporating nitrate to account for potential algal growth reduction due to nitrogen limitation, as well as including ammonium to represent the biogeochemical processes within sea ice more realistically. At any given time, the growth of simulated ice algae is limited by one of the four limiting factors: light, ice melt, silicate, or nitrate. A limitation index for each factor is determined as a non-dimensional index that varies between 0 and 1 as in Lavoie et al. (2005). The ice algal growth rate is then determined by the minimum of the four indices multiplied by the specific growth rate at a given temperature of the ice skeletal layer (A2).

Figure 1

Schematic diagram of the coupled sea ice-ocean biogeochemical model. Circles represent the model state variables: nitrate (NO3), ammonium (NH4), silicate (Si), ice algae (IA), small phytoplankton (P1), large phytoplankton (P2), microzooplankton (Z1), mesozooplankton (Z2), small detritus (D1), large detritus (D2), and biogenic silica (BSi). Sinking variables are bounded by yellow circles. Black and red arrows represent paths of nitrogen and silicon transfers between the variables, respectively: photosynthesis (PH), nitrification (NI), diffusive mixing (DI), flushing (FL), seeding (SE), linear mortality (LM), quadratic mortatlity (QM), remineralization (RE), grazing (GR), ingestion (IN), sloppy feeding (SL, for inefficient grazing that leaves unconsumed but dead prey), and excretion (EX). DOI: https://doi.org/10.1525/elementa.229.f1

To study the sympagic-pelagic ecological interactions at the lower-trophic level, the sea ice biogeochemical model was coupled to a ten-compartment (small and large phytoplankton, microzooplankton, mesozooplankton, small and large detritus, biogenic silica, nitrate, ammonium, and silicate) pelagic biogeochemical module based on Steiner et al. (2006). This module was updated by including mesozooplankton as a prognostic variable and by partitioning detritus into small and large size classes. At the ice-water interface dissolved nutrients are exchanged through molecular diffusion. Ice algae released into the water column are treated similarly as in the coupled model of Lavoie et al. (2009): sloughed ice algae enter either the large phytoplankton pool in which they continue to grow or the large detritus pool in which they sink rapidly as aggregate products. The equations and parameters for the coupled biogeochemical model are provided in Appendix A2.

#### 2.1.3 Experimental design

The 1-D model was applied to simulate ice algae and pelagic primary production within and under the landfast first-year sea ice in Resolute Passage, at a location with a water depth of 141 m. Resolute Passage was chosen for the study site because extensive field research has been conducted in the area (Cota et al., 1987; Lavoie et al., 2005; Papakyriakou and Miller, 2011; Galindo et al., 2014; Brown et al., 2015–1; Geilfus et al., 2015). Specifically, model simulations were conducted for a location representative of the Arctic Ice Covered Ecosystem (Arctic-ICE) field campaign (74.71°N, 95.25°W). This field campaign took place during the spring of 2010 in order to study the physical and biological processes controlling the timing of ice algae and under-ice phytoplankton blooms (Mundy et al., 2014). The model was divided into 10 uniformly-spaced layers for sea ice and 100 layers for the upper 100 m of the water column. With the ultimate goal of implementing the parameterizations considered into coarser-resolution regional or global ocean circulation models, we do not attempt to resolve small-scale under-ice processes finer than 1 m. In order to limit the ultimate computational burden, we compared the 10-layer model to 5- and 2-layer simulations, deciding that the minor differences (1–2%) in output did not justify curtailing the effort at this stage.

The model was integrated for 8 months (1 February – 30 September, 2010) with a timestep of 10 minutes, and forced with Environment Canada’s hourly weather data (including surface air temperature, zonal and meridional wind at 10 m above the sea surface, surface air pressure, relative humidity, cloud cover, and precipitation) collected at the Resolute airport, located within 10 km of the study site. Temperature, salinity, and horizontal velocity fields of the ocean were restored over the entire water column with restoring timescale of 1 day (temperature and salinity) and 10 minutes (horizontal velocity) to the output of a 3-D regional model simulation (NEMO-LIM2) used in Dukhovskoy et al. (2016). We chose to restore the model this often in order to tightly constrain the physical water column properties and thus focus on comparing biogeochemical components of the model. The initial snow and melt pond depths and ice thickness were set to 5, 0, and 55 cm, respectively. The initial concentration of ice algae was set to 1.0 mmol N m–3 (ca. 3.5 mg Chl a m–3; the observed range of C:N:Chl a ratios is described in Appendix A2). The initial concentration of nitrate (silicate) was set to a constant value of 7.2 mmol N m–3 (14.7 mmol Si m–3) throughout the bottom ice and the water column, based on the measurements of these nutrients during the Arctic-ICE 2010 field campaign (Mundy et al., 2014; Galindo et al., 2014). The initial concentrations of ammonium both in the sea ice and the water column were assumed to be small (e.g., Harrison et al., 1990), and hence, set to 0.01 mmol N m–3. Similarly, the initial concentrations of all other pelagic biogeochemical state variables were set to 0.01 mmol N m–3 (mmol Si m–3 for biogenic silica) throughout the water column.

### 2.2 Observations

Observational data used to evaluate the model results include snow and melt pond depths, ice thickness, under-ice PAR, and chlorophyll a (Chl a). Measurements of these variables were conducted during May and June of 2010 as part of the Arctic-ICE field campaign. Observed snow and melt pond depths, ice thickness, and Chl a in the bottom 3 cm of sea ice were sampled at various sites of high, medium, and low snow covers. The mean value of Chl a is therefore an estimate of the site average, as presented in Galindo et al. (2014), and is comparable to a grid cell average in a regional or global model. Concentrations of Chl a in the water column were determined by collecting samples at five depths (2, 5, 10, 25, and 50 m below the sea surface) using 5 L Niskin bottles and following the procedures outlined in Galindo et al. (2014). In situ time series data for daily-mean under-ice (2 m below sea surface) PAR were collected using two independent tethers moored to the sea ice below high (>40 cm prior to snowmelt onset) and low (<20 cm prior to snowmelt onset) snow cover sites (within 4 – 6 m of the CTD casts). Technical details of these PAR measurements are provided in Mundy et al. (2014). In addition to the tether measurements, instantaneous under-ice PAR was estimated by extrapolating the 20 m depth CTD-based PAR measurement to the surface following Frey et al. (2011). Casts of CTD and a biospherical 4 pi sensor were obtained daily through the main sampling hole within a heated tent on the sea ice. Details of the CTD-based under-ice PAR estimates are described in Gale (2014).

## 3 Results

Results are divided into three parts based on the types of model simulations conducted. The first subsection evaluates the performance of the standard run. The second subsection compares the result of the standard run with a simulation that excludes ice algae. The third subsection provides the results of parameter sensitivity experiments. Specific setups of these runs are described in each of these subsections.

### 3.1 Model evaluation

The standard run was conducted with the setup outlined in the previous section (Experimental design) and by applying the Rayleigh distribution for representing the subgrid-scale snow depth variability (see Appendix A1). Abraham et al. (2015) indicated a better fit for the Rayleigh distribution than gamma probability distribution based on observations from the Arctic-ICE 2010 study (not shown).

#### 3.1.1 Snow and melt pond depths and ice thickness

In many previous 1-D model studies, the temporal evolution of snow depth was either prescribed to observed snow depth data (e.g., Lavoie et al., 2005; Pogson et al., 2011; Palmer et al., 2014) or simulated by prescribing a constant snowfall rate (Flato and Brown, 1996; Abraham et al., 2015). In this study, snow depth was simulated by prescribing a variable snowfall rate based on observed precipitation data. The simulated and observed time series of snow and melt pond depths are shown in Figure 2a. The simulated snow depth increased occasionally as a result of snowfall events until the maximum depth (ca. 20 cm) was reached by mid-May. In the standard run, the simulated snow started melting toward the end of May and completely vanished within 3 weeks. Snowmelt resulted in the formation of melt ponds which reached a maximum depth of 5 cm shortly after the snow disappearance. In comparison with the field measurements presented in Figure 2a, the timing of melt events was simulated reasonably with the distributed snow case.

Figure 2

Simulated and observed snow depth, melt-pond depth, and ice thickness. Time series of (a) simulated daily-mean snow (solid line) and melt pond (dashed line) depths, observed snow/melt pond depth (circles), and (b) simulated daily-mean (line) and observed (circles) ice thickness. Circles represent the site-average values with one standard deviations indicated by vertical bars. DOI: https://doi.org/10.1525/elementa.229.f2

Figure 2b shows the simulated and observed time series of ice thickness. In the standard run, simulated ice grew gradually to a maximum thickness of about 150 cm by early June and then started melting following the initial snowmelt. In the standard case, the distributed snow parameterization represents snow-free areas, which allows the ice to start melting before all the snow has disappeared. The simulated ice vanished completely in early July after which the sea surface remained ice-free until late September. The simulated ice thickness agreed well with the observations throughout the sampling period (Figure 2b), whereas the ice break up in the simulation occurred a week earlier than in the observations (Galindo et al., 2014). This difference could be attributed to dynamic processes of sea ice (e.g., wind-driven ridging and rafting) which are not accounted for in our 1-D model.

#### 3.1.2 Surface area fractions and under-ice PAR

Simulation of the light penetration through snow and sea ice is crucial for simulating a reasonable ice algal bloom, as the initial phase of the bloom is typically limited by light (Gosselin et al., 1990; Lavoie et al., 2005; Leu et al., 2015). During the melt period, surface area fractions of simulated snow, melt ponds, and bare ice undergo changes that affect the amount of light reaching the ice base as indicated in Figure 3. In the standard simulation, the surface of the simulated ice was fully snow-covered prior to the snowmelt onset. Consequently, the simulated daily-mean under-ice PAR during this period was less than 1 µmol photons m–2 s–1. This value is lower than either of the tether measurements, but in good agreement with most of the CTD-based estimates. In the model, about 10% of the snow surface was replaced with melt ponds due to snowmelt during the first week of June, resulting in an increase of the daily-mean under-ice PAR to more than 1 µmol photons m–2 s–1. This value is comparable to the tether measurements at high snow cover station, as well as to the CTD-based estimates. By June 9, the surface area coverage of simulated melt ponds extended to 30% (the maximum value prescribed by the model). Further areal loss of simulated snow resulted in an emergence of bare ice, which covered 70% of the ice surface following the snow disappearance. The pulsed effect in melt pond area in mid-June (Figure 3a) reflects daily signals associated with daytime melting and overnight freezing (causing surface bare ice). The simulated under-ice PAR during this period exceeded 10 µmol photons m–2 s–1 (Figure 3b), which is comparable to both the tether and the CTD-based observations. As expected, the simulated gridbox-mean under-ice PAR was quantitatively closer to the CTD-based (site-average) estimates than the tether (point) measurements. Furthermore, the standard simulation successfully reproduced the smooth seasonal transition of under-ice PAR that is evident in the tether measurements during the melt period.

Figure 3

Simulated snow, melt-pond depth, and bare ice area, and simulated and observed PAR. Time series of (a) surface area fraction of simulated snow (red), melt ponds (green), and bare ice (blue) and (b) simulated daily-mean (line) and observed (circles) under-ice PAR during the Arctic-ICE 2010 study period. In (b), the units for the simulated PAR values were converted from W m–2 to µmol photons m–2 s–1 by a conversion factor of 4.56 following Lavoie et al. (2005). Vertical bars associated with the solid line represent the diurnal range of simulated under-ice PAR. Red and blue circles represent the daily-mean values measured using tethers deployed over high (HSC) and low (LSC) snow cover sites, respectively. Yellow circles are the instantaneous values based on CTD casts (CTD). DOI: https://doi.org/10.1525/elementa.229.f3

#### 3.1.3 Sea ice ecosystem

Figure 4 shows the simulated time series of sea ice ecosystem variables. The standard run simulated an ice algal bloom that is comparable to the observations in terms of both the magnitude and timing of the bloom (Figure 4a). In the following, we discuss the dynamics of simulated sea ice ecosystem by partitioning into growth and decline phases.

Figure 4

Simulated and observed ice algal biomass, nutrients, growth limitations, and simulated sympagic and pelagic production. Time series of (a) simulated (line) and observed (circles) Chl a concentrations in the bottom 3 cm of the sea ice, (b) simulated nitrate (solid black), ammonium (dashed black) and silicate (red) concentrations in the bottom 3 cm of sea ice, (c) simulated daily-mean growth limitation index for light (yellow), nitrogen (black), silicate (red), and ice melting (green), and (d) primary production rates of simulated ice algae (solid line) and phytoplankton (dashed line). In (a), circles represent the site-average values with one standard deviations indicated by vertical bars. DOI: https://doi.org/10.1525/elementa.229.f4

The growth phase of simulated ice algal bloom lasted from late March to mid-May, while the bloom decline phase is from mid-May to late June. During the growth phase of the ice algal bloom, the simulated ice algal biomass in the standard run increased up to 1050 mg Chl a m–3 (Figure 4a). This maximum value in the bloom is within the range of observed values during the peak of the ice algal biomass (800 – 1300 mg Chl a m–3). Note that this wide range in the observed peak is due to sampling over different snow depth conditions, and that the model succeeded in simulating a bloom that falls near the center of the observed range. Up until the end of April, concentrations of simulated nitrate and silicate in the ice decreased rapidly due to uptake by ice algae, while the simulated ammonium concentration increased as a result of remineralization of dead ice algal cells (Figure 4b). During this time, the ice algal growth rate declined slightly even though nutrients are not yet limiting, likely due to the quadratic term in the parameterization of mortality. Consequently, both nitrate and silicate concentrations recovered slightly until they were drawn down further by ice algae during their bloom peak in mid-May. The ice algal growth was generally light-limited during the growth phase (Figure 4c), except for a day in the beginning of May when the nitrate concentration reached nearly 0.5 mmol m–3 (Figure 4b).

At the peak of the ice algal bloom, simulated nutrients became extremely low, nearly 0 mmol m–3 for nitrate and ammonium and 1 mmol m–3 for silicate (Figure 4b). Consequently, the ice algal growth became nitrogen-limited following the peak (Figure 4c), and remained so until the end of the bloom in late June (Figure 4a). The simulated range of nitrate concentration (0 – 8 mmol m–3; Figure 4b) matches with the observed range reported in Galindo et al. (2014). In contrast, the simulated range of ammonium concentration (0 – 0.05 mmol N m–3) is much smaller than the range typically observed in the bottom ice (e.g., Vancoppenolle et al., 2013). This discrepancy is most likely due to the fact that much of the ammonium found in the bottom ice is trapped in the ice matrix and therefore not accessible to ice algae residing in the brine phase of the ice (Vancoppenolle et al., 2013). The model simulates the remaining fraction of ammonium available for ice algae which is low in abundance due to rapid turnover of ammonium production and removal processes. Figure 4d presents the time series of depth-integrated production rates by simulated ice algae and phytoplankton (i.e., sum of P1 and P2). The production rate of simulated ice algae was around 0.1 g C m–2 d–1 during its bloom peak in mid-May. The time-integrated production by ice algae and phytoplankton over the simulation period was about 4 and 60 g C m–2, respectively. Hence, the primary production by simulated ice algae accounted for 6% of the entire sea ice and water column primary production. This fraction is within the range of the observational and model estimates for first-year Arctic sea ice (2 – 33%; Legendre et al., 1992; Gosselin et al., 1997; Lavoie et al., 2009).

#### 3.1.4 Pelagic ecosystem

Figure 5 shows the comparison of simulated and observed time series of chlorophyll a concentrations in the upper 80 m of the water column. In mid-June, the model simulated an under-ice phytoplankton bloom in the upper 10 m of the water column (Figure 5a). This bloom was dominated by large phytoplankton (Figure S1b), and reached a peak concentration of 13 mg Chl a m–3 in late June. The timing, magnitude, and vertical extent of the simulated under-ice phytoplankton bloom are consistent with the observed bloom (Figure 5b), which was also dominated by large cells (Mundy et al., 2014). The simulated bloom migrated downward and formed a subsurface chlorophyll maximum of 18 mg Chl a m–3 at 15 – 40 m between late June and early July. During the ice-free period, increased light penetration allowed the deepening of the simulated subsurface chlorophyll maximum to a depth of 75 m where it maintained fairly large concentrations (above 6 mg Chl a m–3) until the end of August. The formation and subsequent deepening of a deep chlorophyll-maximum is a typical feature in the Arctic where surface nutrients are low (the chlorophyll maximum typically follows the nitricline). No direct observations are available for this particular time period near Resolute to evaluate the deepening of the subsurface chlorophyll maximum simulated by the model. However, observations taken during the last decade in the Beaufort Sea and Canadian Archipelago show the subsurface chlorophyll maxima with depths ranging from 35 and close to 100 m depending on time and location measured (Tremblay et al., 2008; Carmack et al., 2010; McLaughlin and Carmack 2010; Carmack and McLaughlin, 2011) which is also represented in model results (Steiner et al., 2015). The chlorophyll maximum in the Chukchi Sea tends to be much shallower (Brown et al., 2015–2), while the deepest maxima have been observed in the Beaufort Sea. A maximum depth of 75 m for the deep chlorophyll maximum in the Canadian Arctic Archipelago is within the range of observations. The daily production rates corresponding to the under-ice phytoplankton bloom (1.2 g C m–2 d–1) and the subsurface chlorophyll maximum (up to 1.6 g C m–2 d–1) simulated by the model (Figure 4d) are comparable to the observed rates in Resolute Passage (1.1 g C m–2 d–1; Mundy et al., 2014) and in the Beaufort Sea (1.4 g C m–2 d–1; Mundy et al., 2009), respectively.

Figure 5

Simulated and observed Chl a concentration. Time series of (a) simulated and (b) observed Chl a concentrations in the upper 80 m of the water column. DOI: https://doi.org/10.1525/elementa.229.f5

Figure 6a–c illustrates the temporal evolution of simulated dissolved nutrients in the upper 80 m of the water column. Prior to the development of the under-ice phytoplankton bloom in mid-June (Figure 5a), the concentrations of simulated nitrate (Figure 6a) and silicate in the upper 15 m (Figure 6c) were reduced as a result of the uptake by ice algae. In contrast to nitrate and silicate, concentrations of simulated ammonium increased slightly below the nitracline due to the remineralization of dead ice algal cells released into the water column (Figure 6b). In late June, these nutrients were drawn down by large phytoplankton, and decreased to <1 mmol m–3 (nitrate; Figure 6a), <0.04 mmol m–3 (ammonium; Figure 6b), and <4 mmol m–3 (silicate; Figure 6c) in the upper 10 m of the water column. These values of simulated nitrate and silicate concentrations are close to the values (0.2 mmol m–3 for nitrate+nitrite and 2.8 mmol m–3 for silicate) reported at the end of the sampling period (21 June) in Resolute Passage (Mundy et al., 2014). The concentrations of simulated nutrients remained below these levels until the end of the simulation period (Figure 6a–c) because large detritus, which consists of dead cells of ice algae and large phytoplankton and fecal pellets, sank quickly (50 m d–1 as specified in the model, following Lavoie et al., 2009) into the deeper water column before they could be remineralized in the upper water column. The rapid sinking of large detritus resulted in the accumulation of ammonium at depth below the nitracline in late June onwards (Figure 6b).

Figure 6

Simulated water column concentrations of nutrients and biological uptake and drawdown of nitrate. Simulated time series of (a) nitrate, (b) ammonium, and (c) silicate concentrations in the upper 80 m of the water column (depth of entire water column is 141 m). (d) Simulated time series of cumulative depth-integrated nitrate uptake and drawdown. In (d), areas filled in red represent the cumulative uptake by ice algae integrated over the bottom 3 cm of the ice skeletal layer, areas filled in blue represent the cumulative uptake by phytoplankton (P1 and P2) integrated over the upper 80 m of the water column, and the black line represents the cumulative amount of nitrate drawn down from the upper 80 m of the water column. Note that the sum of the two uptake terms (red+blue) does not balance with the drawdown during the ice-free period; the mismatch represents the uptake of nitrate entrained from the layer below 80 m. DOI: https://doi.org/10.1525/elementa.229.f6

To demonstrate that the ice algal uptake and the nutrient removal in the water column are balanced in the model, the time series of depth-integrated (3 cm) cumulative nitrate uptake by ice algae is displayed with the depth-integrated cumulative nitrate drawdown and total uptake by phytoplankton in the upper 80 m of the water column (Figure 6d). Clearly, the total amount of nitrate consumed by ice algae is equivalent to the amount removed from the water column until the onset of the pelagic bloom in mid-June. The result demonstrates an important role of ice algae in reducing the ambient nutrients in the upper water column. This important aspect of sympagic-pelagic ecological coupling will be examined further in a later section. The decreasing trend of simulated nitrate in the water column during May and June (Figure 6a) is generally in good agreement with the observed nitrogen time series in the ice and underlying water column as reported in Galindo et al. (2014).

### 3.2 Sympagic-pelagic ecosystem coupling

In order to assess the impact of the simulated ice algal bloom on the underlying pelagic ecosystem, we conducted an additional simulation that turned off the ice algal bloom (referred to as the exclusion run). This scenario was established by setting the initial biomass of ice algae to zero, while all other setups are identical to the standard run. Hence, the difference in the results between the standard and the exclusion runs represents the impact of ice algae on the pelagic ecosystem.

Figure 7 displays the comparison of the two runs in terms of Chl a concentrations in the upper 50 m of the water column. The differences between the two runs are most evident in late June, which correspond to the under-ice bloom in the upper 10 m of the water column (Figure 7c). Both the timing and magnitude of the bloom were affected by the presence/absence of ice algae. When ice algae were excluded from the simulation (Figure 7b), the onset of the under-ice bloom was delayed by a few days. This delay in the bloom onset is due to the lack of seeding by ice algae in the exclusion run (Hayashida et al., 2017). Despite the delay in the development of the under-ice bloom, the exclusion run simulated a higher peak in Chl a (with a concentration difference of about 7 mg Chl a m–3) than the standard run. The enhanced peak in the exclusion run is due to the absence of nutrient drawdown by ice algae (Figure 8), which makes a concentration difference of about 3 mmol N m–3 in the upper 10 m of the water column at the onset of the under-ice bloom. (It is not due to the absence of light-shading by the ice algae, as the pelagic bloom does not begin in the standard run until after the ice algal bloom has ended). The effects of ice algae in the pelagic ecosystem appear to be relatively small below the upper 10 m of the water column, as there is no substantial difference in either Chl a or nitrate concentrations between the standard and the exclusion runs.

Figure 7

Water columm Chl a concentration when ice algae are present, absent, and the difference. Simulated phytoplankton bloom in the upper 50 m of the water column when ice algae are present (a), absent (b), and the difference (c). Phytoplankton are sum of large and small (P1 and P2) groups. DOI: https://doi.org/10.1525/elementa.229.f7

Figure 8

Water column nitrate when ice algae are present, absent, and the difference. Simulated NO3 concentration in the upper 50 m of the water column when ice algae are present (a), absent (b), and the difference (c). DOI: https://doi.org/10.1525/elementa.229.f8

#### 3.2.1 Sinking rate of large detritus

In the model, large detritus (D2) represents the non-living particulate matter originating mainly from ice algal and large phytoplankton cells. The simulated large detritus is assumed to sink at a constant rate (wd2; Appendix A2) which is presumably faster than the sinking rate of small detritus, another form of detritus considered in the model. In the standard run, a sinking rate of 50 m d–1 was prescribed for large detritus following Lavoie et al. (2009). However, observations of this rate span a range of values. Onodera et al. (2015) observed sinking rates from 37 to more than 85 m d–1 for diatoms in the western Arctic Ocean. Higher and lower rates have also been measured, with sinking rates well over 100 m d–1 among Antarctic ice algal aggregates (Sibert et al., 2010) and near 20 m d–1 in lab tests with the common Arctic ice algae diatom Nitzschia frigida (Aumack and Juhl, 2015).

In this sensitivity analysis, we assessed the simulated phytoplankton response to a change in the fast sinking rate. Runs with a slower sinking rate do not show much difference in the pelagic ecosystem until the sinking rate is lowered below a threshold of approximately 10 m d–1. Above this threshold, large detritus is effectively removed from the euphotic layer and transported to depth before it can be remineralized (Figure 9a and b). Below that threshold, e.g., at wd2 = 5 m d–1, a secondary sub-surface bloom, comprised of small phytoplankton (P1), forms after the first bloom (Figure 9c). This secondary bloom results from an increased supply of nitrogen. The remineralization rate from biogenic silica is an order of magnitude slower than that from D2 (0.01 d–1, and 0.3 d–1, respectively, Table S2), and hence the second bloom does not allow for silicate-dependent large phytoplankton.

Figure 9

Phytoplankton in the water column with fast-sinking detritus. Simulated phytoplankton in the upper 50 m of the water column, with fast-sinking detritus (D2) set at 50 m d–1(a), 15 m d–1(b), and 5 m d–1(c). First bloom is dominated by large phytoplankton (P2, diatoms) and the later bloom in (c) is dominated by small phytoplankton (P1, flagellates). DOI: https://doi.org/10.1525/elementa.229.f9

### 3.3 Sensitivity analyses for ice algae

Given the influence of simulated ice algae on the underlying pelagic ecosystem, it is of utmost interest to investigate the physical and biological controls on the simulated ice algal bloom (and subsequently on the underlying ecosystem). The shape of these control functions is set via parameter values which are often not measured directly, but inferred from concentrations of observed variables that are also not well constrained. Sensitivity analyses focus on parameters that represent key uncertainties in the observational record. By varying each parameter over the range of observed (or estimated if not constrained by observations) uncertainty and determining which parameters have the strongest impact on properties of the simulated ice algal bloom, we can identify which observations would be most beneficial to improve our understanding of the system.

The growth of the ice algal bloom is dependent on both physical and biogeochemical processes. In the simulated ice algal bloom, several key parameters determine the strength of these influences. In the standard simulation, parameters controlling the onset, growth, maximum biomass, and termination of the modelled ice algae have been adjusted to result in good agreement with observations. In this section, key parameters associated with over-wintering (pre-bloom) ice algal biomass, mortality, photosynthetic sensitivity, and nutrient limitation, are varied independently in order to determine the sensitivity of the simulated bloom.

The experiments testing photosynthetic efficiency (not shown) demonstrated that increasing photosynthetic efficiency does not increase the maximum biomass substantially, because of nutrient limitation. The experiments varying the ratio of intracellular silicate to nitrogen (also not shown) indicated that increasing the intracellular ratio Si:N by ~20% was enough for the ice algal growth to become silicate-limited instead of nitrogen-limited.

#### 3.3.1 Pre-bloom algal biomass in the ice

In previous work, the pre-bloom ice algal biomass in simulations has been estimated based on water column measurements during ice formation (Steiner et al., 2016), or from early bloom measurements (Lavoie et al., 2005; Pogson et al., 2011). It is possible that processes involved in ice formation can preferentially pick up certain marine particles, such as algal cells, and ice algal biomass concentrations up to 2 orders of magnitude higher than the underlying water column have been observed in sea ice in fall and winter (Różańska et al., 2008; Niemi et al., 2011).

The year-to-year variability in the amount of ice algae before the bloom may have a strong effect on the timing of the bloom onset. The timing of the onset of the simulated ice algal bloom (defined as when the biomass surpasses 100 mg Chl a m–3) depends on the pre-bloom, or over-wintering, concentration (Figure 10b). The pre-bloom concentration is implemented in the model as a minimum ice algal biomass. In the standard run, the pre-bloom concentration was set at 1 mmol N m–3 (or 3.533 mg Chl a m–3) to match the observed bloom onset. This value is approximately 20% of the value used in Lavoie et al. (2009) of 0.5 mg Chl a m–2 (16.7 mg Chl a m–3, assuming a 3 cm ice algal layer). This value was set for the more productive Beaufort Sea, but is an order of magnitude larger than that observed by Niemi et al. (2011) (0.1 mg Chl a m–3 for first-year ice in the Beaufort Sea). Because of the large difference between these estimates of the pre-bloom concentration, the simulation was run with pre-bloom concentrations at 200%, 150%, 50%, and 10% of the standard value in order to test the importance of this value on the onset and maximum concentration of the ice algal bloom. With pre-bloom concentrations of 50% (200%) of the standard value, the subsequent ice algal blooms are slightly later (earlier), with the biomass reaching 100 mg Chl a m–3 approximately 4 days later (earlier). It is evident that a multiplicative change in the pre-bloom biomass results in an additive offset in the time needed to reach a specified biomass (consistent with exponential growth through the earlier parts of the bloom).

Figure 10

Snow and ice thickness and ice algal biomass, varying pre-bloom biomass. Snow and ice thickness (cm) and ice algal biomass during sensitivity analyses of the simulated ice algal bloom to variation of pre-bloom biomass (b), prescribed at 10, 2, 1/2, and 1/10 times that in the standard simulation (solid black line). DOI: https://doi.org/10.1525/elementa.229.f10

Modelled ice algal blooms for runs with the pre-bloom ice algal biomass values an order of magnitude larger or smaller than the standard value (Figure 10b) indicate that, at higher pre-bloom biomass, the bloom occurs earlier, but the maximum biomass is not much greater than in the standard run because the growth is terminated by nutrient limitation. When the pre-bloom ice algal biomass is one-tenth of that in the standard run, the timing of the bloom onset is delayed and the bloom levels off (at the time of the maximum biomass in the standard run). This is because the NO3 limitation in that time period is approximately 0.2 day–1 (not shown). In an idealized 12 hour day, and no other limitation, the daily averaged minimum limitation would be half of that (0.1 day–1), which is roughly equal to the grazing rate.

These results are in agreement with those of Jin et al. (2006), who found that doubling the initial ice algal biomass does not affect the maximum biomass of the bloom, and results in an onset 3 – 5 days earlier. In addition, Pogson et al. (2011) find that using the observed low initial biomass under high snow cover causes underestimation of the simulated maximum biomass when compared to the observations.

#### 3.3.2 Mortality

The mortality rate for marine algae is commonly parameterized as some combination of linear and quadratic dependencies on biomass. To our knowledge mortality rates have not been directly measured for ice algae and the contribution of linear and quadratic contributions needs to be tested. Here, the mortality rate for ice algae (M; Appendix A2) is defined as a function of biomass:

$M={m}_{\mathit{\text{lia}}}\mathrm{exp}\left({b}_{\mathit{\text{ia}}}{\left[T\right]}_{\mathit{\text{ia}}}\right)+{m}_{\mathit{\text{qia}}}\left[\mathit{\text{IA}}\right]$
(1)

where bia, [T]ia, and [IA] represent the temperature sensitivity coefficient, temperature in the ice skeletal layer, and ice algal biomass, respectively (see the Appendix for details). mlia represents the rate constant for the temperature-dependent linear mortality and mqia is the rate constant for the quadratic mortality. The linear term represents ice algal biomass-independent processes, in which a specified fraction of the population is lost per unit time. Lavoie et al. (2005) defined this term as the grazing rate on ice algae, and prescribed it at 10% of the growth rate. The quadratic term represents crowding effects, in that the fraction of biomass lost per unit time increases with higher biomass. (Although the quadratic formulation is a commonly used approach in representing the crowding effect of large phytoplankton cells (i.e., diatoms) in marine ecosystem models (e.g., Steiner et al., 2006; Aumont et al., 2015), we do not implement it in the case of small plankton because they do not reach high enough densities). Based on the model tuning to match observations, mlia and mqia are respectively set to 0.03 d–1 and 0.00015 d–1 in the standard run. As the simulated bloom grows, the population will have a quasi-exponential growth if the linear contribution to mortality varies slowly with time, and the biomass is small enough that the quadratic contribution to mortality is small.

Figure 11b presents the standard run along with multiple runs in which the linear and quadratic mortality parameters have been increased or decreased. As expected, when both parameters are increased (decreased), the simulated ice algae has a lower (higher) maximum biomass than the standard run. When the two are changed in opposite directions, the magnitude of the maximum biomass does not vary substantially, but the onset timing is earlier or later. In Figure 11c, the red box from Figure 11b is enlarged in order to show when the simulated ice algal bloom crosses the 100 mg Chl a m–3 threshold. With a 25% decrease (increase) to this parameter, the bloom reaches the 100 mg Chl a m–3 threshold 2 days earlier (later).

Figure 11

Snow and ice thickness, ice algal biomass varying mortality function, and onset of the bloom. Snow and ice thickness (cm) (a) and ice algal biomass (mg Chl a m–3) differing linear and quadratic dependencies on mortality (b). The black solid line in (b) is the standard run, the dashed red (blue) line is the simulated bloom with both linear and quadratic dependencies decreased (increased) by 25%. The solid colored lines are for blooms with linear and quadratic dependencies changed in opposite directions, e.g., increased for linear and decreased for quadratic. The onset of the bloom in the red box in (b) is expanded in (c). DOI: https://doi.org/10.1525/elementa.229.f11

Dashed lines correspond to simulations in which the linear and quadratic mortality parameters have been changed in the same way. These different runs cross the 100 mg Chl a m–3 threshold at almost the same time, indicating that the bloom onset (when ice algal biomass is small) is relatively insensitive to the quadratic mortality dependency. Therefore, these two mortality parameters can be adjusted to best fit observations for the timing of the bloom onset and magnitude of maximum biomass.

## 4 Discussion

The recent model study by Abraham et al. (2015) showed that grid-cell mean simulations of light and heat fluxes through sea ice could be improved by parameterizing the subgrid-scale snow depth variability, relative to simulations with spatially uniform snow depth distribution. These authors pointed out the need to examine biological responses to this new parameterization. In the first part of the present study, we investigated the impact of the new parameterization on simulated ice algae. The results indicate an improvement in simulating the ice algal bloom especially during the melt period, owing to an improvement in simulating the gradual increase in light availability to the ice algae. However, in this study, we are unable to further assess the impact of the new light parameterization on earlier stages of the bloom because the observed time series of ice algal biomass are confined mostly to the decline phase of the bloom. Measurements focusing on ice algal biomass during the onset and early growth of blooms are needed for assessing this impact.

As discussed in Arrigo (2014), the presence of ice algae affects several important processes in the underlying water column ecosystem. However, it is logistically difficult to isolate the contribution of ice algae from that of phytoplankton in terms of observed nutrient drawdown and biomass production. It is similarly difficult to observationally assess the seeding of the phytoplankton bloom by ice algae. Hence, process models become important tools to address questions like: What if ice algae were excluded from a given environment? In particular, the absence of advective processes in 1-D models allows focus on the in situ sympagic-pelagic ecosystem coupling. The present analysis demonstrated some of the influences of ice algae on the pelagic ecosystem. The results indicate that both the timing and magnitude of the simulated under-ice phytoplankton bloom are affected by the presence of ice algae. The timing of the bloom is affected due to seeding as a result of ice algal flushing, whereas the magnitude is affected due to the nutrient drawdown by the earlier ice algal bloom. These impacts of ice algae further influence other important biogeochemical processes, such as the production of dimethylsulfide (Hayashida et al., 2017). Previous model studies also indicated the timing and magnitude of the ice-associated pelagic bloom as an important response to ice algal seeding Jin et al. (2007); Tedesco et al. (2012). However Jin et al. (2007) highlighted the importance of stratification on the response suggesting that sudden mixing events following ice melt would disrupt the ice-associated pelagic bloom. More quantitative estimates for the effects of ice algae on the underlying ecosystem can be achieved by conducting simulations (including the exclusion run) in a full 3-D model using the parameterizations considered in this study. Deal et al. (2011) and Jin et al. (2012) 3-D model applications highlight both high regional variability as well as the seasonal importance of ice algal primary production.

The model applies several simplified assumptions due to lack of observations in the ice. For instance, the simulated ice algal nitrogen uptake preference (${p}_{\mathit{\text{no}}3}^{\mathit{\text{ia}}}$ in Equation 27) is constant throughout the simulation. However, Harrison et al. (1990) observed that nitrogen utilization by ice algal communities of Barrow Strait shift from a nitrate- to an ammonium-dominated metabolism. In addition, the nitrification rate (NH4 to NO3) in sea ice is set to a constant rate, even though bacteria in the ice that facilitate the process (Fripiat et al., 2014) may experience variations due to environmental fluctuations.

In both the modelled ice and water column, nutrient depletion due to phytoplankton uptake leads to near-zero concentrations in the limiting nutrient. Observations of post-bloom nutrient concentrations in an area with little horizontal transport could allow assessment of this result. The influence of horizontal advection on the nutrient drawdown below the ice could be assessed in 3-D model simulations.

The model results indicate that a combination of linear and quadratic mortality terms is required to adequately represent the development and decline of the ice algal bloom. The application of a quadratic mortality term implies a larger specific mortality at higher ice algal concentrations, representing lysis due to viral infection and other overcrowding processes that occur at higher ice algal concentrations. Additional field observations during the height of the bloom could help to constrain this term.

In the standard simulation, the growth of ice algae was initially limited by light, and then by nutrients (nitrate) during the peak and the decline of the bloom, which is consistent with the findings of previous studies (Mundy et al., 2014). The simulated under-ice bloom was similar to the observed bloom in terms of the magnitude, timing, and the species composition (dominated by diatoms, Galindo et al. (2014)). During the ice-free period, the simulated under-ice bloom was succeeded by the formation of a subsurface chlorophyll maximum. While this is a common feature in low-nutrient Arctic waters, observations are lacking for this particular time and location. It is possible that high tidal and/or horizontal mixing could prevent a deep chlorophyll maximum from developing in particular regions.

The parameters were adjusted to this specific dataset (particular year, particular place). Applications for different years and locations, and subsequent implementation in a 3-D model, will indicate if some retuning may be necessary. A need for retuning would hint at processes that are incompletely understood and indicate whether further measurements to constrain the process are required.

## 5 Conclusions

This 1-D study is intended as a step in the development of a 3-D model, one of a growing number that incorporate biogeochemical processes in order to represent the sympagic ecosystem and its coupling to the underlying pelagic ecosystem.

In order to establish a set of parameterizations which can be transferred into a 3-D regional Arctic model which couples sea-ice, ocean and associated ecosystems, this 1-D model study investigates the physical and biological controls on sympagic and pelagic primary production using observations from Resolute Passage. Results of the standard simulation, including a snow distribution function allowing for a slow evolution towards bare ice and melt ponds, were generally in good agreement with the variability of snow/melt pond depths, ice thickness, under-ice PAR, and bottom-ice and seawater Chl a observed during the melt season in 2010. The simulated ice algal and under-ice phytoplankton blooms in the standard run were in reasonable agreement with the observations in terms of timing and magnitude.

Several findings can be taken from the sensitivity analyses. (1) Ice algal growth limits subsequent pelagic biomass in the upper water column by removing nutrients and limiting their availability to the phytoplankton, with a decrease of ~50% of the maximum phytoplankton concentration in the upper 10 m in the standard run relative to the run without ice algae. (2) Photosynthetic sensitivity and pre-bloom biomass determine the onset timing of the ice algal bloom. (3) The maximum biomass is relatively insensitive to the pre-bloom ice algal biomass. (4) A combination of linear and quadratic parameterizations of mortality rate is required to adequately simulate the evolution of the ice algal bloom, indicating that processes associated with each of these functional forms are occurring within the ice algal bloom phase. And (5), a large detrital (D2) sinking rate greater than a threshold of ~10 m d–1 effectively strips the upper water column of the potential to regenerate the limiting nutrient after the bloom by transporting it to depth. For this scenario a deep chlorophyll maximum develops, as is characteristic for low nutrient Arctic waters. A D2 sinking rate slower than this threshold allows for a subsequent subsurface P1 bloom due to availability of remineralized ammonium (from detritus) after the initial (P2 dominated) pelagic bloom.

Measurements needed to better constrain the simulated ice algal bloom include ice algal concentration in the winter, in situ mortality rate, and sinking rates for ice algal aggregates. This 1-D study is part of two subsequent 1-D studies, implementing sulfur (dimethyl sulfide, or DMS) and inorganic carbon cycles. The work in all three of these studies will be used as a basis for the implementation of ice algae, DMS, and carbon cycles into a 3-D coupled ice-ocean biogeochemical regional model of the Canadian marine Arctic.

## 6 Appendix

### 6.1 A1. Parameterizations for subgrid-scale snow depth distribution and light penetration through snow, sea ice, and melt ponds

To improve light and heat flux estimates through sea ice in regional and global models, Abraham et al. (2015) applied two kinds of one-parameter probability density functions for describing subgrid-scale snow depth variability: Rayleigh and gamma distributions. In this study, the Rayleigh distribution is used in model simulations since Abraham et al. (2015) indicated a better fit with observed snow depth evolution. The probability density function for the Rayleigh distributed snow (pdf(h)) is defined as:

$\begin{array}{l}\mathit{\text{pdf}}\left(h\right)=\frac{\pi h}{2{h}_{s}^{2}}\mathrm{exp}\left(-\frac{\pi {h}^{2}}{4{h}_{s}^{2}}\right),\hfill \end{array}$
(2)

where h is the snow depth (m) and hs represents the gridbox-mean snow depth (m), which is simulated by the sea ice model. The light transmission through snow and sea ice is described by the Beer-Lambert law, which is defined in a generalized form as:

$\begin{array}{l}I\left(z\right)={I}_{0}\mathrm{exp}\left(-\kappa z\right),\hfill \end{array}$
(3)

where I(z) is the radiation at depth z (m) below the surface (W m–2) and κ is the extinction coefficient of the medium (m–1). I0 represents the amount of incident light that penetrates into the snow/ice/melt ponds surface (W m–2), which is defined in a generalized form as:

${I}_{0}=\mathit{\text{SWR}}\left(1-\alpha \right){i}_{0},$
(4)

where SWR is the incident shortwave radiation (W m–2), α is the surface albedo (dimensionless), and i0 is the transmissivity coefficient (dimensionless). For a fully snow-covered surface of non-uniform snow depth, the gridbox-mean light intensity at the snow base is obtained by averaging the Beer-Lambert law over all snow depths weighted by the relative probabilities:

${\overline{I}}_{s}=\mathit{\text{SWR}}\left(1-{\alpha }_{s}\right){i}_{0,s}{\int }_{0}^{\infty }pdf\left(h\right)\mathrm{exp}\left(-{\kappa }_{s}h\right)\text{d}h,$
(5)

where αs, i0,s, and κs respectively represent the albedo, the transmissivity coefficient, and the extinction coefficient for snow.

During melt periods, the ice surface may have different covers, such as snow, bare ice, and melt ponds. To account for different surface conditions in the parameterization of the gridbox-mean light intensity at the ice surface, Abraham et al. (2015) accounted for these different surface conditions within a grid box by introducing surface area fractions of snow (As), bare ice (Ai), and melt ponds (Am), such that:

$\begin{array}{l}{A}_{s}+{A}_{i}+{A}_{m}=1.\hfill \end{array}$
(6)

The parameterizations for As, Ai, and Am are described in Abraham et al. (2015). The gridbox-mean light intensity at the ice surface (Ītop) is then defined as a sum of the incident light that has: 1) penetrated through snow; 2) reached the bare ice; and 3) penetrated through melt ponds. Hence, Ītop is given by:

${\overline{I}}_{\mathit{\text{top}}}={A}_{s}{\overline{I}}_{s}+{A}_{i}{I}_{0,i}+{A}_{m}{I}_{0,m}\mathrm{exp}\left(-{\kappa }_{m}{h}_{m}\right)$
(7)

where Īs is Equation 5. κm and hm are the extinction coefficient and the depth of melt ponds, respectively. I0,i and I0,m respectively represent the amounts of incident light that penetrates through the ice and melt ponds surface:

${I}_{0,i}=\mathit{\text{SWR}}\left(1-{\alpha }_{i}\right){i}_{0,i}$
(8)
${I}_{0,m}=\mathit{\text{SWR}}\left(1-{\alpha }_{m}\right){i}_{0,m}$
(9)

where αi and αm are the albedos and i0,i and i0,m are the transmissivity coefficients for sea ice and melt ponds, respectively. The optical parameters used in this study are listed in Table S1. Note that different values for the extinction coefficients and albedos are set between the freezing and melting phases of snow and sea ice. To allow a smooth transition between the values under the freezing and melting phases, the extinction coefficients and albedos of snow and sea are defined based on Abraham et al. (2015) as:

${\kappa }_{s}=\frac{{\kappa }_{s,f}+{\kappa }_{s,m}}{2}+\frac{{\kappa }_{s,m}-{\kappa }_{s,f}}{2}\mathrm{tanh}\left(T-{273.15}^{°}\text{C}\right)$
(10)
${\kappa }_{i}=\frac{{\kappa }_{i,f}+{\kappa }_{i,m}}{2}+\frac{{\kappa }_{i,m}-{\kappa }_{i,f}}{2}\mathrm{tanh}\left(T-{273.15}^{°}\text{C}\right)$
(11)
${\alpha }_{s}=\frac{{\alpha }_{s,f}+{\alpha }_{s,m}}{2}+\frac{{\alpha }_{s,m}-{\alpha }_{s,f}}{2}\mathrm{tanh}\left(T-{273.15}^{°}\text{C}\right)$
(12)
${\alpha }_{i}=\frac{{\alpha }_{i,f}+{\alpha }_{i,m}}{2}+\frac{{\alpha }_{i,m}-{\alpha }_{i,f}}{2}\mathrm{tanh}\left(T-{273.15}^{°}\text{C}\right)$
(13)

Following Zeebe et al. (1996), it is assumed that only PAR penetrates into the ice interior, while the radiation outside of PAR bands is absorbed by the uppermost layer of snow, bare ice, or melt ponds. Therefore, the gridbox-mean PAR at the ice base $\left({\overline{PAR}}_{bot}\right)$ is defined as:

${\overline{\mathit{\text{PAR}}}}_{\mathit{\text{bot}}}={\overline{I}}_{\mathit{\text{top}}}\mathrm{exp}\left(-{\kappa }_{i}{h}_{i}\right)$
(14)

where hi is the sea ice thickness. Finally, the gridbox-mean PAR in the water column under the ice $\left({\overline{\mathit{\text{PAR}}}}_{\mathit{\text{ui}}}\left(z\right)\right)$ is first attenuated by ice algae before it reaches the uppermost layer of the water column, and is further reduced as it penetrates through each model layer due to absorption and scattering by seawater itself, as well as by phytoplankton and detritus:

${\overline{\mathit{\text{PAR}}}}_{\mathit{\text{ui}}}\left(z\right)={\overline{\mathit{\text{PAR}}}}_{\mathit{\text{bot}}}\mathrm{exp}\left(-{\kappa }_{\mathit{\text{ia}}}\left[\mathit{\text{IA}}\right]{h}_{\mathit{\text{ia}}}-{\kappa }_{\mathit{\text{ui}}}{z}_{\mathit{\text{sw}}}\right)$
(15)

where κia is the extinction coefficient of ice algae, [IA] is the ice algal biomass, and hia is the thickness of the ice skeletal layer, which are all defined in the next section. zsw is the depth of seawater under the ice. κui is the total extinction coefficient in the water column defined as:

${\kappa }_{ui}={\kappa }_{sw}+{\kappa }_{pd}{\int }_{0}^{{z}_{sw}}\left(P1+P2+D1+D2\right)\text{d}z$
(16)

where κsw and κpd are the extinction coefficients for seawater and for both phytoplankton and detritus. κsw is computed by the GOTM model assuming the Jerlov type I (Burchard et al., 2006). The concentrations of phytoplankton (P1, P2) and detritus (D1, and D2) at the given model layer are defined in the next section.

### 6.2 A2. Ecosystem model equations

The coupled sea ice-ocean biogeochemical model consists of 14 state variables (Figure 1 and Table S1). Nitrogen is used as the currency for the model state variables other than [Si]ia, [Si], and [BSi], which are expressed in silicon units. For comparison with observations, conversion from nitrogen (N) to Chl a and carbon (C) are required. For small phytoplankton (P1), the Redfield C:N ratio of 106:16 (mol:mol) following Redfield et al. (1963) and a fixed C:Chl a ratio of 50:1 (wt:wt) following Lavoie et al. (2009) are used. For ice algae and large phytoplankton (IA and P2), a fixed C:N ratio of 106:12 (mol:mol) following Palmer et al. (2014) and a fixed C:Chl a ratio of 28:1 (wt:wt) following Lavoie et al. (2009) were used. In Lavoie et al. (2009), the C:Chl a ratio of 28:1 (wt:wt) was used for ice algae, while 50:1 (wt:wt) was used for phytoplankton (corresponding to P2 in this study) assuming that ice algae was more acclimated to low light conditions. However, in this study, prescribing different C: Chl a ratios between IA and P2 would result in the violation of mass conservation of chlorophyll a in P2 because IA and P2 are coupled through the flushing of ice algae to enter the P2 pool, and also because chlrophyll a is not modeled explicitly. We represent the photoacclimation of IA relative to P2 by prescribing different values of photosynthetic parameters. The parameters for the ecosystem model are listed in Table S3.

#### 6.2.1 Sea ice component

The sea ice biogeochemical model consists of four state variables: ice algae, nitrate, ammonium, and silicate. The equation for ice algae is:

$\frac{\partial }{\partial t}\left[\mathit{\text{IA}}\right]={\mu }_{\mathit{\text{ia}}}\left[\mathit{\text{IA}}\right]-M\left[\mathit{\text{IA}}\right]-\left(\frac{\left[{R}_{\mathit{\text{ice}}}\right]+{f}_{\mathit{\text{pond}}}{r}_{\mathit{\text{pond}}}}{{z}_{\mathit{\text{ia}}}}\right)\left[\mathit{\text{IA}}\right],$
(17)

where [IA] is the ice algal biomass in the bottom 3 cm of the skeletal layer (mmol N m–3). The first term in Equation 17 represents the growth of ice algae due to primary production. The rate of ice algal growth (µia) depends on the ambient temperature of the bottom ice ([T]ia in °C) following Eppley (1972), as well as on the minimum value among the four limitation functions (Lnit, Lsil, Llig, and Lice):

${\mu }_{\mathit{\text{ia}}}={\mu }_{\mathit{\text{ia}}}^{\mathit{\text{max}}}\mathrm{ln}\left(2\right)\mathrm{exp}\left({b}_{\mathit{\text{ia}}}{\left[T\right]}_{\mathit{\text{ia}}}\right)\mathrm{min}\left({L}_{\mathit{\text{nit}}},{L}_{\mathit{\text{sil}}},{L}_{\mathit{\text{lig}}},{L}_{\mathit{\text{ice}}}\right),$
(18)

where ${\mu }_{\mathit{\text{ia}}}^{\mathit{\text{max}}}$ represents the specific growth rate and bia is the temperature sensitivity coefficient for the ice algal growth. Following Lavoie et al. (2005), ${\mu }_{\mathit{\text{ia}}}^{\mathit{\text{max}}}$ and bia are set to 0.85 d–1 and 0.0633°C–1. The nitrogen (Lnit) and silicon (Lsil) limitation functions are based on the Monod formulation (Monod, 1949):

${L}_{\mathit{\text{nit}}}=\frac{{\left[N{O}_{3}\right]}_{\mathit{\text{ia}}}+{\left[N{H}_{4}\right]}_{\mathit{\text{ia}}}}{{k}_{n}+{\left[N{O}_{3}\right]}_{\mathit{\text{ia}}}+{\left[N{H}_{4}\right]}_{\mathit{\text{ia}}}}$
(19)
${L}_{\mathit{\text{sil}}}=\frac{{\left[Si\right]}_{\mathit{\text{ia}}}}{{k}_{s}+{\left[Si\right]}_{\mathit{\text{ia}}}}$
(20)

where [NO3]ia, [NH4]ia, and [Si]ia respectively represent the concentrations of nitrate (mmol N m–3), ammonium (mmol N m–3), and silicate (mmol Si m–3) in the bottom 3 cm of the ice skeletal layer. kn and ks are the half-saturation constants for nitrogen and silicon uptake, which are set to 1 mmol N m–3 and 4 mmol Si m–3 following Lavoie et al. (2005), respectively. The light limitation function (Llig) is formulated following Lavoie et al. (2005):

${L}_{\mathit{\text{lig}}}=\mathrm{tanh}\left(\frac{{\alpha }^{\mathit{\text{ia}}}}{{P}_{m}^{\mathit{\text{ia}}}}{\left[\mathit{\text{PAR}}\right]}_{\mathit{\text{ia}}}\right)$
(21)

where αia is the photosynthetic efficiency (mg C (mg Chl a)–1 h–1 (µmol photons m–2 s–1)–1), ${P}_{m}^{\mathit{\text{ia}}}$ is the maximum photosynthetic rate (mg C mg Chl a–1 h–1), and [PAR]ia is the PAR (in µmol photons m–2 s–1) penetrating the ice algae skeletal layer, which is computed by the sea ice thermodynamic model. Lavoie et al. (2005) used various values (0.1–0.45) for the ratio of αia and ${P}_{m}^{\mathit{\text{ia}}}$ depending on snow cover. We set the ratio of these two photosynthetic parameters $\left({\alpha }^{\mathit{\text{ia}}}/{P}_{m}^{\mathit{\text{ia}}}\right)$ to 0.44, based on the results of our sensitivity studies.

The ice growth/melt limitation function (Lice) is formulated similarly to Lavoie et al. (2005):

${L}_{\mathit{\text{ice}}}=\mathrm{max}\left(1-\frac{\left[{R}_{\mathit{\text{ice}}}\right]}{{C}_{m}},0\right)$
(22)

where [Rice] is the ice growth/melt rate at the ice base (m s–1) computed by the physical model (see Flato and Brown (1996) for details) and Cm is the critical ice growth/melt rate, which is set to 0.015 m d–1 following Lavoie et al. (2005).

The second term in Equation 17 represents the loss due to mortality. The mortality rate of ice algae (M) includes both linear and quadratic dependence on the ice algal biomass:

$M={m}_{\mathit{\text{lia}}}\mathrm{exp}\left({b}_{\mathit{\text{ia}}}{\left[T\right]}_{\mathit{\text{ia}}}\right)+{m}_{\mathit{\text{qia}}}\left[\mathit{\text{IA}}\right]$
(23)

where mlia is the rate constant for the temperature-dependent linear mortality and mqia is the rate constant for the quadratic mortality. We set mlia and mqia to 0.03 d–1 and 0.00015 d–1, based on the results of our sensitivity studies.

The last term in Equation 17 represents the loss due to flushing. There are four sources of meltwater that contributes to flushing: ice meltwater originating from the 1) top, 2) interior, and 3) bottom of the sea ice; and 4) snow meltwater that drains out of melt ponds. The sum of the ice melt rates (Rice) and the areal fraction of melt ponds (fpond) are computed by the physical model (see Flato and Brown (1996) and Abraham et al. (2015) for details). rpond represents the drainage rate of meltwater in the melt ponds, which is set to 0.0175 m d–1 following Taylor and Feltham (2004). Finally, zia is the thickness of the ice skeletal layer, which is set to 0.03 m in this study.

The equations for nitrate ([NO3]ia) and ammonium ([NH4]ia) in the ice skeletal layer (mmol N m–3) are defined as:

$\begin{array}{ll}\frac{\partial }{\partial t}{\left[N{O}_{3}\right]}_{ia}=\hfill & -{\mu }_{ia}\left[IA\right]{p}_{no3}^{ia}\frac{{\left[N{O}_{3}\right]}_{ia}}{{\left[N{O}_{3}\right]}_{ia}+{\left[N{H}_{4}\right]}_{ia}}\hfill \\ \hfill & +{r}_{nit}{\left[N{H}_{4}\right]}_{ia}-\frac{{f}_{pond}{r}_{pond}}{{z}_{ia}}{\left[N{O}_{3}\right]}_{ia}\hfill \\ \hfill & +\frac{D}{{h}_{\nu }}\left(\frac{\left[N{O}_{3}\right]-{\left[N{O}_{3}\right]}_{ia}}{{z}_{ia}}\right)\hfill \end{array}$
(24)
$\begin{array}{ll}\frac{\partial }{\partial t}{\left[N{H}_{4}\right]}_{ia}=\hfill & -{\mu }_{ia}\left[IA\right]\left(1-{p}_{no3}^{ia}\right)\frac{{\left[N{O}_{3}\right]}_{ia}}{{\left[N{O}_{3}\right]}_{ia}+{\left[N{H}_{4}\right]}_{ia}}\hfill \\ \hfill & -{r}_{nit}{\left[N{H}_{4}\right]}_{ia}-\frac{{f}_{pond}{r}_{pond}}{{z}_{ia}}{\left[N{H}_{4}\right]}_{ia}\hfill \\ \hfill & +\frac{D}{{h}_{\nu }}\left(\frac{\left[N{H}_{4}\right]-{\left[N{H}_{4}\right]}_{ia}}{{z}_{ia}}\right)\hfill \\ \hfill & +{f}_{lia}{m}_{lia}exp\left({b}_{ia}{\left[T\right]}_{ia}\right)\left[IA\right]\hfill \end{array}$
(25)

The first term in Equation 25 and Equation 26 represents the loss due to uptake by ice algae. To discriminate the uptake of nitrate from that of ammonium, the relative preference index (${p}_{\mathit{\text{no}}3}^{\mathit{\text{ia}}}$) based on Denman (2003) is defined as:

$\begin{array}{l}{p}_{\mathit{\text{no}}3}^{\mathit{\text{ia}}}=\frac{{\nu }_{n}}{{\nu }_{n}+{\left[N{H}_{4}\right]}_{\mathit{\text{ia}}}},\hfill \end{array}$
(26)

where νn represents the half-saturation constant for preferential uptake of nitrate, which is set to 0.2 mmol N m–3 following Denman (2003). The second term in Equation 25 and Equation 26 represents the nitrification, which is a source for nitrate and a sink for ammonium. In this study, a first order reaction is assumed for this process in the sea ice with a rate constant of 0.01 d–1. The third term in Equation 25 and Equation 26 represents the loss due to flushing, which is formulated in a similar manner as the flushing of ice algae. The fourth term in Equation 25 and Equation 26 represents the diffusive exchange at the ice-ocean interface, where [NO3] and [NH4] are the concentrations (mmol N m–3) of nitrate and ammonium in the uppermost layer of the water column. D is the molecular diffusion coefficient for dissolved nutrients at the ice-water interface, which is set to 4.7 × 10–8 m2 s–1 (Rebreanu et al., 2008). The exchange rate at the ice-water interface depends on the molecular sublayer thickness (hν) defined as:

${h}_{\nu }=\frac{\nu }{{u}_{\tau }},$
(27)

where ν is the kinematic viscosity of seawater, which is set to 1.85 × 10–6 m2 s–1 (Lavoie et al., 2005). uτ represents the friction velocity, which is defined as:

${u}_{\tau }={\left(\frac{|\tau |}{{\rho }_{o}}\right)}^{1/2},$
(28)

where τ is the ice-ocean stress and ρ0 is the seawater density. Following Goosse and Fichefet (1999), τ is defined as:

$\tau ={\rho }_{o}{C}_{\mathit{\text{io}}}\text{|}{u}_{i}-{u}_{o}\text{|}\left({u}_{i}-{u}_{o}\right),$
(29)

where Cio is the drag coefficient at the ice-ocean interface, which is set to 0.0054 (Shirasawa and Ingram, 1997). ui and uo are the magnitudes of horizontal velocity fields of ice and surface water, respectively. Substituting Equations 29 and 30 into Equation 28 gives:

${h}_{\nu }=\frac{\nu }{\text{|}{u}_{i}-{u}_{o}\text{|}}{C}_{\mathit{\text{io}}}^{-\frac{1}{2}},$
(30)

where uo is computed by the physical ocean model, and we assume ui = 0 for landfast ice considered in this study.

The fifth term in Equation 26 represents the remineralization of dead ice algal cells. A direct transfer from ice algae to ammonium is implemented here due of lack of representation of detritus pool in the sea ice. It is assumed that only a fraction (flia) of linear mortality enters the ammonium pool, while the remaining fraction is released into the water column and enters the large detritus pool. In this study, flia is set to 0.3.

The equation for silicate ([Si]ia) in the ice skeletal layer (mmol Si m–3) is given by:

$\begin{array}{l}\frac{\partial }{\partial t}{\left[Si\right]}_{\mathit{\text{ia}}}=-{q}_{\mathit{\text{si}}2n}{\mu }_{\mathit{\text{ia}}}\left[\mathit{\text{IA}}\right]-\frac{{f}_{\mathit{\text{pond}}}{r}_{\mathit{\text{pond}}}}{{z}_{\mathit{\text{ia}}}}{\left[Si\right]}_{\mathit{\text{ia}}}\\ +\frac{D}{{h}_{\nu }}\left(\frac{\left[Si\right]-{\left[Si\right]}_{\mathit{\text{ia}}}}{{z}_{\mathit{\text{ia}}}}\right)\end{array}$
(31)

where the first term represents the uptake by ice algae, the second term represents the flushing, and the third term represents the diffusive mixing at the ice-ocean interface. The uptake rate is converted from nitrogen to silicon units by assuming a fixed nitrogen-to-silicon intracellular ratio (qsi2n) of 1.7 mol Si:mol N based on Mundy et al. (2014). [Si] represents the concentration of silicate in the uppermost layer of the water column (mmol Si m–3).

#### 6.2.2 Oceanic component

The ocean biogeochemical model consists of ten state variables: small and large classes of phytoplankton, zooplankton, and detritus, biogenic silica, nitrate, ammonium, and silicate. The equations for small and large phytoplankton are:

$\begin{array}{l}\frac{\partial }{\partial t}\left[P1\right]={\mu }_{p1}\left[P1\right]-{\Gamma }_{z1}\frac{\left[P1\right]}{{F}_{z1}}\left[Z1\right]-{M}_{\mathit{\text{lp}}1}\left[P1\right]\\ +\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[P1\right]}{\partial z}\right)\end{array}$
(32)
$\begin{array}{l}\frac{\partial }{\partial t}\left[P2\right]={\mu }_{p2}\left[P2\right]-{\Gamma }_{z2}\frac{\left[P2\right]}{{F}_{z2}}\left[Z2\right]-{M}_{\mathit{\text{lp}}2}\left[P2\right]\\ +\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[P2\right]}{\partial z}\right)-{m}_{\mathit{\text{qp}}2}{\left[P2\right]}^{2}\\ +{f}_{\mathit{\text{seed}}}\left(\left[{R}_{\mathit{\text{ice}}}\right]+{r}_{\mathit{\text{pond}}}\right)\left[\mathit{\text{IA}}\right]\frac{{z}_{\mathit{\text{ia}}}}{{z}_{\mathit{\text{oc}}}}{\delta }_{z,{z}_{1}}\end{array}$
(33)

where [P1] and [P2] respectively represent the biomass of small and large phytoplankton (mmol N m–3). The first term in Equations 33 and 34 represent the growth due to primary production. The growth rates of phytoplankton (µp1 and µp2) depend on the ambient temperature, nutrient and light conditions:

${\mu }_{p1}={\mu }_{p1}^{\mathit{\text{max}}}{Q}_{10p}^{\left(\left[T\right]-10\right)/10}\mathrm{min}\left({L}_{\mathit{\text{nit}}}^{\mathit{\text{phy}}},{L}_{\mathit{\text{lig}}}^{p1}\right)$
(34)
${\mu }_{p2}={\mu }_{p2}^{\mathit{\text{max}}}{Q}_{10p}^{\left(\left[T\right]-10\right)/10}\mathrm{min}\left({L}_{\mathit{\text{nit}}}^{\mathit{\text{phy}}},{L}_{\mathit{\text{lig}}}^{p2},{L}_{\mathit{\text{sil}}}^{p2}\right)$
(35)

where ${\mu }_{p1}^{\mathit{\text{max}}}$ and ${\mu }_{p2}^{\mathit{\text{max}}}$ represent the maximum specific growth rates, which are set to 0.5 and 2.0 d–1, respectively (Steiner et al., 2006). The temperature dependence is based on the Q10 formulation with a Q10 factor (Q10p) of 1.55 (Suzuki and Takahashi, 1995). [T] is the ambient seawater temperature (°C), which is simulated by the physical model. The nutrient limitation functions are based on the Monod function (Monod, 1949):

${L}_{\mathit{\text{nit}}}^{\mathit{\text{phy}}}=\frac{\left[N{O}_{3}\right]+\left[N{H}_{4}\right]}{{k}_{n}+\left(\left[N{O}_{3}\right]+\left[N{H}_{4}\right]\right)}$
(36)
${L}_{\mathit{\text{sil}}}^{p2}=\frac{\left[Si\right]}{{k}_{\mathit{\text{si}}}+\left[Si\right]}$
(37)

where the same half-saturation constants are used between the sea ice and ocean biogeochemical modules. The light limitation functions are defined as:

${L}_{\mathit{\text{lig}}}^{p1}=1-\mathrm{exp}\left(-{\alpha }^{p1}/{P}_{m}^{p1}\left[\mathit{\text{PAR}}\right]\right)$
(38)
${L}_{\mathit{\text{lig}}}^{p2}=1-\mathrm{exp}\left(-{\alpha }^{p2}/{P}_{m}^{p2}\left[\mathit{\text{PAR}}\right]\right)$
(39)

where αp1 and αp2 represent the photosynthetic efficiency and ${P}_{m}^{p1}$ and ${P}_{m}^{p2}$ represent the maximum photosynthetic rate. In this study, the ratios of photosynthetic parameters (i.e., ${\alpha }^{p\text{1}}/{P}_{m}^{p1}$ and ${\alpha }^{p\text{2}}/{P}_{m}^{p2}$) are set to 0.07 for both phytoplankton groups. [PAR] represents the photosynthetically active radiation in the water column (µmol photons m–2 s–1) which is defined in the previous section. The second term in Equations 33 and 34 represent the loss due to grazing by zooplankton. The rates of small and large zooplankton (Γz1 and Γz2) are defined as:

${\Gamma }_{z1}={\gamma }_{z1}{Q}_{10z}^{\left(\left[T\right]-10\right)/10}\frac{{F}_{z1}^{2}}{{k}_{z1}^{2}+{F}_{z1}^{2}}$
(40)
${\Gamma }_{z2}={\gamma }_{z2}{Q}_{10z}^{\left(\left[T\right]-10\right)/10}\frac{{F}_{z2}^{2}}{{k}_{z2}^{2}+{F}_{z2}^{2}}$
(41)

where γz1 and γz2 represent the maximum specific grazing rates, which are set to 1 d–1 for both zooplankton groups in this study. The Q10 factor for zooplankton (Q10z) is set to 2.14 (Buitenhuis et al., 2006). kz1 and kz2 represent the half-saturation constants for grazing, which are set to 1 mmol N m–3 for both zooplankton groups in this study. Fz1 and Fz2 respectively represent the available food sources for small and large zooplankton, which are defined as:

${F}_{z1}=\left[P1\right]+{f}_{d1}\left[D1\right]$
(42)
${F}_{z2}=\left[P2\right]+{f}_{d2}\left[D2\right]+{f}_{z1}\left[Z1\right]$
(43)

where fd1 represents the grazing preference of small zooplankton on small detitus. fd2 and fz2 represent the grazing preferences of large zooplankton on small and large detritus, respectively. fd1, fd2, and fz1 are all set to 0.5 following Steiner et al. (2006).

The third terms in Equations 33 and 34 represent the linear mortality of small and large phytoplankton. Assuming a temperature depnedence following the same Q10 formulation, the rates of linear mortality are given by:

${M}_{\mathit{\text{lp}}1}={m}_{\mathit{\text{lp}}1}{Q}_{10p}^{\left(\left[T\right]-10\right)/10}$
(44)
${M}_{\mathit{\text{lp}}2}={m}_{\mathit{\text{lp}}2}{Q}_{10p}^{\left(\left[T\right]-10\right)/10}$
(45)

where mlp1 and mlp2 represent the rate constants for linear mortality of small and large phytoplankton, respectively. In this study, both of these rates are set to 0.03 d–1.

The fourth terms in Equations 33 and 34 represent the vertical mixing of small and large phytoplankton at a given model layer with those in the adjacent (upper and/or lower) layers. Kz represents the vertical eddy diffusivity coefficient, which is calculated by the physical ocean model.

The fifth term in Equation 34 represents the quadratic mortality of large phytoplankton, which depends on a rate constant, mqp2. In this study, mqp2 is set to 0.05 d–1 (mmol N m–3)–1. The sixth term in Equation 34 represents the source for large phytoplankton due to seeding of ice algae. Following Lavoie et al. (2009), it is assumed that a certain fraction (fseed) of ice algae flushed into the water column enters the large phytoplankton pool, while the remaining fraction (1 – fseed) enters the large detritus pool. In this study, fseed is set to 0.1, assuming that 10% of ice algae flushing results in large phytoplankton seeding. hoc is the thickness of the uppermost layer of the ocean model (i.e., 1 m in this study), and ${\delta }_{z,{z}_{1}}$ is the Kronecker’s delta which equals 1 in the uppermost layer of the ocean (z1), while it is set to 0 elsewhere:

${\delta }_{z,{z}_{1}}=\left\{\begin{array}{ll}0\hfill & \text{if} z\ne {z}_{1}\hfill \\ 1\hfill & \text{if} z={z}_{1}\hfill \end{array}$
(46)

The equations for small and large zooplankton are defined as:

$\begin{array}{l}\frac{\partial }{\partial t}\left[Z1\right]={a}_{1}{\Gamma }_{z1}\left[Z1\right]-{M}_{\mathit{\text{lz}}1}\left[Z1\right]-{\Gamma }_{z2}\frac{{f}_{z1}\left[Z1\right]}{{F}_{z2}}\left[Z2\right]\\ +\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[Z1\right]}{\partial z}\right)\end{array}$
(47)
$\begin{array}{l}\frac{\partial }{\partial t}\left[Z2\right]={a}_{2}{\Gamma }_{z2}\left[Z2\right]-{M}_{\mathit{\text{lz}}2}\left[Z2\right]-{m}_{\mathit{\text{qz}}2}{\left[Z2\right]}^{2}\\ +\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[Z2\right]}{\partial z}\right)\end{array}$
(48)

where [Z1] and [Z2] respectively represent the biomass of small and zooplankton (mmol N m–3).

The first terms in Equations 48 and 49 represent the sources due to ingestion, which are proportional to the zooplankton grazing. a1 and a2 respectively represent the assimilated fractions of grazing by small and large zooplankton, both of which are set to 0.7 following Lavoie et al. (2009).

The second terms in Equations 48 and 49 represent the sinks due to linear mortality. The rates of zooplankton linear mortality (Mlz1 and Mlz2) are parameterized in the same way as those of phytoplankton:

${M}_{\mathit{\text{lz}}1}={m}_{\mathit{\text{lz}}1}{Q}_{10z}^{\left(\left[T\right]-10\right)/10}$
(49)
${M}_{\mathit{\text{lz}}2}={m}_{\mathit{\text{lz}}2}{Q}_{10z}^{\left(\left[T\right]-10\right)/10}$
(50)

where the rate constants for linear mortality of small and large zooplankton (mlz1 and mlz2) are both set to 0.03 d–1 (Lavoie et al., 2009).

The third term in Equation 48 represents the loss of small zooplankton due to (carnivorous) grazing by large zooplankton. The third term in Equation 49 represents the loss of large zooplankton due to quadratic mortality which is parameterized in the same way as the quadratic mortality of large phytoplankton with a rate constant (mqz2) of 0.1 d–1 (mmol N m–3)–1 in this study.

The fourth terms in Equations 48 and 49 respectively represent the vertical mixing of small and large zooplankton. The equations for small and large detritus, as well as biogenic silica are given by:

$\begin{array}{l}\frac{\partial }{\partial t}\left[D1\right]=\left(1-{a}_{1}\right){\Gamma }_{z1}\left[Z1\right]+{M}_{\mathit{\text{lp}}1}\left[P1\right]\\ -{\Gamma }_{z1}\frac{{f}_{d1}\left[D1\right]}{{F}_{z1}}\left[Z1\right]-{R}_{d1}\left[D1\right]\\ -{w}_{d1}\frac{\partial \left[D1\right]}{\partial z}+\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[D1\right]}{\partial z}\right)\end{array}$
(51)
$\begin{array}{l}\frac{\partial }{\partial t}\left[D2\right]=\left(1-{a}_{2}\right){\Gamma }_{z2}\left[Z2\right]+{M}_{\mathit{\text{lp}}2}\left[{P}_{2}\right]\\ -{\Gamma }_{z2}\frac{{f}_{d2}\left[D2\right]}{{F}_{z2}}\left[Z2\right]-{R}_{d2}\left[D2\right]\\ -{w}_{d2}\frac{\partial \left[D2\right]}{\partial z}+\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[D2\right]}{\partial z}\right)\\ +{m}_{\mathit{\text{qp}}2}{\left[{P}_{2}\right]}^{2}\\ +\left(1-{f}_{\mathit{\text{seed}}}\right)\left(\left[{R}_{\mathit{\text{ice}}}\right]+{r}_{\mathit{\text{pond}}}\right)\left[\mathit{\text{IA}}\right]\frac{{z}_{\mathit{\text{ia}}}}{{z}_{\mathit{\text{oc}}}}{\delta }_{z,{z}_{1}}\\ +\left(1-{f}_{\mathit{\text{lia}}}\right){m}_{\mathit{\text{lia}}}\mathit{\text{exp}}\left({b}_{\mathit{\text{ia}}}{\left[T\right]}_{\mathit{\text{ia}}}\right)\left[\mathit{\text{IA}}\right]\frac{{z}_{\mathit{\text{ia}}}}{{z}_{\mathit{\text{oc}}}}{\delta }_{z,{z}_{1}}\\ +{m}_{\mathit{\text{qia}}}\left[\mathit{\text{IA}}\right]\frac{{z}_{\mathit{\text{ia}}}}{{z}_{\mathit{\text{oc}}}}{\delta }_{z,z0}\end{array}$
(52)
$\begin{array}{l}\frac{\partial }{\partial t}\left[BSi\right]={q}_{\mathit{\text{si}}2n}\left(\left(1-{a}_{2}\right){\Gamma }_{z2}\left[Z2\right]+{M}_{\mathit{\text{lp}}2}\left[{P}_{2}\right]-{\Gamma }_{z2}\frac{{f}_{d2}\left[D2\right]}{{F}_{z2}}\left[Z2\right]\right)\\ -{R}_{\mathit{\text{bsi}}}\left[BSi\right]-{w}_{d2}\frac{\partial \left[BSi\right]}{\partial z}+\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[BSi\right]}{\partial z}\right)\\ +{q}_{\mathit{\text{si}}2n}{m}_{\mathit{\text{qp}}2}{\left[{P}_{2}\right]}^{2}\\ +{q}_{\mathit{\text{si}}2n}\left(1-{f}_{\mathit{\text{seed}}}\right)\left(\left[{R}_{\mathit{\text{ice}}}\right]+{r}_{\mathit{\text{pond}}}\right)\left[\mathit{\text{IA}}\right]\frac{{z}_{\mathit{\text{ia}}}}{{z}_{\mathit{\text{oc}}}}{\delta }_{z,{z}_{1}}\\ +{q}_{\mathit{\text{si}}2n}\left(1-{f}_{\mathit{\text{lia}}}\right){m}_{\mathit{\text{lia}}}\mathit{\text{exp}}\left({b}_{\mathit{\text{ia}}}{\left[T\right]}_{\mathit{\text{ia}}}\right)\left[\mathit{\text{IA}}\right]\frac{{z}_{\mathit{\text{ia}}}}{{z}_{\mathit{\text{oc}}}}{\delta }_{z,{z}_{1}}\\ +{q}_{\mathit{\text{si}}2n}{m}_{\mathit{\text{qia}}}\left[\mathit{\text{IA}}\right]\frac{{z}_{\mathit{\text{ia}}}}{{z}_{\mathit{\text{oc}}}}{\delta }_{z,z0}\end{array}$
(53)

In each of Equations 52, 53, and 54, the first term represents the unassimilated fraction of grazing that enter the detrital pool, the second term represents the source for detritus due to phytoplankton linear mortality, the third term represents the loss of detritus due to grazing, and the fourth term represents the remineralization of detrital materials. The rates of remineralization of small and large detritus and biogenic silica (Rd1, Rd2, Rbsi) depend on seawater temperature following the Q10 formulation:

${R}_{d1}={r}_{d1}{Q}_{10b}^{\left(\left[T\right]-10\right)/10}$
(54)
${R}_{d2}={r}_{d2}{Q}_{10b}^{\left(\left[T\right]-10\right)/10}$
(55)
${R}_{\mathit{\text{bsi}}}={r}_{\mathit{\text{bsi}}}{Q}_{10b}^{\left(\left[T\right]-10\right)/10}$
(56)

where rd1, rd2, and rbsi represent the rate constants for remineralization, which are respectively set to 0.03 d–1 (this study), 0.3 d–1 (Lavoie et al., 2009), and 0.01 d–1 (Steiner et al., 2006).

The fifth terms in Equations 52, 53, and 54 represent the loss due to sinking. It is assumed that detritus sink at a fixed rate; a rate of 1 m d–1 is given to small detritus, while a faster sinking rate of 50 m d–1 is prescribed for large detritus and biogenic silica. Both of these rates are taken from Lavoie et al. (2009). Detritus reaching the seafloor is accumulated in the deepest layer of the water column. The sixth terms in Equations 52, 53, and 54 represent the vertical mixing of respective detrital materials.

The seventh terms in Equations 52 and 53 represent the sources for small and large detritus due to quadratic mortality of large phytoplankton. The last three terms in these equations represent the influxes of detrital materials from the ice skeletal layer due to flushing (the ninth terms), linear mortality (the tenth terms), and quadratic mortality (the eleventh terms) of ice algae. Note that the equation for biogenic silica is identical to the equation for large detritus except that the two state variables have different remineralization rates and that those terms expressed in nitrogen units are converted to silicon units by assuming a fixed intracellular nitrogen-to-silicon ratio of 1.7 mol N:mol Si (Mundy et al., 2014).

The equations for dissolved nutrients in the water column are defined as:

$\begin{array}{l}\frac{\partial }{\partial t}\left[N{O}_{3}\right]=-\left({\mu }_{p1}\left[P1\right]+{\mu }_{p2}\left[{P}_{2}\right]\right){p}_{\mathit{\text{no}}3}\frac{\left[N{O}_{3}\right]}{\left[N{O}_{3}\right]+\left[N{H}_{4}\right]}\\ +\eta \left[N{H}_{4}\right]-\frac{D}{{h}_{\nu }}\left(\left[N{O}_{3}\right]-{\left[N{O}_{3}\right]}_{\mathit{\text{ia}}}\right){\delta }_{z,{z}_{1}}\\ +\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[N{O}_{3}\right]}{\partial z}\right)\end{array}$
(57)
$\begin{array}{l}\frac{\partial }{\partial t}\left[N{H}_{4}\right]=-\left({\mu }_{p1}\left[P1\right]+{\mu }_{p2}\left[{P}_{2}\right]\right)\\ \left(1-{p}_{\mathit{\text{no}}3}\frac{\left[N{O}_{3}\right]}{\left[N{O}_{3}\right]+\left[N{H}_{4}\right]}\right)-\eta \left[N{H}_{4}\right]\\ +{M}_{\mathit{\text{lz}}1}\left[Z1\right]+{M}_{\mathit{\text{lz}}2}\left[Z2\right]+{R}_{d1}\left[D1\right]\\ +{R}_{d2}\left[D2\right]-\frac{D}{{h}_{\nu }}\left(\left[N{H}_{4}\right]-{\left[N{H}_{4}\right]}_{\mathit{\text{ia}}}\right){\delta }_{z,{z}_{1}}\\ +\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[N{H}_{4}\right]}{\partial z}\right)\end{array}$
(58)
$\begin{array}{l}\frac{\partial }{\partial t}\left[Si\right]=-{q}_{\mathit{\text{si}}2n}{\mu }_{p2}\left[{P}_{2}\right]+{R}_{\mathit{\text{bsi}}}\left[BSi\right]\\ -\frac{D}{{h}_{\nu }}\left(\left[Si\right]-{\left[Si\right]}_{\mathit{\text{ia}}}\right){\delta }_{z,{z}_{1}}+\frac{\partial }{\partial z}\left({K}_{z}\frac{\partial \left[Si\right]}{\partial z}\right)\end{array}$
(59)

where [NO3], [NH4], and [Si] represent the concentrations of nitrate, ammonium, and silicate, respectively.

The first terms in Equations 58 and 59 represent the uptake by both small and large phytoplankton. The first term in Equation 60 represents the uptake by large phytoplankton only, assuming an insignificant silicon uptake by small phytoplankton. Similarly to the nitrogen uptake by ice algae, the relative preference index for nitrate uptake is formulated as a function of ammonium concentration following Denman (2003):

${p}_{\mathit{\text{no}}3}=\frac{{v}_{n}}{{v}_{n}+\left[N{H}_{4}\right]}$
(60)

The second terms in Equations 58 and 59 represent the nitrification, which is a source for nitrate and a sink for ammonium. The rate of nitrification depends on temperature (Kawamiya et al., 1995) and is assumed to be photo-inhibited (Aumont et al., 2015):

$\begin{array}{l}\eta =\frac{{r}_{\mathit{\text{nit}}}\mathrm{exp}\left({b}_{\mathit{\text{nit}}}\left[T\right]\right)}{1+\left[\mathit{\text{PAR}}\right]}\hfill \end{array}$
(61)

where the rate constant (rnit) and the temperature sensitivity coefficient (bnit) for nitrification are set to 0.03 d–1 and 0.0693 (°C)–1 following Kawamiya et al. (1995).

The third and the fourth terms in Equation 59 respectively represent the sources for ammonium due to linear mortality of small and large zooplankton. The fifth and the six terms in the same equation represent the sources due to remineralization of small and large detritus, respectively.

The second term in Equation 60 represents the source for silicate due to remineralization of biogenic silica. The second last terms in Equations 58, 59, and 60 represent the diffusive exchange of nutrients between the uppermost layer of the water column and the ice skeletal layer. The last terms in these equations represent the vertical mixing of dissolved nutrients between the model layers in the water column.

Lastly, for each of the living organisms (IA, P1, P2, Z1, and Z2), all the sink terms in the equation are neglected (i.e., set to zero) when its biomass is below its initial concentration in order to maintain its overwintering level.

## Acknowledgements

This is a contribution to the Scientific Committee on Oceanic Research (SCOR) Working group on Biogeochemical Exchange Processes at Sea Ice Interfaces (BEPSII), the Network on Climate and Aerosols: Addressing Key Uncertainties in Remote Canadian Environments (NETCARE), ArcticNet, Québec-Océan, Arctic Science Partnership (ASP), and the Canada Excellence Research Chair unit at the Centre for Earth Observation Science. The atmospheric forcing data were provided by Environment Canada.

We would also like to thank the reviewers for their feedback and advice.

## Funding information

EM and HH were funded via the NETCARE and ArcticNet projects. NS and DL were supported by Fisheries and Oceans Canada. VG, MG, AM, and CJM were supported by funding from Natural Sciences and Engineering Research Council of Canada (NSERC), Fonds de recherche du Québec Nature et Technologies (FRQNT), Canada Economic Development and Polar Continental Shelf Program (PCSP) of Natural Resources Canada.

## Competing interests

The authors have no competing interests to declare.

## Author contributions

• Contributed to conception and design: EM, HH, NS, AM
• Contributed to model development: EM, HH, NS, AM, DL
• Contributed to model experiments: EM, HH
• Contributed to acquisition of restoration data: XH
• Contributed to acquisition of field data: MB, MGa, VG, MGo, CJM
• Contributed to analysis and interpretation of data: EM, HH, NS, AM
• Drafted and/or revised the article: All
• Approved the submitted version for publication: All

## References

1. Abraham C, Steiner N, Monahan A and Michel C 2015. Effects of subgrid-scale snow thickness variability on radiative transfer in sea ice. Journal of Geophysical Research: Oceans 120(8): 5597–5614, DOI: http://dx.doi.org/10.1002/2015JC010741 ISSN 2169-9291.

2. Arrigo KR 2014. Sea ice ecosystems. Annual Review of Marine Science 6(1): 439–467, DOI: http://dx.doi.org/10.1146/annurev-marine-010213-135103

3. Arrigo KR and Sullivan CW 1994. A high resolution bio-optical model of microalgal growth: Tests using sea-ice algal community time-series data. Limnology and Oceanography 39(3): 609–631, DOI: http://dx.doi.org/10.4319/lo.1994.39.3.0609 ISSN 1939-5590.

4. Aumack CF and Juhl AR 2015. Light and nutrient effects on the settling characteristics of the sea ice diatom. Nitzschia frigida. Limnology and Oceanography 60(3): 765–776, DOI: http://dx.doi.org/10.1002/lno.10054 ISSN 1939-5590.

5. Aumont O, Ethé C, Tagliabue A, Bopp L and Gehlen M 2015. PISCES-v2: an ocean biogeochemical model for carbon and ecosystem studies. Geoscientific Model Development 8(8): 2465–2513, DOI: http://dx.doi.org/10.5194/gmd-8-2465-2015 ISSN 1991-9603.

6. Barber DG Hop H Mundy CJ Else B Dmitrenko IA et al. 2015. Selected physical, biological and biogeochemical implications of a rapidly changing Arctic Marginal Ice Zone. Progress in Oceanography 139: 122–150, DOI: http://dx.doi.org/10.1016/j.pocean.2015.09.003 ISSN 0079-6611.

7. Brown KA Miller LA Mundy CJ Papakyriakou T Francois R et al. 2015–1. Inorganic carbon system dynamics in landfast Arctic sea ice during the early-melt period: Spring Sea Ice Inorganic Carbon System. Journal of Geophysical Research: Oceans 120(5): 3542–3566, DOI: http://dx.doi.org/10.1002/2014JC010620 ISSN 2169-9275.

8. Brown ZW Lowry KE Palmer MA van Dijken GL Mills MM et al. 2015–2. Characterizing the subsurface chlorophyll a maximum in the Chukchi Sea and Canada Basin. Deep Sea Research Part II: Topical Studies in Oceanography 118(Part A): 88–104, DOI: http://dx.doi.org/10.1016/j.dsr2.2015.02.010 ISSN 0967-0645.

9. Bruggeman J and Bolding K 2014. A general framework for aquatic biogeochemical models. Environmental Modelling Software, DOI: http://dx.doi.org/10.1016/j.envsoft.2014.04.002 ISSN 1364-8152.

10. Buitenhuis E Le Quéré C Aumont O Beaugrand G Bunker A et al. 2006. Biogeochemical fluxes through mesozooplankton. Global Biogeochemical Cycles 20(2): GB2003.DOI: http://dx.doi.org/10.1029/2005GB002511 ISSN 1944-9224.

11. Burchard H Bolding K Kühn W Meister A Neumann T et al. 2006. Description of a flexible and extendable physical biogeochemical model system for the water column. Journal of Marine Systems 61(3–4): 180–211, DOI: http://dx.doi.org/10.1016/j.jmarsys.2005.04.011 ISSN 0924-7963.

12. Carmack E and McLaughlin F 2011. Towards recognition of physical and geochemical change in Subarctic and Arctic Seas. Progress in Oceanography 90(1–4): 90–104, DOI: http://dx.doi.org/10.1016/j.pocean.2011.02.007 ISSN 0079-6611.

13. Carmack EC, McLaughlin FA, Vagle S, Melling H and Williams WJ 2010. Structures and property distributions in the three oceans surrounding Canada in 2007: A basis for a long-term ocean climate monitoring strategy. Atmosphere-Ocean 48(4): 211–224, DOI: http://dx.doi.org/10.3137/OC324.2010 ISSN 0705-5900.

14. Cota GF Prinsenberg SJ Bennett EB Loder JW Lewis MR et al. 1987. Nutrient fluxes during extended blooms of Arctic ice algae. Journal of Geophysical Research: Oceans 92(C2): 1951–1962, DOI: http://dx.doi.org/10.1029/JC092iC02p01951 ISSN 2156-2202.

15. Deal C Jin M Elliott S Hunke E Maltrud M et al. 2011. Large-scale modeling of primary production and ice algal biomass within arctic sea ice in 1992. Journal of Geophysical Research 116(C7)DOI: http://dx.doi.org/10.1029/2010JC006409 ISSN 0148-0227.

16. Denman KL 2003. Modelling planktonic ecosystems, Modelling planktonic ecosystems: parameterizing complexity. Progress in Oceanography 57(3–4): 429–452, DOI: http://dx.doi.org/10.1016/S0079-6611(03)00109-5 ISSN 0079-6611.

17. Dukhovskoy DS Myers PG Platov G Timmermans ML Curry B et al. 2016. Greenland freshwater pathways in the sub-Arctic Seas from model experiments with passive tracers. Journal of Geophysical Research: Oceans 121(1): 877–907, DOI: http://dx.doi.org/10.1002/2015JC011290 ISSN 2169-9291.

18. Eppley RW 1972. Temperature and phytoplankton growth in the sea. Fish. Bull 70(4): 1063–1085.

19. Flato GM and Brown RD 1996. Variability and climate sensitivity of landfast Arctic sea ice. Journal of Geophysical Research: Oceans 101(C11): 25767–25777, DOI: http://dx.doi.org/10.1029/96JC02431 ISSN 2156-2202.

20. Frey KE, Perovich DK and Light B 2011. The spatial distribution of solar radiation under a melting Arctic sea ice cover. Geophysical Research Letters 38(22): L22501.DOI: http://dx.doi.org/10.1029/2011GL049421 ISSN 1944-8007.

21. Fripiat F Sigman DM Fawcett SE Rafter PA Weigand MA et al. 2014. New insights into sea ice nitrogen biogeochemical dynamics from the nitrogen isotopes. Global Biogeochemical Cycles 28(2): 2013.DOI: http://dx.doi.org/10.1002/2013GB004729 GB004729. ISSN 1944-9224.

22. Gale M 2014. On The Biophysical Factors That Control Under-Ice Phytoplankton Bloom Onset in the Central Canadian Archipelago. Ph.D. thesis,

23. Galindo V Levasseur M Mundy C Gosselin M Tremblay JE et al. 2014. Biological and physical processes influencing sea ice, under-ice algae, and dimethylsulfoniopropionate during spring in the Canadian Arctic Archipelago. Journal of Geophysical Research: Oceans, DOI: http://dx.doi.org/10.1002/2013JC009497 ISSN 2169-9291.

24. Geilfus NX Galley RJ Crabeck O Papakyriakou T Landy J et al. 2015. Inorganic carbon dynamics of melt-pond-covered first-year sea ice in the Canadian Arctic. Biogeosciences 12(6): 2047–2061, DOI: http://dx.doi.org/10.5194/bg-12-2047-2015 ISSN 1726-4189.

25. Goosse H and Fichefet T 1999. Importance of ice-ocean interactions for the global ocean circulation: A model study. Journal of Geophysical Research: Oceans 104(C10): 23337–23355, DOI: http://dx.doi.org/10.1029/1999JC900215 ISSN 2156-2202.

26. Gosselin M, Legendre L, Therriault JC and Demers S 1990. Light and nutrient limitation of sea-ice microalgae (Hudson Bay, Canadian Arctic). Journal of Phycology 26(2): 220–232, DOI: http://dx.doi.org/10.1111/j.0022-3646.1990.00220.x ISSN 1529-8817.

27. Gosselin M, Levasseur M, Wheeler PA, Horner RA and Booth BC 1997. New measurements of phytoplankton and ice algal production in the Arctic Ocean. Deep Sea Research Part II: Topical Studies in Oceanography 44(8): 1623–1644, DOI: http://dx.doi.org/10.1016/S0967-0645(97)00054-4

28. Grenfell TC and Maykut GA 1977. The optical properties of ice and snow in the Arctic Basin. Journal of Glaciology 18: 445–463. ISSN 0022-1430.

29. Harrison WG, Cota GF and Smith REH 1990. Nitrogen utilization in ice algal communities of Barrow Strait, Northwest Territories, Canada. Marine Ecology Progress Series 67: 275–283, DOI: http://dx.doi.org/10.3354/meps067275

30. Hayashida H, Steiner N, Monahan A, Galindo V, Lizotte M and Levasseur M. 2017. Implications of sea-ice biogeochemistry for oceanic production and emissions of dimethyl sulfide in the Arctic. Biogeosciences 14: 3129–3155, DOI: http://dx.doi.org/10.5194/bg-14-3129-2017

31. Jin M Deal C Lee SH Elliott S Hunke E et al. 2012. Investigation of Arctic sea ice and ocean primary production for the period 1992–2007 using a 3-D global ice–ocean ecosystem model. Deep Sea Research Part II: Topical Studies in Oceanography 81–84: 28–35, DOI: http://dx.doi.org/10.1016/j.dsr2.2011.06.003 ISSN 0967-0645.

32. Jin M Deal C Wang J Alexander V Gradinger R et al. 2007. Ice-associated phytoplankton blooms in the southeastern Bering Sea. Geophysical Research Letters 34(6): L06612.DOI: http://dx.doi.org/10.1029/2006GL028849 ISSN 1944-8007.

33. Jin M Deal CJ Wang J Shin KH Tanaka N et al. 2006. Controls of the landfast ice–ocean ecosystem offshore Barrow, Alaska. Annals of Glaciology 44(1): 63–72, DOI: http://dx.doi.org/10.3189/172756406781811709

34. Kawamiya M, Kishi MJ, Yamanaka Y and Suginohara N 1995. An ecological-physical coupled model applied to Station Papa. Journal of Oceanography 51(6): 635–664, DOI: http://dx.doi.org/10.1007/BF02235457

35. Kirst GO and Wiencke C 1995. Ecophysiology of polar algae. Journal of Phycology 31(2): 181–199, DOI: http://dx.doi.org/10.1111/j.0022-3646.1995.00181.x

36. Krembs C, Eicken H and Deming JW 2011. Exopolymer alteration of physical properties of sea ice and implications for ice habitability and biogeochemistry in a warmer Arctic. Proceedings of the National Academy of Sciences 108(9): 3653–3658, DOI: http://dx.doi.org/10.1073/pnas.1100701108 ISSN 0027-8424, 1091-6490.

37. Lavoie D, Denman K and Michel C 2005. Modeling ice algal growth and decline in a seasonally ice-covered region of the Arctic (Resolute Passage, Canadian Archipelago). Journal of Geophysical Research 110(C11)DOI: http://dx.doi.org/10.1029/2005JC002922 ISSN 0148-0227.

38. Lavoie D, Macdonald RW and Denman KL 2009. Primary productivity and export fluxes on the Canadian shelf of the Beaufort Sea, Primary productivity and export fluxes on the Canadian shelf of the Beaufort Sea: A modelling study. Journal of Marine Systems 75(1–2): 17–32, DOI: http://dx.doi.org/10.1016/j.jmarsys.2008.07.007 ISSN 0924-7963.

39. Legendre L Ackley SF Dieckmann GS Gulliksen B Horner R et al. 1992. Ecology of sea ice biota. Polar Biology 12(3–4): 429–444, DOI: http://dx.doi.org/10.1007/BF00243114 ISSN 0722-4060, 1432-2056.

40. Leu E Mundy CJ Assmy P Campbell K Gabrielsen TM et al. 2015. Arctic spring awakening – Steering principles behind the phenology of vernal ice algal blooms. Progress in Oceanography 139: 151–170, DOI: http://dx.doi.org/10.1016/j.pocean.2015.07.012 ISSN 0079-6611.

41. Light B, Grenfell TC and Perovich DK 2008. Transmission and absorption of solar radiation by Arctic sea ice during the melt season. Journal of Geophysical Research 113(C3)DOI: http://dx.doi.org/10.1029/2006JC003977 ISSN 0148-0227.

42. Lindsay R and Schweiger A 2015. Arctic sea ice thickness loss determined using subsurface, aircraft, and satellite observations. The Cryosphere 9(1): 269–283, DOI: http://dx.doi.org/10.5194/tc-9-269-2015 ISSN 1994-0424.

43. McDonald S Koulis T Ehn J Campbell K Gosselin M et al. 2015. A functional regression model for predicting optical depth and estimating attenuation coefficients in sea-ice covers near Resolute Passage, Canada. Annals of Glaciology 56(69): 147–154, DOI: http://dx.doi.org/10.3189/2015AoG69A004 ISSN 0260-3055.

44. McLaughlin FA and Carmack EC 2010. Deepening of the nutricline and chlorophyll maximum in the Canada Basin interior, 2003–2009. Geophysical Research Letters 37(24): L24602.DOI: http://dx.doi.org/10.1029/2010GL045459 ISSN 1944-8007.

45. Monod J 1949. The growth of bacterial cultures. Annual Review of Microbiology 3(1): 371–394, DOI: http://dx.doi.org/10.1146/annurev.mi.03.100149.002103

46. Mundy CJ Gosselin M Ehn J Gratton Y Rossnagel A et al. 2009. Contribution of under-ice primary production to an ice-edge upwelling phytoplankton bloom in the Canadian Beaufort Sea. Geophysical Research Letters 36(17): L17601.DOI: http://dx.doi.org/10.1029/2009GL038837 ISSN 1944-8007.

47. Mundy CJ Gosselin M Gratton Y Brown K Galindo V et al. 2014. Role of environmental factors on phytoplankton bloom initiation under landfast sea ice in Resolute Passage, Canada. Marine Ecology Progress Series 497: 39–49, DOI: http://dx.doi.org/10.3354/meps10587 ISSN 0171-8630, 1616-1599.

48. Niemi A, Michel C, Hille K and Poulin M 2011. Protist assemblages in winter sea ice, Protist assemblages in winter sea ice: setting the stage for the spring ice algal bloom. Polar Biology 34(12): 1803–1817, DOI: http://dx.doi.org/10.1007/s00300-011-1059-1 ISSN 0722-4060, 1432-2056.

49. Onodera J, Watanabe E, Harada N and Honda MC 2015. Diatom flux reflects water-mass conditions on the southern Northwind Abyssal Plain, Arctic Ocean. Biogeosciences 12(5): 1373–1385, DOI: http://dx.doi.org/10.5194/bg-12-1373-2015 ISSN 1726-4189.

50. Palmer MA, Saenz BT and Arrigo KR 2014. Impacts of sea ice retreat, thinning, and melt-pond proliferation on the summer phytoplankton bloom in the Chukchi Sea, Arctic Ocean. Deep Sea Research Part II: Topical Studies in Oceanography 105: 85–104, DOI: http://dx.doi.org/10.1016/j.dsr2.2014.03.016 ISSN 0967-0645.

51. Papakyriakou T and Miller L 2011. Springtime CO2 exchange over seasonal sea ice in the Canadian Arctic Archipelago. Annals of Glaciology 52(57): 215–224, DOI: http://dx.doi.org/10.3189/172756411795931534

52. Perovich DK, Grenfell TC, Light B and Hobbs PV 2002. Seasonal evolution of the albedo of multiyear Arctic sea ice. Journal of Geophysical Research: Oceans 107(C10): 8044.DOI: http://dx.doi.org/10.1029/2000JC000438 ISSN 2156-2202.

53. Pogson L, Tremblay B, Lavoie D, Michel C and Vancoppenolle M 2011. Development and validation of a one-dimensional snow-ice algae model against observations in Resolute Passage, Canadian Arctic Archipelago. Journal of Geophysical Research 116(C4)DOI: http://dx.doi.org/10.1029/2010JC006119 ISSN 0148-0227.

54. Rebreanu L, Vanderborght JP and Chou L 2008. The diffusion coefficient of dissolved silica revisited. Marine Chemistry 112(3–4): 230–233, DOI: http://dx.doi.org/10.1016/j.marchem.2008.08.004 ISSN 0304-4203.

55. Redfield AC, Ketchum BH and Richards FA 1963. The influence of organisms on the composition of sea-water In: The composition of sea-water: comparative and descriptive oceanography, The sea, ideas and observations on progress in the study of the seas. New York: Interscience Pub, 2pp. 26–77.

56. Riedel A, Michel C and Gosselin M 2006. Seasonal study of sea-ice exopolymeric substances on the Mackenzie shelf: implications for transport of sea-ice bacteria and algae. Aquatic Microbial Ecology 45: 195–206, DOI: http://dx.doi.org/10.3354/ame045195 ISSN 0948-3055, 1616-1564.

57. Różańska M, Poulin M and Gosselin M 2008. Protist entrapment in newly formed sea ice in the Coastal Arctic Ocean. Journal of Marine Systems 74(3–4): 887–901, DOI: http://dx.doi.org/10.1016/j.jmarsys.2007.11.009 ISSN 0924-7963.

58. Shirasawa K and Ingram RG 1997. Currents and turbulent fluxes under the first-year sea ice in Resolute Passage, Northwest Territories, Canada. Journal of Marine Systems 11(1–2): 21–32, DOI: http://dx.doi.org/10.1016/S0924-7963(96)00024-3 ISSN 0924-7963.

59. Sibert V Zakardjian B Saucier F Gosselin M Starr M et al. 2010. Spatial and temporal variability of ice algal production in a 3D ice-ocean model of the Hudson Bay, Hudson Strait and Foxe Basin system: Modelling ice algal production in the Hudson Bay system. Polar Research 29(3): 353–378, DOI: http://dx.doi.org/10.1111/j.1751-8369.2010.00184.x ISSN 0800-0395.

60. Smith REH 1988. Abundance and production of ice algae in Resolute Passage, Canadian Arctic. Mar Ecol Prog Ser 48: 251–263, DOI: http://dx.doi.org/10.3354/meps048251

61. Smith RE, Harrison WG, Harris LR and Herman AW 1990. Vertical fine structure of particulate matter and nutrients in sea ice of the high Arctic. Canadian Journal of Fisheries and Aquatic Sciences 47(7): 1348–1355, DOI: http://dx.doi.org/10.1139/f90-154

62. Steiner N Deal C Lannuzel D Lavoie D Massonnet F et al. 2016. What sea-ice biogeochemical modellers need from observers. Elementa: Science of the Anthropocene 4: 000084.DOI: http://dx.doi.org/10.12952/journal.elementa.000084 ISSN 2325-1026.

63. Steiner N, Denman K, McFarlane N and Solheim L 2006. Simulating the coupling between atmosphere–ocean processes and the planktonic ecosystem during SERIES. Deep Sea Research Part II: Topical Studies in Oceanography 53(20–22): 2434–2454, DOI: http://dx.doi.org/10.1016/j.dsr2.2006.05.030 ISSN 0967-0645.

64. Steiner NS Sou T Deal C Jackson JM Jin M et al. 2015. The future of the subsurface chlorophyll-a maximum in the Canada Basin – a model intercomparison. Journal of Geophysical Research: Oceans, : n/a–n/a, DOI: http://dx.doi.org/10.1002/2015JC011232 ISSN 2169-9275.

65. Suzuki Y and Takahashi M 1995. Growth responses of several diatom species isolated from various environments to temperature. Journal of Phycology 31(6): 880–888, DOI: http://dx.doi.org/10.1111/j.0022-3646.1995.00880.x ISSN 1529-8817.

66. Taylor PD and Feltham DL 2004. A model of melt pond evolution on sea ice. Journal of Geophysical Research: Oceans 109(C12): 007.DOI: http://dx.doi.org/10.1029/2004JC002361 ISSN 2156-2202.

67. Tedesco L, Vichi M and Thomas DN 2012. Process studies on the ecological coupling between sea ice algae and phytoplankton. Ecological Modelling 226: 120–138, DOI: http://dx.doi.org/10.1016/j.ecolmodel.2011.11.011 ISSN 0304-3800.

68. Tremblay JÉ Simpson K Martin J Miller L Gratton Y et al. 2008. Vertical stability and the annual dynamics of nutrients and chlorophyll fluorescence in the coastal, southeast Beaufort Sea. Journal of Geophysical Research: Oceans 113(C7): S90.DOI: http://dx.doi.org/10.1029/2007JC004547 ISSN 2156-2202.

69. Vancoppenolle M Goosse H de Montety A Fichefet T Tremblay B et al. 2010. Modeling brine and nutrient dynamics in Antarctic sea ice: The case of dissolved silica. Journal of Geophysical Research 115(C2)DOI: http://dx.doi.org/10.1029/2009JC005369 ISSN 0148-0227.

70. Vancoppenolle M Meiners KM Michel C Bopp L Brabant F et al. 2013. Role of sea ice in global biogeochemical cycles, Role of sea ice in global biogeochemical cycles: emerging views and challenges. Quaternary Science Reviews 79: 207–230, DOI: http://dx.doi.org/10.1016/j.quascirev.2013.04.011 ISSN 0277-3791.

71. Vaughan D Comiso J Allison I Carrasco J Kaser G et al. 2013. Observations: Cryosphere In: Stocker, T Qin, D Plattner, GK Tignor, M Allen, S et al. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge, United Kingdom and New York, NY, USA: Cambridge University Press, pp. 317–382. 978-1-107-66182-0

72. Zeebe RE, Eicken H, Robinson DH, Wolf-Gladrow D and Dieckmann GS 1996. Modeling the heating and melting of sea ice through light absorption by microalgae. Journal of Geophysical Research: Oceans 101(C1): 1163–1181, DOI: http://dx.doi.org/10.1029/95JC02687 ISSN 2156-2202.