Benthic silicon cycling in the Arctic Barents Sea: a reaction–transport model study

Over recent decades the highest rates of water column warming and sea ice loss across the Arctic Ocean have been observed in the Barents Sea. These physical changes have resulted in rapid ecosystem adjustments, manifesting as a northward migration of temperate phytoplankton species at the expense of silica-based diatoms. These changes will potentially alter the composition of phytodetritus deposited at the seafloor, which acts as a biogeochemical reactor and is pivotal in the recycling of key nutrients, such as silicon (Si). To appreciate the sensitivity of the Barents Sea benthic system to the observed changes in surface primary production, there is a need to better understand this benthic–pelagic coupling. Stable Si isotopic compositions of sediment pore waters and the solid phase from three stations in the Barents Sea reveal a coupling of the iron (Fe) and Si cycles, the contemporaneous dissolution of lithogenic silicate minerals (LSi) alongside biogenic silica (BSi), and the potential for the reprecipitation of dissolved silicic acid (DSi) as authigenic clay minerals (AuSi). However, as reaction rates cannot be quantified from observational data alone, a mechanistic understanding of which factors control these processes is missing. Here, we employ reaction–transport modelling together with observational data to disentangle the reaction pathways controlling the cycling of Si within the seafloor. Processes such as the dissolution of BSi are active on multiple timescales, ranging from weeks to hundreds of years, which we are able to examine through steady state and transient model runs.

Steady state simulations show that 60 % to 98 % of the sediment pore water DSi pool may be sourced from the dissolution of LSi, while the isotopic composition is also strongly influenced by the desorption of Si from metal oxides, most likely Fe (oxyhydr)oxides (FeSi), as they reductively dissolve. Further, our model simulations indicate that between 2.9 % and 37 % of the DSi released into sediment pore waters is subsequently removed by a process that has a fractionation factor of approximately −2 ‰, most likely representing reprecipitation as AuSi. These observations are significant as the dissolution of LSi represents a source of new Si to the ocean DSi pool and precipitation of AuSi an additional sink, which could address imbalances in the current regional ocean Si budget. Lastly, transient modelling suggests that at least one-third of the total annual benthic DSi flux could be sourced from the dissolution of more reactive, diatom-derived BSi deposited after the surface water bloom at the marginal ice zone. This benthic–pelagic coupling will be subject to change with the continued northward migration of Atlantic phytoplankton species, the northward retreat of the marginal ice zone and the observed decline in the DSi inventory of the subpolar North Atlantic Ocean over the last 3 decades.

1 Introduction

Diatoms are photosynthesising algae that take up dissolved silicic acid (DSi) from seawater to build silica-based frustules (termed “biogenic silica” (BSi) or “opal”), which are then recycled or reworked in transition to and within the seafloor. The seafloor acts as a biogeochemical reactor, generating a benthic return flux of DSi across the pan-Arctic region that is estimated to equal the input from all Arctic rivers (März et al., 2015). These recycling and reworking processes are therefore important for the regional silicon (Si) budget and for fuelling subsequent blooms, where seafloor-derived nutrients are able to be advected into the photic zone. Typically, Barents Sea phytoplankton spring blooms are dominated by diatoms (Wassmann et al., 1999; Orkney et al., 2020). However, temperate flagellate species are becoming more dominant in the Eurasian Basin of the Arctic Ocean and are expected to become the resident bloom formers in the region (Neukermans et al., 2018; Orkney et al., 2020; Oziel et al., 2020; Ingvaldsen et al., 2021). This shift in species composition is thought to be driven in part by an expansion of the Atlantic Water realm (“Atlantification”) (Fig. 1). Furthermore, nutrient concentrations in Atlantic Water flowing into the Barents Sea have declined over the last 3 decades and are forecast to do so throughout the 21st century (Neukermans et al., 2018, and references therein). Crucially, a much more significant drop in DSi concentrations has been observed relative to nitrate (Rey, 2012; Hátún et al., 2017), creating less favourable conditions for diatom growth (Neukermans et al., 2018).

This shift in phytoplankton community composition is predicted to reduce the export efficiency of phytodetritus, with potentially significant implications for pelagic–benthic coupling and thus the Si cycle (Fadeev et al., 2021; Wiedmann et al., 2020). Observations from long-term sediment trap data show that carbon export and aggregate sinking rates are 2-fold higher underneath diatom-rich blooms in seasonally sea-ice-covered areas of the Fram Strait, compared with that in Phaeocystis pouchetii-dominated blooms in the ice-free region (Fadeev et al., 2021). A similar contrast was observed in carbon export fluxes measured using short-term sediment trap deployments north of Svalbard (Dybwad et al., 2021). It is estimated that 40 %–96 % of surface ocean primary production is exported to the seafloor in the Barents Sea (Cochrane et al., 2009, and references therein), while the export efficiency of net primary production out of the euphotic zone in the central gyres is typically < 10 % (Turner, 2015, and references therein).

Figure 1 Map of the Barents Sea shelf, including ChAOS stations, the main hydrographic features and water masses from Ward et al. (2022) (PF – oceanic polar front, AW – Atlantic Water, ArW – Arctic Water, NCCW – Norwegian Coastal Current Water, BSW – Barents Sea Water, BSO – Barents Sea Opening, BSX – Barents Sea Exit). Dashed red path refers to a subsurface water mass (Lien et al., 2013). High-resolution bathymetry data taken from the General Bathymetric Chart of the Oceans (GEBCO) (Jakobsson et al., 2012).

Given the changes forecast in the pelagic–benthic coupling of Si in the Arctic, it is important to understand the baseline benthic biogeochemical system in order to anticipate the implications of further perturbations. Based on Si isotopic data from various reactive sedimentary pools and the sediment pore water dissolved phase from the Barents Sea seafloor, Ward et al. (2022) hypothesised that the Si cycle is isotopically coupled to the redox cycling of metal oxides, most likely solid-phase Fe (oxyhydr)oxides. The reductive dissolution of Fe (oxyhydr)oxides and release of adsorbed Si (FeSi) are thought to drive marked shifts in the isotopic composition of the Barents Sea sediment pore water DSi pool towards lower values. Further, Ward et al. (2022) propose that sediment pore water undersaturation drives the contemporaneous dissolution of lithogenic silicate minerals (LSi) alongside BSi, some of which are reprecipitated as authigenic clay minerals (AuSi), representing a sink of isotopically light Si to the regional Si budget. Finally, Ward et al. (2022) propose that seasonal pelagic phytoplankton blooms generate stark peaks in pore water DSi that dissipate on the order of weeks to months. However, to fully understand the early diagenetic cycling of Si within the seafloor of the Barents Sea, we must be able to quantify the relative contribution of LSi and BSi to the DSi pool, as well as establish whether AuSi precipitation removes a significant portion of that pool. Here we employ steady state reaction–transport modelling to reconstruct the benthic cycling of Si in the Barents Sea, informed by our dataset of solid- and dissolved-phase Si isotopic compositions (Ward et al., 2022) to test these hypotheses. Such techniques allow for the disentangling and quantification of the aforementioned early diagenetic reactions (Geilert et al., 2020a; Ehlert et al., 2016a; Cassarino et al., 2020), as well as the return benthic flux of DSi to the overlying bottom water. Furthermore, reaction–transport modelling allows for the quantification of processes on much shorter timescales; thus we use transient model runs to validate the hypothesis that the pulsed deposition of bloom-derived BSi can perturb the benthic Si cycle. We then quantify the bloom-derived BSi contribution to the total annual benthic DSi flux, the deposition of which is subject to the anticipated shifts in community compositions of pelagic primary producers across the Arctic Ocean.

Understanding the key aims presented here not only is important for anticipating the biogeochemical response of the Barents Sea seafloor to physical, chemical and biological changes in the surface ocean but also has implications for the pan-Arctic Si budget. Currently there are disparities in the isotopic and mass balances of the Arctic Ocean Si budget, with Torres-Valdés et al. (2013) concluding that the Arctic Ocean is a slight net exporter of Si. Furthermore, a recent isotopic assessment identified the need for an additional benthic sink of light Si to close the Si budget (Brzezinski et al., 2021). However, current understanding is limited by a lack of direct observations from major gateways, including the Barents Sea (Brzezinski et al., 2021). By coupling observational data with reaction–transport modelling, we are able to construct a balanced Si budget for the Barents Sea (Sect. 3.5), contributing to the data gaps that currently limit our understanding of pan-Arctic Ocean Si cycling.

2 Oceanographic setting, materials and methods

2.1 Oceanographic setting

The Barents Sea is one of seven shelf seas encircling the central Arctic Ocean and lies on the main inflow route for Atlantic Water. Oceanic circulation is driven by regional cyclonic atmospheric circulation and constrained by areas of prominent bathymetry (Fig. 1) (Smedsrud et al., 2013). Atlantic Water is fed in through the Barents Sea Opening between mainland Norway and Bear Island. This water mass then flows northwards, where it is met by colder, fresher Arctic Water, infiltrating the Barents Sea from the northern openings (Oziel et al., 2016). The oceanic polar front delineates these two water masses, the geographic position of which is tightly constrained in the western basin by the bathymetry but is less well defined in the east (Barton et al., 2018; Oziel et al., 2016). The heat content of the Atlantic Water-dominated region south of the polar front maintains a sea-ice-free state year-round, whereas the northern Arctic Water realm is seasonally sea-ice-covered, with a September minimum and a March/April maximum (Årthun et al., 2012; Faust et al., 2021).

The Barents Sea winter sea ice extent has been in decline since circa 1850 (Shapiro et al., 2003), but from 1998 the rate of retreat has become the most rapid observed on any Arctic shelf (Oziel et al., 2016; Årthun et al., 2012). Current forecasts suggest the Barents Sea will become the first year-round, sea-ice-free Arctic shelf by 2075 (± 28 years) (Onarheim and Årthun, 2017). The atmospheric and water column warming driving this sea ice retreat is a result of both anthropogenic and natural processes, with recent Atlantification arising from a northward expansion of the Atlantic Water realm (Årthun et al., 2012) and a reduction in sea ice import to the northern Barents Sea. The impact of these changes is an increase in upward heat fluxes, which inhibits sea ice formation (Lind et al., 2018).

The dynamic nature of the Barents Sea with respect to the physical oceanographic characteristics is reflected biologically in the ecosystems of the two main hydrographic realms. Annual primary production is estimated to range from 70 to 200 g C m⁻², with lower values found in the northern Arctic Water realm, where a deep meltwater-formed pycnocline limits nutrient replenishment through wind-induced mixing (Sakshaug, 1997; Wassmann et al., 1999). However, the most distinct peaks in the rates of primary production are found in the marginal ice zone (MIZ) (reaching 1.5–2.5 g C m⁻² d⁻¹; Hodal and Kristiansen, 2008; Titov, 1995) (Fig. 2), which forms in spring/early summer as sea ice melts and retreats northwards, stratifying the water column and stabilising the nutrient-rich photic zone (Wassmann et al., 2006; Reigstad et al., 2002; Olli et al., 2002; Krause et al., 2018; Wassmann et al., 1999; Vernet et al., 1998; Wassmann and Reigstad, 2011). The phytoplankton communities of the Barents Sea in proximity to the polar front and MIZ tend to be dominated by pelagic and ice-associated diatom species, as well as the prymnesiophyte P. pouchetii (Syvertsen, 1991; Wassmann et al., 1999; Degerlund and Eilertsen, 2010; Makarevich et al., 2022).

2.2 Reaction–transport modelling

2.2.1 General approach

We use the Biogeochemical Reaction Network Simulator (BRNS) to disentangle the interplay of chemical and physical processes involved in the early diagenetic cycling of Si at stations B13, B14 and B15 of the Changing Arctic Ocean Seafloor (ChAOS) project in the Barents Sea (Figs. 1 and 2, Table 1). These stations span the main hydrographic features (polar front) and realms (Atlantic Water and Arctic Water) of the Barents Sea. BRNS is an adaptive simulation environment suitable for large, mixed kinetic–equilibrium reaction networks (Regnier et al., 2002, 2003; Aguilera et al., 2005), which is based on a vertically resolved mass conservation equation (Eq. 1) (Boudreau, 1997), simulating concentration changes for solid and dissolved species (i) in porous media at each depth interval and time step.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\mathit{\delta }\mathit{\sigma }{C}_i}{\mathit{\delta }t}}& ={\displaystyle \frac{\mathit{\delta }}{\mathit{\delta }z}}\left(D_{\mathrm{bio}}\mathit{\sigma }{\displaystyle \frac{\mathit{\delta }{C}_i}{\mathit{\delta }z}}+D_i\mathit{\sigma }{\displaystyle \frac{\mathit{\delta }{C}_i}{\mathit{\delta }z}}\right)-{\displaystyle \frac{\mathit{\delta }\mathit{\sigma }\mathit{\omega }{C}_i}{\mathit{\delta }z}}\\ & +{\mathit{\alpha }}_i\mathit{\sigma }\left(C_i\left(\mathrm{0}\right)-C_i\right)+\sum _j{\mathit{\lambda }}_i^j{R}_j,\end{array}\end{array}$$

where ω, C_i, t and z represent the sedimentation rate, concentration of species i, time and depth respectively. The porosity term σ is given as $\mathit{\sigma }=(\mathrm{1}-\mathit{\phi })$ for solid species and σ=φ for dissolved species, where φ is sediment porosity (Table S2). This term ensures that the respective concentrations represent the amount or mass per unit volume of sediment pore water or solids as required (Boudreau, 1997). D_bio is the bioturbation coefficient (cm² yr⁻¹) and was determined experimentally alongside this study (Solan et al., 2020). D_i (cm² yr⁻¹) is the effective molecular diffusion coefficient (D_i=0 for solids), and α_i (yr⁻¹) represents the bioirrigation rate (α_i=0 for solids). R_j represents the rate of each reaction (j), and ${\mathit{\lambda }}_i^j$ is its stoichiometric coefficient. A full description of the model can be found in Sect. S2 of the Supplement.

Figure 2 Changing Arctic Ocean Seafloor (ChAOS) project summary. Chlorophyll α represents the peak value measured at each station during JR16006 CTD casts; data are available at https://www.bodc.ac.uk/data/published_data_library/catalogue/10.5285/89a3a6b8-7223-0b9c-e053-6c86abc0f15d/ (last access: 14 April 2022). Benthic nutrient flux magnitudes are for DSi measured in this study and Ward et al. (2022), as well as ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ from Freitas et al. (2020). The red box is a schematic summary of the main processes involved in the early diagenetic cycling of Si in the Barents Sea, derived from the results of the steady state model simulations in this study. Sea ice extent represents an approximation of the conditions at the time of sampling in 2017 and 2019. BSi reactivities were determined by the steady state model simulations. Please see Supplement Sect. S4 for a description of sediment pigment extraction methods. BSi – biogenic silica; LSi – lithogenic silica; AuSi – authigenic clay minerals; DSi – dissolved silicic acid.

2.2.2 Steady state reaction–transport modelling

Steady state modelling was employed to reproduce Si isotopic observational data in order to quantify the reaction rates of key processes involved in the cycling of Si within the seafloor. The version of the Si BRNS model employed here is adapted from Cassarino et al. (2020), which largely follows the approach of Ehlert et al. (2016a) and assumes a steady state. To ensure a steady state was achieved in the baseline simulations, the applied run time was dependent upon the sedimentation rate (0.05–0.06 cm yr⁻¹; Zaborska et al., 2008; Faust et al., 2020) and core length so as to allow for at least two full deposition cycles (∼ 2500 years for a 50 cm Barents Sea core). The implemented reaction network accounts for a pool of pore water DSi, sourced by a dissolving BSi phase, from which Si can be incorporated into authigenic clay minerals (AuSi) as they precipitate. The kinetic rate law for the dissolution of BSi follows Eq. (2) (Hurd, 1972), where k_diss is the reaction rate constant (yr⁻¹) and BSi_sol is the solubility of BSi (mol cm⁻³), implying that the rate of dissolution is proportional to the saturation state. The rate of BSi dissolution is allowed to decrease exponentially downcore in order to account for a reduction in reactivity due to BSi maturation and interaction with dissolved Al, as well as the preferential dissolution of more reactive material at shallower depths (Rickert, 2000; Van Cappellen and Qiu, 1997b; Rabouille et al., 1997; Dixit et al., 2001). The rate constants for BSi dissolution (k_diss, Eq. 2) were constrained using the solid-phase BSi content measurements (Fig. 3).

Equation (2) represents a simplification of the reaction rate law, which in reality is influenced by processes not incorporated into the model, such as surface area, temperature, pH, pressure and salinity. It is possible in some circumstances for the dissolution rate to deviate from the linear rate law (Van Cappellen et al., 2002); however, it is generally accepted that the dissolution of BSi is predominantly driven thermodynamically by the degree of undersaturation, leading to the linear rate law implemented in this study (Van Cappellen et al., 2002; Rimstidt and Barnes, 1980; Van Cappellen and Qiu, 1997b; Loucaides et al., 2012).

$$\begin{array}{}\text{(2)}& R_{\mathrm{d}b}=k_{\mathrm{diss}}\cdot \left[\mathrm{BSi}\right]\cdot \left(\mathrm{1.0}-{\displaystyle \frac{\left[\mathrm{DSi}\right]}{{\mathrm{BSi}}_{\mathrm{sol}}}}\right)\end{array}$$

The precipitation of AuSi was modelled through Eq. (3), where k_precip is the precipitation rate constant (Ehlert et al., 2016a). This rate law assumes that the reaction will proceed, providing the concentration of DSi is greater than the solubility of the AuSi (AuSi_sol). The rate is thus proportional to the degree of pore water DSi oversaturation (Ehlert et al., 2016a). We assume a value of 50 µM for AuSi_sol at all three stations (Lerman et al., 1975; Hurd, 1973) (Table S2). As with BSi dissolution, the rate of AuSi precipitation was allowed to decrease exponentially with depth, compatible with the hypothesis that the majority of AuSi precipitation occurs in the upper portion of marine sediment cores. Here, DSi can more easily precipitate in the presence of more readily available dissolved Al, the concentration of which is typically higher in the upper reaches of shelf sediments, sourced from the dissolution of reactive LSi (e.g. feldspar and gibbsite) contemporaneously to that of BSi (Aller, 2014; Rabouille et al., 1997; Van Beusekom et al., 1997; Ehlert et al., 2016a).

$$\begin{array}{}\text{(3)}& R_{\mathrm{p}}=k_{\mathrm{precip}}\cdot \left({\displaystyle \frac{\left[\mathrm{DSi}\right]}{{\mathrm{AuSi}}_{\mathrm{sol}}}}-\mathrm{1.0}\right)\text{if}\left[\mathrm{DSi}\right]>{\mathrm{AuSi}}_{\mathrm{sol}}\end{array}$$

In addition to the dissolution of BSi and precipitation of AuSi accounted for in previous early diagenetic modelling studies of the benthic Si cycle (Ehlert et al., 2016a; Cassarino et al., 2020), we incorporate the dissolution of LSi, which is thought to be an important oceanic source of numerous elements, including Si (Geilert et al., 2020a; Tréguer et al., 1995; Jeandel et al., 2011; Fabre et al., 2019; Ehlert et al., 2016b; Jeandel and Oelkers, 2015; Pickering et al., 2020; Morin et al., 2015). Here we assume that the dissolution of LSi is predominantly driven by the degree of undersaturation (Eq. 4), although as with the dissolution of BSi, LSi dissolution is a complex reaction and sensitive to processes that are not included in the model, including the potential for being catalysed by microbes (Vandevivere et al., 1994; Vorhies and Gaines, 2009; Liu et al., 2017). The undersaturation of Si minerals is known to include most primary and secondary silicates; thus dissolution extends beyond BSi in marine sediments (Isson and Planavsky, 2018, and references therein). Indeed, a suite of experiments have shown that primary silicates and clay minerals can rapidly release Si when placed in DSi-undersaturated seawater and take up Si in DSi-enriched waters (Siever, 1968; Mackenzie et al., 1967; Lerman et al., 1975; Hurd et al., 1979; Fanning and Schink, 1969; Mackenzie and Garrels, 1965; Gruber et al., 2019; Pickering, 2020). Lerman et al. (1975) determined in one such experiment that the dissolution of eight clay minerals could be described by a first-order reaction rate law driven by the saturation state, consistent with that applied here. Further, our assumption that LSi dissolution is driven by the degree of undersaturation is consistent with the suggestion that low bottom water DSi concentrations of the North Atlantic Ocean could allow for the dissolution of silicate minerals and thus account for high benthic DSi flux magnitudes in areas almost devoid of BSi (< 1 wt %) (Tréguer et al., 1995). A value of ∼ 100 µM was used for the solubility of LSi (LSi_sol) at all three stations, consistent with observations during multiple dissolution experiments of common silicate minerals in seawater (Table S3) as well as the estimated solubility of amorphous silica in high-detrital-component estuarine sediments (Kemp et al., 2021).

$$\begin{array}{}\text{(4)}& R_{\mathrm{d}l}=k_{\mathrm{LSidiss}}\cdot \left[\mathrm{LSi}\right]\cdot \left(\mathrm{1.0}-{\displaystyle \frac{\left[\mathrm{DSi}\right]}{{\mathrm{LSi}}_{\mathrm{sol}}}}\right)\end{array}$$

The desorption of Si from solid Fe (oxyhydr)oxide phases under anoxic conditions was simulated using a simple reaction rate constant (k_FeSi), representing the rate of desorption. The value assigned to k_FeSi was calculated during the modelling exercise, and no assumed amount of FeSi was included in the upper boundary conditions. This parameter likely represents a significant simplification; however the exact process pertaining to the adsorption of Si onto Fe (oxyhydr)oxides is unclear and requires further study (Geilert et al., 2020a). Step functions were included in the FeSi reactions in the model to simulate the desorption of this phase at specific depth intervals, representing the Fe redox boundaries identified in Ward et al. (2022). The step functions act as a cut-off mechanism, either setting reaction rates to zero or activating them at specific depths. A full description of the model, including all boundary conditions and how isotopic fractionation was imposed in the AuSi precipitation and FeSi desorption reactions, can be found in Sect. S2 of the Supplement.

Our estimates for all reaction rate constants in the steady state simulations (k_diss, k_precip, k_LSidiss and k_FeSi) were not based on published values and were model-derived (Table S2). These values were constrained by ensuring the best fit of the observational data with the simulated solid-phase BSi content and pore water DSi concentration and isotopic compositions, which were obtained by minimising the root-mean-square error (RMSE) between simulated and measured values (Table S1). Despite being model-derived, k_diss values (0.0055–0.074 yr⁻¹) are found to lie within the published range for marine sediment BSi (Sect. 3.2). After the best-fit scenarios were established for each station, a sensitivity experiment was carried out by sequentially setting each reaction rate constant to zero in order to assess the importance of each process to the model fit (Fig. 3).

2.2.3 Processing the simulated data

Depth-integrated rates (R) of a given reaction (j) were calculated across the model domain using Eq. (5), for the best-fit simulation data of each station.

$$\begin{array}{}\text{(5)}& R_j\left(x\right)=\sum _j\underset{x=\mathrm{0}}{\overset{x=L}{\int }}{R}_j\mathrm{d}x,\end{array}$$

where L is the model domain length and dx denotes the given depth interval.

The deposition flux of BSi at the sediment–water interface (SWI) (J_BSi,in) was then calculated based on Eq. (6), which states that J_BSi,in equates to the sum of the flux of BSi out of the sediment (assumed to equate to the integrated rate of BSi dissolution, R_db) and the BSi burial flux (J_BSi,bur) (Burdige, 2006; Freitas et al., 2021). J_BSi,bur was estimated at the base of the model domain (50.4 cm), following Eq. (7) on mass accumulation (Varkouhi and Wells, 2020, and references therein), which is controlled solely by advection. The sedimentation rate at depth (ω_z) was corrected for compaction following Eq. (8) (Berner, 1980). J_BSi,bur calculations assume a sediment wet bulk density of 1.7 g cm⁻³, consistent with previously assumed values for the Arctic seabed (Brzezinski et al., 2021; Backman et al., 2004) and that measured in clay-rich sediments of the Barents Sea (Orsi and Dunn, 1991).

$$\begin{array}{}\text{(6)}& {\displaystyle }{J}_{\mathrm{BSi},\mathrm{in}}=R_{\mathrm{d}b}+J_{\mathrm{BSi},\mathrm{bur}}\text{(7)}& {\displaystyle }{J}_{\mathrm{BSi},\mathrm{bur}}=(\mathrm{1}-{\mathit{\phi }}_z)\cdot {\mathit{\omega }}_z\cdot [\mathrm{BSi}{]}_z\cdot {\mathit{\rho }}_z\text{(8)}& {\displaystyle }{\mathit{\omega }}_z={\mathit{\omega }}_{\mathrm{0}}\cdot (\mathrm{1}-{\mathit{\phi }}_{\mathrm{0}})/(\mathrm{1}-{\mathit{\phi }}_z)\end{array}$$

We are also able to use the model simulation output to determine the total benthic flux of DSi at the SWI (J_tot), which has multiple constituent parts that contribute to the benthic flux magnitude. Following Eq. (9) (Freitas et al., 2020), we calculate J_tot and thus the relative contributions from bioturbation (J_bioturb), bioirrigation (J_bioirr), advection (J_adv) and molecular diffusion (J_diff) to complement the calculated J_diff estimates and core-incubation-derived J_tot of Ward et al. (2022).

$$\begin{array}{}\text{(9)}& J_{\mathrm{tot}}=J_{\mathrm{bioturb}}+J_{\mathrm{bioirr}}+J_{\mathrm{adv}}+J_{\mathrm{diff}}\end{array}$$

2.2.4 Transient reaction–transport modelling

The influence of seasonality in pelagic primary production on the benthic Si cycle of the Barents Sea has been inferred through interpretation of pore water DSi depth profiles at station B14 (Ward et al., 2022). However, BRNS assumes a steady state and therefore cannot resolve seasonal biogeochemical dynamics without modification to allow certain boundary conditions to become time-dependent, enabling their activation and deactivation on a temporal scale. Here we use transient model runs to test the hypothesis that the pulsed deposition of bloom-derived BSi can rapidly perturb the benthic DSi pool, which is then able to recover on the order of weeks to months.

The steady state baseline simulations at station B14 represent a data–model best fit of the 2018 observational data, wherein we do not observe transient peaks in sediment pore water DSi concentrations (Fig. 4). This steady state scenario was used as the initial conditions for the transient simulations, which were run for 1 simulation year, producing output data at weekly time intervals. We simulate the phytoplankton spring bloom event by incorporating a step function to temporarily increase the BSi deposition flux and reactivity, simulating the effects of a 1- to 3-week spring bloom in late spring/early summer on the delivery of BSi to the seafloor. Durations of 1 and 3 weeks are thought to represent the typical length of a Barents Sea MIZ and spring bloom respectively (Sakshaug, 1997; Dalpadado et al., 2020) (Fig. 4). All boundary conditions were kept constant, with the exception of those related to the bloom-derived BSi pool.

During the 1- to 3-week time interval, the bloom-derived material was deposited at a rate equivalent to an increase in the steady state BSi deposition flux of 10- to 30-fold (Fig. 4). The background BSi deposition flux magnitude at station B14 is a model-derived parameter, constrained in the steady state simulations with the measured sediment BSi content and pore water DSi depth profiles. A 30-fold increase on the background deposition flux is similar to observations of a 10 d post-bloom diatom mass sinking event in the subpolar North Atlantic down to 750 m depth (26-fold increase) (Rynearson et al., 2013). The Barents Sea covers a relatively shallow continental shelf (average depth 230 m), and intense physical mixing at the polar front has been shown to enhance rates of vertical organic carbon flux at depth, close to station B14 (Wassmann and Olli, 2004). We could therefore anticipate an even greater increase in BSi depositional flux under bloom conditions than the maximum value assumed here. The reactivity of the bloom-derived material (k_dissbloom) ranged from 5 to 20 yr⁻¹, which is within the reactivity range of fresh pelagic BSi (3 to 100 yr⁻¹; Ragueneau et al., 2000; Nelson and Brzezinski, 1997) (Fig. 4). Each of these three boundary conditions (length of the bloom period, k_dissbloom and the deposition flux) was varied across multiple simulations within the constraints of published values to assess the influence of each parameter on the size and longevity of the sediment pore water DSi peak. The boundary conditions were determined to test whether peaks in pore water DSi concentrations are able to form within 1.5 months and dissipate within 3 months of the bloom, as proposed by Ward et al. (2022) and evidenced by contrasting sea ice conditions across the three cruises (Sect. 3.2; Fig. 4).

2.3 Observational data

As many reactions responsible for the biogeochemical cycling of Si between the solid and dissolved phases fractionate the isotopes of Si (²⁸Si, ²⁹Si, ³⁰Si) relative to each other, we are able to use stable Si isotopes as a tool to trace these pathways. All solid-phase, core top and sediment pore water samples collected for Si isotopic analysis were collected over three summers in the Barents Sea (30^∘ E transect spanning 74 to 81^∘ N) aboard the RRS James Clark Ross (2017, 2018, 2019). Dissolved-phase pore and core top water DSi concentration measurements were determined on board using a Lachat QuikChem 8500 flow injection analyser with an accuracy of 2.8 %, defined using certified reference materials (CRMs; Kanso Technos Co., Ltd.). Stable Si isotopic compositions of the samples were determined at the University of Bristol in the Bristol Isotope Group laboratory. Isotopic compositions are expressed in δ³⁰Si notation (per mille, ‰), relative to the international Si standard NBS-28 (Eq. 10). A full description of field methods, as well as Si isotopic and concentration data of the solid- and dissolved-phase reconstructed using BRNS, is provided in Ward et al. (2022).

$$\begin{array}{}\text{(10)}& {\mathit{\delta }}^{\mathrm{30}}\mathrm{Si}=\left({\displaystyle \frac{{\left({}^{\mathrm{30}}\mathrm{Si}/{}^{\mathrm{28}}\mathrm{Si}\right)}_{\mathrm{sample}}}{{\left({}^{\mathrm{30}}\mathrm{Si}/{}^{\mathrm{28}}\mathrm{Si}\right)}_{\text{NBS-28}}}}-\mathrm{1}\right)\cdot \mathrm{1000}\end{array}$$

Model upper boundary conditions for the dissolved phase were based on the DSi concentration and Si isotopic composition of the core top waters (+1.64 ± 0.19 (n=5), +1.46 ± 0.15 (n=3) and +1.69 ± 0.18 ‰ (n=6) at stations B13, B14 and B15 respectively). The Si isotopic compositions of the solid phases (BSi, LSi, FeSi) implemented in the model were determined based on a series of sequential digestion experiments carried out on surface sediment samples of the three aforementioned stations (Ward et al., 2022). The digestion sequence (0.1 M HCl, 0.1 M Na₂CO₃, 4 M NaOH) activates operationally defined reactive pools of Si, including Si-HCl (Si associated with metal oxide coatings on BSi, assumed here to represent the FeSi pool), Si-Alk (authigenically altered and unaltered BSi) and Si-NaOH (LSi and refractory BSi) (Pickering et al., 2020). The composition of the FeSi (−2.88 ± 0.17 ‰, n=20) and LSi (−0.89 ± 0.16 ‰, n=18) pools in Barents Sea sediment leachates was within long-term reproducibility of diatomite standard measurements (2σ ± 0.14 ‰, n=116) across the three stations, whereas the BSi phase was found to vary spatially, with an isotopic composition of +0.82 ± 0.16 ‰ (n=14) at station B15 and +1.43 ± 0.14 ‰ (n=8) and +1.50 ± 0.19 ‰ (n=7) at stations B13 and B14 respectively (Ward et al., 2022).

Table 1 Sampling station information averaged across the three Changing Arctic Ocean Seafloor (ChAOS) cruises.

3 Results and discussion

Gaining a better mechanistic understanding of benthic biogeochemical Si cycling in the Barents Sea is important to anticipate the effects of further natural and anthropogenically driven environmental perturbations, as well as to begin to address key knowledge gaps that currently limit our understanding of the pan-Arctic Ocean Si budget (Brzezinski et al., 2021). Ward et al. (2022) propose that LSi is dissolving alongside BSi in Barents Sea sediments, that part of the DSi pool is reprecipitated as AuSi, that the benthic Si and Fe redox cycles are coupled, and that seasonal variability in the deposition of BSi is expressed in the sediment pore water DSi pool. Here we reconstruct the benthic Si cycle of the Barents Sea using a reaction–transport model to further investigate and disentangle the interplay of processes that combined to produce our observational dataset and test these hypotheses. This approach allows for the quantification of reaction rates and fractionation factors, as well as of deposition and benthic flux magnitudes, which are used to inform a balanced Barents Sea Si budget that could have implications for the pan-Arctic Ocean budget.

3.1 What can reaction–transport modelling reveal about the controls on the background, steady state benthic Si cycle?

Coupling a reaction–transport model with δ³⁰Si values measured in the solid and dissolved phases represents a powerful tool with which to trace early diagenetic reactions, given that multiple reaction pathways have been shown to fractionate Si isotopes. One early diagenetic process that fractionates isotopes of Si is the formation of clay minerals, which preferentially takes up the lighter isotope from the dissolved phase, leaving residual waters relatively isotopically heavy in composition. A Si isotopic fractionation factor (³⁰ϵ) associated with AuSi formation has yet to be thoroughly established; however Ehlert et al. (2016a) and Geilert et al. (2020a) modelled a ³⁰ϵ of −2 ‰ for marine AuSi formation. A ³⁰ϵ of this magnitude is consistent with the formation of clay minerals in riverine and terrestrial settings (−1.8 ‰ to −2.2 ‰) (Hughes et al., 2013; Ziegler et al., 2005a, b; Opfergelt and Delmelle, 2012) and similar to that observed in the adsorption of Si onto Fe (oxyhydr)oxide minerals (³⁰ϵ of −0.7 ‰ to −1.6 ‰) (Zheng et al., 2016; Delstanche et al., 2009; Wang et al., 2019; Opfergelt et al., 2009). However, the magnitude of ³⁰ϵ associated with AuSi precipitation can reach up to −3 ‰ in deep-sea settings (Geilert et al., 2020b), likely depending on pore water properties (pH, temperature, salinity, saturation states). This relatively high ³⁰ϵ is also consistent with repetitive clay mineral dissolution–reprecipitation cycles (Opfergelt and Delmelle, 2012). Similarly, ³⁰ϵ during Si adsorption onto Fe (oxyhydr)oxides is thought to increase with mineral crystallinity (³⁰ϵ of −1.06 ‰ for ferrihydrite and −1.59 ‰ for goethite) (Delstanche et al., 2009). Isotopic fractionation during dissolution is less well constrained, with previous work suggesting that BSi dissolution could induce a slight fractionation that enriches the DSi pool in the lighter isotope (³⁰ϵ of −0.55 ‰ to −0.86 ‰; Demarest et al., 2009; Sun et al., 2014) or occur without isotopic fractionation (Wetzel et al., 2014; Egan et al., 2012). Here we assume the latter and thus impose a ³⁰ϵ of 0 ‰ for the dissolution of BSi and LSi, which is consistent with similar, previous reaction–transport model studies (Geilert et al., 2020a; Ehlert et al., 2016a).

3.1.1 LSi dissolution and AuSi precipitation

Model results show that the sediment pore water DSi profiles cannot be reproduced by the dissolution of BSi alone at all three stations (B13, B14 and B15) (dash-dotted blue lines, Fig. 3). At station B13 the simulations suggest that while there is sufficient DSi released to reproduce the asymptotic DSi concentration due to the higher BSi content at depth compared to B15, the rate of release in the upper sediment layers is not consistent with that in the pore water DSi concentration profiles downcore of the SWI in the observational data (Fig. 3). This observation is in contrast to B15, where the simulated asymptotic DSi concentration is just 23 µM when BSi is the only source of DSi, consistent with the measured BSi content profiles that suggest a cessation in BSi dissolution by the middle of the sediment cores (∼ 15 cm, asymptotic BSi content of ∼ 0.2 wt %) (Fig. 3). Because of the continued increase in DSi with depth at station B13 (solid grey line), partly driven by the elevated BSi content in the mid-core relative to B15, a relatively slow rate of AuSi precipitation is required at depth to take up the excess DSi and reproduce the observed asymptotic DSi value. Generally it is assumed that AuSi precipitation is concentrated in near-SWI sediments (0–5 cm), where the concentration of other essential solutes (Al, Fe, Mg²⁺, K⁺, Li⁺, F⁻) is generally highest, sourced from Fe and Al (oxyhydr)oxides and reactive LSi (Ehlert et al., 2016a; Van Cappellen and Qiu, 1997a; Mackin and Aller, 1984; Aller, 2014). However, the uptake of DSi through AuSi precipitation has previously been inferred in terrigenous-dominated shelf sediments of the Arctic Ocean at > 50 cm depth (März et al., 2015).

Figure 3 Dashed red lines show best-fit steady state model simulations for stations B13 (a, b, c) and B15 (d, e, f) (see Fig. S2 for station B14), which require BSi and LSi dissolution, AuSi precipitation, and FeSi desorption (see Table S2 for all boundary conditions). Additional simulations are based on a series of sensitivity experiments carried out to assess the importance of each reaction pathway. Left to right: sediment pore water DSi concentration, sediment pore water DSi δ³⁰Si, solid-phase BSi content. Vertical dashed black lines represent the core top water δ³⁰Si from 2017; open shapes show observational pore water data (Ward et al., 2022). Error bars on δ³⁰Si data are based on long-term reproducibility, derived from repeat measurements of the diatomite standard (2σ ± 0.14, n=116), unless the error derived from analytical replicates was greater.

Due to the discrepancies between observational and simulated DSi pore water concentration data, we incorporated a LSi phase into our model, which dissolves according to the presumed degree of undersaturation. Without this additional phase (when k_LSidiss is set to zero), model simulations show that the rate of DSi release is insufficient to reconstruct the observational DSi data (Fig. 3). Implementing LSi dissolution in conjunction with BSi produced the best data–model fit (dashed red lines). The dissolution of LSi has been inferred in a similar study of marine sediment pore waters of the Guaymas Basin (Geilert et al., 2020a), as well as in beach and ocean margin sediments (Fabre et al., 2019; Ehlert et al., 2016b). Indeed, Morin et al. (2015) report Si dissolution rates for basaltic glass particles in seawater that exceed those of diatoms (Pickering et al., 2020, and references therein).

The inference based on the pore water DSi concentration profiles that an additional phase, most likely LSi, is dissolving into Barents Sea pore waters is supported by the simulated isotopic composition of the pore water phase at stations B13 and B15. Without the dissolution of the LSi phase, the δ³⁰Si of the pore waters represents a mixture of the composition of the core top water and the BSi phase composition (+1.44 ‰ and +1.04 ‰ at B13 (solid grey line) and B15 (dash-dotted blue line) respectively) (Fig. 3). In this simulation scenario at station B15, the integrated rate of AuSi precipitation is zero as the concentration of pore water DSi has not surpassed the imposed AuSi solubility so cannot influence the sediment pore water δ³⁰Si. Therefore, the model set-up without LSi dissolution cannot reproduce the intricacies of the downcore isotopic profile, as the lack of DSi released results in an insufficient concentration to allow for the precipitation of AuSi; thus the relative shift from isotopically lighter to heavier compositions between 0.5 and 2.5 cm cannot be resolved. The downcore shift to heavier isotopic compositions between 0.5 and 2.5 cm is thought to be caused by AuSi formation as the pore water DSi concentration crosses the saturation of the AuSi phase, facilitating its precipitation. Autochthonous benthic diatoms have been found in the Barents Sea at a maximum depth of 245 m (Druzhkova et al., 2018), which could fractionate the DSi pool during uptake with a ³⁰ϵ of ∼ −1.1 ‰ (De La Rocha et al., 1997). However, they are in very low abundance and so unlikely to contribute significantly to sediment pore water DSi uptake or isotopic fractionation in the uppermost sediment layers.

Figure 3 indicates that the observational data from stations B13 and B15 can be reproduced by a model that assumes a steady state, which is less the case for the data from station B14, particularly for the 2017 and 2019 profiles (Fig. S2). This observation suggests that the reaction–transport model will need to resolve transient, non-steady state dynamics in order to better represent the observational data, which will be discussed further in Sect. 3.2.

These model observations are consistent with a mass balance calculation using the isotopic compositions of the 0.5 cm pore water sample, as well as the BSi and LSi leachate samples at stations B13 and B14, which indicate a contemporaneous release of both phases. Our model findings therefore support the hypothesis of Ward et al. (2022) that LSi minerals are likely to be dissolving in the upper few centimetres of the Barents Sea seafloor. The depth-integrated reaction rates of the best-fit steady state simulations suggest that between 60 % and 98 % of the DSi released into the sediment pore water from the solid phase is sourced from the dissolution of LSi (Table 2). This range was determined by calculating the depth-integrated rate of LSi dissolution across the model domain and three stations as a proportion of the total integrated rate of DSi input from the three simulated sources (BSi dissolution, LSi dissolution and desorption from metal oxides). The predominance of LSi over BSi dissolution is consistent with the observation that Barents Sea sediments consists of ∼ 96 % terrigenous material (Ward et al., 2022), which is compatible with previous work showing that clay mineral assemblages in the Barents Sea are dominated by terrestrial signals from Svalbard and northern Scandinavia (Vogt and Knies, 2009). The Si isotopic composition of the LSi phase in surface sediments of stations B13, B14 and B15 (−0.89 ± 0.16 ‰; Ward et al., 2022) is also closer to secondary clay minerals (−2.95 ‰ to −0.16 ‰; Opfergelt and Delmelle, 2012, and references therein) than primary silicates of the crust and mantle (∼ 0 ‰ and −0.34 ‰ respectively; Opfergelt and Delmelle, 2012, and references therein). This isotopic composition is not surprising, given the predominance of the clay and silt size fraction in these sediment cores (87 % < 63 µm) (Faust et al., 2020). However, δ³⁰Si measured in the Si-NaOH pool could also represent a combination of reactive primary and secondary silicates, ranging from ∼ 0 ‰ to −2.95 ‰.

While previous mass balance calculations and the model-derived integrated reaction rates presented here agree that both LSi and BSi dissolution contribute to the sediment pore water DSi pool, the magnitudes of the LSi dissolution contribution vary significantly. Ward et al. (2022) suggest that just 14 % and 13 % of the DSi pool in the 0.5 cm pore water interval are sourced from the dissolution of LSi at stations B13 and B14 respectively, while at station B15 it is inferred that the δ³⁰Si can be resolved through BSi dissolution alone. This is in contrast to the 84 %, 60 % and 98 % contributions calculated here (Table 2). There are multiple contributing factors to this discrepancy. Firstly, previous mass balance calculations are based on one depth interval, whereas the estimates presented here are derived from depth-integrated dissolution rates of the entire 50.4 cm model domain. Furthermore, reaction–transport modelling has revealed that this contradiction is likely born of the assumption that the BSi pool at all three stations is sufficient to fuel the pore water DSi stock. Dissolution dynamics were not taken into account in the simple mass balance calculation of Ward et al. (2022); however here we have shown that the composition of Barents Sea surface sediments, which are almost devoid of BSi (0.26 wt %–0.52 wt %, or 92–185 µmol g⁻¹ dry wt), cannot reproduce the rate of pore water DSi build-up with depth from the SWI and can only support an asymptotic DSi concentration of 23 µM at station B15 (dash-dotted blue line, Fig. 3). This observation implies that the additional assumption of Ward et al. (2022) that the 0.5 cm pore water δ³⁰Si value is not impacted by AuSi precipitation could be invalid. If the 0.5 cm pore water interval was directly or indirectly influenced by AuSi precipitation, such an assumption would lead to an underestimation of the LSi contribution as the δ³⁰Si value would be isotopically heavier than if it were derived solely from dissolving solid phases mixing with trapped core top water.

Previously it has been suggested that the quantification of AuSi precipitation rates in marine sediments is not critical in order to fully understand the early diagenetic cycling of Si as reverse weathering typically represents a diagenetic solid-phase conversion from BSi to AuSi, via the dissolved phase (DeMaster, 2019). Model simulations reveal that 37 %, 2.9 % and 13.8 % of the DSi released across the 50 cm model domain at stations B13, B14 and B15 respectively are taken out of solution in the formation of AuSi. This observation is consistent with a similar, previous study of the Peruvian margin upwelling region, which determined that 24 % of DSi released from the dissolution of BSi was reprecipitated as AuSi at a rate of 1.53 mmol Si m⁻² d⁻¹ (Ehlert et al., 2016a). AuSi precipitation rates in sediments of the Amazon Delta on the other hand can reach 7.7 mmol Si m⁻² d⁻¹ (Michalopoulos and Aller, 2004). Reprecipitation of the DSi pool within the relatively shallow cores studied here will inhibit its exchange with overlying bottom waters. Therefore, in this context, AuSi precipitation can be considered a sink term for the regional Si budget (Ward et al., 2022). In a setting where AuSi precipitation occurs at such a depth where pore water DSi exchange with bottom waters is not possible, this reaction pathway could instead be considered an early diagenetic solid-phase conversion (from BSi and/or LSi to AuSi), as opposed to a true sink term (Frings et al., 2016; DeMaster, 2019).

Table 2 Depth-integrated reaction rates of the steady state, best-fit simulations across the upper 50 cm of sediment (Fig. 3), as well as calculated BSi deposition rates and model-derived DSi benthic fluxes. All magnitudes are given in units of mmol Si m⁻² d⁻¹.

3.1.2 Evidence for coupling of the benthic Fe and Si cycles in the Barents Sea

Ward et al. (2022) suggest that the benthic Fe and Si cycles are coupled in the Barents Sea, evidenced by a contemporaneous increase in pore water Fe concentrations with an enrichment in the lighter Si isotope of the DSi pool at all three stations. Model simulations support this hypothesis by demonstrating that the Barents Sea DSi pore water profiles can be reconstructed when applying the dissolution of both a BSi and a LSi phase; however under the model scenario where the desorption of FeSi is inhibited (k_FeSi=0), the δ³⁰Si pore water profiles are inconsistent with the observational data (solid black lines, Fig. 3). With the dissolution of the LSi phase implemented at both stations B13 and B15, it is possible to resolve the δ³⁰Si pore water profiles in the upper 2.5 and 8.5 cm respectively. However, below these depths the simulated profiles have isotopically heavier compositions than the observational data. Release of an isotopically light phase at specific depth intervals (beginning at 1.5 cm at B13 and 10 cm at B15) results in a simulated δ³⁰Si profile within the range of the observational data (Fig. 3). This isotopic shift to lower pore water δ³⁰Si is interpreted to represent the desorption of Si from solid Fe (oxyhydr)oxide phases, wherein the depth intervals of the same isotopic shifts correspond to the depths at which Fe is released into the pore waters (Faust et al., 2021; Ward et al., 2022). This increase in pore water Fe also occurs at similar depth intervals to decreases in pore water dissolved O₂ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations, consistent with a transition from oxic to anoxic conditions (Freitas et al., 2020; Ward et al., 2022).

As discussed above, the Si-HCl reactive Si pool is isotopically light and thought to be associated with metal oxide coatings on BSi (Pickering et al., 2020). The δ³⁰Si composition of the Si-HCl phase is assumed here to represent the composition of the phase desorbing across the Fe redox boundaries. When using a composition of −2.88 ‰ for the FeSi phase, simulation scenarios have identical DSi concentration profiles whether k_FeSi is active or set to zero (dashed red and solid black lines, Fig. 3). This similarity is in contrast to the δ³⁰Si pore water profiles, which cannot be reproduced without release of this isotopically light phase at the redox boundaries of all three modelled stations.

This discrepancy between the sediment pore water DSi and δ³⁰Si profiles when the FeSi desorption is active and inactive may explain why the influence of FeSi desorption is so apparent in Barents Sea sediment cores and more ambiguous in similar, previous studies of the Si cycle in lower-latitude marine sediments. The preferential adsorption of ²⁸Si onto Fe (oxyhydr)oxides and the subsequent dissolution or formation of these minerals have been used to interpret both heavy and light δ³⁰Si marine sediment pore water signals in previous work (Ehlert et al., 2016a; Geilert et al., 2020a). In addition, while not inducing a clear signal in the δ³⁰Si, redox cycling of Fe was highlighted as a potential regulating factor in the release of DSi into pore waters of the Greenland margin (Ng et al., 2020). Ng et al. (2020) hypothesised that the reductive dissolution of Fe mineral coatings increased the reactivity of the BSi pool, hence the elevated DSi concentrations found in cores with increased pore water Fe (Ng et al., 2020). In the Barents Sea, FeSi desorption across sedimentary redox boundaries is thought to be so prominent in the δ³⁰Si data because the asymptotic concentration of pore water DSi (∼ 100 µM) is much lower than that in the aforementioned studies (∼ 350–900 µM). The low sediment pore water DSi concentration allows for the direct detection of this process, whereas in previous studies the influence of Fe on the benthic Si cycle is inferred either through elevated DSi and Fe concentrations (Ng et al., 2020) or by depositional context, for example in cores sampled from systems with an abundance of reactive Fe (e.g. hydrothermal vent systems (Fe sulfides) – Geilert et al., 2020a, or the Peruvian oxygen minimum zone – Ehlert et al., 2016a).

At stations B13 and B14, as well as to a lesser extent at B15, there is an increase in the sediment pore water DSi concentration downcore of the middle to the base of the profiles. This feature can be reproduced in the model simulations with the desorption of FeSi from the respective redox boundary depths (Fig. S1); however the required ratio of k_FeSi for ²⁸Si and ³⁰Si suggests an isotopic composition of the FeSi phase of just −1.0 ‰ to −1.5 ‰. This isotopic composition is heavier than that measured in the Si-HCl pool (−2.88 ‰), likely reflecting a complexity in the desorption process not captured by the model. Nevertheless, both scenarios support the release of an isotopically light phase at depth, most likely sourced from the Fe redox cycle.

In summary, model simulations somewhat support the hypothesis, based on observational data, that the Barents Sea benthic Si cycle is influenced by the Fe redox system. Model results suggest that the influence of the Fe redox cycle is relatively unimportant for the magnitude of the pore water DSi pool, which appears to be controlled by release from the BSi and LSi phases. However, the coupling of these element cycles is evidenced in the pore water δ³⁰Si data, indicating that the Fe cycle is important for the isotopic budget within the seafloor. We suggest the influence of FeSi desorption is detectable in the Barents Sea pore water δ³⁰Si data, due to the relatively low pore water DSi concentrations and the distinctly isotopically light nature of the FeSi phase. These findings indicate that FeSi desorption should be considered when interpreting downcore δ³⁰Si trends, especially in low-DSi-concentration settings.

3.2 What can transient reaction–transport modelling reveal about non-steady state dynamics in the benthic Si cycle?

Observational pore water DSi concentration data from station B14 suggest that a pulsed increase in the deposition of reactive phytodetritus to the seafloor, derived from phytoplankton blooms, can drive transient peaks in pore water DSi of up to ∼ 300 µM (Ward et al., 2022). This non-steady state dynamic is evidenced in the 3 consecutive years of pore water DSi concentration data collected during the summers of 2017 to 2019, which show that in 2017 and 2019, when the MIZ was above station B14 just 1.5 months prior to sediment coring, a sediment pore water DSi peak is present. This characteristic is in contrast to 2018 when station B14 had been sea-ice-free for 3 months prior to sampling (Downes et al., 2021) and no peak in DSi concentration was observed in the pore water nutrient data (Fig. S6). This observation may indicate that the sediment pore water DSi peaks form under the MIZ, which supports the formation of phytoplankton blooms in late spring/early summer, supplying fresher BSi to the benthos relative to the background BSi pool. Dissolution of this fresher BSi could then fuel subsurface peaks in pore water DSi concentration that dissipate between 1.5 and 3 months after formation, driven by the increased concentration gradient across the SWI that would enhance the rate of molecular diffusion. Additional model simulations were carried out on the baseline, steady state, best-fit model scenario of station B14 to assess this hypothesis.

Results of the transient simulations show that it is possible for the deposition of fresh, bloom-derived BSi to reproduce the observed peaks in pore water DSi concentration in 2017 and 2019 (Fig. 4). Calculated rates of background BSi deposition across the three stations (0.17 to 1.71 mmol Si m⁻² d⁻¹, Table 2) are similar to BSi export rates measured in short-term and moored sediment traps in Kongsfjorden, Svalbard at 100 m depth and the eastern Fram Strait at 180–280 m depth (0.2–1.3 mmol Si m⁻² d⁻¹) (Lalande et al., 2013, 2016). A simulated 3-week, 10-fold increase in this BSi depositional flux at a much higher reactivity (k_dissbloom 20 yr⁻¹) than the background value (k_diss 0.074 yr⁻¹) derived from the 2018 baseline simulation, results in a DSi peak consistent with the magnitude of that in the observational data after 1.5 months (dashed red lines, Fig. 4). This simulated peak in DSi concentration is then able to dissipate by 3 months after bloom initiation. Similarly, with a shorter bloom (1 week), compatible with the typical length of an ice edge bloom in the Barents Sea (Dalpadado et al., 2020) and a 30-fold increase in BSi depositional flux of k_dissbloom 15 yr⁻¹, the DSi peak is able to form and dissipate on a time frame similar to that of the former scenario (solid black lines, Fig. 4). The generated peak in DSi concentration must be able to disperse after 3 months if the timing of core sampling relative to MIZ retreat is valid as the explanation for the lack of pore water DSi peak observed in the 2018 data.

Figure 4 Results of the transient simulations under different conditions. Multiplier refers to the increase in the deposition flux (J_BSi,in) of the bloom-derived BSi, which was assigned a reactivity (k_dissbloom in yr⁻¹) and deposited over either 1 or 3 weeks. Note that the “Bloom Initiated” panel represents the first day of the simulated bloom; therefore all the modelled scenarios overlap and reflect the background steady state simulation. The panels “+1.5 months” and “+3 months” refer to the time elapsed since the bloom was initiated, which are presented to demonstrate how the DSi peak evolves over time. Ward et al. (2022) suggested that DSi peaks could form within 1.5 months and dissipate after 3 months, as evidenced by contrasting sea ice conditions across the three cruises.

The implemented k_dissbloom values (5–20 yr⁻¹) are higher than those typically observed in marine sediment cores. k_diss values of 0.98–1.38 yr⁻¹ have been observed at 1 to 20 cm depth in sediment cores collected from the Porcupine Abyssal Plain (mean water depth of 4850 m), considered high for BSi in deep marine sediments (Ragueneau et al., 2001), while values of up to 6.8 yr⁻¹ have been measured in much shallower sediments of Jiaozhou Bay in the Yellow Sea (Wu et al., 2015). k_diss values are highly sensitive to species assemblage and temperature, generally decreasing from the surface ocean towards the seafloor (Natori et al., 2006; Rickert et al., 2002; Rickert, 2000; Roubeix et al., 2008; Kamatani, 1982; Tréguer et al., 1989; Ragueneau et al., 2000). However, k_diss values in fresh diatoms are thought to range from 3 to ∼ 100 yr⁻¹ (Ragueneau et al., 2000; Nelson and Brzezinski, 1997), and previous experiments have measured dissolution rate constants of 27 yr⁻¹ and above (Rickert, 2000; Kamatani and Riley, 1979) at lower temperatures (2–4 ^∘C) in seawater, consistent with the values employed in this study. Station B14 is located beneath the polar front, which is a location of intense physical mixing due to the interleaving of multiple water masses (Barton et al., 2018), which has been shown to drive enhanced depositional fluxes of particulate organic carbon relative to that measured both north and south of the frontal zone (Wassmann and Olli, 2004). Furthermore, data collected from sediment traps deployed to the north and north-west of Svalbard have uncovered an approximately 2-fold-higher vertical carbon export flux from diatomaceous aggregates formed in seasonally sea-ice-covered regions, compared with aggregates from P. pouchetii blooms (Fadeev et al., 2021; Dybwad et al., 2021). Therefore, it is considered here that given the shallow depth of the Barents Sea seafloor and the location of station B14 beneath the polar front, fresh and reactive BSi formed in MIZ blooms could be efficiently ferried to the seafloor.

The k_dissbloom of the bloom-derived BSi is much higher than that required to reproduce the solid-phase BSi content profiles of the background system (k_diss) (Fig. 3). At stations B13, B14 and B15, the implemented k_diss values in the steady state simulations were 0.0055, 0.074 and 0.0105 yr⁻¹ respectively (Table S2, Fig. 2), which are not dissimilar to those estimated for deep-sea sediments from DSi pore water profile-fitting procedures and flow-through reactor experiments (0.006–0.44 yr⁻¹) (McManus et al., 1995; Rabouille et al., 1997; Ragueneau et al., 2001; Rickert, 2000). The inverse of the dissolution rate constant provides an estimate of the residence time or mean lifetime of the given pool of BSi (McManus et al., 1995), which suggests that the less reactive pool of BSi has a mean lifetime of 182, 13.5 and 95 years at stations B13, B14 and B15 respectively. These mean lifetimes are too long to influence the Si cycle on a seasonal scale, which must be < 1 yr⁻¹ to do so, as is the case for organic matter (Burdige, 2006). The model-derived estimates of k_dissbloom on the other hand would suggest a mean lifetime of approximately 20 d for the fresh BSi.

Therefore, this work suggests that there are at least two types of BSi in Barents Sea sediments, one less reactive pool that dissolves at a slower rate and one fresher, bloom-derived pool that is able to perturb the sediment pore water DSi stock on a seasonal timescale. This conclusion is compatible with findings from a study of the equatorial Pacific region (McManus et al., 1995) and observations from the Arabian Sea, indicating that the bulk sediment BSi content should not be treated as a single pool of uniform reactivity but should instead be separated into reactive and unreactive fractions (Rickert, 2000; Schink et al., 1975). Consistent with these conclusions, previous experiments have demonstrated that dissolution patterns of some diatom frustules can be best described by two k_diss values, an order of magnitude apart (Boutorh et al., 2016; Moriceau et al., 2009). These results indicate the presence of two phases of BSi within diatom frustules, denoting a potential physiological basis for the differentiation in reactivity of seafloor BSi.

3.3 What is the simulated benthic DSi flux, and how important is the contribution of bloom-derived BSi dissolution to the annual flux?

The simulated J_diff magnitudes (0.11–0.27 mmol Si m⁻² d⁻¹) that contribute to J_tot (Table 2) are within error of previously calculated J_diff values (0.10–0.37 mmol Si m⁻² d⁻¹) for these Barents Sea sediment cores (Ward et al., 2022). Thus, our model-derived benthic DSi fluxes are well within the range of a compilation of pan-Arctic shelf benthic DSi fluxes (−0.03 to +6.2 mmol Si m⁻² d⁻¹) (Bourgeois et al., 2017). Based on previous calculations and the simulation results presented here, we estimate that the mean DSi J_diff benthic flux magnitude for the Barents Sea is +0.23 (±0.11 1σ) mmol Si m⁻² d⁻¹, ranging from +0.08 to +0.54 mmol Si m⁻² d⁻¹.

J_tot at all stations is dominated by the molecular diffusive component (J_diff) (76 %–85 %), in agreement with simulated estimates of phosphate fluxes at the same stations (Freitas et al., 2020) (Fig. 5). J_tot at station B13 has the highest contribution from bioturbation (6 %), consistent with the highest experimentally determined bioturbation diffusion coefficient of the Barents Sea stations (Solan et al., 2020). The advective component of J_tot is negligible at all stations, while the bioirrigation element represents the greatest source of uncertainty in the simulated flux magnitudes as this parameter was not constrained in parallel to this work; thus a global value was assumed (Thullner et al., 2009). At station B15, the model-derived J_tot (+0.25 mmol Si m⁻² d⁻¹) is greater than J_BSi,in (0.17 mmol Si m⁻² d⁻¹) (Table 2). This observation points to the release of an additional source of DSi to the dissolving BSi, which is compatible with the hypothesis that LSi is being released into the pore water dissolved phase.

DSi benthic flux magnitudes were also calculated for the transient simulations carried out on station B14 to quantify the influence of fresh bloom-derived BSi. Dissolution of the fresher BSi has an immediate and significant effect on the benthic flux, doubling the steady state background value of +0.34 to +0.66 mmol Si m⁻² d⁻¹ within 1 week, which peaks 2 weeks after bloom material deposition at +1.84 mmol Si m⁻² d⁻¹, representing a 5-fold increase (an additional +1.5 mmol Si m⁻² d⁻¹) in the background steady state benthic flux magnitude. The prominent DSi peak then dissipates and becomes largely undetectable after 3 months (Fig. 4). The average DSi benthic flux at the SWI over the 12-week period is +1.07 mmol Si m⁻² d⁻¹, indicating that the bloom-derived BSi releases an additional +0.73 mmol Si m⁻² d⁻¹ to the overlying bottom water. The steady state and transient model simulations therefore suggest that the background benthic flux of Si from the benthos is +124 mmol Si m⁻² yr⁻¹ at station B14, while the additional contribution over the 12-week period sourced from fresh BSi dissolution is +61.3 mmol Si m⁻² (based on a rate of +0.73 mmol Si m⁻² d⁻¹). This estimate suggests that a minimum of 33 % of the total annual benthic flux of DSi discharging from the seafloor at station B14 is sourced from the deposition of fresh BSi during the 1-week MIZ bloom.

The contribution of bloom-derived BSi dissolution to the annual benthic DSi flux magnitude reported here (an additional +0.73 mmol Si m⁻² d⁻¹ over the 3 months) for station B14 is greater than that in Ward et al. (2022) (an additional +0.23 mmol Si m⁻² d⁻¹ over the same time interval), although the proportion is consistent across the two estimates (approximately one-third of the total annual benthic DSi flux). In part, this is due to the simulated fluxes incorporating the contribution from bioirrigation and bioturbation (J_tot). When using only the simulated J_diff component, an additional +0.46 mmol Si m⁻² d⁻¹ is estimated to be sourced from the bloom-derived BSi, which is more consistent with the observational data calculations. However, this disparity is also due to the nature of the simulated flux calculation. The model-derived benthic flux magnitudes are calculated at the SWI, whereas previous J_diff estimates are based on observational data of a much lower resolution, with the concentration gradient determined based on the DSi concentration in the core top water and in the sediment pore water at 0.5 cm depth. Furthermore, the simulated benthic flux estimates are based on a mean value derived from a weekly temporal resolution, which is not accessible in the observational data. However, both estimates can be used to draw a range of possible contributions from the bloom-derived BSi, and although there is a disparity in the benthic flux magnitude, both methodologies suggest that at least one-third of the annual DSi benthic flux at station B14 is sourced from the dissolution of BSi deposited after a short MIZ bloom.

Figure 5 Benthic DSi flux magnitudes at the SWI calculated from steady state simulations, including contributions from molecular diffusion, bioirrigation and bioturbation. B14 TP refers to the B14 DSi flux magnitude at the peak of the transient simulations (×30 BSi depositional flux, 15 yr⁻¹ k_dissbloom, 1-week bloom duration).

3.4 How much BSi is buried in the long term in the Barents Sea?

Traditionally, burial efficiencies of BSi were not included within Si budgets of the Arctic Ocean due in part to a low mean BSi content (< 5 wt %), as well as to a low estimated sedimentation rate (a few mm kyr⁻¹) (Brzezinski et al., 2021; März et al., 2015). However, more recently Arctic Ocean seafloor BSi burial efficiencies are being re-evaluated due to revised models of sediment accumulation rates (Brzezinski et al., 2021). Asymptotic BSi content profiles approach 0.2 wt % at 14.5 cm below the SWI and remain constant with depth at all three Barents Sea stations (Figs. 3, S2), suggesting that some BSi is buried and able to accumulate in the seabed. Assuming a sedimentation rate of 0.05–0.06 cm yr⁻¹ (Zaborska et al., 2008), BSi particulates deposited at the SWI would approach 14.5 cm depth within approximately 240–300 years due to sediment accumulation, where BSi dissolution appears to cease. This timescale is of a similar magnitude to the calculated mean lifetime based on k_diss from the model simulations at station B13 (182 years), indicating that some BSi from this pool could accumulate at depth. However, the steady state BSi pools at stations B14 and B15 have mean lifetimes of just 95 and 13.5 years respectively, theoretically precluding burial of this pool below 14.5 cm depth. The observation that BSi is present at ∼ 40 cm depth at all three stations with a concentration of ∼ 0.2 wt % indicates that there may be three fractions of BSi constituting the bulk BSi phase: fresh bloom-derived BSi represented by k_dissbloom, the less reactive background BSi (k_diss) and a non-reactive pool that contributes to the 0.2 wt % BSi at depth alongside any surviving component from the background pool. A mean BSi burial rate of 0.024 mmol Si m⁻² d⁻¹ is estimated here for the Barents Sea (Eq. 7, Table 2), corresponding to 0.012 Tmol Si yr⁻¹ for the whole shelf, assuming an area of 1.4×10⁶ km².

Global seafloor BSi burial efficiencies (calculated as the BSi burial flux divided by the deposition rate) range from 1 %–97 % (Westacott et al., 2021; Frings, 2017; Ragueneau et al., 2001, 2009; DeMaster, 2001; DeMaster et al., 1996; Liu et al., 2005), averaging ∼ 11 % (Tréguer et al., 2021; Frings, 2017). Typically, BSi burial efficiencies are much higher in coastal and shelf settings and underneath polar fronts, due in part to higher rates of sedimentation (DeMaster, 2001). Indeed, studies of the Southern Ocean (Indian sector and the Ross Sea) and the Peruvian margin uncovered a 4- to 30-fold increase in BSi burial efficiency at the continental shelf stations (22 %–58 %), relative to the open-ocean counterparts (2 %–17 %) (Dale et al., 2021; DeMaster, 2001; Ragueneau et al., 2001). BSi burial efficiencies estimated here for three Barents Sea stations (1.4 %–12 %, Table 2) are within the range of published values and similar to the global mean at stations B13 and B15 but low relative to other continental shelves. Our estimated BSi burial efficiencies are based on the same sedimentation rates employed in the model (Table S2), which are similar to the Barents Sea mean of 0.07 cm yr⁻¹ (Zaborska et al., 2008). However, Barents Sea sedimentation rates of up to 0.21 cm yr⁻¹ have been estimated since the last glacial period (Faust et al., 2021), which would significantly increase the estimated BSi burial efficiencies (5 %–33 %) (Eq. 7).

While the BSi burial efficiency calculated at station B14 is also within previously published values, it is much lower than the Barents Sea stations to the north and south. Solid-phase BSi contents in the surface intervals were determined from samples collected during the third cruise in summer 2019. As discussed, the pore water DSi profiles at B14 from this cruise are thought to be influenced by the dissolution of bloom-derived BSi, which would account for the elevated BSi content in the surface sediment relative to that at stations B13 and B15. The higher BSi content has likely resulted in an elevated estimation of J_BSi,in and thus a reduced burial efficiency, which would also explain why the model-derived estimate of the contribution of LSi dissolution to the DSi released from the solid phase is lower at station B14 (60 %), due to an elevated R_db, which is enhanced by the deposition of fresher BSi. In addition, this would also explain why the estimated mean lifetime of the BSi from the steady state simulations at station B14 is much shorter than for the other two sites (13.5 years vs. 95–182 years). If the measured surface BSi content at station B14 is indeed influenced by residual bloom-derived BSi, this would result in an overestimate of the background BSi reactivity in the model.

Relatively low burial efficiencies of BSi in sediments underneath oxygen-depleted bottom waters (7 %–12 %), similar in magnitude to those calculated here for the Barents Sea, have previously been attributed to low rates of bioturbation, resulting in less efficient export of BSi towards more saturated pore waters (Dale et al., 2021). Bioturbation coefficients were determined experimentally for the Barents Sea stations (2–6 cm⁻² yr⁻¹) (Table S2) and are much lower than might be expected, based on an empirical global relationship with water depth (∼ 24 cm⁻² yr⁻¹) (Middelburg et al., 1997). Furthermore, the impacts of the low rates of macrofaunal mixing on BSi burial efficiency are likely exacerbated by slow rates of sediment accumulation in the Barents Sea. High rates of BSi burial in the Bohai Sea (60 %) and Yellow Sea (42 %) are thought to be driven by high sediment accumulation rates (Liu et al., 2002), which as with bioturbation is much lower in the Barents Sea than might be expected based on an empirical global relationship with water depth (0.55 cm yr⁻¹) (Middelburg et al., 1997). The combination of low rates of macrofaunal mixing and sediment accumulation may therefore be the cause of the lower BSi burial efficiencies observed here relative to other continental shelves.

3.5 What are the implications of this work for the Arctic Ocean Si budget?

Brzezinski et al. (2021) uncovered an imbalance in the Arctic Ocean Si budget after carrying out an assessment using Si isotopes. The δ³⁰Si values measured in the main Arctic Ocean water mass inflow and outflows are similar (∼ +1.70 ‰; Giesbrecht, 2019; Liguori et al., 2020; Brzezinski et al., 2021), suggesting that the cycling of Si within the Arctic Ocean has little net effect on δ³⁰Si. Brzezinski et al. (2021) conclude that given the relatively isotopically light input from fluvial sources (+1.30 ± 0.3 ‰; Sun et al., 2018), as well as that emerging from seafloor sediments (+1.16 ± 0.11 ‰; Ward et al., 2022), balance must be maintained through the burial of isotopically light Si by BSi. However, a mass balance showed that the isotopically light inputs to the system are only partially offset by the burial of BSi (0.16–0.30 Tmol Si yr⁻¹), assumed to have a δ³⁰Si of +1.16 ± 0.10 ‰ (Brzezinski et al., 2021). There must therefore be an additional sink of isotopically light Si if the Arctic Ocean Si isotope budget is to maintain balance (Fig. 6). The absence of direct isotopic observations from some of the major gateways, including the Barents Sea shelf over which most of the DSi sourced from the Atlantic Ocean flows (Torres-Valdés et al., 2013), as well as a lack of data for the isotopic composition of BSi in Arctic Ocean sediments, must be addressed to confirm the mechanisms proposed by Brzezinski et al. (2021). Si isotopes measured in the weak alkaline leachate (0.1 M Na₂CO₃, δ³⁰Si_Alk) extracted from surface sediment sequential digestion experiments and measurements of δ³⁰Si in core top waters (Ward et al., 2022), coupled with reaction–transport modelling (this study) for stations B13, B14 and B15, contribute to the Arctic Ocean Si isotope dataset and help to fill these knowledge gaps.

Figure 6 The Arctic Ocean Si budget (Brzezinski et al., 2021) (left) and a proposed Si budget for the Barents Sea (right), including a benthic flux recategorisation (i.e. contributions from BSi and LSi) and AuSi burial. Boxes include flux magnitudes given in Tmol Si yr⁻¹ (top values) and the flux δ³⁰Si in per mille (‰; italicised bottom values). “In” and “Out” refer to the Si fluxes, discounting the water mass inflow and outflow (i.e. in (blue) includes rivers + LSi; out (red) includes AuSi + BSi). Grey boxes and arrows represent internal cycling. See Supplement Sect. 3 for further information on how the Barents Sea Si budget was calculated.

δ³⁰Si_Alk in the Barents Sea ranges from +0.82 ± 0.16 ‰ at station B15 to +1.50 ± 0.19 ‰ at B14 (Ward et al., 2022). The Na₂CO₃ leachate activates an operationally defined reactive pool of Si, thought to be associated with authigenically altered and unaltered BSi (Pickering et al., 2020). Molar $\mathrm{Al}/\mathrm{Si}$ ratios of the Na₂CO₃ support this concept, which fall within the range of that expected of BSi (Ward et al., 2022). δ³⁰Si values measured in the core top waters at the Atlantic Water station (B13, +1.64 ± 0.19 ‰) and Arctic Water station (B15, +1.69 ± 0.18 ‰) are similar to the composition of the main Arctic Ocean inflow and outflow water masses (∼ 1.7 ‰) (Giesbrecht, 2019; Brzezinski et al., 2021; Liguori et al., 2020) and heavier than that measured in the BSi deposited at the seafloor. Therefore, assuming the composition of the BSi below the BSi dissolution zone within the seafloor is similar to that at the SWI, Barents Sea sediments represent a sink of ²⁸Si, relative to the composition of the inflow waters. However, δ³⁰Si_Alk at stations B13 and B14 is still isotopically heavier than the Arctic Ocean riverine input (+1.30 ± 0.3 ‰), as well as the assumed composition of the BSi buried across the Arctic seabed (+1.16 ± 0.10 ‰) in the Si budget (Brzezinski et al., 2021) (Fig. 6). Through reaction–transport modelling we have estimated that between 2 % and 40 % of the sediment pore water DSi pool is sourced from the dissolution of BSi. Moreover, 2.9 %–37 % of the total amount of DSi released is reprecipitated as AuSi (Table 2). AuSi preferentially takes up the lighter isotope in the Barents Sea with a ³⁰ϵ of −2.0 ‰ to −2.3 ‰ (Table S2), thereby enhancing the preservation of BSi and further enriching the solid phase in the lighter isotope. Clays formed during weathering have a δ³⁰Si composition ranging from −2.95 ‰ to −0.16 ‰ (Opfergelt and Delmelle, 2012). The burial of AuSi alongside BSi could therefore account for some of the isotopic imbalance.

Our reaction–transport model study has also highlighted the important contribution of LSi dissolution to the sediment pore water DSi pool (60 %–98 %). If our findings on the dissolution of LSi are consistent across other Arctic shelves, a portion of the benthic DSi flux cannot be defined as internal cycling of Si and should be recategorised as an additional input to that from the major ocean gateways and discharge from rivers. It is currently estimated that the benthic flux of DSi across the whole Arctic Ocean seafloor is ∼ 0.39 Tmol Si yr⁻¹ (März et al., 2015); therefore between 0.23 and 0.38 Tmol Si yr⁻¹ may represent an input of Si rather than a recycling term. This recategorisation could account for the additional Si inputs required to close the Si budget, as currently 32 %–47 % (or 0.21–0.38 Tmol Si yr⁻¹) of the estimated net Si output is unaccounted for (Brzezinski et al., 2021). The addition of AuSi as an output to resolve the isotopic imbalance would offset to some extent the release of DSi from LSi, while the LSi input compounds the isotopic imbalance identified by Brzezinski et al. (2021). However, here we show that through AuSi precipitation acting as an additional sink for ²⁸Si, both mass and isotopic balance can be attained in our proposed Si budget for the Barents Sea (Fig. 6, Supplement Sect. S3). Future work should look to assess whether similar relationships exist between the dissolution of LSi and precipitation of AuSi on other Arctic Ocean shelves if we are to use these mechanisms to balance the pan-Arctic Si budget. Additional work could include empirical assessments, such as through batch or flow-through reactor experiments, to study dissolution (both BSi and LSi) and precipitation (AuSi) kinetics in Arctic Ocean sediments and to further examine and better constrain the relationship between benthic Si and Fe redox cycles.

4 Conclusions

In this study we quantify and disentangle the processes involved in the early diagenetic cycling of Si in the Arctic Barents Sea seafloor by reproducing Si isotopic and DSi concentration data from the solid and dissolved phase in a reaction–transport model (Figs. 2 and 3). Baseline simulations are able to reproduce the observational data well; however we have also shown that the benthic Si cycle is responsive on the order of days to the delivery of fresh BSi. Therefore, while the transient disturbances appear to be short-lived, future work should look to incorporate these processes into the baseline simulations.

Baseline model simulations also reveal that a significant proportion of the Si released from the solid phase within Barents Sea surface sediments is sourced from the dissolution of LSi (60 %–98 %) on account of the low BSi contents (0.26 wt %–0.52 wt %). Furthermore, we demonstrate that without the influence of the Fe redox cycle, which results in the release of Si adsorbed onto solid Fe (oxyhydr)oxides under anoxic conditions, the observed isotopic composition of the pore water DSi pool cannot be reconciled. Both the LSi and FeSi sources are depleted in the heavier isotope (−0.89 ‰ and −2.88 ‰ respectively), as demonstrated in a sequential digestion experiment (Ward et al., 2022), consistent with the observation that sediments of the Barents Sea represent a source of light DSi to the overlying bottom waters (Brzezinski et al., 2021). Of the DSi sourced from BSi, LSi and FeSi, we show that between 2.9 % and 37 % is reprecipitated as AuSi. Coupled with the observation that a significant proportion of the sediment pore water DSi pool is sourced from the dissolution of LSi, this hypothesis is significant for the regional Si budget. The dissolution of LSi represents a source of “new” Si to the ocean DSi pool, and the precipitation of AuSi inhibits exchange of pore water DSi with overlying bottom waters and therefore represents a sink term. These observations could require a recategorisation of a portion of the benthic flux in the Arctic Ocean Si budget, which is currently defined as a recycling term, as well as the inclusion of an additional Si sink. If LSi dissolution and AuSi precipitation are not exclusive to the Barents Sea shelf, the additional input and isotopically light output could account for both the isotopic imbalance and the remaining proportion of net Si outflow that is currently unaccounted for (Brzezinski et al., 2021).

Model simulations also highlight a dichotomy in the cycling of Si in the Barents Sea seafloor, which is hypothesised to occur on at least two timescales. Observational data at stations B13 and B15 can be reproduced by assuming a steady state dynamic, thus representing a background system, which is controlled by the release of Si into the DSi pool from LSi and the reprecipitation of DSi as AuSi. However, sampling across 3 years at station B14 has uncovered Si cycling on a much shorter timescale, controlled by the deposition of fresh phytodetritus. In this transient dynamic, the release of DSi is controlled by the dissolution of more reactive BSi. The processes occurring on the former steady state time frame will likely remain largely unaltered with further Atlantification of the Barents Sea due to the mineralogical control on DSi release, whereas the latter transient system is reliant upon the seasonal delivery of fresh BSi, which is subject to change as the community compositions of the MIZ and spring phytoplankton blooms shift to favour temperate Atlantic flagellate species (Neukermans et al., 2018; Orkney et al., 2020) or diatoms with lower silica content than polar species (Lomas et al., 2019). Furthermore, we have shown that the benthic DSi flux magnitude can increase 5-fold after a simulated 1-week bloom, which is calculated here to contribute a minimum of one-third of the total annual flux of DSi from the seafloor at station B14. Any perturbation in the delivery of bloom-derived, relatively reactive BSi to the seafloor could therefore be detrimental to the total annual supply of DSi from Barents Sea sediments.