Assessing the effective settling of mineral particles in the ocean with application to ocean-based carbon-dioxide removal

Ocean alkalinity enhancement (OAE), a potential approach for atmospheric carbon dioxide removal (CDR), can involve introducing milled mineral particles into the ocean to promote carbon dioxide uptake. The effectiveness of this method relies on particles remaining in the ocean mixed layer while dissolution takes place, which depends on particle settling rates. Conventionally, particle settling rates are assessed using the Stokes settling velocity in stagnant conditions. However, recent numerical modeling reveals that in dynamic, stratified ocean environments, sediment vertical transport can be up to an order of magnitude faster than Stokes settling because of two types of fluid instabilities that can take place at the mixed layer base. Here, we estimate effective settling velocities in the presence of these instabilities and assess the implications for the efficacy of this particular OAE approach for CDR. The new effective settling rate estimates are sufficiently rapid that there is negligible particle dissolution before particles settle out of the mixed layer. This result is independent of initial particle size for the range of sizes considered here. Findings underscore the importance of considering ocean dynamics and stratification in assessing particle settling rates and provide valuable insights for optimizing OAE applications in diverse marine settings.


Background and motivation
It is projected that global warming exceeding 1.5 • C will likely occur in the coming decades (IPCC 2018).To prevent the Earth's climate from surpassing this critical threshold, substantial reductions in carbon emissions and effective carbon dioxide removal (CDR) strategies need to be implemented now.Ocean alkalinity enhancement (OAE) is a potential CDR strategy that can involve adding natural minerals (e.g.olivine and basalt) or industrially produced substances (e.g.magnesium oxide and calcium oxide) to the surface of the ocean.When the substance dissolves in the surface ocean, the ocean pH increases and the system equilibrates to result in a flux of CO 2 from the atmosphere to the ocean.The dissolution results in the storage of CO 2 in the ocean as stable dissolved bicarbonate and carbonate (Bach et al 2019, Fakhraee et al 2023).For the CDR to be effective, the substance added must remain in the surface layer of the ocean, which is in contact with the atmosphere.
In this paper, we describe fluid dynamics that are relevant for the application of OAE where it is essential to understand the behaviour of mineral particles in a dynamic marine environment.In particular, we develop the results of a recent numerical modeling study (Yang et al 2023) to establish bounds on how long particles remain in the ocean surface layer in context with dissolution timescales to provide a refined understanding of the potential effectiveness of this type of OAE.We recognize that there are significant obstacles to scaling-up this CDR strategy, such as the substantial energy requirements to produce the particles (e.g.Hangx and Spiers 2009), and the potentially negative side effects on ocean ecosystems (e.g.Bach et al 2019); discussion of these is beyond our focus here.
The simplest description of gravitational settling of an individual particle in a fluid is the Stokes settling velocity.This describes the terminal velocity under gravity g of a sphere of density ρ p and diameter D, given by where µ is the dynamic viscosity and ρ is the density of the fluid.We take parameters corresponding to the ocean (µ = 10 −3 kg m −1 s −1 and ρ = 1000 kg m −3 ), as well as ρ p = 3200 kg m −3 , which corresponds to the density of olivine (Köhler et al 2013).Further, we take particle diameters in the range D = 4 − 40 µm (this is within the size range considered for OAE applications; Köhler et al 2013, Renforth andKruger 2013).
The Stokes settling velocity in addition to some estimate of ocean mixed-layer depth have been used to quantify how quickly particles will settle out of the surface ocean for the purposes of evaluating potential carbon draw down (Köhler et al 2013, Renforth andKruger 2013).This calculation, however, presumes a stagnant, unstratified water column.In reality, any marine environment in which scaled-up OAE experiments would take place might include stratification, ocean currents (e.g.vertical velocity gradients, or shear), waves and tides.We analyze the numerical simulations of Yang et al (2023) that examine how the former two features, stratification and velocity shear, modify particle settling from the canonical Stokes settling.We assume a uniform size distribution of particles, and that the concentration is initially well-mixed in the surface ocean.In practice, a well-mixed condition might be achieved by implementing an OAE experiment via injection of a sediment-laden flow, similar to an estuarine river inflow (e.g.Lu et al 2022).The appropriateness of the well-mixed assumption may depend on the specific method employed (e.g.spreading rock dust over the ocean surface vs. injecting a particle-laden flow), which would need to be assessed for any given approach.In addition, we do not take into consideration the effects of particle flocculation (i.e.particles bonding with each other to form larger clusters), although this may be an important factor depending on the particle types and sizes (e.g.Rouhnia andStrom 2015, 2017).

Particle settling in a dynamic marine environment
We consider two common and distinct fluid flow types (instabilities) that can influence the amount of time that particles remain in the surface mixed layer of the ocean (the basic flows are schematized in figure 1).The first can take place in an initially stratified upper ocean (e.g. the ocean mixed layer separated from the deeper waters by stratification across its base) to which a concentration of particles is added; as particles sink below the base of the upper layer, an unstable density profile can develop.The result can be the development of the Rayleigh-Taylor (RT) instability, or the sudden onset of dense sedimentladen sinking plumes.This instability is observed over a wide range of particle concentrations in numerical simulations (e.g.Yu et al 2014, Burns andMeiburg 2015), laboratory experiments (e.g.Blanchette andBush 2005, Davarpanah Jazi andWells 2016), and field observations (e.g.Warrick et al 2008, Hill andLintern 2021).The RT instability, manifest as plumes with characteristic mushroom caps (figure 1(a)), can play a significant role in transport of particles out of the surface mixed layer such that the instability gives rise to a downward movement of particles that is much faster than the Stokes settling velocity of individual particles (e.g.Bradley 1965, Yu et al 2014).
Horizontal velocity shear (e.g.arising from currents in the surface mixed layer overlying a relatively quiescent deeper ocean) is a common feature of the global oceans and sets up conditions that could give rise to another type of fluid instability, the Kelvin-Helmholtz instability (KH), characterized by roll-up billows (figure 1(b)) that can overturn a density interface (Thorpe 1971, Smyth andMoum 2012), and can be observed at the base of the ocean mixed layer (e.g.Kaminski et al 2021).As for the RT instability, the KH instability also influences how particles settle.Furthermore, the interaction of the two fluid instabilities (e.g.see the schematic depiction in figure 1(c)) (has important consequences to the effective particle settling rates (Yang et al 2023).
The general structure of the types of fluid flows described above may be seen in representative numerical simulations of a two-layer system that mimics an ocean surface layer overlying the deeper ocean (figures 2(a)-(d)); the simulation experiments are described in detail by Yang et al (2023).In these experiments, the density stratification is set up as a salt concentration, with the upper fluid layer being fresh and the lower fluid layer initialized to have a salinity contributing to a density difference of 20 kg m −3 .This is representative of a layer of freshwater entering an estuary, for example; the main results described here are insensitive to the precise value of this initial stratification.The upper layer is initialized with a particle concentration of 15 g l −1 , which yields an initial density difference between the upper and lower layers (including the salinity contribution) of 10 kg m −3 .For useful context with OAE, this can be related to some capacity for atmospheric CO 2 draw down.It is estimated that the dissolution of 1 tonne of magnesium olivine has the potential to absorb 1.25 t of CO 2 (Hangx and Spiers 2009).Given a concentration of olivine particles of 15 g l −1 distributed over a 1 m thick ocean layer of area 100 km square, the maximum potential for carbon draw down would be about 0.2 GtCO 2 (cf anthropogenic CO 2 emissions of 36.8GtCO 2 yr −1 ).A particle concentration of 15 g l −1 is comparable to olivine concentrations In the absence of velocity shear (figures 2(a) and (b)), one can see the early onset of the RT instability in the simulation (figure 2(a)) manifest as multiple rising and sinking finger-like plumes at the base of the surface layer.The sinking plumes are sediment-laden and transport particles from the upper to the lower layer.At later times, the plumes grow larger and form heads with mushroom caps (figure 2(b)).Note that the time for the onset of the RT instability is on the order of 1 min.In this time, the vertical distance that particles sink purely via the Stokes settling velocity (i.e.equation ( 1)) is only about 1 cm for the particles considered here.Therefore, the RT instability dominates the particle vertical transport.In the example simulation where there is velocity shear (figures 2(c) and (d)), the onset of the KH instability is seen at early times as the presence of a small roll-up billow structure at the interface between the two fluid layers (figure 2(c)).At this stage, there is negligible particle transport out of the surface layer.At later times, both the RT and KH instabilities are observed in the sediment concentration field, with the RT instability appearing to dominate in the form of tilted sedimentladen plumes in the lower layer (figure 2(d)).At the later stages of both experiments (figure 2(b) and (d)), one can visualize the rapid transport of particles out of the surface layer, carried by these plumes, and an effective settling of particles that is much faster than the Stokes settling.

The fate of particles in the ocean mixed layer
A useful metric for particle transport for the flows described above is the ratio between an effective settling velocity and the Stokes settling velocity of individual particles (U eff /U st ).The effective settling velocity is calculated as the vertical sediment flux normalized by the sediment concentration (e.g.Yu et al 2014, Burns andMeiburg 2015).Sediment fluxes increase rapidly as soon as RT plumes form, and then plateau when the plumes are fully developed; Yang et al (2023) quantify this stabilized effective settling velocity for a range of parameters.We note that the initial strength of the stratification has only a minor influence on the effective settling velocity ratio, with U eff /U st being most sensitive to the particle properties (see Yang et al 2023, for details).
In the absence of shear, U eff , enhanced by the RT instability, is up to 100 times faster than the Stokes settling for the range of particle sizes considered here (figure 2(e)).The ratio, U eff /U st , decreases with an increase in the size of particles.For a typical mixed layer depth of ∼150 m, particles of diameter D = 4−40 µm remain in the mixed layer on the order of 24 h to 6 h (compared to 3 months to 22 h for Stokes settling).This suggests that the residence time (a day or less) is relatively insensitive to particle size when particles are transported via plumes arising from the RT instability.
The presence of velocity shear modifies the structure of the RT plumes (figure 2(d)), and therefore the effective settling of particles, while the effective settling velocity remains faster than the Stokes settling for the range of particle sizes considered (figure 2(e), gray curve).For small particle sizes, the KH instability grows rapidly and suppresses the RT instability; the effective settling velocity is reduced compared to the case without shear.An interesting feature of flows with larger particles, and one that illustrates the complexity of particle settling in the real world, is that the RT instability inhibits the KH instability and the effective settling velocity is enhanced simply due to the additional energy input by shear (figure 2(e), gray curve); see Yang et al (2023) for further discussion.
A key question of relevance to CDR is the time that particles remain in the ocean surface mixed layer relative to their dissolution timescales.Renforth and Kruger (2013) applied a shrinking core model to quantify the fractional change in diameter of magnesium oxide particles (assumed to be spherical) in the ocean, considering typical parameters for the dissolution of MgO in seawater.After some time t, the remaining fraction of reacted material is expressed as where a value of W r = 10 −9 kg m −2 s −1 is taken for the dissolution rate.Given the significant uncertainty in W r (e.g.Fuhr et al 2022, Fakhraee et al 2023), and the simplicity of the model, the result should only be considered as an order of magnitude estimate.We take the mixed-layer depth to be 150 m and consider particles settling at the Stokes settling velocity, with a time-dependent particle size due to dissolution, and compute the fraction of a particle remaining at the time it leaves the mixed layer (figure 2(f)).
The model indicates that particles initially less than about 10 µm in diameter dissolve entirely before settling out of the mixed layer.Conversely, for particles with initial diameters larger than about 15 µm settling under the Stokes settling velocity, there is effectively negligible dissolution by the time particles leave the mixed layer (figure 2(f), blue curve).This becomes the case for the full range of initial particle sizes considered (4-40 µm) if particles settle at the enhanced effective settling velocity described here, having significant implications to the effectiveness of this type of OAE.
There are other fluid-dynamical complexities not considered in the dissolution model (equation (2)).For example, Amann et al (2022) found that larger particles can have proportionately higher weathering rates compared to smaller particles; this is because fluid can more easily move through the larger and interconnected pores of the larger particles.More complex dissolution models are required to incorporate such processes which are also a function of the material.

Implications
The broad results described here highlight that assessing the viability of this type of OAE must account for both the shear and stratification properties of the ocean to estimate an effective settling velocity because the Stokes settling rate may not apply.Results further suggest the optimal regions for field tests of this form of OAE would be where the mixed layer is moving faster than the deeper ocean (and this holds in large parts of the world's oceans); this would yield particles having a longer residence time in the mixed layer than in the absence of shear.In addition, smaller particles are preferable as they allow for the growth of the KH instability (that suppresses vertical particle fluxes) as opposed to large particles that favor the dominance of the RT instability.
The occurrence of the convective instabilities described here is contingent upon specific conditions.In a quiescent marine setting, convective instabilities have been shown to occur when particle concentrations are dilute (volume fractions < 1.5%) and particle diameters are on the order of 1−100 µm (see e.g.Warrick et al 2008, Yu et al 2014, Davarpanah Jazi and Wells 2016).Small particles (< 1 µm) can remain in suspension because Brownian motion can counteract gravity and the convective instability never develops (e.g.Cosgrove 2010, Monroy et al 2017).In a dynamic marine setting, convective instabilities may be impeded depending on the strength of the stratification relative to the strength of the horizontal velocity shear (the ratio of these quantities is the Richardson number).Flows characterized by Richardson numbers smaller than 1/4 are prone to KH instabilities.This regime is frequently observed throughout the global oceans (e.g.Smyth 2020, Hughes et al 2021, Kaminski et al 2021).
A next step for global and regional models (that generally do not explicitly capture the fluid instabilities described here) would be to formulate parameterizations for effective settling of particles that are functions of stratification, shear and particle size.This would allow for an exploration of optimal particle sizes depending on the dynamics at play in different ocean settings, for example.This will require further study to quantify the influence of the shear strength on the effective settling velocity of particles.Other relevant factors for future research include examining the influence of tidal flows, surface waves, and particle flocculation to improve our understanding of particle transport for OAE applications in a variety of dynamic marine settings.

Figure 1 .
Figure 1.Schematic diagram depicting sediment-laden fluid flows in a dynamic marine setting such stratified upper ocean.In a field implementation of OAE, milled particles are put into the ocean at the surface and ultimately are transported out of the surface mixed layer by fluid flows that depend on the strength of the stratification, the particle size and density, and the difference in velocity (i.e.shear) between the mixed layer and the underlying water.Schematics on the left depict the basic structure of two types of flows: (a) the Rayleigh-Taylor instability that can arise in the presence of stratification and (b) the Kelvin-Helmholtz instability that can arise in the presence of shear.In general, both types of flow can be at play (c) and influence how long milled particles remain near the surface in any given marine setting.

Figure 2 .
Figure 2. Numerical simulation results showing the evolution of the particle concentration in the absence of background velocity shear (a), (b), and in the presence of velocity shear (c), (d).Early representative times for each experiment are shown in the top row (a), (c), and later times in the bottom row (b), (d).In the absence of background shear, the Rayleigh-Taylor instability is responsible for the downward transport of particles.In the presence of shear, both the Kelvin-Helmholtz and Rayleigh-Taylor instability appear in the flow.(e) Ratio of effective particle settling velocity to the corresponding Stokes settling velocity (U eff /Ust) as a function of particle diameter (D; corresponding Ust is shown on the top axis) with and without background velocity shear (gray and black curves, respectively).(f) The remaining fraction χ of a dissolving particle at the time it leaves the mixed layer.