Estimates of carbon sequestration potential in an expanding Arctic fjord (Hornsund, Svalbard) affected by dark plumes of glacial meltwater

In polar regions, glaciers are retreating onto land, gradually widening ice-free coastal waters, which are known to act as new sinks of atmospheric carbon. However, the increasing delivery of inorganic suspended particulate matter (iSPM) with meltwater might significantly impact their capacity to contribute to carbon sequestration. Here, we present an analysis of satellite, meteorological, and SPM data as well as results of a coupled physical–biogeochemical model (1D GOTM-ECOSMO-E2E-Polar) with a newly implemented iSPM group to show the impact of iSPM on the ecosystem dynamics in a warming polar fjord (Hornsund, European Arctic) with numerous shallow-grounded marine-terminating glaciers. Our results indicate that with a longer melt season (9 d per decade, 1979–2022), the loss of sea ice cover (44 d per decade, 1982–2021) and the formation of new marine habitats after the retreat of marine-terminating glaciers (around 100 km² in 1976–2022, a 38 % increase in the total area), glacial meltwater has transported increasing loads of iSPM from land (3.7 g m⁻³ per decade, reconstructed for 1979–2022). The simulated light limitation induced by the iSPM input delayed and decreased the peaks in phytoplankton, zooplankton, and macrobenthos. The newly ice-free areas still markedly contributed to plankton primary and secondary production and carbon burial in sediments (5.1, 2.0, and 0.9 Gg C yr⁻¹, respectively, on average for 2005–2009 in the iSPM scenario). However, these values would have been 5.0, 2.1, and 0.1 Gg C yr⁻¹ higher, respectively, without the iSPM input. Since carbon burial was the least affected by iSPM (a decrease of around 16 %, in comparison to 50 % for plankton primary and secondary production), the impact of marine ice loss and enhanced land–ocean connectivity should be investigated further in the context of carbon fluxes in expanding polar fjords.

1 Introduction

Organic carbon burial in marine sediments represents the dominant natural pathway toward long-term sequestration and hence plays a key role in controlling atmospheric O₂ and CO₂ concentrations (Berner, 1982; Hedges and Keil, 1995). While important carbon sinks at coastal wetlands (mangrove forests, salt marshes, and seagrass beds) are declining globally (Duarte et al., 2005; Howard et al., 2014), new marine habitats are opening up in the Arctic and West Antarctic due to retreating glaciers and giant iceberg calving (Ficetola et al., 2021). Within these coastal ecosystems, CO₂ drawdown by phytoplankton and ice algae is supported by nutrient input from land, intensifying the cascade from carbon capture into storage and burial in sediments (Ardyna and Arrigo, 2020; Arrigo et al., 2008; Wadham et al., 2019). Due to the high sedimentation rates, emerging and expanding fjords play an important role as efficient carbon burial hot spots (Bianchi et al., 2020; Cui et al., 2022; Smith et al., 2015). Thus, the loss of marine ice in polar coastal waters might to some extent compensate for decreasing coastal carbon sinks elsewhere (Barnes, 2017; Peck et al., 2010; Zwerschke et al., 2022).

Despite recent increases in primary and secondary production due to, among other factors, the earlier break-up of seasonal sea ice, the polar regions' potential for long-term carbon burial in sediments is ultimately limited by multifarious mechanisms. Changes in the duration and composition of ice algae blooms weaken the sympagic–benthic coupling, consequently leaving more biomass that can be utilized and dispersed in the pelagic system (Fadeev et al., 2021; Lalande et al., 2019; Riser et al., 2008). Thus, warming induces the maturation of polar fjords, i.e. the transition to a more complex but effective pelagic food web consuming most of the available organic matter, and less carbon is deposited at the bottom (Węsławski et al., 2017; Zaborska et al., 2018). Furthermore, the delivery of inorganic suspended particulate matter (iSPM) with glacial meltwater dims the underwater light later in the productive season (summer and autumn) (Szeligowska et al., 2022) and results in a significant reduction in phytoplankton and phytobenthic biomass (Blain et al., 2021; Deregibus et al., 2016; Holt et al., 2016) that influences the carbon burial potential in glacial bays.

Further warming will likely exacerbate sediment inputs through increasing precipitation, storm activity (enhancing erosion), glacial melt, permafrost thaw, and sea-level rise (Syvitski et al., 2005, 2022). Moreover, in situ observations and numerical simulations from Arctic fjords suggest that after marine-terminating glaciers retreat onto land, subglacial discharge and nutrient upwelling cease, enhancing surface stratification, weakening vertical mixing, and thereby reducing the productivity in coastal zones (Hopwood et al., 2018; Meire et al., 2017). However, our understanding of the rapid transformations of polar marine ecosystems under climatic stressors remains insufficient due to, among other reasons, the scarcity of long-term standardized monitoring data (Schofield et al., 2010). While numerical models have been essential for filling knowledge gaps related to the mechanisms of nutrient supply with meltwater (Castelao et al., 2019; Oliver et al., 2020) and substantial effort has been put into incorporating modules representing biogeochemistry in sea ice (Steiner et al., 2016), only a few of them resolve inorganic particulate matter dynamics in glacier-fed basins (Neder et al., 2022). So far, these models do not typically represent the impact of the delivery of terrigenous material on biological production and carbon budgets.

This study aimed to assess the gains and losses in plankton primary and secondary production and carbon burial due to the transformations of European Arctic coastal waters. We investigated Hornsund (Svalbard, West Spitsbergen) as a model high-latitude fjord since it is among the best-studied fjords in the Arctic and represents an area of rapid regional warming, with many bays affected by the recession of glaciers. Thus, here we (1) map the extent of emerging habitat after the retreat of marine-terminating glaciers and (2) simulate how the ecosystem dynamics and carbon sequestration are affected by sediment discharge in these bays using a 1D coupled physical–biogeochemical model (GOTM-ECOSMO-E2E-Polar) with a newly implemented iSPM group. We present the results of our simulations for 2005–2009, i.e. a period with an exceptionally strong warming signal (Muckenhuber et al., 2016; Promińska et al., 2017), in the context of multidecadal (1976–2022) changes in the physical environment to discuss the potential of newly ice-free areas to act as emerging carbon sinks and their role in the global carbon cycle.

2 Methods

2.1 Study area

Hornsund is a glaciomarine fjord in Svalbard, with inner basins affected by glacial outflow (Fig. 1) (Błaszczyk et al., 2019). Since the strong polar front formed by the West Spitsbergen Current (saline and warm Atlantic Water) and the Sørkapp Current (cold and relatively fresh Arctic Water) reduces the advection of Atlantic Water into Hornsund in comparison to other West Spitsbergen fjords (Promińska et al., 2017), it is considered a less mature, highly productive cold-water fjord with an Arctic-type resident biota and relatively high sequestration of organic carbon (Węsławski et al., 2017; Zaborska et al., 2018). Characterized by dynamic paraglacial coastal systems with high sediment mobility, Brepollen is the most extensive bay in Hornsund, where > 85 km of new shoreline was formed in the last century after an ice retreat (Strzelecki et al., 2020). The area is known for exhibiting one of the fastest retreat rates of marine-terminating glaciers in the Svalbard archipelago; it has accelerated this century, reaching around 3 km² yr⁻¹ in 2001–2010 (Błaszczyk et al., 2013). Importantly, the ice bridge between Brepollen and Hambergbukta (Fig. 1b; currently < 5 km wide) is predicted to break up in the coming decades (Grabiec et al., 2018; Osika et al., 2022), thus reopening a direct connection to the Barents Sea and changing the hydrodynamic conditions for biological production and carbon burial by either stronger sea ice or Atlantic Water advection.

Figure 1 (a) The Svalbard archipelago, showing the monthly mean sea ice concentration > 15 % in March 2005 (plain colour) and 2006 (dotted). Land and glacier extents were downloaded from https://geodata.npolar.no/, last access: 10 September 2023. The red frame indicates the location of Hornsund. (b) Newly ice-free areas in Hornsund which have opened since 1976 (blue lines – the positions of the glacier front in 1976, 2006, 2010, and 2022) and the width of the ice bridge between Brepollen and Hambergbukta. The dots indicate modelled stations (1–20, red) and three data points each for the sea and ice surface temperature (SST, pink) and sea ice concentration (SIC; blue). The red frame indicates the location of Hansbukta. (c) Long-term SPM and sediment flux monitoring stations in Hansbukta (M1–M5, red dots) and Polish Polar Station Hornsund (PPS). The Landsat 8 satellite image (4 April 2020) was downloaded from https://glovis.usgs.gov/app, last access: 27 January 2023.

2.2 Datasets

2.2.1 The area and volume of newly ice-free marine habitats

Summertime Landsat images of Hornsund were downloaded from https://glovis.usgs.gov/app, last access: 3 October 2022 (Table S1 in the Supplement). Only cloud-free images with no sea ice cover (from July to early September) were used. When present, the 4–3–2 and 3–2–1 spectral bands (Landsat 8 and Landsat 1–7, respectively) were used to prepare RGB composites, and a panchromatic band (8) was used to enhance the resolution. Newly ice-free areas were manually delineated, with the position of the glacier front in 1976 used as a reference since this was the first year for which Landsat images taken during the summer were available (Fig. 1b). The same person (MS) repeated the procedure three times for each year to test the repeatability of manual digitization. The standard deviation was up to 0.23 km². Importantly, the fronts of marine-terminating glaciers undergo seasonal fluctuations, which might increase the uncertainty (Błaszczyk et al., 2021, 2023). However, here we narrowed the analysis to the main melt season (from July to early September). Marine habitat volume was calculated based on digitized area and bathymetry data from Hornsund (grid size 100 m) (Moskalik et al., 2014) using the zonal statistics method in ArcGIS Pro 2.8.0. Marine habitat volume was calculated until 2010, since the bathymetric data were not available for the glacial bays that emerged after 2010.

2.2.2 Sea and ice surface temperature and sea ice concentration

Sea and ice surface temperature (SST) and sea ice concentration (SIC) were extracted from the Arctic Sea and Ice Surface Temperature dataset (L4, 5 km, daily). These data were provided by the Danish Meteorological Institute and the MyOcean regional data assembly centre and were created using multisensor satellite surface temperature observations. Since the dataset did not cover all the fjord, data were extracted from points in the outer/central parts (Fig. 1, Table S2), assuming that they reflect the SST/SIC conditions in the inner fjord. This assumption is supported by previous studies (Arntsen et al., 2019; Błaszczyk et al., 2021; Sutherland et al., 2013) and the fact that this analysis focuses on relative changes in melt season intensity rather than absolute values. The data were extracted for three adjacent cells (Table S2) and averaged. Sea-ice-free (SIF) days was defined as the fraction of the year with SIC < 15 %. The monthly mean extents of SIC > 15 % in March 2005 and 2006 are shown in Fig. 1a. The sum of all daily SST values > 0 °C (the positive degree day, PDD SST) was calculated for each year (annual) and each melt season (summertime, June–August) as a proxy for submarine melt potential (Hock, 2005; Rignot et al., 2008).

2.2.3 Air temperature and precipitation

Air temperature and precipitation datasets from Polish Polar Station Hornsund (PPS, Fig. 1c) were downloaded from Wawrzyniak and Osuch (2020) (1979–2018, https://doi.org/10.1594/PANGAEA.909042) and from SIOS – Svalbard Integrated Arctic Earth Observing System (Wawrzyniak and Osuch, 2020, 2018–2022). The daily average air temperature (AT) was used to calculate the sum of all daily AT values > 0 °C (PDD AT) for each year (annual) and each melt season (summertime, June–August) as a proxy for surface melt potential (Hock, 2005; Rignot et al., 2008). Annual and summertime (June–August) precipitation were calculated by summing the daily precipitation measurements (mm). The start of the melt season was defined as the start of the first period of 6 consecutive days with AT > 0 °C; similarly, the end of the melt season was defined as the first of 6 consecutive days with AT < 0 °C (modified from Błaszczyk et al., 2021). This takes into account the delays in meltwater and particulate matter delivery to the fjord, and the 6 d window was shown to be well correlated with sediment flux (in this study and D'Angelo et al., 2018). Melt season duration was calculated as the number of days between the end and the start of the melt season and provided as a fraction of the year.

2.2.4 Suspended particulate matter, sediment flux, and salinity

Datasets for suspended particulate matter (SPM), sediment flux, and salinity collected in Hansbukta (2015–2021) at long-term monitoring stations (Fig. 1c, Table S2) were downloaded from https://dataportal.igf.edu.pl/group/longhorn, last access: 4 November 2022 (see the detailed description in Moskalik et al., 2018). In this study, sediment flux data from sediment traps deployed for 1 day at 5, 10, 15, and 20 m depths were considered for analysis. Inorganic and organic SPM concentrations (g m⁻³) and the sediment flux (g m⁻² d⁻¹) were calculated based on the total SPM or sediment flux and the loss on ignition (Moskalik et al., 2018). The integrated iSPM concentration in the water column was averaged for the main melt season (June–August 2016–2021). The annual variability in the SPM levels was visualized using a kernel density estimate (KDE) plot prepared based on SPM data (2015–2021) from discrete depths. The sinking rate of SPM (m d⁻¹) was calculated by dividing the sediment flux by the SPM concentration sampled at the corresponding depth layers (Mugford and Dowdeswell, 2011). Sediment flux and salinity datasets were used for model parametrization (see Sect. 2.3.1), whereas the SPM dataset was used for model assessment (see Sect. 2.4).

2.3 Numerical model

2.3.1 General setting

To study the dynamics of the West Spitsbergen coastal waters affected by the iSPM input, numerical experiments were designed using the polar version of the biogeochemical ECOSystem Model (ECOSMO-E2E-Polar) coupled with the General Ocean Turbulence Model (GOTM) (Burchard et al., 1999). ECOSMO-E2E-Polar represents the three main nutrient cycles (nitrogen, phosphorus, and silica) in the pelagic and sympagic systems, three functional groups of primary producers (ice algae, diatoms, and flagellates), two zooplankton groups (micro- and meso-), one macrobenthos group, and chlorophyll a as a prognostic variable, allowing a flexible chlorophyll-to-carbon ratio. The development of ECOSMO has been fully described by Benkort et al. (2020), Daewel et al. (2019), Daewel and Schrum (2013), and Yumruktepe et al. (2022). We extended the biogeochemical model to include iSPM in the model formulation (Fig. 2). The model was built with the Fortran-based Framework for Aquatic Biogeochemical Models (FABM) (Bruggeman and Bolding, 2014) to facilitate coupling with the physical model. In this first application, a 1D numerical framework was used, and physical processes in the water column were calculated by GOTM. Simulation resolved the profiles of the velocity, temperature, salinity, turbulent mixing, and transport of ecosystem state variables in 20 vertical layers with surface zooming of 1.5 and bottom zooming of 0.1. This approach neglects horizontal transport and considers vertical exchange processes only. This allows feasible parameterization, verification, and sensitivity tests to study processes with low computational effort, but it hinders the model's skill in representing advection and upwelling. However, Atlantic Water advection is considered to be limited in Hornsund in comparison to other West Spitsbergen fjords, in particular in the inner bays, whereas upwelling is most important up to a distance of 500 m from the glacier front (Pasculli et al., 2020).

Figure 2 Schematic diagram of the biogeochemical model ECOSMO-E2E-Polar with the new iSPM group. The three systems (benthic, pelagic, and sympagic) are included. Chl_d and Chl_f stand for phytoplankton-type-specific chlorophyll a content for diatoms and flagellates, respectively.

The model was implemented at 20 stations located within the newly ice-free areas in Hornsund (Fig. 1b, Table 1). The simulations were run from the beginning of 2005 to the end of 2009. Input data from 2005–2009 were averaged, repeated five times, and used as a spin-up to allow the model to reach equilibrium under the applied forcing. Temperature and salinity vertical forcing were used from the 3D hydrodynamic numerical model of Hornsund (HRM; Jakacki et al., 2017), which represents nine sources of freshwater input from the Hornsund drainage basin, including all components (ablation, precipitation, snow, and rivers). Sea ice thickness and concentration were extracted from the S800 model simulation at the closest grid cells (Albretsen et al., 2017); these are the same data as used for the HRM model. The sea ice input was smoothed using a 30 d rolling average, as the 1D setup does not represent advection and high variability of the thickness and concentration affects the performance of the ice algae module. The atmospheric conditions were prescribed from meteorological monitoring at Polish Polar Station Hornsund: i.e. the air temperature (2 m above the surface), eastward (u) and northward (v) wind speeds, cloudiness, relative humidity, and pressure. The model was run with a 30 min time step, and the daily average was saved as an output. Two sets of scenarios were performed to evaluate the gains in carbon sequestration potential due to the retreat of marine-terminating glaciers and to evaluate the losses due to the iSPM discharge by evaluating plankton primary and secondary production and carbon burial. The SPM scenario included the iSPM input prescribed to the model according to Eq. (1), and the noSPM scenario was a control run without iSPM input.

Table 1 Sources of the input data and the modelling setup.

2.3.2 ECOSMO developments

Two state variables were added to the ECOSMO-E2E-Polar model framework to account for the iSPM and the iSPM sediment pool (sed_iSPM) (Fig. 2). The input of iSPM prescribed to the model was calculated based on the inorganic sediment flux (iSF) and its relationship with air temperature and salinity (Fig. S1 in the Supplement), as developed from field data, which takes the form of Eq. (1):

$$\begin{array}{}\text{(1)}& \mathrm{iSF}={\mathrm{10}}^{\mathrm{0.04}\cdot \mathrm{6}\mathrm{accPDD}\mathrm{AT}+\mathrm{0.174}\cdot \left(\text{refS-meanS}\right)+\mathrm{0.815}},\end{array}$$

where 6accPDD AT is the accumulated daily air temperature for positive degree days for a 6 d window (°C), refS is a reference salinity for Atlantic Water (34.9) (Moskalik et al., 2018), and meanS is the mean salinity above the sediment trap. The ordinary least squares (OLS) function (in the statsmodels library in Python) was used to generate a linear model for iSF estimates. The root-mean-square deviation (RMSD) was used to measure the differences between the observed iSF and linear model. The iSF calculated for each depth layer was prescribed to the model as a daily input (C_iSPMinpput) in mg m⁻³. While some runoff data are available for Hornsund (Van Pelt et al., 2019; Błaszczyk et al., 2019), it was not feasible to parametrize the iSPM input based on the meltwater discharge here due to the structure of the hydrodynamic model (1D in contrast to 3D) and the lack of data on sediment loads in glacial plumes. However, in this study, salinity depended on the discharge provided in the 3D hydrodynamic model that was a source of input data (HRM; Jakacki et al., 2017). Therefore, we used the salinity as a proxy for the inorganic sediment.

The state variables in ECOSMO (all the state variables are listed in Table S3) are solved using prognostic equations in the form of Eq. (2):

$$\begin{array}{}\text{(2)}& C_t+\left(w_{\mathrm{d}}\right)C_z={\left(A_{\mathrm{v}}{C}_z\right)}_z+R_C,\end{array}$$

with $C_x=\frac{\mathrm{d}C}{\mathrm{d}x}$, where x represents either time (t) or depth (z). This equation includes vertical turbulent subscale diffusion, sinking rates, and chemical and biological interactions. The vertical turbulent sub-scale diffusion coefficient (A_v) is estimated by the hydrodynamic core of ECOSMO. The sinking rate (w_d) is a constant that is only non-zero for detritus, opal, and iSPM. The sinking rate (w_d) that allowed the dynamics of iSPM to be properly represented in the model was 0.8 m d⁻¹, which is in the lower range of the sinking rates observed in the field (all parameters are listed in Table 2). Chemical and biological interactions are employed in the interaction term R_C, which is different for each variable (C), based on relevant processes.

Table 2 List of parameters, corresponding description, and units used in the model.

The rate of change in the iSPM concentration (C_t term) is calculated via Eq. (3):

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{iSPM}}}{\mathrm{d}t}}=C_{\mathrm{iSPM}}+C_{\mathrm{iSPMinput}}.\end{array}$$

The interaction term R_C is calculated via Eq. (4):

$$\begin{array}{}\text{(4)}& R_{\mathrm{iSPM}}={\left[\left({\mathit{\lambda }}_{\mathrm{s}\mathrm{2}\mathrm{d}}{C}_{\mathrm{sediSPM}}-{\mathit{\lambda }}_{\mathrm{d}\mathrm{2}\mathrm{s}}{C}_{\mathrm{iSPM}}\right)/\mathrm{d}z\right]}_{z=\mathrm{bottom}}.\end{array}$$

iSPM enters a new sediment pool with a sedimentation rate (λ_d2s) of 3.5 m d⁻¹ if the bottom stress is <τ_crit or is resuspended with a resuspension rate (λ_s2d) of 26 d⁻¹ if the bottom stress is >τ_crit. The critical bottom shear stress (τ_crit) was set to 0.07 N m⁻², which is in the range reported by Wölfl et al. (2014).

As sed_iSPM exchanges occur locally at the bottom and the group is not exposed to mechanical displacement, Eq. (2) can be simplified to Eq. (5):

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{sediSPM}}}{\mathrm{d}t}}={\left[R_{\mathrm{sediSPM}}\right]}_{z=\mathrm{bottom}}.\end{array}$$

The interaction term R_C is calculated via Eq. (6):

$$\begin{array}{}\text{(6)}& R_{\mathrm{sediSPM}}=-{\mathit{\lambda }}_{\mathrm{s}\mathrm{2}\mathrm{d}}{C}_{\mathrm{sediSPM}}+{\mathit{\lambda }}_{\mathrm{d}\mathrm{2}\mathrm{s}}{C}_{\mathrm{iSPM}}.\end{array}$$

As the iSPM has an impact on light penetration, the photosynthetically active radiation in the water column had to be updated and is calculated via Eq. (7):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}I\left(x,y,z,t\right)=& {\displaystyle \frac{I_{\mathrm{s}}\left(x,y\right)}{\mathrm{2}}}\exp \\ & (-k_{\mathrm{w}}z\\ & -k_{\mathrm{Chl}}{\int }_z^{\mathrm{0}}{\sum }_{j=\mathrm{1}}^{\mathrm{2}}{\mathrm{Chl}}_{Pj}\partial z\\ & -k_{\mathrm{iSPM}}{\int }_z^{\mathrm{0}}{C}_{\mathrm{iSPM}}\partial z\\ & -k_{\mathrm{DOM}}{\int }_z^{\mathrm{0}}{C}_{\mathrm{DOM}}\partial z),\end{array}\end{array}$$

where I_s(xy) is the short-wave radiation (W m⁻²) at the surface, x and y identify the horizontal grid points, z is the water depth in meters, and k_x are extinction coefficients (Table 2).

In Hornsund, most of the variability in the optical properties in the summers of 2009 and 2010 was attributed to particles of mineral origin (Sagan and Darecki, 2018), and thus the input of organic particles with meltwater was considered negligible here. The attenuation coefficient for iSPM measured in another polar fjord in Greenland (0.13 m² g⁻¹) (Lund-Hansen et al., 2010) was high compared to other published values: 0.07 m² g⁻¹ (Christian and Sheng, 2003), 0.06 m² g⁻¹ (Pfannkuche and Schmidt, 2003), and 0.065 m² g⁻¹ (Oliver et al., 2020). Thus, here, a light extinction coefficient (k_iSPM) of 0.065 m² g⁻¹ was prescribed for the model; this gave reasonable results in terms of the light limitation and is in the range of field measurements.

The light limitation also depends on the plankton photosynthesis efficiency parameter (a). Here, it was increased to 0.04 (W m⁻²)⁻¹, which is within the range reported for Arctic coastal and shelf waters (Van De Poll et al., 2018; Stuart et al., 2000; Strom et al., 2016; Platt et al., 1982) and is in line with previous studies showing that fjord plankton communities are adapted to low light (Simo-Matchim et al., 2016; Holding et al., 2019). The light limitation is calculated as in Eq. (8):

$$\begin{array}{}\text{(8)}& \mathit{\alpha }\left(I\right)=\tanh \left(a\right)I\left(x,y,z,t\right).\end{array}$$

In this 1D setup, we do not simulate fish due to their migration, which reduces the uncertainty of the current simulations. Thus, the macrobenthos loss term only consists of excretion (ε_MBC_MB) and natural mortality (m_MBC_MB), as shown in Eq. (9):

$$\begin{array}{}\text{(9)}& R_{{\mathrm{MB}}_{\mathrm{loss}}}={\mathit{\epsilon }}_{\mathrm{MB}}{C}_{\mathrm{MB}}+m_{\mathrm{MB}}{C}_{\mathrm{MB}}.\end{array}$$

Similarly, the reaction terms for zooplankton, detritus, and dissolved organic matter (DOM) were changed accordingly to remove fish grazing (Daewel et al., 2019). Predation mortality from the fish functional group was accounted for by increasing macrobenthos natural mortality (m_MB) to 0.03 d⁻¹.

We do not provide the nutrient input with meltwater due to the lack of data for parametrization and to disentangle it from the effect of the iSPM discharge; thus, the burial rates in the carbon and nitrogen sediment pool (sedCN; Eq. 10) and in the silicate sediment pool (sedSi; Eq. 11) are set to 0 to prevent decreasing nutrient concentrations over the simulation time. For the full description of the equations, the reader is referred to Daewel and Schrum (2013). We speculate that the bias introduced by not providing the nutrient input is relatively low, considering the characteristics of the discharge (see Sect. 4.5).

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}{R}_{\mathrm{sedCN}}& ={\mathit{\lambda }}_{\mathrm{d}\mathrm{2}\mathrm{s}}{C}_D-{\mathit{\lambda }}_{\mathrm{s}\mathrm{2}\mathrm{d}}{C}_{\mathrm{sedCN}}\\ & -\mathit{\theta }\left({\mathrm{O}}_{\mathrm{2}}\right)\mathrm{2}{\mathit{\epsilon }}_{\mathrm{sedCN}}\left(T\right)C_{\mathrm{sedCN}}\\ & -\mathit{\theta }\left(-{\mathrm{O}}_{\mathrm{2}}\right){\mathit{\epsilon }}_{\mathrm{sedCNdenit}}\left(T\right)C_{\mathrm{sedCN}}-{\mathit{\delta }}_{\mathrm{bur}}{C}_{\mathrm{sedCN}},\end{array}\text{(11)}& {\displaystyle }{R}_{\mathrm{sedSi}}={\mathit{\lambda }}_{\mathrm{d}\mathrm{2}\mathrm{s}}{C}_{\mathrm{opal}}-{\mathit{\lambda }}_{\mathrm{s}\mathrm{2}\mathrm{d}}{C}_{\mathrm{sedSi}}-{\mathit{\delta }}_{\mathrm{bur}}{C}_{\mathrm{sedSi}}.\end{array}$$

The carbon burial potential (CB; Eq. 12) was calculated as a 70 % burial efficiency of the carbon and nitrogen sediment accumulation rate, as previously reported for Hornsund (Koziorowska et al., 2018):

$$\begin{array}{}\text{(12)}& \mathrm{CB}={\mathit{\eta }}_{\mathrm{bur}}{R}_{\mathrm{sedCN}}.\end{array}$$

2.3.3 Model assessment

Satellite data products of suspended particulate matter were not available for the glacial bays, and analysis of the long-term trends and model validation (2005–2009) were not possible here. Thus, we performed the model assessment based on the available field data from 2015–2021. The summertime mean iSPM concentration for past conditions was reconstructed based on measurements of the iSPM concentration and PDD AT, and it showed a high correlation with the simulated iSPM concentration parametrized based on the complementary dataset of sediment flux measurements (Fig. S2, R² = 0.928, p = 0.009). The iSPM concentrations at modelled station 2 (HH1) in 2006 and 2009 were also compared with the iSPM field data at monitoring stations M4 (H1_09) and M5 (H1_11) from 2019 (Fig. S3), which represented the closest environmental conditions (PDD SST, PDD AT, melt season duration, and precipitation in Fig. 3) to those in simulation period. Results showed that the model realistically simulated the seasonal pattern and vertical distribution of the iSPM (rho > 0.74, p<0.001 for Spearman's correlation; see Table S4). Despite the fact that the iSPM input was parametrized for Hansbukta (which was the only bay with sufficient data and is the most studied bay in Hornsund), and the iSPM load and discharge can differ between glaciers, the spatial patterns from measurements of iSPM at the surface conducted in all of Hornsund in summer 2017 were in line with the simulation results (Fig. S4). The literature data (Table S5) on the concentrations of all the nutrients and functional groups showed that the model performed well when compared to the current knowledge of the West Spitsbergen fjords.

Figure 3 Long-term trends in the newly ice-free marine habitat and the melt potential in Hornsund: area (km²) and volume (km³) of the newly ice-free marine habitat (assessed for summers between 1976 and 2022), sea-ice-free days (fraction of the year, 1982–2020), accumulated positive degree days for the sea surface temperature and air (PDD SST and PDD AT, °C, 1979–2022), inorganic SPM concentration reconstructed from 6 years of monitoring (g m⁻³, 1979–2022), melt season duration (fraction of the year, 1979–2022), and precipitation (mm, 1979–2020). ^* p<0.05, ${}^{**}$ p<0.001 for the modified Mann–Kendall test. Grey shading indicates the modelling period (2005–2009).

2.4 Data analysis and visualization

The maps and satellite images were generated and processed in ArcGIS Pro 2.8.0. The plots were prepared in Python 3.7 (Van Rossum and Drake, 2009) using Matplotlib 3.1.1 (Caswell et al., 2019), Pandas 1.0.5 (Mckinney, 2010; Reback et al., 2020), and Seaborn 0.11.1, and they were arranged in Inkscape 0.92.4.

The Hamed and Rao modified Mann–Kendall (mMK) test was used to determine whether a trend existed in time series data (SIF, PDD SST, PDD AT, iSPM, and precipitation) with a significance level of 0.05 or 0.001 (Python library pymannkendall 1.4.2).

For each modelled station and each scenario, the 5-year (2005–2009) averages of SIC and SIT in May, the mean summertime integrated iSPM, and the rates of phytoplankton primary production (phyPP), zooplankton secondary production (zooSP), and carbon burial (CB) were calculated. Then, the average values of the phyPP, zooSP, and CB rates for all 20 stations were multiplied by the average newly ice-free area between 2006 and 2010 (64.21 km²). The resulting phyPP, zooSP, and CB under the SPM scenario were considered as gains in carbon sequestration potential due to the retreat of marine-terminating glaciers, whereas the differences between the noSPM and SPM scenarios were considered as losses due to the iSPM discharge with meltwater.

The influence of the iSPM discharge on the ecosystem dynamics was exemplified by presenting the biomasses of ice algae (IA) and macrobenthos (MB) as well as the biomasses of phytoplankton (PHY) and zooplankton (ZOO) and the silicate, phosphorus, nitrogen, and light limitation index (SIL, PLI, NLI, and LLI) integrated across the whole water column at three modelled stations (2, 9, and 14) that were comparable due to their similar depths (42.45–49.55 m) but presented low, intermediate, and high levels of summertime iSPM input. Also, 2 years with contrasting sea ice conditions (2008 and 2009) were displayed.

3 Results

3.1 Newly ice-free marine habitats

The area of newly ice-free coastal waters due to the retreat of marine-terminating glaciers in Hornsund increased by ∼ 99.4 km² between the summers of 1976 and 2022 (an increase in total area of around 38 %; see Figs. 1 and 3), whereas the volume gained by 2010 was ∼ 3.3 km³. The trends were linear ($y=\mathrm{2.1406}x-\mathrm{4231.2}$; R²=0.995 for the area and $y=\mathrm{0.097}x-\mathrm{191.38}$; R²=0.984 for the volume; t-test p<0.001), with rates of ∼ 21.4 km² per decade and 1.0 km³ per decade (Fig. 3). While advances of glacier fronts due to surge events were observed for some marine-terminating glaciers in Hornsund, these did not influence the overall increasing trends. Along with the glacial retreat, the number of SIF days (the fraction of the year with SIC < 15 %) increased significantly (∼ 0.1 per decade, i.e. by around 44 d; p<0.001, mMK test). Despite high interannual variability, the central part of Hornsund has become mostly devoid of sea ice since 2006, but there still is seasonal sea ice cover in the newly formed glacial bays (Fig. 6a).

3.2 Melt and SPM discharge potential

The annual sum of daily SST > 0 °C (PDD SST) showed no significant trend in outer Hornsund due to strong variability between years (p>0.1, mMK test) (Fig. 3), but it significantly increased for the summer months (June–August; 46.8 °C per decade; p<0.05, mMK test). The annual and summertime sums of positive daily air temperatures (PDD ATs) as well as the annual precipitation showed significant increases (60.5, 31.4, and 56.0 mm per decade, respectively; p<0.001, mMK test). Melt season duration increased significantly (p<0.001 mMK test) at a rate of ∼ 9 d per decade (2.5 % of the year). At the beginning of the measurements, the melt season started in June and ended in late September to mid-October, whereas it can currently start as early as February and ends mostly in October (Fig. S5).

The 6-year monitoring dataset of summertime SPM concentration in Hansbukta (Fig. 1c) was not sufficient to show long-term trends. However, average integrated iSPM levels were correlated with both the annual sum of PDD AT ($y=\mathrm{0.061}x-\mathrm{19.549}$, R²=0.68; p<0.05, t test) and the summertime sum of PDD AT (June–August) ($y=\mathrm{0.221}x-\mathrm{75.047}$, R²=0.78; p<0.05, t test). Even though the correlation was stronger for the summertime PDD AT, the estimates displayed numerous negative values. However, the annual sum of PDD AT allowed a coarse reconstruction of past conditions and revealed significant increases in iSPM concentration (3.7 g m⁻³ per decade in 1979–2022; p<0.001, mMK test). Importantly, within the modelled time range (2005–2009; Fig. 3, grey shading), both iSPM estimates gave similar results in 2006 and 2009 (8.6 and 12.0; 8.1 and 9.8 g m⁻³, respectively).

3.3 SPM dynamics

The concentration of iSPM varied between seasons, with the highest levels occurring in July–October (up to 150 g m⁻³) and the lowest between November and May (up to 50 g m⁻³), whereas the highest levels of organic SPM were observed between April and June (up to 20 g m⁻³) (Fig. 4a). Sediment flux observed for iSPM ranged between 1–6648 g m⁻² d⁻¹, while it was 0.9–333 g m⁻² d⁻¹ for organic SPM (Fig. 4b). The sinking rate of iSPM ranged between 0.6–265 m d⁻¹ (mean 25.3, median 12.2 m d⁻¹) (Fig. 4c), while the sinking rate of organic SPM was 1 order of magnitude lower, with a range of 0.3–28.9 m d⁻¹ (mean 2.8, median 1.7 m d⁻¹). The sediment flux of iSPM, which represents the temporary dynamics of the iSPM input, was dependent on the accumulated daily air temperature for positive degree days for a 6 d window (6accPDD AT) and mean salinity in the layer above (R²=0.662; p<0.001, t test; Fig. 4d). Within the range of frequently observed values of 6accPDD AT (0–40 °C) and salinity (30–35), the estimated inorganic sediment flux could reach 1860 g m⁻² d⁻¹ (Fig. S1). Importantly, the regression model (Eq. 1) performed well for an inorganic sediment flux of < 2000 g m⁻² d⁻¹ (RMSD = 290.1 g m⁻² d⁻¹), which corresponded to 95 % of the dataset, and mostly underestimated the highest inorganic sediment flux values (RMSD = 823.3 g m⁻² d⁻¹ for all the dataset).

Figure 4 SPM dynamics at monitoring stations in Hansbukta in 2016–2021. (a) Kernel density estimates of SPM concentration (g m⁻³; inorganic on the left and organic on the right). Colours indicate the distribution between months (high to low). (b) Inorganic (blue) and organic (yellow) sediment fluxes (g m⁻² d⁻¹; the dashed grey line is the mean value). (c) Inorganic and organic matter sinking rate (m d⁻¹; the dashed grey line is the mean value). (d) Inorganic and organic sediment fluxes (SF; dots, log scale) and the accumulated daily air temperature for positive degree days for a 6 d window (°C; 6accPDD AT; line).

Figure 5 Spatial patterns in (a) average sea ice thickness (SIT, cm) and concentration (SIC, %) in May and the summertime integrated inorganic SPM concentration (iSPM, g m⁻³), (b) plankton primary production (phyPP), (c) plankton secondary production (zooSP), and (d) carbon burial (CB) (g C m⁻² yr⁻¹), with carbon gains (in blue) due to the retreat of marine-terminating glaciers (in green; SPM scenario) and losses due to the inorganic SPM discharge with meltwater (in pink; noSPM−SPM scenario, Gg C yr⁻¹, average for 2005–2009) also shown. The lines indicate the extent of newly ice-free areas in 2006. The ecosystem dynamics at stations 2, 9, and 14 are presented in Fig. 6. Landsat 8 satellite images (from 14 May 2022 for a and from 4 August 2020 for b, c, d) were downloaded from https://glovis.usgs.gov/app, last access: 29 March 2023.

3.4 Spatial patterns of sea ice, iSPM, plankton production, and carbon burial

The mean SIT and SIC in May were highest in the southern and inner parts of Hornsund (the 5-year average was up to 8 cm and 19.3 %, respectively) and lowest in the northern and outer parts (5-year average was 3 cm and 8.7 %, respectively) (Fig. 5a). The mean summertime integrated iSPM concentration was highest in the inner glacial bay (modelled station 14; 5-year average: 164.4 g m⁻³), where rates of plankton primary and secondary production and carbon burial were lowest (5-year average: 11.0, 1.5, and 5.5 g C m⁻² yr⁻¹, respectively; Fig. 5b, c, d). At other stations, the mean summertime integrated iSPM concentration was in the range between 2.1–7.1 g m⁻³, which allowed phyPP to reach rates between 66.3–100.7 g C m⁻² yr⁻¹ (versus 131.3–171.2 g C m⁻² yr⁻¹ under the noSPM scenario), whereas rates of zooSP were between 17.7–47.2 g C m⁻² yr⁻¹ (versus 48.7–75.7 g C m⁻² yr⁻¹ under the noSPM scenario), and the CB rate was in the range of 6.0–17.7 g C m⁻² yr⁻¹ (versus 6.5–23.0 g C m⁻² yr⁻¹ under noSPM). In the simulation period (2005–2009), the newly ice-free areas in Hornsund substantially contributed to phyPP, zooSP, and CB (on average 5.1, 2.0, and 0.9 Gg C yr⁻¹, respectively; see the gains in carbon sequestration potential in the SPM scenario in Fig. 5, which are shown in green). However, the potential was hindered by 5.0, 2.1, and 0.1 Gg C yr⁻¹, respectively, due to the iSPM input (see the losses due to the difference between the noSPM and SPM scenarios in Fig. 5, which are shown in pink). Thus, without the release of mineral particles, plankton primary and secondary production could have been around 2 times higher (10.1 and 4.1 Gg C yr⁻¹, respectively, under the noSPM scenario), whereas carbon burial was less affected by iSPM input (1.0 Gg C yr⁻¹ under the noSPM scenario – around 16.5 % higher than carbon burial under the SPM scenario).

3.5 Ecosystem dynamics

The ecosystem dynamics relating to the sea ice and iSPM in the newly ice-free areas are presented for three modelled stations with a low, intermediate, and high influence of iSPM (stations 9, 2, and 14; Fig. 6a, b, c, respectively) in 2 contrasting years (cold 2008 and warm 2009). The sea ice thickness and concentration (SIT and SIC) were lower in 2008 than in 2009 and in the outer than in the inner glacial bay. Thus, only in 2009 did the ice algae bloom reach a biomass of 0.16 g C m⁻² in inner Hornsund (stations 9 and 14 in Fig. 6a and c, respectively), and the presence of sea ice (up to 0.5 m) delayed the phytoplankton bloom by around 10 d. Under a low or intermediate influence of iSPM (Fig. 6a, b), the light limitation index was slightly lowered before the main melt season (March to early June), whereas a significant effect of light limitation due to the iSPM input started around late June (up to 24 and 6 g m⁻³ at stations 9 and 2 in Fig. 6a and b, respectively). Due to the worsened underwater light conditions, the peaks in the spring and summer phytoplankton blooms were delayed by around 10–14 d and the summer peak reached a lower biomass (∼ 0.4–0.5 g C m⁻³ under the SPM scenario and 0.7 g C m⁻³ under the noSPM scenario), which also affected the zooplankton (the peak was delayed by ∼ 9 d and there was 0.1–0.2 g C m⁻³ less biomass) and macrobenthos (∼ 10 g C m⁻² less biomass). At the highest levels of iSPM (up to 500 g m⁻³ at station 14; Fig. 6c), strong light limitation started early in March. Thus, phytoplankton, zooplankton, and macrobenthos reached very low biomasses (< 0.2 g C m⁻³, < 0.05 g C m⁻², and < 5 g C m⁻², respectively). The iSPM-related delays in phytoplankton blooms led to delays in silicate limitation and increases in the ice algae biomass in spring, particularly at the station with the highest levels of iSPM (the difference between the SPM and noSPM scenarios was up to 0.01 g C m⁻²).

Figure 6 Changes in ecosystem dynamics due to the iSPM input in 2008 and 2009 at three modelled stations: (a) station 9 (low iSPM influence), (b) station 2 (medium iSPM influence), and (c) station 14 (high iSPM influence). Line plots show the sea ice thickness (SIT; m), the biomass of ice algae (IA; g C m⁻² or mg C m⁻²; black), the integrated silicate (grey), phosphorus (red), and nitrogen (blue) and the light limitation index (SLI, PLI, NLI, and LLI, respectively; all dimensionless), the integrated biomasses of phytoplankton and zooplankton (PHY and ZOO, respectively; g C m⁻³), and the macrobenthos biomass (MB; g C m⁻²). Full line – SPM scenario; dashed line – noSPM scenario. SLI and LLI values equal to 1 indicate that phytoplankton is not limited by either silicate or light. Colour plots indicate where SIC < 15 % (blue: open water), differences in the ice algae, phytoplankton, zooplankton, and macrobenthos biomasses between the SPM and noSPM scenarios (SPM−noSPM), and the integrated inorganic SPM concentration in the SPM scenario (iSPM; g m⁻³). Brown arrows indicate the start of the melt season (30 May 2008 and 3 June 2009), and black arrows indicate delays in the peak abundances of phytoplankton and zooplankton. Note the different scales (^*b for panel b and ^*c for panel c).

4 Discussion

4.1 Newly ice-free marine habitats

We report significant increases in new marine habitat area (∼ 100 km²) and volume (> 3.3 km³) between 1976 and 2022 in Hornsund (Figs. 1b, 3, 7a) due to the retreat of marine-terminating glaciers. These results are in line with cryosphere studies in West Spitsbergen fjords (Błaszczyk et al., 2021, 2023; Grabiec et al., 2018; Strzelecki et al., 2020) and in polar regions in general (Kochtitzky et al., 2022; Pfeffer et al., 2014). In the coastal Arctic and Antarctic, glaciers and ice sheets have lost mass due to increased submarine (basal) melting and iceberg calving (dos Santos et al., 2021; Błaszczyk et al., 2013, 2023), and a doubling of ice mass loss in Svalbard by 2100 was predicted (Geyman et al., 2022). The retreat of many marine-terminating glaciers has already produced newly ice-free areas, and some of them have receded onto land (Błaszczyk et al., 2013; Jerosch et al., 2019; Kochtitzky et al., 2022). Recently, the rapid loss of numerous glaciers was related to both external forcing such as increases in atmospheric and oceanic temperatures and a lack of sea ice buttressing or internal dynamics such as surges (Błaszczyk et al., 2013, 2023; Strzelecki et al., 2020). Here, we show increasing trends in the length of the melt season (∼ 9 d per decade) and the sum of PDD (46.79 °C per decade for summer PDD SST and 60.54 °C per decade and 31.43 °C per decade for annual and summer PDD AT, respectively; Fig. 3). While the melting potential is rising, the annual runoff in Svalbard is expected to increase till 2060; then it will likely decrease towards 2100 due to the reduction in glacier storage as the glaciers shrink (Bliss et al., 2014; Van Pelt et al., 2021; Nowak et al., 2021).

Furthermore, we report a significant loss in sea ice duration in central Hornsund (∼ 44 d per decade) (Fig. 3). However, as glacial retreat opens new coastal areas, it also increases the potential for winter sea ice formation in the more protected inner bays (Fig. 5a). In contrast to the glaciers where mass loss cannot be stopped nor reversed once induced, sea ice was shown to be more responsive to variations in both ocean and air temperatures (Muckenhuber et al., 2016). Thus, there can still be land-fast ice (sea ice attached to the coastline) that covers the inner parts of West Spitsbergen fjords for a limited time during winter and spring. Moreover, the ice bridge in inner Hornsund (Figs. 1b, 7a) is predicted to vanish in the coming decades (2030–2055) (Grabiec et al., 2018; Osika et al., 2022), which will transform Hornsund from a fjord into a strait, enabling sea ice advection from the Barents Sea. However, the loss of the ice bridge could also result in the increased presence of warm Atlantic Water in the area and, therefore, further sea ice loss. These seemingly contrasting predictions highlight the importance of continuous evaluation of the changing Hornsund environment and its potential as a model area for studies on regime shifts.

4.2 SPM dynamics

Based on the coarse reconstruction and modelling results presented in this study, we suggest that the Hornsund bays have been under the strong influence of dark glacial plumes since the beginning of the simulation period (2005) (Figs. 3, 5a, b, 6). In this study, the reconstructed iSPM concentration increased after 2013, and further rises are expected (3.7 g m⁻³ per decade integrated for the water column in summer). We show that air temperature variability – specifically, the accumulated daily air temperature above the melting point for a 6 d window (6accPDD AT), which takes into account the delays in meltwater discharge – modulates the iSPM flux (Figs. 4d, S1), as was also suggested for a glacial bay in Kongsfjorden, another West Spitsbergen fjord (D'Angelo et al., 2018). Recent studies also indicated that sediment production and fluxes to the coastal zones in the polar regions have increased due to higher air temperatures (Overeem et al., 2017; Szczuciński et al., 2009). Thus, it is anticipated that even central fjords will receive high inputs of mineral particles in the future, as turbid glacial plumes will spread farther from the source (Fig. 7a) (Castelao et al., 2019; Kanna et al., 2018; Hudson et al., 2014), thereby extending the influence of meltwater discharge onto the shelf and considerably affecting marine systems downstream (Meire et al., 2015, 2017; Milner et al., 2017). The relationship between melting potential and sediment input might differ between catchments; in particular, it could change after glaciers retreat onto land.

Figure 7 (a) 3D representations of the inner Hornsund bay (Brepollen) in the summers of 1976 and 2022. Landsat satellite images (from 18 July 1976 and 15 August 2022) were downloaded from https://glovis.usgs.gov/app, last access: 3 October 2022. Digital elevation model data were downloaded from https://arcticdem.apps.pgc.umn.edu/, last access: 17 August 2023. (b) Schematic representation of the positive (blue arrows) and negative (red arrows) feedback mechanisms influencing biological production and carbon burial in the Arctic fjords.

The iSPM discharge was most extensive during summer (Fig. 4), although it could also be observed during autumn and winter, when it is intensified by tidal resuspension, resulting in a relatively high concentration of organic and inorganic suspended particles (Moskalik et al., 2018). In the future, more days with open-water conditions (no sea ice), which can increase wave action and particle removal from the beaches and tidal flats, as well as a longer melt season (Figs. 3, S5) could potentially lead to iSPM affecting a substantial part of the productive season, including not only summer and autumn but also spring. Here, we show high variability of the iSPM dynamics, with sinking rates between 0.6–265.9 m d⁻¹ and sediment fluxes between 1.0–6647.7 g m⁻² d⁻¹ (Fig. 4b, c), which should be investigated further in the context of the driving mechanisms, such as flocculation (Moskalik et al., 2018).

4.3 Ecosystem dynamics

Observational data and previous modelling studies have shown that the continuing retreat of marine-terminating glaciers will negatively affect planktic and benthic communities, especially in enclosed shallow bays such as Brepollen (Fig. 7a) (Neder et al., 2022; Szeligowska et al., 2021, 2022; Torsvik et al., 2019). Indeed, we observed decreases in phytoplankton, zooplankton, and macrobenthos biomasses and delays in their peak occurrences close to the glacial fronts (by around 10–14 d as compared to the noSPM scenario; Fig. 5). These decreases were related to the input of particulate matter from land, which, even at a relatively low concentration in spring, can affect phytoplankton due to light attenuation (Fig. 7b). Under the SPM scenario, plankton primary production rates reached 66.3–100.7 g C m⁻² yr⁻¹ with a mean summertime integrated iSPM concentration of 2.1–7.1 g m⁻³, whereas it was around 2 times higher in the noSPM scenario (131.3–171.2 g C m⁻² yr⁻¹). Both ranges are comparable with the field measurements in inner and outer Hornsund and other West Spitsbergen fjords (Hodal et al., 2012; Iversen and Seuthe, 2011; Piwosz et al., 2009; Vonnahme et al., 2021) (Table S5, Fig. S6). Sea ice algae biomass was extremely low in most years (< 12 mg C m⁻², except for 2009, when it was up to 160 mg C m⁻²) due to thin ice (< 50 cm) that disappeared before the main productive season. Ice algae did not seem to be negatively affected by iSPM and, as modelling results suggest, their biomass was slightly higher in the SPM scenario than in the noSPM scenario. Importantly, we suggest that sea ice loss leading to the earlier onset of spring pelagic production might become a compensation mechanism for the higher iSPM input in summer (Fig. 7b).

The modelled carbon burial rate was within the reported values (Koziorowska et al., 2018; Kuliński et al., 2014; Zaborska et al., 2018) and constituted around 10 %–20 % of the primary production, which is also in line with the current observations for polar and sub-polar fjords (Włodarska-Kowalczuk et al., 2019). Unfortunately, no field data for the assessment of plankton secondary production rates exist from this region. However, the values simulated here fell between the plankton primary production and carbon burial rates, as expected. Plankton secondary production was reduced due to the decreased food base (17.7–47.2 g C m⁻² yr⁻¹ versus 48.7–75.7 g C m⁻² yr⁻¹ under the SPM and noSPM scenarios, respectively). According to our simulations, carbon burial was the least affected by iSPM (6.0–17.7 g C m⁻² yr⁻¹ versus 6.5–23.0 g C m⁻² yr⁻¹ under the SPM and noSPM scenarios, respectively). Since the burial of accumulated material depends on the vertical flux of the organic matter originating from phytoplankton and zooplankton, food intake by benthic fauna, and rates of benthic mineralization, we hypothesize that the changes in the phytoplankton bloom timing might have shifted the carbon pathway from the zooplankton and macrobenthos pool to carbon burial in sediments, and thus carbon burial was still relatively high in the SPM scenario (∼ 16 % lower than under the noSPM scenario). Only the extremely high levels of iSPM (a mean summertime integrated iSPM concentration of 164.4 g m⁻³), which can be observed directly inside the turbid plumes, resulted in an almost complete absence of phyto- and zooplankton and macrobenthos and relatively low plankton production rates (11.0 and 1.5 g C m⁻² yr⁻¹ for phyPP and zooSP, respectively) but still considerable burial rates (5.5 g C m⁻² yr⁻¹). Thus, we speculate that sediment discharge to polar coastal zones might result in less complex food webs constituted by species that are better adapted to high iSPM concentrations and sedimentation rates, as shown for Antarctic benthos (Clark et al., 2013; Krzemińska and Kukliński, 2018; Sahade et al., 2015). This could reduce the biomass that is utilized in the pelagic and benthic system, leading to higher carbon burial in sediments (Fig. 7b).

4.4 Carbon gains and losses

Marine sediments in polar fjords have recently been recognized as efficient organic carbon sinks and incorporated into global carbon burial estimates (Bianchi et al., 2020; Cui et al., 2022; Smith et al., 2015), highlighting their societal importance as a climate-regulating ecosystem service (Barnes et al., 2021; Bax et al., 2021). They might become more efficient at the capture and long-term carbon storage due to high sedimentation rates and their restrictive nature compared to more open coastal environments, particularly with the expansion of shallow and isolated bays and increased land–ocean connectivity (Fig. 7b) (Smith et al., 2015). Here, we show that newly ice-free areas in Hornsund (∼ 64 km² between 2006–2010) markedly contributed to plankton primary (5.1 Gg C yr⁻¹) and secondary (2.0 Gg C yr⁻¹) production and carbon burial (0.9 Gg C yr⁻¹) (green in Fig. 6, indicating a carbon gain under the SPM scenario). This carbon burial constitutes only a small fraction of the globally estimated rates for the seafloor (2.9 × 10⁴–1.6 × 10⁵ Gg C yr⁻¹; Bauer et al., 2013; Cai, 2011; Hedges and Keil, 1995). However, emerging marine habitats could gain more relevance considering that the organic carbon burial efficiency in fjords is 2 times higher than the global ocean average (Smith et al., 2015) and recognizing the scale of marine ice loss across the Arctic and Antarctic. Due to the anticipated negative effects of glacier ice loss (Hunter, 2022), here we show that part of the potential gains in carbon sequestration related to the newly ice-free areas turns into losses of plankton primary (−5.0 Gg C yr⁻¹) and secondary (−2.1 Gg C yr⁻¹) production and burial (−0.1 Gg C yr⁻¹) under the SPM scenario (Fig. 6, red). Without the delivery of mineral particles from land, plankton primary and secondary production could be around 2 times higher (10.1 and 4.1 Gg C yr⁻¹ under the noSPM scenario, respectively), whereas carbon burial is less affected by iSPM input (1.0 Gg C yr⁻¹ under the noSPM scenario). Importantly, the carbon burial efficiency is highly variable and differs between fjords (Koziorowska et al., 2018), thus limiting direct generalizations.

4.5 Current limitations and future perspectives

While the coupled physical–biogeochemical model with the newly implemented iSPM input performed well according to our assessment, the field data for model parametrization and validation were not available for the simulated period (2005–2009), whereas remotely sensed products for iSPM concentration did not cover the inner fjords and were frequently limited by clouds. Despite that, reconstructions of previous conditions and an assessment based on two complementary datasets (SPM and sediment flux) collected in recent years (2015–2021) suggest that the simulated spatial and temporal dynamics of both inorganic and organic SPM are rather realistic and in line with the current knowledge of the West Spitsbergen marine ecosystem. While it should be considered that this reconstruction was based on a few years of measurements, which might limit its robustness, particularly towards the beginning of the AT measurements, the correlation with annual PDD AT seems to yield reasonable estimates. A recent multi-year study (2010–2016) of another West Spitsbergen fjord (Kongsfjorden) also indicated the relationship between particle fluxes and air temperature above the melting point (D'Angelo et al., 2018). Importantly, satellite data products calibrated for the glacial bays should become available (Klein et al., 2021; Walch et al., 2022) to verify the long-term trends in the iSPM discharge revealed in this study.

In this study, we used the meteorological forcing from observations performed at the Polish Polar Station located in the outer part of Hornsund for all the modelled stations, since there is no long-term weather monitoring in the inner fjord. A previous study showed that in summer, the air temperature in the inner fjord was lower by 0.6–1 °C than values reported for the Polish Polar Station, and the highest difference was observed during winter (around 2 °C) (Araźny et al., 2018). While the proper atmospheric representation is crucial and, in general, the spatial variations could affect the result, the differences in daily temperature (AT) and precipitation between the inner and outer fjord according to atmospheric fields derived from ERA-Interim reanalysis (Dee et al., 2011) were relatively low (Fig. S7). This could be related to the fact that Hornsund is a small fjord and its opening is mostly influenced by the Sørkapp Current, which transports Arctic Water from the Barents Sea. The polar front that exists there reduces the advection of warm Atlantic Water into Hornsund. Thus, the entire area retains Arctic character, in contrast to other West Spitsbergen fjords (Promińska et al., 2017; Cisek et al., 2017).

So far, only a few numerical models have been implemented in the polar and subpolar regions to study the dynamics of the SPM input from land. 3D models have indicated the areas with long residence times and high accumulation rates of iSPM (Neder et al., 2022) and considered a light limitation that led to a shallowing of the photic zone within the dark plumes (Le Fouest et al., 2010; Marín et al., 2013; Møller et al., 2023). Moreover, 2D models have been developed to simulate the sedimentation induced by ice-rafted debris (Mugford and Dowdeswell, 2010) and by glacial meltwater plumes emerging from marine-terminating glaciers (Dowdeswell et al., 2015; Mugford and Dowdeswell, 2011). However, they differ from each other and from our study in the parametrization and approach due to the various data used, processes represented, and numerical models available for the respective regions. To the best of our knowledge, this is one of the first attempts to implement the influence of iSPM in the coupled physical–biogeochemical model in polar coastal zones. While this 1D approach does not represent upwelling or spatial fluctuations in the glacial plumes (e.g. those implemented for Hansbukta in 2D; De Andrés et al., 2018, 2021) or the flocculation of the particles (Dowdeswell et al., 2015; Mugford and Dowdeswell, 2011), it is a first step in addressing the technical challenges related to the coupling between the sympagic, pelagic, and benthic systems and their responses to glacial discharge and retreat.

Even though Hornsund is amongst the most-studied Svalbard fjords, our study was limited to a 5-year period due to the lack of long-term input data for temperature and salinity (Jakacki et al., 2017; Torsvik et al., 2019), as most of the hydrodynamic models do not simulate coastal zones with sufficient horizontal resolution and do not consider changes in glacial bay extent. It should also be considered that the sea ice concentration and thickness were extracted from the closest data points available, and thus the sea ice conditions might have been different in the glacial bays. However, smoothing the data for more stable model runs could have resulted in more accurate forcing. Also, the advection of Atlantic Water is not represented in this 1D setup; however, due to the strong boundary in the form of the polar front and sills, most of the primary and secondary production in the glacial bays of Hornsund is assumed to be local, in contrast to other West Spitsbergen fjords, which experience high advection of plankton (Basedow et al., 2004; Gluchowska et al., 2016).

Here, we disentangled the effects of the iSPM input from those of other factors such as organic matter and nutrient delivery with meltwater. The influence of terrestrial organic matter on light attenuation was assumed to be negligible in Hornsund for the time covered by the simulation (Sagan and Darecki, 2018). Despite that, the release of large amounts of petrogenic organic carbon that has been isolated for millennia under the ice has recently emerged as an important component of carbon burial in fjords, and the fluxes of this carbon – as well as the transformations of it by microorganisms, leading to greenhouse gas emissions – should be better constrained for future model development (Fig. 7b) (Ruben et al., 2023). Moreover, several modelling and field studies in Arctic coastal waters have shown that the upwelling effect of submarine plumes and nutrient fluxes with meltwater supports primary production in the glacial bays and on the shelf (Castelao et al., 2019; Luo et al., 2016; McGovern et al., 2020; Oliver et al., 2020). However, the net effect depends on the lithology, the subglacial discharge rate, and the depth of the glacier grounding line as well as the seasonal dynamics of coastal currents, winds, and eddy activity, and it was not possible to represent these properly in this study. Studies in deep Greenland fjords indicate that macronutrients were primarily supplied to the surface waters by mixing and not through transport from land with glacial meltwater, as this was shown to have a relatively low nutrient load (Hopwood et al., 2020). However, the Svalbard fjords are relatively shallow, and thus the upwelling pump might not be as efficient as for Greenland fjords or the shallower, nutrient-deficient waters might be transported (Hopwood et al., 2018). Furthermore, while macronutrient concentrations can be higher in the Arctic rivers, most of the discharge in Hornsund comes from marine-terminating glaciers (Błaszczyk et al., 2019). Also, rivers were shown to deliver nutrients mostly in August (McGovern et al., 2020), when phytoplankton are already limited due to the light attenuation by iSPM. Even though nutrient input was not provided per se, setting the nutrient burial rate to 0 allowed nutrients that would otherwise be excluded to be kept in the system, and it could, to some degree, compensate for the lack of nutrient input with meltwater. Thus, the overall bias introduced by not providing nutrient input in our simulations might be relatively low.

The ecosystem dynamics are a result of the combined interaction of, inter alia, the dynamic coastline, hydrographic and sea ice conditions, nutrients, and sediment discharge, and thus this interdisciplinary work adds to the current understanding of the complex influence of glaciers on marine productivity and carbon fluxes (Hopwood et al., 2020). The presented numerical framework allows us to disentangle the effects of various processes and permits efficient hypothesis testing. Despite inherent weaknesses, it provides reliable results that are comparable with field measurements. The limitations of this study could be readily addressed by the further development and implementation of high-resolution general circulation models in polar regions (Szeligowska et al., 2024) and by coupling them with biogeochemical modules such as those presented here. Thus, skilful 3D fine-scale ecosystem models could arise from such work in the future.

5 Conclusions

In this study, we used Hornsund as a model high-latitude fjord that is particularly sensitive to a changing climate. We presented the accumulated effects of interactions between the atmosphere, ocean, cryosphere, and the dynamic coastline and explored how these affect the carbon sequestration potential. By combining the results of numerical modelling, remote sensing, and in situ observations, we provided a broad view of the periglacial environment and a framework for future simulations of ecosystem dynamics affected by terrigenous matter input with meltwater. Relatively well-studied areas adjacent to rapidly retreating marine-terminating glaciers in Hornsund are representative of similar coastal environments with shallow grounding line depths and, therefore, shed light on the formation and development of new marine habitats on not only a local but also a regional scale. Here, we show that despite the negative influence of iSPM input, the loss of marine ice in polar regions can be expected to ultimately lead to higher net productivity and the emergence of carbon sinks due to the formation of newly ice-free areas. Thus, glacial retreat and the terrigenous matter input should be implemented in current ocean models applied to such coastal systems to resolve carbon fluxes more accurately. However, the intertwined complexity of changes in High-Arctic coastal zones complicates the estimation of net effects on carbon burial in sediments. Considerable uncertainties remain, in particular those related to the petrogenic organic carbon release. Here, we also highlight the importance of maintaining long-term observations and implementing the FAIR (findability, accessibility, interoperability, and reusability) principles in data infrastructures to improve our understanding of the evolution of deglaciating coasts and subsequent influences on the marine ecosystem, which is one of the research priorities in the context of climate change impacts on polar regions.