Heterogeneous future Arctic Ocean primary productivity changes projected in CMIP6

The Arctic Ocean is experiencing profound environmental changes due to climate change, with Net Primary Production (NPP) broadly projected to increase this century. This study analyzes NPP trends and their drivers across pan-Arctic and sub-regional scales throughout the 21st century, comparing Coupled Model Intercomparison Project Phase 6 (CMIP6) and Phase 5 (CMIP5) projections to assess how model generations differ. Using a multi-model approach, we assess projections for different Phytoplankton Functional Types (PFTs), diatoms and nanophytoplankton, and examine the role of physical and biogeochemical constraints including light, nutrient, and temperature limitations. Our results reveal that Arctic Ocean NPP increases are primarily driven by reduced sea ice cover, leading to longer ice-free seasons in the expanding seasonal ice zone. However, NPP changes exhibit pronounced spatial heterogeneity, with strong increases in Arctic inflow shelf regions, tempered by decreases in Baffin Bay and Nordic Seas. These differences are due to the varying balance between physical and biogeochemical NPP constraints across the Arctic Ocean. The multi-model mean Arctic Ocean NPP increase is four times larger in CMIP6 than in CMIP5, under comparable radiative forcing, with a three times higher uncertainty at the end of the century. This difference is attributed to higher baseline nutrient levels in CMIP6, combined with more pronounced sea ice loss and greater warming than in CMIP5. Key aspects to better simulate future Arctic Ocean NPP remain the representation of present-day nutrient levels, light transmission through sea ice and reduced model uncertainty in climate sensitivity.

1 Introduction

Phytoplankton are marine microorganisms that serve as primary producers in ocean ecosystems, forming the foundation of marine food webs (Vincent and Laybourn-Parry, 2008). Their development depends on biomass accumulation and growth rates, which are controlled by environmental factors such as light availability, nutrient concentrations and temperature. Marine Net Primary Production (NPP) represents a fundamental indicator of phytoplankton activity, quantifying net carbon fixation after losses due to cellular respiration and maintenance processes are accounted for. As an indicator of primary producer activity, NPP drives the biological ocean carbon pump and therefore plays a key role in long-term ocean carbon sequestration (Sarmiento, 2013).

The dependence of NPP on environmental conditions makes it sensitive to climate change. Global NPP has declined by −2.1 % per decade from 1998 to 2015 according to satellite observations (Gregg and Rousseaux, 2019). Model hindcast simulations indicate NPP has declined by −6.5 % since 196) and is projected to continue to decline throughout the 21st century (Bopp et al., 2013). This decrease is particularly apparent in tropical and mid-latitude regions where vertical mixing and nutrient supply decline due to enhanced upper-ocean stratification (Doney, 2006). In contrast, at high latitudes, and in particular in the Arctic Ocean where NPP is primarily light limited, observations suggest a recent increase in NPP, rising by 57 % between 1998 and 2018 (Lewis et al., 2020). However, there is disagreement on the sign of the projected changes in these regions (Laufkötter et al., 2015).

The Arctic Ocean is particularly sensitive to the effects of climate change (Kwiatkowski et al., 2020; Laufkötter et al., 2015; Vancoppenolle et al., 2013) and is projected to experience warming nearly four times the global average (Rantanen et al., 2022). This warming results in ocean freshening through increased freshwater input from multiple sources, including melting glaciers and sea ice, enhanced precipitation, and greater river inflow (McCrystall et al., 2021; Shu et al., 2018). Moreover, climate change also impacts global ocean circulation patterns such as the Atlantic Meridional Overturning Circulation (AMOC), which is an important regulator of the climate system, and which could in turn have an impact on the Arctic Ocean and thus on NPP (Weijer et al., 2020). In polar regions, stratification is primarily driven by salinity gradients rather than temperature differences, making this freshening particularly impactful for ocean mixing and circulation patterns (Timmermans and Marshall, 2020).

The Arctic Ocean phytoplankton community consists of diatoms and smaller taxa including coccolithophores, prymnesiophytes, flagellates, and picoeukaryotes (Ardyna and Arrigo, 2020). Climate change is projected to drive a compositional shift toward smaller phytoplankton (Bopp et al., 2005), fundamentally altering marine food web structure and efficiency. This size-structure transition extends food chain length and reduces trophic transfer efficiency, consequently reducing carbon export. The shift particularly weakens export through two mechanisms: smaller phytoplankton are more efficiently recycled within the microbial loop, while larger phytoplankton like diatoms contribute disproportionately to carbon export via gravitational sinking. These structural modifications have cascading effects on both ecosystem functioning and carbon sequestration capacity (Grebmeier et al., 2010; Ward et al., 2012). Furthermore, phytoplankton community restructuring propagates through upper trophic levels, ultimately impacting commercial fisheries and marine food security (Ardyna and Arrigo, 2020; Hegseth and Sundfjord, 2008; Neukermans et al., 2018). The implications of global warming on the Arctic Ocean extend beyond primary productivity and encompass Arctic Ocean biogeochemical cycles, especially the carbon pump and the regional contribution to climate regulation (Sarmiento, 2013).

Arctic Ocean NPP has increased by > 50 % over recent decades, according to satellite observations (Arrigo and van Dijken, 2015; Lewis et al., 2020). This increase is largely attributed to the decrease in sea ice coverage and thickness, as well as associated changes in sea ice scape, which enhances light penetration into the water column and fosters higher phytoplankton productivity (Lannuzel et al., 2020). Additionally, riverine fluxes and coastal erosion supply a substantial amount of nutrients, supporting approximately one-third of Arctic Ocean productivity (Terhaar et al., 2021).

The Arctic Ocean is projected to experience a continued increase in NPP under climate change (Tagliabue et al., 2021). Future sea ice retreat will further extend the open-water season duration, enhance light availability and stimulate phytoplankton growth. However, the melting of sea ice will also introduce large volumes of freshwater, which may further stratify the water column, reducing nutrient transport from deeper waters and potentially modifying NPP dynamics (Pabi et al., 2008; Popova et al., 2010) and the distribution of the phytoplankton in the water column (Steiner et al., 2016).

CMIP5 Arctic Ocean NPP projections exhibited considerable model divergence with even the sign of future Arctic Ocean NPP anomalies uncertain (Vancoppenolle et al., 2013). However, there is model agreement on the driving mechanisms: decreasing sea ice extent (and thus increasing light availability) and declining nitrate concentrations, almost reaching an oligotrophy onset in 2100. Increased temperatures could also contribute to enhanced NPP, however the analysis of CMIP5 projections did not consider thermal sensitivity as a mechanism influencing NPP (Vancoppenolle et al., 2013). There is also a high uncertainty amongst CMIP5 models concerning the deepening of the future Subsurface Chlorophyll Maximum (SCM) (Steiner et al., 2016).

Despite the recognized importance of environmental drivers in controlling Arctic phytoplankton productivity, several critical knowledge gaps remain. First, it is unclear how projections of key environmental factors (light availability, nutrient concentrations, and temperature) have evolved between CMIP5 and CMIP6 model generations, and whether improved model physics in CMIP6 has led to more consistent or divergent projections. Second, the relative importance of these drivers in CMIP6 varies spatially across the Arctic Ocean, but a systematic analysis of these spatial patterns is lacking. Third, different PFTs respond differently to environmental changes, yet how these differential responses are captured across regions remains poorly understood. Finally, while previous studies have identified sign inconsistencies in CMIP5 NPP projections (Vancoppenolle et al., 2013), the mechanisms underlying the substantially larger NPP increases projected by CMIP6 models have not been systematically analyzed in the Arctic Ocean (Tagliabue et al., 2021). These gaps limit our ability to assess the reliability of future Arctic Ocean productivity projections and their implications for marine ecosystems and biogeochemical cycles.

To address this knowledge gap, we investigate the environmental drivers of phytoplankton growth using both CMIP6 and CMIP5 models. Our analysis aims to understand how projections and associated uncertainties have changed between the two model generations over the 21st century, examine spatial variations in these projections across the Arctic Ocean, and assess how different phytoplankton functional types respond to these environmental changes.

2 Material and Methods

2.1 Model selection

We used output from the Coupled Model Intercomparison Project Phase 6 (Eyring et al., 2016) downloaded from the Earth System Grid Federation (ESGF) servers (Table 1). Each ESM member is identified by a nomenclature composed of letters and various numbers, such as rWiXpYfZ, where W indicates the realization number (different for each run), X denotes the initialization method, Y refers to the set of physical parameterizations used, and Z represents the forcing configuration number. For each model, we prioritized selecting the ensemble member r1i1p1f1 when available. If r1i1p1f1 was not available, the ensemble member most similar to it, such as r2i1p1f1, r1i1p1f2, or the closest alternative was chosen. The historical period spans 1850–2014, while the SSP scenarios cover 2015–2100. Five SSP scenarios were selected: SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP5-3.4 and SSP5-8.5 (Meinshausen et al., 2020).

Table 1 Description of CMIP6 models used in this study, along with the associated available variables and simulations. The first group includes models used to project environmental conditions (NPP, NO₃, SIC, and SST). The second group includes models used for both environmental projections and analysis of phytoplankton growth limitation terms (N_lim, L_lim and T_f).

To compare with the previous generation of models, from the Coupled Model Intercomparison Project Phase 5 (Taylor et al., 2012), we downloaded output from ESGF servers (Table 2) for the scenario RCP8.5 (van Vuuren et al., 2013). As with CMIP6, we prioritized r1i1f1 when possible. The CMIP5 historical period spans from 1850 to 2005 while the RCP8.5 scenario covers the period 2006–2100.

Table 2 CMIP5 models used in this study, along with the available variables and simulations.

Although additional scenarios were briefly examined, the detailed analysis is restricted to SSP5-8.5 and RCP8.5 (as in Vancoppenolle et al., 2013), since the primary objective of this study is to investigate the dominant mechanisms and processes underlying future changes rather than to provide precise scenario-dependent projections.

2.2 Models diagnostic

2.2.1 Study region and interpolation methods

The Arctic Ocean was defined as the waters above 66.5° N latitude and divided into 10 different regions, following (Arrigo and van Dijken, 2011). Both model output and observations were interpolated onto a regular latitude-longitude grid with a resolution of 360×180 using the CDO (Climate Data Operators) software with the remapdis interpolation method (Schulzweida, 2023).

2.2.2 NPP diagnostics

Model outputs are not always available for all variables of interest, so we categorized the models into two groups based on data availability (Table 1), corresponding to the two sections of our study: In the first section, we examined the prognostic environmental variables that may influence NPP, namely NO₃ concentrations, sea ice concentrations and sea surface temperatures (respectively named in CMIP6: NO₃, sic and tos) and in the second section we investigated the light, nutrient and thermal limitation terms that directly impact phytoplankton growth rates. The first section uses outputs from Group 1 and Group 2 models, while the second section focuses on models from Group 2, which provide direct limitation term diagnostics. These limitation terms are derived from the variables mentioned above and represent key factors influencing phytoplankton growth rate (μ) for a given phytoplankton type:

$$\begin{array}{}\text{(1)}& \mathit{\mu }={\mathit{\mu }}_{max}\times T_{\mathrm{f}}\times L_{lim}\times N_{lim}\end{array}$$

where μ_max is the maximum growth rate, L_lim and N_lim are the limitation terms for light and nutrient respectively, and T_f represents the temperature function.

2.2.3 Nutrients, SST, Sea ice and limitation factors

Nitrate (NO₃) was selected as the limiting nutrient, as the Arctic Ocean is primarily nitrogen-limited, and ammonium (NH₄) is not available for all models (Browning and Moore, 2023; Tremblay and Gagnon, 2009). Two different approaches were used to calculate nitrate concentrations. When data were available, nitrate concentrations were integrated over the mixed layer (MLD). However, as mixed-layer depth and vertically-resolved NO₃ data were not available for all scenarios (especially for CMIP5), some analyses were conducted using surface nitrate concentrations, i.e., from the upperocean model layer. An oligotrophy threshold was set at 1.6 mmol m⁻³, i.e., at the half-saturation constant of diatoms for NO₃ uptake (Sarthou et al., 2005). Additional analysis is provided in the appendix, where NO₃ is integrated over the upper 100 m of the water column (see Fig. A1).

As a diagnostic of ice coverage, we used monthly mean sea ice concentration (SIC) outputs. From this, we derived the September sea ice area (SSIA), the mean September sea ice concentration (SSIC), and the mean ice-free season duration (IFSD). The IFSD was defined as the number of months with sea ice concentration lower than 15 % (Lebrun et al., 2019).

We analyzed growth rates for the two PFTs available in model output: diatoms and miscellaneous phytoplankton. The direct limitation terms provided by the models include: the light limitation term for miscellaneous phytoplankton ($L_{lim}^{\mathrm{misc}}$, named limirrmisc in CMIP6 terminology), the light limitation term for diatoms ($L_{lim}^{\mathrm{diat}}$, limirrdiat), the nutrient limitation term for miscellaneous phytoplankton ($N_{lim}^{\mathrm{misc}}$, limnmisc) and the nutrient limitation term for diatoms ($N_{lim}^{\mathrm{diat}}$, limndiat). These limitation terms can vary between 0 and 1, with lower values indicative of greater constraint on phytoplankton growth rates.

The temperature function (T_f) is not directly provided by the models, but it is calculated using the Eppley function ($T_{\mathrm{f}}=e^{K_{\mathrm{Eppley}}^T}$, where K_Eppley=0.063 °C⁻¹) except for UKESM1-0-LL which calculates the temperature function as T_f=1.006^T, based on vertically gridded potential temperature (3Dtemp, thetao). Monthly T_f was computed for each depth level and weighted according to the biomass of each PFT (phymisc or phydiat) to obtain PFT-specific depth integrated temperature functions for each model ($T_{\mathrm{f}}^{\mathrm{misc}}$ and $T_{\mathrm{f}}^{\mathrm{diat}}$).

2.3 Data used for model evaluation

2.3.1 Primary Production

The primary production observational data product used in this study provides recalculated depth-integrated NPP values based on satellite-derived chlorophyll concentrations, sea surface temperature, and sea ice cover (Lewis et al., 2020). This dataset spans 1998–2018 and offers satellite-based estimates in the Arctic Ocean. It was regridded onto a 1°×1° spatial resolution. Only NPP data from 1 March to 30 September were considered for model evaluation, due to the lack of ocean color observations during the dark season.

2.3.2 Nitrate (NO₃)

To evaluate nitrate concentrations, we used a climatology derived from in situ hydrographic observations from the World Ocean Atlas 2018 (Garcia et al., 2019; Randelhoff et al., 2020). This dataset provides an average nitrate concentration from 1978 to 2017 on a 1°×1° grid. Observed nitrate concentrations were extracted over the upper 100 m of the ocean. Depending on analysis, depth-averaging was performed consistently with model output, as detailed in Sect. 2.2.3. However, nitrate data in the Arctic Ocean remain sparse, particularly in summer, and may be subject to significant biases (Gibson et al., 2020).

2.3.3 Sea ice concentration

Passive microwave satellite sea ice concentration fields were taken from the OSI-SAF Global sea ice concentration data record v3 (EUMETSAT OSI SAF, 2022). This dataset provides daily averaged sea ice cover (in percentage) at 25 km spatial resolution, covering the period 1979–2022.

2.3.4 Sea Surface Temperature (SST)

The sea surface temperature observational dataset is derived from the Climate Change Initiative (CCI) satellite product, covering the period 1982–2019, providing daily analyses at a high spatial resolution of 0.05° (Embury and Good, 2024). This dataset provides a long-term, high-quality observational record of SST, allowing for comparison with model outputs.

3 Results

3.1 Model evaluation

Simulated ranges of NPP, MLD, NO₃ concentration, SSIA and SST encompass data-based estimates over the observational period (Fig. 1). While the simulated trends in Arctic Ocean warming and declining SSIA are broadly consistent with observations, simulated increases in NPP are low biased. The multi-model mean NPP is slightly higher than observed, over the 1998–2018 period over which all models consistently simulate an NPP increase. Observed NPP increases from 256 in 1998 to 391 Tg C yr⁻¹ in 2018, increasing by 135 Tg C yr⁻¹ in 21 years. The increase in NPP in models is lower, increasing from 354 to 394 Tg C yr⁻¹, increasing by 40 Tg C yr⁻¹ in 21 years (Fig. 1a, Table 3). The multi-model annual mean NO₃ concentration for the period 1979–2005 is 6.90 mmol m⁻³, which is higher than the observations of 4.67 mmol m⁻³ for the same period. However, this discrepancy is not unexpected, as Arctic Ocean NO₃ observations contain significant uncertainties due to limited accessibility beneath sea ice cover (Fig. 1b). Models and observations agree on a loss of sea ice coverage, although the simulated rate of decrease is lower than observed (Fig. 1c, Table 3), confirming previous studies (Notz and Community, 2020). Similarly, observed Arctic Ocean warming is ∼ 50 % higher than that simulated over the same period (Fig. 1d, Table 3).

Figure 1 Projections of Arctic Ocean (a) annual mean depth-integrated NPP, (b) NO₃ concentration in the mixed layer depth (MLD), (c) September Sea Ice Area (SSIA), and (d) annual mean SST, derived from individual CMIP6 models (Group 1) using historical and SSP5-8.5 simulations over the period 1900–2100. In panel (b), three metrics are displayed as multi-model mean (black lines): the annual maximum (February, upper line), annual mean (middle line), and annual minimum (July, lower line); individual model trajectories are shown only for the annual mean. The horizontal grey dashed line indicates the NO₃ concentration threshold associated with oligotrophic conditions. In all panels, colored lines represent individual model simulations, the thick red line represents the observations, and the thick black line denotes the multi-model mean.

Table 3 Observed and simulated NPP, SSIA and SST values throughout the observational period. Nitrate concentrations are not given as an Arctic Ocean time series is not available.

The drivers of these projected anomalies are changing in the same direction across all Arctic subregions: NO₃ and September sea ice decrease and SST increases for all models and regions, consistent with earlier CMIP-based assessments (Vancoppenolle et al., 2013) (see Figs. A1, A2, A3, A4).

Several additional figures on individual model and regional analysis are provided in the appendix.

The spatial distribution of key variables across the Arctic Ocean is generally well represented by the multi-model mean and comparable to observations despite some regional discrepancies that reflect both model limitations and observational uncertainties (Fig. 2).

Figure 2 NPP, NO₃ (integrated over the upper 100 m of the water column), IFSD, and SST, over the recent past, as represented by the multi-model mean (left), observational data products (middle) and the multi-model standard deviation (right). The averaging period corresponds to the availability of observations prior to 2015. Model outputs are combined from historical simulations for years prior to 2015 and SSP5-8.5 thereafter. For NPP, the March to September average is shown.

March–September NPP resembles satellite-derived estimates, with the highest values in the Nordic Seas, the Bering Strait, and along the Siberian shelf, while the Arctic Basin exhibits the lowest NPP values (Fig. 2a and b). The multi-model mean IFSD and SST spatial distributions also closely align with observations. The warmest regions (Nordic Seas, Barents Sea, and Chukchi Sea) correspond to areas with the least sea ice (Fig. 2g, h, j and k). While the multi-model mean represents low values for coastal NO₃ concentrations, in the Chukchi Sea and along the eastern Arctic Basin, it may overestimate nitrate concentrations in the Arctic Basin, though observations in this region remain limited. The coastal shelves of the Kara, Laptev and East Siberian Seas show low NO₃ values – of less than 1.5 mmol m⁻³ in both models and observations (Fig. 2d and e).

3.2 Projections of NPP and its environmental drivers

Across all SSPs, Arctic Ocean NPP and SST are projected to increase over this century, while nitrate concentration and IFSD decline (Fig. 3). However, the magnitude of change in each variable depends on the scenario, with generally larger anomalies and associated uncertainty under higher radiative forcing. An exception to this is SSP5-3.4, for which there is lower model uncertainty, possibly due to the smaller model ensemble (Fig. 6).

Figure 3 CMIP6 projected pan-Arctic Ocean anomalies of NPP (a), surface NO₃ (b), SSIA (c) and SST (d) over the period 1900–2100 relative to 1850–1899, for the different emission scenarios. Vertical bars represent multi-model means and standard deviation in 2100. The last vertical bar and triangle represent the range and mean value for the CMIP5 RCP8.5 scenario in 2100. The numbers within parentheses in the legend give the number of available models. Additional figures representing the timeseries for subregions are shown in the appendix.

Projected Arctic Ocean NPP and nitrate inter-model uncertainty is greater than inter-scenario uncertainty by the end of the century (Fig. 3a and b). Contrary to the other variables, IFSD exhibits reduced uncertainty under higher radiative forcing scenarios. This reflects the strong model agreement that the Arctic Ocean will be ice-free for most months by the end of the century (Fig. 3c). Surface ocean warming shows significant sensitivity across scenarios, with the greatest warming (multi-model mean +4.78 °C in 2081–2100), simulated in SSP5-8.5 (Fig. 6d).

Projected twenty-first century NPP increases are highly divergent between CMIP6 and CMIP5 scenarios of comparable high radiative forcing (Fig. 3a). For RCP8.5, the relative anomaly in NPP reaches 11.5±15.7 % at the end of the century. In contrast, SSP5-8.5 projects a much larger increase of 46.4±53.7 %. This highlights greater CMIP6 model agreement of increasing future NPP alongside a more than 3-fold increase in associated uncertainty (Fig. 3a). The projected twenty-first century Arctic Ocean NO₃ surface concentration decline is higher in SSP5-8.5 (4.38±1.98 mmol m⁻³) than in RCP8.5 (3.09±1.45 mmol m⁻³, Fig. 3b). The increase in IFSD is greater in CMIP6 than in CMIP5 (7.8 months vs. 4.2 months), indicating a stronger reduction in sea-ice duration in CMIP6 compared to CMIP5. However, when expressed as anomalies relative to the reference period, the magnitude of change is slightly lower in CMIP6 than in CMIP5 (Fig. 3c). The surface ocean warming is two times higher in SSP5-8.5 (4.78±2.09 °C) than in RCP8.5 (2.32±0.75 °C) with higher associated uncertainty (Fig. 3d).

Multi-model mean declines in NO₃ and increases in IFSD and SST occur across the Arctic Ocean domain, however NPP exhibits both regional-scale increases and decreases (Fig. 4). NPP increases are greatest in the northern Barents Sea, along the east coast of Greenland, and in the Chukchi Sea, reaching more than 20 g C m⁻² yr⁻¹. In contrast, there are NPP decreases along the coast of the Kara Sea, in the North Atlantic, and in Baffin Bay (Fig. 4a). NO₃ concentration decreases across the entire basin, with the most pronounced reductions on the Siberian shelf and off the East Siberian Shelf, where declines can exceed 4 mmol m⁻³ (Fig. 4b). The region experiencing the highest increase in the ice-free season duration is the northern Barents Sea, east of Svalbard Archipelago. In the central Arctic Basin, sea ice decline is also strong, with a loss of more than 6 months ice coverage (Fig. 4c). The Barents and Greenland Seas are the regions most affected by rising SST, with increases of up to 6 °C (Fig. 4d).

Figure 4 CMIP6 multi-model mean anomalies of NPP (a), NO₃ (b), IFSD (c) and SST (d). Anomalies are for 2081–2100 of SSP5-8.5 relative to 1995–2014 of the historical simulation. Hatched areas indicate where < 80 % of the models agree on the sign of the change.

Most models project consistent NPP increases under SSP5-8.5 except over the shelves (where anomalies are less pronounced) and in the central Arctic Ocean. For other variables, such as nitrate concentration, IFSD and SST, models agree on the sign of future change across the entire basin (Fig. 4).

By the end of the century, most CMIP6 models simulate a transition toward oligotrophic conditions, with an increasing fraction of ensemble members exhibiting oligotrophy and longer seasonal persistence of these conditions (Fig. 5). The seasonal cycle of nutrient concentration indicates that during the period 1995–2014, only seven CMIP6 models reach oligotrophy during the productive period, generally from mid-May to September (CanESM5, MPI-ESM1-2-LR, MPI-ESM1-2-HR, GFDL-ESM4, CESM2, CESM2-WACCM, and IPSL-CM6A-LR). An exception occurs for IPSL-CM6A-LR, which remains oligotrophic until November, and for CanESM5, which is oligotrophic throughout the year. By the end of the century, nine models reach oligotrophy (all except MPI-ESM1-2-HR and CNRM-ESM2-1), with duration extending from early May to February. However, the multi-model mean is influenced by models with higher nutrient concentrations and reaches oligotrophy only during July and August.

Figure 5 Seasonal cycle of nitrate concentration during the two last decades of the historical simulation (1995–2014) (left) and the last two decades of SSP5-8.5 (2081–2100) (right). Each model is represented in a different color while the multi-model mean is shown in black.

The multi-model mean SCM deepens throughout the century across the Arctic Ocean, while the inter-model standard deviation increases (Fig. 6). During the period 1995–2014, the SCM is below 25 m over most of the Arctic Ocean, whereas it increases to up to 50 m by the end of the century, particularly in the Beaufort Sea and the central Arctic. A slight deepening is projected along the coastal regions of the Chukchi, Siberian, Laptev, and Kara seas.

Figure 6 Depth of subsurface chlorophyll maximum (SCM) in July for the periods 1995–2014 (a) and 2081–2100 of SSP5-8.5 (b), and the difference between these two periods (c) (top row). Standard deviation of the July maximum chlorophyll depth across models for 1995–2014 (d), 2081–2100 of SSP5-8.5 (e), and the difference in standard deviation between these two periods (f) (bottom row).

Models largely disagree on the depth of the SCM, particularly in future projections where deepening is observed. In some regions, the uncertainty is nearly as high as the SCM depth, possibly indicating that while some models predict a deep SCM, others do not, consistent with findings from CMIP5 models (Steiner et al., 2016).

Across all Arctic sub-regions, multi-model mean NPP, NO₃, SSIC and SST are relatively stable until ∼ 1980, after which all basins begin to exhibit varying levels of surface ocean warming alongside NO₃ and SSIC declines (Fig. 7). The sign of regional-scale anomalies in NO₃, SSIC and SST is consistent across the CMIP6 ensemble, however model uncertainty varies across regions (Figs. 7 and 8).

Figure 7 CMIP6 multi-model mean basin-scale evolution of NPP (green), NO₃ (pink), SSIC (blue) and SST (orange), over the period 1900–2100 of the historical and SSP5-8.5 simulations, relative to 1850–1899.

Figure 8 SSP5-8.5 ensemble anomalies of NPP (green), NO₃ (pink), SSIC (blue) and SST (orange) for each Arctic basin and the entire Arctic Ocean domain. Each boxplot represents the multi-model mean, interquartile range and ensemble outliers in 2081–2100 relative to 1995–2014.

The pan-Arctic Ocean ΔNPP is relatively low compared to some other regions, with an ensemble range of −3.62 to 18.3 g C m⁻² yr⁻¹ for the entire Arctic Ocean and −2.52 to 91.3 g C m⁻² yr⁻¹ for the Chukchi Sea, where the increase and associated uncertainty is the highest. The decrease in NO₃ concentration is relatively consistent across Arctic sub-regions, ranging from approximately −2.38 in the Laptev Sea to −4.12 mmol m⁻³ in the Arctic Basin. SST increases across all subregions, with an average increase of approximately 3.66 °C across the entire Arctic Ocean. The highest multi-model ΔSST uncertainty is in the Barents Sea, with a model spread of 6.95 °C, while the lowest uncertainty occurs in Baffin Bay, with a spread of 3.42 °C (Fig. 8).

Baffin Bay differs from other regions as the only area where the multi-model mean NPP is projected to decline under SSP5-8.5. During the historical period (1995–2014), this subregion already has low annual mean sea ice coverage, therefore although it becomes ice-free by the end of the century under SSP5-8.5, the change in coverage remains small. Baffin Bay is also projected to experience the lowest multi-model mean surface ocean warming (2.53 °C) compared to other Arctic Ocean regions, with low uncertainty (Fig. 8).

3.3 Phytoplankton growth limitation terms

As phytoplankton growth is the product of multiple limiting factors (Eq. 1), the analysis of their relative anomalies throughout the 21st century is necessary to assess their respective contributions to phytoplankton growth rate perturbations and the projected increase in Arctic Ocean NPP. Diatoms exhibit a smaller increase in NPP compared to miscellaneous phytoplankton, with respective increases of approximately 11 % and 52 % relative to 1995–2014. However, the uncertainty associated with both groups is high and increases throughout SSP5-8.5 (Fig. 9a). The greatest driver of relative increases in multi-model mean phytoplankton growth rates is a reduction in light limitation (Fig. 9). Alongside this, warming directly enhances growth rates via increases in T_f although this is offset by enhanced nutrient limitation. The L_lim and T_f increase is consistent between PFTs (Fig. 9b, c). However, the ΔL_lim uncertainty increases sharply, exceeding ±30 % by 2100 for both PFTs while ΔT_f and ΔN_lim uncertainty is between ±10 % and ±20 % by 2100 (Fig. 9b, c and d). Both PFTs are increasingly nutrient limited in the Arctic Ocean over the 21st century, but diatoms are more impacted by nutrient concentration declines (Fig. 9b). These trends are consistent across models (see Fig. A5).

Figure 9 Multi-model mean SSP5-8.5 pan-Arctic Ocean anomalies of diatom (blue) and miscellaneous phytoplankton (orange) NPP (a), nutrient limitation (b), light limitation (c) and temperature function (d). Anomalies are relative to 1995–2014 values with the multi-model standard deviations shaded.

To better understand NPP uncertainties, their causes and their distribution across the Arctic, we represented the multi-model standard deviation of SSP5-8.5 anomalies in NPP and growth limitation terms (Fig. 10). Diatoms and miscellaneous phytoplankton show comparable Arctic Ocean patterns of multi-model uncertainty associated with SSP5-8.5. Relative ΔNPP uncertainty is particularly pronounced in the central Arctic Ocean (Fig. 10a and e). This is coincident with the location of greatest ΔL_lim uncertainty which reaches > 210 % over the century in this region (Fig. 10c and g). In contrast, ΔT_f and ΔN_lim exhibit much lower uncertainty, remaining below 30 % across the entire Arctic Ocean (Fig. 10b, d, f and h). Given the very low NPP values in the Central Arctic Ocean, it is not surprising that the regions of greatest relative NPP uncertainty (central Arctic Ocean) differ from those exhibiting the largest absolute NPP increases (shelf regions, Fig. 4a).

Figure 10 Regional differences in the multi-model standard deviation of anomalies in NPP (a, e), nutrient limitation (b, f), light limitation (c, g) and temperature function (d, h) in 2081–2100 of SSP5-8.5 relative to 1995–2014 of the historical. The upper row represents diatoms and the lower row miscellaneous phytoplankton.

In the two decades 1995–2014, diatoms are not the dominant contributor to total NPP across the majority of the Arctic Ocean with the exception of the coastal East Siberian Sea where diatom NPP is particularly important (Fig. 11a). Under SSP5-8.5, miscellaneous phytoplankton dominance is enhanced, notably in the central Arctic Ocean where they contribute to up to 70 % of total NPP. In contrast, diatom NPP remains dominant in the coastal East Siberian Sea and exhibits increasing dominance in the Kara and Laptev Seas (Fig. 11b).

Figure 11 Multi-model mean fraction of total NPP due to diatoms in 1995–2014 (a), and 2081–2100 under SSP5-8.5 (b). Change in % of total NPP due to diatoms (c). Panels (d) and (e) show the anomalies in miscellaneous phytoplankton NPP (d), and diatom NPP (e) in 2081–2100 compared to 1995–2014.

The NPP increase is the greatest in Chukchi and Barents seas for both PFTs but highest for miscellaneous phytoplankton (Fig. 11d, e). It confirms that the increase in total NPP in the pan-Arctic Ocean is driven mainly by enhanced productivity of miscellaneous phytoplankton, rather than an increase in diatom productivity (Fig. 4). NPP of both PFTs, and in particular diatom NPP, declines in the Baffin Bay, as well as in regions of the Nordic Seas and on the West Siberian shelf (Fig. 11d, e).

In the reference period (1995–2014), N_lim is lower for diatoms than for miscellaneous phytoplankton (Fig. 12a and b). This is likely indicative of diatoms typically having higher nutrient half-saturation constants than nanophytoplankton in most ocean biogeochemical models (e.g. Laufkötter et al., 2015). It could also be due to differences in the limiting nutrient but this is unlikely given the agreement between the spatial distributions of N_lim and ΔN_lim for each phytoplankton type (Fig. 12m and n). In the Canadian Basin and in Baffin Bay, diatom N_lim decreases by over 40 %, with changes in diatom growth rate dominated by N_lim anomalies. For miscellaneous phytoplankton, nutrient limitation is also enhanced but to a lesser extent, with a maximum N_lim decrease below 20 % (Fig. 12m and n).

Figure 12 CMIP6 multi-model mean nutrient limitation, light limitation and temperature function for each PFT in 1995–2014 of the historical simulation (top), 2081–2100 of SSP5-8.5 (middle) and anomalies between these periods (bottom).

For both PFTs, L_lim and T_f exhibit similar regional distributions during the reference period (1995–2014). L_lim is highest in the Nordic and Barents Seas, reaching 0.4, while T_f is also highest in these regions with values exceeding 1.8. Conversely, both L_lim and T_f show their lowest values in the Arctic Basin (Fig. 12c, d, e, and f). The temporal evolution of these terms is consistent between PFTs. L_lim increases most substantially in the central Arctic, with an increase of approximately 180 % relative to the reference period (Fig. 12o and p) demonstrating strong alleviation of light limitation, while T_f increases are most pronounced in the Kara Sea, Barents Sea, and Nordic Seas, exceeding 50 % relative to the reference period (Fig. 12q and r).

In SSP5-8.5 simulations, increasing L_lim (reduced light limitation) is the principal driver of multi-model mean increases in phytoplankton growth rates for both diatoms and miscellaneous phytoplankton across the Arctic Ocean, particularly in the central Arctic where this is coincident with NPP increases (Fig. 13). In the Nordic, Barents, and Kara Seas, phytoplankton growth rate changes of both PFTs are mainly driven by increasing T_f (thermal enhancement) and generally consistent with NPP increases, except off the Norwegian coast. In contrast, in Baffin Bay, declining N_lim (greater nutrient limitation) is the dominant contributor to growth rate changes for diatoms and is consistent with declining diatom NPP (see Fig. A6, A7, A8, A9). For miscellaneous phytoplankton however, increasing L_lim is the dominant contributor to growth rate changes. This is indicative of enhanced growth rates but miscellaneous phytoplankton NPP is projected to decline in Baffin Bay (Fig. 13b). This suggests that other mechanisms, such as enhanced zooplankton grazing rates (Rohr et al., 2023) may contribute to projected NPP declines in this region.

Figure 13 Dominant driver of phytoplankton growth rate changes (Δμ) in SSP5-8.5 for diatoms (a) and miscellaneous phytoplankton (b). Colors denote whether L_lim (yellow), T_f (orange) or N_lim (pink) is the largest contributors to Δμ. NPP anomalies are negative in stippled areas and positive elsewhere.

4 Discussion

Projected twenty-first century increases in Arctic Ocean NPP are higher and more variable in CMIP6 than in CMIP5. The multi-model mean NPP increase is four times larger under comparable high radiative forcing in CMIP6 than in CMIP5 (respectively 46.4 % and 11.5 % compared to the historical period, see Fig. 6a), with uncertainty at the end of the century three times higher in CMIP6 than CMIP5, consistent with previous studies (Tagliabue et al., 2021). The CMIP6 and CMIP5 ensembles also exhibit differing temporal evolution of NPP over the 21st century. All CMIP6 models exhibit monotonic increases in Arctic Ocean NPP, whereas CMIP5 models exhibit more diverse responses, with some models showing an initial increase in NPP followed by a subsequent decrease due to the emergence of oligotrophic conditions (Vancoppenolle et al., 2013). In absolute values, the most substantial NPP increases in CMIP6 are observed in the inflow shelf regions, which experience the greatest warming, sustained nutrient levels and limited change in light availability. This is consistent with the trend in the current observations (Lewis et al., 2020). The model uncertainty associated with CMIP6 Arctic Ocean NPP projections is higher than the scenario uncertainty (Fig. 3). This result is consistent with global NPP projections, which indicates finely balanced limitations on NPP, which often compensate one another (Bopp et al., 2013; Kwiatkowski et al., 2020).

The direct influence and balance among climate-driven drivers of phytoplankton growth rates, namely temperature, light and nutrient availability, shifts substantially in CMIP6 compared to CMIP5. The stronger Arctic Ocean surface warming projected in CMIP6 and greater associated uncertainty (Fig. 6) relates to higher and more variable climate sensitivity (Zelinka et al., 2020) with similar Arctic amplification simulated across both ensembles (Hahn et al., 2021). This results in greater thermally-driven increases in phytoplankton growth rates than previously simulated in CMIP5 (Laufkötter et al., 2015; Nakamura and Oka, 2019). The more pronounced Arctic Ocean sea ice loss in CMIP6 (Notz and Community, 2020) increases light availability in the upper ocean, thereby further acting to enhance phytoplankton growth rates and associated NPP. In addition to this intensified sea ice loss, a number of CMIP6 models include an improved representation of sea ice light attenuation, resulting in generally higher under-ice light transmission (Lebrun et al., 2023). Although CMIP6 simulations exhibit a larger decrease in upper Arctic Ocean nitrate concentrations over the 21st century (Fig. 6b), oligotrophic conditions are reached later and less frequently compared to CMIP5 models (Vancoppenolle et al., 2013), due to higher simulated mean state nutrient concentrations.

Beyond changes in mean nutrient concentration, CMIP6 projections point to a reorganization of the seasonal dynamics of nutrient limitation throughout the century in the Arctic Ocean. The lengthening of the oligotrophic period in CMIP6 projections suggests a progressive extension of nutrient limitation beyond the peak bloom phase, potentially affecting both the magnitude and the phenology of phytoplankton production (Manizza et al., 2023). The intensity of the bloom is constrained by the seasonal depletion of nutrients, despite enhanced light availability and warmer temperatures, while a delayed replenishment in autumn and winter could limit the recovery of nutrient inventories prior to the following growing season. Contrary to expectations, nutrient replenishment is not occurring in autumn despite the expansion of open water areas. Under usual conditions, the reduced sea ice cover in this season would increase fetch, which could allow stronger storm-driven mixing and entrain nutrient-rich deep waters into the surface layer (Ardyna et al., 2014). Our analysis suggests that if his phenomenon is occurring, it is inefficient.

Although NPP is projected to increase in most regions of the Arctic Ocean, many models simulate local declines, particularly in Baffin Bay and in the Nordic Seas. These regions are characterized by rising temperatures, reduced nutrient availability and little or no sea ice in present-day conditions. As such, the potential for changes in light supply is low in these areas and phytoplankton growth shows little to no increase due to increased light availability. While higher temperatures do stimulate growth rates, the limiting effect of declining nutrient availability dominates, resulting in overall NPP declines in Baffin Bay and Nordic Seas. The projected decline in NPP in Baffin Bay is consistent with recent observations of decreasing productivity in this region (Ardyna and Arrigo, 2020).

In addition to this projected NPP increase, in CMIP6 models, the SCM deepens in CMIP6, consistent with CMIP5 (Steiner et al., 2016). However, the uncertainty associated with this deepening is high and the modeled historical SCM depth is still inconsistent with observational estimates.

Arctic Ocean phytoplankton communities are projected to shift toward smaller size classes, with miscellaneous phytoplankton increasingly dominating primary production throughout the 21st century (Fig. 11). This change in phytoplankton community structure aligns with previous projections from ESMs (Bopp et al., 2005; Fu et al., 2016), and is driven by differential nutrient affinities among PFTs. Miscellaneous phytoplankton demonstrate greater competitive advantage than diatoms under the increasingly low-nutrient conditions (Ward et al., 2012). In most of the Arctic Basin where NPP is increasing, productivity exhibits substantial gains, with miscellaneous phytoplankton driving most of this enhancement. This production is facilitated by newly ice-free waters that provide increased light availability and extended growing seasons. In contrast, in Baffin Bay and Nordic Seas, the decline of NPP is driven by a pronounced decrease in diatom productivity. While smaller miscellaneous phytoplankton can efficiently exploit low-nutrient environments and capitalize on higher temperatures, larger diatoms require higher nutrient concentrations to maintain growth (Marinov et al., 2010). These contrasting patterns suggest a fundamental shift in Arctic marine ecosystems, with implications for higher trophic levels and thus carbon pump efficiency (Grebmeier et al., 2010; Ward et al., 2012).

In particular, a transition toward smaller phytoplankton cells is expected to reduce the efficiency of energy transfer to zooplankton, especially large lipid-rich copepods such as Calanus glacialis and Calanus hyperboreus, whose development depends on the timing and magnitude of diatom-dominated blooms for feeding and reproduction (Leu et al., 2011; Wassmann, 2011). Such temporal mismatches can cascade through the food web, affecting the growth, recruitment and distribution of key predators including fish, seabirds and marine mammals (Haug et al., 2017). Moreover, although export production is not assessed in this study, reduced export of large, fast-sinking diatom aggregates may weaken pelagic–benthic coupling, ultimately altering benthic communities and their role in carbon remineralization (Wassmann et al., 2006).

This study highlights the complexity of projecting future NPP in the Arctic Ocean under anthropogenic pressure. While Arctic Ocean NPP is projected to increase throughout the 21st century, substantial uncertainties persist alongside pronounced spatial heterogeneity, with diverse impacts. The projected shift in phytoplankton community structure toward smaller taxa may exacerbate the weakening of the Arctic Ocean carbon sink, which is already projected to decline due to surface ocean warming that reduces CO₂ solubility and enhances stratification (Oziel et al., 2025). Concurrently, these ecosystem changes could cause unclear impacts to local communities, whose economic livelihoods, food security, public health, and cultural identity depend heavily on marine resources (Malik and Ford, 2025; Mudryk et al., 2021). The magnitude and spatial distribution of future NPP changes remain highly uncertain due to differences in key limiting factors, biogeochemical parameterizations and their impact on phytoplankton community structure. In particular, the limited representation of sea-ice algae in most CMIP6 models might affect the phenology of pan-Arctic NPP as well as its vertical distribution and associated carbon export. Specifically, uncertainties in light penetration through varying sea ice thickness and optical properties, coupled with incomplete understanding of nutrient availability and cycling beneath ice-covered waters, represent critical knowledge gaps. Improving future Arctic Ocean NPP simulations requires enhanced representation of present-day nutrient levels, more accurate light transmission parameterizations through sea ice, and reduced climate sensitivity uncertainties across model ensembles.

5 Conclusion

The increase in Arctic Ocean NPP projected by CMIP6 models for the end of the 21st century under high anthropogenic emissions is about four times higher than in CMIP5. This is mostly due to the higher climate sensitivity of CMIP6 models, leading to greater warming and earlier September sea-ice loss and a longer open-water season, together with a lower limitation by nutrients. CMIP6 models also incidentally simulate higher nutrient levels than their CMIP5 counterparts. The NPP increase is accompanied by a shift toward smaller phytoplankton species, driven by projected nutrient decline over the 21st century – despite CMIP6 models simulating higher absolute nutrient levels than CMIP5. This community shift has implications for higher trophic levels, biogeochemical cycling, and carbon export. The implications for ecosystem services and indigenous communities that depend on marine resources remain largely unassessed.

Appendix A: Analyses of individual models

Figure A1 Projected regional annual mean NPP during the period 1900–2100 (historical and SSP5-8.5). Multi-model mean is shown in black.

Figure A2 Projected regional annual mean MLD NO₃ concentration during the period 1900–2100 (historical and SSP5-8.5). Multi-model mean is shown in black.

Figure A3 Projected regional September sea ice concentration during the period 1900–2100 (historical and SSP5-8.5). Multi-model mean is shown in black.

Figure A4 Projected regional SST during the period 1900–2100 (historical and SSP5-8.5). Multi-model mean is shown in black.

Figure A5 Evolution of SSP5-8.5 NPP (a), N_lim (b), L_lim (c) and T_f (d) for each model and each PFT, relative to the period 1995–2014 averaged over the pan-Arctic Ocean.

Figure A6 Projected NPP relative to 1995–2014 during the period 2015–2100 of SSP5-8.5 in each subregion of the Arctic Ocean. Diatoms are represented in blue and miscellaneous phytoplankton in orange. Shading represents the multi-model standard deviation.

Figure A7 Projected nutrient limitation factor relative to 1995–2014 during the period 2015–2100 of SSP5-8.5 in each subregion of the Arctic Ocean. Diatoms are represented in blue and miscellaneous phytoplankton in orange. Shading represents the multi-model standard deviation.

Figure A8 Projected light limitation factor relative to 1995–2014 during the period 2015–2100 of SSP5-8.5 in each subregion of the Arctic Ocean. Diatoms are represented in blue and miscellaneous phytoplankton in orange. Shading represents the multi-model standard deviation.

Figure A9 Projected temperature function relative to 1995–2014 during the period 2015–2100 of SSP5-8.5 in each subregion of the Arctic Ocean. Diatoms are represented in blue and miscellaneous phytoplankton in orange. Shading represents the multi-model standard deviation.