Research ArticleCLIMATOLOGY

Sea-ice algal phenology in a warmer Arctic

See allHide authors and affiliations

Science Advances  08 May 2019:
Vol. 5, no. 5, eaav4830
DOI: 10.1126/sciadv.aav4830


The Arctic sea-ice decline is among the most emblematic manifestations of climate change and is occurring before we understand its ecological consequences. We investigated future changes in algal productivity combining a biogeochemical model for sympagic algae with sea-ice drivers from an ensemble of 18 CMIP5 climate models. Model projections indicate quasi-linear physical changes along latitudes but markedly nonlinear response of sympagic algae, with distinct latitudinal patterns. While snow cover thinning explains the advancement of algal blooms below 66°N, narrowing of the biological time windows yields small changes in the 66°N to 74°N band, and shifting of the ice seasons toward more favorable photoperiods drives the increase in algal production above 74°N. These diverse latitudinal responses indicate that the impact of declining sea ice on Arctic sympagic production is both large and complex, with consequent trophic and phenological cascades expected in the rest of the food web.


Sea-ice algae form a large fraction of sea-ice (sympagic) biomass (1). Before the onset of the seasonal phytoplankton bloom, sea ice provides critical habitat for marine algae and the upper trophic levels for whom sea-ice algae are their sole food-source during this period. Together with sub-ice phytoplankton, sea-ice algae represent the foundation of ecological interactions in the sea-ice biome (2). Algal growth in sea ice is limited by the strong seasonality of light, nutrients, and brine connectivity and, in first-year ice, by the relatively narrow time window during which sea ice exists. At high latitudes, light conditions in the sea-ice habitat are subject to extreme seasonal changes. Light availability is largely regulated by the photoperiod, which depends on latitude and time of the year, and by the albedo and light attenuation of different sea-ice surfaces (e.g., snow-covered ice, bare ice, and ponded ice). Nutrients in sea ice are replenished through brine connectivity with the underlying water column and through local regeneration processes. Brine connectivity depends on brine permeability, a function of both sea-ice temperature and salinity. The warmer and less salty the ice, the greater its permeability and thus its habitable space.

We expect future higher sea-ice surface temperature (3) and thus warmer and more permeable sea ice. However, steeper temperature gradients through thinner ice/snow, together with the expected general freshening of the Arctic Ocean (4, 5), might have a counterbalancing effect on sea-ice permeability. Large uncertainties in future nutrient dynamics (i.e., how mixed layer depth might change and be sustained) limit our understanding of the future Arctic Ocean ecosystem (6, 7).

In sea ice, the phenology of the narrow time window for biological growth is constrained by the incidence of suitable light conditions and by the length of the ice season (1, 810). While sea-ice freeze-up (breakup) is overall expected to occur later (earlier) in the season at all latitudes (3), there are more uncertainties in the future timing of sympagic algal blooms (1113), with snow changes considered to be the single most important predictor of future sea-ice algal phenology (14). A few conceptual and qualitative models have been proposed to describe the anticipated response of the seasonal ice zone biota to a warmer climate (2, 1519), and some more quantitative works have used climate model projections and earth system models to analyze the effects of a fresher, ice-free Arctic on phytoplankton primary production (6, 7). Melt onset trends have shown considerable advances on pan-Arctic scale with marked regional differences and the southerly regions melting earlier (20). The strongest trend since 1979 has occurred in the Barents Sea and East Greenland Sea (25 and 30 days, respectively), while less rapid trend is observed in the Baffin Bay, Kara Sea, and Hudson Bay (22, 15, and 10 days, respectively), and the weakest trend is found in the Central Arctic Sea (5 days). There is increasing consensus that these advances in melt onset, together with thinner sea ice and increased stratification, have led to conditions generally more favorable for sub-ice pelagic blooms in the last 30 years (21), suggesting that the advancement of the phytoplankton bloom peak has made the annual season of primary productivity overall longer as well (22).

However, less is known of the conditions that lead to sea-ice algal blooms, which are also dependent on other factors such as brine connectivity. Changes in sea-ice nutrient concentrations are expected to affect the overall amplitude (8, 9, 23) but are less likely to influence the timing of the growth of sea-ice biota (14, 24). Thus, sympagic and pelagic dynamics might not be directly related (25), and the onset of seasonal primary production in the Arctic marine food web may occur in subsequent waves made of sub-ice and open-water phytoplankton blooms (22) on top of sea-ice algal blooms. These pulses of primary production trigger higher trophic levels, from zooplankton to fish to top predators such as whales and polar bears. These phenological cascades are expected to be stronger in the Arctic Ocean where there is tight coupling between abiotic conditions and the timing of primary production and where the food chains are short (22).

Climate model projections are the only available tools that can be used to quantitatively assess potential long-term changes in sea-ice biogeochemical dynamics in a warmer Arctic. Here, we use a near-mechanistic modeling approach that combines sea-ice drivers obtained from an ensemble of 18 CMIP5 [Climate Model Intercomparison Project Phase 5 (26)] climate models with a state-of-the-art sea-ice biogeochemical model (24) to provide a quantitative overview of the effect of different physical drivers on the potential phenological changes in pan-Arctic first-year ice primary production. We focus our study on first-year ice, which is becoming the dominant type of sea ice in the Arctic Ocean (3, 20). In this work, we consider the Arctic first-year ice as a single pan-Arctic physico-biogeochemical unit. This is motivated by the inherent constraints of climate projections and by the sparse availability of recurring studies on seasonal sea-ice biogeochemistry. We seek to elucidate the first-order impacts of large-scale, climate-driven ice changes on sea-ice algal production rather than considering the more local and subregional components, e.g., as done in (27) to explain the coupling between sea ice and tundra in the Svalbard Archipelago since it would require observational efforts that are not yet available on a pan-Arctic scale.


Changes in sea-ice physical drivers of biological responses

Daily and monthly simulations of physical properties of sea ice were used to derive empirical probability density functions of sea-ice drivers in latitudinal bands of 2°, based on 18 CMIP5 climate models (table S1) for recent historical (1961–2005) and business-as-usual representative concentration pathways [RCP8.5, 2061–2100 (28)] scenarios (see Materials and Methods for details). This method is superior to the use of zonal averages when dealing with variable events such as the day of freeze-up and breakup of sea ice, and it has been successfully used in climate model analyses (29, 30). It also allows us to derive statistical descriptors of the extremes of the distribution, without any assumption on the shape of the distribution itself. Models were selected according to the availability of high-frequency daily data at the time of analysis. Sea-ice drivers considered were date of freeze-up, date of breakup, maximum sea-ice thickness, first-year ice extent, first month of snow cover, last month of snow cover, maximum snow depth, and minimum near-surface atmospheric temperature, as well as the dates at which the maxima and minima occurred. The climate data were obtained from the same CMIP5 dataset used in earlier works on future Arctic changes under the RCP8.5 scenario [e.g., (31)]. However, here, we considered minimum air temperature rather than mean air temperature, focused on first-year ice only rather than on all Arctic sea ice, and explored daily data of sea-ice drivers rather than monthly averages, providing an original contribution to the analysis of CMIP5 model outputs.

We initially verified the ensemble means of first-year ice freeze-up and breakup days and the area extent for the period 1979–2005 from the 18 single distributions against remote sensing data (Fig. 1). The observations of sea-ice concentration for current climate conditions were obtained from passive microwave data (32). The sea-ice drivers of phenology (freeze-ups and breakups) from the climate model ensemble distributions compare well with the distributions computed from the satellite data (Fig. 1, A and B). In particular, the simulated medians of the days of freeze-up and breakup match the observations at all latitudes [root mean square differences (RMSDs) are 4 and 11 days, respectively]. The spread of the modeled distribution with respect to observations is similar for the freeze-up (Fig. 1A) but narrower for the modeled breakup, indicating less variability in CMIP5 models (Fig. 1B). Less agreement is found when comparing the medians of observed and modeled sea-ice extent (RMSD = 0.2 × 106 km2; Fig. 1C), likely because of the known biases of CMIP5 models to simulate this sea-ice driver (33, 34). The difference is particularly large in the latitudinal band between 70°N and 78°N. Our statistical approach demonstrates acceptable performances of the climate models, especially in terms of making the sea-ice drivers of phenology usable for the analysis of changes in seasonal sea-ice and their impact on ice algal phenology. The larger spread in sea-ice extent will instead be taken into consideration and discussed when integrating our latitudinal-based results on a pan-Arctic scale.

Fig. 1 Empirical probability density functions of the estimated climatic drivers.

(A) Freeze-up and (B) breakup dates and (C) sea-ice area of Arctic first-year ice as a function of latitude are shown from passive microwave data (32) (shading and green lines) and from the ensemble mean (blue lines) of CMIP5 data (table S1). The gray shading indicates the normalized empirical probability density value for each variable derived from the satellite data over the period 1979–2005. A darker color shows the more frequent days of freeze-up and breakup and the more typical first-year ice area at each latitude. The 50th (solid lines), 10th, and 90th (dashed lines) percentiles from the satellite and model ensemble distributions are shown.

At all latitudes (Fig. 2), the projected physical drivers of sea ice trend toward thinner sea ice (approximately from 0.3 m thinner in the low latitudes to 1.3 m thinner at high latitudes), warmer near-surface minimum atmospheric temperature (up to 15°C warmer above 70°N), later freeze-ups (average, 1.5 months), earlier breakups (average, 1 month), thinner snow depth (0.05 to 0.15 m), and generally shorter ice seasons. The variation of the freeze-up medians is larger in the future scenario, highlighting the intermodel variance in simulated summer sea-ice losses over the 21st century (35). Despite the large projected increase in Arctic precipitation over the century (3), our analysis confirms recent work on the general reduction of snowfall, with rainfall becoming a more dominant form of precipitation in a warmer Arctic (36). Snow depth is halved over first-year ice at high latitudes (Fig. 2G), with a similar response to that observed on multiyear ice (37, 38).

Fig. 2 Latitudinal changes in the sea-ice drivers.

Medians of the empirical probability density functions from the 18 models (thin lines) and ensemble mean (thick lines) in the historical (1961–2005, blue lines) and RCP8.5 scenario (2061–2100, red lines) simulations for (A) date of freeze-up, (B) date of breakup, (C) maximum sea-ice thickness, (D) first-year ice extent, (E) first month of snow cover, (F) last month of snow cover, (G) maximum snow depth, and (H) minimum near-surface atmospheric temperature.

The use of empirical probability density functions shows that the models have good skills in simulating the latitudinal variation of the observable sea-ice drivers (Fig. 1), although it shows the occurrence of first-year ice above 84°N in 12 of the 18 models, which is not seen in the satellite distribution. Another bias that was not previously reported, e.g., by Notz (33), is the underestimation of the areal extent above 70°N (Fig. 2D). This is due to the large variation in the area of first-year ice between 70°N and 80°N, which causes this latitudinal band to be characterized by the largest increase of seasonal sea-ice extent by the end of the century (Fig. 2D).

We used the ensemble means in Fig. 2 to construct representative pan-Arctic sea-ice seasons at 2° latitudinal resolution for both time slices in the form of continuous forcing functions (see Materials and Methods for details). The physical forcings were then used to drive the sea-ice biogeochemical model, simulating the permeability and biogeochemistry of the dynamic bottom permeable sea-ice layer and the development of environmental factors for sea-ice algal growth (24). The biogeochemical model simulated the sympagic system, focusing on one class of sea-ice algae resembling Arctic pennate diatoms, which are the dominant community in bottom sea ice (39) and are characterized by efficient light acclimation while limited by nutrients found in the bottom permeable layer (see Materials and Methods for details). While the initial composition of the sympagic community can be highly heterogeneous [e.g., (40)], mirroring that of the ocean underneath, the blooming phase is usually dominated by pennate diatom species after natural selection of the most adapted taxa to survive and grow when suitable environmental conditions return [e.g., (39)].

Response of sea-ice algae to habitat changes

Assuming unchanged nutrient availability, Arctic sympagic gross primary production (GPP) in first-year ice is projected to overall increase in the future (i.e., from 2.47 to 4.02 Tg C y−1, that is, 63%; Fig. 3 and table S2), indicating an increased uptake of carbon dioxide and thus an overall increase of the carbon biological pump. However, the response of the biological community to quasi-linear latitudinal changes in sea-ice physical drivers of Fig. 2 is markedly nonlinear.

Fig. 3 Latitudinal changes in sea-ice annual GPP.

The percentage changes in the magnitude of sea-ice GPP are shown as a function of latitude and differences in the timing of the bloom peak. This indicator is the difference between the RCP8.5 scenario (2061–2100) and the historical simulations (1961–2005), scaled by the historical value. (A) Percentage change in annual sea-ice GPP per unit area (red bubbles), i.e., integrated temporally along the ice season length. (B) Percentage change in annual sea-ice GPP integrated both temporally along the ice season length and spatially over the ice extent of each simulation period (yellow bubbles). The area of the bubbles is proportional to the percentage change.

Despite similar latitudinal advances in projected breakup dates, the changes in timing of sea-ice algal blooms are different and grouped into two main clusters. Below 66°N, where the breakup date is projected to occur about 25 days earlier in the future climate (Fig. 2B), the bloom peak develops between 45 and 70 days earlier than in the historical climate (Fig. 3A). Above 66°N, the advance of the algal bloom is more in line with the advancing of the breakup dates (Fig. 2B), i.e., between 15 and 40 days earlier in the future projections (Fig. 3A). The number of clusters increases to three when we factor in the relative changes in GPP. For these clusters, low-latitude seasonal sea ice is expected to be up to 150% more productive (i.e., increasing from 0.48 to 1.20 g C m−2 y−1 at 65°N). Smaller and negative variations are expected between 68°N and 74°N (i.e., between −25 and 11% relative change), while the largest increase is expected above 74°N, with up to more than 550% increase in GPP at 77°N (i.e., from 0.05 to 0.33 g C m−2 y−1; Fig. 3A and table S2).

With inclusion of the projected latitudinal changes in ice extent (Fig. 2D), the relative increase in GPP per unit area predicted below 74°N diminishes because of the large loss of seasonal sea ice, while above 74°N, it becomes more important because of first-year ice largely expanding at the expense of multiyear ice (Fig. 3B and table S2). The combination of area-specific production increase and first-year ice expansion is expected to give rise to a sea-ice habitat that is more than 2500% more productive than in the recent climate for algae located northward of 80°N (i.e., from 7.48 to 199.48 Gg C y−1 at 83°N). Thus, the anticipated increase in first-year ice algal production would make this region productively equivalent to what is currently simulated 10° more to the south (table S2), which is a non-negligible change. In light of the CMIP5 model bias in sea-ice extent (particularly large between 70°N and 78°N; Fig. 1C), we anticipate the relative increase in GPP found at these latitudes (Fig. 3B) to be even underestimated.

We searched for mechanisms linking the changes in physical drivers to the phenological differences illustrated in Fig. 3. We specifically investigated alterations in the timing of the peak bloom, driven by latitudinal variation in light availability. The available light in turn is a function of sea-ice and snow thicknesses, length of the biological time window, and photoperiod. We present here one latitude example for each cluster: 58°N to 66°N (lower panels), 66°N to 74°N (mid panels), and 74°N to 84°N (upper panels), representing respectively low-, mid-, and high-latitude pan-Arctic sea-ice ecosystems for historical (left) and RCP8.5 (right) scenario simulations (Fig. 4).

Fig. 4 Drivers and mechanisms of sea-ice algal phenology.

Simulation results at three representative latitudes (61°N, 71°N, and 81°N) for the historical period (1961–2005, left) and RCP8.5 scenario (2061–2100, right), showing sea-ice season length, sea-ice/snow thickness (black), downward solar radiation [from the ERA-Interim climatology (56), used as the reference between historical climate scenario; red], algal chlorophyll anomaly (green), nonpermeable sea ice (gray), permeable and potentially productive sea ice [color gradient from pale blue to green is proportional to normalized GPP (0 to 1), with pale blue representing 0 and green representing 1] Dotted lines highlight changes in snow thickness (black) and biological time window (green). The chlorophyll anomaly (green) is an indicator of algal bloom phenology (22) and is computed as the difference between the concentration and its SD over the sea-ice season (fig. S1).

The first cluster (Fig. 4, 61°N) is typical of regions such as the Bering Sea, the Sea of Okhotsk, the Hudson Bay, or the land-fast ice conditions found in Greenland fjords (41) and the Baltic Sea (42). According to our modeling results, below the Arctic Circle, sea ice is projected to become much thinner and entirely permeable. The thinner snow depth, expected to be less than 0.1 m, leads to a long productive season, with a large and sharp asymmetrical bloom that is characteristic of little light limitation (see also fig. S1C).

The second cluster (Fig. 4, 71°N) includes basins such as the Beaufort Sea and the Mackenzie shelf in the Canadian Arctic (25), the Chukchi Sea, the Baffin Bay, and the Kara Sea, where sea-ice production under current climate conditions is limited by thick ice and snow despite the favorable photoperiod. The latitudes of up to 74°N are currently characterized by one of the largest extents of first-year ice (Fig. 1), and these areas are projected to double by the end of the century (Fig. 2D). In the RCP8.5 scenario simulation at 71°N, both sea ice and snow are expected to thin largely, although not enough to allow sufficient growth within the narrower biological time window due to earlier breakups. The resulting magnitude of sea-ice production is expected to mimic current climate conditions, but the phenology is nonetheless altered toward earlier blooms (Fig. 3 and fig. S1B).

The Barents Sea deserves specific mention because it is mostly located within this cluster’s range of latitudes, but it is characterized by sea-ice conditions more similar to that of the first cluster (Fig. 4, 61°N). In the Barents Sea, we expect an increase in GPP similar to that seen at the lowest latitudes (Fig. 3A) because of significant sea ice and snow thinning (Fig. 4, 61°N). However, we expect this to be modulated by the sharper photoperiod typical of these latitudes. This might reduce the length of the ice algal growing season in a way that is similar to the second cluster (Fig. 4, 71°N).

The third cluster (Fig. 4, 81°N) is typical of first-year ice found in the high Arctic. An example of this cluster is the offshore Lincoln Sea north of Greenland. The area is currently characterized by first-year ice that is generally less productive than other shelf regions, with algal biomass comparable to multiyear ice conditions (43). At these high latitudes, sea-ice blooms are limited by low brine permeability and thick ice and snow. The climate model projections indicate a marked thinning of the sea ice to levels similar to the current-day 71°N cluster (Fig. 2C). At these highest latitudes, the key factor responsible for the largest relative increase in GPP in first-year ice (Fig. 3) is the advance of the sea-ice season toward a more favorable photoperiod, which would give rise to blooms that are expected to be larger than in the adjacent southern cluster (see also fig. S1A).


Dealing with uncertainties and the impacts on the Arctic marine food web

In spite of wide variability, the projections of the climate models consistently show that the thinning and retreat of Arctic sea ice will continue throughout this century, although at different paces in different scenarios (3). Despite the models’ limitations (33), the proposed statistical approach allowed us to extract reliable descriptors of current sea-ice conditions from coarse resolution climate models (Fig. 1), which were used to provide a quantitative description of the impact of future sea-ice projections on the phenology and magnitude of sea-ice GPP (Figs. 3 and 4).

The use of 2° latitudinal bands, driven by the coarseness of climate model data, limits our ability to resolve the exact geographic boundaries of the sea-ice algal response. Our statistical methodology to include daily variability as a way to increase climate model skills in reproducing the seasonal patterns of sea-ice drivers cannot overcome the reduced regional variability, as indicated by the narrower spread of the model ensemble percentiles in Fig. 1. This can only be achieved with models at higher spatial resolution and more detailed sea-ice physical processes. We acknowledge the potential compromises that this approach presents with regard to ignoring variability in the local climatic conditions for the different sea-ice regions.

The anticipated phenological changes in GPP in a warmer Arctic indicate an overall early shift at all latitudes but large heterogeneity in the combination of bloom timing and intensity (Fig. 3). These findings stress the nonlinearity of the response of the biological community to quasi-linear latitudinal changes in physical drivers, highlighting the complexity of the sympagic system. Important trade-offs are at play below 74°N; the expected increases in production are limited by the diminished seasonal areas, while blooms in areas with expanding first-year ice are limited by narrow growth windows. These are also the latitudes currently characterized by most of the coastal Arctic sea ice and where the terrestrial inputs of nutrients are most relevant.

How can we incorporate the observational and climate model uncertainties in a more comprehensive overview of the potential impacts on the Arctic marine food web? The hybrid numerical model can also provide insights on the spread of algal response to a range of uncertain factors. We focused on the role of the unconstrained changes in future nutrient distributions (6, 7) and the wide variability of climate model projections of sea-ice season duration and snow depth (Fig. 2, A, B, and G) by considering more extreme values.

It is known that the final algal biomass is also determined by winter nutrient availability in the ocean (10, 14). The available measurements of winter dissolved silicate in the Arctic Ocean point to great heterogeneity among latitudinal bands (44). The effects of this latitudinal variation of nutrients on sea-ice production have been tested with a sensitivity analysis to the 10th and 90th percentiles and comparison with the median value (Fig. 5 and table S3), which was used in the presented standard simulations of Fig. 3. Overall, the different ranges of nutrients give values of sea-ice GPP and Chlorophyll-a (Chl-a) falling within the reported ranges in the Arctic Ocean (tables S2 and S3) (45), although within their lower limits. This is due to the fact that our results are representative of wide latitudinal bands, thus flattening the influence of possible hot spots of biological production that can be found in the literature. Little nutrient variability is currently observed below 66°N, and this has a negligible effect on current and projected sea-ice GPP (Fig. 5, 61°N; nutrient values are given in the legend). The largest observed nutrient heterogeneity is found between 66°N and 74°N (Fig. 5, 71°N), which is an indication of the diversity of oceanic and coastal conditions at this latitudinal band. The simulated algal biomass is directly proportional to nutrient availability (14, 24), which is an expected response given that the vernal development of sea-ice algae occurs when all nutrient replenishing mechanisms are at their optimal conditions. Above 74°N, the observed spread in nutrient distribution is between the one found in two other cases and has an expected intermediate effect on the resulted GPP (Fig. 5, 81°N).

Fig. 5 Changes in GPP associated with nutrient heterogeneity along latitudinal bands.

Simulation results at three representative latitudes (61°N, 71°N, and 81°N) for the historical period (1961–2005) and RCP8.5 scenario (2061–2100), showing sea-ice GPP under different levels of availability of dissolved silicate according to the calculated 10th, 50th, and 90th percentiles of measured winter silicates at a depth of 10 m between 1948 and 2000, retrieved from the Hydrochemical Atlas of the Arctic Ocean (44).

The presence of local upwelling and nutrient pulses generated by coastal features are likely to further enhance the productivity of the sea ice. South of the Arctic Circle, these events will be nevertheless confined to the strongly reduced areas of seasonal sea ice (Fig. 2D), while to the north, the compensation between areal expansion and early breakup would eventually maintain the current conditions. Given the wide spread of the tested nutrient concentrations, opposite speculations can be made, considering the projected enhancement in stratification that would result in lower or suppressed levels of nutrients (7) and thus reduced sea-ice primary production.

This uncertainty needs to be combined with the spread of the climate models and their limitations in deriving robust estimates of the projected trends in sea-ice seasonality (4648). In particular, CMIP5 models fail to capture the recent steepening trend in sea-ice extent (46), and the large spread has been linked to the model’s mean sea-ice thickness of multiyear ice, which is difficult to constrain given the limited observations (35, 48). Although this issue may not be as relevant for first-year ice projections—which are the focus of our study—we need to consider the possibility that, even within the chosen most pessimistic scenario (i.e., the RCP8.5), the ensemble medians would be too conservative.

More extreme scenarios can be assessed by considering the ensemble means of the single model 10th and 90th percentiles, which are shown in Fig. 6, A to D for a selected set of sea-ice drivers (see table S3 for a complete overview). Two extreme cases have been selected: a model realization characterized by the shortest sea-ice seasons along all latitudes (warm extreme; see Materials and Methods for details) and a drastic reduction to bare ice conditions (Fig. 6C). The latter case aims to analyze the impact that uncertainties associated with the modeled snow depth data (49, 50) may have on algal phenology and production. The results give a quantitative overview of the anticipated shifts in the algal bloom features as the projected sea-ice season shortens with the scenarios. While the bloom peak under current sea-ice conditions is evenly distributed across latitudes (Fig. 6E, circles), in the warm extreme case, it shows a rather linear response to shorter ice seasons (Fig. 6E, diamonds), as envisaged in conceptual models (17, 18). The higher latitudes will have more pronounced, earlier blooms as the ice season becomes dramatically shorter and algal production in an extremely warm climate is ultimately constrained by the duration of the ice season at all latitudes. In the scenario driven by ensemble medians (Fig. 6E, stars), ice algal phenology shows nonlinear interactions between sea-ice conditions and algal growth (Figs. 3 and 4).

This response is sensitive to the snow depth simulated by climate models, as indicated in the scenario experiment where the snow-covered surface is forced with prescribed bare ice; the bloom would advance by at least 50 days at all latitudes (Fig. 6E, squares and fig. S3). We notice that the GPP of the high latitude clusters identified in Fig. 3 would be considerably larger (fig. S3). The smallest differences would be expected at the lowest latitudes, since there is a comparable decrease of maximum snow depth (<0.1 m) in the RCP8.5 snow-covered scenario (Fig. 4, 61°N, right panel). At the mid-latitudes, the little increase seen for the snow-covered RCP8.5 scenario (Fig. 4, 71°N, right panel) would be replaced by a more intense increase because of the snow removal.

Fig. 6 Uncertainties in the future projections of sea-ice drivers and simulated responses of algal phenology.

Spread of the sea-ice driver projections and extreme cases for (A) ice season duration, (B) maximum sea-ice thickness, (C) maximum snow depth, (D) minimum near-surface air temperature. The lines are the ensemble means of empirical density distribution medians (50th percentile as in Fig. 2) for the RCP8.5 scenario (2061–2100) simulations. The shaded areas indicate the interval between the ensemble means of the 10th and 90th percentiles from the same distributions. The purple lines show the sea-ice drivers derived from the extreme warm scenario. The bare ice scenario is identical to the ensemble mean except for the snow depth set to zero [gray line in (C)]. (E) The resulting ice algae phenology, expressed as the day of bloom peak (y axis) under different scenarios (historical, ensemble RCP8.5, ensemble RCP8.5 for bare ice experiment, and RCP8.5 of the warmest extreme among models) and as a function of sea-ice season duration (x axis) for all the simulations and latitudinal bands.

These results, albeit constrained by the uncertainties discussed above, offer quantitative support to other paradigms linked to the loss of less-productive multiyear ice and the mismatch between the timing of primary and secondary production (18, 19). Our analysis indicates a doubling of first-year ice at the expense of multiyear ice at higher latitudes (Fig. 2D). The notable increase of annual production in first-year ice (Fig. 3) would make the Arctic sea ice at these latitudes a more productive system (Fig. 4, 81°N, right panel). At the same time, observations from the central Arctic (43, 51) suggest that thick ice may in fact be already more productive than initially thought; therefore, the projected thinner and more permeable multiyear ice might enhance sea-ice production even further.

While the seasonal timing of the sea-ice algal bloom is critical to the successful reproduction of copepod grazers, the timing of the phytoplankton bloom is critical to the survival of copepod offspring (52). Disruption of the seasonality of the algal blooms by the sea-ice changes analyzed here (Figs. 2 and 6E) has already created mismatches for the timing of zooplankton production, which has consequences for higher trophic levels (17, 22, 51). These copepods are highly adapted to feed on ice algae, and their overwintering survival is phenologically synchronized with the availability of this food source. We have shown here that earlier breakups, associated with altered photoperiods and thinner sea ice and snow, can shorten the length of the biological time window at mid-latitudes, thus setting an upper limit for ice algal production, which is not the case at low and high latitudes (Figs. 3 and 4). This might affect timing and amplitude of zooplankton production, which in turn might affect their main predator, Arctic cod, up to the top predators at the highest trophic levels. Because of the diverse latitudinal response of sea-ice algae (Fig. 3) to changes in climatic drivers (Fig. 2), we confirm the possibility of phenological uncoupling with secondary and tertiary consumers as described by Post (22).

The recent shifting paradigm toward a high level of wintering activity of plankton populations (53) also suggests that we should revisit our way of considering the winter Arctic marine ecosystem. The projected increase and advanced sea-ice GPP found in our simulations (Fig. 3) may also be an earlier potential pulse of energy for the active winter plankton communities in the low latitudes of the Arctic, although the large reduction in first-year ice extent would likely limit this effect.

Since the Arctic marine food web is short, poorly diverse, and seasonally driven by limited pulses of energy, the projected latitudinal changes in ice algal phenology have the potential of trophic cascades and phenological shifts at higher trophic levels. The phenological synchrony across trophic levels is expected to be driven by the pace of advancing sea-ice breakups (22) and consequent shift in algal bloom timing (this study). Thus, as supported by our results in Fig. 3, we might expect the strongest phenological cascades (i.e., shifts in timing) at the lowest latitudes where algal blooms are projected to occur significantly earlier, while the largest trophic cascades (i.e., changes in abundance at adjacent trophic levels) at the highest latitudes where the largest relative increase in biomass is projected.

Sea ice is inhabited also by other algal groups besides the dominant pennate diatoms (40). We thus may expect the potential development of new acclimation and/or adaptation patterns or the shift in dominance of other competitive algal groups in a warmer Arctic. The modeling of these ecological dynamics would require an increase in model complexity beyond that of the present study, which was constrained by the current limited knowledge on the ecological and physiological traits of other algal groups (39). To resolve the local features driving algal changes, we urgently need organized multidisciplinary studies of sea-ice physics, biogeochemistry, and associated biota on multiple temporal and spatial scales to achieve the proper level of details that models would need to make reliable future projections of algal phenological and trophic changes at those scales.


Study design

Climate model data from the CMIP5 (26) were obtained from the Earth System Grid Federation (ESGF) data distribution center. The years 1961–2005 were used to represent historical climate conditions, and the years 2061–2100 from the RCP8.5 (28) simulation were used as an example of a warming business-as-usual scenario. The models (listed in table S1) were selected according to the availability of daily data of sea-ice thickness and sea-ice concentration at the time of the analysis. We also did not exclude models based on performances because our method works on empirical probability density functions and not on zonal or monthly averages; therefore, we sought to preserve any possible interannual or regional variability that each model can elicit. In assessing the projections for an ice-free Arctic, Liu et al. (54) retained the models with the simulated September sea-ice extent falling within 20% of the observations. Twenty of the 30 models they examined satisfied that requirement, and 16 of the 18 models we considered are among these 20 [the additional two models we used were not available at the time of writing of Liu et al. (54)]. All original observations and model data grids were preserved, and stored daily data were used when available without any further spatial or temporal averaging. The observations of sea-ice concentration for current climate conditions (1979–2005) were obtained from Nimbus-7 SMMR and DMSP SSM/I-SSMIS passive microwave remote sensing data (32).

The indicators for the freeze-up and breakup dates as a function of latitude were obtained by computing the empirical probability density function of the days at which a chosen threshold of 15% sea-ice cover was passed, counting all model grid points in the interval of 2° latitude. A positive transition indicated a freeze-up event, while a negative transition indicated a breakup event. The method is robust to the choice of the sea-ice variable because it gives the same results if sea-ice thickness with a threshold of 0.1 m is used. It is also not sensitive to the subjective choice of the thresholds as long as they are considered representative of low sea-ice conditions; a variation of ±50% around the chosen values did not lead to substantial changes in the overall statistics. To avoid counting successive events of melting and refreezing, it was decided to only consider the last melting event before the first freezing of the next sea-ice season and the last freezing event before the seasonal melting. During data collection, the algorithm also collected the area of the ocean (within the latitudinal band) that was characterized by first-year ice. The same algorithm was applied to the historical climate and future scenario simulations, as well as to satellite data for validation. Not all variables were available for all the models in the ESGF servers because of the recent overhauling (table S1). Given the use of empirical probability density functions, the medians and other statistical descriptors were nevertheless computed over a large number of data points. The CMIP5 daily variables considered were near-surface atmospheric temperature and sea-ice thickness, while snow depth was available from CMIP5 models on a monthly basis. All values not corresponding to a sea-ice–covered ocean were masked. The analysis was performed for first-year ice only, with 2° latitudinal bands from 58°N to 84°N. The latitudinal empirical probability density functions and the relative statistical descriptors (10th, 50th, and 90th percentiles) were computed for freeze-up day, breakup day, maximum sea-ice thickness, day of maximum sea-ice thickness, maximum snow depth, month of maximum snow depth, minimum near-surface temperature, and day of minimum near-surface temperature (see fig. S4 for an example). If the corresponding first day of sea ice fell within the same (or previous) month of the snow season starting, then we used that day as the starting day of the snow season and the sea-ice season. Otherwise, we used the first day of the following month for the snow season. For ascribing the last day of snow cover, if the snow season showed to end within the same month when sea ice vanished, then we chose the first day of the corresponding month, otherwise, we used the 15th of the previous month. The day of maximum snow depth was instead assumed to be the midmonth day of the corresponding month.

The ensemble mean of the median values of the empirical probability density functions was used to construct an analytically idealized ice season (i.e., snow and sea-ice dynamics) for the historical and future climate at every band of latitudes as proposed by Tedesco and Vichi (24). Sea-ice thickness was prescribed to grow with a cubic function from the first freezing event (i.e., below freezing temperature and the sea-ice thickness >0.1 m) until the day when sea-ice thickness reached its seasonal maximum. This method takes into account the reduced growth due to sea-ice thickening and aging. After the seasonal maximum, sea ice was modeled to start melting with a cubic function. Snow was prescribed to accumulate linearly on sea ice from the first snow event (freezing temperature and snow depth of >0 m) until its maximum and then to melt linearly until the last day of snow cover (snow depth of >0 m). The surface temperature was modeled to decrease linearly from the freezing temperature to its minimum (dependent on latitude) and then to increase linearly until it reached the freezing point, which was then kept until all sea ice and snow were melted completely.

The ensemble means of the 10th and 90th percentiles of the single model distributions were used to analyze the impact of extreme case scenarios (Fig. 6). It is physically inconsistent to combine the extreme drivers to create a scenario case where, for instance, the freeze-up occurs very early in the year, the snow depth is minimal, and the minimum air temperature is the lowest. The adopted solution consisted of selecting one model realization that represented the combination of the driver values in the tails of the distributions in a physically consistent manner. We selected the model with the shortest duration of the sea-ice season because all the other drivers were then generally located around or beyond the ensemble mean of the percentiles (Fig. 6A to D).

At all latitudes, surface seawater salinity was fixed at 31 (practical salinity unit) and wind speed was fixed at 5 m/s, while sea-ice bulk salinity (Sbu,fy) was computed according to the analytical function provided by (55) for Arctic first-year ice, which is dependent only on the normalized ice depth (0 ≤ z ≤ 1) and is based on a fitted salinity profile derived from a complex parameterization of sea-ice bulk salinitySbu,fy(z)=za+bz+c(1)where a = 1.0964, b = −1.0552, and c = 4.41272.

Incoming solar radiation was retrieved for both scenarios and at each considered latitude from the ERA-interim European Centre for Medium-Range Weather Forecasts climatological (1979–2010) downward surface solar radiation (56). The incoming radiation was converted to visible wavelengths (400 to 700 nm) and then to photosynthetically active radiation (PAR) and, lastly, was corrected for surface albedo and light extinction through snow and sea ice with a Semtner-like sea-ice thermodynamic model (57). Because of the lack of information from CMIP5 models on sea-ice halodynamics, we limited our investigation to the biological production of the bottom sea-ice habitable space, operationally defined as the sea-ice biologically active layer (BAL) (57). This layer comprises the permeable part of sea ice that extends from the bottom-most sea ice upward until the brines volume reaches the 5% limit. The BAL fraction can be calculated given its boundaries. At the lower sea-ice interface with the ocean, seawater is assumed to be at freezing point, and brine volume is considered to be at 100%, i.e., all medium is liquid. At the upper interface, the BAL is delimited by the brines volume being 5% at the maximum. Given the computation of sea-ice bulk salinity after (56), brine salinity and brine/sea-ice temperature were calculated from thermal equilibrium. Given the BAL sea-ice temperature, the dynamics of the BAL thickness were reconstructed with a cubic function, as done for sea-ice thickness. The computed physical variables needed by the ocean component of the biogeochemical model were seawater temperature, seawater salinity, wind speed, and PAR. The computed physical variables needed by the sea-ice component of the biogeochemical model were BAL thickness, temperature, bulk salinity, brine salinity, brine volume, and PAR.

The biogeochemical model set-up

The sea-ice biogeochemical model (24) allowed flexible stoichiometry in the algal composition and physiological responses (58); in this specific case, it also featured a nonlimiting amount of dissolved nitrogen and phosphorous and realistic historical concentrations of dissolved silicates (SiO4), assuming no intracellular storage and a Monod-like uptake rate (59) as modeled in phytoplankton. This choice was made (i) to limit the uncertainties due to the little available knowledge on sea-ice nutrient co-limitation; (ii) because of previous sensitivity studies on the role of nutrients on algal bloom timing (24); (iii) because silica is an essential macronutrient for the generally dominant functional group of algae found in the Arctic bottom sea-ice habitat, which is pennate diatoms (14, 39). Pennate diatoms were characterized by three basic constituents: carbon, Chl-a, and silica. Note that the internal stoichiometric ratios were allowed to vary, particularly the Chl:C ratio, which is an emulator of light acclimation capability. Other functional groups included detritus and gases, for nine state variables in total. The simulated sea-ice biological processes were photosynthesis, respiration, mortality/excretion, and nutrient uptake. The sea-ice biogeochemical model was coupled to pelagic biogeochemistry in a simple slab ocean, representing an averaged mixed layer depth under sea ice (15 m for all latitudes and both scenarios). The ocean model was initialized in both scenarios with the zonal winter median (November to May) concentrations of dissolved silicates in the first 10 m of the Arctic Ocean between 1948 and 2000 (44) for every latitudinal band. 1 mg C/m3 of diatoms in seawater for all latitudes and in both scenarios were used as the typical background winter value. The sea-ice model did not need initialization since the fluxes at the interface, which are a function of sea-ice growth/melt velocities and thickness of the BAL, controlled the exchanges of dissolved and particulate matter between sea ice and seawater (58). With this setup, we investigated the relative changes in future GPP, which, assuming unchanged nutrient dynamics in the future, are independent on the background nutrient conditions used, which are highly uncertain (6, 7). This operative strategy thus removes the uncertainties associated with current and future nutrient dynamics, while it increases our confidence toward reliable estimates of future sea-ice GPP.

The biogeochemical prognostic variables presented are sea-ice algal Chl-a and GPP in the brine pockets and their vertical integral over the BAL, both temporally integrated over the duration of the sea-ice season. An indicator of phenology, the anomaly index of chlorophyll concentration in the brine volume has been developed (24), where positive values represent the period of bloom activity in sea ice. This index was computed as the difference between the amplitude of the concentration and its SD over the production period.


Supplementary material for this article is available at

Table S1. List of the 18 CMIP5 models and of the related sea-ice variables used in this study.

Table S2. Temporally and spatially integrated GPP of Arctic first-year ice along latitudes and scenarios.

Table S3. Maximum sea-ice algal Chl-a of bottom Arctic first-year ice along latitudes and scenarios.

Fig. S1. Sea-ice algae Chl-a anomalies.

Fig. S2. Sea-ice GPP.

Fig. S3. Latitudinal changes in the bloom phenology and sea-ice annual GPP in the case of RCP8.5 bare ice scenario.

Fig. S4. Empirical probability density functions of sea-ice drivers for the recent historical (1961–2005) and scenario simulations (2061–2100) for an example CMIP5 model (GFDL-CM3).

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (table S1) for producing and making the model output available. We also acknowledge the use of the BFM system ( We are thankful to J. Field for proofreading an earlier version of this work and to two anonymous reviewers for valuable comments. L.T. thanks colleagues from the BEPSII (Biogeochemical Exchange Processes at the Sea-Ice Interfaces) WG for valuable discussion on sea-ice matters during the past years and BEPSII funding agencies (CliC/IASC/SOLAS/SCAR) for making this network possible. M.V. acknowledges the availability of the graphical package. Funding: M.V. acknowledges funding from DST/NRF through the South African National Antarctic Programme. Author contributions: L.T. developed the original idea that led to this paper. M.V. and E.S. collected and analyzed the climate model simulations. L.T. prepared the combined physical forcings and performed and analyzed the biogeochemical model simulations. L.T. and M.V. wrote the main paper, with input from E.S. All authors discussed the results and implications and participated in the manuscript improvement at all stages. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All climate model output data from the CMIP5 are publicly available at and at other designated ESGF data centers. Given that some data used in this study are currently not found on ESGF servers, they are made available from the corresponding author upon request, as well as the scripts to process them. The National Snow and Ice Data Center sea-ice concentration data and the ERA-Interim reanalysis data are publicly available through the respective web portals. The sea-ice biogeochemistry model is available on the BFM website The empirical probability density functions for all the models are available as further additional supplementary material at All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article