Effects of sea ice model characteristics on Arctic marine primary production

The Arctic environment is rapidly changing and drastic declines in sea ice extent and thickness have been observed over the past decades. The thinning of the ice pack affects its dynamics and makes it more susceptible to form fractal features such as leads and ridges at scales of meters to kilometres, as well as larger openings known as polynyas of tens to tens of thousands of square kilometres. Observations show that their presence influences polar phytoplankton communities and productivity, pointing at the importance of representing such sea ice features in biogeochemical ocean modelling. The following study investigates how Arctic primary production and carbon export simulated with the biogeochemical model HYCOM-ECOSMO are affected by different sea ice models. Three simulations are available; one coupled to the Los Alamos sea ice model (CICE) with an elastic-viscous-plastic rheology. Two other simulations were run with the next generation sea ice model (neXtSIM) with a brittle rheology that better represents features such as leads and ridges. The neXtSIM(NOlight) simulation was run with the standard neXtSIM code and does not allow for light flux through the ice. In the neXtSIM(light) simulation, the light flux was parameterised to go through the young ice class with the objective to investigate the effect of refreezing leads on phytoplankton growth. Sea ice properties and key biogeochemical fields were analysed for the Chukchi Sea. We found that the choice of sea ice model rheology caused differences in the spring-bloom timing via sea ice concentration (SIC) and associated light availability. This further affected the timing of spring carbon export. The light flux through ice in CICE and neXtSIM(light) increased the light availability particularly during late spring and summer, and affected SIC in neXtSIM(light) compared to neXtSIM(NOlight). These findings suggest that the choice of sea ice model affects simulated primary production in HYCOM-ECOSMO, and that a better representation of mechanical damage in sea ice models is relevant for primary production and carbon export estimates in Arctic Ocean biogeochemical models.

Under climate change, the Arctic environment is observed to warm much faster than the rest of the globe, a phenomenon known as Arctic amplification. The most drastic estimates suggest that the region warms four times faster than the rest of the world, which strongly affects the terrestrial and marine environment (Rantanen et al., 2022). The past 47 years of available satellite observations show a long-term negative trend in sea ice cover, and Arctic sea ice extent hit several record lows over this period, with the lowest maximum extent on record in March 2025 (Meier et al., 2025). Average winter sea ice thickness shows a decrease from ∼3.6 m in the 1980s to currently ∼2 m, largely due to the loss of thicker multi-year ice, which survives at least one summer season (Kwok, 2018).

Primary production (PP) is limited by photosynthetically active radiation (PAR), which is the spectrum of short-wave radiation that can be utilised as source of energy by phytoplankton (Mõttus et al., 2012) and the availability of nutrients. In the Arctic Ocean, sea ice conditions significantly regulate the light environment and influence upper-ocean stratification, vertical mixing and thus nutrient supply (Arrigo and Van Dijken, 2011; Rudels and Carmack, 2022). Once light becomes available in the water column with the end of the polar night at the beginning of spring, and the mixed layer starts to shoal, primary production rates increase exponentially and peak in late spring. This is referred to as the spring-bloom, which is currently mostly driven by the growth of diatoms (Ardyna and Arrigo, 2020). Mixed layer shoaling determines the distribution and therefore light exposure of phytoplankton in the water column. Sverdrup (1953) hypothesised that the spring-bloom initiates once the mixed layer shoaling transports phytoplankton above a critical depth between 30 to 300 m, where the light exposure is sufficient for phytoplankton growth to exceed losses. However, it has been shown that the distribution of grazers and concentrations of pathogens can vary significantly within the water column due to vertical mixing and dilution, which further effects the concentrations of phytoplankton in the euphotic zone during the spring-bloom onset (Behrenfeld, 2010; Hoppe et al., 2024). While the availability of light triggers the start of the growth season, or may even limit growth underneath the pack ice (Clement Kinney et al., 2023), nutrient inventories can control the overall magnitude of primary production by setting the trophic status of an area (Tremblay et al., 2015). Nitrogen has been observed and simulated to be the limiting nutrient in the group of macro-nutrients in the Arctic Ocean (Clement Kinney et al., 2023; Tremblay and Gagnon, 2009). The commonly cited pathways for the replenishment of macro-nutrients (hereafter referred to as nutrients) in the euphotic zone are advection, vertical mixing and upwelling, as well as in-situ remineralisation by heterotrophic organisms (Tremblay and Gagnon, 2009; Tremblay et al., 2015; Terhaar et al., 2021). Furthermore, Terhaar et al. (2021) estimate that 29 to 51% of current annual Arctic net primary production could be sustained by terrigenous nitrogen fluxes from rivers and coastal erosion.

In-situ observations of Arctic Ocean primary production, mainly available after 1975, show a wide range of values between 17.5 and 5,260 mgC m -2 d -1 (Lee et al., 2016), and remote sensing-based estimates report large differences in production rates between Arctic regions. For example, Arrigo and Van Dijken (2011) estimated mean annual PP during the period of 1998 to 2009 at 96.9 gC m -2 in the Chukchi Sea, values over 100 gC m -2 for Greenland, Barents, Kara and Siberian seas, and ∼70 gC m -2 of mean annual PP in Beaufort Sea and Baffin Bay, including the most northern latitudes of each region into their estimates. However, other studies have shown that the inflow shelves are hot spots of primary production, while the Arctic central basin’s annual productivity is very low due to the effects of persistent sea ice on light, and low nutrient levels (Codispoti et al., 2013; Hill et al., 2013; Sakshaug, 2004). Together with the decline in sea ice, satellite observations show that the Arctic Ocean net primary production has increased by 57% since 1998 (Lewis et al. 2020), and 19.2% between 2003 and 2024 (Frey et al. 2024). This trend is mainly attributed to increases in open water area and duration from sea ice decline, and increases in Chlorophyll a concentrations that must be sustained by increasing nutrient availability (Lewis et al., 2020). The open water duration is leading to an earlier spring-bloom onset and longer growth season (Ardyna and Arrigo, 2020). The estimated rates of change in PP vary across the Arctic, and even a decline in PP has been documented (Arrigo and Van Dijken, 2011; Lewis et al., 2020). Negative trends were detected mostly on the interior and outflow shelves, and have been associated with increasingly nutrient-depleted outflow waters that supported increasing PP ‘upstream’, or increasing light and nutrient limitation due to larger inflow volumes of turbid and nutrient-poor river runoff (Arrigo and Van Dijken, 2015; Lewis et al., 2020). How Arctic Ocean PP will continue to change throughout the 21 st century remains uncertain. Some models are projecting an increase of net PP mainly associated with the decrease in sea ice extent and the related increase in light availability (Champiot-Bayard et al., 2025). Productivity may further be enhanced by higher influx of terrigenous nutrients from increasing river flows or new sources such as the melt of glaciers (Oziel et al., 2025; Terhaar et al., 2021). In contrast, enhanced meltwater-driven stratification could potentially inhibit the replenishment of nutrients from deeper water layers (Von Appen et al., 2021).

A more seasonal ice pack that is increasingly made up of first-year ice (Kwok, 2018) may experience earlier spring melt (Bliss and Anderson, 2018) and could strongly affect the timing and rate of biological productivity. Additionally, a shift to more mobile (Rampal et al., 2009), thinner ice results in dynamical changes to the consolidated pack ice. The sea ice becomes more susceptible to form fractal features such as leads and ridges at scales of meters to kilometres, as well as larger openings of tens to tens of thousands of square kilometres known as polynyas (Rampal et al., 2009). These openings may significantly change the light environment in the water column. It was previously observed that light availability in a pack ice area of high lead fraction was sufficient to trigger and sustain an under-ice phytoplankton bloom, due to light transmission through thin ice with thin snow cover in refreezing leads (Assmy et al., 2017). Furthermore, brine rejection from refreezing of leads during winter was shown to increase upper ocean salinity and vertical mixing (Morison et al., 1992; Smith and Morison, 1998), while wind forcing can shape upwelling and downwelling dynamics inside leads (Bourgault et al., 2020; Lee et al., 2025) that may affect nutrient availability. Furthermore, the downward transport of organic detritus is enhanced by vertical mixing, which thereby contributes to carbon sequestration (i.e. carbon export) (Oziel et al., 2025). At the same time, leads may influence the timing of mixed layer shoaling through increased meltwater influx into the upper ocean. Vivier et al. (2016) found that under mild atmospheric conditions, lead openings triggered an earlier onset of basal melt compared to surface melt in correspondence with higher water temperatures likely caused by the higher solar radiation influx through the openings. This may affect the distribution and therefore light exposure of phytoplankton, with the potential to alter the dynamics and magnitude of the spring-bloom (Hoppe et al., 2024). Considering the mechanisms described above, it is assumed that these pack ice features could strongly influence polar phytoplankton communities and productivity away from the ice edge, and could eventually alter carbon export dynamics.

Predicting future changes in Arctic marine primary production and carbon export in a changing climate requires ocean models that realistically represent sea ice characteristics and their effects on phytoplankton growth. Fractal pack ice features have been observed to influence phytoplankton growth through their effects on light availability, and are known to affect upper-ocean stratification. Yet, resolving these structures explicitly requires high model resolution, which increases computational costs. This raises a central question: to what extent does representing such fine-scale sea-ice features actually affect the modelling of Arctic primary production and carbon export? To answer this question, we investigated three coupled sea ice-ocean-biogeochemical model simulations with different configurations for the sea ice. Differences between the simulation set-ups are described in the Materials and Methods in section 2, which also introduces the variables that were selected for analysis. We selected the Chukchi Sea for regional analysis mainly because of significant sea ice cover until the late summer, which allows us to investigate whether the representation of sea ice affects the simulation of PP during the spring-bloom. Our results of this analysis are presented in section 3 and discussed in section 4, considering available observations from the Chukchi Sea.

To investigate the effects of sea ice variability on Arctic marine primary production, we analyse three simulations of a sea ice-ocean-biogeochemical model with different choices of the sea ice model and light parameterisation through sea ice. The following subsections detail the model choices and simulation setups.

Model components and configurations

The Hybrid Coordinate Ocean Model (HYCOM) used in this study is a regional ocean general circulation model that introduces hybrid vertical coordinates to better represent the physical conditions in the weakly stratified mixed layer or shallow regions (Bleck, 2002). Its domain covers the North Atlantic from about 40°N and the Arctic Ocean. The model has been configured and tested for the North Atlantic and Arctic oceans (Lisæter et al., 2007; Sakov et al., 2012; Xie et al., 2019), and the coupled system HYCOM-CICE is currently used as the operational Arctic forecast model by the European Copernicus Marine Service (E.U. Copernicus Marine Service, 2015).

ECOSMO II(CHL) is a 3D bio-physical ECoSystemMOdel that simulates hydrodynamics and an intermediate-complexity biology of the lower trophic levels, including nutrient cycles, phytoplankton and zooplankton evolution and detritus pathways. Sediment and water column exchange processes, such as resuspension, are also considered (Daewel and Schrum, 2013; Yumruktepe et al., 2022). The model simulates the trophic interactions between three phytoplankton groups (diatoms, flagellates, coccolithophores) and two zooplankton groups (micro- and meso-zooplankton). It models nutrient pathways for nitrate, ammonia, phosphate and silicate, using a fixed molar Redfield ratio of C:N:Si:P=106:16:16:1. The model also accounts for dissolved organic matter (DOM) and particulate (detrital) organic matter (Yumruktepe et al., 2022). The model is currently used for Arctic Ocean biogeochemistry analysis and forecast by the European Copernicus Marine Service (E.U. Copernicus Marine Service, 2021).

In the first simulation (CICE), the biogeochemical model HYCOM-ECOSMO was coupled to the Los Alamos sea ice model (CICE) (Hunke et al., 2015), that uses an elastic-viscous-plastic rheology (EVP). CICE is run with 5 ice thickness categories, and a salinity that changes over time, so that meltwater going into the ocean has a variable salinity. The flux of shortwave radiation and thus photosynthetically active radiation (PAR) are allowed to penetrate through the ice into the ocean in CICE as a function of the sea ice thickness. Absorption from algal pigments in the ice is parameterised.

Two other simulations were run with HYCOM-ECOSMO coupled to the next generation sea ice model (neXtSIM) (Ólason et al., 2025). neXtSIM is run on a triangular Lagrangian mesh and uses the brittle Bingham-Maxwell (BBM) rheology. Compared to the EVP rheology used in CICE, the BBM rheology better represents small-scale sea ice deformations at the same resolution. The Lagrangian framework is used to maintain a low level of numerical diffusion, to avoid losing the high gradients associated with leads after they have formed. neXtSIM does not have a sea ice thickness distribution per se like CICE, but each element of the mesh can be divided into three categories: open water, young ice, and consolidated ice. When ice begins to form from open water, it is added to the young ice category until it reaches a certain thickness, whereby it is upgraded to the consolidated ice category where it remains until it is melted. Consolidated ice never returns to the young ice category. This is the main difference with the CICE sea ice categories that are purely based on thickness, independently of the history of the ice in each category. The ice salinity in neXtSIM is fixed at 5 psu. More details are described by Ólason et al. (2025).

The difference between the two neXtSIM simulations was the parameterisation of light through ice. neXtSIM(NOlight) is run with the standard neXtSIM code and does not allow for light flux through the ice, regardless of thickness. In reality, in a lead that is refreezing, a thin layer of ice would prevent gas exchange and sensible and latent heat-flux between the ocean and air, but light could still penetrate and support primary production. The second neXtSIM simulation, neXtSIM(light), is designed to investigate the effects of refreezing in leads in early spring by letting all shortwave radiation through any young ice. It therefore represents an upper bound of the potential effects of this process. An overview of the main differences between the sea ice models is presented in Table 1. All simulations were run without accounting for fluxes in sea ice algae or nutrients between the ice and ocean. Nutrient cycling and biological activity are handled solely by HYCOM-ECOSMO.

Simulation set up, initialisation and forcing

For each of the three simulations, the domain covers the Arctic and includes the most northern regions of the Atlantic, which are relevant for simulating the Atlantic inflow through the Fram Strait. In this study, we use a horizontal resolution of 17 km on average. The model bathymetry is interpolated from the Bathymetric Chart of Oceans database (Jakobsson et al., 2024). The ocean model assumes 50 hybrid vertical coordinates, using isopycnals in the stratified open ocean and z-coordinates in the upper ocean mixed layer due to its homogeneous potential density. This method has shown to be suitable for representing the density-driven inflow of Atlantic water and deep water formation in the Arctic (Sakov et al., 2012). The upper ten depth levels are set as z-levels, in order to keep a minimum vertical resolution in areas of very shallow mixed layer depth. This is important to ensure realistic representation of the decrease in biomass production due to the light gradient in the upper ocean (Yumruktepe et al., 2022). Advection and mixing, i.e. the transport of biological state variables, is handled by the ocean model HYCOM via its native tracer-transport routines. ECOSMO is coupled to HYCOM via the open source Fortran-based Framework for Aquatic Biogeochemical Models (Bruggeman and Bolding, 2014). Therefore, they share the same spatial resolution and time-step (10 min).

A long HYCOM-CICE-ECOSMO simulation was used to spin up the ocean, sea ice and biogeochemistry. The ocean physics were initialised in 1977 with lateral boundary nudging of World Ocean Atlas 2018 (WOA2018) climatological temperature and salinity (Locarnini et al., 2019; Zweng et al., 2019). The sea ice was also initialised in 1977, based on the sea ice climatology of 1991–2020 from an Arctic Ocean Physics Reanalysis product (E.U. Copernicus Marine Service, 2020). The spin up of the biogeochemical tracers, i.e. nutrients and carbon, required less time and was therefore initialised in 1987 using World Ocean Atlas 2013 (WOA2013) monthly climatology profiles (Garcia et al., 2014).

The CICE simulation for this study was run from 1993, using restarts from the long HYCOM-CICE-ECOSMO simulation. The two neXtSIM simulations were initialised with ocean and biological fields from the CICE run on the 1 st of January 2006. Initial sea ice conditions for the neXtSIM simulations were taken from a long neXtSIM-HYCOM simulation run that was initialised in 1979 with Pan-Arctic Ice Ocean Modeling and Assimilation System (Zhang and Rothrock, 2003) and WOA2018. All simulations were run with the same lateral boundaries using nesting conditions for physics from the GLORYS12 V1 reanalysis (daily outputs at 1/12° horizontal resolution and 50 vertical levels (E.U. Copernicus Marine Service, 2018)), including sea surface height (SSH), temperature, salinity and horizontal ocean currents. Values were relaxed to WOA2018 monthly climatology. Nutrients and carbon at the boundaries continued to be relaxed toward WOA2013 monthly climatology. The analysis period started in 2007 to allow for a one-year spin-up period for the sea ice and surface ocean in the neXtSIM-simulations.

The simulations were forced with atmospheric data from the European Centre for Medium-Range Weather Forecasts’ (ECMWF) ERA5 reanalysis (Hersbach et al., 2020). The CICE-simulation was forced with 6-hourly atmospheric fields, whereas both neXtSIM-simulations were run with the default 1-hourly forcing, which is currently implemented to better catch the response of sea ice to wind. The atmospheric forcing fields include 10 m winds, 2 m air temperature, dew-point temperature at 2 m, atmospheric pressure, total precipitation and long- and shortwave radiation. River forcing is simulated with monthly climatologies derived from a combination of AHYPE3 and E-HYPE v3.1. Riverine nutrients fluxes (nitrate, phosphate, silicate) were collected from Global NEWS 2 dataset (Mayorga et al., 2010; Seitzinger et al., 2010) and simply scaled by the runoff climatology mentioned above for each river (Yumruktepe et al., 2022). Thus, inter-annual variability in river runoff and nutrient fluxes is not accounted for. An overview of the set up for all three simulations is provided in Table 1.

Study region and selected variables

We selected the Chukchi Sea for regional analysis as an important inflow region of the Arctic Ocean. Its location in the Arctic domain is shown on the map presented in Figure 1. The simulations show generally good agreement of both magnitude and seasonal variability in sea ice area when compared to OSISAF (EUMETSAT OSI SAF, 2017) (Supplementary Figure S1), and indicate that much of the region is usually covered with sea ice until late summer. All models tend to overestimate sea ice volume in the latter part of the winter when compared to CS2SMOS (Ricker et al., 2017), but still simulate a strong seasonal cycle and peak at the right time. There are some notable differences between the three simulations, making it a good candidate to investigate impacts of different sea ice conditions.

We selected the years 2012 and 2013 for in-depth analysis. 2012 is a known low sea ice year in the Arctic, and we expected these low sea ice conditions to trigger responses in GPP. Adding 2013 to our analysis allows us to investigate how effects of sea ice minima may influence the next season. The following variables were extracted for analysis: temperature, salinity, nutrient concentrations (silicate, nitrate and phosphate), mixed layer depth (MLD) as well as gross primary production (GPP) values, concentrations of the phytoplankton functional groups diatoms, flagellates and coccolithophores, concentrations of zooplankton (sum of micro- and mesozooplankton), carbon export, photosynthetically active radiation (PAR), sea ice concentration (SIC), sea ice thickness (SIT), eastward and northward sea ice velocities. All biogeochemical variables are available as the grid cell mean. Primary production can only be extracted as GPP, because the model does not have an explicit term for respiration (Yumruktepe et al., 2022). It computes GPP rates separately for each phytoplankton type, as a function of a default maximum growth rate, limiting factors for light and nutrients, and phytoplankton biomass. GPP was vertically integrated over 500 m. Carbon export is calculated as the product of detrital biomass and sinking speed and is thus a combination of mixing, vertical advection, and gravitational sinking. Values were analysed at depth of 1 m above the seabed, to also capture carbon export in areas of shallow bathymetry. Carbon export values will therefore include resuspended matter and are thus treated as an overestimation of the carbon that ends up in the sediments.

We analysed daily time series of spatially averaged variables for 2012 to 2013, to understand how the differences in sea ice models and resulting sea ice parameters affect biogeochemical conditions, and to identify the variables with the strongest effects on integrated GPP and carbon export. The visual comparison of the monthly mean sea ice conditions and biogeochemical conditions during spring-bloom (March to June) across the Chukchi Sea allowed us to further understand which factors caused the differences in GPP and carbon export in the models.

Sea ice conditions

All three model simulations in our study exhibit similar seasonal changes of SIC and approach a maximum SIC of 1 (100%) during winter (Figure 2). Maximum concentrations in CICE are similar in 2012 and 2013, while the neXtSIM simulations show slightly higher winter SIC in 2012 compared to 2013. Minimum summer concentrations are lowest in the neXtSIM(light) simulation in both years, and summer SIC are lower in all simulations in 2012 (∼0.1) compared to 2013 (∼0.3). The neXtSIM(light) simulation is the first to show a decrease in SIC in April 2012 and March 2013, but the decline in CICE is faster during June and July compared to both neXtSIM simulations.

SIC in neXtSIM is spatially more heterogeneous and smoother in CICE (Figure 3). In CICE, SIC decreases faster along the coastlines compared to the rest of the region until May in both years. Both neXtSIM simulations show a decrease in SIC around the Bering Strait in addition to a decrease along the Alaskan coast. This decrease is spatially larger in neXtSIM(light) compared to neXtSIM(NOlight). In April 2012, a larger area of SIC < 0.9 starts to appear along the Alaskan coast in all simulations, which in the following year is not present in CICE, and less prominent in both neXtSIM simulations (Supplementary Figure S4). In the northern latitudes of the Chukchi Sea, SIC < 0.9 in both neXtSIM simulations as of April indicates the presence of lead features in the pack ice. SIC in March and April remain higher in CICE compared to both neXtSIM simulations in the areas along the coastline and Bering Strait (red areas in Supplementary Figure S10). However, in May and June of both years, SIC in CICE is lower than in neXtSIM along the Alaskan coast and in the centre of the region. SIC in neXtSIM(NOlight) is higher compared to neXtSIM(light) across most of the region until May in both years.

SIT is thickest in CICE throughout 2012 and 2013, and thinnest in the neXtSIM(light) simulation (Supplementary Figure S2). Maximum SIT values are reached between April and May in all simulations. CICE simulates notably larger maximum SIT in 2012 compared to 2013 (∼3.5 m vs. ∼2.75 m), while maximum SIT in both neXtSIM simulations are only slightly thinner in 2013 compared to 2012 (between ∼2.3 m and 2.0 m). Summer minimum SIT of both neXtSIM simulations approach 0 m (2012) and ∼0.25 m (2013), while CICE remains at ∼0.25 m (2012) and ∼0.75 m (2013). CICE shows a significantly stronger thickness decrease in 2012 compared to both neXtSIM simulations, while the decrease is similar in 2013 in all simulations. Both sea ice models simulate large areas of thick sea ice (≥2 m) between March and June in both years (Supplementary Figure S3). Areas of thinner sea ice are mostly located along the Alaskan coast during those months. Particularly in both neXtSIM simulations, these areas of thinner ice correspond to where SIC decline commences (Figure 3). SIT is more heterogeneous in the higher latitudes of the region in both neXtSIM simulations compared to CICE.

Gross primary production

The onset of the spring-bloom appears earlier in the neXtSIM(light) simulation compared to the other two simulations in both years (Figure 4). Production starts earlier in 2013 compared to 2012 in all simulations. The highest production rates are reached during the months of May and June in all simulations and both years, slightly later in neXtSIM(NOlight) compared to both other simulations in 2012. A large difference between peak production rates appears between the CICE and the neXtSIM(light) simulation in both years. The maximum production rate in CICE is more than twice the production rate in neXtSIM(light) in 2012, and more than triple in 2013. Similarly, the neXtSIM(NOlight) simulation shows values for integrated GPP in 2013 that remain below 200 mgC m -2 d -1, while the peak rates in CICE reach over 400 mgC m -2 d -1. The annual mean GPP rates for Chukchi Sea are higher in 2012 compared to 2013 in all simulations (Table 2). CICE simulates the highest and neXtSIM(light) the lowest mean annual rates for GPP. Over the full simulation period of 2007 to 2014 (Supplementary Figure S8), neXtSIM(light) always has the earliest onset of primary production. Peak production rates are highest in CICE in all years (except in 2012), followed by neXtSIM(NOlight).

The spring-bloom first appears in March along the Alaskan coast in all simulations, and continues to spread across most of the region until June (Figure 5 for 2012). This spread is smoother in CICE compared to both neXtSIM simulations. The largest production values during the spring-bloom onset in March and April of 2012 and 2013 are found in neXtSIM(light), followed by CICE. The spatial extent of the spring-bloom onset in neXtSIM(light) is much larger in March 2013 (Supplementary Figure S9) compared to March 2012. During May of both years, CICE is the most productive across the entire region in comparison to both neXtSIM simulations, and productivity in neXtSIM(NOlight) is more confined to the coastal areas and the inflow region of the Bering Strait. In June 2012, neXtSIM(NOlight) shows the highest GPP values across the Siberian half of the region compared to both other simulations. The following year, productivity in June is higher in CICE and neXtSIM(NOlight) compared to neXtSIM(light), with maximum GPP values showing mostly in the south-central regions of the Chukchi Sea.

Macro-nutrients

Regionally-averaged surface concentrations of nitrate, silicate and phosphate vary greatly between simulations (Figure 6). In 2012 and 2013, surface nitrate concentrations in the neXtSIM(NOlight) simulation are approximately double of the concentration in the neXtSIM(light) simulation. The neXtSIM(light) and CICE simulations approach nitrate concentration values of 0 mmol m -3 in both summers, while minimum concentrations in neXtSIM(NOlight) remain above 0 mmol m -3. Considering these spatial averages, the region appears to be nitrate limited in summer in neXtSIM(light) and CICE. Phosphate surface concentrations follow the same concentration pathways as nitrate, with maximum concentrations over 1.1 mmol m -3 in neXtSIM(NOlight) in January 2012. Silicate surface concentrations do not show a significant replenishment of silicate during autumn 2012 and winter 2013. Throughout 2007 to 2014 (Supplementary Figure S6), winter nitrate concentrations are highest in the neXtSIM(NOlight) simulation, and lowest in the neXtSIM(light) simulation.

The highest surface concentrations of nitrate are simulated in the southwestern areas of the Chukchi Sea in March and April 2012 in all simulations (Figure 7). In 2013 (maps not shown), areas of high nitrate surface concentrations appear further south in CICE and neXtSIM(NOlight) compared to 2012. All simulations show a decrease in nitrate surface concentrations between March and June in both years, and large areas where surface concentrations approach 0 mmol m -3. The neXtSIM(NOlight) simulation shows the largest area of NO 3 > 10 mmol m -3 in all months compared to the other simulations in both years. This is in agreement with this simulation having the highest concentration values in the time series (Figure 6). Spatial distributions of silicate and phosphate (maps not shown) are similar to nitrate concentrations, without showing regions of limitation in either simulation.

Light availability

neXtSIM(light) has the earliest increase in surface light availability in March of both years; this is both earlier and higher in 2013 compared to 2012 (Figure 8). The CICE and the neXtSIM(NOlight) simulations show an increase in PAR starting in April of both years. neXtSIM(light) remains with the highest overall light availability until the beginning of June in both years. Highest maximum values of PAR during summer occur in CICE, while neXtSIM(NOlight) shows the lowest values of PAR throughout both years. Maximum PAR values in 2012 are higher in all simulations compared to 2013.

In all three simulations, PAR first reaches the sea surface earliest in the season along the Alaskan coast, and becomes available in the water column across the entire region until June with varying intensity (Figure 9). Additionally, both neXtSIM simulations show elevated PAR values around the Bering Strait as of March, and in the northern latitudes of the Chukchi Sea. The highest values of PAR at the sea surface in March to May show along the Alaskan coast. The spread of light availability across the region is smoother in CICE compared to both neXtSIM simulations, which show more spatially confined areas of elevated PAR at the sea surface. The neXtSIM(light) simulation clearly shows the largest areas of elevated PAR at the sea surface in March and April compared to both other simulations. This reflects the earlier light availability already observed in the times series in Figure 8. CICE simulates higher PAR compared to neXtSIM(NOlight) across most of the region during all months in 2012 and 2013.

Mixed layer depth

Both neXtSIM simulations simulate a deeper winter mixed layer in January to April in 2012 compared to CICE (Figure 10). The maximum depth of mixing in 2013 is similar in all simulations (∼35 m), but is reached later in CICE (April) compared to both neXtSIM simulations (January). The neXtSIM(light) simulation shows slightly deeper winter mixing at the beginning of 2012 and 2013, and a ∼5 m deeper mixed layer in December 2013 compared to the neXtSIM(NOlight) simulation. CICE simulates shoaling of the winter mixed layer in January to February of both years.

Spring shoaling starts earlier in both neXtSIM simulations in both years, and is slightly stronger in the neXtSIM(light) simulation in the end of April 2012 compared to neXtSIM(NOlight). Average summer MLD is similar in all simulations, between 5 to 10 m, but minimum depth is reached earlier in both neXtSIM simulations compared to CICE. Fall deepening starts earlier in both neXtSIM simulations compared to CICE in 2012, but has similar timing in 2013, when CICE also shows very strong temporal variability.

Across the southern extent of the Chukchi Sea, deep winter mixing mainly follows the shallow bathymetry of the region. In 2012, both neXtSIM simulations furthermore simulate deeper mixing in the Beaufort Gyre and in the northern latitudes compared to CICE (Figure 11). In 2013, CICE also simulates deeper winter mixing in the northern latitudes (maps not shown). The earlier spring shoaling of the mixed layer in both neXtSIM simulations shows in the southern latitudes in May of both years.

Carbon export

Peak export rates at 1 m above the seabed in 2012 are highest in CICE, reaching nearly 0.72 gC m -2 d -1 (60 mmol m -2 d -1) at the end of September, compared to approximately half this rate in both neXtSIM simulations (Figure 12). Peak rates are lower and more similar between simulations in 2013, around 0.24 to 0.36 gC m -2 d -1 (20 to 30 mmol m -2 d -1). During the months of January to April, both neXtSIM simulations show visible export rates, especially in January to March 2013. Mean annual carbon export rates in the Chukchi Sea for each target year are shown in Supplementary Table S1. Note that export values at 1 m above the seabed include resuspension. Due to the lack of a separate term for resuspended materials in HYCOM-ECOSMO, it is not possible to distinguish between freshly and previously deposited material. Carbon export values in this study are thus likely an overestimate of the total carbon export.

The winter carbon export in both neXtSIM simulations is mostly simulated along the edges of the Beaufort Gyre and in proximity to the Bering Strait (Figure 13). The spatial distribution of spring export is similar in all simulations and shows first along the southern extent of the Alaskan coast. Carbon is exported across a larger area in neXtSIM(light) in April compared to the two other simulations. Except in the northern latitudes, carbon export spreads across the entire region until June in all simulations.

In this section, we discuss how the differences in sea ice in the simulations affect light availability, and how the different rheologies and light parameterisations are amplifying these differences. In this context we explore how the choice of sea ice model caused the differences in primary production in HYCOM-ECOSMO, mainly the different onset times of the spring-bloom and the very low peak production rates during summer in neXtSIM. Similarly, we investigate the differences in carbon export flux at 1 m above the seabed and how they relate to GPP and sea ice characteristics in each model simulation.

Differences in sea ice and atmospheric radiation influx

Sea ice plays a key role in mediating the influx of solar radiation into the Arctic Ocean and thereby affects the availability of PAR (Arrigo and Van Dijken, 2011). Additionally, the density of openings in the pack ice together with the ice thickness regulate the heat influx into the upper ocean and can therefore influence the melt and refreezing of sea ice and thus the SIC (Maslanik et al., 2007; Perovich et al., 2001; Perovich and Polashenski, 2012).

The combined effects of fractal features and light flux through the young sea ice class made neXtSIM(light) have the earliest increase in PAR (Figure 8) with light values that stay above the other simulations until the end of May in both years. CICE showed the highest values of PAR during summer, while neXtSIM(NOlight) had lowest light availability in both years.

The effects of sea ice differences on light availability further showed when comparing maps of SIC (Supplementary Figure S10 and Supplementary Figure S11) and PAR (Figure 9 and Supplementary Figure S7). In March to May of 2012, and in April and May of 2013 the areas of lower SIC in neXtSIM(NOlight) compared to CICE co-located well with areas that showed higher light availability in the neXtSIM(NOlight) simulation, mainly along the Alaskan coast and at the Bering Strait. This indicates that during the spring months, the lower SIC and related elevated open water area (neXtSIM(NOlight)) had a strong influence on light availability compared to the light flux through ice in CICE. Ice thickness and snow cover depth have been shown to significantly control the light conditions under the sea ice in the Arctic (Stroeve et al., 2021), which may prevent significant light flux through the sea ice in CICE in early spring, as CICE showed thickest sea ice throughout 2012 and 2013 compared to the neXtSIM simulations (Supplementary Figure S2). Furthermore, the brittle rheology performs well in simulating fractal features when sea ice is thick and consolidated, characteristics that are often found during winter and early spring in regions away from the ice edge. This was also visible in the northern latitudes of the study region, which showed distinct areas of lower SIC in both neXtSIM simulations (Figure 3). All of this suggests that the brittle rheology of the sea ice model neXtSIM leads to more openings in the sea ice during early spring, and that the related effects on SIC and light availability are especially relevant for primary production modelling during the time of the spring-bloom onset. That CICE simulates higher PAR in March 2013 in areas where the simulation has higher SIC compared to neXtSIM(NOlight) may be related to the lower SIT values in 2013 compared to 2012 (Supplementary Figure S2), allowing for more light flux through the thinner sea ice in CICE that compensates for less fractal features compared to neXtSIM.

The change in light parameterisation in neXtSIM(light) compared to neXtSIM(NOlight) led to earlier increases in PAR, particularly in 2013 (Figure 8). It also resulted in slight differences in the sea ice variables. neXtSIM(light) showed lower SIT values in 2012 and 2013 (Supplementary Figure S2), and a stronger decrease in SIC during spring and summer in both years (see Figure 2). The differences in SIC were particularly pronounced in areas where the SIC decline commenced along the Alaskan coast (Supplementary Figure S10). The sea ice cover in neXtSIM(light) was clearly thinner between April and June 2012 in the same areas of SIC decrease along the Alaskan coast (Supplementary Figure S3). Refreezing leads have previously been observed to enhance sea ice melt due to the lower albedo of the thinner ice cover inside the leads, which allows for higher energy influx into the upper ocean (Perovich et al., 2001). Together with the energy influx through open leads, this may have enhanced the melt of the surrounding sea ice in nextSIM(light), and may have prevented subsequent refreezing of neighbouring leads, thus causing the lower SIT and SIC values. The daily time series of sea surface temperatures (SST, Supplementary Figure S5) supports that assumption by showing an earlier increase in SST and higher summer SST in both years in neXtSIM(light) compared to neXtSIM(NOlight).

We found the largest values of PAR (Figure 8), SIC (Figure 2) and SST (Supplementary Figure S5) in CICE during both summers, which shows that the effects of light flux through ice in CICE are amplified once the sea ice melt advances. This corresponds to the findings of Stroeve et al. (2021), who found that the light availability under sea ice at the end of the melt season was mainly controlled by SIT. However, with the current light flux parameterisation in neXtSIM(light), which allows only for light flux through the young ice class, the decrease of SIT towards summer will not change incoming radiation significantly, as most ice belongs to the consolidated ice class during the melt season. Therefore, we hypothesise that the advantages of the brittle sea ice features for light and energy influx through open and refreezing leads in neXtSIM over the more sophisticated thermodynamic simulation of light flux through varying SIT in CICE is outweighed once SIT falls below a certain value. Determining such thresholds is out of scope for this study, but our findings point at the importance of simulating both fractal features and light flux through varying thicknesses of sea ice to account for the evolution of conditioning factors for light availability under the sea ice throughout spring and summer.

Differences in primary production rates and timing

Reported annual primary production values for the Chukchi Sea vary based on the selection of the regional extent, observation technique and whether under-ice phytoplankton blooms are considered. In-situ observations report annual rates of PP for the Chukchi Sea ranging between 20 and 400 gC m -2 (Sakshaug, 2004; Lee et al., 2007; Codispoti et al., 2013). Based on satellite-observations from 1998 to 2009, Arrigo and Van Dijken (2011) estimate an annual mean PP of ∼98 gC m -2 for the Chukchi Sea (including the northern latitudes of the central basin) and identify a significant positive trend in total annual primary production over the observational period, similar to Lewis et al. (2020) who report a positive trend in PP over the past decade. Using a biogeochemical model that also simulates ice algae growth, Payne et al. (2021) modelled phytoplankton dynamics for the period of 1988 to 2018, during which under-ice blooms had frequently been reported in the Chukchi Sea. Values of annual net primary production range from ∼80 to ∼95 gC m -2 depending on whether production was dominated by under-ice algae or marginal ice zone phytoplankton blooms, or a mixed regime. Our coupled models simulated mean annual GPP rates and maximum annual rates (Table 2) within the reported ranges of in-situ observations for the Chukchi Sea, while our mean annual GPP rates are of the same order of magnitude compared to the satellite-based estimates by Arrigo and Van Dijken (2011) and model-based estimates by Payne et al. (2021).

The most notable differences between simulations in GPP (Figure 4) were the significantly earlier spring-bloom onset in neXtSIM(light) in both years and the much lower production peaks during summers of 2012 and 2013 in comparison to CICE, as well as the much lower production rates in summer 2013 in neXtSIM(NOlight).

As discussed above, the differences in sea ice models led to an earlier increase in PAR in neXtSIM(light) (Figure 8), which corresponds to the earlier onset of the spring-bloom in both years. This is expected as light availability triggers the spring-bloom (Tremblay and Gagnon, 2009). Concerning the effects of mixed layer shoaling on the spring-bloom onset, the comparison of daily time series for MLD (Figure 10) and GPP (Figure 4) showed that the spring-bloom started earlier than MLD shoaling in CICE. In addition, the monthly maps of MLD (Figure 11) compared to the maps of GPP (Figure 5) did not show an exclusive co-location of MLD shoaling and primary production onset in this simulation. In neXtSIM(NOlight), the onset of the spring-bloom roughly aligned with the beginning of the MLD shoaling in both years. The earlier onset of the spring-bloom in neXtSIM(light) appears while the model continues to simulate deeper winter mixing in the region. However, a small drop in MLD can be observed right before the spring-bloom onset, which may together with the early PAR availability trigger the earlier spring-bloom onset. Due to the different behaviours in our model simulations we assume that mixed layer shoaling is not a relevant factor that triggers the spring-bloom onset in the Chukchi Sea in our simulations. This may be due to the overall shallow bathymetry in much of the region, with an average depth of 50 m (Woodgate et al. 2005), which is close to the top end of the range of 30 to 300 m for the critical depth horizon based on the theory of Sverdrup (1953).

Nutrient availability and consequences for primary production

The availability of nutrients and associated replenishment during winter mixing is known to condition the total amount of biomass production, with nitrate as the limiting macro-nutrient in the Arctic Ocean (Clement Kinney et al., 2023; Tremblay and Gagnon, 2009). This is reflected in the daily time series of nutrient surface concentrations (Figure 6), which showed that regionally averaged nitrate concentration values approached 0 mmol m -3 in July of both years in the CICE and neXtSIM(light) simulation. It is also reflected in the regional maps of nutrient concentrations, which did not show visible areas of silicate or phosphate limitation (maps not shown), but increasing areas of nitrate surface concentrations approaching 0 mmol m -3 until June in both years and all simulations (Figure 7). This points to the spatial distribution of nitrate concentrations in the region as a potential cause of the lower biomass production in neXtSIM(NOlight) in 2013 and neXtSIM(light) in both years.

Nitrate surface concentrations appeared overall lower in the Alaskan sector of the Chukchi Sea compared to the Siberian sector (Figure 7). This may be related to circulation patterns of nutrient-rich Pacific inflow water through the western Bering Strait (Codispoti et al., 2013), which has recently been shown to flow along the Siberian Coast and recirculate back around Wrangel Island (location marked in Figure 7) into the Chukchi Sea, contributing to the high nutrient concentrations east of the island (Clement Kinney et al., 2022). Overall nitrate concentrations are lower in neXtSIM(light) compared to both other simulations. Where PAR (Figure 9) becomes available and triggers the spring-bloom onset along the Alaskan coast and at the Bering Strait in neXtSIM(light) in March 2012, nitrate surface concentrations already approach 0 mmol m -3 in the following month. This suggests that nitrate availability is insufficient to sustain substantial primary production in areas with light availability. In contrast, the CICE and neXtSIM(NOlight) simulation can utilise their much higher nitrate concentrations to sustain PP later in spring and summer. The same combination of factors leading to lower production rates in neXtSIM(light) was observed in 2013.

The significantly lower GPP rates in neXtSIM(NOlight) in 2013 cannot be explained with overall low nitrate surface concentrations, as the simulation showed the highest concentration values out of all three simulations (see times series in Figure 7 and maps in Figure 6). Concentration maxima appear in proximity to the Bering Strait and in the more south-central areas of the Chukchi Sea, where light availability is low due to higher SIC (and no radiation flux through ice) between March and May (Supplementary Figure S4). This suggests that the growth conditions during early spring were more inhibited by light compared to the previous year and compared to the CICE simulation.

Nitrate concentration differences are largest between the two neXtSIM simulations (Figure 6), while their set up is the same except for the parameterisation of light flux. The low nitrate concentrations, which are likely the cause of the lower production rates in neXtSIM(light) may be caused by differences in vertical mixing and replenishment during winter, or differences in horizontal advection (Tremblay and Gagnon, 2009; Tremblay et al., 2015). As we did not see a significant difference in the winter MLD between the neXtSIM simulations (Figure 10), it is unlikely that the low nitrate concentrations in neXtSIM(light) are related to differences in vertical mixing. The horizontal advection of nutrients through the Bering Strait may vary between simulations. On the one hand, differences in SIC and related open water area can affect sea ice mobility, with a potential influence on vertical and horizontal velocities. On the other hand, nutrient advection may be different due to consumption in the source regions. The earlier and higher levels of light availability in neXtSIM(light) and related changes in GPP may have also caused an earlier and stronger depletion of nutrients through primary production south of the Bering Strait, thus leading to a downstream lack of nitrate in the Chukchi Sea. This assumption is furthermore supported by the overall lower nutrient availability in 2013 compared to 2012 (Figure 6), considering that the low sea ice year 2012 may have led to higher biomass production and related depletion of nutrients elsewhere. The analysis of lateral nutrient fluxes is beyond the scope of the paper, but the available results indicate that the representation and analysis of biogeochemical processes beyond the Arctic gateways are important to further study the effects of sea ice characteristics on Arctic primary production.

Carbon export

Observations of particulate organic Carbon (POC) export in the Chukchi Sea by Moran et al. (2005) vary from ∼3 to ∼10 mmol m -2 d -1 between spring and summer, and Lepore et al. (2007) report an average POC export flux at 50 m during spring and summer of 20 mmol m -2 d -1. Maximum carbon export rates at 1 m above the seabed in our simulations (Figure 12) exceed these reported values. Furthermore, the mean annual export values (Supplementary Table S1) in our simulations are higher than their mean annual GPP values (Table 2), which likely reflects resuspension in the shallow shelf regions counted towards the carbon flux at 1 m above seabed. neXtSIM(light) showed an earlier increase of carbon flux in spring compared to the other simulations (Figure 12), which aligns with the earlier spring-bloom onset in this simulation (Figure 4). This suggests that the differences in sea ice models and their effects on spring-bloom onset are also affecting the timing of spring carbon export. Both neXtSIM simulations showed carbon export during January to March in both years, which is likely related to winter-time resuspension. Within the scope of our study and the analysed parameters, it was not possible to identify which factors led to the carbon export during winter in neXtSIM, and what caused the discrepancies between total biomass and total export flux in the simulations. Further investigation involving more hydrographic and biological parameters would be needed to explain the mechanisms of carbon export in the Chukchi Sea in the available simulations.

Our comparison of the HYCOM-CICE-ECOSMO with both HYCOM-neXtSIM-ECOSMO simulations showed that the difference in sea ice model rheology caused differences in the spring-bloom timing via SIC and associated light availability. The low nitrate availability in neXtSIM(light) seems to have caused the low maximum GPP rates compared to the other simulations during both years. The low maximum production rates in neXtSIM(NOlight) in 2013 however are associated with unfavourable light conditions during early spring due to high SIC. The effects of the brittle rheology and its influence on PAR in the neXtSIM simulations appeared most influential during early spring when SIC was higher, while the more sophisticated thermodynamics and related light flux through ice in CICE became more important for light availability during summer. The higher energy influx in neXtSIM(light) compared to neXtSIM(NOlight), due to light penetrating through the young ice in refreezing leads, may have contributed to lower SIC and SIT, with an additional influence on light availability and spring-bloom timing. Together with the potential upstream effects of nutrient depletion, this points at the importance of a good representation of both light flux through ice and sea ice dynamics for estimating total biomass production in the Chukchi Sea and other sea ice covered regions. The sea ice affected the timing of the carbon export, which corresponded to the onset of the spring-bloom in all simulations. We recommend to look at a region with deeper bathymetry to analyse carbon export rates in the water column of the open ocean to reduce the effect of resuspension.

We recommend a detailed analysis of regions with high SIC throughout the year to determine whether the effects on light availability from fractal features (neXtSIM) are outweighed by the sophisticated light flux through sea ice (CICE) once a certain SIT is surpassed. Snow depth should be considered during such analysis to account for the changing albedo of sea ice with snow cover, particularly during the beginning of the melt season. Furthermore, we recommend to analyse flow dynamics, currents and lateral nutrient transport with respect to sea ice differences to better understand the dynamics of nutrient distribution and advection into the Chukchi Sea, as well as effects on carbon export and resuspension. In addition, we propose to carry out an analysis of biogeochemical variables across neighbouring regions, to determine differences in upstream nutrient depletion.

The results of this study show that the differences in sea ice characteristics caused by differences in sea ice model rheologies influence primary production modelling in the biogeochemical ocean model HYCOM-ECOSMO. This highlights the importance of employing sea ice models that simulate both the sea ice cover with its fractal features and the light flux through ice. Considering the anticipated changes in sea ice towards a thinner and more mobile ice cover, improving the representation of both processes in sea ice models will improve simulations of future primary production and carbon export, and the already observed changes in bloom dynamics such as an earlier spring-bloom onset and under-ice blooms in the Arctic Ocean.