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

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.

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

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

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

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.

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.

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

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.

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

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

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

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

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

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

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.

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 (w_d2; 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 w_d2 = 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.

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.

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

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 NO₃ 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.

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:

where b_ia, [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). m_lia represents the rate constant for the temperature-dependent linear mortality and m_qia 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, m_lia and m_qia 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).

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.

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 (⁠$pno3ia$ 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 (NH₄ to NO₃) 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.

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.

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:

where h is the snow depth (m) and h_s 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:

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). I₀ 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:

where SWR is the incident shortwave radiation (W m^–2), α is the surface albedo (dimensionless), and i₀ 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:

where α_s, i_0,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 (A_s), bare ice (A_i), and melt ponds (A_m), such that:

The parameterizations for A_s, A_i, and A_m 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:

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

where α_i and α_m are the albedos and i_0,i and i_0,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:

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 $(PAR¯bot)$ is defined as:

where h_i is the sea ice thickness. Finally, the gridbox-mean PAR in the water column under the ice $(PAR¯ui(z))$ 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:

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

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.

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.

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

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 (L_nit, L_sil, L_lig, and L_ice):

where $μiamax$ represents the specific growth rate and b_ia is the temperature sensitivity coefficient for the ice algal growth. Following Lavoie et al. (2005), $μiamax$ and b_ia are set to 0.85 d^–1 and 0.0633°C^–1. The nitrogen (L_nit) and silicon (L_sil) limitation functions are based on the Monod formulation (Monod, 1949):

where [NO₃]_ia, [NH₄]_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. k_n and k_s 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 (L_lig) is formulated following Lavoie et al. (2005):

where α^ia is the photosynthetic efficiency (mg C (mg Chl a)^–1 h^–1 (µmol photons m^–2 s^–1)^–1), $Pmia$ 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 $Pmia$ depending on snow cover. We set the ratio of these two photosynthetic parameters $(αia/Pmia)$ to 0.44, based on the results of our sensitivity studies.

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

where [R_ice] 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 C_m 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:

where m_lia is the rate constant for the temperature-dependent linear mortality and m_qia is the rate constant for the quadratic mortality. We set m_lia and m_qia 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 (R_ice) and the areal fraction of melt ponds (f_pond) are computed by the physical model (see Flato and Brown (1996) and Abraham et al. (2015) for details). r_pond 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, z_ia is the thickness of the ice skeletal layer, which is set to 0.03 m in this study.

The equations for nitrate ([NO₃]_ia) and ammonium ([NH₄]_ia) in the ice skeletal layer (mmol N m^–3) are defined as:

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 (⁠$pno3ia$⁠) based on Denman (2003) is defined as:

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 [NO₃] and [NH₄] 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 m² s^–1 (Rebreanu et al., 2008). The exchange rate at the ice-water interface depends on the molecular sublayer thickness (h_ν) defined as:

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

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

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

where u_o is computed by the physical ocean model, and we assume u_i = 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 (f_lia) 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, f_lia is set to 0.3.

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

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

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:

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:

where $μp1max$ and $μp2max$ 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 (Q_10p) 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):

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

where α^p1 and α^p2 represent the photosynthetic efficiency and $Pmp1$ and $Pmp2$ represent the maximum photosynthetic rate. In this study, the ratios of photosynthetic parameters (i.e., $αp1/Pmp1$ and $αp2/Pmp2$⁠) 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:

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 (Q_10z) is set to 2.14 (Buitenhuis et al., 2006). k_z1 and k_z2 represent the half-saturation constants for grazing, which are set to 1 mmol N m^–3 for both zooplankton groups in this study. F_z1 and F_z2 respectively represent the available food sources for small and large zooplankton, which are defined as:

where f_d1 represents the grazing preference of small zooplankton on small detitus. f_d2 and f_z2 represent the grazing preferences of large zooplankton on small and large detritus, respectively. f_d1, f_d2, and f_z1 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:

where m_lp1 and m_lp2 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. K_z 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, m_qp2. In this study, m_qp2 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 (f_seed) of ice algae flushed into the water column enters the large phytoplankton pool, while the remaining fraction (1 – f_seed) enters the large detritus pool. In this study, f_seed is set to 0.1, assuming that 10% of ice algae flushing results in large phytoplankton seeding. h_oc is the thickness of the uppermost layer of the ocean model (i.e., 1 m in this study), and $δz,z1$ is the Kronecker’s delta which equals 1 in the uppermost layer of the ocean (z₁), while it is set to 0 elsewhere:

The equations for small and large zooplankton are defined as:

where [Z₁] and [Z₂] 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. a₁ and a₂ 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 (M_lz1 and M_lz2) are parameterized in the same way as those of phytoplankton:

where the rate constants for linear mortality of small and large zooplankton (m_lz1 and m_lz2) 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 (m_qz2) 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:

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 (R_d1, R_d2, R_bsi) depend on seawater temperature following the Q10 formulation:

where r_d1, r_d2, and r_bsi 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:

where [NO₃], [NH₄], 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):

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

where the rate constant (r_nit) and the temperature sensitivity coefficient (b_nit) 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.

• Environmental forcing data: http://climate.weather.gc.ca/

• The source code and the configuration that can be used to reproduce the model results of the present study are available via https://doi.org/10.5281/zenodo.825274.

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.

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.

The authors have no competing interests to declare.

• 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