Towards circularity in shellfish aquaculture: stimulating mussel shell dissolution in marine sediments

Ocean alkalinity enhancement (OAE) is receiving considerable attention as a carbon dioxide (CO₂) removal strategy, and novel approaches to increase the total alkalinity (A_T) of the surface ocean are being explored. In bivalve aquaculture, calcification during shell growth consumes A_T, thus leading to CO₂ emissions. After consumption, shells are typically landfilled or incinerated, which can generate additional CO₂ emissions. Here, we investigate how the CO₂ footprint of shellfish aquaculture can be reduced by using bivalve shells as a resource for mineral-based OAE. The idea is to grind the calcium carbonate (CaCO₃) shells to increase the reactive surface area and distribute them into permeable, oxygen-rich sediments, where their dissolution produces A_T that could then compensate the CO₂ emitted during calcification. To evaluate this concept, we conducted microcosm incubations of sediments amended with crushed mussel shells (∼ 8 wt %), and monitored the sediment geochemistry and sediment-water exchange over 24 weeks. Control sediments exhibited low and constant CaCO₃ dissolution rates (R_diss = 0.9 ± 0.5 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and A_T fluxes (F_AT = 3.2 ± 1.1 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). In contrast, shell-amended sediments showed markedly higher R_diss and F_AT values, which exhibited a transient response modulated by oxygen and organic matter availability. Initially, shell dissolution was restricted, most likely by oxygen availability due to the intense mineralization of shell-associated organic matter. Subsequently, following gradual sediment reoxygenation, dissolution rates increased, reaching a maximum R_diss of 22.7 ± 2.6 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ after 9 weeks, corresponding to a measured F_AT of 43.0 ± 6.0 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. After that, CaCO₃ dissolution rates declined as organic matter availability decreased, approaching control rates, with a slightly elevated R_diss of 2.2 ± 1.1 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ at the end of the experiment. After 24 weeks, ∼ 6 % of the initial shell mass had dissolved, and extrapolation of the dissolution rate at the end of the experiment suggests that complete dissolution would require ∼ 38 years. Our results suggest that organic matter availability limits CaCO₃ dissolution in the permeable sediment investigated. This constraint, however, can be alleviated by targeting environments with high organic matter deposition for in-situ applications, such as sediments beneath mussel farms, thereby promoting mussel aquaculture circularity.

1 Introduction

Climate stabilization requires direct carbon dioxide removal (CDR) from the atmosphere alongside conventional mitigation strategies that prevent new emissions (Gasser et al., 2015; IPCC, 2023; Sanderson et al., 2016). As a consequence, CDR research has gained considerable traction over the last few years (Lück et al., 2025; Minx et al., 2024). Marine CDR options are receiving increased attention due to the large CO₂ buffering capacity of the ocean, which has taken up about 26 % of the global anthropogenic emissions since 1850 (Friedlingstein et al., 2025). The amount of CO₂ that can be stored in the ocean is largely controlled by total alkalinity (A_T), which is a measure of the excess of proton acceptors (bases) over proton donors (acids) in solution (Dickson, 1981; Zeebe and Wolf-Gladrow, 2001). An increase in A_T causes the carbonate system to respond by dissociating dissolved CO₂ to bicarbonate (${\mathrm{HCO}}_{\mathrm{3}}^{-}$) and carbonate (${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$), resulting in an increased uptake of atmospheric CO₂. Ocean alkalinity enhancement (OAE) represents one of the proposed marine CDR methods, and aims to establish higher oceanic CO₂ uptake by artificially increasing A_T levels in the surface ocean (Oschlies et al., 2023; Renforth and Henderson, 2017).

In recent years, a number of OAE approaches have been investigated in which minerals (silicates and carbonates) are added to coastal sediments (Fuhr et al., 2022, 2023, 2024, 2025). These minerals subsequently undergo a slow weathering process, which releases A_T. Up until now, most attention has gone towards purposely mined minerals, such as the silicate mineral olivine (e.g., Flipkens et al., 2023; Geerts et al., 2025; Hangx and Spiers, 2009; Meysman and Montserrat, 2017; Montserrat et al., 2017). Yet, mining operations require energy, generate environmental impacts, and are challenging to rapidly scale, as new mining sites require long permitting trajectories. Therefore, an alternative option is to consider minerals for OAE that are sourced from side and waste streams from existing industrial activities (Bullock et al., 2021). One such activity is bivalve aquaculture. During shell formation, bivalves remove dissolved inorganic carbon (DIC) from the water to form calcium carbonate (CaCO₃). Yet, like any form of carbonate formation, this process also removes A_T, as the precipitation of 1 mol of CaCO₃ consumes 2 mol of A_T and 1 mol of DIC:

$$\begin{array}{}\text{(1)}& {\mathrm{Ca}}^{\mathrm{2}+}+\mathrm{2}{\mathrm{HCO}}_{\mathrm{3}}^{-}\leftrightarrow {\mathrm{CaCO}}_{\mathrm{3}}+{\mathrm{CO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\end{array}$$

Figure 1 (a) Thermodynamic equilibrium plot of alkalinity (A_T) as a function of dissolved inorganic carbon (DIC) superimposed on contours of partial pressure of CO₂ (pCO₂) at 15 °C and a salinity of 35. Black dot: initial seawater concentration. Calcification (vector 1) reduces A_T and DIC in a 2:1 ratio, which increases pCO₂ and leads to CO₂ degassing to the atmosphere. Carbonate dissolution (vector 2) increases A_T and DIC in a 2:1 ratio, which decreases pCO₂ and leads to CO₂ uptake from the atmosphere. (b) Conceptual scheme illustrating the life cycle of bivalve shells from aquaculture (using mussels as example). Bivalves remove A_T and dissolved inorganic carbon (DIC) from seawater to form calcium carbonate shells (CaCO₃), a process that releases CO₂. After harvest and consumption, shell waste is typically incinerated (with associated CO₂ emission and landfill of CaO-containing ashes) or directly landfilled. The dashed arrows depict the circular route investigated here. Crushed shells are spread onto shallow sediments, where CaCO₃ dissolution generates A_T and leads to atmospheric CO₂ uptake. The overall process offsets the CO₂ released during shell growth.

Overall, carbonate formation results in an increase in the partial pressure of CO₂ (pCO₂) in seawater (Fig. 1a; Frankignoulle et al., 1994, 1995). In shallow coastal environments, where bivalve aquaculture typically takes place, rapid equilibration of the elevated pCO₂ with the atmosphere will result in outgassing of most of the CO₂ produced during CaCO₃ formation. As such, bivalve production induces a release of CO₂ into the atmosphere (Pernet et al., 2025).

With a current global production of over 17 million t of bivalves per year (FAO, 2024, 2025), bivalve aquaculture generates considerable CO₂ emissions associated with calcification (∼ 4 million  t CO₂ yr⁻¹). Life cycle assessments (LCAs) of bivalve aquaculture typically do not account for CO₂ emissions from calcification (see Ray et al., 2018 and references therein), thus resulting in an underestimation of its climate impact. When taken into account, calcification may account for up to 90 % of the total emissions (Ray et al., 2018; Thrane, 2004). Moreover, shell waste is typically disposed of as general waste, either through landfilling or incineration. While landfilling retains the CaCO₃ mineral phase, incineration can result in additional CO₂ emissions because CaCO₃ thermally decomposes into calcium oxide (CaO, quicklime) and CO₂ at high temperatures (Fig. 1b):

$$\begin{array}{}\text{(2)}& {\mathrm{CaCO}}_{\mathrm{3}}\to \mathrm{CaO}+{\mathrm{CO}}_{\mathrm{2}}\end{array}$$

Following incineration, the resulting ash may be disposed of in landfill sites. Under these conditions, some CaO may subsequently react with CO₂ to reform CaCO₃; however, the extent of this carbonation remains uncertain and depends strongly on landfill conditions.

In this study, we explore an alternative destination for the CaCO₃-containing shells after consumption. The core idea is to offset CO₂ emissions associated with calcification (Eq. 1) by reintroducing mussel shells into the marine environment and stimulating shell dissolution (Fig. 1b). The key question is whether the mussel shells will effectively dissolve, and if so, under which conditions and over which timeframes. The precipitation and dissolution of CaCO₃ are controlled by the CaCO₃ saturation state of seawater (Zeebe and Wolf-Gladrow, 2001) and the reactive surface area of the CaCO₃ minerals. Since surface waters are typically oversaturated with respect to CaCO₃, dissolution generally does not occur in the water column (Kleypas, 2011; Morse and Mackenzie, 1990). However, CaCO₃ dissolution is possible in the pore water of marine sediments under specific conditions, in a process referred to as metabolic carbonate dissolution, in which CO₂ generated by aerobic mineralization of organic matter builds up in the pore water until undersaturation is achieved and CaCO₃ starts to dissolve (Aller, 1982; Archer et al., 1989; Emerson and Bender, 1981; Lunstrum and Berelson, 2022; Morse et al., 1985; Subhas et al., 2022):

$$\begin{array}{}\text{(3)}& {\mathrm{CaCO}}_{\mathrm{3}}+{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{Ca}}^{\mathrm{2}+}+\mathrm{2}{\mathrm{HCO}}_{\mathrm{3}}^{-}\end{array}$$

The above reaction equation also succinctly summarizes the conditions that promote metabolic CaCO₃ dissolution: (1) permeable sediments that enable deep oxygen penetration; (2) sediments sufficiently rich in organic matter, which drives metabolic dissolution; and (3) low ambient CaCO₃ concentration, thereby limiting buffering from background CaCO₃ (Goossens et al., 2026). Dissolution efficiency could further be improved by grinding shells to a finer size to increase their reactive surface area and enhance their intrinsic dissolution rate. Likewise, the sediments should be located within shallow waters with a well-mixed water column as to allow rapid CO₂ equilibration with the atmosphere.

Here, we assessed the potential for enhanced sedimentary dissolution of CaCO₃ shells derived from the blue mussel (Mytilus edulis), a species that dominates shellfish aquaculture in Belgium and the Netherlands (FAO, 2025). Previous studies on biogenic shell dissolution (see Gazeau et al., 2013 for an overview) have considered a range of taxa, including oysters (Waldbusser et al., 2011; Welladsen et al., 2010), clams (Green et al., 2004, 2009; McClintock et al., 2009), limpets (McClintock et al., 2009), and snails (Nienhuis et al., 2010), but only a limited number focused on sedimentary environments (Green et al., 2004, 2009). Furthermore, relatively little attention has been given to mussel shells; the few studies that address their dissolution primarily focus on ocean acidification experiments using chemically altered (sea)water conditions (Carlson et al., 2025; Cubillas et al., 2005; Ericson and Ragg, 2022; Gazeau et al., 2007; Melzner et al., 2011). Here, we address this gap by performing laboratory flux incubations with permeable, CaCO₃-poor coastal sediments amended with crushed blue mussel shells.

Figure 2 (a) Map of the Belgian part of the North Sea (BPNS), showing the sediment collection site (Bruggen en Wegen) and the mussel farming site (Westdiep Sea Farm). An overview map of the North Sea is provided with the location of the BPNS in red. (b) Experimental setup, consisting of 3 control chambers and 3 treatment chambers amended with crushed mussel shells. Chamber and sediment dimensions are displayed in cm. (c) Chamber top lids contained 3 sampling ports: one was used for an oxygen probe, and the other two for subsampling of the overlying water. Arrows indicate water flow during subsampling: subsamples were taken from one port, while replacement water was added through the other port. (d) Stirring disk rotation induced pore-water exchange through a radial pressure gradient. Arrows indicate pore-water flow through the sediment column.

2 Materials and Methods

2.1 Sediment collection and incubation setup

Sediment was collected from a shallow site in the southern North Sea along the Belgian coast (station “Bruggen en Wegen Oostende”, 51.28300° N 2.92133° E; Fig. 2a) on 20 March 2023. The location is a dredge deposit site characterized by permeable sediment. Sediment was collected using a Van Veen grab sampler, transported to the laboratory, homogenized and stored in the dark at stable room temperature (18–22 °C) with 20 cm of overlying seawater until further use. While the sediment was not sieved prior to incubation, visual inspection during homogenization and sediment sampling did not reveal any presence of macrofauna. Air pumps were installed to keep the overlying seawater oxygenated. A subsample of sediment was taken for analysis of grain size distribution and porosity. Natural seawater was used as overlying water in the sediment incubations. This seawater was collected from the Eastern Scheldt (The Netherlands) and filtered through a sand filter (0.5 mm) to remove macrofauna. Blue mussels (Mytilus edulis) were sourced from the supermarket and originated from a hanging culture at the Westdiep Sea Farm, a mussel farming site ∼ 5 km off the Belgian coast at 15 m water depth in the southern North Sea (51.16975° N, 2.63104° E; Fig. 2a). Mussel tissue was manually removed, shells were rinsed with fresh water and dried overnight in a drying oven at 60 °C. Dried shells were ground with an ultra-centrifugal mill (ZM 200, Retsch) and the fraction that passed a 0.25 mm sieve was retained. The potential influence of fine shell particles on sediment permeability therefore cannot be completely ruled out. Ground shell samples were collected for analysis of organic and inorganic carbon concentration.

The sediment microcosms were adapted for the closed incubation of permeable sediments following a well-known flux chamber design (Huettel and Gust, 1992). The flux chambers consisted of polymethyl methacrylate (PMMA) core barrels (19 cm inner diameter, 30 cm height) fitted with a top and bottom lid (Fig. 2b and c). Top lids were equipped with a stirring disk (14 cm diameter), which mixes the overlying seawater and also creates a radial pressure gradient that drives advective pore-water exchange in the permeable sediment (Fig. 2d; Huettel and Gust, 1992). Three sampling ports were installed in the top lids (Fig. 2c): one was used for a fiber-optical oxygen probe (FireSting OXROB10, PyroScience), and the other two for subsampling of the overlying seawater during the incubation.

The experimental setup consisted of 6 flux chambers: 3 replicate control chambers were filled solely with natural sediment (∼ 13 cm), while 3 replicate treatment chambers were filled with a base layer of natural sediment (∼ 11 cm) and a top layer of natural sediment (∼ 2 cm) mixed with 100 g of ground shells, to achieve a concentration of ∼ 8 wt % of added shell material in the mixed top layer. Solid-phase samples for the analysis of total, organic, and inorganic carbon (and its stable isotope composition, i.e. δ¹³C values) were collected from the surface sediment in the control and treatment chambers. Oxygenated filtered seawater (∼ 4.2 L) was added to the chambers at the start of the incubation. Top lids were installed and chambers were placed in a closed container filled with water to ensure a dark and stable environment.

2.2 Solute flux incubation

Consecutive flux sessions were performed to quantify the uptake or release of solutes across the sediment-water interface by monitoring their accumulation in the overlying water. A total of 11 flux incubation sessions were conducted over the course of 24 weeks (each flux session lasted ∼ 10 d). Before each flux session, the overlying water in the flux chamber was replaced in order to start a new monitoring period of solute accumulation. Each flux session consisted of two parts: a short “closed” part (∼ 3 d) and a longer “open” part (∼ 7 d), as further detailed below. Throughout the session, the rotating disks were set at 80 rounds per minute (rpm) to stimulate advective pore-water exchange with the overlying water while not disturbing the sediment (i.e., no resuspension). The total time of a single incubation was adjusted to prevent excessive solute concentration build-up (i.e., to avoid saturation effects and secondary mineral precipitation). Before each session, oxygen optodes were calibrated using a two-point calibration with seawater at 0 % (saturated with Na₂SO₃) and 100 % (bubbled with air pumps) O₂ saturation.

In the first “closed” part of the incubation session, atmospheric exchange was prevented by closing the gas-tight top lids. The aim was to obtain O₂ and DIC fluxes, which necessitate isolation from the atmosphere. Closed incubations lasted until the O₂ saturation had dropped from ∼ 100 % at the start to ∼ 60 % to prevent excessive disturbance of the balance between aerobic and anaerobic processes in the sediment (which may influence CaCO₃ dissolution). Water samples (25 mL) were taken from one sampling port, while simultaneously replacing the sampled volume with stock seawater through the other port (Fig. 2c). Sampling was done at variable time intervals (5 time points), with samples collected after every ∼ 8 % decrease in O₂ saturation. Oxygen concentrations were continuously monitored at one-second intervals with the optode. Samples for dissolved inorganic carbon (DIC, 12 mL) and the stable isotope composition of DIC (δ¹³C-DIC, 6 mL) were filtered (Acrodisc Syringe Filter, 0.8/0.2 µm) and collected in glass exetainers, poisoned with 10 µL HgCl₂, and stored upside down in the dark at 4 °C.

After the closed incubation, chambers were opened and allowed to reoxygenate, and salinity was measured in the overlying water. In the subsequent “open” incubation, atmospheric exchange was enabled by opening the top lids and placing small bars between the chambers and the lids to keep the overlying water oxygenated. Water samples (60 mL) were taken directly from the overlying water at fixed time intervals (∼ 24 h in between) with 5–6 time points in total. At every time point, salinity was measured in the overlying water (as to track potential evaporation). Samples for A_T (50 mL) and nutrients (9 mL; ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$) were filtered through a dual-layer filter (Acrodisc Syringe Filter, 0.8/0.2 µm), collected in plastic vials, and stored unpoisoned in the dark at 4 and −20 °C, respectively. Both A_T and DIC samples were analyzed within the same week of collection or, at most, during the following week to ensure accurate monitoring of fluxes and minimize potential changes in sample composition during storage.

2.3 Pore-water sampling

At the end of the final flux session, 2 small core liners were inserted into the sediment and recovered for pore-water analysis. Cores were sectioned in an anaerobic glove box at 0.5 cm intervals for the 0–2 cm depth range, 1 cm intervals for the 2–6 cm depth range, and 2 cm intervals for the remaining depth. Pore water was extracted from sediment slices using rhizon samplers. Pore-water samples from one core liner were collected for Inductively Coupled Plasma Optical Emission Spectroscopy (ICP-OES) analysis (Ca²⁺; 0.5 mL), stored in plastic tubes, diluted 10 times with 0.7 M HNO₃, and stored in the dark at 4 °C. Pore-water samples from the second core liner were collected for DIC and δ¹³C-DIC (only samples ≥ 1 mL were retained), poisoned with 10 µL HgCl₂, and stored upside down in exetainers in the dark at 4 °C.

2.4 Analytical methods

Grain size distribution was determined using a laser diffraction particle size analyzer (Malvern 2000, Malvern Panalytical) with a precision of < 3 % relative standard deviation (RSD) for D₅₀, and < 5 % for D₁₀ and D₉₀. Porosity was derived from weight loss after freeze-drying, accounting for salt content in the pore water, which was assumed to be equal to the salt content in the overlying water.

Solid-phase samples were analyzed for total carbon and organic carbon content as well as δ¹³C values through combustion using an elemental analyzer (EA 1110, CE Instruments), connected via a Conflo IV Universal Continuous Flow Interface (Thermo Fisher Scientific) to an Isotope Ratio Mass Spectrometer (IRMS; Thermo Delta V Advantage, Thermo Fisher Scientific). Calibration was done using certified reference material IAEA-600 (δ¹³C = −27.77 ‰; carbon content = 49.49 %) and in-house standards Leucine (δ¹³C = −13.73 ‰; carbon content = 54.72 %) and Tuna (δ¹³C = −17.96 ‰; carbon content = 45.24 %), which were calibrated against IAEA-600. Values of δ¹³C are expressed relative to the international standard Vienna PeeDee Belemnite (VPDB), with a precision for carbon content of < 2 % (RSD) and < 0.08 ‰ (1 SD) for δ¹³C. To determine organic carbon content and δ¹³C, subsamples were weighed into Ag cups and acidified with HCl (10 %, 40 µL or until no further reaction was observed) to remove inorganic carbon. Solid-phase samples were analyzed for inorganic carbon content and associated δ¹³C composition using a Gasbench II coupled to an IRMS (Thermo Delta Plus XP, Thermo Fisher Scientific). Calibration was done using certified reference materials LSVEC (δ¹³C = −46.60 ‰) and NBS19 (δ¹³C = 1.95 ‰) and in-house standards Merck (δ¹³C = −9.65 ‰) and Fluka (δ¹³C = +2.36 ‰), which were calibrated against LSVEC and NBS19. Values were expressed relative to VPDB, with a precision of < 0.1 ‰ (1 SD).

Nutrients were analyzed using a continuous flow analyzer (SAN++, Skalar Analytical). Calibration was done using stock solutions made with NH₄Cl and NaNO₃, verified with diluted standard quality control solutions for ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (Merck Millipore, reference 1.19812.0500) and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Merck Millipore, reference 1.19811.0500), respectively. The precision (RSD) was 5.9 % (low range) and 0.8 % (high range) for ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and 1.5 % (low range) and 1.4 % (high range) for ${\mathrm{NO}}_{\mathrm{3}}^{-}$. Dissolved inorganic carbon was analyzed using a DIC analyzer (AS-C6L, Apollo SciTech) coupled to a trace gas analyzer (LI-7810, LICOR). Calibration was done using Dickson certified reference material (batch 209; DIC = 2060.05 µmol kg⁻¹) with a precision (RSD) of 0.16 %. Alkalinity was analyzed using a titrator (Titrando 888, Metrohm) with a 1 mL buret. Calibration was done using Dickson certified reference material (batch 209; A_T = 2210.40 µmol kg⁻¹) with a precision (RSD) of 0.08 %. Dissolved pore-water Ca²⁺ was analyzed using an ICP-OES elemental analyzer (AVIO500; Perkin Elmer), using indium as internal standard and the precision was generally < 2.5 % (RSD).

2.5 Solute flux calculations

Sediment-water flux calculations were performed using the R-based software FLUXER, which also provides robust data quality checks (Hylén and van de Velde, 2025). Solute fluxes (F, $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, normalized to sediment surface area) were determined from the linear regression of overlying water concentrations over time.

$$\begin{array}{}\text{(4)}& F={\displaystyle \frac{\mathrm{1}}{A}}\cdot V_{\text{ow}}\cdot {\displaystyle \frac{\mathrm{d}{C}_{\text{ow}}}{\mathrm{d}t}}=H_{\text{ow}}\cdot {\displaystyle \frac{\mathrm{d}{C}_{\text{ow}}}{\mathrm{d}t}}\end{array}$$

V_ow is the volume of the overlying water, A is the surface area of the sediment, and C_ow is the concentration of the solute in the overlying water. The term $V_{\text{ow}}/A$ can be replaced by the height H_ow of the overlying water, which was measured once during each session (as the mean of 4 different locations along the circumference of the chamber). Solute concentrations were corrected for dilution resulting from sample volume replacement in both types of incubations before flux calculations. Daily measurement of salinity during open incubations indicated no effect of evaporation.

FLUXER assumes simple linear regression by default, but can switch to a quadratic regression model when non-linear trends are detected, for example due to saturation effects from dissolved compound buildup (Devol et al., 1997; Forja and Gómez-Parra, 1998). To evaluate whether a linear of quadratic regression model is appropriate, FLUXER uses the corrected Akaike Information Criterion (cAIC; Hurvich and Tsai, 1989). The regression model with the lowest cAIC was selected, and a difference of > 2 was considered to be meaningful (Burnham and Anderson, 2004).

For each model, FLUXER displayed five diagnostic graphs to ensure good model fit and assess whether assumptions for linear regression were fulfilled (Montgomery et al., 2012): (1) residuals versus fitted values, which provided information about data linearity; (2) scale-location plot, which provided information about homoscedasticity; (3) residuals versus time, which provided information about the randomness of errors; (4) normal Q-Q (quantile-quantile) plot, which provided information about the error distribution; and (5) influence plot, which helped to identify points that had a high influence on the model (i.e., outliers) using studentized deleted residuals (SDR), hat values and Cook's distance values as indicators. We refer the reader to the theoretical background section of FLUXER for a detailed discussion of how to interpret the diagnostic graphs (Hylén and van de Velde, 2025).

2.6 CaCO₃ dissolution rates

Carbonate dissolution rates were calculated using the Keeling plot method (Keeling, 1958, 1961; Pataki et al., 2003), which is based on an isotope mass balance. In closed incubations, the DIC concentration (µmol L⁻¹) of the overlying water (DIC_ow) is the sum of the background DIC concentration at the start of the incubation (DIC_b) and the source DIC concentration (DIC_s), which is added to the overlying water by sedimentary processes:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathrm{DIC}}_{\mathrm{ow}}={\mathrm{DIC}}_{\mathrm{b}}+{\mathrm{DIC}}_{\mathrm{s}}\text{(6)}& {\displaystyle }\begin{array}{rl}\left({\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}{\mathrm{DIC}}_{\mathrm{ow}}\right)\cdot {\mathrm{DIC}}_{\mathrm{ow}}=& \left({\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}{\mathrm{DIC}}_{\mathrm{b}}\right)\\ & \cdot {\mathrm{DIC}}_{\mathrm{b}}\left({\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}{\mathrm{DIC}}_{\mathrm{s}}\right)\cdot {\mathrm{DIC}}_{\mathrm{s}}\end{array}\end{array}$$

Combining Eqs. (5) and (6) gives a first-order equation in the measured values of δ¹³C-DIC (‰) of overlying water (δ¹³C-DIC_ow) over $\mathrm{1}/{\mathrm{DIC}}_{\mathrm{ow}}$.

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}{\mathrm{DIC}}_{\mathrm{ow}}=& {\mathrm{DIC}}_{\mathrm{b}}\cdot \left({\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}{\mathrm{DIC}}_{\mathrm{b}}-{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}{\mathrm{DIC}}_{\mathrm{s}}\right)\\ & \cdot \left({\displaystyle \frac{\mathrm{1}}{{\mathrm{DIC}}_{\mathrm{ow}}}}\right)+{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}{\mathrm{DIC}}_{\mathrm{s}}\end{array}\end{array}$$

The isotope ratios of source and background can be obtained by a linear regression. The intercept yields δ¹³C-DIC_s. Since there are two sedimentary processes that generate DIC, i.e., organic matter mineralization and CaCO₃ dissolution, that each have a distinct isotopic signature, the relative contribution of each process to the total DIC can be derived.

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{f}_{\text{OM}}+f_{\text{carb}}=\mathrm{1}\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}{\mathrm{DIC}}_{\mathrm{s}}=f_{\text{carb}}\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{PIC}}+f_{\text{OM}}\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{POC}}\end{array}$$

The fraction of DIC generated by carbonate dissolution hence becomes:

$$\begin{array}{}\text{(10)}& f_{\text{carb}}={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{-}{\mathrm{DIC}}_{\mathrm{s}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{POC}}}{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{PIC}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{POC}}}}\end{array}$$

In these equations, f_carb and f_OM are the fractions of total DIC derived from CaCO₃ dissolution and organic matter mineralization, respectively, while δ¹³C_PIC and δ¹³C_POC are the isotopic compositions of CaCO₃ and organic matter. At saturation (Ω = 1), the theoretical maximum f_carb is 0.5 (50 %), corresponding to fully efficient metabolic carbonate dissolution, where half of the total DIC pool originates from organic matter mineralization and the other half from CaCO₃ dissolution (Eq. 3). Finally, the CaCO₃ dissolution rate (R_diss, $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, normalized to sediment surface area) and metabolic CaCO₃ dissolution efficiency (MCDE, %) can be calculated from f_carb and the measured DIC flux (F_DIC) via:

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{R}_{\text{diss}}=f_{\text{carb}}\cdot F_{\mathrm{DIC}}\text{(12)}& {\displaystyle }& {\displaystyle }\text{ MCDE}=\mathrm{2}\cdot f_{\text{carb}}\cdot \mathrm{100}\end{array}$$

Note that our experimental setup precludes any external methane sources.

Table 1 Particulate organic carbon (POC), particulate inorganic carbon (PIC), and calcium carbonate (CaCO₃) concentrations (assuming all PIC is CaCO₃) of the sediments in the control and treatment chambers, as well as the added mussel shells. The δ¹³C isotopic compositions are also provided. Values for treatment chambers are based on the top layer (2 cm) composition (sediment mixed with crushed mussel shells). Uncertainties of the chamber values are standard deviations based on triplicate chambers. Uncertainties of the mussel shell values are standard deviations based on 4 samples.

3 Results

3.1 Sediment properties

The natural sediment had a median grain size of 262 µm and porosity of 0.45 ± 0.03 (n=3) and can be classified as medium sand according to the Wentworth scale. Consistent with a permeable sediment, the organic carbon content was low (0.08 %), while the background level of inorganic carbon was also moderate (0.95 %). Particulate organic carbon (POC) and particulate inorganic carbon (PIC) concentrations of the different substrates used in the experiments (including the added mussel shells) and their corresponding δ¹³C values are summarized in Table 1.

Figure 3 (a) Oxygen fluxes ($F_{{\mathrm{O}}_{\mathrm{2}}}$). (b) Dissolved inorganic carbon fluxes (F_DIC). (c) Alkalinity fluxes (F_AT). Positive fluxes indicate fluxes into the overlying water (i.e., net production by the sediment), negative fluxes indicate fluxes into the sediment (i.e., net consumption by the sediment). Means are plotted with error bars representing standard deviations based on three replicate chambers. No O₂ and DIC data were available for week 7, as only open incubations were performed during that week. (d) 1:1 plot, comparing F_DIC and $F_{{\mathrm{O}}_{\mathrm{2}}}$. (e) 1:1 plot, comparing F_AT and F_DIC.

3.2 Sediment-water exchange fluxes

Most incubations showed linear changes in overlying water solute concentrations over time (Figs. S1–S5 in the Supplement). An exception was observed during week 1, where A_T accumulation in the treatment chambers was nonlinear, with the accumulation rate decreasing halfway through the incubation (Fig. S3). Similarly, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations in the treatment chambers displayed nonlinear trends during week 1 (Figs. S4 and S5): concentrations of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ remained constant in the beginning of the incubation, after which they increased, while ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations increased at the beginning and decreased at the end.

Fluxes of O₂, DIC and A_T were calculated by regression of overlying water concentrations with time, and showed good replication between replicates at each time point (Fig. 3a–c). Fluxes showed a clear difference between control ($F_{{\mathrm{O}}_{\mathrm{2}}}$ = −5.4 ± 1.2 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, F_DIC = 7.4 ± 1.0 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, F_AT = 3.2 ± 1.1 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; mean flux ± standard deviation over all replicate chambers and flux sessions) and treatment chambers ($F_{{\mathrm{O}}_{\mathrm{2}}}$ = −16.8 ± 4.3 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, F_DIC = 32.0 ± 9.4 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, F_AT = 25.1 ± 13.5 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; Table S1 in the Supplement). While fluxes in the control chambers showed little variation over time, the treatment chambers displayed a clear temporal pattern. Both O₂ (Fig. 3a) and DIC (Fig. 3b) fluxes were initially markedly higher than those of the control chambers, and gradually decreased over the course of the experiment, approaching fluxes of the controls towards the end. Fluxes of A_T first increased with time until a maximum value of 48.0 ± 6.4 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ was reached at week 7, after which they decreased and approached the fluxes of the controls towards the end (Fig. 3c). Fluxes of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in the control chambers were negligible throughout all flux sessions (Table S1). In the treatment chambers, elevated ${\mathrm{NH}}_{\mathrm{4}}^{+}$ effluxes were only observed during the first part of week 1 ($F_{{\mathrm{NH}}_{\mathrm{4}}^{+}}$ = 5.1 ± 1.5 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), after which they dropped below the detection limit (Table S1). A positive ${\mathrm{NO}}_{\mathrm{3}}^{-}$ efflux was observed in the second part of week 1 ($F_{{\mathrm{NO}}_{\mathrm{3}}^{+}}$ = 5.8 ± 1.0 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; Table S1), and in weeks 3 and 5, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ fluxes were still observed, although they were substantially lower (∼ 1 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Finally, flux scatter plots revealed a clear additional DIC source in the treatment chambers beyond mineralization-derived DIC (Fig. 3d), as well as a proportionally greater increase in A_T relative to DIC (Fig. 3e).

Figure 4 (a) Mean fraction of total DIC derived from CaCO₃ dissolution (f_carb) as determined using the Keeling plot method. Values are expressed as percentages. (b) Corresponding mean CaCO₃ dissolution rates (R_diss). Values of R_diss were only calculated for positive values of f_carb, as discussed in the main text. Error bars are standard deviations based on three replicate chambers.

3.3 CaCO₃ dissolution rates

Figure 4a and b and Table S2 show the derived f_carb and corresponding R_diss values across the experiment, while the corresponding MCDE values are provided in Table S2 (see Fig. S6 for all Keeling plots and Table S2 for all model intercepts and p values). A temporal pattern was observed in the treatment chambers over the course of the experiment, with f_carb and R_diss increasing over time, reaching a maximum of f_carb = 59.5 ± 6.5 % (R_diss = 22.7 ± 2.6 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) in week 9, followed by a gradual decrease towards a final f_carb of 17.3 ± 12.7 % (R_diss = 2.2 ± 1.1 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). In the control chambers, f_carb remained relatively stable over time, although some variability was observed both among replicate chambers and across incubation weeks. As the Keeling plot method combines measurements of both DIC and δ¹³C-DIC, each with associated uncertainties, such variability is expected, particularly in control chambers, where concentration and isotopic composition changes over time are relatively small. The occurrence of negative f_carb values in some control chambers indicated that the regression intercept fell outside the isotopic range defined by the CaCO₃ and organic matter end-members. This has no geochemical meaning, other than suggesting that little to no CaCO₃ dissolution is occurring, and underscores the uncertainty that the Keeling plot approach inherently incorporates when variations in DIC and/or δ¹³C-DIC are relatively modest. For this reason, R_diss was only calculated for chambers with positive f_carb values (Fig. 4b).

3.4 Pore-water geochemistry

We tracked the color changes in the sediment over time by taking photographs throughout the experimental period (Fig. S7). These images reveal a markedly different color evolution in control versus treatment, thus suggesting a different development of the geochemistry. In the controls, the sediment remained light-grey colored, and the only color change observed was the development of a thin (∼ 5 mm) orange layer just below the sediment-water interface (Fig. S7). This feature suggests remobilization of ferrous iron in deeper layers, which then precipitates out as iron (hydr)oxide upon upward diffusion and contact with O₂. In contrast, in the treatments, there was a far more notable color change, as the entire layer with mussel shell fragments became black over the course of the first 3 weeks, after which it gradually became grey again as the experiment progressed (Fig. S7).

Figure 5 Pore-water profiles. (a) Calcium (Ca²⁺). (b) Dissolved inorganic carbon (DIC). (c) Isotopic DIC composition (δ¹³C-DIC). (d) 1:1 plot, comparing Ca²⁺ and DIC concentrations. Pore-water samples were only taken once at the end of the experiment, explaining the buildup of solutes in deeper anoxic layers, as further discussed in the main text.

Figure 5 displays the depth profiles of pore-water concentrations (Ca²⁺, DIC) and isotopic composition (δ¹³C-DIC) recorded at the end of the experiment. Higher pore-water Ca²⁺ concentrations were observed in the treatment chambers as compared to the control chambers (Fig. 5a). In the control chambers, pore-water Ca²⁺ concentrations remained constant with depth, and were similar to the concentration of ∼ 9.5 mM of the overlying seawater. In contrast, in the treatment chambers, dissolved Ca²⁺ markedly increased with depth, reaching 12–13 mM at 12 cm depth. Pore-water DIC profiles showed an increasing trend with sediment depth for both control and treatment chambers, but with a different shape and higher DIC concentrations observed in the treatment chambers (Fig. 5b). Profiles of δ¹³C-DIC showed a decreasing trend in δ¹³C values with sediment depth in the control chambers. In the treatment chambers, however, δ¹³C values decreased in the top layer but then remained relatively stable below a depth of ∼ 1 cm (Fig. 5c). Pore-water 1:1 plots of Ca²⁺ and DIC concentrations revealed a clear Ca²⁺ source in the treatment chambers, alongside higher DIC production (Fig. 5d).

4 Discussion

4.1 Mussel shell dissolution provides a significant additional A_T source

Metabolic CaCO₃ dissolution is thought to occur more efficiently in oxygenated sediments, where the aerobic mineralization of organic matter produces the CO₂ that drives the pore water towards undersaturation (Archer et al., 1989; Emerson and Bender, 1981; Lunstrum and Berelson, 2022; Morse et al., 1985). To meet this condition, permeable sediment (medium sand; porosity 0.45) was used here in the experiment. While the natural sediment employed does contain a sizeable baseline amount of PIC (0.95 ± 0.20 %; Table 1), low A_T fluxes were observed in the control chambers (Fig. 3c). Likewise, the Keeling plot analysis indicated low rates of metabolic carbonate dissolution (Fig. 4b). The mean CaCO₃ dissolution rate of 0.9 ± 0.5 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ as observed in the controls, corresponds to an A_T flux 1.8 ± 1.0 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, thus providing ∼ 56 % of the total measured A_T flux of 3.2 ± 1.1 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. The low organic carbon content of the natural sediment (0.08 ± 0.02 wt %; Table 1) resulted in low mineralization rates, as indicated by the low O₂ consumption rates (Fig. 3a), which are a proxy for total sedimentary organic matter mineralization (Glud, 2008; Jørgensen et al., 2022). The low CO₂ production from mineralization thus constrained metabolic CaCO₃ dissolution in the natural sediment.

Figure 6 Shell-amended chamber flux comparison of total A_T, which was measured directly in the overlying water, and CaCO₃ A_T, which was derived using the Keeling plot analysis. (a) Mean weekly fluxes with standard deviations based on three replicate chambers. (b) 1:1 plot, comparing the individual weekly total and CaCO₃-derived A_T fluxes per chamber.

In the treatment chambers, more favorable conditions for metabolic dissolution were established, as the added mussel shells still had a considerable amount of associated organic matter (2.20 wt %; Table 1), which stimulated organic matter mineralization, as reflected in the O₂, DIC and A_T fluxes (Fig. 3a–c). First, elevated O₂ fluxes in the treatment chambers suggest high mineralization rates driven by the organic matter associated with the shells, and thus an increased metabolic dissolution potential. Second, elevated DIC fluxes in the treatment chambers confirmed increased metabolic dissolution driven by the more intense mineralization (Fig. 3d). The Keeling plot analysis demonstrates there is a clear additional DIC source coming from CaCO₃ dissolution in the treatment chambers (Fig. 4b), which is further supported by the higher Ca²⁺ and DIC pore-water concentrations observed in the treatment chambers (Fig. 5a, b, and d). It should be noted that pore-water samples were only taken once at the end of the experiment (24 weeks), which allowed for the downward buildup of these solutes over time. While enhanced CaCO₃ dissolution only occurred in the shell-amended top layer in the treatment chambers, downward diffusion likely resulted in the elevated concentrations in the deeper anoxic layers. Likewise, the overlying water was replaced after every incubation, which enabled upward diffusive and advective transport of Ca²⁺ into the overlying water, thus resulting in a lower buildup in the top layer. Third, high A_T fluxes in the treatment chambers (Fig. 3c and e) also provide support for substantial CaCO₃ dissolution, as 2 mol of A_T are produced per mol of CaCO₃ dissolved (Eq. 1). In addition, comparison of the A_T fluxes attributed to CaCO₃ dissolution (as derived from the Keeling plot analysis) with the total A_T fluxes (as measured in the overlying water) revealed a strong linear correlation for nearly all flux sessions (Fig. 6a and b), demonstrating that CaCO₃ dissolution was the dominant A_T source in the treatment chambers. A small number of flux sessions (weeks 3, 11, and 13) deviated from this relationship, indicating the presence of additional A_T-generating or A_T-consuming processes. These processes are discussed in the following section.

4.2 Factors controlling metabolic CaCO₃ dissolution

The observed geochemistry in the treatment chambers carries a clear signature of metabolic CaCO₃ dissolution. Given that metabolic CaCO₃ dissolution is constrained by the availability of oxygen, organic matter, and CaCO₃ (Eq. 3), an important question is which of these three factors acts as the primary limiting factor. Metabolic carbonate dissolution efficiencies (MCDE, Eq. 12) offer a useful metric to assess this issue. Efficiencies < 100 % can be due to different processes: (1) the CaCO₃ pool is either too small or not reactive enough for dissolution to fully utilize the potential of mineralization, (2) the pore water receives A_T from anaerobic mineralization processes, which suppresses the undersaturation with respect to CaCO₃ (Soetaert et al., 2007), or (3) the acid (i.e. CO₂) produced by aerobic mineralization is flushed out of the pore water by advective transport before it can drive the pore water towards undersaturation (Goossens et al., 2026). We now discuss these different factors in more detail.

4.2.1 Limitation by CaCO₃ availability and reactivity

To assess whether CaCO₃ availability limited dissolution in the treatment chambers, the total amount of CaCO₃ dissolved over the course of the experiment was calculated by integrating the R_diss values obtained from the Keeling plot analysis over their respective incubation periods and the sediment surface area (Fig. 4b). This resulted in a total of 6.1 ± 1.2 g of CaCO₃ dissolved, corresponding to 6.1 ± 1.2 % of the initial shell mass that was added. These results suggest that CaCO₃ availability most likely did not limit dissolution.

In addition to CaCO₃ content, CaCO₃ dissolution efficiency can also be limited by reactivity. When pore waters transition from over- to undersaturated conditions, carbonate minerals dissolve more readily at high-energy surface sites such as steps and kinks, which are more reactive than smoother surfaces (Adkins et al., 2021). In this experiment, mussel shells were freshly crushed, resulting in more high-energy surface sites relative to the natural CaCO₃ present in the sediment. This suggests that CaCO₃ reactivity was unlikely to be limiting carbonate dissolution in the shell-amended chambers.

Figure 7 Conceptual figure illustrating the three distinct phases observed after mussel shell amendment. Black dots represent mean metabolic CaCO₃ dissolution efficiencies (MCDE) determined using the Keeling plot method. During phase 1, MCDE is limited by oxygen (O₂) availability. During phase 2, CaCO₃ dissolution is limited by organic matter (OM) availability while occurring at maximum efficiency. During phase 3, dissolution efficiency is limited by OM availability, and both dissolution efficiencies and rates gradually approach those observed in sediments without added shells. In a real-world shell amendment scenario, however, metabolic CaCO₃ dissolution would be continuously fueled by the supply of fresh OM from the water column, potentially driving the system back towards maximum efficiency.

4.2.2 Limitation by anaerobic mineralization and pore-water flushing

Undersaturation of pore waters with respect to CaCO₃ is promoted by oxic conditions, which drive CO₂ production through aerobic mineralization. In coarse sediments, this is stimulated by advection, which allows deeper oxygen penetration (Huettel and Gust, 1992) and increases the zone in which CaCO₃ undersaturation can develop. The degree of undersaturation is further controlled by the pore-water residence time (PRT), which is the mean time that seawater remains within the sediment before being replaced by new overlying water. If the PRT is too short, acid produced by aerobic mineralization is flushed out before it can drive the pore water towards undersaturation. If the PRT is too long, A_T from CaCO₃ dissolution will accumulate and drive the pore water back towards saturation (Goossens et al., 2026).

Below the oxic zone, anaerobic mineralization generates A_T instead of acidity, along with reduced metabolites (Soetaert et al., 2007). However, these reduced compounds are typically transported upward into the oxic zone, where their reoxidation produces acidity that lowers the CaCO₃ saturation state of the pore water. This means that metabolic carbonate dissolution can occur directly via aerobic mineralization and indirectly via anaerobic mineralization coupled to reoxidation in the oxic zone (Dale et al., 2024). Iron sulfide precipitation, however, represents an exception to this process.

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathrm{2}{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+\mathrm{2}{\mathrm{H}}^{+}\to \mathrm{2}{\mathrm{CO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{S}+\mathrm{2}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\text{(14)}& {\displaystyle }{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+\mathrm{4}\mathrm{FeOOH}+\mathrm{8}{\mathrm{H}}^{+}\to {\mathrm{CO}}_{\mathrm{2}}+\mathrm{4}{\mathrm{Fe}}^{\mathrm{2}+}+\mathrm{7}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\end{array}$$

$$\begin{array}{}\text{(15)}& \begin{array}{rl}\mathrm{15}{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+\mathrm{8}{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}& +\mathrm{4}\mathrm{FeOOH}+\mathrm{16}{\mathrm{H}}^{+}\\ & \to \mathrm{15}{\mathrm{CO}}_{\mathrm{2}}+\mathrm{4}{\mathrm{FeS}}_{\mathrm{2}}+\mathrm{25}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\end{array}\end{array}$$

During sulfate reduction, hydrogen sulfide (H₂S) is produced (Eq. 13). When H₂S reacts with reactive iron generated through iron reduction (Eq. 14), iron sulfides can precipitate (Eq. 15, shown for pyrite formation). This prevents the upward transport of H₂S and its subsequent reoxidation, resulting in a net source of A_T (Hu and Cai, 2011; Middelburg et al., 2020). Under such conditions, the efficiency of metabolic carbonate dissolution is reduced since a potential source of acidity from mineralization is removed.

Ultimately, the interplay between aerobic/anaerobic mineralization and pore-water exchange dictates whether pore waters become under- or oversaturated, thus controlling the MCDE. Throughout the experiment, three distinct phases can be distinguished (Fig. 7): phase one (weeks 1–9) was characterized by initially low MCDE values that gradually increased over time; phase two (weeks 9–15) was characterized by maximum MCDE values (∼ 100 %); and phase three (weeks 15–24) was characterized by decreasing MCDE values over time.

During the first phase of the experiment (weeks 1–9), intense mineralization took place, likely driven by the fresh organic matter associated with the mussel shells. This process was immediately visibly apparent, as during the first week, the top layer of the shell-amended sediment (where the mussel shell fragments had been added) turned black, indicative of sediment anoxia and potential iron sulfide formation (Fig. S7). The oxygen consumption in the treatment chambers was high, suggesting that the oxic zone was likely strongly diminished. Week 1 was also the only week during which we observed ${\mathrm{NH}}_{\mathrm{4}}^{+}$ accumulation in the overlying water, suggesting that nitrification was inhibited, most likely by a lack of oxygen. In all ensuing flux sessions, the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux from the sediment remained negligible (Table S1). In the latter part of week 1, oxic conditions apparently reestablished, as nitrification initiated, leading to a depletion of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (Fig. S5) and a buildup of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Fig. S4). Nitrification consumes A_T by releasing protons, which likely contributed to the observed nonlinear A_T accumulation in week 1. Anaerobic mineralization generates A_T, which elevates the CaCO₃ saturation state in the pore water and thus likely contributed to the low MCDE (∼ 45 %) observed in week 1.

In the rest of the first phase (after week 1), the mineralization remained high (as indicated by high O₂ and DIC fluxes; Fig. 3a and b), but the pore waters recovered from the initial “mineralization burst” and became gradually more oxygenated, stimulating metabolic CaCO₃ dissolution (Eq. 3). As a result, the MCDE increased with time, until reaching its maximum value in week 9.

During the second phase (weeks 9–15), carbonate dissolution occurred at maximum efficiency in the shell layer. The δ¹³C-DIC pore-water profiles taken at the end of the experiment indicated that organic matter (δ¹³C = −21.8 ‰) and CaCO₃ (δ¹³C = −0.7 ‰) contributed approximately equally to the pore-water DIC pool, as reflected by the constant δ¹³C-DIC value of ∼ −10 ‰ in the deeper sediment layers (Fig. 5c). While this end-point profile does not reveal when these processes occurred, it indicates that, overall, metabolic dissolution contributed substantially to the accumulated pore-water DIC. Given that R_diss values were highest during the second phase, this integrated signal is expected to be dominated by this period, during which the dissolution efficiency was highest (MCDE ≈ 100 %; Fig. 7). Maximum efficiencies also indicate that R_min in the shell layer was sufficiently large to overwhelm any potential effect of iron sulfide formation below the oxic layer on f_carb.

In weeks 9–11, the MCDE even exceeded 100 % (Fig. 7). Theoretically, this is not possible when only metabolic dissolution is at play, but it can occur when an additional acid source is present. The observed MCDE > 100 % are hence either related to uncertainty in their estimation, or they can also result from the transient reoxidation of a pool of iron sulfides, a process which generates acid and consumes A_T (Rimstidt and Vaughan, 2003). After week 9, the intensity of mineralization decreased, which may have caused an increased and deeper oxygenation, and hence, a reoxidation of previously accumulated iron sulfides, consistent with the observed offset between total A_T fluxes and A_T fluxes attributed to CaCO₃ dissolution (Fig. 6). The observed continuous transition of the sediment color from black to light grey is also consistent with this interpretation, as it suggests potential reoxidation of iron sulfides (Fig. S7).

During the third phase (weeks 15–24), MCDE decreased with time (Fig. 7). At this point, the most labile fraction of organic matter had been mineralized, leading to reduced mineralization rates and decreased O₂ consumption rates (Fig. 3a). At this point, CaCO₃ dissolution likely became limited by organic matter availability, and dissolution rates decreased accordingly (Fig. 4b). The decreasing R_min could not generate sufficient acidity to drive the pore water fully towards undersaturation before being flushed, thereby reducing R_diss, f_carb and MCDE. Eventually, the treatment chambers appeared to evolve towards the control rates, with R_diss remaining slightly higher than in the control chambers (Fig. 4b).

At the end of the experiment, CaCO₃ dissolution appeared limited by organic matter availability. We hence speculate that the addition of new fresh organic matter could again increase the dissolution rate and associated MCDE (Fig. 7). The latter aspect was not examined in the current experiment but provides a testable hypothesis for future studies.

4.3 Enhanced CaCO₃ dissolution in the global aquaculture context

The flux incubations performed here provide an estimate for the areal CaCO₃ dissolution rate and A_T flux induced by enhanced mussel shell dissolution. These values can be used to arrive at a first-order estimate of the potential of bivalve dissolution approaches within an OAE context. When extrapolating our results to the global aquaculture context, we need to consider both CO₂ emissions associated with mussel shell aquaculture as well as CO₂ sequestration obtained from a shell amendment scenario.

Global yearly CO₂ emissions (${\mathrm{m}}_{{\mathrm{CO}}_{\mathrm{2}}}$; t CO₂ yr⁻¹) associated with calcification in mussel aquaculture can be calculated as:

$$\begin{array}{}\text{(16)}& m_{{\mathrm{CO}}_{\mathrm{2}}}=m_{\text{mussel}}\cdot x_{\text{shell}}\cdot x_{{\mathrm{CaCO}}_{\mathrm{3}}}\cdot \mathit{\psi }\cdot {\displaystyle \frac{M_{{\mathrm{CO}}_{\mathrm{2}}}}{M_{{\mathrm{CaCO}}_{\mathrm{3}}}}}\end{array}$$

In this, m_mussel represents the total mass of harvested mussels from aquaculture (t yr⁻¹), x_shell denotes the shell fraction of the total mussel mass (∼ 29 %; Adler et al., 2024), $x_{{\mathrm{CaCO}}_{\mathrm{3}}}$ represents the CaCO₃ fraction of shells (92.65 %; Table 1), ψ is the molar ratio of CO₂ released versus CaCO₃ precipitated (Frankignoulle, 1994; Frankignoulle et al., 1994), and $M_{{\mathrm{CO}}_{\mathrm{2}}}$ and $M_{{\mathrm{CaCO}}_{\mathrm{3}}}$ are the molar masses of CO₂ (44.01 g mol⁻¹) and CaCO₃ (100.01 g mol⁻¹), respectively. Calculation of ψ was done using the CRAN:seacarb package in R software (Gattuso et al., 2024), with input parameters representative of surface seawater conditions: DIC = 2100 µmol kg⁻¹, salinity = 35, temperature = 15 °C, and atmospheric pCO₂ = 427 µatm (Lan et al., 2023). This yielded a value of ψ = 0.72. Blue mussel aquaculture has a global production of approximately 132 kt yr⁻¹ (FAO, 2025), and so, the CO₂ emissions linked to calcification are ∼ 11.2 kt CO₂ yr⁻¹. Using the results obtained here, 6.1 ± 1.2 % of the initial shell mass was dissolved after 24 weeks, and so the corresponding CO₂ compensation over the same time period would be 6.1 % of the calcification-related emissions incurred during the growth phase (representing 686 ± 135 t CO₂ globally). It should be noted, however, that at the end of the experimental period, the CaCO₃ dissolution rate and corresponding A_T flux in the treatment chambers were still higher than those in the control chambers, thus indicating that enhanced shell dissolution would still continue for a longer period. Assuming the dissolution rate stays constant in time, it would take ∼ 38 years to completely dissolve the added mussel shell fragments, if we use the dissolution rate measured at the end of the experiment. Alternatively, if we use the average dissolution rate over the experiment, the dissolution timescale reduces to ∼ 7 years. However, the extent to which crushed shells remain exposed to oxic conditions over these timescales depends on sediment accumulation rates, which typically range between 0.1 and 1 cm yr⁻¹ in shelf environments (Restreppo et al., 2020), and oxygen penetration depths. Assuming an oxygen penetration depth of 2 cm, corresponding to the thickness of the amended sediment layer in our experiment, this suggests that crushed shells may remain exposed to oxic conditions for ∼ 2–20 years. Once buried below the oxic layer, precipitation may become favored again in oversaturated pore waters, thereby reducing the potential for continued enhanced mussel shell dissolution. Bioturbation, however, can substantially increase the depth of oxygen penetration through bio-irrigation (Aller, 1980; Meysman et al., 2006), and extend the residence time of crushed shells within the oxygenated zone. Considering an average bioturbated layer depth of 10 cm (Boudreau, 1998), the residence time of crushed shells within this zone increases to ∼ 10–100 years (Kuderer and Middelburg, 2024). Therefore, in bioturbated sediments, crushed shells may remain exposed to conditions favorable for metabolic dissolution long enough to approach complete dissolution.

Furthermore, our results clearly indicate that carbonate dissolution was limited by organic matter availability from week ∼ 15 onward in our experimental setup, leading to lowered mineralization rates and decreased metabolic CaCO₃ dissolution. Such a depletion of fresh organic matter is a direct consequence of the incubation procedure. In real-world application scenarios, fresh organic matter would be continuously supplied from the water column, especially during spring blooms. This would fuel metabolic CaCO₃ dissolution, as conceptually illustrated in Fig. 7, and hence reduce the dissolution time of bivalve shell material. Consequently, dissolution under field conditions may proceed considerably faster than suggested by the estimates derived from our laboratory experiment, increasing the likelihood that complete dissolution occurs before permanent burial, even in sediments with limited bioturbation. As such, the potential of mussel-shell amendment as a strategy to compensate for aquaculture CO₂ emissions could be greater than suggested by our laboratory experiment. This interpretation is supported by a recent study, which simulated an application scenario of two types of biominerals to coastal sediments by means of a reactive transport model. These simulations showed that enhanced A_T production from shell dissolution could be sustained for over 30 years before transitioning back to steady state (Biçe et al., 2025). Sensitivity analysis also confirmed the limiting effect of lowered mineralization rates on CaCO₃ dissolution rates, supporting our findings that shell amendment would be most effective in environments with high organic matter degradation (Biçe et al., 2025).

Enhanced bivalve shell dissolution is most promising in CaCO₃-poor permeable sediments residing in shallow areas with high local primary production. Globally, ∼ 70 % of the continental shelf sediments are classified as sandy (i.e. permeable; Burdige, 2007), while ∼ 88 % is classified as CaCO₃-poor (Smith and Mackenzie, 2016), thus showing that the coastal environment has a large areal potential for OAE applications based on enhanced bivalve shell dissolution. Bivalve aquaculture produces approximately 12 million  t of shell material per year (FAO, 2025), which at the concentrations employed here (8 % addition) would require ∼ 5000 km² of new application area per year. This area only represents 0.04 % of the CaCO₃-poor, sandy sediments found on the global continental shelf, and hence, the space for application does not appear to be limiting.

Still, the selection of an application site for enhanced bivalve shell dissolution requires judicious consideration. Sandy sediments are generally low in organic matter content, and so their potential for metabolic CaCO₃ dissolution is highly dependent on a continuous local input of organic matter (de Beer et al., 2005; Goossens et al., 2026). Therefore, sandy sites with a suitably increased flux of organic matter towards the seabed seem preferential. When the organic matter input is too low, metabolic carbonate dissolution becomes limited by organic matter. Likewise, when organic matter inputs are too high, metabolic carbonate dissolution can become oxygen limited. Intensified mineralization can deplete sediment oxygen and create sulfidic conditions, which may remove the presence of bioturbating infauna (Christensen et al., 2003), and promotes A_T production by iron sulfide formation. Under such circumstances, CaCO₃ dissolution shifts from being organic-matter-limited to being oxygen-limited. Finally, hydrodynamic conditions will determine the fate of added shell material. Resuspension events are common in shallow coastal environments, but their frequency and intensity depend strongly on local wave and current forcing (Guillén et al., 2002; Miles et al., 2013). Resuspension may have contrasting effects on shell dissolution: while physical disturbance can increase oxygen exposure and abrasion of shelf particles, it also increases exposure to oversaturated water and may redistribute particles away from the application area into environments with less favorable conditions for dissolution. Therefore, local sediment transport dynamics should be considered when selecting suitable application sites.

One potential application scenario could be to redistribute shell waste in or near mussel aquaculture environments, which are known to have induced high organic matter deposition rates (Christensen et al., 2003; Hargrave et al., 2008; Lavoie et al., 2024). This would also create direct compensation for the calcifying-induced CO₂ emissions that are locally generated (i.e. within-project offsetting). Moreover, it could potentially also improve aquaculture conditions by releasing A_T and increasing the seawater pH, although this effect will strongly depend on the residence time of the water at the site. Previous studies have investigated organic matter accumulation beneath mussel farms and its effect on the local sediment geochemistry, showing that the magnitude of the impact depends on multiple factors. While the mussel density determines the total organic matter flux from the water column, the distance from the farm and the local hydrodynamics control how much of this material is deposited within the sediment (e.g., Christensen et al., 2003; Cranford et al., 2007; Hargrave et al., 2008; Lacoste et al., 2018; Lavoie et al., 2024; McKindsey et al., 2011; Wilson and Vopel, 2015). This emphasizes the need to carefully characterize sediment conditions prior to mussel shell waste application. Optimal sites could be areas with low hydrodynamic forcing that are located near but not directly beneath mussel farms, allowing elevated organic matter deposition without inducing anoxia, or areas with high hydrodynamic forcing directly beneath mussel farms, where organic matter fluxes are high but sufficiently dispersed to prevent highly localized deposition.

5 Summary and outlook

In this study, we investigated the dissolution of mussel shell waste in permeable sediments as a novel type of OAE approach, which aims for enhanced circularity and reduced emissions in bivalve aquaculture operation. While mussels (and bivalves in general) already have a low CO₂ footprint per gram of protein compared to other food sources, this CO₂ footprint could be lowered by dissolving the shells in marine conditions. Permeable sediment amended with crushed mussel shells exhibited significantly higher dissolution rates, driven by aerobic mineralization linked to organic matter associated with the shells. A transient dissolution pattern was observed, where initial oxygen limitation was followed by organic matter limitation. After 24 weeks of sediment incubation, 6.1 ± 1.2 % of the initial shell mass had dissolved, which upon extrapolation provides a dissolution timescale of ∼ 38 years. Deploying mussel shell waste in natural sediments with high and continuous organic matter deposition, however, could increase the compensation potential and shorten the dissolution timescale. Potential application areas include mussel farms, which would allow for direct in-project CO₂ compensation, with the potential of improving aquaculture conditions. Careful characterization of the sediment geochemistry prior to mussel shell waste application is crucial to ensure optimal conditions for CaCO₃ dissolution.

While enhanced mussel shell dissolution seems feasible, it should be noted that the results presented here are derived from a controlled laboratory experiment. In contrast, natural sediments experience seasonal variation in temperature and organic matter deposition, are subject to variable hydrodynamic forcing (including storms), and experience infaunal activity such as bioturbation and bio-irrigation, all of which can influence oxygen availability or benthic mineralization rates. Bridging the gap between laboratory and field conditions will require experiments conducted under more realistic settings, such as large-scale mesocosms and eventually field trials. To date, OAE studies are heavily inclined towards modeling approaches. While models are crucial for exploring large-scale scenarios, their predictive capacity ultimately depends on experimental data, which is currently still largely lacking. In this study, we tested a single mussel shell loading in a single sediment type under controlled conditions. It remains unknown whether higher shell concentrations could yield a higher CO₂ compensation potential, or how this potential might change in coarser sediments with deeper oxygen penetration. Future research exploring a range of shell waste types, loadings and sediment grain sizes in realistic field conditions is needed to identify the optimal conditions for enhanced shell dissolution. Finally, future assessments should integrate a full life-cycle perspective on shell reuse pathways, including emissions associated with shell collection, processing, transport, and application, to determine the net CO₂ benefit of enhanced shell dissolution under realistic deployment scenarios.