Long-term mesoscale variability of modelled sea-ice primary production in the northern Baltic Sea

The Baltic Sea features seasonal sea ice structurally similar to that of the polar oceans (Vihma and Haapala 2009). However, the large riverine inputs to the Baltic Sea and limited water exchange with the ocean result in low-salinity seawater and, consequently, low sea-ice brine volumes. Brine salinity together with ice temperature determine the permeable fraction of sea ice, i.e., the habitable space for the biota, which is small in the Baltic Sea (Granskog et al. 2006) compared with the more saline polar oceans. In the northernmost areas just south of the Arctic Circle, the sunlight hours are few throughout most of the ice-covered season, indicating that light is a major limiting factor for sea-ice algal growth (Haecky et al. 1998). During non light-limiting periods, significant correlations between photosynthetic capacity, chlorophyll a (Chl-a) and inorganic phosphorus (P) point to actively growing algal communities and to P availability as a crucial limiting factor of photosynthesis in the lowermost sea-ice layers (Kaartokallio et al., 2007). Kuparinen et al. (2007) compiled several sea-ice datasets of biogeochemical properties measured from the Baltic Sea between 1996 and 2006 and concluded that Baltic sea ice is a highly active ecosystem that may contribute to annual nutrient and carbon (C) cycles, with complex controlling mechanisms not yet fully understood. The sea-ice bottom communities of undeformed ice were considered the most biologically active, with presumed decreasing habitable space from the southernmost to the northernmost basins, due to decreasing sea-ice bulk salinity.

The Gulf of Bothnia is the northernmost arm of the Baltic Sea; it comprises the southern Bothnian Sea and the northern Bothnian Bay basins (Figure 1). Several rivers flow into the Gulf, and a salinity gradient is present from north to south. In the south, the water is the normal brackish water of the Baltic Sea, with salinity about 6 (practical salinity unit, henceforth unitless), while in the Bothnian Bay the water is nearly fresh (Leppäranta and Myrberg, 2009). The Bothnian Bay is approximately 36 800 km² in area, with average and maximum depths of 41 m and 146 m, respectively. The mean water circulation is cyclonic with mesoscale features, including coastal jets (Leppäranta and Myrberg, 2009). During average winters, sea-ice formation starts as early as October, and the sea ice melts completely by early June at the latest. Undeformed sea ice in the Bothnian Bay has low Chl-a concentrations and low bacterial turnover rates, even when high concentrations of allochthonous dissolved organic matter (DOM) have been reported (Stedmon et al., 2007). This suggests that bacterial growth is controlled mostly by the algal production of the most labile fraction of DOM (Kaartokallio, 2001; Meiners et al., 2002; Kaartokallio, 2004; Kaartokallio et al., 2007; Kuparinen et al., 2007).

Here, a new sea-ice development within the BFM (Biogeochemical Flux Model, http://www.bfm-community.eu) community (i.e., open-source) model in offline coupling configuration with NEMO (Nucleus for European Modelling of the Oceans, http://www.nemo-ocean.eu) is described and was applied throughout the Bothnian Bay. Selected output variables (such as sea-ice thickness, h_i) of a NEMO (v. 3.3) implementation in the Baltic Sea and North Sea (NEMO-Nordic, Hordoir et al., 2015) were used as input to a previously validated version of the BFM for sea ice, also applied to the Baltic Sea (BFM-SI, Tedesco et al., 2010), which computed the biogeochemical dynamics. Numerical simulations of sea-ice physical and biogeochemical properties were performed for the period 1991–2007.

We first briefly review the main features of the NEMO-Nordic and BFM-SI. We present the coupling strategy between the models and its assumptions, followed by a description of the model setups. We show the spatial and temporal variability in the simulated inorganic macronutrients (referred to here as “nutrients”), sea-ice algal Chl-a and C contents in the basin. We then determined the controlling factors of sea-ice primary production among light, nutrients and sea-ice permeability, i.e., the factors related to the habitable space that affect the fluxes of organic and inorganic matter at the ice-water interface. Finally, we address the potential consequences of the changes in ecosystem functioning suggested by Eilola et al. (2013) under climate change scenarios of Baltic sea-ice thinning and retreating (Meier et al., 2004).

A Baltic-North Sea configuration (NEMO-Nordic) of the three-dimensional ocean model NEMO v. 3.3. developed by Madec and the NEMO team (2016) coupled with the Louvain-la-Neuve Sea-Ice Model (LIM3) (Vancoppenolle et al., 2012) was used to generate the prescribed ice conditions (such as h_i) used for the BFM-SI. NEMO solves the primitive Navier-Stokes equations and is based on a nonlinear free surface ocean general circulation model (OGCM) on a C-grid. The LIM3 model is a representation of both thermodynamic and dynamic sea-ice processes. The thermodynamics calculate the melting and freezing rates to achieve a local balance of heat and water fluxes. The dynamics calculate the ice velocity by balancing forcing from wind stress, ocean drag and sea-surface tilt with internal ice stresses. The LIM3 also includes an explicit five ice thickness category distribution. Heat, freshwater fluxes and surface stresses are computed from the atmospheric state and modified by the ice model at every fifth time step of the ocean model.

The BFM-SI is a comprehensive stoichiometric biomass-based model directly derived from the pelagic BFM (Vichi et al., 2007, 2015). Previously, the model was used to describe Arctic and Baltic sea-ice biota in one-dimensional studies (Tedesco et al., 2010, 2012; Tedesco and Vichi, 2014; Thomas et al., 2017). Here, the BFM-SI is used for the first time in a basin-scale application. The first implementations of the BFM-SI were devoted to process studies of the sea-ice ecosystem, with special focus on sea-ice algae (Tedesco and Vichi, 2010; Tedesco et al., 2010) and coupling with phytoplankton (Tedesco et al., 2012; Thomas et al., 2017). In an advanced stoichiometrically flexible framework, the model describes the inorganic nutrients nitrate+nitrite (NO_x), ammonium (NH₄), phosphate (PO₄), silicate (SiO₄), two functional groups of algae (i.e., adapted, resembling mostly fast-growing, light-adapted and high-nutrient-demanding diatoms, and survivors, resembling mostly slow-growing, light-acclimated and low-nutrient-tolerant flagellates and others), particulate organic matter (POM), DOM, and gases such as carbon dioxide (CO₂) and oxygen (O₂), for a total of 22 state variables. Two other living functional groups, sea-ice bacteria and fauna, are also represented in the BFM-SI, but are not discussed here, though their contributions to in situ nutrient regeneration processes are recognized. All sea-ice variables are exchanged at the ice-ocean interface with the underlying seawater, while only inorganic nutrients, organic matter and gases have boundary fluxes at the ice-atmosphere interface. Nutrient supply for algal growth comes not only from the mixed layer, but also from snow deposition through brine drainage for surface communities and from in situ regeneration processes. At the ice-ocean interface, the entrapment of dissolved matter follows the same partitioning of salt in sea ice, while the release to the water column is dependent only on the sea-ice melting rate and concentration (as for particulate matter). The entrapment of particulate matter is assumed to be only a function of the seawater concentration, the sea-ice growth rate and the habitable space. The 22 state variables are listed in Table 1, and the interactions between the various functional groups and with the environment are presented in Figure 2. The biogeochemical processes involved in algal net primary production (NPP) include photosynthesis, respiration, exudation and lysis. The nutrient requirements are regulated by optimum, minimum and maximum uptakes, and by multi-nutrient stress. Dynamic Chl:C ratios allow algae to acclimate to less favourable light conditions than their optima (both lower and higher). For a more detailed description of the biogeochemical dynamics of the BFM-SI and for the model equations, the reader is referred to the BFM-SI manual (Tedesco and Vichi, 2010) and to Tedesco et al. (2010).

The spatial resolution of the NEMO-Nordic (Hordoir et al., 2015) is approximately 2 nautical miles, comprising 56 ocean vertical levels from 0 m to 700 m, with a 3-m resolution at the top level, increased layer thickness towards greater depths and 10 sea-ice vertical layers. Numerical simulations in this study were performed and forced for the period 1991–2007 by the Swedish Meteorological and Hydrological Institute (SMHI) reanalysis product EURO4M (Dahlgren et al., 2016; Landelius et al., 2016), with a horizontal resolution of 22 km. The forcing consisted of 3-hourly wind, longwave and shortwave radiation, air temperature, humidity and precipitation readings. Selected averaged physical components were saved daily. Because the NEMO-Nordic prescribed constant sea-ice salinity in the study area, we used the model of Tedesco et al. (2010) to calculate the salinity-dependent variables: bulk sea-ice salinity, brine salinity and volume, given the sea-ice temperature and sea-ice thickness from NEMO. The model outputs used as the boundary conditions for the BFM-SI included ocean temperature and salinity, wind speed, surface sea-ice and under-ice irradiance, sea-ice thickness, growth/melting rates, sea-ice temperature, sea-ice bulk salinity, brine volume and brine salinity. Because the BFM-SI is able to compute the response of the biological community to instantaneous irradiance, daily average irradiance or daylight average with explicit photoperiod, we forced the BFM-SI with daily average irradiance from NEMO and then converted to photosynthetically active radiation (PAR) by applying a factor of 0.5. Following the analyses of Kuparinen et al. (2007), who identified the lowermost permeable layers of undeformed ice in the Bothnian Bay as biologically active, we used the approach of Tedesco et al. (2010) to model the habitable space. This space represents the sea-ice permeable fraction, i.e., that part of the ice that extends from the bottom upwards, with a brine volume larger than 5%, operationally defined as the sea-ice biologically active layer (BAL) by Tedesco et al. (2010).

The BFM-SI was run in coupled configuration with its pelagic counterpart the BFM, thus with an interactive biogeochemical ocean. Due to our interest in the sea-ice season only, we simplified the standard pelagic BFM of Vichi et al. (2007) to include only those functional groups that are also found in the sea-ice environment. These include inorganic nutrients, diatoms and flagellates, bacterioplankton, and single groups of micro- and mesozooplankton (although the latter not allowed to enter the sea ice), as was done previously (Tedesco and Vichi, 2010; Tedesco et al., 2010; 2012; Thomas et al., 2017). The coupling between the pelagic and sea-ice models required no sea-ice initialization and ensured full mass conservation because the initial concentrations in the sea ice were regulated by the fluxes at the sea ice-ocean interface, as detailed in Tedesco and Vichi (2010); Tedesco et al. (2010); Tedesco and Vichi (2014).

The BFM-SI runs were initiated on the first day of ice formation and ended on the day of complete melting at each of the 3106 grid points and for each ice season between 1991 and 2007. The contribution of the surrounding rivers to the overall distribution of nutrients and organic matter was taken into account implicitly, as the pelagic model was initialized every ice season with available wintertime nutrient and Chl-a data measured within the study area (see Figure 1 for the station locations and Table 2 for the spatial and temporal distributions of the measured variables). The monitoring data were taken from the International Council for the Exploration of the Sea (ICES) Dataset on Ocean Hydrography, Copenhagen, 2011. Additional Chl-a data were taken from the Swedish environmental monitoring and the UMF (Ume Marine Sciences Centre) database dBotnia. The nutrient samples were taken mostly during November and December; if both were available we used those from November. The Chl-a measurements were taken mostly during January. The nutrient data in the upper 15 m of the water column were averaged to represent the mean value for the typical mixed-layer depth and then optimally interpolated throughout the area modelled. Due to the limited number of observations of Chl-a winter data, the average value was often used for the entire area.

The observed dates of freezing, breakup and maximum h_i retrieved from three ice stations in the Bothnian Bay (Röyttä, Hailuoto and Märäskär, see Figure 1 for station location) were compared with the NEMO-Nordic results extrapolated from the corresponding grid point for 1991–2007 in Figure 3 and summarized in Table 3. Overall, the model showed favourable agreement with the observations, tending to slightly delay the start and end of the ice seasons and to overestimate the maximum h_i. The standard deviations for all the variables considered, both modelled and observed, showed wide ranges at all three locations.

At Röyttä (Figure 3, left panels), the model tended to advance the start of the ice season by approx. 4.5 days and to delay the end of the season by approx. 8 days, which led to an overall 12-day-longer simulated than observed length of the ice season. The modelled maximum average h_i (0.79 m) was 0.03 m larger than that observed (0.76 m). The derived trend lines (i.e., linear fits over the simulated period) from Röyttä showed varying patterns. Both the modelled and observed freezing dates were consistently delayed over time, with a significant (p-value 0.000) correlation coefficient of 0.8896. However, the observed date of breakup was anticipated over time, while the modelled date showed almost no trend or a slight delay. Finally, there was no clear trend in the observed maximum h_i, while the model showed a sharp decrease in the maximum h_i from 1991 to 2007, the only trend statistically significant for Röyttä (p-value 0.0202).

At Hailuoto (Figure 3, centre panels), the modelled ice season began approx. 10 days later than that observed and ended being delayed by 8 days, leading to only a 1.6-day shorter simulated length than in the observations. The model resulted in overestimation of the maximum h_i by approx. 0.012 m at this station. The trend lines showed a pattern similar to that for Röyttä. The beginning of the ice season was delayed over time (correlation coefficient 0.7437, p-value 0.0015), and there was less evident anticipation of the end of the season (correlation coefficient 0.6772, p-value 0.0040). There was also widening divergence between the non-trend observed in maximum h_i versus the consistently decreasing h_i in the modelled trend, which was also the only statistically significant trend for Hailuoto (p-value 0.0189).

Comparison between the modelled and observed physical properties at Märäskär (Figure 3, right panels) showed patterns very similar to those at Hailuoto. The modelled beginning of the ice season was delayed by about 9 days, and the modelled date of breakup was delayed by about 7 days. The observed versus modelled length of the ice season differed on average by only 1.9 days. At Märäskär, the model also resulted in slight overestimation of the maximum h_i (0.51 m versus the observed 0.45 m). The trend lines showed patterns similar to those at the other stations for the dates of breakup (no trend) and the maximum h_i (no observed trend versus decreasing modelled trend, although the values were not significant, p-value 0.3802). There were, instead, contrasting trends in the freezing dates, which were slightly anticipated, while the model resulted in an almost consistent delay over time. Neither trend line was significant (p-value 0.5966 and 0.1441), but both were significantly correlated (coefficient 0.5207, p-value 0.0386).

The NEMO-Nordic model based on NEMO v. 3.3. includes no representation of snow, which may have a two-fold counteracting effect. Firstly, the snow cover effectively attenuates the solar radiation such that neglecting the snow cover would lead to enhanced transmittance through the ice, eventually advancing the end of the ice seasons, which however was not the case at any of the ice stations we analysed (Figure 3 and Table 3). Secondly, due to the insulation properties of the snow cover, the modelled h_i was expected to result in slight overestimation of the thermodynamic growth of sea ice, which was in fact the case on average, although the overall overestimation was significantly small (3–12 cm, Figure 3 and Table 3).

To account for the potential overestimation of available light, we tested the feasibility of using the NEMO under-ice irradiance as representative of the sea-ice BAL irradiance. To our knowledge, the only published under-ice measurements of irradiance within the study area are those of Haecky and Andersson (1999). The observed h_i ranged between 0.35 m and 0.7 m (Figure 4a), while 0.06–0.13 m of snow were observed between mid-January and late April. As the modelled sea ice was thinner and melted more quickly than the observed ice (Figure 4a), the modelled under-ice irradiance was mostly larger than the irradiance measured (Figure 4b). However, in January during the reported snow cover period, the modelled under-ice irradiance was already smaller than that observed.

We further compared the measured h_i and transmittance obtained by Ehn et al. (2004) during March and April 2000 in a southern coastal ice station with our modelled results by extracting the corresponding model grid point (Figure 4c, d). The h_i observed ranged between 0.28 m in March and 0.1 m in April. During the sampling period, the sea ice was mostly snow-free, with only some snow patches at the surface. However, further analysis revealed that the snow melted during warm periods prior to sampling and percolated through the ice. The transmittance measured in the PAR region averaged about 30% throughout the ice season and in creased to about 70% when the observed h_i was reduced to 0.1 m. Here, the model resulted in underestimation of the observed h_i, and the sea-ice season terminated earlier than that observed (Figure 4c). Consequently, the modelled transmittance was generally higher (Figure 4d), but closer to that observed when the difference between the modelled and observed h_i was also smaller.

We also compared our results with the phytoplankton light saturation indices measured by Rintala et al. (2010). The light saturation index refers to the time when the algal community ceases to be light-limited. The observed sea ice was mostly pack ice, with h_i ranging between 0.1 m and 0.63 m, and was snow-covered, with some snow-ice structural features. We compared our modelled sea-ice BAL irradiance (Figure 4f) with the measured light saturation indices from Rintala et al. (2010). We observed that the index values were always higher than our modelled sea-ice irradiance in March, when sea ice usually melts, indicating that the modelled light availability for the ice was in a realistic range. We observed favourable matching between the measured and modelled h_i values (Figure 4e), highlighting the better performance of the model for pack ice.

Finally, we compared the maximum observed snow thickness at Röyttä (0.31 m), Hailuoto (0.27 m) and Märäskär (0.22 m) with the snow depth at the same sites when the sea ice had reached its maximum h_i (0.08 m, 0.06 m and 0.05 m, respectively, Table 3). We observed that the latter readings were, on average, less than the critical values of snow depth for algal growth, which are between 0.1 m and 0.2 m (Lavoie et al. 2005; Mundy et al., 2007; Lange et al., 2015). It must be noted that the ocean heat fluxes of the Baltic Sea are extremely small, due to the mostly permanent stratification, and thus the heat contribution from below is very limited (Launiainen and Cheng, 1998). This typical feature implies that the sea ice of the Baltic Sea primarily melts from the top and only after the snow cover has already melted. The model results showed that at the time of the bloom initiation, i.e., after maximum h_i, the observed snow depth was already under the critical limiting threshold for algal light limitation. We conclude that although snow was not represented in this first implementation of the NEMO-Nordic, substituting under-ice irradiance for sea-ice BAL irradiance may result in an acceptable approximation for light availability within the sea-ice BAL.

Published measurements of sea-ice biogeochemical properties are limited and scattered worldwide. Large-scale observations, such as those provided by remote sensing, are likewise not available for sea-ice biogeochemical properties. This scarcity of observations challenges testing and validation/verification of sea-ice biogeochemical models (Steiner et al., 2016). One way to overcome this problem is to test and validate/verify models at specific study sites, where high-temporal resolution data are available, and then extend the use of the model to larger areas, comparing the results with the measurements from the study area and within the time period investigated (Arrigo et al., 1997; Deal et al. 2011; Saenz & Arrigo, 2014). Here, we used the latter approach, extending the applicability of the already validated model against weekly observations collected in the Baltic Sea (Tedesco and Vichi, 2010) to the entire Bothnian Bay. We compared our model results with available observations collected within the bay during the simulated period (1991–2007). Seven sea-ice biogeochemical datasets collected in the Bothnian Bay were published within the timeframe of our model simulations. Four datasets (Haecky and Andersson, 1999; Kaartokallio 2001; Meiners et al. 2002; Vähä, 2005) were already compiled in Kuparinen et al. (2007). To these, Steffens et al. (2006); Rintala et al. (2010); Piiparinen et al. (2010) are added here. A summary of the physical, chemical and biological properties of these datasets, which were compared with the model simulations, is given in Table 4. To compare the results from the modelled sea-ice BAL, we took values from the bottommost section of the measured ice cores, which ranged from 1 cm to approximately 20 cm at the maximum. Only in one case did we have no section depth information, so we compared the model results with measurements averaged over the entire ice core. Overall, these observations indicated maximum h_i ranging between 0.3 m and 0.7 m and snow cover ranging from 0.06 m to 0.3 m. The biogeochemical data showed wide ranges of values for all the properties, apparently as a function of ice types (land-fast versus deformed) and location (coastal versus open sea), where the coastal ice is more nutrient-enriched and the deformed ice more productive.

Comparison between the model results and measurements was done by extracting the corresponding modelled variables from the grid point corresponding to the station location (or the closest grid point in same cases; see Table 4 for details) and throughout the measured period. In case more than one station was sampled, the entire domain containing the station locations was extracted. Ranges of values for h_i, nutrients and Chl-a were then compared between the measured and modelled variables (Table 4). In general, there was favourable agreement between the modelled and observed h_i values, although the former tended to fall towards the lower end of the range observed. Similar agreement was found between the observed and modelled PO₄ values, while less agreement was shown between the modelled and measured NO_x and SiO₄ values. However, both the NO_x and SiO₄ values were highly variable in time and space, perhaps reflecting the fact that P is the most limiting nutrient, as typically reported for the Bay (e.g., HELCOM, 2009). The modelled Chl-a values were mostly smaller than those observed, which probably resulted from most of the measurements being made on thicker and more productive coastal sea ice. In fact, the NEMO-Nordic lacks a land-fast ice parameterization and may tend to result in underestimation of h_i in coastal areas, if the model allows sea ice to drift from the coast. None of the monitoring stations used for model initialization was very coastal, which may have led to some underestimation of nutrient availability along the coastline. The only favourable agreement for all of the variables was in the comparison with the measurements presented in Rintala et al. (2010). This case was the only one in which the model results were compared with measurements taken on a spatial scale larger than our model grid. As the timeframe of the observations was very small (i.e., 7 days), we conclude that the model performs better on a limited temporal than spatial scale compared with observations collected at point stations.

The average physical and biogeochemical properties of sea ice in the Bothnian Bay within the period of simulations (1991–2007, Figure 5) are described here, with the aim of analysing their spatial variability within the basin. Henceforth, the term maximum refers to maximum concentration, while the term peak refers to the corresponding day of the year of maximum concentration. Unless otherwise specified, all the results presented refer to the biogeochemical properties within the sea-ice BAL, i.e., the two-dimensional state variables of the biogeochemical model presented in Table 1 were all integrated over the BAL thickness.

There was a clear positive gradient from south to north in the length of the ice season, breakup dates, maximum sea-ice irradiance, maximum h_i and average thickness of the BAL. On average, the ice season in the Bothnian Bay lasted from at least 110 days in the southern parts of the basin to as many as 170 days on the northernmost coast (Figure 5a). The average sea-ice breakup dates were between late March in the south and early May in the northeastern-most part of the bay (Figure 5b), and a similar pattern can be seen in the maximum sea-ice irradiance, which varied between 240 and 360 µE m^–2 s^–1 (Figure 5c). The maximum h_i varied between 0.3 m in the southern basin and 0.9 m in the northernmost basin, with thicker ice towards the coastal areas, indicating ice drifting in the middle of the basin and deformation towards the coast (Figure 5d), as reported previously (e.g., Uotila, 2001). The average thickness of the sea-ice BAL varied from a minimum of 0.02 m in the southern basin to a maximum of 0.075 m in the northern basin, pointing to a rather small amount of habitable space for sea-ice algae (Figure 5e), compared with that of the more saline ice-covered oceans.

The maximum amounts of PO₄ and NO_x (Figure 5g, h, respectively), which are generally found when the maximum h_i is reached and prior to algal uptake and growth, were seen along the coast (up to 0.35 and 21.0 mmol m^–3, respectively). The sea-ice N:P ratios (Figure 5i) indicated that the sea ice was characterized by P-limitation, similar to its parent waters, in contrast to the rest of the Baltic Sea, which is generally considered nitrogen (N)-limited (HELCOM, 2009).

The BFM-SI describes DOM as produced by sea-ice organisms via exudation and lysis of the cytoplasm and POM as produced by organisms via lysis of the structural parts of the cells (Tedesco and Vichi, 2010; Tedesco et al., 2010). The accumulation of dissolved organic carbon (DOC, Figure 5j, up to 75 mg C m^–3) was much higher than that of particulate organic carbon (POC, Figure 5k, up to 3 mg C m^–3). This difference indicated that most of the organic pool was excreted via algal activity and/or from nutrient-stress excretion, rather than released from algae as lysis of the structural part of their cells. The DOC pool represented almost 50% of the average annual NPP (up to 155 mg C m^–3 y^–1, Figure 5l), pointing to a low level of production in the bay. This evaluation was also confirmed by the very low Chl-a values found in the sea ice, with average maximum concentrations of less than 1.2 mg Chl-a m^–3 (Figure 5f). Examination of the spatial variability revealed clearly similar trends in the maps of DOC, Chl-a and NPP. The largest fraction of NPP, as well as the highest Chl-a values, were found along the coast and were more pronounced on the northern coast and in the central basin, especially the southern part. This pattern was also observed in the DOC spatial variability and partly in the POC distribution.

The interannual variability in NPP between 1991 and 2007 is presented in Figure 6. Most of the ice seasons were characterized by low sea-ice NPP, between 60 and 160 mg C m^–3 y^–1, while the productivity of two ice seasons was exceptionally low (1997–1998 and 1998–1999, with maximum concentrations of less than 80 mg C m^–3 y^–1, Figure 6g, h, respectively). Only one ice season was exceptionally more productive, 2005–2006, with NPP rates as high as 720 mg C m^–3 y^–1 (Figure 6o).

In the rest of this section, the potential limiting factors of sea-ice algal growth are investigated by comparing the simulated physical and biogeochemical sea-ice properties between the most productive year (2005–2006) and the least productive year (1998–1999). In particular, the ice season length, habitable space and available light were analysed as potential physical limiting factors, and nutrient concentrations and initial conditions as potential biogeochemical limiting factors.

In 1998–1999, the ice season was long (Figure 7e), with sea ice already forming in early November in the northern parts of the basin and by late December at the latest in the rest of the bay (Figure 7a). The maximum h_i ranged between 0.3 m in the south and 1.2 m in the north (Figure 7g) and was reached between mid-January and late March at the latest (Figure 7i). The latest peak in h_i, reached in the western part of the basin, may have been indicative of some ridging. The sea ice melted completely in early April at the earliest in the southern parts and western coast, and by early May in the remaining areas (Figure 7c). The sea-ice BAL was thicker than average, ranging between 0.03 m and 0.10 m in the northernmost areas (Figure 7k). The maximum sea-ice PAR was reached just prior to ice breakup and was indicative of the maximum light levels received by the sea-ice algal community. In 1998–1999, the BAL maximum PAR showed the highest values in those northern areas where the dates of ice breakup were also the latest of the season (Figure 7m, c, respectively), illustrating that long ice seasons (or late breakups) expose the sea ice to more sunlight hours.

During the ice season in 2005–2006, the sea ice formed significantly later (i.e., between mid-December in the northern bay and mid-January in the southern bay, Figure 7b). The length of the ice season was significantly shorter (by at least 1 month throughout the bay, Figure 7f). However, the ice season ended relatively later, by at least 1–2 weeks in most of the bay (Figure 7d). The sea ice was generally thinner, with maximum h_i ranging between 0.1 m and 0.7 m (Figure 7h), which was reached again between mid-January and late March at the latest (Figure 7j). The sea-ice BAL was on average thinner (0.01–0.05 m, Figure 7l). The maximum sea-ice PAR levels were generally higher (Figure 7n), due both to thinner ice, but also to later sea-ice breakups (i.e., more sunlight hours). Thus, in this season, the habitable space appeared to be more limiting and the light less limiting than in 1998–1999.

The maximum nutrient concentrations in the sea ice were similar in both years (Figure 8a–d), and the highest biomass in 2005–2006 was attributed to sea-ice diatoms (Figure 8f), although there were small contributions from other types of algae as well (Figure 8h). Interestingly, during the most productive year the initial concentration of PO₄ in seawater was lower (Figure 8j) and the concentration of Chl-a significantly higher (Figure 8n). Additionally, the dates of maximum biomass in both years (Figure 8q, r) were generally reached earlier than the dates of maximum h_i (Figure 7i, j). These results led to the hypothesis that some of the biomass accumulated in sea ice during 2005–2006 was driven by physical processes (e.g., entrapment during ice growth) rather than by biological activity. This effect could also be seen in the positive correlation between sea-ice average annual NPP (Figure 6) and initial concentrations of seawater Chl-a (Figure 9) throughout the simulated period.

Investigation of the limiting factors on algal growth is continued here, where the seasonal variability in the physical and biogeochemical properties in the least and most productive years are compared (Figure 10). Some sea-ice PAR levels showed wide variability only in December and January (Figure 10e, f), as some parts of the basin were already ice-covered while others were not. However, all of the other properties showed wide variability throughout the ice seasons in both years, pointing to a very dynamic system, despite its low productivity compared with that of typical model values reported in the Arctic and Antarctic (e.g., Deal et al., 2011; Saenz and Arrigo, 2014, respectively). In both years, the nutrients (Figure 10g–j) accumulated until the sea ice reached its maximum h_i (Figure 10a, b) and were then consumed during algal growth. A modest contribution to the levels of PO₄ and NO_x was likely provided by regeneration towards the very end of the ice season in 2005–2006 (Figure 10h, j, respectively). While the diatom Chl-a and C contents showed similar patterns during the ice season in 2005–2006, they were less evident during the least productive season (Figure 10l, k, respectively). Low Chl-a values were registered throughout the season, while the biomasses of diatoms and surviving algae were more variable, especially early and late in the season, potentially indicating photoacclimation at first and photoinhibition towards the end of the season.

We presented here the first long-term mesoscale numerical simulations of sea-ice biogeochemistry in the Bothnian Bay, northern Baltic Sea, with a sea-ice biogeochemical model of relatively high complexity. The methodology was based on prescribed physical properties computed by a regional ocean model and used as inputs to the biogeochemical model. It is a feasible method for sea-ice areas that do not feature widely significant ridging. This feasibility was shown in our case given the rare occurrence of thicker-than-average sea ice during the simulated period (Figure 11). We showed an application of this methodology, based on the coupling of the NEMO-Nordic model (Hordoir et al., 2015) with the BFM-SI model (Tedesco and Vichi, 2010; Tedesco et al., 2010, 2012).

The modelled spatial and temporal (inter- and intra-annual) variabilities throughout the Bothnian Bay during the ice seasons in 1991–2007 were investigated. The model was highly sensitive to the seawater conditions that were derived from the early winter measurements, highlighting the importance of extensive winter monitoring programmes to better constrain future modelling efforts. The Bothnian Bay is a low-productivity basin (e.g., Meiners et al., 2002), but biologically active throughout the ice season. In addition to light being a major limiting factor for sea-ice algal growth (Haecky et al., 1998), short ice seasons, such as those that characterize the Bothnian Bay, apparently provide less available time for biomass build-up. Early breakups mean fewer sunlight hours, thus even more light limitation, as previously suggested by Tedesco and Vichi (2014) while investigating the ecological consequences of projected shorter ice seasons on Arctic sea-ice microbiology.

In considering a larger portion of the Baltic Sea from the southern Gulf of Finland to the northernmost Bothnian Bay and corresponding decreasing gradients in sea-ice bulk salinities and brine volumes, Kuparinen et al. (2007) speculated that the sea-ice biologically active fraction also decreased along the same gradient. In considering only the Bothnian Bay, our modelling results showed that the sea-ice permeable fraction instead displayed an increasing gradient. This finding indicates that, on a smaller spatial scale than the entire Baltic sea-ice domain, the habitable space in the ice is affected less by the slight variability in sea-ice bulk salinity than by the thickness of the ice.

The ice seasons in 2000–2001, 2001–2002, 2004–2005, and 2005–2006 were characterized by higher biological production prior to the sea ice reaching its maximum h_i (compare results in Figure 11 and 12). This pattern points to the relevance of the initial pelagic conditions for the entire dynamics of the following ice season, as also suggested by e.g. Kaartokallio (2004). During the other ice seasons examined, the highest productivity was reached after the onset of ice melting, as also observed in e.g. Haecky and Andersson (1999), thus indicating active algal growth in the ice. We then hypothesize that seawater conditions prior to ice formation regulate the initial nutrient and biomass that are found in sea ice in the early stages of ice formation. This biological community then maintains its activity throughout the entire ice season, as observed by Kaartokallio (2004) in the southern Gulf of Finland, and functions as a seeding population for the blooming of ice algae once the light no longer limits their growth.

In all, ice-algal production in the Bothnian Bay appears to be limited by several linked factors. In the southern basin, the smaller amount of habitable space, shorter ice seasons, earlier dates of ice breakup and fewer hours of sunlight are more pronounced limiting factors than in the northern basin. Offshore, the nutrients seem to be more limiting and are less available than in the coastal zones, as expected in the study area. Finally, the seawater conditions prior to ice formation vary widely over an interannual basis and strongly determine the following sea-ice biogeochemical dynamics and their potential for biological production.

A change towards milder ice winters has already been detected on longer timescales (Vihma and Haapala, 2009; Uotila et al., 2015) and can partly be seen here for both h_i and length of the ice season (Figure 13a, b, respectively). Analysis of the temporal evolution of sea-ice NPP showed no significant trend over time (Figure 6), but did demonstrate rather wide interannual variability (Figure 13). Both the average annual NPP and yearly cumulative NPP throughout the domain showed no direct correlation with h_i or with length of the ice season (Figure 13a, b, respectively). This result points to the complex responses of the ice-associated ecosystems to more monotonic trends in climatic stressors, such as increasing temperatures, thinning of the ice and changes in sea-ice extent and ice cover duration.

Scenario simulations suggest that the annual maximum sea-ice extent will decrease between 50% and 80% by the end of the century in the Baltic Sea (Meier et al., 2004). Due to its low level of production, sea ice is not a major contributor to the overall production in the Bothnian Bay. However, the disappearance of sea ice will expose the pelagic community to new light regimes and water-mixing patterns, as suggested by Eilola et al. (2013). These authors concluded that spring blooms in the Baltic Sea will start and end as many as 20–30 days earlier in the Bothnian Bay, with more intense production in areas that were previously ice-covered. Our modelling results and previous works (Tedesco et al., 2012; Tedesco and Vichi, 2014) suggest that light regimes together with duration of the ice season and breakup dates are crucial regulating factors for ice algal growth. This complexity may also apply to phytoplankton that will be exposed to different light levels; i.e., earlier blooms may not necessarily be more intense, because the sunlight hours will be fewer when the ice seasons end earlier. This question can be addressed only by models, such as the BFM-SI, that embed a dynamic Chl-a:C ratio and thus can reproduce photoacclimation and photoinhibition processes, as found previously (Tedesco and Vichi, 2010; Tedesco et al., 2010, 2012) and also discussed here.

The model results were satisfactorily compared with the physical and biogeochemical measurements reported from the study area (Table 4). It must be stressed that comparing point data stations with coarse grid numerical simulations of sea-ice biogeochemistry has several shortcomings. Sea-ice physical and biogeochemical properties show spatial variabilities much smaller than our model grid size, which also applies to Baltic sea ice (Granskog et al., 2005; Steffens et al., 2006). Most of the measurements were made on landfast ice, which can be thicker than pack ice and is more productive, due to the vicinity of the coast and higher nutrient supplies. Our model lacks parameterization of landfast ice, which can lead to thinner sea ice than that observed near the coasts in cases where sea ice drifts away from the coast. These results point to the importance of including landfast ice parameterization in regional ocean models, such as the NEMO-Nordic. Including landfast ice would lead to better model performances in terms of h_i, light transmittance, habitable space and resulting biogeochemical properties in coastal areas, where most of the sea-ice physical and biogeochemical observations are collected. The distribution of the winter pelagic monitoring stations was rather uneven (Figure 1) compared with our model grid, and none of the stations was located near the coastline. Not surprisingly, the best agreement between the modelled and observed sea-ice quantities was with those measured by Rintala et al. (2010), the only published sea-ice dataset from non-coastal ice. It used a spatial scale larger than our model resolution, highlighting the importance of consistently comparing model results with observations from a similar scale.

On a global scale, observation-based estimates of annual primary production in sea ice have been as much as 23 g C m^–2 y^–1 in the Canadian Arctic and 38 g C m^–2 y^–1 in the Antarctic sea ice of McMurdo Sound. In comparison to pelagic waters, these levels accounted for 3% to 28% of the annual primary production in ice-covered polar waters (Arrigo et al., 2010). We are not aware of any estimates available on sea-ice annual NPP in the Baltic Sea region besides the one presented here. Based on weekly to biweekly observations at a coastal ice station only, Haecky and Andersson (1999) are the only authors that have estimated annual sea-ice primary production in the Baltic Sea. Their results of approx. 0.11 g C m^–2 y^–1 during a 3.5-month season constituted less than 1% of the total annual primary production and less than 10% of the primary production during the ice-covered season. In addition to h_i and under-ice irradiance (Figure 4a, b), we also compared our NPP modelling results from the same time period and extrapolated from the closest grid point to the location of the ice station of Haecky and Andersson (1999) (Figure 14). The comparison was made after integrating their estimates (mg C m^–2 d^–1) over the observed h_i (m), to compare them in terms of units with our model estimates (mg C m^–2 d^–1) along the BAL thickness (m). Apart from one of the last measurements that showed high rates of production in the ice while our model results showed sea ice already almost completely melted (Figure 4a), we observed the same order of magnitude and generally favourable agreement between the modelled and observed NPP values. This finding supported early estimates by Haecky and Andersson (1999) and indicated that our basin-scale estimates (Figure 6) are reasonable. Due to the scarcity of data and substantial uncertainties, drawing general conclusions from large-scale or mesoscale numerical simulations of sea-ice biogeochemistry must be done with caution and be considered only as best estimates on temporal and spatial scales which are otherwise impossible to observe by sea-ice sampling.

All of the data discussed in this study are included in the tables and figures of this manuscript. The sea-ice biogeochemical model BFM-SI used in this work is part of the BFM package and is written in FORTRAN90. The BFM-SI model code is available for download from the BFM web site (http://www.bfm-community.eu). After registration, go to http://bfm-community.eu/download and download the version including the sea ice. Follow the instructions for setting up, compiling and running the model. All of the simulations presented in this work are stored in local archives and can be made available upon request.

This work was written under the auspices of BEPSII – Biogeochemical Exchange Processes at the Sea-Ice Interfaces – Clic/IASC/SCAR/SOLAS WG. The Swedish Meteorological and Hydrological Institute (SMHI) is acknowledged for providing the NEMO-Nordic model and the atmospheric forcing used in this study. The monitoring data were taken from the International Council for the Exploration of the Sea (ICES) Dataset on Ocean Hydrography, Copenhagen, 2011. Additional Chl-a data were taken from the Swedish environmental monitoring and the UMF (Umeå Marine Sciences Centre) database dBotnia (http://www.umf.umu.se/english/environmental-monitoring-and-environmental-analysis/the-database-dbotnia). The authors acknowledge the use of the BFM system (http://bfm-community.eu/). The authors wish to thank the editor Jody Deming, the associate editor Lisa Miller, and three anonymous reviewers for constructive comments on an earlier version of this work.

LT and EM acknowledge support from the Maj and Tor Nessling Foundation, project number 2012170. LT also knowledges support from the Estonian Science Foundation, contract number MJD62. JH was funded by the Academy of Finland SmartSea project no. 292985.

The authors have no competing interests to declare.

• Contributed to conception and design: LT with contribution from HK, JH

• Contributed to acquisition of data: EM, BWA

• Contributed to analysis and interpretation of data: LT, JH, HK

• Drafted and/or revised the article: LT with contribution from EM, JH, HK

• Approved the submitted version for publication: LT, EM, BWA, JH, HK