Against the backdrop of increasing atmospheric carbon dioxide (CO2) concentrations, a fast development towards net-zero emissions is urgently needed1,2,3,4. For achieving this goal, it is necessary to develop negative emissions technologies to account for residual emissions that cannot be avoided5,6. Aside from established methods, such as carbon capture and storage, the enhanced weathering of rocks in the ocean is gaining increasingly more attention, as the process consumes CO2 via ocean alkalinity enhancement (OAE)7,8,9.

Apart from highly reactive hydroxides, calcite (a calcium carbonate morphotype) and dunite (a mafic silicate rock that mainly consists of olivine10), have been suggested as possible materials for OAE8,11,12,13,14,15,16. Both materials consume protons upon dissolution, increasing total alkalinity (TA) and subsequently leading to an uptake of CO2 in the water column and storage in the form of dissolved bicarbonate17. However, the abiotic dissolution of these materials in the water column is either chemically limited due to oversaturation in the case of calcium carbonates or, in the case of olivine and dunite, not feasible due to slow dissolution kinetics18,19.

A more promising method to achieve OAE, using calcite or dunite, is enhanced benthic weathering (EBW). In this approach, the minerals are distributed directly onto marine sediments, whereupon dissolution is furthered as they are either exposed to grain collision in high-energy environments, or the “benthic weathering engine”12,20,21. The latter describes a combination of benthic processes where particulate organic matter is respired by microorganisms to CO2, thus lowering the pH of the sedimentary porewaters. Bioturbation and bioirrigation enable oxygen to penetrate deeper into sediments, further stimulating CO2 production and acidification. Besides the microbes catalyzing these chemical reactions, invertebrates in- and digest sediment particles, where abrasion additionally contributes to mineral dissolution12,20. In the following, the combination of these processes will be referred to here as bioactivity.

The growing number of publications on marine carbon dioxide removal over the last few years clearly indicates that science is rapidly advancing in this field4. Still, the majority of these studies are based on models, which use mineral dissolution rates derived from laboratory experiments that do not reflect natural conditions, which can alter the kinetics of mineral dissolution tremendously 21,22,23,24. It is thus necessary to close the knowledge gap that exists between these assumed weathering rates and those that are actually valid in natural sediments, to make informed estimates on the efficiency and potential costs of EBW. This study is the first, to the best of our knowledge, to combine data from recent laboratory studies on EBW, where calcite and dunite were added to organic-rich sediments of the Baltic Sea23,24, with results from long-term experiments (September 2022–September 2023) in benthocosms. The latter are large plastic containers (~0.8 m2) that are partly filled with sediments, with constant flow-through of Baltic Sea bottom water (Fig. 1)25. The sediments were amended with calcite and dunite in triplicate25. Three additional benthocosms were left unamended to serve as controls. These experiments elucidate the actual weathering behavior of calcite and dunite under conditions as close to the natural system as possible. Weathering rates are subsequently implemented in a 4-box numerical model to simulate the addition of calcite and dunite to a 25 km2 application area in Eckernförde Bay (25 m water depth), a seasonally stratified and seasonally (~late July – ~ early October) hypoxic to anoxic inlet in the southwestern Baltic Sea (Fig. 1). The model simulates the seasonal hydrological and geochemical properties of this study area and calculates the additional CO2 uptake from the atmosphere by EBW as well as the amount of calcite and dunite that needs to be added to achieve this uptake. Additionally, for both materials, efficiency and cost calculations were carried out to understand the overall efficiency and cost per tonne of sequestered CO2 for each source material. Furthermore, the efficiencies and uptake potentials are then extrapolated to other marine systems with similar geochemical properties as the Baltic Sea to assess the global potential of EBW as a climate change mitigation measure.

Fig. 1: Map of the southwestern Baltic Sea. Blue shading indicates the bathymetry.
figure 1

Red areas denote the distribution of muddy sediments47. The adjacent sediment consists of gravel and sand with low organic carbon content. Yellow, red, and green lines describe the approximate spatial extent of the boxes considered for the box-model (Supplementary Fig. S1), where the yellow area is the area of application. The insert (upper right) shows a schematic sketch of the benthocosms used for the long-term experiment. Modified after Riebesell et al. 25.

Full size image

Results and discussion

Observed dissolution of calcite and dunite

The TA fluxes measured in mini-benthic chambers that were placed inside benthocosms (see methods and Riebesell et al.25), follow a clear pattern for the different treatments (Fig. 2). Within two weeks after the addition of alkaline minerals (0.44 g m−2 with a mean grain radius of 5.58 µm for dunite and 0.36 g m−2 with a mean grain radius of 19.39 µm for calcite), the TA fluxes increased in both treatment experiments from a background of 2.27 ± 0.00 mol m−2 a−1 to 10.60 ± 2.26 mol m−2 a−1 (calcite) and 6.06 ± 1.31 mol m−2 a−1 (dunite). In contrast to dunite, the dissolution of calcite is strongly dependent on the saturation state of the aqueous phase26,27. The fact that calcite dissolved despite oversaturated bottom waters after two weeks (Fig. 2b), indicates that the material must have been worked into the sediment, where undersaturated pore waters or the benthic weathering engine might have fostered the dissolution12,23. Over the first 224 days of the experiment (until May 2023) the fluxes closely followed the change in temperature under well-oxygenated conditions (Fig. 2a, c). When oxygen concentrations in bottom waters declined (Fig. 2c), this temperature correlation weakened as TA fluxes in calcite treatments stagnated or slightly decreased. This is likely a result of reduced bioactivity and less oxidation in upper sediments, both of which would lead to lower dissolution rates12,28. The undersaturation in porewaters in the upper 6 cm of sediment, expressed as a term (Disspot) that indicates the extent to which calcite can dissolve in porewaters23 (see “Methods”), appears to be less critical over the same period (Fig. 2d). This is emphasized by the fact that in both other experiments (dunite treated and control), high DissPot values did not apparently induce measurable natural calcite dissolution. Furthermore, the dissolution rates calculated based on the TA fluxes are more temperature sensitive compared to previously published values for both materials27,29(see Supplementary Table S1). This suggests that the intra-sedimentary calcite dissolution is driven by a variety of oxygen- and temperature-dependent processes. Most likely, bioactivity, which depends on both factors, pushes abiotic mineral dissolution to values higher than ascribable to the effect of ambient pore water saturation state. In a previous study dissolution of calcite in organic rich sediments in the same study area was simulated with a reaction transport model that accounts for bioturbation, but not local effects such as digestion by invertebrates and undersaturation along the walls of worm barrows28. As a result, the dissolution rates and subsequent TA fluxes to bottom waters were ~4 fold smaller, than observed in the benthocosm experiment in this study. This clearly substantiates the principle working function of the “benthic weathering engine”12. However, the macrofaunal community might not represent the natural system, as firstly the water that transports larvae into the system was pumped from relatively shallow depth (~15 m) and secondly the benthocosms were operated on shore, resulting in lower pressure compared to the sediment source area in the Eckernförde Bay at 25 m water depth.

Fig. 2: Geochemical data obtained from the benthocosm experiments.
figure 2

a TA fluxes (in mol m−2 a−1, bars) and temperature (in °C, curve), bΩCal, c dissolved oxygen (O2 in μmol l−1) and d Disspot (in cm) values in each treatment as a function of time. Values at t = 0 denote data obtained directly before the addition of alkaline materials. “Cal” denotes the calcite treatment, “Dun” the dunite treatment and “C” the untreated controls. ΩCal is the saturation state of bottom waters with respect to calcite whereas Disspot indicates the extent of calcite undersaturation in porewaters23. Error bars are reported as SD of the three replicas for each treatment. Several TA fluxes in the control and dunite-treated experiments were below detection limit (2.27 mol m−2 a−1) and are therefore not shown. Sampling months are given at the top of the figure.

Full size image

The benthocosms were run ashore, with bottom waters pumped in from the Kiel Fjord (Fig. 1). This inlet is shallower than Eckernförde Bay, which is why the bottom water oxygenation state is more sensitive to wind driven up- or down-welling and mixing. Hence, the bottom waters in the benthocosms did not follow a full annual cycle, but were oxic at the beginning of the experiment, until the hypoxic phase began after ~280 days. Consequently, the material had already been worked into the sediment and fluxes in both treatments were not as high as observed in previous studies under low-oxygen conditions, where the material remained on top of the sediment24. The TA fluxes induced by calcite addition during the oxic period were ~2 fold higher compared to dunite-treated experiments (Fig. 2a). For the latter, through February and March, TA fluxes were below the detection limit. This was also observed in the control experiments from December to May Overall, the TA fluxes induced by the addition of calcite and dunite are in good agreement with findings of prior studies on weathering on sediments from the same study area23,24, given the difference in reactive surfaces of the added material. This, and the fact that the treatment with calcite and dunite was the only difference with respect to the controls, allows for the assumption that the increased TA flux compared to the control experiment resulted from the dissolution of the added alkaline minerals—despite the fact that corresponding Ca2+ and H4SiO4 fluxes (from calcite and dunite weathering) could not be resolved against the high background concentrations.

Modeled impact of the reactive surface area on the CO2 uptake efficiency

For the assessment of additional CO2 uptake induced by calcite and dunite addition to a 25 km2 application site in Eckernförde Bay, the dissolution rates based on benthic chamber measurements in the benthocosm experiments and previous incubation experiments were used as inputs to a box model23,24. The model was run for 35 years with an annual mineral addition at the beginning of seasonal hypoxic to anoxic phase, so that the added mineral would remain exposed to corrosive bottom waters as long as possible. The model considers two layers in the water column for Eckernförde Bay and the adjacent Kiel Bight (Fig. 1, Supplementary Fig. S1) as well as dissolution of added alkaline materials on and within the sediment. Note that the used model is a strong simplification of the natural system. Detailed hydrological features such as mixing rates and water exchange with the adjacent sea lack precision, and the results must therefore be interpreted with care. Still, they can be assumed to be a reasonable approximation.

The model shows that after 30 years, for both dunite and calcite addition, higher average annual CO2 uptake can be achieved via either increased addition of material or a higher reactive surface, implying a smaller grain size (Fig. 3a, b). Note that for dunite the geometrical surface (AGEO) was used, which is the theoretical surface if all grains were perfectly spherical, while for calcite the measured Brunauer–Emmett–Teller (BET)30 surface is reported, which is the actual surface of non-spherical grains (ABET). The latter was ~4.55 times larger than the AGEO for calcite grains. The ratio of CO2 uptake to alkaline mineral addition is non-linear. For dunite, the additional uptake induced by higher additions increases continuously for all grain sizes (Fig. 3a). For calcite, very low additions (up to ~2500 t yr−1) lead to a fast-increasing uptake, regardless of the grain size. Afterward, the additional uptake per tonne of added material is smaller, linear and more dependent on the reactive surface area (Fig. 3b). The uptake efficiency, calculated as the total uptake in tonnes of CO2 per tonne of added material indicates that for both materials a higher mean annual uptake is achieved at the expense of lower uptake efficiency (Fig. 3c, d). This further requires a higher usage of material, that is, the amount of material needed to sequester 1 t of CO2.

Fig. 3: Modeled CO2 uptake and efficiency.
figure 3

Mean additional CO2 uptake (above natural background) as a function of mineral addition in tonnes per year (a,b) and the CO2 uptake efficiency (c, d) for 30 years of dunite addition (a,c) and calcite addition (b,d) over an application area of 25 km2. The width of the lines indicates the geometrical (dunite) and BET (calcite) surface area, respectively. Blue-filled symbols (a,b) indicate the amount of added material, up to which experimentally derived dissolution kinetics are valid. Red circles indicate the maximum addition that can be applied before reprecipitation of calcite might compromise CO2 uptake (see text). Color codes of the symbols in (c) and (d) represent the amount of material added in each model run.

Full size image

Whilst the dependence of the efficiency of dunite weathering on the reactive surface and the amount of added material follows the expected trends, calcite dissolution shows a surprising behavior. The CO2 uptake efficiency does not change significantly for low additions, and is also not affected by the BET surface of the added material up to an addition of ~1 mol m−2 yr−1. This pattern can be explained by the fact that in the model, calcite is added at the beginning of the hypoxic phase, which occurs regularly in study area due to seasonal stratification of the water body31,32,33,34. During these phases, bioturbation is practically absent, as laboratory incubations under such conditions have shown24. Thus, the added substrate is unlikely to be directly worked into the sediment, as observed in the benthocosm experiment, but instead remains on the surface and dissolves at a rate depending on the bottom water saturation state24 (Eq. 3, Methods). This way, for low additions in the model, the entire amount of added calcite has dissolved in corrosive bottom waters, before the beginning of the oxic phase and a subsequent downward transport of the substrate into the sediment by bioturbation. If more material is added with a large reactive surface, calcite dissolves quickly, and bottom waters in the model become oversaturated after a short time, which inhibits further dissolution. For smaller BET surface areas, the bottom waters remain undersaturated for a longer period of time, but the grains do not dissolve entirely before the end of the anoxic phase due the small surface area. Hence, the amount of calcite that dissolves before it is worked into the sediment depends on the amount of added mineral, rather than the physical properties of the added material; at least for the grain size range considered in this study.

These findings imply that the grain size range of added calcite material can be rather large in practice, which facilitates easier production and lowers costs. The precise amount of calcite that can dissolve in bottom waters before they reach supersaturation is dependent on the remineralization rate of particulate organic carbon in bottom waters and surface sediments, which is the major source for CO2 in bottom waters, and thus on the rate of primary production. Moreover, the water exchange rate with the adjacent sea and the upwelling rate determine the dilution of weathering products in bottom waters, and consequently the saturation state and the amount of calcite that can dissolve. In the model, mixing in the two larger Kiel Bight boxes reintroduces alkaline water into the smaller bottom water box (Supplementary Fig. S1). This can lead to an overestimation of ΩCal values and thus diminish calcite dissolution. For this reason, a hydrologically more accurate model will be needed to assess the precise uptake threshold up to which EBW in the study area is highly efficient.

After the uptake threshold of ~1000 t yr−1, which equals a yearly addition of 0.75–1 mol m−2 over the entire study area, dissolution takes place mainly within the sediment, as incremently more added calcite is worked into the sediment before it can fully dissolve in undersaturated bottom waters. Since the model does not consider changes in the saturation of pore waters induced by mineral dissolution, the rate does not change with larger calcite content in the sediment and thus the efficiency remains almost constant. The dissolution rates in sediments, however, may be higher if microbial activity leads to higher Disspot values induced by extremely low pH values as observed in previous studies23,35,36.

The initial addition of calcite and dunite in the benthocosm experiment was 22 mol m−2. After 1 year, the material was mixed into the upper ~3.3 cm of sediment (Supplementary Fig. S2). Hence, the dissolution rates calculated from the benthocosm experiments are valid for an amount of ~22 mol of calcite in the upper 3.3 cm of sediment, which subsequently constitutes a threshold up to which calculations by the model can be assumed as reliable (blue symbols, Fig. 3a, b). A second threshold is the maximum saturation state with respect to calcite, above which the spontaneous precipitation of CaCO3 can occur diminishing the uptake efficiency due to the release of CO2. This threshold depends on saturation state, as well as on equilibration with CO2 and the used material15,22,37. Since carbonate precipitation was observed at ΩCal = ~6–8 in previous studies22,38 in the presence of olivine, a value of 6.1 in the water column was defined as a threshold for both material additions, despite the fact that recent studies about the stability of alkalinity observed precipitation of CaCO3 only at values of ΩCal > 1237,39. Red circles in Fig. 3a, b define the amount of added mineral that lead to a transgression of this threshold. A third threshold, that could not be defined in the frame of this study is the effect of additionality, e.g. whether the addition of alkaline minerals leads to a reduced natural TA release40. This will demand dedicated studies to ensure an effective application of EBW. This is also true for the thresholds emerging from a potential ecosystem impact, which will need to be thoroughly studied and understood, before any possible large-scale application scenario in the field.

Identification of a candidate material for EBW based on cost estimates

For the estimation of overall efficiency and costs for the modeled application scenario, the single steps from the production of a respective calcite and dunite material to the actual spreading process in the application area were analyzed for their expenditure and their respective CO2-footprint (Supplementary Tables S8–S10). As to be expected, the costs per tonne of sequestered CO2 via addition of calcite and dunite (Fig. 4) follow the inverted uptake efficiency (Fig. 3c, d), which explains the constant costs for calcite application up to ~1000 t annual uptake. For dunite, the lower energy demand for the production of coarser grains did not compensate for the lower uptake efficiency, which is why the total costs per tonne of CO2 increase with smaller grain size over the grain size range assumed for the model (~5 µm–~30 µm). Smaller grain sizes were not considered as these might remain in suspension and would thus not be suitable for EBW19. For both materials the costs for production and grinding exceeded the transport costs by a factor of ~5 for the coarsest grainsize (Supplementary Tables S8–S10), the implication being that proximity to the application site or the option for long-distance transport via ship are key to extrapolate these pricing estimates to other application sites.

Fig. 4: Costs for CO2 sequestration via enhanced benthic weathering.
figure 4

Mean costs including production, comminution, transport, and spreading per tonne of sequestered CO2 for dunite (a) and calcite (b) for a 30-year application period. Color codes indicate the annual addition of material. Line width indicates the geometrical surface (olivine) and the corresponding BET surface (calcite), respectively. Note different range of x-axis in panels a and b.

Full size image

In sum, the low CO2 uptake efficiency, combined with the high energy demand for grinding exclude dunite as a potential material for efficient EBW in the investigated study area.

For calcite, these constraints do not apply. The softer material requires less energy for comminution, which lowers the costs and CO2 footprint compared to dunite. Combined with the higher uptake efficiency (Fig. 3c) this leads to prices between 81.6–99.5 € tCO2−1 for a total CO2 uptake of up to ~1000 t over the entire application area (Fig. 4b). For higher, but less efficient total uptake, the costs vary between 172€ tCO2−1 (ABET = 3.0 m2 g−1) and 462€ tCO2−1 (ABET = 0.5 m2 g−1). For this uptake range, the lower uptake efficiency of the coarser material is—as for dunite—not outbalanced by the lower production costs and CO2 footprint. Furthermore, the additional uptake for all treatments increased over time for all treatments (Supplementary Figs. S3, S4). Hence, an application over a longer time period would increase the overall uptake efficiency and lower the costs per tonne of sequestered CO2., Especially below the uptake of 1000 t, the prices for calcite based EBW would be highly competitive when compared to other technical CDR methods such as direct air capture (DACCS) and bioenergy with carbon capture and storage (BECCS)2,41. Still, the exact potential for EBW and the precise threshold for high uptake efficiencies need further research before a final comparison can be made.

Enhanced benthic weathering in a global context

The results of this study clearly show that mineral dissolution can lead to considerable marine CDR, as long as the alkalized waters are in contact with the atmosphere to actually take up CO2. These new findings have implications for the perception of EBW as a climate change mitigation measure. The Baltic Sea alone comprises 133,000 km2 of muddy sediments42. Assuming a similar weathering behavior for calcite as in Eckernförde Bay, this would imply an uptake between ~4.5 and 30 Mt yr−1. The latter, though, with lower efficiency.

Despite the fact that the study area is a rather small inlet, the chemical features bear some similarities to large-scale upwelling systems where nutrient rich waters foster high rates of primary production and organic carbon rain rates. The export of organic matter in such region leads to high organic carbon contents in sediments that are similar to those in the study area Eckernförde Bay20,23,43,44,45. The degradation of this organic carbon combined with the upwelling of pCO2-rich waters alters the bottom water chemistry in a way that makes these large systems comparable to the inlet investigated in this study. Thus, the presented results may well be extrapolated to a globally relevant scale as the total area of eastern boundary upwelling systems alone accounts for ~1.04\(\cdot\)106 km2 46. When including the highly productive East China Sea continental margin, the potential application area could be up to ~2.0\(\cdot\)106 km2 plus several smaller upwelling regions46. Considering a higher, but less efficient uptake of 320 t km−2 yr−1 and also including smaller potential application sites, annual global potential CO2 uptake might reach several gigatonnes. Although it is unlikely that EBW will reach a potential to serve as a sole sink for global non-abatable CO2 emissions, it may well be an efficient way for countries bordering upwelling systems or organic-rich shelf sediments to offset non-point-sourced residual emissions or earn carbon credits.

Methods

Experimental set-up

Between 22 and 24 February 2022, a total of nine polyethylene boxes with a base area of 111 cm × 71 cm were filled with ~20 cm of organic-rich sediment from Boknis Eck (Eckernförde Bay, south-western Baltic Sea) using a sediment grab. Beforehand, a side wall of each benthocosms was drilled with sealable holes to insert rhizons for porewater extraction and subsequent analysis. The depth of the supernatant bottom water in each box was ~30 cm. Subsequent to recovering the sediment, the boxes were placed on a quay in Kiel Fjord adjacent to GEOMAR. In groups of three, the boxes were connected to a header tank that ensured constant flow rates. From these header tanks, ~10 m3 d−1 of bottom water (~15 m depth) was pumped into the nine benthocosms via a dispenser to avoid resuspension of the fine-grained sediment (Fig. 1). The benthocosms were then left for 6 months to equilibrate and to allow for natural chemical pore water gradients to develop.

On 27 September 2022, 22 mol m−2 of limestone (99.8% calcite, German Limestone Association) and dunite (90% olivine, SibelcoTM) powder was distributed evenly on the sediment surface of six benthocosms with three replicates for each treatment. The other three benthocosms were left as unamended controls. An overview of the precise treatments is given in Supplementary Table S11.

Sampling and measurements

Benthic fluxes in each benthocosm were measured with sealed polycarbonate mini-benthic chambers with a base area of 100 cm2 and a total volume of 400 cm3. The chambers were first flipped upside down to fill with bottom water and then placed carefully on the sediment. A ~ 0.5 cm rim around the chamber penetrated into the sediment to assure that the chamber water was isolated from the surrounding water. Directly before deployment, a bottom water sample was taken to determine initial conditions. After 3 h, a sample from inside the chamber was taken via a PVC tube that was attached to the chamber in the middle of the side wall. The first 15 ml of water was discarded. Afterward, a 20 ml sample was filtered using 0.2 µm regenerated cellulose acetate syringe filters and refrigerated in polyethylene (PE) vials without headspace for a maximum of 2 days for further measurements. Total alkalinity (TA) was determined via titration using 0.02 N HCl48 as described in previous studies49. A 3 ml aliquot was acidified with 30 µl suprapure HNO3 for cation analysis via inductively coupled plasma optical emission spectrometry (ICP-OES, Varian 720-ES).

Before the addition of substrates, pore water samples were taken using 10 cm long 0.15 µm standard rhizon samplers attached to 20 ml PE syringes. After recovery of the pore waters, the syringes were sealed, refrigerated, and transported to the laboratory for further processing. In an N2-filled glove bag, the samples were transferred to 20 ml PE vials. TA and cations were determined the same way as for the benthic chamber samples.

Micro-profiling of oxygen and pH in each benthocosm was conducted with UnisenseTM micro-sensors Ox-100 and pH-100, mounted on a motorized UnisenseTM micromanipulator. To maximize measuring resolution for oxygen profiles the step size was adjusted to 100 µm. Oxygen sensors were calibrated according to UnisenseTM guidelines. For pH profiles, the resolution was lowered to 300 µm to minimize profiling time. In addition to technical buffers (pH 4 and 7), the pH sensors were calibrated with TRIS-buffers50,51. pH values in this study are reported on the total scale.

Samples for analysis of solid components over the upper 6 cm were taken using cut-off syringes (2.5 cm diameter). The samples were obtained directly after addition of the alkaline substrates to determine conditions at the starting point of the experiment and at the end of the experiment. Sediment samples were freeze-dried, ground, and analyzed via flash combustion using the EuroEA 3000 element analyzer (EuroVector, Pavia, Italy) to determine total carbon (TC), total organic carbon (TOC), total nitrogen (TN), and total sulfur (TS). The TIC content was calculated by subtracting the TOC value from the TC measurement. To assess the accuracy of the analytical method, method blanks and two reference standards were employed: 2.5-Bis(5-tert-butyl-2-benzo-oxazol-2-yl)thiophene (HEKAteckTM), and an internal sediment standard.

Flux calculations from the mini-benthic chambers

Fluxes of solutes, Fsol, were determined from the initial and final concentrations with the following equation:

$${{{{\rm{F}}}}}_{{{{\rm{sol}}}}}=\frac{{{{\rm{d}}}}{{{{\rm{C}}}}}_{{{{\rm{sol}}}}}}{{{{\rm{dt}}}}}\times {{{{\rm{H}}}}}_{{{{\rm{Chamber}}}}}$$
(1)

where dCsol/dt is concentration change over time and HChamber is the height of the overlying water.

Calculation of Disspot values

To assess the impact on porewater saturation state on calcite dissolution, the Disspot values (cm) were calculated following Fuhr et al. (2023) by integrating ΩCal values over the region where pore waters were undersaturated as

$${{{{\rm{Diss}}}}}_{{{\rm{pot}}}}=\int ^{{{{\rm{z}}}}({{\Omega }} < 1)}_{{{{\rm{z}}}}({{\Omega }} > 1)}(1-{{\Omega }}_{{{\rm{cal}}}}({{{\rm{z}}}})) \, \, {{{\rm{dz}}}}$$
(2)

where Ωcal(z) is the saturation state at a certain depth, z (cm), in the sediment given. z(Ω < 1) and z(Ω > 1) denote depths where Ωcal values fall below 1 and above 1, respectively.

ΩCal values were calculated based on the pH profiles at the respective time step and the porewater TA as well as pore water Ca concentrations sampled at the beginning of the experiment, based on seawater constants following Zeebe and Wolf-Gladrow52. To reach the resolution of the pH profile (300 µm), TA and Ca concentrations were interpolated linearly. These data are available online (see data availability statement) and in Supplementary Tables S12 and S13.

Box model set-up

To assess the impact of calcite and dunite addition to an application area of 25 km2, a numerical box model was set up using WolframTM MathematicaTM v12. The model considers four boxes, which represent a two-layered water body including the mineral application area in Eckernförde Bay and the adjacent Kiel Bight (Fig. 1, Supplementary Fig. S1). Box sizes were chosen accordingly. The water depth was defined as 25 m. A boundary between the upper and lower boxes was defined at 15 m depth as this is the depth where the halocline usually occurs during summer33. The Kiel Bight box was considered to be 100 times larger than the Eckernförde Bay box. The area of the bottom box of Kiel Bight was 34% smaller than the surface box in order to represent the area of Kiel Bight that is covered with sandy sediments47,53,54. Thus, the benthic fluxes contributing from sandy sediments55 were not considered for the model. The model considers primary production (PP) that transports reduced dissolved inorganic carbon (DIC) in the surface layer and exports it to the lower box as particulate organic carbon (POC), where POC is then remineralized to DIC. Furthermore, the model accounts for water mass mixing between the upper and lower boxes, as well as water mass exchange between the two surface boxes. To account for a quasi-upwelling that is induced by a saltwater wedge that invades and replaces bottom water in spring34, the model includes an upwelling term (Supplementary Fig. S1). The upwelling flux enters the Kiel Bight bottom box and flushes the entire system before it exits via the Kiel Bight surface box. Hence, the Kiel Bight bottom box serves to pre-form the water that circulates through Eckernförde Bay. Accordingly, the model was tuned to produce naturally prevailing conditions in both bottom boxes31,33. Hydrological parameters were estimated based on a 3-D numerical model34. Primary production was qualitatively assessed based on chlorophyll measurements33,56, normalized, and then tuned in such a way that the bottom water properties matched observed values31.

The upper two boxes can exchange CO2 with the atmosphere, which allows quantifying the natural (without addition of alkaline materials) as well as the enhanced (induced by addition of alkaline materials) CO2 flux across the sea-atmosphere interface. The lower Eckernförde Bay box includes DIC and TA fluxes from the sediment that derive from natural CaCO3 dissolution and degradation of sedimentary POC. Natural background fluxes were determined via the unamended benthocosm experiment. The sediment model consists of three equally sized layers that describe dissolution of added minerals within the upper 10 cm of the sediment (not shown in the schematic in Supplementary Fig. S1). Due to the long simulation period of 35 years, the TA and DIC fluxes from the sediment were assumed to directly enter the bottom water without delay via molecular diffusion.

The model solves the carbonate system in each of the boxes based on DIC and TA via a set of coupled mass balance equations (ordinary differential equations) using finite differences and the method-of-lines approach as implemented in the partial differential equation solver. For simplicity, other chemical processes such as oxidation and higher-level food chains are not considered in the model. Initial conditions for carbonate system properties were inferred from long-term observations at Boknis Eck31,33. Salinity, temperature, and oxygen concentrations were implemented as time-dependent parameters for upper and lower boxes. Salinity and temperature were used to calculate the solubility and dissociation constants for the different carbonate species52. For bottom water oxygen concentrations, the observed long-term decreasing trend of −1.4 µmol l−1 was considered33.

To ensure full equilibration, the model was spun-up for five years before calcite or dunite were added in July, at the beginning of the anoxic phase, and then every July for another 30 years (Supplementary Table S2). The added materials were assumed to stay on top of the sediment until bottom water oxygen concentrations reached 100 µmol l−1, which is in line with observations from the benthocosm experiment and previous studies23,24. Above this threshold, the material is allowed to be transported into the upper 3.3 cm of sediment by bioturbation. This depth stems from previous reports32 including own observations23. The transport was parameterized using a first-order rate constant that was chosen such that 80% of the material was worked into the sediment after ~3 weeks and 99% after 7 weeks. This is based on observations from the benthocosm experiment and previous studies23,24. For the transport into the deeper sediment (i.e., the third layer) the model assumed burial using a sedimentation rate of ~3.3 mm yr−1 57,58. Below 10 cm, no further dissolution was assumed as porewaters below this depth are oversaturated with respect to calcite23,24 and bioturbation does usually not reach deeper than 10 cm59. The combination of sedimentation rate and depth of bioturbation implies that the residence time of added material in the top 10 cm is ~30 years, which defines the chosen model simulation period.

For each box, an ordinary differential equation was formulated as

$$\frac{{{{\rm{d}}}}{{{{\rm{C}}}}}_{{{{\rm{i}}}}}}{{{{\rm{dt}}}}}={{{{\rm{T}}}}}_{{{{\rm{upw}}}},{{{\rm{i}}}}}+{{{{\rm{T}}}}}_{{{{\rm{mix}}}},{{{\rm{i}}}}}+{{{{\rm{T}}}}}_{{{{\rm{ex}}}},{{{\rm{i}}}}}+{{{{\rm{F}}}}}_{{{{\rm{CO}}}}2}+{\sum }_{{{{\rm{i}}}}}{{{{\rm{R}}}}}_{{{{\rm{i}}}}}$$
(3)

where i is DIC or TA and Tupw, Tmix and Tex describe the water transport between adjacent boxes due to upwelling, mixing and lateral exchange, respectively (Supplementary Fig. S1) and R defines changes due to formation/remineralization of POC and dissolution of minerals. FCO2 describes the CO2 exchange with the atmosphere for DIC only.

Dissolution of calcite

Calcite dissolution in undersaturated bottom waters (ΩCal< 1) was simulated as follows:

$${{{\rm{CaC}}}}{{{{\rm{O}}}}}_{3}+{{{{\rm{H}}}}}_{2}{{{\rm{O}}}}+{{{\rm{C}}}}{{{{\rm{O}}}}}_{2}\to {{{{\rm{Ca}}}}}^{2+}+2{{{\rm{HC}}}}{{{{\rm{O}}}}}_{3}^{-}$$
(4)

The rate for this dissolution process, \({{{{\rm{r}}}}}_{{{{\rm{CalDiss}}}}}\), was defined as:

$${{{{\rm{r}}}}}_{{{{\rm{CalDiss}}}}}({{{\rm{t}}}})={{{{\rm{Cal}}}}}_{{{{\rm{top}}}}}({{{\rm{t}}}})\times {{{{\rm{A}}}}}_{{{{\rm{BET}}}}}\times {{{{\rm{k}}}}}_{{{{\rm{CalDiss}}}}}\times {(1-{{{\Omega }}}_{{{{\rm{cal}}}}})}^{{{{{\rm{n}}}}}_{{{{\rm{CalDiss}}}}}}$$
(5)

where Caltop is the amount of calcite on the sediment in mol, ABET is the reactive surface of the calcite grains in cm mol−1, kCalDiss is a first-order reaction constant in mol cm−1 d−1, ΩCal is the saturation state with respect to calcite and nCalDiss is the reaction order. kCalDiss and nCalDiss were dependent on ΩCal following previous studies60. ΩCal was defined as

$${{{\Omega }}}_{{{{\rm{Cal}}}}}=\frac{\left[{{{{\rm{Ca}}}}}^{2+}\right]\times \left[{{{\rm{C}}}}{{{{\rm{O}}}}}_{3}^{2-}\right]}{{{{{\rm{K}}}}}_{{{{\rm{sp}}}}}}$$
(6)

[Ca2+], \([{{{{\rm{CO}}}}}_{3}^{2-}]\) and Ksp are the porewater concentrations of calcium and carbonate ions, and the solubility product of calcite, respectively52.

Calcite dissolution in the three sediment boxes (rCallDissSed) was estimated based on the enhanced TA fluxes measured during the benthocosm experiments (Supplementary Table S4) following:

$${{{{\rm{r}}}}}_{{{{\rm{CallDissSed}}}}}({{{\rm{t}}}},{{{\rm{z}}}})={{{{{\rm{k}}}}}_{{{{\rm{depth}}}}}({{{\rm{z}}}})\times {{{\rm{Cal}}}}}_{{{{\rm{Sed}}}}\; {{{\rm{i}}}}}({{{\rm{t}}}})\times {{{{\rm{A}}}}}_{{{{\rm{BET}}}}}\times {{{{\rm{rT}}}}}_{{{{\rm{Cal}}}}}({{{\rm{t}}}})\times {{{{\rm{r}}}}}_{{{{\rm{mean}}}}}$$
(7)

where z denotes the average depth of each of the three sediment boxes, CalSed i (t) is the total amount of calcite in the respective sediment box in mol, ABET is the reactive surface, rTCal (t) is the temperature sensitivity (Supplementary Eq. S26) and rmean is the mean calculated dissolution rate (Supplementary Eq. S25). Since the added material only mixes into the upper three centimeters of sediment after one year in the benthocosms, the bulk dissolution kinetics for sedimentary dissolution of calcite is only valid for the upper sediment box. For deeper sediment layers, an exponential decrease of dissolution was assumed and dissolution was calculated by multiplying Eq. 5 with a pre-factor kdepth(z):

$${{{{\rm{k}}}}}_{{{{\rm{depth}}}}}\left({{{\rm{z}}}}\right)={{{{\rm{e}}}}}^{\left(-0.538{{{\rm{z}}}}+0.778\right)}$$
(8)

where z is the depth in the sediment. kdepth was chosen in a way that it is equal to 1 when averaged over the upper 3.3 cm and 0.01 at 10 cm depth, which constitutes the depth at which no further bioturbation and thus dissolution is assumed.

Olivine dissolution

Olivine dissolution was defined as:

$$\begin{array}{c}{\left({{{{\rm{Mg}}}}}_{0.93},{{{{\rm{Fe}}}}}_{0.07}\right)}_{2}{{{\rm{Si}}}}{{{{\rm{O}}}}}_{4}+{4.07{{{\rm{H}}}}}_{2}{{{\rm{O}}}}+0.035{{{{\rm{O}}}}}_{2}+3.72{{{\rm{C}}}}{{{{\rm{O}}}}}_{2}\to \\ 1.86{{{{\rm{Mg}}}}}^{2+}+{{{{\rm{H}}}}}_{4}{{{{\rm{SiO}}}}}_{4}+0.14{{{\rm{Fe}}}}{({{{\rm{OH}}}})}_{3}+3.72{{{{\rm{HC}}}}{{{{\rm{O}}}}}_{3}}^{-}\end{array}$$
(9)

Assuming that olivine consists of 93% forsterite and 7% fayalite23,24, the dissolution rate can be parameterized following literature values27 and adjusted based on sediment incubation experiments for anoxic to hypoxic conditions24 as:

$${{\mathrm{log}}}\left({{{{\rm{r}}}}}_{{{{\rm{ol}}}}}({{{\rm{t}}}})\right)=4.07-0.256\times {{{\rm{pH}}}}-\frac{3465}{{{{\rm{T}}}}}$$
(10)

where rol is the dissolution rate given in mol m−2 s−1, T is temperature in K and pH is the pH in bottom waters. Dissolution in the sediment was assumed following the same approach as for calcite, with mean dissolution rates calculated from benthocosm experiments with dunite addition.

pH simulation

TA and DIC were defined as follows:

$${{{\rm{TA}}}}=[{{{\rm{HC}}}}{{{{\rm{O}}}}}_{3}^{-}]+2[{{{{\rm{CO}}}}}_{3}^{2-}]+[{{{\rm{B}}}}{({{{\rm{OH}}}})}_{4}^{-}]+[{{{{\rm{OH}}}}}^{-}]-[{{{{\rm{H}}}}}^{+}]$$
(11)
$${{{\rm{DIC}}}}=[{{{{\rm{CO}}}}}_{2}]+[{{{\rm{HC}}}}{{{{\rm{O}}}}}_{3}^{-}]+[{{{{\rm{CO}}}}}_{3}^{2-}]$$
(12)

Each term on the right-hand side of Eq. (11) was replaced for an expression based on the total concentration of each weak acid, leading to:

$${{{\rm{TA}}}} = \frac{{{{{\rm{K}}}}}_{1} \cdot \left[{{{{\rm{H}}}}}^{+}\right] \cdot [{{{\rm{DIC}}}}]}{{\left[{{{{\rm{H}}}}}^{+}\right]}^{2}+\left[{{{{\rm{H}}}}}^{+}\right] \cdot {{{{\rm{K}}}}}_{1}+{{{{\rm{K}}}}}_{1} \cdot {{{{\rm{K}}}}}_{2}}+\frac{2 \cdot {{{{\rm{K}}}}}_{1} \cdot {{{{\rm{K}}}}}_{2} \cdot [{{{\rm{DIC}}}}]}{{\left[{{{{\rm{H}}}}}^{+}\right]}^{2}+\left[{{{{\rm{H}}}}}^{+}\right] \cdot {{{{\rm{K}}}}}_{1}+{{{{\rm{K}}}}}_{1} \cdot {{{{\rm{K}}}}}_{2}} \\ +\frac{{{{{\rm{K}}}}}_{{{{\rm{B}}}}} \cdot [{{{\rm{TB}}}}]}{{{{{\rm{K}}}}}_{{{{\rm{B}}}}}+\left[{{{{\rm{H}}}}}^{+}\right]}+\frac{{{{{\rm{K}}}}}_{{{{\rm{W}}}}}}{\left[{{{{\rm{H}}}}}^{+}\right]}-\left[{{{{\rm{H}}}}}^{+}\right]$$
(13)

The calculation of protons was calculated at each time step by solving Eq. 13 as described by Zeebe and Wolf-Gladrow16. Stoichiometric equilibrium coefficients were calculated based on temperature and salinity following Zeebe and Wolf-Gladrow52.

The carbonate ion concentration in each box Eq. (6) was then calculated as:

$${{{\rm{C}}}}{{{{\rm{O}}}}}_{3}^{2-}=\frac{{{{{\rm{K}}}}}_{1} \cdot {{{{\rm{K}}}}}_{2}}{{\left[{{{{\rm{H}}}}}^{+}\right]}^{2}+\left[{{{{\rm{H}}}}}^{+}\right] \cdot {{{{\rm{K}}}}}_{1}+{{{{\rm{K}}}}}_{1} \cdot {{{{\rm{K}}}}}_{2}} \cdot {{{\rm{DIC}}}}$$
(14)

and the partial pressure of CO2 we calculated as:

$${{{{\rm{pCO}}}}}_{2}=\frac{\left[{{{{\rm{CO}}}}}_{2}\right]}{{{{{\rm{K}}}}}_{0}}$$
(15)

where K0 is the solubility constant for CO2, and where

$$\left[{{{\rm{C}}}}{{{{\rm{O}}}}}_{2}\right]=\frac{\left[{{{\rm{HC}}}}{{{{\rm{O}}}}}_{3}^{-}\right] \cdot \left[{{{{\rm{H}}}}}^{+}\right]}{{{{{\rm{K}}}}}_{1}}$$
(16)

Initial conditions for TA and DIC were chosen based on values reported by Melzner et al.31 and reported in Supplementary Table 3.

Additional CO2 uptake induced by addition of alkaline substrates

The model considers equilibration of CO2 with the atmosphere to estimate the additional uptake after addition of alkaline substrates. For this purpose, the model was run for 35 years without addition of alkaline materials to obtain the natural baseline. Subsequently, the natural CO2 uptake was subtracted from the uptake with addition of calcite and dunite, respectively.

The air-sea gas flux (Eq. 3) from surface box i, \({{{{\rm{F}}}}}_{{{{\rm{CO}}}}2,{{{\rm{i}}}}}\), was defined as

$${{{{\rm{F}}}}}_{{{{\rm{CO}}}}2,{{{\rm{i}}}}}={{{{\rm{k}}}}}_{{{{\rm{CO}}}}2}\times {{{{\rm{K}}}}}_{0}\times ({{{\rm{pC}}}}{{{{\rm{O}}}}}_{2{{{\rm{atm}}}}}-{{{\rm{pC}}}}{{{{\rm{O}}}}}_{2,{{{\rm{i}}}}})$$
(17)

where kCO2 is the gas transfer velocity for CO2, and pCO2,atm and pCO2,i are the partial pressures of CO2 in the atmosphere and the surface water of box i, respectively. A value of 419 µatm was used as the initial condition for pCO2,atm. An annual increase of 2.4 µatm was assumed61.

The gas transfer velocity normalized for CO2 at 20 °C and zero salinity (k(600), cmh−1) was determined via in-situ measurements in coastal Baltic Sea in June and September 2018 (Dobashi and Ho, in preparation) and parameterized as:

$${{{\rm{k}}}}(600)=1.26\times {{{{\rm{u}}}}}_{10{{{\rm{m}}}}}({{{\rm{for}}}}\,{{{{\rm{u}}}}}_{10{{{\rm{m}}}}} < {7{{{\rm{ms}}}}}^{-1})$$
(18)

and

$${{{\rm{k}}}}\left(600\right)=8.83\times({{{\rm{for}}}}\,{{{{\rm{u}}}}}_{10{{{\rm{m}}}}}\ge {7{{{\rm{ms}}}}}^{-1})$$
(19)

where u10m is the wind speed at 10 m height62. Subsequently, the k(600) values were corrected for in-situ salinity and temperature as

$${{{{\rm{k}}}}}_{{{{\rm{CO}}}}2}={{{\rm{k}}}}(600){\times \left(\frac{{{{{\rm{Sc}}}}}_{{{{\rm{CO}}}}2}}{600}\right)}^{-\frac{1}{2}}$$
(20)

where ScCO2 is the Schmidt number (i.e., kinematic viscosity of water divided by diffusion

coefficient of gas in water) for CO2, which was calculated following Dobiashi & Ho63

$${{{{\rm{Sc}}}}}_{{{{\rm{CO}}}}2} = 1914.2+5.6330\times {{{{\rm{Sal}}}}}_{{{{\rm{up}}}}}-\left(123.18+0.32481\times {{{{\rm{Sal}}}}}_{{{{\rm{up}}}}}\right)\times {{{{\rm{T}}}}}_{{{{\rm{up}}}}}\\ +\left(4.2040+0.010208\times {{{{\rm{Sal}}}}}_{{{{\rm{up}}}}}\right)\times {{{{\rm{T}}}}}_{{{{\rm{up}}}}}^{2}\\ -(0.079322+0.00018296\times {{{{\rm{Sal}}}}}_{{{{\rm{up}}}}})\times {{{{\rm{T}}}}}_{{{{\rm{up}}}}}^{3}\\ +(0.00062463+1.3992\times {10}^{-6}* {{{{\rm{Sal}}}}}_{{{{\rm{up}}}}})\times {{{{\rm{T}}}}}_{{{{\rm{up}}}}}^{4}$$
(21)

where Salup is salinity in surface waters and Tup is surface temperature in °C.

A detailed description of the model, including all model parameters and boundary conditions is provided in the supplement.

Efficiency and cost assessment

CO2 emissions from production, transport and spreading of calcite and dunite were calculated. For dunite, this assumes production in Åarheim, Norway, a 1400 km transport by ship to Kiel harbor, milling (at a theoretical facility in Kiel harbor), and spreading by ship in Eckernförde Bay (~150 km by a small vessel). For calcite, we considered production and milling in a quarry 300 km south of Kiel, transport to Kiel by train, and subsequent spreading in the application area. Further impacts from mining, such as altered land use were not considered.

Based on Strefler et al.64 the production costs per tonne of dunite were calculated excluding grinding and transport. Additionally, energy demand (Egr) in kWh for grinding based on the grainsize (gr) in µm was calculated following64:

$${{{{\rm{E}}}}}_{{{{\rm{gr}}}}}=6.62\times {{{{\rm{gr}}}}}^{-1.162}$$
(22)

Costs for grinding were then estimated using the current industry electricity price (0.15€ kWh−1). For these calculations, we assumed a near-shore grinding facility, negating the need for further transport from the harbor to the mill and back.

The CO2 emissions for the production of raw dunite were added to emissions for grinding estimated based on the energy demand (Eq. 22) and calculated for the energy mix in Germany in 202065.

For calcite, the production costs and CO2 emissions for the different grain sizes were provided by the German Limestone Association. To these, the costs and emissions for transport and spreading were added. A precise description of the assumptions made for cost and emission calculations can be found in the supplement. Finally, the combined emissions for each material were subtracted from the uptake and the final usage of material multiplied with the total costs per tonne of material.

In order to derive comparable efficiency and cost estimates, the geometrical surfaces of dunite grains were calculated for grain sizes that correspond to the BET surfaces assumed for calcite addition. For this purpose, the ratio of ABET to AGEO for calcite was calculated and assumed the same ratio for all BET surfaces. The uncertainty inherent in this assumption needs to be accepted. The AGeo was then calculated accordingly.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.