CSIB v1 (Canadian Sea-ice Biogeochemistry): a sea-ice biogeochemical model for the NEMO community ocean modelling framework

Process-based numerical models are a useful tool for studying marine ecosystems and associated biogeochemical processes in ice-covered regions where observations are scarce. To this end, CSIB v1 (Canadian Sea-ice Biogeochemistry version 1), a new sea-ice biogeochemical model, has been developed and embedded into the Nucleus for European Modelling of the Ocean (NEMO) modelling system. This model consists of a three-compartment (ice algae, nitrate, and ammonium) sea-ice ecosystem and a two-compartment (dimethylsulfoniopropionate and dimethylsulfide) sea-ice sulfur cycle which are coupled to pelagic ecosystem and sulfur-cycle models at the sea-ice–ocean interface. In addition to biological and chemical sources and sinks, the model simulates the horizontal transport of biogeochemical state variables within sea ice through a one-way coupling to a dynamic-thermodynamic sea-ice model (LIM2; the Louvain-la-Neuve Sea Ice Model version 2). The model results for 1979 (after a decadal spin-up) are presented and compared to observations and previous model studies for a brief discussion on the model performance. Furthermore, this paper provides discussion on technical aspects of implementing the sea-ice biogeochemistry and assesses the model sensitivity to (1) the temporal resolution of the snowfall forcing data, (2) the representation of light penetration through snow, (3) the horizontal transport of sea-ice biogeochemical state variables, and (4) light attenuation by ice algae. The sea-ice biogeochemical model has been developed within the generic framework of NEMO to facilitate its use within different configurations and domains, and can be adapted for use with other NEMO-based sub-models such as LIM3 (the Louvain-la-Neuve Sea Ice Model version 3) and PISCES (Pelagic Interactions Scheme for Carbon and Ecosystem Studies).

1 Introduction

Biogeochemical processes at the sea-ice–ocean interface play an active role in polar marine ecosystems and global cycling of important chemical elements and compounds. For example, microalgae that colonize the base of sea ice in spring can have a strong influence on primary production of phytoplankton through light attenuation, nutrient drawdown, and seeding as well as on secondary production by providing a food source for grazers (Arrigo, 2014). Furthermore, these ecological processes regulate the production and removal of greenhouse gases (e.g., carbon dioxide and nitrous oxide) and other climatically important gases (e.g., dimethylsulfide) in ice-covered regions, and the exchange of these gases with the overlying atmosphere (Vancoppenolle et al., 2013).

However, our current understanding of many of these processes remains limited due to both logistical and technical challenges for field observations (Miller et al., 2015). Process-based numerical models representing sea-ice biogeochemistry can both fill gaps between sparse measurements and aid in the interpretation of these measurements. Furthermore, these models can be used in systematic intercomparisons that can build confidence in our understanding of polar marine science such as has been done for pelagic ecosystem models (e.g., Popova et al., 2012; Jin et al., 2015).

Dupont (2012)Jin et al. (2012)Watanabe et al. (2015)Castellani et al. (2017)

Table 1 Comparison of pan-Arctic 3-D sea-ice biogeochemical model configurations developed in various frameworks. dx: the horizontal resolution; dzo: the vertical resolution of the uppermost water column; dzi: the vertical extent of the biologically active layer at ice base; i₀ (snow surface): the fraction of incoming shortwave radiation that penetrates through the snow surface; Shading: attenuation of light by ice algae; Runoff: river discharge of nitrate.

^∗ Confirmed through personal communication with the lead author (2018).

Although considerable effort has been invested in developing process-based numerical models for sea-ice biogeochemistry over the last 3 decades following the pioneering work of Arrigo et al. (1991), most of these models were applied in one-dimensional (1-D) frameworks, and the results are therefore limited to particular locations (see Vancoppenolle and Tedesco, 2016). Only a few of these models have been applied in three-dimensional (3-D) framework coupled to either a regional or global sea-ice–ocean general circulation model (see Table 1 for a list of 3-D model configurations developed for pan-Arctic studies). More efforts toward developing such 3-D sea-ice biogeochemical models are needed to better understand the large-scale variability in biogeochemical processes within sea ice and their role in underlying pelagic and benthic ecosystems.

In this study, we present CSIB v1 (Canadian Sea-ice Biogeochemistry version 1), a new sea-ice biogeochemical model implemented into the Nucleus for European Modelling of the Ocean (NEMO), a state-of-the-art modelling framework for oceanographic research (https://www.nemo-ocean.eu/, last access: 13 May 2019). To the best of our knowledge, Tedesco et al. (2017) is the only previous study in which a sea-ice biogeochemical model has been coupled to NEMO. However, the coupling is done in an offline configuration in that study. An important advance of the present study is that the model is written within the NEMO code to allow in-line coupling (i.e., physical dynamics are computed simultaneously with biogeochemistry) and the computation of horizontal transport of sea-ice biogeochemical state variables associated with sea-ice drift. These implementations allow for more realistic simulation of sea-ice biogeochemistry and intercomparison of process-based numerical ice algae models. The main objectives of the present study are to describe the development of the coupled model in a pan-Arctic configuration (Sect. 2), present the basic features of the simulation (Sect. 3), and assess the model sensitivity to modifications of parameters and parameterizations (Sect. 4). Key findings of the present study are summarized in Sect. 5. We note that this study is intended as a model description paper, and the analysis focuses on results for the year 1979, corresponding to the end of a decadal model spin-up. The analysis of the simulation beyond 1979, in which more observational data are available for evaluation (Hayashida, 2018a), is planned to be published as a journal article separately.

2 Model description and setup

The fundamental constituents of NEMO are the following three sub-models: ocean physics, sea-ice physics, and ocean biogeochemistry. In the present study, we adopted version 3.4 of NEMO (NEMO v3.4; Madec, 2008) and developed within it an additional sub-model, sea-ice biogeochemistry. Technical details on the code structure of the model developed in this study are provided in Appendix A for those who are interested in using the newly added sea-ice biogeochemical model.

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

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

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.2.2 Correction to the fractionation of under-ice PAR

The second modification to CanOE was made to the PAR fraction of incident solar radiation. PAR is the shortwave radiation in the 400–700 nm wavelength range, which is available for photosynthesis. In CanOE, PAR is 43 % of the downwelling shortwave radiation reaching the sea surface, a well established estimate for PAR in open water (e.g., Morel, 1988). However, this assumption underestimates PAR reaching the sea surface under sea ice. The shortwave radiation penetrating through snow and ice is almost entirely PAR, as radiation outside of the 400–700 nm range is absorbed by the snow and ice (e.g., Zeebe et al., 1996). Thus, we have set the fraction of the downwelling shortwave radiation to unity (instead of 43 %) when computing the sea-surface PAR under sea ice.

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

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial X}{\partial t}}=-\mathrm{\nabla }\cdot \left(\overrightarrow{U}X\right)+D{\mathrm{\nabla }}^{\mathrm{2}}X+\mathrm{SMS}\left(X\right),\end{array}$$

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.3.1 Horizontal transport

Horizontal transport of sea-ice biogeochemical state variables is computed simultaneously and in the same way as the sea-ice physical properties of LIM2 (i.e., snow and sea-ice volume, heat content, and areal coverage). Specifically, it is done by solving the advection–diffusion equation for sea ice. Advection (by which we refer to the transport associated with the resolved motion of sea ice) is computed using the scheme of Prather (1986). Diffusion, on the other hand, represents transport by unresolved motions (random component of sea-ice motion analogous to turbulence in fluids; Thorndike, 1986; Rampal et al., 2009, 2016), and is often tuned to improve numerical stability. Diffusion is computed within the ice pack by evaluating the second-order diffusive operator using the Crank–Nicholson scheme (Crank and Nicolson, 1996), while it is set to zero at the ice edge. The horizontal diffusion coefficient (D) is set to 5 m² s⁻¹, as suggested by Vancoppenolle et al. (2012). The impacts of horizontal transport of sea ice on modelled biogeochemical state variables are discussed in Sect. 4.3.

2.3.2 New ice formation

The bottom 3 cm of newly formed ice is assumed to contain the same concentrations of biogeochemical state variables as those in the underlying water column. Thus, the concentration (X) of any sea-ice biogeochemical state variable is updated as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }X={\displaystyle \frac{{\mathrm{SIC}}_{t-\mathrm{1}}}{{\mathrm{SIC}}_t}}{X}^{\ast }+{\displaystyle \frac{{\mathrm{SIC}}_t-{\mathrm{SIC}}_{t-\mathrm{1}}}{{\mathrm{SIC}}_t}}{X}_{\mathrm{ui}},\end{array}$$

where SIC_t−1 and SIC_t respectively denote the sea-ice concentrations in the previous and current time step. X^∗ denotes the concentration of the sea-ice biogeochemical state variables (meltwater equivalent) after the computation of advection and diffusion but prior to the computation of biological and chemical sources and sinks. X_ui denotes the concentration of the biogeochemical state variable in the uppermost layer of the water column under the ice. Equation (3) neglects the density difference between sea ice and seawater, and therefore violates mass conservation. However, this simplification has a negligible effect on ocean biogeochemistry given the relatively thin sea-ice biogeochemical layer (Hayashida et al., 2017). For ice algae only, a minimum biomass is set at 10 mmol C m⁻³ in order to mimic reasonable overwintering biomass. This threshold is derived based on the observed range of ice algal biomass in young sea ice (Garrison et al., 1983) and by assuming a fixed carbon-to-chlorophyll ice algal cell quota (Mortenson et al., 2017).

2.3.3 Biological and chemical sources and sinks

The biological and chemical processes represented as sources and sinks of the sea-ice biogeochemical state variables are described in detail in Mortenson et al. (2017) and Hayashida et al. (2017). For the ecosystem component, these processes include photosynthesis, mortality, and remineralization of dead organic matter. The growth rate of ice algae is dependent on ambient temperature (of underlying seawater), PAR, and nutrient concentrations (nitrate and ammonium). Note that the growth rate dependence on ice melt considered in Mortenson et al. (2017) has been neglected in the present study because (1) our preliminary results indicated that ice algal blooms were generally insensitive to it, (2) the parameterization for ice melt limitation was applied for a specific location and might not be appropriate for other locations, and (3) the parameterization lacks observational evidence.

In addition to the computation of biological and chemical sources and sinks, processes relevant to the ice–ocean fluxes are computed, including (1) turbulent molecular diffusive exchange of nutrients, (2) release of all state variables into the water column due to basal ablation, and (3) flushing of these variables by flow of water through the ice from rainfall and surface melting (including flooding due to negative freeboard). For process (1), the effects of turbulence are approximated by parameterizing the molecular sublayer as a function of friction velocity, and molecular diffusion is calculated using the observed diffusion coefficient of dissolved silica measured in seawater at 2 ^∘C (Rebreanu et al., 2008).

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

Figure 4 Comparison 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.

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 2 List of model experiments.

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

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 3 List of selected model parameters in the NEMO namelists.

2.5.3 Output

The output of the model experiments was saved as annual means for the first 10 years (1969–1978) and 5-day means for the final year (1979). Ice (snow) volume was defined as the sum of the product of grid-cell-mean ice thickness (snow depth) and the grid-cell area. Ice extent was defined as the areal sum of all grid cells with an ice concentration of at least 0.15. Primary productivity of ice algae and phytoplankton was quantified in terms of depth-integrated (bottom 3 cm of sea ice and upper 90 m of water column, respectively) net primary productivity (NPP). Ice algal NPP is assumed to equal the growth term in the model equation (Mortenson et al., 2017), as the specific growth rate associated with that term is derived from Eppley (1972). This rate is a measure of particulate production, which is considered to provide values closer to NPP than gross primary productivity (GPP) (e.g., Sakshaug et al., 1997; Hashimoto et al., 2005). Thus, the loss due to respiration is implicitly included in the growth term in the model equation for ice algae. On the other hand, CanOE has an explicit representation of respirational loss, and so phytoplankton NPP is defined as the growth minus the respiratory cost of biosynthesis (Christian et al., 2019; see Appendix A3 of Hayashida, 2018a). For the analysis of NPP and PAR only, any grid cell whose ice concentration is 0.15 or greater was considered “under-ice” following Zhang et al. (2010). To investigate the interannual variability in pan-Arctic primary productivity, the ice algal NPP, phytoplankton NPP, and under-ice NPP were integrated annually and horizontally to derive respective pan-Arctic annual quantities. The term pan-Arctic is defined here as the region north of the Arctic Circle (66.5^∘ N). The pan-Arctic mean refers to an area-weighted average over the region north of the Arctic Circle. This areal restriction allows for a consistent comparison to some previous studies (e.g., Legendre et al., 1992; Jin et al., 2012).

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 Reference simulation (EXP0)

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

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

3.1.2 Spatial variability

Figure 7 shows the spatial variability in modelled March- and September-mean ice thickness fields in EXP0 and PIOMAS. The extent of modelled Arctic sea ice can be inferred from the location of the ice edge, defined here as the 0.15 contour of ice concentration (Fig. 7a, b, d, e). Overall, the locations of the ice edge within our model domain are similar between EXP0 and PIOMAS for both March and September. Beyond the model domain, the ice coverage in March extends to Hudson Bay and the Sea of Okhotsk in PIOMAS (Fig. 7b). The March-mean ice thickness distribution in EXP0 includes a band of >5 m thick ice along the coast of the Canadian Arctic Archipelago extending east to Greenland, and a region of relatively thick ice (∼4 m) in the Arctic Basin north of the East Siberian Sea (Fig. 7a). The band is also present in PIOMAS, although it is restricted to the north of Greenland (Fig. 7b). The thick-ice region in the Arctic Basin north of the East Siberian Sea, on the other hand, is absent in PIOMAS. Besides these particular regions, EXP0 generally simulates thicker ice than PIOMAS in the Greenland Sea and various shelf regions (Fig. 7c). On the other hand, EXP0 simulates thinner ice than PIOMAS on the Canadian Polar Shelf and in the Chukchi Sea, the Barents Sea, the Kara Sea, and an area near the North Pole (Fig. 7c). Overall, the mean absolute difference in the ice thickness distribution over the NAA domain is 0.43 m (30 %). Note that the difference is still calculated even if the ice is absent by considering a thickness of 0 m. Also, note that even though PIOMAS assimilates data, it is still a model product, and therefore the difference is not a definite measure of accuracy.

In September, the most notable features in the ice thickness distribution are the presence of thinner (<2 m) ice in an area near the North Pole in EXP0 (Fig. 7d) and thicker (>5 m) ice along the coast of Siberia in PIOMAS (Fig. 7e). While the latter feature is not discussed in any of the literature on PIOMAS, it seems unrealistic considering that it is thicker in September than in March and it is thicker than the multi-year ice present along the band north of the Canadian Arctic Archipelago and Greenland. Both of these features constitute negative ice thickness anomalies in the model relative to PIOMAS (Fig. 7f). The difference is also negative and large (∼3 m) on the Canadian Polar Shelf; this could be due to the fact that the horizontal resolution of PIOMAS (∼22 km; Zhang et al., 2010) is too coarse to resolve the circulation through these relatively narrow channels, resulting in the simulation of first-year ice that is too thick in this region at this particular time of the year. The ice thickness is greater in EXP0 than in PIOMAS in the Arctic Basin, part of the East Siberian Sea, and the Laptev Sea as well as along the eastern coast of Greenland. The mean absolute difference is 0.31 m (38 %), similar to the March comparison.

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

3.2.2 Spatial variability

Figure 8 shows the spatial distribution of annual-mean snow depth, surface seawater nitrate concentration, ice algal NPP, and phytoplankton NPP for 1979. The largest values of ice algal annual NPP (>10 g C m⁻² yr⁻¹) are present on the Canadian Polar Shelf and in the coastal regions of Baffin Bay, the Chukchi Sea, the East Siberian Sea, and the Kara Sea (Fig. 8c). All of these regions have relatively thin snow (less than 0.1 m; Fig. 8a), demonstrating the control of light on ice algal growth. In contrast, the nutrient control on ice algal production is less pronounced; although high ice algal NPP usually coincides with high surface seawater nitrate, it is also present in a few areas where the nitrate levels are relatively low (Baffin Bay and Chukchi Sea; Fig. 8b). Overall, ice algal production is mostly confined to shelf regions (excluding the Barents Sea), consistent with previous model studies (Deal et al., 2011; Dupont, 2012; Jin et al., 2012, 2018).

There are a few noteworthy similarities and differences in the spatial variability in modelled ice algal annual production between the present study and previous model studies. All studies show a moderate-to-high level of ice algal production in Baffin Bay. In contrast, disagreement in the ice algal production is found along the eastern coast of Greenland and in the Bering Sea; the values along the eastern coast of Greenland are moderate (5–10 g C m⁻² yr⁻¹) in Deal et al. (2011) and Jin et al. (2012), while they are low (less than 5 g C m⁻² yr⁻¹) in Dupont (2012), Jin et al. (2018), and the present study. Similarly, although the Bering Sea is identified as a region of high ice algal production by Deal et al. (2011), Jin et al. (2012), and Jin et al. (2018), Dupont (2012) and the present study simulate low ice algal production in this region. A possible explanation for the lower ice algal production in this region in the latter studies is an insufficient nutrient supply from the Pacific boundary as discussed in Dupont (2012). Lastly, the recent study by Jin et al. (2018) finds the Sea of Okhotsk to be a region of elevated ice algal annual production, which we are unable to assess in the present study due to the limited model domain. We also note the difference in the temporal coverage of simulations among these studies, which can explain some of the differences in these results.

The modelled phytoplankton annual NPP is high (>100 g C m⁻² yr⁻¹) in the Atlantic and the Pacific sectors with little-to-no ice cover, moderate (50–100 g C m⁻² yr⁻¹) in the shelf seas along the North American and the Eurasian continents, and low (<50 g C m⁻² yr⁻¹) in the interior of the Arctic Ocean (Fig. 8d). These findings are in quantitative agreement with the results of five different models and satellite-based estimates (Fig. 1 of Popova et al., 2012).

3.2.3 Interannual variability

The modelled pan-Arctic ice algal annual NPP in EXP0 ranges from 10.5 to 18.2 Tg C yr⁻¹ for the period 1970–1979, excluding the initial spin-up year. While this value is on the lower end of the range of observation-based NPP estimates (9–73 Tg C yr⁻¹; Legendre et al., 1992), it is close to the decadal mean of the annual NPP (10.1 Tg C yr⁻¹ for 1998–2007) simulated by Jin et al. (2012). The pan-Arctic estimates by Legendre et al. (1992) are quite speculative as they are based on the integration of a single ice algal production value over a specified ice area (discussed in detail in Deal et al., 2011). The close agreement between the two model-based estimates suggests that the lower end of the observation-based estimates is more plausible than their upper end. Although the upper end accounts for contribution from mat and strand communities that are not represented in our model (representing them would require additional state variable(s) and/or new formulation with more parameters), their contribution to the pan-Arctic production should be small as their spatial distribution is generally localized (e.g., Assmy et al., 2013). Direct comparisons with the results of Deal et al. (2011) and Dupont (2012) are not possible because the reported values in those studies include contributions from below the Arctic Circle. The modelled pan-Arctic phytoplankton annual NPP in EXP0 ranges from 378 to 465 Tg C yr⁻¹, which is in line with the observation-based estimate (>329 Tg C yr⁻¹; “Total High Arctic”) of Sakshaug (2004) and the satellite-based estimate (419 Tg C yr⁻¹ for 1998–2006) of Pabi et al. (2008), but below the model-based estimate (627 Tg C yr⁻¹ for 1998–2006) of Jin et al. (2012).

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

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

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 Sensitivity experiments (EXP1–5)

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

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

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

Figure 14 Model 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.

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

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

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.

5 Conclusions

In the present study, we have developed a sea-ice biogeochemistry model which is coupled to NEMO. A number of modifications to the sea-ice physical model used in the standard distribution of NEMO (LIM2), to the ocean biogeochemical model (CanOE), and to the existing pan-Arctic configuration (NAA) were necessary to properly simulate the physical and biogeochemical processes in ice-covered regions. Results of the reference simulation (EXP0) were discussed and compared with previous studies, with a focus on the year 1979; more thorough evaluation of the model performance over the recent decades is planned for future studies. Adopting a high vertical resolution in the upper water column was found to be necessary to properly represent the effects of a meltwater lens on surface nutrients and the formation of a subsurface chlorophyll maximum. Furthermore, the vertical resolution was shown to have an effect on the magnitude of the modelled surface seawater DMS concentration (∼10 % annually and up to ∼20 % seasonally), which in turn influences DMS emissions. Results of the sensitivity experiments demonstrated that LIM2 requires high-frequency (daily) snowfall forcing data to simulate realistic snow depth (EXP1 and 2); the assumption of no light penetration through snow in LIM2 is unrealistic for simulating an adequate ice algal bloom (EXP3); horizontal transport of sea ice contributes to an enhancement of the pan-Arctic ice algal annual NPP by 16 % (EXP4); and attenuation of light by ice algae has local influence on under-ice NPP but is negligible when estimating larger-scale quantities (e.g., pan-Arctic under-ice annual NPP) (EXP5). While we believe that these findings would be qualitatively similar in other years, it would be worthwhile to quantify their interannual variability. The modifications to LIM2, CanOE, and NAA adopted in the present study are also applicable to other sub-models and configurations of NEMO (e.g., LIM3, PISCES, ORCA) as the code structures are similar, and therefore can be incorporated into future pan-Arctic biogeochemical studies. The sea-ice biogeochemical model developed in the present study has been embedded into NEMO in a generic way (see Appendix A), and can therefore be easily coupled to the aforementioned sub-models. To our knowledge, such a development has not been done previously within NEMO. Further sensitivity experiments and observational constraints are needed to refine the important parameters (e.g., i₀) for sea-ice biogeochemistry.

Appendix A: Implementation of ocean sulfur cycle and sea-ice biogeochemistry into the NEMO source code

Figure A1 shows the structure of the NEMO v3.4 source-code directory (NEMO), which includes the following subdirectories (sub-models): OPA_SRC (OPA), LIM_SRC_2 (LIM2), and TOP_SRC (ocean biogeochemistry). The directory TOP_SRC contains two subdirectories: PISCES and MY_TRC. In this study, the directory PISCES contains the source code of CanOE, as CanOE has been developed using the code structure of the PISCES ocean biogeochemical model. The other directory, MY_TRC, consists of a list of generic modules that can be modified by end users to add their own biogeochemical models; we introduced an ocean sulfur cycle and sea-ice biogeochemistry into this interface. Furthermore, we modified a few modules in the directories LIM_SRC_2 and PISCES for the implementation of sea-ice biogeochemistry into the NEMO modelling system (Table A1).

Numerically, the tendencies for the sea-ice biogeochemical state variables are computed at each time step as follows: first, the concentrations of all state variables from the previous time step are transferred from the module trcsms_my_trc.F90 to the module limtrp_2.F90 to compute the advective and diffusive tendencies. The updated concentrations are transferred back to the module trcsms_my_trc.F90 within which the biological and chemical sources and sinks as well as the ice–ocean fluxes of these state variables are computed.

In NEMO, user-specific modules built within MY_TRC are designed to be activated by defining the C preprocessor (CPP) key key_my_trc. As such, we assigned CPP keys for each component of the newly developed models, which can be activated as needed (Table A2).

Table A1 A list of NEMO modules modified to add ocean sulfur cycle and sea-ice biogeochemistry.

Figure A1 File tree diagram of the OPA–LIM2–CanOE configuration of NEMO v3.4. The modules listed in the diagram (^*.F90) have been modified in order to implement ocean sulfur cycle and sea-ice biogeochemistry into the present configuration.

Table A2 A list of CPP keys created in the present study.

Appendix B: Interannual variability during spin-up

The annual-mean time series of modelled snow and ice volumes, ice extent, seawater nitrate, and ice algal and phytoplankton biomass over the 11 years of EXP0 are shown in Fig. B1. This time period can be considered sufficiently long enough to spin up some of these quantities, while others may require additional time to spin-up. However, none of these quantities reach a steady state in the current setup as the model was driven by interannual surface and lateral boundary conditions. The aim of the present analysis is to examine the temporal variability starting from the initial year and compare with findings of previous model studies. Presenting the results from the initial year is often neglected in the literature, but can be useful for future studies.

The annual-mean modelled snow volume stabilizes around 0.8×10³ km³ after an initial drop of about 0.1×10³ km³ from year 1969 to 1970 (Fig. B1a), indicating a spin-up period of a year or so. In contrast, the annual-mean modelled ice volume variations show an initial reduction during 1969–1971 followed by an overall increase during 1973–1979. The relatively short duration of this simulation does not allow us to distinguish between trends and slow interannual variability, so we cannot determine if the ice volume has spun up based solely on this analysis; this will be addressed in a follow up study. A previous pan-Arctic regional model study of Watanabe (2013) shows a spin-up period of 10 years for modelled ice volume based on a simulation using a fixed annual cycle of atmospheric forcing and restoring of temperature and salinity.

Figure B1 Time series of annual-mean modelled (a) snow and ice volumes, (b) ice extent and depth-integrated (90 m) seawater nitrate concentration, and (c) depth-integrated (3 cm) ice algal NPP and depth-integrated (90 m) phytoplankton NPP in EXP0. The depth-integrated quantities represent averages over the entire model domain.

Modelled ice extent shows a decrease in the first 6 years followed by a stabilization in the last 5 years, suggesting that this quantity spun up at year 1975 (Fig. B1b). This spin-up time is similar to that found in the pan-Arctic model study of Jin et al. (2012), in which their modelled ice area and extent became comparable to the observations after the first 6 years of simulation.

Annual-mean modelled seawater nitrate concentration integrated over the upper 90 m of the water column shows both increases and decreases during the 11 years (Fig. B1b), although the size of the fluctuation (∼20 mmol N m⁻²) is small relative to its mean state (∼490 mmol N m⁻²). Similarly to ice volume, a longer simulation would be needed to distinguish between trends and interannual variability in the modelled nitrate concentration. A previous pan-Arctic model study of Dupont (2012) indicated a spin-up period of at least a decade for nitrate in the upper 100 m water column for the model domain he considered. The modelled primary producers (ice algae and phytoplankton) appear to have spun up within a year of the model simulation, as their annual primary production fluctuates around a steady mean following the first year (Fig. B1c).

Appendix C: Subsurface chlorophyll a maxima

The model simulates subsurface chlorophyll a maxima (SCMs) at various locations in the Arctic. Figure C1 shows an example for the Chukchi Sea where the SCM is formed at depth in mid-June. This SCM depth (>20 m) is comparable to observations (e.g., Brown et al., 2015). A detailed discussion on SCM simulated by CanOE will be provided in a future study (Steiner et al., 2019).

Figure C1 Vertical time series of modelled 5 d mean chlorophyll a concentration in the upper 50 m of the water column at the grid cell corresponding to 170^∘ W and 70^∘ N in the Chukchi Sea.

Appendix D: Parameter tuning to increase snow accumulation

Here we demonstrate that the lack of snow accumulation using the daily climatology for snowfall dataset can be resolved by tuning the model parameter nn_fsbc, which defines the frequency of the computation of surface boundary conditions and sea-ice physics relative to that of ocean dynamics. Figure D1 compares the annual-mean modelled snow depths for year 1970 of EXP2 with those of the simulations that varied nn_fsbc from the default value of 1 (i.e., the time step for surface boundary condition and sea-ice physics is identical to the ocean time step) to 5 and 10 (i.e., surface boundary condition and sea-ice physics are computed at every 5 and 10 ocean time steps, respectively). We find that setting nn_fsbc to 5 or 10 increased the modelled snow depth quite remarkably (Fig. D1). This high sensitivity to the choice of nn_fsbc is somewhat unexpected given that the tested range (1–10 time steps, equivalent to 20–200 min) is far less than the temporal resolution of the CORE-II dataset. A more detailed analysis of the model sensitivity to nn_fsbc is outside the scope of this study. Furthermore, we note that the tuning of this parameter without known constraints is quite arbitrary and might have other implications for modelled dynamics. As discussed in the main text, the usage of high-frequency atmospheric forcing dataset is therefore recommended to realistically simulate snow accumulation.

Figure D1 Sensitivity of modelled snow depth to the parameter nn_fsbc, which defines the frequency of the computation of surface boundary conditions and sea-ice physics relative to that of ocean dynamics. Spatial distribution of annual-mean modelled snow depth for 1970 when nn_fsbc is set to (a) 1, (b) 5, and (c) 10.