On the Role of Biogeochemical Coupling Between Sympagic and Pelagic Ecosystem Compartments for Primary and Secondary Production in the Barents Sea

Primary production in the Arctic marine system is principally due to pelagic phytoplankton. In addition, sea-ice algae also make a contribution and play an important role in food web dynamics. A proper representation of sea-ice algae phenology and the linkage with the pelagic and benthic systems is needed, so as to better understand the ecosystem response to warming and shrinking ice cover. Here we describe the extension of the biogeochemical model ECOSMO II to include a sympagic system in the model formulation, illustrated by implementation in the Barents Sea. The new sympagic system formulation includes four nutrients (NO3, NH4, PO4, and SiO2), one functional group for sea-ice algae and one detritus pool, and exchanges with the surface ocean layer. We investigated the effects of linkage between the three systems (sympagic, pelagic, and benthic) on the ecosystem dynamic; the contribution of the ice algae to total primary production; and how the changes in ice coverage will affect the lower trophic level Arctic food-web dynamics. To solve the scientific and technical challenges related to the coupling, the model was implemented in a 1D application of the General Ocean Turbulence Model (GOTM). Results showed that the model simulated the seasonal pattern of the sympagic components realistically when compared to the current knowledge of the Barents Sea. Our results show that the sympagic system influences the timing and the amplitude of the pelagic primary and secondary production in the water column. We also demonstrated that sea-ice algae production leads to seeding of pelagic diatoms and an enhancement of the zooplankton production. Finally, we used the model to explain how the interaction between zooplankton and ice algae can control the pelagic primary production in the Barents Sea.

Primary production in the Arctic marine system is principally due to pelagic phytoplankton. In addition, sea-ice algae also make a contribution and play an important role in food web dynamics. A proper representation of sea-ice algae phenology and the linkage with the pelagic and benthic systems is needed, so as to better understand the ecosystem response to warming and shrinking ice cover. Here we describe the extension of the biogeochemical model ECOSMO II to include a sympagic system in the model formulation, illustrated by implementation in the Barents Sea. The new sympagic system formulation includes four nutrients (NO 3 , NH 4 , PO 4 , and SiO 2 ), one functional group for sea-ice algae and one detritus pool, and exchanges with the surface ocean layer. We investigated the effects of linkage between the three systems (sympagic, pelagic, and benthic) on the ecosystem dynamic; the contribution of the ice algae to total primary production; and how the changes in ice coverage will affect the lower trophic level Arctic food-web dynamics. To solve the scientific and technical challenges related to the coupling, the model was implemented in a 1D application of the General Ocean Turbulence Model (GOTM). Results showed that the model simulated the seasonal pattern of the sympagic components realistically when compared to the current knowledge of the Barents Sea. Our results show that the sympagic system influences the timing and the amplitude of the pelagic primary and secondary production in the water column. We also demonstrated that sea-ice algae production leads to seeding of pelagic diatoms and an enhancement of the zooplankton production. Finally, we used the model to explain how the interaction between zooplankton and ice algae can control the pelagic primary production in the Barents Sea. Keywords: biogeochemical model, primary and secondary production, Barents Sea, zooplankton grazing impact, sea-ice algae

INTRODUCTION
The Arctic Ocean is one of the environments experiencing the most severe changes due to ongoing climate change (Anisimov and Fitzharris, 2001;Bennett et al., 2015), exemplified by the extent of summer sea ice retreat (Stroeve and Notz, 2018). The current rate of decrease in summer ice cover is estimated to be 13.5% (Peng and Meier, 2018) per decade, with commensurate changes in ice thickness (Lindsay and Schweiger, 2015). Forecasts predict an Arctic Ocean completely free of ice within 30 years (Wang and Overland, 2012) and possibly by 2030 (Diebold and Rudebusch, 2019). These factors all contribute to an earlier exposure of the ocean surface to the atmosphere in spring and ice-free conditions over a longer season (Stroeve et al., 2012;Meier, 2017). All of these changes impact on the physical and ecological characteristic of the Arctic Ocean (Darnis et al., 2012), and its role in the air-sea interface interaction (Cottier et al., 2017;McPhee, 2017).
Sea ice is an essential component of Arctic ecosystems, acting both as a physical barrier to light, wind and dispersion (Clark et al., 2017;Bouchard et al., 2018), and as a substrate for living organisms. The presence of ice modifies and influences exchanges between the atmosphere and the water column such as nutrient deposition, wind stress, water mixing, and light availability, which in turn affects primary production. Moreover, sea ice represents a habitat for a rich microorganism community, mainly composed of microalgae (Thomas, 2017). Ice algae communities play an important role in Arctic ecosystems by contributing to primary production. Sympagic production contributes 3-25%, and in some regions such as the central Arctic, more than 50% of the total primary production (Legendre et al., 1992;Gosselin et al., 1997;Kohlbach et al., 2016). Ice algae also represent a valuable resource for zooplankton grazers, such as copepods (Runge and Ingram, 1991;Leu et al., 2011), amphipods (Werner, 1997), or krill (Stretch et al., 1988;Kohlbach et al., 2017a). In the Arctic shelves, ice algae can also be a key resource early in the season for benthic organisms when they settle to the seabed after ice melt (Søreide et al., 2013;Schollmeier et al., 2018). Ice algae are important for the phenology of pelagic and benthic organisms which rely on this additional energy source for reproduction and growth in spring (McMahon et al., 2006;Søreide et al., 2010), and have a key role in the benthic communities' structure and dynamic (Kohlbach et al., 2019).
An especially productive Arctic region is the Barents Sea (Sakshaug et al., 2009). In the Barents Sea, sea-ice algae can locally represent around 20% of the total primary production (Hegseth, 1998). However, winter sea ice extent in the Barents Sea has substantially decreased especially since around 1990 (Årthun et al., 2012). Impacts of the sea ice decrease on the primary and secondary production phenology and amplitude remain poorly understood in the Arctic and in the Barents Sea. Some studies already highlight a potential mismatch between ice algae production and zooplankton production and suggest important consequences for the reproduction of copepods (Leu et al., 2011;Dezutter et al., 2019). For example, the reproduction of Calanus glacialis has been found to coincide with the spring ice algae bloom while hatching of eggs coincides with the phytoplankton bloom (Søreide et al., 2010). In the Barents Sea around 80% of zooplanktonic biomass consists of copepod species, of which most, such as C. glacialis and Calanus finmarchicus, can actively feed on sea-ice algae attached to the ice Ingram, 1988, 1991), as well as on ice algae released to the water column during the melting process (Michel et al., 1996). While living ice algae cells are important for the zooplankton, dead cells are of relevance for benthic organisms. Few ice algae can survive in the water column and most of them are rapidly exported to the sea floor after ice melt where they contribute to benthic production (Ambrose et al., 2005;Juul-Pedersen et al., 2008;Boetius et al., 2013). Rich in polyunsaturated fatty acids, ice algae represent a valuable food for benthos and play a key role in the detritivorous feed web (McMahon et al., 2006;Sun et al., 2007). For example, Søreide et al. (2013) showed that benthos principally uses particulate organic matter originating from ice algae in the Svalbard region and highlighted that changes in the ice algae production can affect the pelagic and benthic systems. In this context it is essential to thoroughly understand sea-ice algae phenology and its linkage with the pelagic and benthic ecosystems when aiming to identify future changes in the Arctic ecosystem.
Because the Arctic is an extreme environment, in situ observations are not readily available and are often limited to the spring-summer season. Under these circumstances numerical models become a tool of choice for studies covering large spatial and temporal scales. In the last 20 years, there have been important improvements in the modeling of sea ice biogeochemistry (Steiner et al., 2016;Vancoppenolle and Tedesco, 2017). First models were developed and parametrized in a 1D framework focusing on specific locations (Arrigo et al., 1993;Lavoie et al., 2005) and particular processes, for example, the complex brine dynamic and nutrient exchange across the ice-ocean interface (Vancoppenolle et al., 2010). These 1D studies were supported by improvements in understanding of the biogeochemistry and dynamic of the Arctic ecosystem (Yakubov et al., 2019). The studies highlighted instances where ice algae could seed pelagic phytoplankton growth and modify the timing of pelagic primary production (Jin et al., 2006(Jin et al., , 2007Tedesco and Vichi, 2010), and how ice algae could control pelagic phytoplankton biomass by competing for nutrients in the surface water (Mortenson et al., 2017). Subsequently, sea ice modeling has been extended to a 3D framework simulating largescale primary production in the whole Arctic basin (Deal et al., 2011) or more regional area such as the Hudson Bay System (Sibert et al., 2010). In spite of these recent improvements in the biogeochemical modeling, the linkages between the sympagic, pelagic and benthic systems are still poorly represented.
In order to improve knowledge and understanding of the linkage between the sympagic, pelagic and benthic ecosystems, we present here an enhanced version of the biogeochemical model ECOSMO II (Daewel and Schrum, 2013) implemented for the Barents Sea, to include a fully coupled sympagic system. With recent developments of the ECOSMO-E2E (Daewel et al., 2019), the model will then later allow us to explore ecosystem dynamic including predators such as fishes and macrobenthos. These two components (fish and macrobenthos) are dependent on high-quality sea-ice algae food source (Sun et al., 2009;Søreide et al., 2013) and should be linked to the sympagic system to accurately represent Arctic ecosystem dynamic. We used the model in a 1D application to comparatively assess three characteristic regions of the Barents Sea and investigated (1) the importance of linkages between the pelagic and the sympagic systems for ecosystem dynamics; (2) the relative contribution of the ice algae on the total primary and secondary production; and (3) how variability in ice coverage affects Arctic foodweb dynamics.

MATERIALS AND METHODS
We present here a new version of the ECOSMO II (Daewel et al., 2019) model including a new sympagic (sea-ice biogeochemistry) module and its application to the Barents Sea region. The biogeochemical model was built with the Fortran-based Framework for Aquatic Biogeochemical Models (FABM) (Bruggeman and Bolding, 2014) in order to facilitate coupling to the physical model described below (section Physical model and forcing data). Here we work in a 1D numerical framework, which allows model parameterization, verification, sensitivity tests to study process level with a low computation effort without add complex hydrodynamic processes. Therefore, this approach neglects horizontal transport and only take into account vertical exchange processes.

Physical Model and Forcing Data
Physical processes in the water column are calculated by the 1D General Ocean Turbulence Model (GOTM) (Burchard et al., 1999). The simulation resolves the profiles of velocities, temperature, salinity, turbulent mixing, and transport of ecosystem state variables in 50 vertical layers. The physical environment affects the ecosystem dynamics, while an ecosystem feedback to the physics, e.g., heating or albedo is switched off. Photosynthetically available radiation is calculated based on the ecosystem tracers in water and ice as described below. At the surface, the model is forced with atmospheric conditions from the NCEP-R2 reanalysis [NCEP-DOE AMIP-II Reanalysis (R-2); Kanamitsu et al., 2002], and sea ice thickness, sea ice concentration and snow thickness are prescribed from a hindcast simulation with the HYCOM (HYbrid Coordinate Ocean Model) (Bleck, 2002) using a setup based on a model version described in Samuelsen et al. (2015). Tidal currents are prescribed M2 and S2 components of the horizontal velocity in order to parameterize tidal mixing and a spring-neap cycle. Tidal variability of the sea surface elevation is neglected for the present study. Physical conditions are initialized from the World Ocean Atlas v2013 (Boyer et al., 2013), which is also used for a filtered restoring of salinities to the climatological conditions in order to parameterize horizontal fluxes of fresh water. The ecosystem model dynamics are coupled with the FABM framework (Bruggeman and Bolding, 2014) into the 1D model, integrated with a third-order Runge-Kutta ODE solver (Schober et al., 2014), and transported in the vertical due to sinking, vertical advection, and turbulent mixing. The 1D model framework provides a consistent scheme for the analysis of ice-ocean ecosystem exchange fluxes, while horizontal triggers are included only through nudging of salinities and forcing with ice conditions.

Biogeochemical Sea Ice Processes
The ECOSMO II model (Schrum et al., 2006;Daewel and Schrum, 2013), including pelagic and benthic biological and geochemical components, has been adapted here for the Arctic ecosystem. One of the most important characteristics of the Arctic ecosystem is the presence of a sea ice cover and the associated microbial communities, forming a sympagic system (Thomas, 2017). To make the existing ecosystem model applicable to the Arctic, a sympagic system model is required that allows online coupling to the existing model for the pelagic and benthic systems. To achieve this, we added six state variables to the original ECOSMO II model framework, accounting for: one sea-ice algae functional group (IA), four nutrients variables (phosphate PO 4 , nitrate NO 3 , ammonium NH 4 , and silicate SiO 2 ), and one detritus variable (ID) (Figure 1).
Arctic sea ice supports a highly diverse microbial community (Arrigo, 2014;Thomas, 2017). Many species of algae bloom in succession throughout the annual season (Hegseth, 1992;Ambrose et al., 2005). Several studies have emphasized the dominance of diatoms (Werner et al., 2007), especially Nitzschia frigida and Melosira arctica (Horner and Alexander, 1972;Gosselin et al., 1997). However, the ECOSMO II model represents functional groups instead of species, and the entire ice algae community approach is represented by a single state variable. Following the assumption that diatoms are dominant in the system, the group processes are parametrized for diatom species.
Similar to the pelagic environment, ice algae dynamics are influenced by various abiotic and biotic factors such as temperature, salinity, light, nutrients, and consumption by microorganisms (Arrigo et al., 1993;Werner, 1997;Arrigo, 2014;Galindo et al., 2017;Oziel et al., 2019). The rate of change in the ice algae biomass is calculated as follows: where IA is the sea-ice algae biomass in mgC m −2 . The first term in Eq. 1 represents the gross primary production of ice algae, the second and third term account for the loss of ice algae and the fourth term parametrizes the exchange between ice and water column. The ice algae growth σ depends on the ambient temperature and salinity of brine in the cracks and channels in the ice and can be limited by light and nutrient. We assume a single limiting factor approach rather than multiplicative limitation, such that the specific growth rate is given by: Where σ max is the maximum growth rate ( Table 1). f T and f S are respectively, the temperature and salinity dependence of the sea-ice algae growth, and N lim and L lim account for nutrient and light limitation respectively. Laboratory experiments have shown that algal growth is influenced by these two factors especially when salinities reach high values in the brine (Arrigo and Sullivan, 1992). According to Vancoppenolle and Tedesco (2017) temperature and salinity dependencies are expressed as following: The terms T and S are, respectively, the ambient temperature in • C and salinity, r g is the temperature sensitivity coefficient ( Table 1). The availability of light and nutrients are the two main factors limiting the ice algae growth with different seasonal relevance (during the polar night in winter for light, and during the summer season for nutrients). The ice algae growth is therefore limited by the four nutrients (NO 3 , NH 4 , PO 4 , and SiO 2 ) following a Monod-formulation (Monod, 1949). Even though nutrient stoichiometry inside the ice algae is known to be variable, we assume a fixed Redfield stoichiometry (Vancoppenolle and Tedesco, 2017;Selz et al., 2018).
where N i is the concentration of a specific nutrient (NO 3 , NH 4 , PO 4 , and SiO 2 ) in the sea ice in mmol m −2 , and K s (i) is the halfsaturation constant of the nutrient. Ks values for each nutrient are taken from Sarthou et al. (2005) (see Table 1). Light limitation follows a hyperbolic tangent function as described in previous studies (Lavoie et al., 2005;Tian, 2006;Castellani et al., 2017): where PAR (photosynthetically active radiation) is the light available for ice algal photosynthesis, α p is the photosynthetic efficiency and P m is the maximum photosynthetic rate. Values and units for α and P m ( Table 1) are taken from Castellani et al. (2017). The second term in Eq. 1 correspond to the mortality loss term in mgC m −2 . Loss of ice algae includes losses due to lysis, exudation and respiration and is given as a constant loss rate, rM ( Table 1; Lavoie et al., 2005;Vancoppenolle and Tedesco, 2017).
The third term in Eq. 1 refers to the loss of ice algae by zooplankton grazing. Little is known about the grazing on ice algae, but observational evidence highlights the importance of carbon pathways from ice algae to zooplankton in Arctic ecosystems (Michel et al., 1996). Even though ice algae can be grazed inside the ice matrix (Gradinger et al., 1999), Michel et al. (2002) concludes that the ice grazing constitutes only a minor proportion of the ice algae loss. Therefore we consider grazing by pelagic grazers only following previous modeling approaches assuming a constant uptake rate by pelagic grazers, rG ( Table 1; Lavoie et al., 2005;Tedesco, 2014).
This simplified consideration of zooplankton grazing is justifiable as ice algae concentrations are generally relatively low and, during the ice algae bloom, ice algae constitute almost the only food source for zooplankton in the model and prey selectivity can be neglected. The fourth term in Eq. 1 refers to the exchange at the sea ice/ocean interface and it depends on rate of change of sea ice thickness dH i dt .This term is positive during sea-ice growth and negative during sea ice melt.
where Dia is diatom concentration in mgC m −3 , and ρ ice and ρ water the density of ice and sea water respectively. Even though the sympagic and the pelagic diatom communities feature some similarities they are different in structure. It has been shown that the main source of the ice algae communities stem from the water column, the benthic system and the sea ice itself, but the origin, the amount and the seasonality of the ice algae origin is still unclear . It seems, however, sensible to consider that not all the algae entrapped during sea ice formation are adapted to survive inside the sea ice. Therefore, we advance the assumption that the quantity of algae entering in the sympagic system is controlled by a factor f dia (Table1), corresponding to the fraction of algae which continue to develop after being entrapped in the ice. As pelagic algae are not the only source of ice algae in the sea ice, we fixed a minimum threshold of ice algae in the ice (0.01 mgC m −2 ), to avoid a complete depletion during winter and to decrease the delay of the ice algae bloom at spring. The remaining part does not survive and contributes directly to the sea-ice detritus pool. Because of the lack of observations and knowledge, we arbitrary set this value at 0.5, after a sensitivity analysis (see section "Results: Role of Sympagic System: Implementation, Verification and Sensitivity"). The rate of change in the ice detritus biomass is calculated as follows: where ID is the sea-ice detritus biomass in mgC m −2 . The first term in Eq. 10 represents the loss of ice algae, the second term account for the remineralization process, with a constant remineralization rate ε D ( Table 2). The third term parametrizes the exchange between ice and water column. As for ice algae, the exchange at the sea ice/ocean interface and it depends on rate of change of sea ice thickness dH i dt .This term is positive during sea-ice growth and negative during sea ice melt.

Nutrient Dynamics
Inside sea ice, biological communities exist within a complex network of brine channels, extending into the ice layer, with brine volume varying according to the thermohaline characteristics of the ice. Brine volumes are at maximum in the first few centimeters from the ice-ocean interface (Arrigo, 2014;Thomas, 2017), and this is mainly where biological activity is concentrated (Tedesco, 2014). In this present study, we follow an approach suggested in earlier modeling studies (Arrigo et al., 1993;Duarte et al., 2015). Thus, instead of considering a complex brine network dynamic, we assume a bottom layer with biological activity with about 5 cm thickness. The relative brine volume e in this layer is assumed to be 50% (e = 0.5) as suggested by Duarte et al. (2015).
The ice-ocean interface is a dynamic area where exchanges of salt, organic, and inorganic material occur between the two systems (sympagic and pelagic), determining the salinity and nutrient, detritus, and algae concentrations in the sea ice and the water body below. Exchanges occur due to convection, diffusion, ice formation, and ice melting processes (Vancoppenolle et al., 2010;Arrigo, 2014). In the absence of biological activity, it is customary to consider the dynamics of the majority of macroelements in sea ice following the same dynamical changes as salinity. Meaning that, during the formation of sea ice crystalline lattice impurities and salt are not entrapped in the structure and nutrients and other components are as well rejected in the brines network where they will exchange with the water underneath.
The bulk formulation N i bk for each of the aforementioned types of nutrients i is characterized as follows (Belém, 2002;Vancoppenolle et al., 2010): where N i br is the nutrient content for nutrient i in the brine (mmol m −2 ) and e is the brine volume fraction ( Table 2). Ammonia is an exception to this rule, being entrapped in the crystalline lattice of the ice and expelled to the brine channels (Weeks, 2010;Vancoppenolle et al., 2013). Therefore, the brine concentration of ammonia, unlike the other nutrients, is equal to the bulk mass of ammonia. The rate of change of the respective nutrients in the brine is described as Duarte et al. (2015): ∂N br The first terms in Eqs 13-16 represent the nutrient transport at the interface due to the diffusion mechanism.
where K wi is the molecular diffusion coefficient (Jin et al., 2006; Table 2) and dN i br /dz is the vertical gradient of the respective nutrient at the ice-ocean interface. The second terms G ice_flux in Eqs 13-16 are either the entrapment of dissolved matter in the ice due to the congelation of water process or the release of nutrients during the ice melt. This entrapment/melting acts on the nutrient bulk. It is defined following the relationship defined as follows (Tedesco and Vichi, 2014): where N i w is the water concentration of nutrient i (mmol m −3 ). For ammonium, during the congelation process the quantity is not divided by the brine fraction due to the entrapment into the crystalline lattice as well in the brine. The third terms of the Eqs 13-16 are the uptake of nutrients by ice algae. Here p NO3 in Eqs 13 and 14 differentiate the uptake of nitrate from that of ammonium based on Mortenson et al. (2017): where v n is the half-saturation constant for preferential uptake of nitrate ( Table 2). The fourth terms in Eqs 13 and 14 are the nitrification rate, which is fixed as constant here ( Table 2). The last terms in Eqs 14-16 represent the remineralization process from the detritus pool in sea ice. The remineralization rate is constant ( Table 2).

Biogeochemical Pelagic Modifications
The ocean biogeochemical model ECOSMO includes 13 state variables describing three nutrient cycles (NO 3 , NH 4 , SiO 2 , and PO 4 ), two phytoplankton groups (diatoms, flagellates), two zooplankton groups (mesozooplankton and microzooplankton), oxygen, biogenic opal, detritus (fast and slow), and dissolved organic matter. In the here used setup, the benthic system is included in the model formulation consisting of three state variables, of which each corresponds to a respective nutrient pool in the sediment (nitrogen, silicate, and phosphate). Model equations and parameterization are described in detail elsewhere (Daewel and Schrum, 2013;Daewel et al., 2019). To study the Arctic system and the sympagic-pelagic-benthic interactions, some changes have been made. These changes are explained in the section below. In particular we split the detritus group into a fast (D f ) and a slow (D s ) sinking class. Similar to diatom entrapment into sea ice, during the melting phase not all the ice algae are transferred to the diatom pool but only a fraction, corresponding to the surviving algae. The same fraction f dia as for the entrapment is used. The remaining part is added to the fast sinking detritus pool. Additionally, the sea-ice detritus is released into the fast sinking detritus pool in the water column during ice melt. The reaction equations for two detritus groups are given as: where, T is the temperature in • C and z the depth in m. C X is biomass of the state variable of the X in mgC m −3 , where X represents P j [phytoplankton group with j = 1, 2 denotes the phytoplankton type (1: diatom, 2: flagellate)]; Z i [zooplankton group with i = 1, 2 denotes the zooplankton type (1: micro-, 2: meso-zooplankton] and SE 1 is the sediment group representing the sediment carbon pool. G i−1,2 is the zooplankton grazing rate. A sensitivity analysis (not included in the results) on the different sinking speed for the fast sinking detritus group showed that a sinking speed of 50 m d −1 leads to the best representation of the ecosystem dynamics. This value has also been chosen in earlier modeling studies by Lavoie et al. (2005) and Mortenson et al. (2017). The parameters values for Eqs 18-23 are given in Table 3. Additionally, we added a temperature dependence to the phytoplankton growth rate following an exponential function as described by Eppley (1972): where, T is the ambient temperature in • C and T lim is the reference temperature (Table 3).

Light
Light is one major factor controlling growth of both sea-ice algae and phytoplankton in the water column. Light that penetrates in and through the ice reaching the sea-ice algal communities depends on the snow and ice thickness following the Beer-Lambert law as Castellani et al. (2017): (27) where I 0 is the shortwave incoming solar radiation, α is the albedo of the surface component ice or snow, C 0 is the fraction of incoming radiation absorbed in the first centimeters of snow or ice layer, k ice and k snow are the ice and snow attenuation coefficient respectively and H ice and H snow are the ice and snow thickness, respectively. All parameter values and units are referred in Table 4.
Light available for phytoplankton in the water column below the ice depends on the snow and ice thickness but also on the ice algae concentration. In a given grid cell, however, the available light for primary production in the water column is additionally determined by what fraction of the area is covered by ice A (Slagstad et al., 2015). Light intensity reaching the water column below ice is expressed as: (28) where Z bot is the thickness of the biologically active bottom sea ice layer. The total light available in a grid cell is then given by:

Numerical Experiments and Sensitivity Analysis
With the 1D model setup we performed numerical experiments to understand primary and secondary production in the Barents Sea. Even though the Barents Sea is one of the most productive areas in the Arctic Ocean it features a high spatial variability in ecosystem productivity with less production in the northern than in the southern areas (Sakshaug, 1997;Dalpadado et al., 2014). To address the spatial dynamics, we ran simulations at three separate locations (Figure 2) that differ in ice dynamics and environmental conditions. The southernmost station was located in an area where the ocean is free of ice all year round We performed three sets of scenarios to evaluate the primary and secondary production in the Barents Sea area at different spatial and temporal scale. The first set of experiments was designed to understand the response of the pelagic production system to implementation of the sympagic ecosystem dynamics in the model. These scenarios were run at MBS location during 1997. In this first set of experiments, we ran the same setup (i) including all the biogeochemical components and interaction between the sympagic and the pelagic systems (fullBGC scenario), (ii) neglecting the biogeochemical interaction between the sympagic and the pelagic system but including the sea ice physics (noBGC scenario), and (iii) neglecting the sea ice completely (noICE scenario), which allowed us to dissociate the role of the sea ice presence itself and the role of the sea ice biogeochemistry on the pelagic system. This last experiment also explores the ice free environment, which is expected in the next decades considering climate induced sea ice retreat (Onarheim and Årthun, 2017).
The second set of scenario experiments described the spatial variability of the primary and secondary production at three stations (NBS, MBS, and SBS) following a north-south gradient.
The third set of scenario experiments addressed the interannual variability of the primary and secondary production at the MBS station. We ran each simulation for a 10 years' time period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999), where the first 2 years were considered as spin-up period and were disregarded for analysis.
Since the Arctic ecosystem is difficult to access and data on parameters relevant for the sympagic ecosystem are sparse, we additionally performed a sensitivity study to investigate the models' response to some of the most uncertain parameters in the ice-ecosystem formulation. The proportion of phytoplankton entrapped in the ice during freezing determines the initialization of the sea-ice algae and thus consequently ice algae biomass but is hardly measurable and thus difficult to identify. Some evidence shows that ice algae stem from different sources (multi-year ice, pelagic system, benthic system; Syvertsen, 1991;Ratkova and Wassmann, 2005;Werner et al., 2007;Olsen et al., 2017). But it remains unclear how much exactly originates from the pelagic system and finally survives inside the ice. On the other hand, the proportion of ice algae released into the water during the sea ice melt contributes to the algal biomass in the water column and can be used by zooplankton as a valuable food source (Michel et al., 1996;Søreide et al., 2006Søreide et al., , 2013. Moreover, some studies highlight the potential of some ice algae species to seed phytoplankton production during the ice melt (Haecky et al., 1998;Galindo et al., 2014;Szymanski and Gradinger, 2016). Through a sensitivity analysis, we evaluated the impact of the amount of phytoplankton entrapped in sea ice on ice algae phenology and, of the amount of ice algae released to the pelagic system on the primary and secondary production. We ran four scenarios for this experiment: 1. control run: 100% of diatoms entrapped is transferred into the ice algae pool and 100% of ice algae is release to the diatom pool; 2. 75% of the diatoms entrapped are transferred to the ice algae pool while the remaining 25% are transferred to the ice detritus pool and, 75% of the ice algae released are transferred in the diatoms pool while the remaining 25% are transferred to the pelagic fast sinking detritus pool; 3. where the proportion of algae:detritus is 50:50%; 4. where the proportion of algae:detritus is 25:75%.
The main assumption of the 1D modeling approach described here is that lateral advection in both the ice and water layers is negligible compared to other fluxes. We used forcing data retrieved from outputs of a 3D hydrodynamic model -which includes advection -to define the seasonality of temperature, vertical mixing, and ice properties. However, the 1D model itself does account for net lateral fluxes of state variables at the given location. This limitation does not invalidate the approach as deployed in this article for exploring hypotheses on the linkage between sympagic-pelagic-benthic systems provided the deployment locations are chosen carefully to avoid areas with strong horizontal gradients in the model variables so that advection is of secondary importance. However, it could become an issue for interpretation of responses to climate change. Climate warming is not only accompanied by loss of sea ice, but also by changes in wind patterns, water mass characteristics and circulation (Skagseth et al., 2020). These could change net horizontal fluxes, especially of ice at our study sites in the Barents Sea (Kimura and Wakatsuchi, 2001;Koenigk et al., 2009;Kwok, 2009). Hence advection could become a significant term in the rate of change of ice algae biomass at model locations. Thus, in the next step of our work we will apply the model in a fully coupled 3D biological-physical model to account for the dynamics nature of the system. This will additionally account for sea ice dynamic forced by the atmospheric conditions, which has clear implications for both sympagic as well as pelagic productivity.

Role of Sympagic System: Implementation, Verification, and Sensitivity
In Figure 3 we show the seasonal variations of estimated ice algae dynamics (fullBGC scenario) and the related external and internal properties. The results show that ice algae growth started slowly at the beginning of March as soon as sufficient light became available (Figure 3E), and prevailed until mid-June when biomass started to decrease ( Figure 3C) controlled by both ice melting (Figures 3A,B) and grazing by pelagic zooplankton. The ice algae growth dynamic was characterized by three blooms. The two first occurred in early and mid-May, both followed by a sharp decrease of biomass due to two melting events, generating a strong ice algae release ( Figure 3B). A third bloom, which was the most important, occurred in mid-June (Figure 3C), during which the biomass increased up to 378 mgC m −2 with a maximum daily ice algae production rate about 89 mgC m −2 d −1 . Our biomass estimates were consistent with data from field measurement in the Barents Sea area (69-620 mgC m −2 ; Tamelander et al., 2009) and in the Arctic region (3-460 mgC m −2 ; Gosselin et al., 1997). The results on daily ice algae production rate, however, appear to be higher than measured production rates in the northern Barents sea (0.16-55 mgC m −2 d −1 ; Hegseth, 1998;McMinn and Hegseth, 2007). On the other hand, our results showed low biomass of ice algae when compared with data from Rózañska et al. (2009) and Fernández-Méndez et al. (2018), who reported ice algae biomasses up to 2250 mgC m −2 , indicating that our model might underestimate ice algae biomass. The abiotic environment (e.g., light, ice, and snow thickness) strongly influences the daily ice algae production rate. These environmental conditions vary strongly in the Barents Sea both spatially and temporally, which can explain these differences between model results and observations, particularly since our environmental forcing is only representative for one specific location. Moreover, we only consider ice algae originating from entrapment of pelagic algae, which neglecting ice algae advection from other area. As this provide a relatively small quantity of algae it impacts the growth of biomass and might explain the slow growth between March and May, and the potentially too low ice algae biomass. Nonetheless, over the large existing range of measured and simulated production rates (from 0.16 to 100 mgC m −2 d −1 ; Clasby et al., 1973;Mortenson et al., 2017) our simulated values are well in the range of previously reported production rates and our daily production rates appear acceptable.
The simulated ice algae bloom dynamic was not nutrient limited in our scenario (Figure 3E), even though the silicate concentration decreased sharply during the bloom event, it never became limiting (Figures 3D,E). Our simulations show two strong melt events during the ice algae bloom period which removed a substantial amount of ice algae biomass from the ice, which also reduced the nutrient requirements in the ice. In addition, the ice bottom layer at the ocean interface is highly dynamic and nutrients can be replenished by convective mixing (Vancoppenolle et al., 2010). This is in compliance with findings from Werner et al. (2007), who showed that ice algae do not appear to be limited by nutrients in the Barents Sea. The Barents Sea is a shallow area and is significantly influenced by winds and tidal mixing (Sundfjord et al., 2007). Hence nutrients are abundant and not generally a limiting factor for ice algae as also shown by Tamelander et al. (2009). The total annual production of the simulated ice algae was about 2.44 gC m −2 , which represented about 3.8% of the total primary production. This fraction is in the lower part of the range of the observations previously reported for the whole Arctic (3-57%; Gosselin et al., 1997). According to Sakshaug (1997) and Hegseth (1998) who reported estimated total annual production respectively around 6 and 7% of the total primary production in the northern Barents Sea, our total production estimation appeared to be slightly underestimated but reasonable.
In order to assess the role of the sympagic system on the overall ecosystem dynamics, we performed three numerical experiments to separate the impacts of the physical and biological forcing. The results are presented in Figure 4 for the fullBGC (Figures 4A,D,G,J), the noBGC (Figures 4B,E,H,K) and the noICE scenario (Figures 4C,F,I,L). Scenario explanations are given the methods section. The results highlight how the sympagic system generates important changes in the whole ecosystem by modifying the phytoplankton and zooplankton phenology.
What we generally can see is that implementing ice biogeochemistry in the model formulation strongly impacts the dynamics of the pelagic ecosystem, both with respect to magnitude and timing of the different components. Most importantly, we found major differences in the dynamics of diatoms and zooplankton. In the fullBGC scenario, when the sympagic system was fully implemented, the model estimated two diatom blooms during the season (Figure 4A), while only one spring bloom was simulated in noBGC. The first bloom occurred in July in both scenarios, partly under ice, when the ice coverage started to decrease, allowing better light conditions in the water column. In case of fullBGC, this bloom was not limited by nutrients but instead controlled by zooplankton grazing (Figures 4G,J). In the noBGC scenario, on the other hand, this first diatom bloom developed until silicate became limiting. The major process leading to this difference was that generally the zooplankton benefits from feeding on ice algae prior to the first diatom bloom, leading to relatively high zooplankton biomasses early in the season ( Figure 4D). Consequently, zooplankton and phytoplankton maximum biomass occurred at the same time (Figures 4A,D), which systematically modified the overall system dynamics. Therefore, silicate remained in the water column allowing a second diatom bloom in August during full open water condition following the ice break-up. Zooplankton grazing on ice algae, in the ice and in the water column just after being released, appears one of the most important links between the pelagic and sympagic component of the ecosystem leading in our results to a general increase in zooplankton biomass. This important link has earlier been reported on the basis of field observations (Stretch et al., 1988;Runge and Ingram, 1991;Werner, 1997;Juul-Pedersen et al., 2008). Observations indicate that the grazing on ice algae is partly used by female copepods to mature their gonads and produce eggs (Søreide et al., 2010;Durbin and Casas, 2014). In this case, there is an increase in biomass due to egg production. But because eggs do not feed directly on phytoplankton during their development and only participate in phytoplankton consumption later in the season after hatching, the control of our first phytoplankton bloom by the zooplankton grazing could potentially be overestimated. More investigations are necessary here to accurately understand the direct role of ice algae on zooplankton dynamics and vice versa. Nonetheless, the interaction between the pelagic zooplankton and ice algae remains an important process as it determines the fitness of these zooplankton species and includes other species that use ice algae such as amphipods and krill.
The flagellate bloom occurred after the ice break-up at the same time as the second diatoms bloom ( Figure 4A) but with a much higher biomass. Flagellates generally dominated the primary production with a total annual production about 42.73 gC m −2 , which represent around two thirds of the total pelagic primary production (61.61 gC m −2 ). Previous studies estimated annual primary production for the Barents Sea to be between 20 and 200 gC m −2 , depending on the area, with estimates around 76 gC m −2 in the central Barents Sea (Luchetta et al., 2000;Sakshaug et al., 2009), which supports our model estimates. The high flagellate biomass results from the fact that the diatoms are mainly controlled by zooplanktonic grazing. This leads to relatively high nutrient concentrations ( Figure 4J) and consequently allows a strong flagellate bloom in the fall. Also, for mesozooplankton a second maximum follows the spring maximum in July and is closely related to the diatom dynamics. Microzooplankton, on the other hand, followed the flagellate's dynamic with a unique peak by the end of summer. Mesozooplankton dominated the secondary production.
In Figure 5, we showed differences between the fullBGC and the noBGC scenario in the vertical biomass profiles of the two detritus groups (Figure 5A: fast and Figure 5B: slow sinking rate). These results indicate a lower biomass of slow sinking detritus during summer (Figure 5B) when the sympagic system is coupled to the pelagic ecosystem, related to a smaller summer diatom bloom in the fullBGC scenario. This was followed by a second strong detritus production in autumn, which was exported to the benthic system in late autumn, and which corresponds to the strong flagellate production. Direct export of ice algae to the benthic system happened at the end of summer when ice melts. We observed a contribution of the fast sinking detritus group in September with a slightly higher export to the benthic system in the fullBGC scenario compared to the noBGC scenario ( Figure 5A). It has already been shown that ice algae contributes to the export of matter to the benthos due to the fast sinking rate and it represents a valuable resource for benthic organisms early in the season (Søreide et al., 2013;Schollmeier et al., 2018). Due to the forcing data we used, the simulated ice break-up occurred late in the season at mid-June and we did not observe a spring contribution of ice algae to the benthos as was showed in previous studies (Boetius et al., 2013;Szymanski and Gradinger, 2016). However, our results pointed out that the ice algae are exported rapidly to the benthos and earlier than the phytoplankton in the productive season (Figure 5). In addition, our results highlight the indirect contribution to detritus export by changes in zooplankton dynamics ( Figure 5B). Because of the feeding process on ice algae, zooplankton developed earlier and participated in an increase in the slow sinking detritus biomass.
The results indicate temporal changes on the export of detritus biomass to the benthic system when sympagic biogeochemistry is included. The sympagic system therefore does not only play a role in pelagic dynamics but also influences both the amplitude and timing of detritus export to the benthic system. Here, we additionally hypothesize that the indirect contribution due to FIGURE 6 | Simulated time series of (A) diatoms/ice algae flux at the ice-ocean interface; (B) production rates of ice algae; (C) production rates of phytoplankton including diatoms (D: solid line) and flagellates (F: dashed line), and (D) production rates of zooplankton including mesozooplankton (Meso: solid line) and microzooplankton (Micro: dashed line). Black color presents result for the 100% experiment, red color for the 75% experiment, green color for the 50% experiment, and purple color the 25% experiment. In the legend (B-D) are indicated the annual production in gC m −2 for ice algae, phytoplankton including diatoms production (D) and flagellates production (F) and zooplankton including microzooplankton (Micro) and mesozooplankton (Meso). the zooplankton feeding on the ice algae is as important for the timing and magnitude of material export to the benthic system.
The noICE scenario shows similar dynamics to the noBGC scenario. The main difference observed is an earlier onset of spring bloom production with about one-month time difference (Figures 4B,C). This is due to the absence of ice cover leading to sufficient light conditions earlier in the year (Figures 4K,L). The other major difference is the increase of the total annual primary and secondary productivity in the noICE scenario. The structure of the ecosystem groups does not change between these two scenarios suggesting that the presence of ice only as a physical barrier does not seem to influence pelagic ecosystem functioning. Differences between fullBGC and noICE scenarios, i.e., earlier primary production and higher primary production, reflect changes expected for the Arctic region regarding an ice decrease. These changes are further discussed in the section Spatial variability in the Barents Sea, environmental drivers and ecosystem dynamics.
Results of the sensitivity analysis (Figures 6A,B) showed that the proportion of phytoplankton entrapped in the ice influenced slightly the total ice algae production with a change of about 0.25 gC m −2 y −1 between the 100 and 25% numerical experiments, which represent a decrease of about 10% for an entrapment proportion decrease of 75%. The sensitivity experiments did not indicate any changes in ice algae production in terms of timing and amplitude of the maximum growth rate. This principally resulted from the minimum threshold implemented in our model to allow a minimum proportion of ice algae to survive during winter in the ice when conditions are not favorable for growth. This assumption was made because the fate and origin of ice algae during winter remains poorly understood. Different hypotheses were suggested previously including the existence of resting cysts or spores, metabolism adaptation to low light, facultative heterotrophy, or origin of algae in the ice from nearby multi-year ice (Syvertsen, 1991;Niemi et al., 2011;Kauko et al., 2018). Even though the total amount of ice primary production changed only slightly, we found that the magnitude of the initial bloom in May was decreased by almost 50% when only 25% of the entrapped algae were used for production.
Although variations in phytoplankton entrapment did not show a prominent effect on the sympagic production, the modeled pelagic system showed a clear response to the amount of ice algae released into the sympagic system during ice melt (Figures 6C,D). First, our results indicate that released ice algae seeded the diatoms blooms, both under ice and during the ice break-up. This is one of the main assumptions made in our approach. The hypotheses about the ice algae seeding are divergent and controversially discussed in the scientific community, because evidence is not clear. Some field samplings highlight that ice algae contribute to the start of the phytoplankton bloom, due to some species overlapping in both the sympagic and pelagic communities (Haecky et al., 1998;Lizotte, 2001;Mangoni et al., 2009). Contradictory field studies have indicated that ice and pelagic algae communities are distinct and that ice algae apparently do not seed phytoplankton blooms (Riebesell et al., 1991;Mundy et al., 2014). Our results showed that, under the condition of ice algae seeding, the proportion of ice algae released into the diatom pool has consequences for the phenology of the diatoms, with an earlier peak biomass the more ice algae are released. This consequently also results in a modified phenology of the secondary production. The maximum difference in peak bloom timing for the two most extreme experiments (100 and 25%) was about 1 week, both for the under-ice bloom and the late summer bloom. The results thereby showed that the peak production occurred later the less important the ice algae seeding is for the overall diatom biomass. This estimated change in diatom phenology systematically modified the primary producers' dynamics. The associated environmental conditions, such as temperature, ice cover, and available light, in the water column change slightly within this 1-week period. This influences the growth rate of the diatoms such that, although the seeding of diatoms was lower, it led to a shorter bloom with higher peak biomass. After the ice break-up flagellates developed quickly, limiting the diatom bloom. Because ice algae seed the diatom blooms it allowed them to bloom earlier at the ice break-up. Therefore, the smaller the ice algae release is, the weaker the second bloom and the greater the flagellate's production. The difference in total flagellate production between the 100 and 25% experiments was an increase of 40.8%. In light of our results and also model results published by Tedesco et al. (2012) and Jin et al. (2007) ice algae seeding appears to be a key process in the Arctic marine ecosystem dynamics and further investigations should be done in order to decrease uncertainties about this process. Especially as Lizotte (2001) highlighted, climate change can alter and modify algae communities and consequently potentially change the contribution of ice algae to the phytoplankton bloom.
The sensitivity analysis again highlights the importance of clarifying the current parameterization of sympagic components and the exchanges with the pelagic components underneath. Despite ice algae representing a small proportion of the total primary production, we showed that its contribution to the phenology of all the pelagic groups is significant. The results also suggest important changes in the Arctic pelagic structure of planktonic communities might be related to ice algae changes, as already suggested before for copepods in an earlier study (Tedesco et al., 2019). Indeed, there is evidences that females of C. finmarchicus and C. glacialis, when barely out of diapause, benefit from feeding on ice algae to maturate their gonads and produce eggs (Hirche and Kosobokova, 2003;Søreide et al., 2010;Leu et al., 2011). Although, the factors of entry and exit in a diapause state remain not well identified, some evidence about a mismatch between the presence of individuals at the surface and algal bloom have already been highlighted, showing significant consequences for fitness and next generations of these species (Dezutter et al., 2019). Diapause behavior is not yet a part of our model formulation and should be integrated in future work to identify more accurately the planktonic phenology in the Arctic system as modeling work of Pourchez (2018) showed for the Beaufort Sea. We previously discussed a possible model artifact of zooplankton controlling the first phytoplankton bloom due to an earlier growing phase feeding on ice algae. If we consider that a significant zooplankton biomass migrates from the bottom to the surface at springtime, it seems realistic that such a grazing control takes place. It could also help to explain the second phytoplankton bloom in autumn, when naturally zooplankton migrates back to the bottom and relieves the grazing pressure (Pourchez, 2018). Identification of these key parameters becomes all the more urgent in the context of rapid climate change in order to properly understand and manage the near future Arctic ecosystem.

Spatial Variability in the Barents Sea, Environmental Drivers and Ecosystem Dynamic
The Barents Sea area exhibits strong seasonal contrasts as well as high spatial variability, as illustrated in Figure 7 regarding the FIGURE 9 | Simulated time series of production rates of zooplankton including mesozooplankton (solid line) and microzooplankton (dashed line). Black color presents result for the central Barents Sea station over the 1997, red color for the north Barents Sea station over the 1997, and green color for the south Barents Sea station over the 1997. In the legend are indicated the annual production in gC m −2 for zooplankton including microzooplankton (Micro) and mesozooplankton (Meso).
light availability at the ocean surface ( Figure 7A) and the sea ice thickness (Figure 7B). This high variability in environmental conditions results in very different ecosystem dynamics and follows mainly a south-north gradient related to temperature and ice cover (Figures 8, 9). For the northernmost station (NBS), which is covered by ice all year round, our results showed the lowest productivity, with a total annual primary production of 2.86 gC m −2 yr −1, while at the central and southernmost location our model estimated 64.05 gC m −2 yr −1 and 96.89 gC m −2 yr −1 , respectively (Figure 8). This was due to the presence of a thick ice coverage in the north leading to low light availability at the iceocean interface. Our results are consistent with field observations provided by Sakshaug et al. (2009), who estimated that primary TABLE 5 | Annual production of ice algae, phytoplankton (diatom + flagellate), and zooplankton (microzooplankton + mesozooplankton) for the central Barents Sea station over the 7 years simulated (1992)(1993)(1994)(1995)(1996)(1997)(1998).

Year
Annual ice algae production (gC m −2 y −1 ) Annual phytoplankton production (gC m −2 y −1 ) Annual zooplankton production (gC m −2 y −1 ) production varied from <20 to 200 gC m −2 yr −1 over the Barents Sea. Even though ice algae production was relatively small in the north when compared to the central Barents Sea area, it was more important for the overall productivity in the region representing around 42% of the total primary production at this location. In contrast, in the central Barents Sea the ice algae accounted only for 3.8% of the overall primary production ( Figure 8A). Light conditions not only determined the magnitude but also the phenology of the productivity resulting in an onset of spring diatom bloom in the beginning of June in the southern Barents Sea, beginning of July in the central Barents Sea, and middle of August in the northern station ( Figure 8B). A second diatom bloom developed only in the central Barents Sea area in the middle of August when ice break-up allowed more favorable light conditions. Flagellates bloomed twice (end of June and beginning of August) in the south, whereas in the central Barents Sea only one bloom occurred in the middle of August. Environmental conditions at the northern location are generally not favorable for primary production and flagellates were outcompeted by diatoms. In contrast, flagellates dominated the primary production in the south and central part of the Barents Sea by around two thirds. For mesozooplankton and microzooplankton the seasonal dynamics followed closely phytoplankton blooms (Figure 9). Thereby, in the south zooplankton were mainly concentrated at the beginning of July following the diatoms and first flagellate blooms and slightly developed at the fall. In the north, only very little mesozooplankton developed in August following the diatom bloom. In the central area the seasonal dynamics were very different featuring two major zooplankton peaks in July and September, respectively. Contrary to the phytoplankton phenologies our results indicated that no phenological shift occurred in secondary production between the southern and the central Barents Sea location. The reason was that the zooplankton also feeds on ice algae allowing an increase in zooplankton biomass even before the pelagic diatom bloom starts and therefore promoting an early peak in zooplankton biomass.
Another interesting emergent result is that although primary production in the south was around 30% higher than in the central Barents Sea, the zooplankton annual total production was not significantly different between these two regions suggesting a more efficient trophic transfer from phytoplankton to zooplankton, again caused by the feeding process of zooplankton on ice algae.
Going from north to south in the Barents Sea may also be indicative of the Barents Sea ecosystem through time. Considering the expected sea ice retreat in the light of climate change during the next decades in the Barents Sea, we can expect a change toward a more southern ecosystem in terms of structure and functioning. This means an increase in primary production for the central Barents Sea and an increase of both primary and secondary productions for the northernmost regions. Previous modeling studies have already highlighted the potential increase of primary production in the Arctic and in the Barents Sea due to climate change and sea ice reduction and support our study results (Arrigo et al., 2008;Holt et al., 2016). Because the Arctic ecosystem structure has a relatively low complexity with few key players (Darnis et al., 2012), changes in the sea ice extent and ecosystem can greatly affect the whole Arctic food-web structure and functioning. As an example, according to the modeling study of Suprenand et al. (2018) important changes in zooplankton have already happened in the Beaufort Sea, and suggest greater changes in the marine ecosystem will occur in the next decades due to the sea ice decrease and temperature increase presently ongoing. Several species directly or indirectly depend on the ice algae, especially female Calanus for the eggs production as already discussed before.

Interannual Variability Pattern of Primary and Secondary Production
An additionally performed multi-year simulation experiment in the central Barents Sea highlights the interannual variability in this area (Sakshaug et al., 2009) with respect to magnitude and seasonality in primary and secondary production ( Table 5). Within the simulated 7 year time period total annual primary production varied in the range of 34.33-95.92 gC m −2 y −1 , with associated response in secondary production, which varies in the range of 13.61-49.76 gC m −2 y −1 . Over these years the ice algae production represented always less than 10% of the total annual primary production also with strong interannual variation. The results indicated that years with long ice seasons, but thin ice thickness were the most productive for the ice algae (e.g., 1998). Our results did not demonstrate a clear and recurring seasonal pattern over the years due to the high variability of the ice cover dynamic. We can highlight the existence of a relationship between the length of the ice season and productivity of the ecosystem, which was already described by Rysgaard et al. (1999), showing a clear linear relationship between the annual pelagic primary production and length of productive open water period. Figure 10 shows clearly how the timing of the ice break-up drives the magnitude of the primary production. It becomes evident that the earlier in the season the ice breaks up, the higher the total annual primary production, which is consistent with findings by Arrigo and van Dijken (2015). This result also implies an increase of the total primary production in the next decade with the decrease of sea ice coverage and thickness and with a shorter period of ice coverage. This trend is already ongoing and observed in the Barents Sea and the other Arctic Seas (Frey et al., 2018) and is expected to be extended northward in future. Our results however suggested that the maximum productivity did not occur in an ice-free year, but rather during a year showing thin ice coverage and an early ice break-up, as it happens in 1995. This finding indicates that primary production in the Barents Sea may not necessarily increase under upcoming icefree conditions but could even decrease especially in areas with seasonal sea ice dynamics. Additional numerical experiments in a 3D context will help to explore this pattern further with respect to regional scales. Previous modeling work from Holt et al. (2016) agreed with the heterogeneous effects of the climate change and ice retreat on the Barents Sea productivity and showed for instance that, in contrast to the northern and north-eastern Barents Sea, in the coastal southern Barents Sea production is expected to decrease.

CONCLUSION
The overarching objective of our modeling approach was to understand the linkage between sympagic and pelagic ecosystem components by implementing a sympagic component to the biogeochemical model ECOSMO. We used the model in 1D at different locations in the Barents Sea ecosystem and analyzed spatial and interannual differences in ecosystem dynamics. The sympagic model component was based on previous developments (Lavoie et al., 2005;Tedesco and Vichi, 2014;Castellani et al., 2017;Vancoppenolle and Tedesco, 2017) and meant to improve the ECOSMO model performance in the Arctic ecosystem. In general, our results showed that the simulated seasonality of the sympagic components were realistic and in line with the current knowledge of the Barents Sea ecosystem. Compared to earlier model applications, our study specifically highlighted the strong interactions between the sympagic and the pelagic ecosystems. Most importantly, we could show that the ice algae play an important role for zooplankton phenology as it allows zooplankton biomass growth prior to the onset of the pelagic spring bloom. This process consequently determines the whole ecosystem structure as the zooplankton impacts the amplitude and the timing of the pelagic primary and secondary production. Additional relevant processes are the timing of the ice break-up and that ice algae production determine the phenologies of phytoplankton through seeding. Therefore, we propose that a fully coupled pelagic-sympagic ecosystem is a prerequisite to simulate the Arctic Ocean ecosystem and its changes correctly.
Exchange processes at the ice-ocean interface remain still poorly understood and our sensitivity analyses showed that the parametrization of these processes is a key component for our understanding of the Arctic ecosystem dynamics. With our approach we also found that the proportion of ice algae released and surviving in the pelagic system had an effect on the phytoplankton and zooplankton dynamics. Consequently, slight changes in the ice dynamic can disturb the whole structure of the present Arctic system and can modify phenologies not only of phytoplankton but also all the others species closely related of the ice algae, such as zooplankton (Werner, 1997;Søreide et al., 2010;Kohlbach et al., 2017a), fish (Kohlbach et al., 2017b;Steiner et al., 2019), or bird (Ramírez et al., 2017;Cusset et al., 2019). As an example, Polar cod (Boreogadus saida) is a central key species of the Arctic ecosystem and it is closely related to the ice algae especially during their larval development feeding on amphipods, which directly grazed on ice algae and sometimes directly feeding on ice algae itself (Gilbert et al., 1992;Kohlbach et al., 2017b). A decrease in the ice coverage will consequently affect the life cycle of Polar cod, by either creating a mismatch situation with the ice algae or changing the diet of larval Polar cod. We should expect changes in the Polar cod distribution in the next decades, with migration northward, and potentially an extension of the capelin population to replace the Polar cod as already suggested by Hop and Gjøsaeter (2013) in the Barents Sea and Steiner et al. (2019) in the Western Canadian Arctic. Moreover, we can expect important change in the abundance and distribution of the predators feeding on both Polar cod and capelin. This example demonstrates, that further numerical experiments including higher trophic level of the marine ecosystem and coupling in a 3D framework will allow us to further explore these hypotheses.
Due to the harsh environmental conditions the Arctic remains a poorly sampled ecosystem leaving many uncertainties in the understanding of its ecosystem functioning. In the context of rapid climate changes, collaborative work becomes even more urgent to fill the knowledge gaps and improve our understanding of the Arctic system in order to establish sustainable management plan in the next decades.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
DB performed and analyzed the biogeochemical model simulations. DB and UD wrote the main manuscript. All authors discussed the results and implications and participated in the manuscript improvement at all stages.

FUNDING
This work was jointly funded by the UK Research and Innovation Natural Environment Research Council (NERC) and the German Federal Ministry of Education and Research (BMBF).