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.