Introduction

In order to achieve stringent climate goals, such as limiting global warming to 1.5 °\(\mathrm{C}\), hundreds of gigatonnes (Gt) of CO2 need to be captured and stored1,2. To achieve this objective, a wide range of carbon dioxide removal (CDR) technologies would have to be developed and deployed at scale2,3. Among them, enhanced rock weathering (ERW) is a competitive option because it combines direct removal of atmospheric CO2 with long-term storage through conversion into aqueous alkalinity or carbonate minerals3,4,5,6,7,8,9,10,11,12,13,14. Theoretical estimates have shown that the weathering of 1 tonne of basic (e.g. basalts) and ultrabasic rocks (e.g. olivine) can remove ~ 0.3 and 0.8 tonne of CO2, respectively8, indicating a considerable overall theoretical potential4,8,15. ERW can also mitigate ocean acidification10, and has an intermediate cost compared to other CDR technologies3,12.

ERW in soils has a series of additional benefits16. It takes advantage of soil chemistry, e.g. elevated local CO2 concentration and the presence of organic ligands from plant and microbial activities, to accelerate the weathering reactions, while not competing with other land uses7,9,17. Application of ground rocks to ~ 50% of global croplands would potentially sequester up to 2 Gt of CO2 per year5. ERW in soils can improve soil quality by reducing soil acidity, which has been traditionally achieved by the practice of liming18,19, while avoiding CO2 emissions associated with lime application6. ERW can also release macro and micro nutrients and thereby reduce the need of soil fertilization with fertilizers, thus improving productivity across agroecosystems12,20,21,22.

However, for ERW to be deployed at a meaningful scale, the primary challenge is to achieve and maintain high weathering rates and thus carbon removal efficiency. In current analyses, weathering rates have been identified as a major source of uncertainty10,12. In situ weathering rates in the soil environment are particularly variable owing to the complex nature of the biogeochemical reactions in the soil matrix and their interactions with the applied rock materials15,23.

Theoretical weathering rates that are based on batch experiments or reaction rate laws can be helpful14,24,25, but one of the caveats is that they do not account for factors such as transport limitation. Transport limitations to reactants and reaction products are a major constraint on achievable reaction rates, resulting in observed weathering rates that are orders of magnitude lower than the theoretical estimates19,26. This discrepancy is particularly relevant for the optimization of particle size of rock materials for field application given that the energy penalty of grinding and comminution to create more surface area is significant27. The complex feedback between characteristics of rock materials, soil matrix, and transport is further emphasized by the observation that experimental weathering rates may increase when application rates decline, as shown for olivine6. Accordingly, the overarching goal of this work is to establish (i) a mechanistic understanding as well as (ii) quantitative approximations of the extent to which both, biogeochemical reactions and transport, are affected by different environmental controls.

While existing dedicated experiments on enhanced rock weathering in soils are few, and more (laboratory and field) experiments are needed to uncover the mechanisms underlying in situ weathering rates, reactive transport modeling provides a valuable tool to complement experimental studies5,28. It has the advantages of expanding the spatial scales and extending beyond the duration of current experimental investigations which are typically carried out on monthly time scales, and thus to better inform field studies and engineering practices. It can also help assess potential environmental risks that can arise from the weathering products (e.g. heavy metal)7,11, and quantify the amount of CO2 that is captured and stored, because reactive transport modeling can readily track the fate of different chemical species in the system. However, modeling studies to date rarely consider reaction networks that reflect the complex soil chemistry.

To explore the major controlling chemical-physical factors of weathering rates in the soil environment, we used a 1D reactive transport model that considers a reaction network capable of capturing microbe-mediated reactions that control the soil chemistry29,30, and different unsaturated hydrological scenarios of the soil environment. There are numerous possibilities regarding the choice of mineralogy for EWR applications. Given our interest in probing the potentials of the ERW practice, we elected to evaluate the fate of a mineral phase that can be considered as particularly effective in producing alkalinity/carbonate rocks: the olivine group mineral forsterite. To constrain model complexity and to increase our ability to detect fundamental relationships, we limited the number of secondary solid phases, aware that this choice constitutes a likely simplification of outcomes in the natural weathering environment. As for system context, our work leverages information (e.g. composition and chemistry) for a silt loam and a sandy loam from a previous study on groundwater recharge in the Central Valley area of California30. These sites also convey practical relevance to our study given the interests of California in carbon neutrality31 and the need to employ different soil types for meaningful carbon removal via ERW.

Methods

Model description

The multi-component multiphase reactive transport code TOUGHREACT32 was used to simulate a 2-m long vertical 1D soil column under a range of conditions. It was assumed that there were no gas pressure or temperature gradients, which was done by solving Richards’ equation (i.e. module EOS9 of TOUGHREACT33). Diffusive transport of the gas phase was considered, as the partial pressure of CO2 in biologically active soil is different from that in the atmosphere. Kinetic mineral reactions were calculated using the Transition State Theory rate law24. The rates of microbe-mediated reactions that are important in the soil environment were tracked by a general rate law based on the Monod equation. Cation exchange was also included following the Gaines-Thomas convention. Details of the model are documented in the SI.

The computational domain was discretized into 36 grid cells. The first grid cell at the top close to the atmosphere is 2 mm, and the cell size increases gradually by ~ 20% consecutively towards the bottom, reaching 100 mm in grid cell 22. The size of grid cell 23–36 was set at 100 mm. At the top of the computational domain, the gas pressure was set equal to the atmospheric pressure. A constant water flux was introduced at the top to simulate infiltration. Water-table conditions were imposed at the bottom by setting the water saturation to one and the total pressure to atmospheric pressure.

Simulation setup

Simulations that mimic soil column experiments with and without the application of silicate rock materials were performed, to evaluate the weathering rate of the applied rock materials and the resulting changes in soil chemistry. The applied material is forsterite (Mg2SiO4), an end member of the olivine group minerals.

The silt loam and the sandy loam, represent a less well-drained and a well-drained soil (see hydraulic properties in Table S1), respectively. Two infiltration rates (200 and 2000 mm/yr) were considered to bracket common scenarios5,34. Together, four water content profiles were generated (Fig. S1).

The two soils are composed of primarily quartz, some montmorillonite, and small fractions of K-feldspar, ferrihydrite and calcite (Table S2). Five secondary minerals, amorphous silica, illite, kaolinite, magnesite and siderite were included. The first three are common products of rock weathering, and magnesite and siderite are potential carbon sinks given the presence of magnesium and iron. Magnesite is the most stable magnesium carbonate but forms slowly under ambient conditions. But the carbonation rate can be significantly enhanced by organic additives and microbes35,36, as can be encountered in the soil environment. Therefore, magnesite was included here instead of other faster forming but metastable MgCO3 to represent a scenario in which carbonation is enhanced. The kinetic coefficients for all the mineral reactions are summarized in Table S337.

In our simulations, the soil has an equivalent specific surface area of 10 m2/g38, and was divided among the primary minerals and secondary minerals based on their respective volume fraction. The total exchange site density was set to 10 meq/100 g soil for both soil types39. Chemical compositions of the initial and infiltrating water are detailed in Table S5.

The solid organic carbon pool in the soil was represented by cellulose. It was modeled to degrade into acetate following the reaction:

$$C{H}_{2}O\left(s\right)\leftrightarrow 0.5C{H}_{3}CO{O}^{-}+0.5{H}^{+}.$$
(1)

This reaction was maintained at equilibrium to provide a constant source of acetate. A network of kinetic reactions (Table S4) was used to simulate the degradation of acetate and thus, indirectly, the kinetic decomposition of cellulose and other important microbe-mediated reactions in the soil. This representation of soil chemistry was initially developed and tested for a site at Rifle Colorado29, and later adopted in the study of California Central Valley30. In our study, the kinetic coefficient for the microbial oxidation of acetate by oxygen (\({k}_{ox}\)), which is the primary reaction (Eq. (2)) that controls the soil CO2 profile, was re-calibrated for our hydrological scenarios to capture key features of CO2-depth profiles in typical field observations40.

$$C{H}_{3}CO{O}^{-}+2{O}_{2}\to 2HC{O}_{3}^{-}+{H}^{+}.$$
(2)

We considered a onetime application of 160 tonnes forsterite per hectare (i.e. ~ 16 kg/m2), which is equivalent to a surficial application of 1 cm of rock materials with a porosity of 50%. This application rate is also comparable to previous numerical and experimental studies5,6,17,41, and is within the range that would be needed if olivine weathering were to be leveraged for addressing phosphorous deficiency21. Two additional application depths (15 and 50 cm) were simulated, with the same total forsterite amount of 160 tonnes per hectare. In these two cases, the applied rock was mixed with the soil matrix uniformly within the application layer. It was assumed that the soil hydraulic properties were not affected by the mixing.

For the four combinations of hydraulic conditions, three rock application scenarios and a baseline case without rock amendments were simulated. For all cases with rock amendments, three forsterite surface area values (0.1, 1, and 10 m2/g, corresponding to grain sizes from several hundreds of microns to microns5,42,43) were compared. For simulations with an application depth of 15 cm and a specific surface area of 10 m2/g for forsterite, two elevated \({k}_{ox}\) values (10× and 100×) were also tested. For a subset of the cases (i.e. forsterite surface area = 10 m2/g), magnesite precipitation was suppressed by reducing its default surface area of 10 m2/g by 1000 times. Simulations were run to cover a total period of 5 years, which has been used as a reasonable timeframe for comparing weathering rates and analyzing carbon removal efficiency19.

Results and discussion

Changes in soil chemistry from forsterite weathering

With the application of forsterite, evident changes were observed in both the mineral phases and the aqueous composition in all cases, compared to the corresponding baseline simulations without forsterite application (Fig. 1). The major driver is forsterite dissolution, which resulted in elevated pH, and silica and Mg concentration in solution:

Figure 1
figure 1

Vertical profiles at year five for an application of forsterite with a specific surface area of 10 m2/g in the top 15 cm for both soil types at both infiltration rates. (a–d) Changes in the volume fractions of key minerals (from left to right: forsterite, magnesite, calcite, amorphous silica), and (e–h) soil water chemistry (from left to right: pH, alkalinity as Total Inorganic Carbon (TIC), DOC as acetate) and CO2 partial pressure [bar].

Full size image
$$M{g}_{2}Si{O}_{4}\left(s\right)+4{H}^{+}\leftrightarrow Si{O}_{2}\left(aq\right)+2M{g}^{2+}+2{H}_{2}O.$$
(3)

This reaction drives the precipitation of magnesite (MgCO3) and amorphous silica (SiO2(am)) (Fig. 1b,d), which were not observed in the baseline simulations. Given the large surface areas and fast reaction rates used, the precipitation occurred primarily in the region where Mg2+ and SiO2(aq) were produced by forsterite dissolution. For calcite, montmorillonite and k-feldspar, both dissolution and precipitation were observed, whereas in the baseline simulations precipitation was absent or in smaller amounts (Fig. S2). Precipitation of illite, kaolinite and gibbsite either did not occur or was reduced because of the higher pH. No siderite precipitation was observed in any of the cases simulated.

Soil water pH was increased by about one unit, e.g. from 7.0 ~ 8.5 to 8.0 ~ 9.5 in the silt loam at the high infiltration rate. The extent of increase is comparable to what was observed in previous experimental studies, e.g. from 7 to 8.5 in Amann et al.26, from 3.55 to 4.69 ~ 5.18 in Dietzen et al.6, and from 5.0 ~ 6.5 to 6.5 ~ 7.5 in Haque et al.22. The increase in pH led to an increase for about one order of magnitude in acetate concentration, which provides a measure of the total dissolved organic carbon (DOC) in our simulations, because it is produced by cellulose degradation that is pH-dependent. Partial pressure of CO2 (\({P}_{C{O}_{2}}\)) shows a clear decreasing trend in the top layer as CO2 was consumed by mineral reactions, e.g. the precipitation of primarily magnesite and to a lesser extent calcite. The bottom section still shows an increasing trend as in the baseline case, but at values that are about one order of magnitude lower. Total inorganic carbon in the solution (TIC) is in general lower because of the transfer of carbon from solution to magnesite and calcite following the carbonation reaction, and is consistent with the lower \({P}_{C{O}_{2}}\) and higher pH.

The impact of soil water content on forsterite weathering

Figure 2 summarizes the total mole change of forsterite, magnesite, and calcite within the soil column by year five for the four water content cases with the largest forsterite surface area, and shows generally more weathering in well-drained soils (i.e. sandy loam) and at the low infiltration rate, i.e. under conditions generating lower water content.

Figure 2
figure 2

Mole changes of (a) forsterite, (b) magnesite, and (c) calcite for the silt loam (SiL) and the sandy loam (SaL) at the infiltration rate of 200 (_L) and 2000 (_H) mm/yr with different application scenarios for a specific surface area of 10m2/g for forsterite. Application depth: green – 1 cm, blue – 15 cm, coral – 50 cm (the total amount of the rock material applied is the same across the three application scenarios).

Full size image

For an application depth of 15 cm, more than 70% of the applied forsterite (~ 120 mol assuming a cross-section area of 1 m2) was weathered in the sandy loam. The weathering rate was slightly higher at the low infiltration rate. In contrast, about 35% and 15% of the applied forsterite was weathered in the silt loam at the low and high infiltration rate, respectively. The average dissolution rates ranged between ~ 5 × 10–12 mol/m2s to ~ 2 × 10–11 mol/m2s, in comparison to the intrinsic kinetic coefficients (2.29 × 10–11 mol/m2s for the neutral reaction mechanism and 1.41 × 10–7 mol/m2s for the acid reaction mechanism, Table S3). The amount of magnesite precipitation tracks forsterite dissolution, being about twice of the amount of forsterite following the stoichiometry of Mg in the two mineral phases. The variability of calcite precipitation across the four hydrological scenarios is similar to that of magnesite, except for the sandy loam cases, in which more precipitation was observed at the high infiltration rate. Calcium consumed by the precipitation reaction was primarily sourced from cation exchange. For instance, in the sandy loam at the low infiltration rate, the fraction of exchange sites occupied by Ca was reduced from ~ 0.4 to 0.2 in the topsoil where calcite precipitation occurred.

In comparison, with an application depth of 50 cm, forsterite weathering and magnesite precipitation show similar variations across soil types and infiltration rates, but at smaller amounts. There was more calcite precipitation, as precipitation also occurred in deeper sections of the soil column.

For the surficial application scenario (1 cm), forsterite weathering was higher in the silt loam at the low infiltration rate, showing complete weathering. In the sandy loam, the low water content and large forsterite surface area and thus fast forsterite dissolution created large local supersaturation with respect to amorphous silica. As a result, with the same simulation setup, amorphous silica precipitation led to a volume fraction change that is unrealistically high. These data points were hence excluded from the analyses.

Faster forsterite weathering under conditions generating a lower water content implies that atmospheric CO2 must be the dominant CO2 source in these simulations. CO2 available for the weathering reactions in the soil environment comes from three sources: CO2 that is dissolved in the infiltrating water, atmospheric CO2 transported into the soil matrix, and CO2 generated by organic matter (OM) decomposition. The first source accounts for a small fraction and is larger at the high infiltration rate. The second source is influenced by gas transport in our simulations, and the third source is controlled by the reaction network and soil chemistry. Given the reaction network considered, more biogenic CO2 should be generated at higher pH. In contrast, gas diffusion—which follows a power law relationship with respect to gas saturation—is faster when water content is lower (as in the sandy loam).

The CO2 profiles that show decreasing partial pressure in the top layer where forsterite was applied suggest a downward gas transport from the atmosphere. Even though the gradients in the lower section indicate an upward CO2 transport, the supply from soil decomposition was smaller than what would be expected based on the baseline simulations because of the elevated soil water pH and thus lower biogenic CO2.

The observation that with the high forsterite surface area, surficial application resulted in the fastest weathering rate is consistent with the conclusion that under these conditions, atmospheric CO2 was the major source driving the weathering reactions.

The impact of soil chemistry on forsterite weathering

The observations above assumed that the functions of plant and microbial activities governing OM decomposition were neither promoted nor suppressed in response to the application of silicate rock materials and the subsequent changes in the soil chemistry, and the kinetic coefficients for the reaction network describing soil decomposition were fixed. However, studies have shown positive responses in functional and species richness of microbes to elemental and pH changes from enhanced rock weathering44,45. To investigate this potential scenario and the impact of the resulting soil chemistry on forsterite weathering, \({k}_{ox}\) was increased for the application depth of 15 cm.

Increasing \({k}_{ox}\) and thus CO2 partial pressure in the soil column resulted in an increase in forsterite weathering across all hydrological scenarios (Fig. 3a). The extent of increase is especially significant for the less well-drained soil (i.e. the silt loam), for which the availability of atmospheric CO2 was more constrained. With a 100× increase in \({k}_{ox}\), the applied forsterite was completely weathered by year five, except for the case of sandy loam at the low infiltration rate. This is likely caused by the fast gas transport under these conditions, which can limit the accumulation of CO2 in the soil column. In the sandy loam cases at the high infiltration rate, increased \({k}_{ox}\) resulted in higher rates of forsterite dissolution. This dependence on water flow indicates that in these cases, leaching of solutes rather than CO2 transport can be important in controlling the weathering rate34. In the silt loam at the low infiltration rate, increasing \({k}_{ox}\) by a factor of 10 resulted in weathering rate that is comparable with that of surficial application, indicating that CO2 from organic matter decomposition can make up for the constrained atmospheric CO2 supply in the deeper soil layers.

Figure 3
figure 3

Mole changes of forsterite for the silt loam (SiL) and the sandy loam (SaL) at the infiltration rate of 200 (_L) and 2000 (_H) mm/yr (a) with different acetate oxidation rates (\({k}_{ox}\)) for a forsterite specific surface area of 10m2/g with an application depth of 15 cm, and with a forsterite specific surface area of (b) 1m2/g (diagonal grid pattern) and (c) 0.1 m2/g (vertical line pattern) with an application depth of 1 cm (green), 15 cm (blue), and 50 cm (blue).

Full size image

The impact of surface area on forsterite weathering

The simulations presented above used a specific surface area of 10 m2/g, which corresponds to very fine grain sizes of the applied forsterite (~ 0.2 µ diameter assuming uniform spheres). This, however, would require considerable energy input (~ 700 MJ/tonne)5, which depending on the carbon intensity of the energy used may not be the optimal operational option. Here, we evaluate how variations in the specific surface area and thus grain size of the applied forsterite affects its weathering in the presence of transport constraint and microbe-mediated reactions. The average forsterite dissolution rates for simulations using different specific surface areas under all four hydrological scenarios ranged between 5 × 10–13 and 2 × 10–11 mol/m2s. They are comparable to previous pot experiments19, which used olivine grains that is composed of > 80% of forsterite with specific surface areas of 0.03–3 m2/g and measured dissolution rates between 4 × 10–13 and 3 × 10–12 mol/m2s.

When the specific surface area was increased from 0.1 to 1 m2/g (~ 2 and 20 µ grain size assuming uniform spheres, respectively), the weathering rates across all cases evaluated increased, but the extent of increase was not always proportional (Fig. 3b,c). For instance, in the sandy loam, weathering rate increased by a factor of 7.8–8.4 at the high infiltration rate, and by a factor of 8.7–9.4 at the low infiltration rate. This rate increase was much more moderate in the silt loam, being 3.4–6 and 4.4–8.3 at the high and low infiltration rate, respectively. Increasing the surface area of forsterite from 1 to 10 m2/g attenuated the weathering rate increase even further, with a factor of 1.2–1.4 in silt loam at the high infiltration rate and a factor of 4.1–5.7 in the sandy loam at the low infiltration rate. This brackets what has been reported in previous studies46.

With a specific surface area of 0.1 m2/g, forsterite weathering did not vary significantly across different application depths. This implies that CO2 transport was not limiting. Forsterite weathering was faster in the silt loam and at the high infiltration rate instead, as water is needed to keep the mineral saturation state low and thus the thermodynamic driving force high. As the specific surface area increased to 1 m2/g, the controls of CO2 transport became evident, showing similar variations across different application depths in the silt loam to the highest specific surface area simulations reported in previous sections. The increase in the relative magnitude of mineral reaction rates and transport rates, i.e. transport becomes increasingly limiting, also explains why the return of increasing surface area diminishes at higher surface area, as quantified above.

These observations highlight the necessity of considering the interplay between reaction kinetics and transport in the selection of operational parameters. At sites where transport of CO2 from the atmosphere is not or less limiting (e.g. the sandy loam), increasing surface area has a higher return in terms of improving weathering rate, while at sites where CO2 availability is already limiting (i.e. the silt loam), the improvement would be more moderate. In those cases, it would make more sense to reconsider the siting or to introduce microbial interventions rather than continuing with increasing the surface area, i.e. reducing grain sizes.

Carbon accounting

Tracking the transport and fate of carbon is critical for carbon accounting and crediting. In our simulations, other than the exchange with the atmospheric CO2 at the top, all other carbon pools and fluxes were explicitly tracked. The simulations used a strict tolerance (10–12) to ensure mass conservation. Therefore, we can monitor the amount of carbon sequestered and the net carbon exchange with the atmosphere (Fig. 4a). Note that the organic carbon pool was highly simplified in our study and the changes in the organic carbon pool were not meant to capture sequestration via soil organic carbon, but to help break down the sources of carbon that was sequestered into the inorganic carbon pool.

Figure 4
figure 4

The amount of inorganic carbon sequestered [mol] in relation to (a) the amount of forsterite weathered [mol] and (b) the contribution from the atmosphere [mol]. The application depth is color coded: green – 1 cm, blue – 15 cm, and orange – 50 cm. The light and dark color corresponds to low and high \(SS{A}_{for}\), respectively. The black edge indicates simulations with suppressed magnesite precipitation. The symbols indicate the hydrological scenarios: filled circle - silt loam at 200 mm/yr, filled asterisk - silt loam at 2000 mm/yr, filled diamond - sandy loam at 200 mm/yr, and filled square - sandy loam at 2000 mm/yr. The symbol size is enlarged by a factor of 2 and 3 for simulations using \(10{k}_{ox}\) and \(10{0k}_{ox}\), respectively.

Full size image

The inorganic carbon (IC) sequestered measures the net change in the IC pool in the soil column and the amount of dissolved IC that enters the groundwater system at the bottom, and was calculated as the difference between the simulations with and without forsterite application. The amount of forsterite weathered and the amount of IC sequestered follows a linear relationship with a 1:2 ratio. This relationship is consistent with the stoichiometry of forsterite dissolution and magnesite precipitation.

$$M{g}_{2}Si{O}_{4}\left(s\right)+2C{O}_{2}\leftrightarrow 2MgC{O}_{3}+Si{O}_{2}\left(aq\right).$$
(4)

However, this ratio is not dependent on magnesite precipitation. Data points from simulations with suppressed magnesite precipitation (points with black edges in Fig. 4) follow the same trend, i.e. under these simulation conditions CO32– remains in the solution as the major form of carbon and as the major anion that balances the positive charges of the cations. With magnesite precipitation suppressed, Mg accumulation in the aqueous phase reduces the thermodynamic driving force for forsterite dissolution. Therefore, understanding the complex kinetics of magnesium carbonate precipitation under ambient conditions in the soil environment will be critical for evaluating ERW efficiency and engineering practices. IC sequestered is also attributed, to a lesser extent, to calcite precipitation and calcite dissolution avoided. The cumulative amount of carbon going into the groundwater was orders of magnitude lower than the amount contained in solid phases in our simulations, whereas under certain conditions CO2 may be stored primarily as HCO3- in solution without precipitation13.

For most simulations, IC sequestered was primarily sourced from atmospheric CO2 (Fig. 4b). This is consistent with the observation that weathering rates in these cases were subject to CO2 transport limitation. In a few simulations with large application depths or with increased acetate oxidation (\({k}_{ox}\)), carbon was sequestered by reducing the biogenic CO2 emission that would otherwise occur without the application of forsterite (as noted by the negative signs). In these simulations, forsterite weathering removed local soil CO2.

The biotic effects, as reflected by elevated \({k}_{ox}\) values led to shifts in the amount of forsterite weathering and the relative contribution of atmospheric CO2, without affecting the relationship between the IC sequestered and the weathered forsterite. However, the relationship may vary if the biotic reactions change the partitioning of carbon between the solid phase and the aqueous phase, which have different stoichiometries. Moreover, future investigations are needed to fully assess the organic carbon sink in the soil and its relation with rock weathering, among which mineral-organic associations are a topic of great importance47.

Implications for engineering design

While complete forsterite weathering (~ 120 mol) can occur within five years, which gives an equivalent carbon removal rate of ~ 2.3 kg CO2/m2/yr for a onetime application of ~ 16 kg/m2, much less carbon was sequestered under unfavorable conditions. Given that CO2 emission associated with crushing and milling of the rock materials can amount to ~ 10–25% of the CO2 sequestration capacity27, scenarios that result in less than 50 mol of IC sequestered after five years can hardly achieve carbon neutrality, let alone negative emissions. These are primarily simulations with the lowest specific surface area for forsterite. This is consistent with the previous observation that using grain size smaller than 10 µm is necessary to achieve significant CO2 removal48.

The weathering rate of a given rock material is controlled by the accessibility of the reactants, primarily the reactive rock surfaces and CO2 availability at the rock surface (Fig. 5). When atmospheric CO2 is the major source, the weathering rate is constrained by the smaller of the chemical reaction rate and the transport rate. Thus, while grinding the rock materials into fine grain sizes can enhance weathering, the extent of the increase is subject to the constraints on CO2 transport. When chemical reactions are relatively slow, increasing surface area by reducing grain size can boost the weathering rate significantly. Once the transport rate of reactants (e.g. CO2) to the rock materials or of the products away from the rock surfaces becomes limiting, increasing surface area has very limited effect on the weathering rate. It should also be noted that, as the grain size decreases, the energy penalty of grinding the particles down to finer sizes increases disproportionately27. Therefore, it is important to be able to identify an optimal grain size that balances the gain in weathering rate and thus carbon removal and the loss in energy due to comminution. The optimal grain size would depend on the crossing point of the chemical reaction rate and the transport rate. As illustrated in our simulations, for sites with high effectiveness of transporting CO2, the increase in weathering rate from increased surface area can be sustained for a wider range of surface area values, and thus the optimal grain size is likely to be smaller. Therefore, siting and operations that lead to effective CO2 transport will be important for ensuring ERW efficiency. For cases when grain size customization is impractical, other operations that relax the transport restrictions of CO2 at the sites should be considered, such as mixing or tillage that may change soil texture and hydraulic properties to improve CO2 transport. While our simulations only considered gas diffusion for CO2 transport, gas flow arising from complex hydrological processes provide another transport mechanism, and thus soil properties (e.g. textural heterogeneity), processes (e.g. transient flow and porosity–permeability change, see SI for porosity change observed in our simulations) and operations influencing gas flow need to be evaluated further.

Figure 5
figure 5

Schematic illustrating chemical and physical properties and processes that affect the efficiency of enhanced rock weathering in soils, and the impacts of engineering operations.

Full size image

Alternatively, high weathering rate can be achieved if local CO2 sources (i.e. CO2 from OM decomposition in the soil) can be tapped into, to relax the constraint of CO2 transport. This can potentially be realized by increasing the application depth to bring the reactive mineral closer to soil CO2. However, when OM decomposition cannot replenish CO2 consumed by the weathering reaction fast enough, a high weathering rate cannot be sustained. In this case, increasing the application depth becomes undesirable as it moves the rock material away from atmospheric CO2. In contrast, when local CO2 supply is ‘guaranteed’ by fast OM decomposition, increasing the application depth helps achieve the high weathering rate. It can also be expected that in these cases, the positive effect of increasing the surface area can be maintained. The observation that the optimal application depth can be highly variable suggests that decisions regarding application by airplane spreading versus slurry application or tillage should be based on site characteristics.

For local supply of soil CO2 to be consistently high, it is necessary that the microbial and plant activities respond positively to the application of rock materials and the resulting changes in soil chemistry. While there are still a lot of uncertainties regarding this topic, existing studies, though few, have shown promise44,45. In addition to providing CO2 for the weathering reactions, the biotic processes are also expected to enhance rock weathering and CO2 sequestration by increasing the availability of organic/chelating ligands that can bond with reaction products and by facilitating carbonation reactions49,50. Future studies that further our understanding of the underlying mechanisms associated with the interactions between silicate rock weathering and biotic factors and processes – such that the effects on weathering rates can be quantitatively evaluated for operation optimization – are in great need. Such studies will also help develop innovative approaches that enhance weathering rates via microbial intervention or other biotic processes.

It should be noted that there are other operational considerations in engineering designs, although not discussed here, warrant in-depth future investigations. One operational parameter that is anticipated to impact the achievable weathering rate in the field is the application rate and frequencies in addition to the application depth. Moreover, in regions where the mining industry is well established, grinding may not be costly, whereas transportation to the sites may account for a larger portion of the energy penalty51. In this case, spatial analyses are needed to optimize the match between suitable rock materials and the application sites.

Conclusions

Our study aimed to elucidate the impacts of coupled biogeochemical reactions and transport processes on the efficiency of enhanced rock weathering (ERW) in soils and quantify the effects of a set of primary environmental and operational controls. Based on a reactive transport model considering both unsaturated soil hydrodynamics and a validated soil reaction network, our results provide important insights for the practices of enhanced silicate weathering as soil amendments.

ERW in soils can be an effective option for carbon sequestration, which removes CO2 either from the atmosphere directly, or by reducing soil related emissions. Its efficiency is promoted by conditions that maintain high CO2 availability, e.g. good connectivity with atmospheric CO2 and/or sufficient biogenic CO2 supply. Therefore, soil texture, composition, and reaction networks are important factors to consider for proper siting. It is also important to co-optimize the site conditions with engineering operations. For instance, the applied rocks should be placed close to the dominant CO2 source (e.g. via optimizing the application depth), and the optimal grain size for efficient CO2 removal is site-specific and reducing grain size should be considered only when other kinetic controls (e.g. of carbonation reaction) and CO2 supply are nonlimiting. We further recommend future research to improve understanding of biotic-abiotic interactions at their native scales, such that these effects can be quantified and properly assessed to ascertain ERW efficiency in soils.