2.1 Ocean and sea-ice physics (OPA–LIM2)

The physical ocean sub-model is the Océan PArallélisé (OPA), which is a free-surface, hydrostatic, primitive-equation model developed for regional and global ocean circulation studies (Madec, 2008). OPA is coupled to the sub-model for sea-ice physics, namely the Louvain-la-Neuve Sea Ice Model (LIM). The present study uses version 2 of LIM (LIM2; Fichefet and Maqueda, 1997; Bouillon et al., 2009), consisting of a three-layer (one for snow and two for ice) dynamic-thermodynamic model.

Figure 1Shortwave radiative transfer through snow and sea ice in the LIM2/LIM3 model. F_sw represents the incoming shortwave radiation, a fraction of which is reflected due to the surface albedo of snow or ice (a). The remaining radiation is either absorbed within the surface scattering layer (SSL) ($(\mathrm{1}-a)(\mathrm{1}-i_{\mathrm{0}})F_{\text{sw}}$; blue arrows) or penetrates below the SSL ((1−a)i₀F_sw; red arrows). (a) When the SSL is 10 cm (i.e., the sum of snow depth and ice thickness is equal to or greater than 10 cm), the radiation penetrating into the snow and/or ice interior attenuates following the Beer–Lambert law and reaches the water column (magenta arrow). (b) When the SSL is less than 10 cm, the penetrating radiation directly reaches the water column. Modified from Fig. 3.4 of Vancoppenolle et al. (2012).

Download

To model ambient light available for ice algae and under-ice phytoplankton properly, we modified the module that computes the shortwave radiative transfer through snow and sea ice as shown schematically in Fig. 1. In this module, the unreflected fraction (1−a) of the incoming shortwave radiation (F_sw) is parameterized as either being absorbed within a thin layer of surface snow and/or ice or penetrating below this layer. This thin layer at the surface is known as the surface scattering layer (SSL; Grenfell and Maykut, 1977), and is defined in the model as the uppermost 10 cm of snow and/or ice column in NEMO v3.4. When the sum of snow depth and ice thickness is less than 10 cm, the SSL equals this total thickness. The penetrating fraction is determined by the coefficient i₀, which is set to zero in the presence of snow in the default configuration of LIM2 following Maykut and Untersteiner (1971). While this assumption of complete blockage of light may be a reasonable approximation for thermodynamic processes of snow and sea ice, it is problematic for modelling sea-ice biogeochemistry. Specifically, the assumption implies that primary producers can not photosynthesize until snow disappears completely, which is inconsistent with the findings of many field observations that measure high algal biomass at the base of snow-covered sea ice (e.g., Leu et al., 2015). Furthermore, i₀ for snow surface has been set to non-zero values in other sea-ice models in the case of thin or melting snow (Flato and Brown, 1996; Abraham et al., 2015). For these reasons, we use a non-zero value of i₀ for snow surface and parameterize the light transmission through the snow column below the specified surface layer following the Beer–Lambert law. The value of i₀ for snow surface was set to 0.15 following the 1-D sea-ice biogeochemical modelling work of Vancoppenolle et al. (2010). The attenuation coefficient of snow was set to 10 m⁻¹, which falls within the observed range for melting and freezing snow (Grenfell and Maykut, 1977). Model sensitivity to i₀ for snow surface is discussed in Sect. 4.2.

2.2 Ocean biogeochemistry (CanOE)

The sub-model for ocean biogeochemistry adopted in the present study is the Canadian Ocean Ecosystem Model (CanOE), developed by the ocean modelling group at the Canadian Centre for Climate Modelling and Analysis (Christian et al., 2019; see Appendix A3 of Hayashida, 2018a). This model has been developed for the latest version of the Canadian Earth System Model (CanESM5), which will be used in the next phase of the Coupled Model Intercomparison Project (CMIP6). CanOE simulates the lower trophic levels of marine ecosystems (nutrients, phytoplankton, zooplankton, and detritus) and biogeochemical cycling of key elements (carbon, oxygen, nitrogen, and iron). This model is built around the basic code structure of the Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) version 2, the default sub-model for ocean biogeochemistry of NEMO (Aumont et al., 2015). CanOE is an updated version of the Canadian Model of Ocean Carbon (CMOC) used in an earlier version of CanESM (CanESM2; Arora et al., 2011) designed to limit the complexity as much as possible, while addressing insufficiencies in CMOC (e.g., single nutrient–phytoplankton–zooplankton–detritus models are either tuned towards low-nutrient or high-nutrient conditions). CanOE is somewhat more efficient than PISCES as a result of having fewer state variables (19 vs. 24) and fewer computationally expensive parameterizations (Christian et al., 2019; Appendix A3 of Hayashida, 2018a).

In the present study, we made two modifications to CanOE. The first modification is the addition of an ocean sulfur-cycle model and the second modification is the parameterization of the photosynthetically active radiation (PAR).

2.2.1 Addition of an ocean sulfur cycle

Figure 2 shows a schematic of CanOE including the ocean sulfur cycle and the sea-ice biogeochemistry. Here, the sulfur cycle is restricted to sources and sinks of two state variables: dissolved dimethylsulfoniopropionate (DMSPd) and dimethylsulfide (DMS). The ocean sulfur cycle is one-way coupled to CanOE as the sources and sinks of DMSPd and DMS depend on the conditions of primary and secondary producers, but not vice versa. The ocean sulfur-cycle model is based on Hayashida et al. (2017) with the following two modifications. First, the cellular DMSP content of modelled phytoplankton is derived from their carbon content as opposed to the chlorophyll content as in Hayashida et al. (2017). This change was made because there are more observation-based estimates of the intracellular DMSP-to-carbon (DMSP to C) ratio than the DMSP-to-chlorophyll a ratio (e.g., Stefels et al., 2007). The DMSP:C ratios for small and large phytoplankton (respectively high and low DMSP producers) are set to 12 and 4 mmol mol⁻¹.

Figure 2Schematic of the CanOE pelagic ecosystem model and associated sea-ice biogeochemistry and pelagic sulfur-cycle models. Black arrows indicate fluxes of carbon (C), nitrogen (N), and iron (Fe) between compartments; blue arrows indicate sources of dissolved dimethylsulfoniopropionate (DMSPd); gray arrows indicate ice–ocean fluxes of nitrate (NO₃), ammonium (NH₄), ice algae (IA) and large phytoplankton (P_L), DMSPd, and dimethylsulfide (DMS). Flows of dissolved oxygen (O₂) are opposite to those of dissolved inorganic carbon (DIC) and are not explicitly illustrated. Detritus (D_S and D_L) and zooplankton (Z_S and Z_L) are denominated in C units but have implicit N and Fe pools according to fixed elemental ratios; phytoplankton (P_S and P_L) have separate state variables for each currency. O₂ and total alkalinity (TA) have their own currencies, but are shown as white here for simplicity; their sources and sinks follow well established stoichiometries relative to those of DIC. Sources and sinks of TA associated with the nitrogen cycle (Wolf-Gladrow et al., 2007) are included but not shown in the figure. The state variables dFe and CaCO₃ represent dissolved iron and calcium carbonate, respectively. The currencies Chl and S represent chlorophyll a and sulfur, respectively.

Download

Also, the parameterization of sea-to-air flux of DMS was modified to account for the non-linear dependence of the flux on the open-water fraction (Loose et al., 2009):

$$\begin{array}{}\text{(1)}& {\displaystyle }F=f_{\mathrm{ow}}^{\mathrm{0.4}}{k}_{\mathrm{dms}}\mathrm{DMS},\end{array}$$

where F is the DMS flux (µmol m⁻² s⁻¹), f_ow is the open-water fraction (–), k_dms is the gas transfer velocity (m s⁻¹), and DMS (nM) is the DMS concentration in the uppermost layer of the water column.

2.3 Sea-ice biogeochemistry

The sub-model for sea-ice biogeochemistry is a modified version of a three-compartment (ice algae, nitrate, and ammonium) ecosystem based on Mortenson et al. (2017) and a two-compartment (DMS and DMSPd) sulfur cycle based on Hayashida et al. (2017).

Sea-ice biogeochemical processes are assumed to take place in a layer of fixed thickness at the ice base. Hence, this bottom-ice biogeochemical layer is not explicitly modelled and does not correspond to one of the two ice layers in LIM. Although algal biomass in ice core samples above this layer can be substantial (e.g., Melnikov et al., 2002; Olsen et al., 2017), resolving vertical distributions of sea-ice biogeochemistry in 3-D models is computationally impractical at present. The governing equation for any sea-ice biogeochemical state variable is

where X denotes the concentration of the state variable, $\overrightarrow{U}$ denotes the horizontal velocity field of sea ice, and D denotes the horizontal eddy diffusion coefficient. The first two terms on the right-hand side of Eq. (2) represent tendencies associated with horizontal motion of sea ice (Sect. 2.3.1). The third term represents biological and chemical source minus sink (SMS). Note that while LIM2 computes the impact of mechanical redistribution (i.e., deformation due to ridging/rafting) on sea-ice physical properties, these processes are neglected in computations of the sea-ice biogeochemical state variables in the present study as the model uses a simple representation of the sea-ice biogeochemical layer as a layer of fixed thickness (3 cm) at the ice base.

2.4 Spatial domain

The model domain is based on the North Atlantic and Arctic (NAA) configuration developed by the ocean modelling group at the University of Alberta (http://knossos.eas.ualberta.ca/anha/model_configuration.php#naa, last access: 13 May 2019). This configuration was built on the curvilinear orthogonal coordinate system of NEMO that has been successfully applied to study the freshwater budget of the Arctic Ocean in present (Hu and Myers, 2013) and future climates (Hu and Myers, 2014), as well as to investigate pelagic ecosystem processes in the Canada Basin (Steiner et al., 2015). The NAA domain includes the Arctic Ocean, the Canadian Arctic Archipelago, the northern Bering Sea, the North Atlantic Ocean, and the Nordic Seas (Fig. 3). The horizontal resolution of the 568×400 grid varies from 10 km along the North American boundary to 14.5 km along the Eurasian boundary. Vertically, the ocean is divided into 46 layers with variable resolution, from approximately 1 m in the uppermost layer to 255 m in the bottommost layer. This vertical resolution is finer than that of the original NAA configuration in the upper layers, while it is slightly coarser in the deeper layers (approximately below the 20th level; Fig. 4). The bathymetry is based on the 1 arcmin global relief data (ETOPO1; Amante and Eakins, 2009) as described by Hu and Myers (2013). For numerical stability, each ocean grid cell is set to have at least seven vertical levels, corresponding to a depth of approximately 20 m.

Figure 3The domain of the North Atlantic and Arctic (NAA) configuration. The colour map represents the horizontal resolution and the contour lines denote the isobaths at 100 m (red), 1000 m (white), and 2000 m (cyan). The thick (thin) solid black lines indicate the locations of Atlantic and Pacific open (North American and Eurasian closed) boundaries.

Download

Figure 4Comparison of the vertical resolution of the ocean model between the original NAA configuration (NAA6, i.e., approximately 6 m in the uppermost layer) and the configuration adopted in the present study (NAA1, i.e., approximately 1 m in the uppermost layer). Note the log scale on the x axis.

Download

2.5 Experiments

We consider six model experiments (Table 2). The first experiment is a reference simulation (EXP0), which is intended as the most realistic among all simulations considered in this study. The 11-year duration of EXP0 is considered sufficient for the spin-up of sea-ice and near-surface pelagic variables based on previous Arctic biogeochemical model studies (e.g., Dupont, 2012; Jin et al., 2012). The results during this spin-up period are presented in Appendix B. The setup of EXP0 is described below. The remaining experiments (EXP1–5) are designed to assess the sensitivity of the model simulations to changes in uncertain forcing data and parameter values.

Table 2List of model experiments.

Download Print Version | Download XLSX

2.5.1 Initial and lateral boundary conditions, runoff, and atmospheric forcing

The ocean was initialized from rest with temperature and salinity fields for January 1969 derived from the Ocean Reanalysis System 4 (ORAS4; Balmaseda et al., 2013). The initial snow depth, ice thickness, and ice concentration were respectively set to 0.1 m, 2.5 m, and 0.95 for grid cells with temperatures within 2 ^∘C of the seawater freezing point. Elsewhere, these values were set to zero. The initial concentrations of nitrate, dissolved inorganic carbon, and total alkalinity were taken from the annual-mean fields of the GLobal Ocean Data Analysis Project version 2 (GLODAP2; Lauvset et al., 2016). The initial concentrations of dissolved oxygen were set to the annual-mean fields from the World Ocean Atlas 2013 Version 2 (WOA13; Garcia et al., 2014). The initial concentration of dissolved iron was set to 0.6 nM in the entire domain (Aumont et al., 2015). Because the model simulation starts at a time of low biological production (i.e., 1 January), the remaining biogeochemical state variables in the ocean were initialized uniformly in space to very low values (e.g., 0.01 mmol C m⁻³ for the carbon contents of phytoplankton, zooplankton, and detritus). The initial concentrations of sea-ice biogeochemical state variables were set to the same values as their respective variables in the uppermost layer of the ocean.

Open boundary conditions were applied by a radiation-relaxation algorithm (Madec, 2008) along the Atlantic and Pacific boundaries of the model domain, while the other two boundaries (along North America and Eurasia) were assumed to be closed (Fig. 3). The boundary temperature, salinity, and zonal and meridional current fields were interpolated from the interannual monthly-mean fields of ORAS4. The open boundary conditions for ocean biogeochemical state variables were the same as their initial conditions. The relaxation timescales were set to 1 d for inflow and 15 d for outflow. These values are identical to those used in Dupont et al. (2015), but differ from the original NAA configuration (Hu and Myers, 2013). Our preliminary experiments suggested that these changes were needed to prevent salinity drift. Because the feature to prescribe the open boundary conditions for the sea-ice prognostic variables was not available in NEMO version 3.4, these were set to zero for the sea-ice prognostic variables of LIM2 as well as the sea-ice biogeochemical state variables; this feature is available in the subsequent version of NEMO (version 3.6).

Figure 5River runoff of freshwater prescribed in the model. (a) Annual cycle of daily discharge and (b) interannual variability in annual discharge integrated over the region north of 60^∘ N. (c) Spatial distribution of annual discharge rate averaged over the period 1969–1979. In panel (a), the error bars indicate the standard deviations over the period 1969–1979. In panel (c), note the log scale on the colour bar, and names of major rivers and countries or regions are labelled.

Download

River discharge of freshwater was derived from the interannual monthly-mean product of Dai and Trenberth (2002). Figure 5 shows the seasonal and interannual variability (panels a and b) and spatial distribution (c) of the total discharge over the pan-Arctic. The river discharge of biogeochemical state variables was neglected due to the lack of adequate data. Additional external supplies of nutrients (dust deposition and sediment mobilization) were neglected due to the lack of reliable data. Partial pressure of carbon dioxide in the atmosphere was derived from the monthly-mean Mauna Loa CO₂ data (https://www.esrl.noaa.gov/gmd/ccgg/trends/data.html, last access: 19 April 2017).

The surface atmospheric conditions used to drive the sea-ice and ocean model simulations were derived from the Drakkar forcing set 5.2 (DFS; Dussin et al., 2016). The DFS dataset has high resolutions in space (0.7^∘) and time (3-hourly for zonal and meridional wind speed at 10 m height as well as air temperature and specific humidity at 2 m height, and daily for incoming shortwave and longwave radiation, total precipitation, and snowfall). It is based on a combination of ERA-40 and ERA-Interim reanalysis products (Uppala et al., 2005; Dee et al., 2011). The original DFS dataset (https://www.drakkar-ocean.eu/forcing-the-ocean, last access: 13 May 2019) has missing data flags which cause a simulation crash in some years. As a substitute, we used a version provided by Clark Pennelly at the University of Alberta (personal communication) which addressed the missing data flag errors without any modifications to the atmospheric data (the only changes were indexing and ordering of latitudinal coordinates and remaining variables). In the original DFS dataset, the daily total precipitation (rain and snow) and snowfall fields prior to 1979 were set to their respective 1979–2012 climatological daily averages due to the lack of adequate observations to construct the dataset for those years individually (Dussin et al., 2016). However, in EXP0, we prescribed the total precipitation and snowfall for 1979 repeatedly for the simulation over the period 1969–1978, while keeping the remaining atmospheric variables the same as the original DFS dataset. This modification was necessary to simulate adequate snow depths (discussed further in Sect. 4.1).

2.5.2 Additional settings

The time step of the model integration was 20 min. Unlike Hu and Myers (2013), no additional treatments for modelled temperature, salinity, or wind-stress fields near the open boundaries were necessary since no obvious drift was apparent in the simulated fields. Table 3 displays some of the model parameters that were modified from their default values in NEMO v3.4. For a complete list of the parameters, readers are referred to the source-code repository (https://doi.org/10.5281/zenodo.1435254). The coefficients for horizontal eddy diffusion for oceanic and sea-ice tracers (rn_aht_0, ahi0, and rn_ahtrc_0) were reduced to keep diffusion relatively small compared to resolved dynamical processes, as recommended by Vancoppenolle et al. (2012). The other two parameters (hiccrit and pstar) were adjusted to improve the fit with the PIOMAS data product (Sect. 2.6) in terms of sea-ice volume and extent for 1979 (Sect. 3.1.1). While parameter hiccrit is very large considering that it is the initial thickness of newly formed ice, it is also known as a tuning parameter for sea-ice extent in LIM (Wang et al., 2010; Roy et al., 2015). Hence, we chose this value while trying to keep it low as much as possible. Lastly, two parameters of CanOE (Tref and chldeg) were adjusted to simulate reasonable annual primary production in the Arctic Ocean (Sect. 3.2). Note that chldeg is a parameter for photo-oxidation designed specifically for the global configuration of CanOE to tune the global chlorophyll distribution, and therefore we set this to zero in our Arctic configuration.

Table 3List of selected model parameters in the NEMO namelists.

Download Print Version | Download XLSX

2.6 PIOMAS data product

The modelled sea-ice properties were evaluated against the output of the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS), which is a regional coupled sea-ice–ocean circulation model that assimilates some observational data (Zhang and Rothrock, 2003; Schweiger et al., 2011). The monthly-mean ice thickness and ice concentration gridded data products of PIOMAS were interpolated onto the NAA grid in order to perform a grid-to-grid comparison across the same domain. Although the PIOMAS data product is considered here as the best presently available, note that it has its own biases that could result in mismatches with our model results.

3.1 Comparison of sea-ice physical properties with PIOMAS in 1979

3.1.1 Seasonal variability

To assess the model performance in simulating sea ice, the seasonal variability of modelled ice volume and extent in EXP0 are compared to PIOMAS (Fig. 6a and b) for 1979 (after a decadal spin-up). In both EXP0 and PIOMAS, the annual maximum in ice volume (extent) takes place in April (March), while both the ice volume and extent are at their annual minima in September. The ice volume is consistently higher in EXP0 than PIOMAS. The difference in the annual-mean ice volume over the NAA domain is 3.9 km³ (17 %). In contrast, the ice extent comparison is much better with the difference of 0.1×10⁶ km² (1 %) in the annual means.

Figure 6Time series of 5 d mean modelled (a) snow and ice volumes, (b) ice extent, and pan-Arctic-mean surface seawater nitrate concentration, and (c) pan-Arctic ice algal and phytoplankton NPP during 1979 in EXP0. The dashed lines in panels (a) and (b) represent the monthly-mean PIOMAS ice volume and extent, respectively.

Download

Figure 7Spatial distributions of monthly-mean ice thickness in EXP0 (a, d) and the PIOMAS product interpolated onto the NAA grid (b, e) and their difference (c, f) for March and September in 1979. The red lines represent the ice edge, defined here as the 0.15 contour of ice concentration.

Download

3.2 Primary productivity of ice algae and phytoplankton

3.2.1 Seasonal variability

Figure 6 shows the seasonal variability in modelled pan-Arctic-mean ice algal NPP and phytoplankton NPP during 1979 (panel c) along with relevant environmental factors (panels a and b). Ice algal NPP starts increasing in early February, peaks in mid-May, sharply declines in late May to early June, and is near zero by late June. The start of the decline of the ice algal NPP coincides with the decline of the ice volume (Fig. 6a) demonstrating that the decline is driven by the release of ice algae as a result of ice melt. The seasonal progression of the ice algal production is similar to Jin et al. (2012). The phytoplankton NPP starts increasing in early March, peaks in early July, and decreases to near zero by the end of October (Fig. 6d). At the peak in phytoplankton NPP, the pan-Arctic-mean surface seawater nitrate concentration is below 1 mmol N m⁻³ and remains so until the end of August (Fig. 6b).

Figure 8Spatial distribution of annual-mean (a) snow depth and (b) surface seawater nitrate concentration, and (c) depth-integrated (bottom 3 cm) ice algal annual NPP and (d) depth-integrated (upper 90 m) phytoplankton annual NPP in 1979 in EXP0. The red and black dashed lines represent the 0.15 contour of monthly-mean ice concentration in March and September, respectively.

Download

3.3 Vertical distribution of salinity, nitrate, chlorophyll a, and DMS in the upper water column

The seasonal variability in pan-Arctic-mean seawater salinity, nitrate, chlorophyll a, and DMS in the upper 15 m of the water column is shown in Fig. 9. During the summer, a prominent freshening of the uppermost layer occurs as a result of ice melt (Fig. 9a). This freshening results in the formation of a thin layer of low-salinity water known as a meltwater lens, which strengthens stratification and reduces mixing with the underlying water column. The formation of the lens coincides with the bloom of modelled phytoplankton, resulting in the depletion of nitrate first in the uppermost model layer and then in the underlying layers (Fig. 9b). Nutrient depletion in the near-surface waters then results in the formation of subsurface chlorophyll a and DMS maxima during the latter half of July (Fig. 9c and d). Note that the meltwater lens and the subsurface maxima are respectively thicker and shallower than those observed by field measurements (e.g., Brown et al., 2015) because of averaging over the pan-Arctic domain. The purpose of this spatial averaging is to quantify the impacts at a larger scale rather than assessing localized effects. Note that the model does simulate reasonable subsurface chlorophyll a maxima compared to observations (Appendix C). On the other hand, the modelled salinity distribution shows a fresh bias compared to the Polar Hydrographic Climatology version 3.0 (PHC3.0; Steele et al., 2001) throughout the year in the upper water column and most substantially (>2 psu) in the surface layer (Fig. S1 in the Supplement).

Figure 9Time series of 5 d mean and pan-Arctic-mean seawater (a) salinity, and concentrations of (b) nitrate, (c) chlorophyll a, and (d) DMS in the upper 15 m of the water column during April–September in 1979 of EXP0.

Download

These ice-associated physical and biogeochemical processes take place within a relatively shallow upper-water column (∼10 m), and would have been impossible to simulate with a model of coarse vertical resolution. It is for this reason that the near-surface vertical resolution of the NAA configuration considered in the present study is finer than that of the original configuration (6 m in the uppermost layer; Hu and Myers, 2013). Although modelling these small-scale processes probably has a negligible effect on bulk quantities such as depth-integrated NPP, it can have an impact on processes at the air–sea or ice–sea interface (e.g., gas fluxes). To illustrate this point, the time series of modelled pan-Arctic-mean seawater DMS concentration in the uppermost layer of the water column (about 1 m) is compared with the concentration averaged over the top four layers (about 12 m) as a proxy for values simulated by a coarse-vertical-resolution model (Fig. 10).

Figure 10(a) Time series of 5 d mean and pan-Arctic-mean seawater DMS concentration (a) in the uppermost layer (∼1 m; black) and averaged over the upper four layers (∼12 m; gray) during April–September in 1979 of EXP0. (b) The percentage difference between the two time series (the 12 m average minus the 1 m average, divided by the 1 m average).

Download

Modelled DMS concentration is higher in the uppermost layer than the 12 m average throughout most of April–September, while it is slightly lower in August (Fig. 10a). The concentration difference is highest (up to about 20 %) in June–July (Fig. 10b). Overall, the annual-mean DMS concentration averaged over the upper 12 m of the water column is 9 % lower than in the upper 1 m. Here, the averaging over a thicker layer results in dilution of the DMS concentration in the uppermost layer represented in the model. Considering that this difference is present primarily during the ice melt period, and therefore that the sea-surface DMS is released into the atmosphere, the modelled sea-to-air DMS flux would be underestimated by a similar amount in the absence of fine vertical resolution in the upper water column.

4.1 Snowfall forcing frequency (EXP1 and 2)

Two sensitivity experiments (EXP1 and EXP2) are performed with the identical setup as EXP0 except for a change to the atmospheric forcing. In EXP1, all the forcing fields are replaced by the CORE-II dataset as in the original NAA configuration (Hu and Myers, 2013). Note that the temporal resolution of the snowfall and total precipitation fields in the CORE-II dataset is monthly. In EXP2, the snowfall and total precipitation fields over the period 1969–1978 are replaced by their respective 1979–2012 daily climatological values as in the original DFS dataset (Dussin et al., 2016). Comparing between EXP0 and EXP1 allows us to assess the impacts of atmospheric forcing (DFS vs. CORE-II), while comparing between EXP1 and EXP2 allows us to assess the impacts of snowfall dataset (daily vs. daily climatology) on modelled snow depth.

Figure 11Model sensitivity to snowfall forcing frequency. Time series of pan-Arctic-mean (a) prescribed snowfall rate of the CORE-II (blue) and DFS (black/red) datasets and (b) modelled annual-mean snow depth in EXP0 (black), EXP1 (blue), and EXP2 (red). Spatial maps of modelled annual-mean snow depth for the period 1970–1978 in (c) EXP0, (d) EXP1, and (e) EXP2. The units for the snowfall rate are converted from kilogram per square metre per second (kg m⁻² s⁻¹) to millimetre per day (mm d⁻¹) using a constant snow density of 330 kg m⁻³, which is the value assumed in LIM2. In EXP0, the snowfall rate for 1979 is repeated annually throughout 1969–1979.

Download

A comparison of the pan-Arctic-mean snowfall rates between the CORE-II and DFS datasets illustrates the differences between them (Fig. 11a). The monthly CORE-II dataset varies from approximately 1 to 2.4 mm d⁻¹ (meltwater equivalent), while the range of the DFS dataset is 3 times as large (from near 0 to about 4.4 mm d⁻¹ for the year 1979) most likely due to the difference in the temporal resolution of the datasets. The lack of high-frequency variability in the DFS daily climatology is evident from the comparison of the DFS dataset between 1969–1978 and 1979. The daily climatology ranges approximately from 0.2 to 2.2 mm d⁻¹, less than half of the range for the individual daily averages for 1979. The annual-mean CORE-II snowfall rate is higher than that of the DFS dataset in all of these years. The annual mean of the DFS daily climatology is slightly greater than that of the individual daily averages for 1979.

Figure 11b shows a comparison of the modelled pan-Arctic annual-mean snow depth time series among EXP0, EXP1, and EXP2. The snow depth is substantially lower in EXP1 and EXP2 than in EXP0 throughout the period 1969–1979, except for 1969 (in which year the snow depth is affected by its initial value). In EXP2, the extremely low snow depth somewhat recovers in 1979. Figure 11c–e shows a spatial comparison of the modelled annual-mean snow depth over the period 1970–1978 (excluding the first and last year of simulation). There is a clear difference in the distribution between EXP0 and the other two experiments; the ice pack is generally covered by a moderate amount of snow (∼0.1 m) in EXP0, while in EXP1 and EXP2, most regions are nearly snow-free. These results of the latter two experiments are inconsistent with the available snow depth climatology indicating the presence of considerably thicker (>0.2 m) annual-mean snow cover over the Arctic Basin (Warren et al., 1999). As a result of these biases, the modified DFS dataset is used as the reference simulation, rather than the CORE-II dataset or the original DFS dataset.

The results of these sensitivity experiments demonstrate the need for high temporal resolution forcing for snowfall datasets in order to realistically simulate snow depth. The lack of snow accumulation in EXP1 and EXP2 is likely due to the mismatch between the timing of snowfall and the state of ice surface (freezing or melting). The latter is also determined by the forcing dataset (air temperature), which has relatively high temporal resolution (3-hourly for DFS and 6-hourly for CORE-II). Ideally, the temporal resolutions of snowfall and air temperature should be identical for consistency, such that snowfall occurs when air temperature is at or below freezing. Such an effort has been made recently (Tsujino et al., 2018), which should resolve the issue illustrated in our sensitivity experiments. Although the usage of a high-frequency atmospheric forcing dataset is desirable, our further sensitivity experiments indicate that the issue can be resolved by tuning another model parameter, which is presented for readers' interest in Appendix D.

4.2 Light penetration through snow column (EXP3)

Figure 12 compares the modelled sea-ice physical and biogeochemical properties in 1979 in EXP0 with those of EXP3, in which i₀ for snow surface is set to the default LIM2 value of zero. The results for modelled snow and ice volume are almost identical between the two experiments (Fig. 12a), indicating a low sensitivity of these physical quantities to the change in i₀ for the snow surface. On the other hand, an appreciable difference in the modelled bottom-ice PAR prior to the melt season in June results in a large difference in the modelled ice algal NPP (Fig. 12b). Ice algal NPP in EXP3 is restricted to snow-free regions, so an increase in the ice algal NPP due to the change in i₀ for snow surface reflects production in snow-covered regions (Fig. 12c). The pan-Arctic ice algal annual NPP of EXP3 is 3.5 Tg C yr⁻¹, only about a quarter of the value obtained in EXP0. This value is much lower than those obtained in previous studies (see Sect. 3.2.3). This result emphasizes the importance of correct representation of the light penetration through snow, and shows that the original LIM2 provides inadequate light for ice algal growth, resulting in insufficient ice algal NPP. Note that the default value of i₀ for the snow surface is also set to zero in LIM3 (Vancoppenolle et al., 2012).

Figure 12Model sensitivity to light penetration through snow. Time series comparison of modelled 5 d mean (a) snow volume (blue) and ice volume (red) and (b) bottom-ice PAR (blue) and ice algal NPP (red) in 1979 between EXP0 (solid) and EXP3 (dashed). (c) Spatial distribution of the difference in the ice algal annual NPP between EXP0 and EXP3.

Download

Previous 3-D sea-ice biogeochemical models differ in their choices of values for i₀ for the snow surface (Table 1). The studies of Dupont (2012) and Jin et al. (2012) set this value to zero, yet their values for simulated ice algal productivity are relatively high. However, these models use special parameterizations for irradiance and light limitation, respectively, which likely result in realistic ice algal primary production values despite the lack of light penetration through snow. Dupont (2012) imposes a minimum lead fraction of 0.01 in any grid cell, supplying enough ambient light for ice algal growth. In Jin et al. (2012), the light limitation parameter (the ratio of light-limited slope and maximal photosynthetic rate; see Table 2 of Jin et al., 2006) is set to a very high value, nearly double the upper limit of the observed range reported in Table 2 of Lavoie et al. (2005). This reduction in light limitation allows the modelled ice algae to grow faster even under low-light conditions upon snow disappearance.

Figure 13Model sensitivity to the horizontal transport of sea-ice biogeochemical state variables. (a) Time series comparison of 5 d mean and pan-Arctic-mean modelled bottom-ice nitrate (blue) and ice algal daily NPP (red) during January–June of 1979 between EXP0 (solid) and EXP4 (dashed). Spatial maps of the annual-mean bottom-ice nitrate in (b) EXP0 and (c) its difference between EXP0 and EXP4, (d) the difference in the ice algal annual NPP between EXP0 and EXP4, and (e) the magnitude of the ice velocity during May.

Download

Figure 14Model sensitivity to shading by ice algae. (a) Time series comparison of modelled pan-Arctic-mean and 5 d mean under-ice PAR (blue) and NPP (red) between EXP0 (solid) and EXP5 (dashed) during 1979. Spatial maps of (b) monthly-mean under-ice PAR in May in EXP0 and (c) its difference from EXP5, (d) the under-ice annual NPP in EXP0, and (e) its difference from EXP5.

Download

Two other regional modelling studies prescribe non-zero values of i₀ for snow surface. Castellani et al. (2017) set i₀ for snow surface to 0.3 based on the measurements over snow-free ice surface (Grenfell and Maykut, 1977). As such, this value (0.3) should be viewed as an overestimate. Similarly, the light penetration through snow is also overestimated in Watanabe et al. (2015), as i₀ for snow surface is effectively unity in their study. Using these higher i₀ for snow surface reduces light limitation for ice algal growth as long as the formulation for light limitation is inversely proportional to irradiance. However, higher i₀ can result in increasing light limitation, e.g., if the formulation assumes some optimal irradiance level with photoinhibition at higher levels, such as done in Watanabe et al. (2015). As such, the overall impact of i₀ on ice algal production depends on the choice of formulation and parameter(s) for the light limitation function.

To the best of the authors' knowledge, no studies have ever reported an observed value for i₀ for snow surface. For a snow-free ice surface, Grenfell and Maykut (1977) report values ranging between 0.18 and 0.63 depending on both the ice type and whether the incoming shortwave radiation is direct or diffuse. Observation-based estimates of i₀ for snow surface would be useful in order to reduce the uncertainty in ice algal and under-ice phytoplankton growth in models.

4.3 Horizontal transport associated with moving sea ice (EXP4)

EXP4 is designed to quantify the contribution of sea-ice advection and diffusion to the overall budget of sea-ice biogeochemical state variables. This analysis can be useful in deciding a location for conducting 1-D simulations in which horizontal processes of sea-ice biogeochemistry are neglected. Specifically, EXP4 is conducted with the identical model formulation as EXP0 except that the advection and diffusion of sea-ice biogeochemical state variables (the effect of horizontal transport associated with moving sea ice) are artificially suppressed (i.e., the first two terms on the right-hand side of Eq. (2) are removed). Note that the advection and diffusion of sea-ice physical state variables are retained in EXP4, so there is no difference in these variables between EXP0 and EXP4.

A time series comparison of the modelled pan-Arctic-mean bottom-ice nitrate and ice algal NPP for 1979 shows that these quantities are always higher in EXP0 than EXP4 (Fig. 13a). The pan-Arctic annual-mean bottom-ice nitrate and the ice algal annual NPP are higher in EXP0 than EXP4 by 2 % and 16 %, respectively. These results indicate that the overall effect of horizontal transport associated with moving sea ice over the pan-Arctic is an increase in these quantities. However, we note that these values could be quite different in other years given the large interannual variability in wind stress fields driving sea-ice drift patterns.

Although the overall effect is an enhancement, the spatial distribution shows regions of local increase and decrease (Fig. 13c–d). The difference in nitrate concentration between the two experiments is relatively high off the west coast of Baffin Island, where the bottom-ice nitrate concentration is relatively high in EXP0 (Fig. 13b), whereas the difference is relatively small on the Canadian Polar Shelf (Fig. 13c). The difference in ice algal NPP is relatively high in regions of high ice algal NPP except for the Canadian Polar Shelf (Fig. 13d), which is a region of relatively slow ice motion (Fig. 13e). One possible explanation for these spatial differences is the horizontal transport of nutrients into regions of high productivity, which results in further increase in ice algal production.

4.4 Shading of ice algae (EXP5)

In EXP5, the shading effect of ice algae on light transfer through the ice is artificially suppressed in order to assess its impact on under-ice NPP. Effectively, this is done by setting the light extinction coefficient for ice algae to zero (Eq. 15 of Mortenson et al., 2017). Note that this modification has no impact on simulated physics. On the pan-Arctic scale, there is almost no effect, as shown in Fig. 14a. The differences in the pan-Arctic- and annual-mean under-ice PAR and the pan-Arctic under-ice annual NPP between EXP0 and EXP5 are 2 % and 1 %, respectively.

Consistent with the patchiness of the ice algal distribution (Fig. 8c), the shading effect is rather localized as shown in Fig. 14b–e. The influence on under-ice PAR is assessed for the month of the ice algal bloom peak (May; Fig. 6c). The spatial distribution of the difference in the under-ice PAR between EXP0 and EXP5 is simply a reflection of ice algal abundance (Fig. 14c). Similarly, a general decrease in the under-ice NPP is found due to shading in the regions of high ice algal production. However, in some regions, shading results in a slight increase in under-ice NPP which is dominated by small phytoplankton, but a slight increase is also seen for large phytoplankton (Fig. 15c). The underlying mechanism for this response of the modelled ecosystem to a perturbation to light is unclear.

Figure 15Spatial maps of under-ice annual NPP by (a) small and (b) large phytoplankton during 1979 in EXP0, and (c, d) their respective differences from EXP5.

Download

The shading effect of ice algae was recently examined in the model study of Castellani et al. (2017). Their results showed that the effect has greater influence at higher latitudes due to low ambient light. Furthermore, they hypothesized that the onset of the under-ice phytoplankton bloom north of 80^∘ N can be delayed by up to 40 d depending on how their modelled under-ice PAR is affected by shading.

It is difficult to directly compare the results of the present study with those of Castellani et al. (2017), primarily due to the difference in the definition of the term under-ice. As described in Sect. 2.5.3, in the present study, a grid cell is considered “under-ice” as long as the ice concentration is 0.15 or above. Because of the high surface albedo and strong light attenuation by snow and ice, the under-ice PAR defined in the present study is therefore dominated by the light through the open-water fraction (see also Dupont, 2012). Consequently, the under-ice NPP is controlled by the light through the open-water fraction and does not show a strong influence by the shading of ice algae. Furthermore, direct comparison is difficult due to the difference in the target year of simulation; Castellani et al. (2017) simulated 2012, while we consider 1979.

Figure 16Effects of ice algal shading on the onset of under-ice phytoplankton bloom. Spatial maps showing the bloom onset (as the day from 1 January) when the ice algal shading is (a) considered and (b) neglected and (c) the difference between the two cases representing the delay due to the shading in 1979 in EXP0. In panel (c), “No bloom” refers to regions in which the bloom was present in panel (b) but not in panel (a). See Sect. 4.4 for the definition of bloom onset.

Download

Nevertheless, to carry out an analysis comparable to that of Castellani et al. (2017), we calculate the onset of under-ice phytoplankton bloom as follows. A bloom onset is defined as the day when bottom-ice PAR exceeds 0.4 W m⁻² and remains above this value at least for 30 d. This threshold for bottom-ice PAR corresponds to the limit for under-ice algal growth considered in Castellani et al. (2017), assuming a unit conversion (from µmol photon m⁻² s⁻¹ to W m⁻²) factor of 1∕4.56 (Lavoie et al., 2005). Figure 16 shows the spatial variability in the under-ice bloom onset based on the definition above. The bloom takes place mostly in seasonally ice-covered regions, while it is absent in most of the pack ice (as indicated by white regions). Unlike Castellani et al. (2017), the under-ice bloom north of 80^∘ N is absent even without the shading effect. The absence of the bloom in our simulation is due to the presence of snow in this region; despite the extremely low quantity (<0.01 m; data not shown), it keeps the light level below the threshold for the bloom to occur. The median value of the onset is on the 155th day (June 6) when the shading is accounted for (Fig. 16a), while it is 10 d earlier without the shading effect (Fig. 16b). Note that these results are based on the analysis of modelled light transmission through snow and ice columns only, and therefore are not comparable to field observations of under-ice blooms triggered by light transmission through open leads (e.g., Assmy et al., 2017).

Figure 16c shows the spatial variability in the delay in the under-ice bloom onset caused by the ice algal shading. The values range from 5 to 275 d; in some places, the bloom is prevented completely. The present study does confirm the finding of Castellani et al. (2017) that the shading effect is spatially variable and can have a strong impact on the phytoplankton bloom timing under the ice with high ice algal biomass. However, given the patchiness of ice algal distribution (mostly confined to shelf regions) and the control of the light through the open-water fraction, the impact of the shading on the pan-Arctic under-ice annual NPP is negligible.

Besides the shading effect, ice algae can contribute to ice melting through light absorption and conversion to heat (Zeebe et al., 1996). However, our model does not incorporate this bio-physical coupling for simplicity. Kauko et al. (2017) estimated that up to 0.8 cm of bottom ice was melted by ice algae and other particles in a week. Zeebe et al. (1996) showed that the amount of melting is highly sensitive to the amount of solar radiation, the ice algal biomass, and its vertical extent. Furthermore, Lavoie et al. (2005) quantified the impact on ice algal loss (release into the water column), and found that the loss is highly sensitive to the amount of solar radiation and ice algal biomass, but not so to the photosynthetic efficiency of ice algae. Future model studies can assess the impact on under-ice blooms, and eventually the combined effect of both the shading and melting by ice algae.