The integrated benthic silicate flux in the Baltic Sea suggests a major land- or freshwater-derived reactive silicon source

Coastal marine environments are hot spots in the global marine silicon (Si) cycle. Dissolved silicate (dSi) is an essential macronutrient for diatoms, which often dominate primary productivity in temporal coastal seas and constitute a key food source for grazers. Even though benthic release of dSi influences the ecology of coastal marine areas, direct rate measurements of dSi mobilisation remain scarce. The Baltic Sea is no exception, and the spatial coverage of benthic dSi flux data is low and limited largely only to regional reports. We report data from 305 individual measurements (mostly in situ) of benthic dSi fluxes conducted in different basins and sediment types of the Baltic Sea during 2001–2021. Using the benthic dSi flux data in combination with literature values, representative average fluxes for various sediment types in the major basins of the Baltic Sea were determined. An areal extrapolation using Geographical Information System (GIS) tools suggests an integrated annual benthic release of 8520 metric kilotonnes (kt) of dSi for the entire Baltic Sea. This benthic release of dSi is about ten times higher than the reported riverine transport of dSi to the Baltic Sea. Furthermore, this benthic load, together with the reported annual burial rate, is more than three times higher than the autochthonous export production of biogenic silica out of the photic zone. The integrated benthic dSi release being substantially larger than that cycled by diatoms may explain the trend of the increasing dSi standing stock in most of the Baltic Sea basins which has been observed since the 1990s. Overall, other major sources of reactive Si (estimated to be 6390 kt yr⁻¹) to the sediment are suggested to exist, such as deposition of river and groundwater derived reactive (dissolvable) particulate amorphous and/or lithogenic Si. Our results strongly suggest that the biogeochemical Baltic Si cycle is more heavily influenced by reactive Si of terrestrial and freshwater origin than previously known.

1 Introduction

Dissolved silicate (predominantly silicic acid, dSi) is an essential macronutrient for diatoms, which often dominate the annual export productivity in temperate seas. The principal source of dSi to the hydrosphere is by weathering of lithogenic silicon (lSi) in bedrock material (Trapp-Müller et al., 2025 and references therein), and riverine input of dSi to coastal marine areas has fundamental influence on sea productivity. Recent studies also emphasize the supply of dSi to the ocean via submarine groundwater discharge (Ehlert et al., 2016b; Cho et al., 2018) and mobilisation from sandy beaches (Ehlert et al., 2016b; Fabre et al., 2019) as important sources (Tréguer et al., 2021). Rivers also transport large amounts of silicon as suspended matter to the sea. This flux of silicon was, however, previously regarded as negligible from an ecological perspective, as most particle-associated forms of silicon were assumed to dissolve too slowly to influence seawater dSi concentrations on ecologically relevant timescales (Tréguer et al., 1995). This view is under revision, and it is now clear that certain types of particulate Si dissolve relatively easily and may therefore be converted to dSi in freshwater watersheds and in seawater. These reactive Si-pools (rSi) may be grouped into amorphous Si (aSi) and lSi. aSi includes biogenic Si (bSi), which is formed by biological assimilation of dSi as well as certain inorganic fractions, for example Si associated with amorphous iron or aluminium particles (Saccone et al., 2007). The potential of lSi to influence marine dSi concentrations is uncertain, but recent studies show that up to 5 % of Si present in rock fragments transported to the sea may be easily dissolved in seawater (Tréguer et al., 2021 and references therein). Thus, it has become increasingly apparent that aSi and lSi are non-negligible fractions in river water and relatively soluble, and that the input of rSi to the world ocean must be revised upwards (Tréguer and De la Rocha, 2013; Tréguer et al., 2021). It is also increasingly recognized that marine dSi inventory is influenced by reverse weathering processes that take place during early sediment diagenesis and lead to the removal of dSi during formation of authigenic clay minerals (Michalopoulos and Aller, 2004; Ehlert et al., 2016a; Ward et al., 2022a; Aller and Wehrmann, 2024; Trapp-Müller et al., 2025). Together these processes indicate that sediment-seawater exchange processes must be revisited for an improved estimate of the marine Si budget.

These relatively recently explored mechanisms and fluxes related to the marine Si cycle have only started to be investigated in the Baltic Sea, a brackish, semi-enclosed sea in northern Europe. The latest estimates (Humborg et al., 2007; Conley et al., 2008; Papush et al., 2009) of the riverine dSi-load to the Baltic Sea amount to approximately 800–860 kt yr⁻¹ for the period 1980–2000. This is a significant reduction of the anthropogenically unperturbed (year 1900) riverine load of about 1300 kt dSi yr⁻¹ (Conley et al., 2008) and has been related to the development of hydroelectric power dams and eutrophication in freshwater watersheds (Humborg et al., 2008). Despite the abovementioned decline in riverine dSi delivery, reports of dSi limitation of diatom growth are few and limited to certain localized coastal areas (Olli et al., 2008), while in the open sea, diatom production has not been impacted (Wasmund et al., 2013). As a result of decreased riverine loading, decreasing trends of dSi concentrations were seen in Baltic waters from the early 1970s to the end of the 1990s (Sandén et al., 1991; Papush and Danielsson, 2006). In more recent years, this decline has ceased, and dSi concentrations are increasing in several Baltic basins (Skjevik et al., 2024).

Sediment-water fluxes of dSi can be another important source of silicon supporting pelagic productivity. However, there are uncertainties associated with the benthic exchange of dSi, both globally (Boynton et al., 2018) and in the Baltic Sea. This is mainly because there are only a limited number of experimental measurements of benthic dSi exchange, which in the Baltic Sea are almost exclusively limited to the Gulf of Finland (e.g., Conley et al., 1997; Almroth et al., 2009; Tallberg et al., 2017), and the Gdansk Basin and adjacent area (e.g., Thoms et al., 2018; Kendzierska et al., 2020; Borawska et al., 2022).

Recently, benthic dSi fluxes at the whole Baltic Sea scale were estimated based on vertical transport coefficients, obtained from water column distributions of radium (Ra) isotopes, and dSi gradients in the water column (McKenzie et al., 2025). The gradient flux method to estimate vertical fluxes is much less effort-consuming than direct flux measurements based on sediment incubations, and it can integrate over much wider bottom areas. It is also a valuable alternative when sediment incubations are difficult to carry out, such as e.g. on stony or rocky bottoms, or stiff sandy sediments.

We incubated sediment to explore spatial variability of benthic dSi fluxes from more than 300 direct flux measurements, primarily conducted in situ by benthic chamber landers, both in coastal and offshore areas of the major basins of the Baltic Sea during 2001–2021. Representative flux values for various sediment surface types in the different basins were generated from the benthic flux data in combination with literature values from similar measurements. Using GIS tools for area calculations, these representative flux values were used to produce an estimate of the integrated annual benthic flux of dSi on a Baltic Sea system scale.

2 Methods

2.1 Study site

The Baltic Sea (Fig. 1) is an inland, brackish sea in northern Europe. Water of higher salinity enters the Baltic Sea via the Danish straits including the Öresund in the south-west and brackish conditions are maintained due to a long water residence time and substantial freshwater input via rivers, primarily in the northern parts of the Baltic Sea.

Figure 1 Map of the Baltic Sea. Stations where flux measurements were conducted are indicated with circles and sites with values from the literature are marked with triangles. Background colours indicate EMODnet seabed substrate type (see legend in figure) and assessed substrate class at sampling points are indicated using the same colours. Dark blue colour on main map shows archipelago zones.

The Baltic Sea is sub-divided into gulfs and basins of which the northernmost Gulf of Bothnia is connected to the largest basin, the Baltic Proper, via relatively shallow archipelago areas (Åland Sea). The Gulf of Bothnia is further sub-divided into the Bothnian Bay and the Bothnian Sea, which are also partly separated by a 25 m deep sill. The hydrography of the Baltic Sea results in a pronounced north-south salinity gradient. Typical summertime bottom water salinities in the Gulf of Bothnia range from about 3.5 or less in the northern Bothnian Bay up to approximately 7 in the southern Bothnian Sea. Similarly, the bottom water salinity increases from 3 in the inner most part of the Gulf of Finland to approximately 9 at the boarder to the Baltic Proper, and up to around 15 in the southernmost basin of the Baltic Proper (the Arkona basin). A permanent halocline is present in the Baltic Sea. The halocline depth in most sub-basins is 60–80 m and the salinity above and below the halocline differs by approximately 2–4 (Leppäranta and Myrberg, 2009).

This paper is based on numerous benthic flux measurements in the Baltic Proper and the gulfs of Bothnia and Finland (Table 1; Fig. 1). The measurements were conducted in all subbasins of these Baltic Sea areas, excluding the Gulf of Riga, the Gdansk Basin and the Arkona Basin, and during all months of the annual cycle except March and December. The majority of measurements were carried out using in situ benthic chamber landers, primarily in open sea areas, but landers were also deployed in coastal areas in the Bothnian Bay and the Northern Baltic Proper. Ex situ measurements were performed in coastal (Himmerfjärden) and open sea sampling points in the Western Gotland Basin. Some of the results were published previously: Stations NWBP, KH104, KH58 (Ekeroth et al., 2016a, b), stations D, E and F (Hall et al., 2017; Hylén et al., 2021), but have not previously been analysed at the whole Baltic Sea scale. We also compiled ex situ dSi flux values from the literature.

Table 1 List of stations where benthic dSi flux measurements were conducted in this study. n is number of flux measurements, T is temperature, S is salinity and O₂ is dissolved oxygen concentration in bottom water. T, S and O₂ were measured by sensors on the benthic landers (Kononets et al., 2021). “–” means O₂ data not available.

* According to EMODnet (2021).

2.2 Benthic lander measurements

249 out of the 305 benthic flux observations were carried out in situ using benthic landers. Two different chamber landers were used for these measurements in the present study – the small and the big Göteborg landers (Kononets et al., 2021). The landers were equipped with two (small lander) or four (big lander) identical, open-bottomed (400 cm²) benthic incubation chambers which penetrate the sediment and thereby incubate the sediment and overlying water mass (see e.g., Ståhl et al., 2004; Kononets et al., 2021). The inner height of the chambers was 35 cm, and the penetration depth of the chambers was on the order of 15 cm, resulting in a ∼ 20 cm water height of the incubated water mass (total water volume ca. 8 L).

Ten syringes were connected to each chamber by ∼ 10 cm long plastic tubes. One of these syringes was used to inject a known volume (ca. 60 mL) of de-ionized water into each chamber a few minutes after the pre-incubation period. The resulting decrease in salinity, measured by conductivity sensors (Aanderaa Data Instruments, Norway) mounted in the chambers, was used to calculate the volume of the incubated water. The remaining nine syringes were used to withdraw water from the chambers at pre-set times during the incubation period. After recovery of the lander, the syringe sample water was filtered (0.45 µm pore size cellulose acetate filters) and were kept cold until spectrophotometric determination of dSi according to Grasshoff et al. (1999) modified for Alpkem SFA (O. I. Analytical Flow Solution IV, model 319529 with a precision of better than 4 %). The lander syringe samples from the 2020 and 2021 expeditions were treated in the same way but dSi was determined by the automated molybdate-blue method (Strickland and Parsons, 1972) with a Smartchem 200, AMS™ discrete analyser with a precision of 4 %. The concentration change versus time for each chamber and deployment was determined by least squares linear regression analysis. The benthic dSi flux was then estimated from the regression slope and chamber volume. When, in some cases, the rate of concentration increase slowed down during incubations, the data points at the end of incubation were not used to calculate the regression slope. The lander measurements and flux evaluation methodology are described in detail in Kononets et al. (2021).

2.3 Ex situ measurements

56 flux measurements were carried out by means of ex situ core incubation experiments. Sediment for these incubations was retrieved using a multiple corer (K.U.M. Umwelt- und Meerestechnik, Kiel). Intact cores with undisturbed sediment surface were subsampled using transparent plastic liners (4.6 cm internal diameter, 30 cm height). The sediment height in the subsampled cores ranged between 12 and 15 cm. The sediment cores were transferred into a 25 L incubation tank that was previously filled with ambient bottom water, situated in a temperature-controlled room kept at bottom water temperature (Table 1). The cores were left uncapped to acclimatize in the dark for 6 to 12 h, when the water overlying sediments was constantly stirred with magnetic bars (∼ 60 rpm.).

A minimum of four cores from each station were closed with caps with stirring on to determine net fluxes of dSi following the procedure described in Bonaglia et al. (2013). The incubation time was set after a preliminary evaluation of the sediment oxygen consumption rates, to avoid oxygen concentrations decreasing below 20 % of the initial value. The sediment cores were incubated for a period of 8 to 12 h. At the beginning and at the end of the incubation, water samples were taken from the incubation tank (n=4) and from each core, respectively. The water samples were immediately filtered (0.45 µm pore size cellulose acetate filters) and analysed for dSi concentrations according to Grasshoff et al. (1999) as described above. Net fluxes across the sediment-water interface were calculated from the difference in dSi concentrations in the water column through the incubation period and the volume of water column.

2.4 Sediment sampling and analysis

Sediment cores were collected with a multiple corer or a GEMAX corer (Kart Oy, Finland) at each station. The cores were sliced in air immediately after collection at 0.5 cm vertical resolution down to a sediment depth of 2 cm. A set of sediment samples were kept in closed preweighed vials in a refrigerator until the end of expeditions, when they were weighed, dried to constant weight and weighed again. Water content was calculated from the weight of water loss and the original wet weight. Another set of sediment samples were frozen on-board, and freeze dried at the end of expeditions. They were then analysed for total carbon or total organic carbon content as described by Nilsson et al. (2019, 2021).

2.5 GIS analysis

Spatial analysis was performed using the software QGIS (QGIS Development Team, 2016) to calculate the areal extent of various seabed surface types within sub-basins of the Baltic Sea. We here used a very similar procedure to the one reported by Hylén et al. (2025) for upscaling of benthic dissolved inorganic phosphorus fluxes in the Baltic Sea. The division of sub-basins was based on HELCOMs definitions (HELCOM, 2022). The sub-basins along the Swedish and Finnish coasts were further divided into coastal and open sea parts. The delimitation of coastal and open sea areas followed definitions of coastal areas by Swedish and Finnish authorities (Swedish Meteorological and Hydrological Institute, SMHI, and Finnish Environment Institute, Syke).

Calculations of the areal extent of various seabed surface types within the sub-basins were based on the EMODnet Seabed Substrate map 1:1M (EMODnet 2021). Five seabed surface type classes are defined in the EMODnet substrate map (after Kaskela et al., 2019) and the areal extent of each class in each subbasin were calculated using the function intersection in QGIS. Prior to this, the layers were reprojected (EPSG:3035–ETRS89-extended/LAEA Europe) and invalid geometries were fixed manually.

2.6 Classification of seabed substrate type at sampling stations

Classification of the seabed substrate type at a sampling station was primarily based on geographical position in relation to the EMODnet seabed substrate map. It should be noted that the EMODnet seabed substrate map was compiled at a 1:1 000 000 scale, using data with a smallest cartographic unit of 4 km² (Kaskela et al., 2019). Therefore, it is not meant to be used for determining the seabed substrate type for individual sampling points or similar small geographical scales. In cases when observations and measurements (water content, total carbon content) of the actual sediment from a sampling station deviated from the seabed substrate type classification according to EMODnet, it was reclassified based on observational data and measurements (Table 2).

Table 2 Sampling station meta data. Coastal stations are located within the archipelago as defined by Swedish and Finnish authorities. EMODnet seabed surface type shows the seabed surface type according to the geographical position. Determined seabed surface type shows the assumed seabed surface type which was used in the upscaling of benthic fluxes. Sampling stations which were assumed to have a different seabed substrate type than indicated by EMODnet are marked in bold. TC shows the average content (% of dry weight) of total carbon in the upper 2 cm of sediment. When only the content of total organic carbon (% of dry weight) was available, it is indicated with TOC in the table. BB = Bothnian Bay, BS = Bothnian Sea, NBP = Northern Baltic Proper, GOF = Gulf of Finland, WGB = Western Gotland Basin, EGB = Eastern Gotland Basin, Bornh. B. = Bornholm Basin, GOR = Gulf of Riga, GB = Gdansk Basin, and AB = Arkona Basin.

* Reported by Fredriksson et al. (2024).

2.7 Calculation of annual benthic dSi loads

The term sub-section is here used for geographical regions of the Baltic Sea used in the upscaling of benthic fluxes. Each sub-section includes one, a part of one, several, or parts of several Baltic Sea sub-basins (see Table 3). In the formation of sub-sections, sub-basins with lack of flux measurement data (such as the Quark and the Åland Sea) were merged with adjoining sub-basins. Internal load estimates were calculated per sub-section, and coastal (archipelagic) and open sea areas were treated individually (Table 3).

Table 3 Sets of station data used for calculation of representative average fluxes for different seabed surface types in sub-sections of the Baltic Sea. Superscripts at station names refer to literature data.

¹ Borawska et al. (2022). ² Gogina et al. (2018). ³ Tallberg et al. (2017). ⁴ Conley et al. (1997). ⁵ Aigars et al. (2015). ⁶ Kendzierska et al. (2020). ⁷ Thoms et al. (2018).

Sub-section specific average fluxes and standard deviations (weighted based on number of measurements per sampling station) per sub-section and seabed substrate type were produced. These fluxes (mmol dSi m⁻² d⁻¹) were multiplied with the total area of the corresponding seabed surface type within each sub-section (m²) and summed yielding load estimates for each individual sub-section (mmol dSi d⁻¹). Areas with unknown or restricted information regarding seabed surface type were assumed to have the same proportional composition as their respective sub-basin. Benthic flux data were lacking for certain combinations of specific seabed surface types and specific sub-sections. In these cases, relevant flux data from one or several (mostly adjacent) sub-sections were used as substitution. In addition to the data from measurements by the authors of the present paper, literature values of benthic dSi fluxes (Conley et al., 1997; Aigars et al., 2015; Tallberg et al., 2017; Gogina et al., 2018; Thoms et al., 2018; Kendzierska et al., 2020; Borawska et al., 2022) were also used in the calculation of sub-section specific average fluxes per seabed surface type. All literature values were based on ex situ incubation measurements using similar methodology as described in Sect. 2.3. Data from McKenzie et al. (2025) are not included so as to restrict our analyses solely to direct flux measurements.

3 Results

3.1 Benthic fluxes

The dataset of benthic dSi fluxes of this study comprises 305 individual flux values from core or chamber incubations. Representative examples of the evolution of dSi concentration versus time during in situ chamber incubations (Fig. S1) are provided in Supplementary Information. All in situ and ex situ flux values obtained in this study are available from Zenodo (Ekeroth et al., 2026, https://doi.org/10.5281/zenodo.18484732). Lander measured dSi fluxes at stations D, E and F in the EGB in 2015, 2016 and 2017 were excluded from this dataset since these fluxes were influenced by the 2014-15 Major Baltic Inflow (Hylén et al., 2021). Similarly, the measured dSi fluxes at station PV1 in the NBP in 2003 and 2005, and at the coastal station KH104 in the NBP in 2013, were excluded from this dataset since these fluxes were influenced by transient oxygenation-deoxygenation events. All but one fluxes were directed out of the sediment. Station-averaged fluxes ranged from 0.3–9.3 mmol m⁻² d⁻¹ (Fig. 2) with an overall average for the Baltic Sea of 3.7 mmol m⁻² d⁻¹, which is in the range of other brackish marine areas (∼ 2–6 mmol m⁻² d⁻¹) such as the Black Sea (Friedl et al., 1998; Friedrich et al., 2002), Chesapeake Bay (Cowan and Boynton, 1996), Mobile Bay (Cowan et al., 1996) and San Francisco Bay (Grenz et al., 2000; Cornwell et al., 2014). The highest fluxes (station averaged flux >7 mmol m⁻² d⁻¹) were found in muddy and relatively deep areas of the Eastern and Western Gotland Basins (EGB stations D and F, and WGB stations T002, T003, T006, W009 and W010), as well as in the deep coastal basin Kanholmsfjärden in the northwestern Baltic proper (station KH104).

Figure 2 (A) Benthic dSi fluxes in mmol m⁻² d⁻¹ (mean values per station ±1 standard deviation). Bars are grouped according to seabed substrate type and Baltic Sea sub-basins. Basin abbreviations are explained above. See Fig. 1 for locations of stations. (C) = coastal station. (B) Benthic dSi fluxes in mmol m⁻² d⁻¹ (mean values per station ±1 standard deviation) in different basins and at different sediment surface types compiled from the literature. The station names also appear in Table 3, where references are given.

3.2 Annual benthic dSi load

The total area and the areal proportion of sediment surface types for each sub-section included in the GIS-analysis is given in Table 4 and classification of sediment surface type for each sampling station is given in Table 2. Mud-muddy sand (type 1) or mixed sediment (type 4) are the dominating seabed substrate types in all sub-sections except the Arkona basin where sand (type 2) is dominating. Relatively large areas in the open Bothnian Bay and in the Gdansk Basin are also covered by sand. Relatively large proportions of the coastal sub-sections have an unknown or unclassified seabed surface type in the EMODnet data (Table 4). However, in absolute terms, compared to the open sea, these areas are very small (2 % of the total area).

Table 4 Total area and coverage (proportion of area) of various seabed substrate types per sub-section in the Baltic Sea used for calculation of internal dSi loads.

The seabed surface type at specific sampling stations generally agreed with the EMODnet map, based on geographical position (Table 2). For some stations, however, the seabed surface type according to the EMODnet map was erroneous, or at least questionable. As described above (Sect. 2.6), such discrepancies were anticipated due to the relatively low resolution of the EMODnet data. In the upscaling of benthic fluxes, the indicative classification of seabed substrate (EMODnet based on geographical position) was changed for ∼ 25 % of the sampling stations (Table 2). Grain size distribution data is not available, so reclassification of seabed surface type was mainly based on measured total carbon content, water content and/or other observations from field work.

Table 3 shows the data used for calculation of average fluxes for various seabed substrate types in each sub-section. Average flux values per sub-section are shown in Fig. 3.

Figure 3 Benthic dSi fluxes in mmol m⁻² d⁻¹ per seabed surface type used for calculation of the integrated annual dSi load in various sub-sections of the Baltic Sea. The bars show weighted average values based on original data in the present study and literature values (see Table 4). Error bars show 95 % confidence intervals of the weighted average values. Lack of error bars is due to lack of information regarding variability (Gulf of Riga) or that n=1 (BS, Quark (coastal). The different sediment surface types are colour coded.

Using the seabed surface- and sub-section-specific average flux values in Fig. 3, the annual integrated benthic dSi flux for each sub-section was calculated (Fig. 4).

Figure 4 Annual benthic dSi fluxes in kt yr⁻¹ in the sub-sections of the Baltic Sea. Please note the logarithmic scale. Error bars show 95 % confidence intervals based on flux uncertainty. Lack of error bars is due to lack of information regarding variability (Gulf of Riga) or that n=1 (BS, Quark (coastal)). The different sediment surface types are colour coded. The actual fluxes per basin and sediment type (kt yr⁻¹) are given in the table under the bars.

Figure 4 (displaying fluxes on a logarithmic scale) clearly shows that coastal areas have an almost negligible annual benthic flux compared to most open sea sections. This is due to their limited areal coverage, as the average fluxes (per m²) in coastal areas are similar to those in open sea areas (Fig. 3). Mud-muddy sand areas (Type 1) have the greatest contribution to the annual benthic flux in most sub-sections (Fig. 4). An exception is the Gulf of Bothnia (Bothnian Bay and Sea, Quark and Åland Sea) where more dSi is released from mixed sediment (Type 4) than from mud-muddy sand (Fig. 4). This is in part due to that mixed sediment cover a relatively large area of the seafloor in this sub-section (58 %, Table 4) compared to most other sub-sections. Also, the load estimate from mixed sediment in GOB is based on a single sampling station (GOB3) and is therefore uncertain.

The estimated total integrated annual benthic dSi-load from all sub-sections amounts to 8520 kt. About 47 % of the load (4000 kt) is released from mud-muddy sand sediment in the EGB and from mixed sediment in the GOB (Fig. 4). The estimated benthic release from mud-muddy sand sediment in the EGB is based on measurements at six stations and is therefore better constrained than the estimate for GOB, as mentioned above.

The annual benthic load per unit area varies between about 6–35 t km⁻² yr⁻¹ between sub-sections (Fig. 5).

Figure 5 Annual benthic dSi fluxes per unit area (t km⁻² yr⁻¹) in the various subsections of the Baltic Sea. Error bars show 95 % confidence intervals of weighted average values.

4 Discussion

4.1 Input of Si to the Baltic Sea

The annual integrated benthic release of dSi in the Baltic Sea system (8520 kt) is equivalent to 78 % of the dSi standing stock in the water column (10970 kt) and is approximately a factor of ten higher than the annual riverine dSi input (which is 822 kt yr⁻¹, Papush et al., 2009). Thus, our benthic load estimate indicates that benthic dSi recycling can renew the dSi water column stock in less than two years, and it appears fundamental for maintaining pelagic dSi concentrations in the Baltic Sea. Previous estimations of the residence time of dSi range between approximately 7–11 years in the whole Baltic Sea (Wulff et al., 1990; Papush et al., 2009), and between a few months to three years in individual basins (Papush et al., 2009). These previous estimates of residence times are based on net losses of dSi (calculated from concentrations differences between inflows and outflows to/from model boxes) in relation to total reservoirs of dSi, but not on internal recirculation processes like in this study.

It should be noted that estimates of the riverine input of dSi (Papush et al., 2009) do not include dissolvable aSi or lSi, which may potentially be converted to dSi with some delay (Raimonet et al., 2023) and thus influence the standing stock. It is also not known to what extent dSi-input mediated by submarine groundwater discharge (Donis et al., 2017), which is demonstrated in the Baltic Sea for DOC (Szymczycha et al., 2014) and nutrients (Kreuzburg et al., 2023), has an influence on the Baltic dSi standing stock. Presently, these key mechanisms in the global marine Si cycle (Tréguer et al., 2021) have only started to be explored in the Baltic Sea. For example, it has been speculated that riverine input of particulate biogenic silica (bSi) may be an important additional source of dSi upon dissolution in the Baltic Sea, particularly from nutrient-rich catchments (c.f., Conley, 1997), such as from the rivers Vistula, Nemunas and Daugava in the Gdansk Basin, Curonian Lagoon and Gulf of Riga, respectively, where some of the highest concentrations of bSi in Baltic Sea sediments are found (Emelyanov, 2014). Lehtimäki et al. (2013) estimated that up to 40 % of rSi transported to the Baltic Sea is released by rivers in the form of bSi, which suggests that external Si sources to the Baltic Sea are severely underestimated, as existing dSi-budgets for the Baltic Sea (Wulff et al., 1990; Papush et al., 2009) do not account for bSi entering the system. There have been no attempts to quantify riverine or groundwater inputs of other reactive forms of amorphous or lithogenic Si to the Baltic Sea.

4.2 Towards a benthic rSi mass balance

We estimate a deposition of rSi on the Baltic seafloor from the sum of the integrated benthic dSi load and bSi burial under the assumption of steady state. Existing bSi burial estimates for the Baltic Sea amount to approximately 610 kt bSi yr⁻¹ (Conley et al., 2008; Papush et al., 2009). These relatively old estimates do not include dSi-losses by reverse weathering (c.f., Michalopoulos and Aller, 2004; Pickering et al., 2020; Aller and Wehrmann, 2024), which on a global scale is the greatest rSi-sink in the continental margin zone (Tréguer et al., 2021). Adding our calculated integrated benthic dSi release (8520 kt yr⁻¹) to the existing burial value, gives a total annual rSi deposition to Baltic Sea sediments of 9130 kt yr⁻¹ and an overall apparent burial efficiency (i.e., bSi burial divided by rSi deposition) of 7 %. However, since burial estimates of other forms of rSi than bSi do not exist, this calculated burial efficiency is likely an underestimate, and it should be interpreted with caution. Using the reported bSi burial rate and our calculated bSi export out of the photic zone (assuming a negligible dissolution during transit through the water column) provide an estimate of the bSi burial efficiency of ($(\mathrm{610}/\mathrm{2740})\times \mathrm{100}$) 22 %. It should be noted that this is an integrated estimate for the entire Baltic Sea system, and that regional variations should be expected. As a comparison, Tallberg et al. (2017) concluded that the burial efficiency for aSi in the Gulf of Finland is lower than 50 %. More studies of burial of rSi, including also other forms of rSi than only bSi, are needed before reliable rSi burial efficiencies can be quantified on a whole Baltic Sea ecosystem scale.

We further assessed whether autochthonous bSi sources alone can sustain the calculated total rSi deposition (9130 kt yr⁻¹). Ostrowska et al. (2022) estimated, based on remote sensing data, average annual primary production rates of 251–298 mg C m⁻² d⁻¹ with an overall average value of 266 mg C m⁻² d⁻¹ for the entire 2010–2019 period in the Baltic Sea. Using the latter value and the total area of our study area (373 124 km², Table 4) yields an average total primary production estimate of ∼ 36 200 kt C yr⁻¹. Assuming that the export production constitute one third of the total primary production (Heiskanen and Tallberg, 1999; Gustafsson et al., 2013), that diatoms are responsible for 75 % of the export production of C in this type of environment (Nelson et al., 1995), a C:Si ratio of marine diatoms of 3.3 per weight (Brzezinski, 1985), and a negligible dissolution of bSi during sinking from the bottom of the photic zone to the sea-floor (which is reasonable since the average depth of the Baltic Sea is relatively shallow, and the sinking rate of diatom debris is quite fast), an approximate estimate for the annual marine diatom bSi sea-floor deposition in the Baltic Sea is $(\mathrm{36}\mathrm{200}\times \mathrm{1}/\mathrm{3}\times \mathrm{0.75})/\mathrm{3.3}=\mathrm{2740}$ kt (Fig. 6). This corresponds to 30 % of our estimate of the annual benthic rSi deposition (9130 kt yr⁻¹). Thus, we suggest that the rSi deposition in the Baltic Sea is sustained not only by marine bSi produced in Baltic surface waters, but primarily by deposition of other forms of aSi and/or lSi of presumably terrestrial origin.

Our results suggest that about 3.1 ($\mathrm{8520}\mathrm{kt}/\mathrm{2740}\mathrm{kt}$) times more dSi is released annually from the Baltic seafloor than is being cycled by diatoms. What may the fate of this “excess” dSi be? Marine environmental monitoring (Skjevik et al., 2024) has shown an increasing trend of the water column dSi standing stock (winter condition) in all basins of the Baltic Sea (no data available for the GOF) from 1994 (from 2000 for the BS) until 2023. The largest increase is observed in the BS. One likely explanation for the increasing dSi standing stock may be the “excess” benthic dSi release suggested by our results. This would then imply that most of the enrichment of dSi in the Baltic basins is rSi of terrestrial origin which is dissolvable and recycled in sediments

The dSi net export from the Baltic Sea to the Kattegat has been estimated to be 225 kt dSi yr⁻¹ (Papush et al., 2009). This estimate is at least 17 years old and the dSi export may today be considerably higher. If this is the case, much of the “excess” dSi released from Baltic sediments is exported to the Kattegat (and likely beyond). In summary, the fate of the “excess” benthic dSi release suggested by this study may be accumulation in the water column of the Baltic basins and export to the Kattegat.

4.3 Mass balance uncertainty

The measured benthic fluxes show spatial variability on each sediment type. For example, on mud-muddy sand, the station averaged dSi flux varies between about 1 and about 10 mmol m⁻² d⁻¹ across the Baltic Sea (Fig. 2A). Mud-muddy sand spans over large areas including transport bottoms, deep accumulation bottoms and this sediment type contains various degrees of sand (up to 90 %) and gravel (up to 5 %) content per dry weight (Kaskela el al., 2019). It also covers large areas of the seafloor in all of the major Baltic Sea basins with various water depths, trophic states (mesotrophic to highly eutrophic) and hence a large range of vertical particle fluxes. So, this sediment type encompasses the large natural spatial variability of the Baltic Sea, which also is reflected in the benthic dSi fluxes.

Our mass balance has a considerable degree of uncertainty. The sum of the uncertainty in annual benthic fluxes, based on the confidence intervals shown in Fig. 4, amounts to 2099 kt yr⁻¹ (i.e., 25 % of the annual benthic load). This can be viewed as a lower estimate of the uncertainty, based on confidence intervals of average benthic flux values. Another source of error is misclassification in the EMODnet substrate map. This source of error is likely relatively small as the composition of the seafloor on a basin-wide geographical scale is rather well known, even though misclassifications are evident on the local scale (Table 2). The largest potential error term is likely related to the upscaling of benthic fluxes for open sea areas of the Gulf of Bothnia (GOB, consisting of the BS and the BB), where benthic flux data is available from only three stations which are all located in deep waters within each respective subbasin. Open sea areas of the GOB make up almost one third of the Baltic Sea (Table 4) and, according to our mass balance, the integrated release of dSi from the GOB similarly makes up approximately one third of the total for the whole Baltic Sea. Thus, the estimated load for the GOB appears reasonable, but errors in average fluxes per sediment type would undoubtedly have an influence on the total mass balance. In particular, the average flux for mixed sediment in the GOB is based on one sampling station (GOB3) with a relatively high average value (2.6 mmol m⁻² d⁻¹, n=6, range 1.9–3.1 mmol m⁻² d⁻¹), compared to average values for mixed sediment in other open sea areas (0.3–1.5 mmol m⁻² d⁻¹). To illustrate the impact of the average flux value for mixed sediment in the GOB (from measurements at station GOB3), the integrated total load of dSi would decrease from 8520 to 6870 kt yr⁻¹, should the lowest value for mixed sediment (0.3 mmol m⁻² d⁻¹) be used instead. The difference (∼ 1650 kt yr⁻¹ or about 20 % of the integrated load estimate) illustrates the maximum error due to under-sampling mixed sediment in the GOB. In combination with the above-mentioned uncertainty estimate based on confidence intervals, an approximation of the total uncertainty is maximum 40 % of the annual benthic load. Still, even with a 40 % lower estimate for the integrated benthic load autochthonous bSi export alone would be insufficient to “fuel” the deposition of reactive Si in our benthic mass balance. Thus, our conclusion that other sources contribute substantially to the deposition of rSi is relatively unaffected by mass balance uncertainties.

One possibility is that an older bSi pool deposited earlier (i.e. a legacy pool) exists in Baltic sediments, dissolves and contributes to our integrated benthic dSi flux. If dissolution of this old legacy bSi pool is important, the need for rSi from land, including groundwater, will be less to match the sum of the integrated benthic dSi flux and the reported bSi burial rate. However, we are not aware of any observations of a legacy bSi pool in Baltic sediments, nor the quantity of it or any change of this quantity. This is due to the fact that published data on solid phase bSi in Baltic sediments are scarce, and on other forms of particulate rSi are largely absent. For this reason, we assume that the contribution to the integrated benthic dSi flux from any legacy bSi pool is quantitatively of less importance.

Consumption of dSi by authigenic clay formation is removing dSi from the porewater, which would make the benthic dSi flux lower. The influence of this process on the benthic dSi flux is included in the measurements we performed (since the fluxes we measure are net fluxes), and thus so also in the integrated annual dSi flux we have presented. However, it is not the aim of this study to specifically quantify this process or to show evidence of its importance in the Baltic Sea.

McKenzie et al. (2025) measured water column dSi gradients and calculated vertical transport coefficients (K_z) from distributions of short-lived Ra isotopes in the water column, which resulted in an integrated, whole-Baltic benthic dSi release of 23 000–35 000 kt yr⁻¹. This first-order estimate is about 3 to 4 times higher than our estimate. However, the K_z reported by McKenzie et al. (2025) spans three orders of magnitude and this high variability in K_z is most likely associated with the strong vertical water density stratification, which is characteristic for the Baltic Sea.

4.4 Potential sources of aSi and lSi to Baltic Sea sediments

The difference between our estimate of total rSi deposition and annual diatom bSi export production amounts to 6390 kt Si yr⁻¹ and may be explained by deposition of allochthonous riverine bSi (including for example phytoliths and freshwater diatoms), other forms of reactive aSi, and lSi that constitute potentially underestimated sources of reactive Si to coastal marine ecosystems (e.g., Ragueneau et al., 2006; Tréguer and De La Rocha, 2013; Tréguer et al., 2021; Hatton et al., 2023). Up until recently, aSi and lSi fractions were considered to have too low dissolution rates to have any significant impact on marine dSi availability on decadal timescales (Papush, 2011; Tréguer and De La Rocha, 2013). Porewater dSi is produced from dissolution of all the different forms of particulate Si that exist in the sediment and is dissolvable (reactive). So how much dSi (contributing to the benthic flux) that originates from dissolution of one form of particulate Si (e.g. lSi) or another (e.g. bSi) is a function of the relative abundance of each of the reactive forms of particulate Si in a given sediment and their respective dissolution rates. However, this cannot be resolved with the benthic flux measurements we performed. Determinations of the Si isotopic composition of the dSi in the benthic chambers may have helped to disentangle the origin of the dSi in the benthic flux, but this study did not include isotopic determinations. Several research findings in recent years have, however, shown that lithogenic forms of Si may have a more significant influence on benthic dSi fluxes than previously thought (Ehlert et al., 2016b; Fabre et al., 2019; Geilert et al., 2020). For example, Ward et al. (2022b) showed that 60 %–98 % of the dSi-pool in the porewater was sourced from dissolution of lSi in Barents Sea sediments. In line with these latter studies and overall, our results suggest that allochthonous bSi and other forms of reactive aSi, as well as lSi, of terrestrial and freshwater origin have a greater influence on the Baltic Sea Si cycle than previously acknowledged.

5 Conclusions

In this study, we used sediment-water fluxes of dSi to quantify the benthic rSi recycling in the Baltic Sea system. Station-averaged fluxes ranged from 0.3–9.3 mmol m⁻² d⁻¹ with an overall average for the Baltic Sea of 3.7 mmol m⁻² d⁻¹. By extrapolating the large data set of benthic dSi fluxes in different subsections and sediment types in the Baltic Sea by GIS upscaling, we estimated an annual integrated benthic load of 8520 kt dSi. This number is about ten times higher than the presently reported riverine dSi load. Assuming steady state and using the reported bSi burial rate, deposition of 9130 kt rSi yr⁻¹ are needed to sustain this flux, which is about 3.3 times higher than the expected autochthonous export production of bSi by diatoms out of the photic zone. Thus, our results indicate a substantial missing source of 6390 kt yr⁻¹ of rSi to sediments of the Baltic Sea. Based on recent literature on other coastal marine areas, riverine and groundwater input of particulate amorphous and/or lithogenic Si would contribute in a major way to the depositional flux of reactive (dissolvable) Si. The integrated benthic dSi flux is substantially larger than the estimated requirements of primary producers and may explain the increasing trend of the water column dSi standing stock in almost all basins of the Baltic Sea that has been reported from 1994 (from 2000 for the BS) through 2023. Taken together, this would then imply that most of the enrichment of dSi in the Baltic basins is rSi of terrestrial origin which is dissolvable and recycled in sediments. It may also be that the export of dSi from the Baltic Sea to the Kattegat and beyond is larger than presently assumed. Although our calculation of the annual integrated benthic dSi load includes uncertainties, and certain seabed substrate types like sand have not been investigated often, and – especially for the GOB with relatively few sampling stations – may be an overestimation, particulate amorphous and/or lithogenic Si from land and freshwater is suggested to have a much more profound influence on the Si cycle of the Baltic Sea than previously known.

Results of this study should facilitate and stimulate the construction of a revised and more complete Si budget for the Baltic Sea. Future studies should in this regard include new estimates of burial of rSi (not only limited to bSi), of export of dSi and rSi to the Kattegat, of delivery of dSi and rSi via both rivers and groundwater, of reactivity of rSi in river water, groundwater and sediment, as well as continued observations of changes of the water column dSi standing stock.