By Lindsey J. Heagy, Thibaut Astic, Joseph Capriotti, John Weis, and Douglas W. Oldenburg
Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, Vancouver, British Columbia
e-mail: [email protected]
Abstract
In addition to reducing anthropogenic emissions of CO2, it is increasingly clear we also need to remove CO2 from the atmosphere in order to avoid some of the worst case scenarios for climate change. Geologic sequestration of CO2 is among the most attractive approaches because of the large global capacity and long-time scales for storage. One mechanism of geologic storage is through carbon mineralization. Some mafic and ultramafic rocks contain minerals that will react with CO2 in a carbonation reaction and convert it to carbonated minerals. This is effectively a permanent CO2 storage mechanism. The geologic question we are faced with is then to locate, delineate and estimate the volume of potentially reactive rocks. Using a synthetic model that emulates a prospective site for carbon mineralization in British Columbia, we simulate and invert gravity and magnetic data to delineate reactive rocks. We begin by inverting each data set independently and introduce a proxy experiment to contend with the challenging problem of choosing an appropriate physical-property threshold to estimate volumes from the recovered model. We use this proxy experiment to estimate thresholds for standard, l2 inversion of the gravity and magnetics, as well as for inversions which use sparse and compact norms. A Petrophysically and Geologically Guided Inversion (PGI) framework is used to construct quasi-geologic models from which volumes can be estimated directly. We apply the PGI framework to the magnetics and gravity data independently. The framework is also used to jointly invert these data and produce a model that is consistent with both data sets. Cumulative volume estimates with depth are informative and can help decide whether in situ or ex situ sequestration might be appropriate. Using each of the inverted models, we estimate cumulative volume of reactive rock as a function of depth.
Introduction
Even if we are successful in reaching net-zero carbon emissions by 2050, the recent Intergovernmental Panel on Climate Change (IPCC) report predicts that global warming will exceed 1.5°C in the 21st century (IPCC 2021). This conclusion highlights that while it is important to transition off of hydrocarbon-based energy sources, it is also critical that we remove CO2 from the atmosphere. Since 1850, over 2000 Gt of anthropogenic CO2 have been emitted (International Energy Agency 2021), and despite the COVID-19 pandemic, energy-related emissions in 2020 were 31.5 GtCO2. Thus, we are faced with the need for negative emissions technologies (NETs), also referred to as carbon dioxide removal (CDR) methods, to remove hundreds to thousands of Gt of CO2 from the atmosphere (IPCC 2021; National Academies of Sciences, Engineering, and Medicine 2019). The 2019 report from the National Academies (National Academies of Sciences, Engineering, and Medicine 2019) provides an overview of 6 categories of approaches to capture and store CO2. Of all of the approaches presented, geologic storage is the most desirable because of both the long time scales of storage as well as the global capacity. Two mechanisms involve geologic storage: sequestration of CO2 in sedimentary formations and carbon mineralization of CO2.
The first mechanism is the more mature of the two and involves the injection of supercritical CO2 into depleted hydrocarbon reservoirs or deep saline aquifers (National Academies of Sciences, Engineering, and Medicine 2019; Kelemen et al. 2019). This approach relies on there being a porous, permeable reservoir overlain by a caprock that prevents CO2 migration out of the reservoir. Major risks for CO2
migration include non-sealed fractures in the caprock and leaky wellbores. Geophysical data, including seismic, InSAR, and electromagnetics can be employed to monitor wellbore integrity, injection, and leakage; there are multiple field examples worldwide (e.g. (Vasco et al. 2010; Ajo-Franklin et al. 2013; Heagy and Oldenburg 2019; Wilt et al. 2020)).
The second mechanism for geologic storage is carbon mineralization; this is the focus of this paper. This approach utilizes magnesium and calcium-bearing minerals in mafic or ultramafic rocks that react with CO2and convert it to solid, stable, carbonate minerals (National Academies of Sciences, Engineering, and Medicine 2019; Seifritz 1990; Lackner et al. 1995). Although the science and technology of carbon mineralization are less mature than geologic storage in sedimentary units, this approach presents several potential advantages. Storage is permanent with minimal risk of leakage (Zhang and DePaolo 2017), and it is estimated that there may be orders of magnitude more global capacity than sequestration in sedimentary units (Kelemen et al. 2019). Broadly speaking, there are two strategies for promoting reaction with CO2: (a) in-situ approaches which require circulating CO2-bearing fluids through subsurface formations where they can react, and (b) ex-situ (or surficial) approaches use minerals that have been brought to the surface, for example in mine-tailings (National Academies of Sciences, Engineering, and Medicine 2019; Kelemen et al. 2019). In some places, reactive rocks are co-located with minerals of economic interest. At these locations, it may be possible for mines to have net-zero, or even net negative, emissions and “carbon mines” that would bring reactive minerals to the surface might become economically viable. For in-situ approaches, there are only two pilot projects that have been conducted: the Carbfix Project in Iceland (Gíslason et al. 2018) and the Wallula Project in Washington State (McGrail et al. 2014). Both of these sites are injections into basalts. To the best of our knowledge, there have not yet been in-situ experiments in ultramafic rocks.
For either an in-situ or ex-situ approach, the first step is to identify the location and volume of reactive rocks. Recent work in (Cutts et al. 2021; Mitchinson et al. 2020) examines the physical properties of ultramafic rocks as they are altered. Of particular interest are rocks that have been altered through serpentinization and can react with CO2 through a carbonation reaction. As ultramafic rocks go through each reaction phase, their magnetic susceptibility and density are altered. This motivates the use of geophysical data to locate these rocks, estimate their depth and volume, and potentially even estimate their reactivity.
In this paper, we will identify some of the opportunities and areas of future research for the use of geophysical data for carbon mineralization. We will focus on ultramafic rocks using the physical property work (Cutts et al. 2021) and the initial inversion of magnetic data over ultramafic sites in BC (Mitchinson et al. 2020) as motivation. The paper is organized as follows. First, we will discuss how the physical properties of ultramafic rocks change as they are altered. We then introduce a representative model and simulate gravity and magnetic data. Next, we will introduce approaches for independently and jointly inverting these data. Finally, we will conclude with a discussion of some of the practical challenges and research opportunities for the use of geophysical data for carbon mineralization. All of the code used to generate figures for this paper is available at https://github.com/simpeg-research/Fast-Times-2021-Carbon-Mineralization.
Physical Properties
In the context of ultramafic rocks, rocks that have been serpentinized are of particular interest for carbon mineralization because they are particularly reactive with CO2. Following the discussion in Cutts et al. 2021, the general formula for a serpentinization reaction can be summarized by:
Rocks that have been serpentinized are then susceptible to reaction with CO2 through carbonation reactions. Brucite is of particular interest for carbon mineralization because it reacts rapidly with CO2 at surface conditions (Kelemen et al. 2019). The general formulae for carbonation reactions are (Hansen et al. 2005):
The key connection point between any question about the subsurface and the use of geophysical data is physical properties. Serpentinization and carbonation reactions change both the density and magnetic susceptibility of rocks. The density and magnetic susceptibility of fresh ultramafic rocks should be similar to that of olivine and pyroxene (3.10-3.30 g/cc and 10-3 SI). Minerals produced through serpentinization (serpentine, brucite and magnetite) are typically less dense than fresh ultramafics. Brucite and serpentine have densities of approximately 2.57 and 2.39 g/cc, respectively. Magnetite also has a larger magnetic susceptibility than unaltered phases. So overall, serpentinized rocks have a lower density and higher magnetic susceptibility than fresh ultramafics. Brucite, serpentine and magnetite are consumed in a carbonation reaction forming magnesite as well as talc and quartz. Magnesite has a density of approximately 3.00 g/cc, and thus carbonated rocks are more dense and have a lower magnetic susceptibility than serpentinites (Cutts et al. 2021; Hansen et al. 2005).
To quantify how physical properties of ultramafic rocks change as a function of alteration, Cutts et al. 2021 made geochemical and petrophysical measurements of over 400 samples of ultramafic rocks from the Canadian Cordillera. In Figure 1, we reproduce portions of Figure 6 in their paper and show physical property measurements as a function of Loss On Ignition (LOI), which is a proxy for alteration. From an LOI of 0% to ~13%, fresh ultramafic rocks undergo a serpentinization reaction where fresh ultramafics are hydrated with the addition of water. Beyond an LOI of ~13%, magnesium and calcium-rich minerals react with CO2in carbonation reactions. Rocks that have been serpentinized but not yet carbonated have the most CO2 uptake capacity and therefore are of the most interest for carbon mineralization. Figure 1 shows that (a) density and (b) magnetic susceptibility vary as a function of LOI. The trend with density is very clear: as fresh ultramafic rocks are serpentinized, the density decreases from 3.2 g/cc to 2.7 g/cc. Magnetic susceptibility is much more variable, but broadly speaking, serpentinized rocks are more susceptible than fresh ultramafic or carbonated rocks. Remanent magnetization is also important to consider in these settings. Though also quite variable, serpentinized rocks can have a substantial remanent magnetization, with many samples in the data collected by Cutts et al. 2021 having Koenigsberger ratios larger than 1 (see figures 6 and 8 in Cutts et al. 2021). In the examples we develop in this paper, we will not include remanence. In practice, it will be important to consider Mitchinson et al. 2020 and inversion approaches that enable us to recover the magnetization vector such as developed in Fournier, Heagy, and Oldenburg 2020.
Figure 1. Physical property measurements of (a) density and (b) magnetic susceptibility of ultramafic rocks as a function of Loss On Ignition from (Cutts et al. 2021). From an LOI of 0% to ~13%, rocks are undergoing serpentinization. Beyond an LOI of 13%, rock undergo carbonation. (c) Cross plot of density and magnetic susceptibility colored by Loss On Ignition (LOI).
Figure 1 (c) shows the same physical property information as a cross-plot of density and magnetic susceptibility with the points colored by LOI. The orange yellow points are rocks that have been serpentinized but not yet significantly carbonated. Cross-plots are generally insightful and often yield well-defined clusters in physical property space that are associated with particular rock types. Although this is not easily visible here, there is a concentration of samples that have a high magnetic susceptibility, low density, and within the desired LOI range for reactive serpentinized rocks.
Motivating Example
For both in-situ or ex-situ approaches to carbon mineralization, the goal of using geophysical data is to locate and delineate potentially reactive rocks. Important factors to characterize include: (a) the depth, which determines whether an in-situ or ex-situ approach should be used, (b) the volume, which influences the CO2 sequestration capacity, and (c) for ultramafic rocks, the degree of serpentinization or carbonation which controls how reactive the rocks are.
To illustrate the role of geophysical data and avenues for future development, we introduce a simple representative model in Figure 2a. This model is motivated by the geologic setting at the Decar site (e.g. see the maps and data in (Mitchinson et al. 2020)). We summarize the physical properties of each unit in Table 1. A serpentinized block with a magnetic susceptibility of 0.15 SI and a density of 2.7 g/cc is embedded in a halfspace background (0 SI, 2.9 g/cc: the background is composed of mafic volcanic and sedimentary rocks). Within that block, there is a region that has been carbonate-altered, it is wider in the north (2km) and tapers to the south (1km). The carbonated region has a magnetic susceptibility of 0.05 SI and a density of 3.0 g/cc.
Table 1. Physical properties of the rock units shown in Figure 2.
Figure 2. (a) Representative model of a serpentinized rock unit that contains a center region that has been carbonated. The physical properties assigned to each unit are shown in Table 1. The total volume of serpentinized rock is 35 km3 and carbonated rock is 15 km3. Panels show the simulated (b) magnetic and (c) gravity data.
We simulate the gravity anomaly and magnetic response (assuming a vertical inducing field) and show those simulated data in Figure 2 b and c, respectively. The data spacing in x and y is 250m. Gaussian random noise of 0.05 mGal and 1 nT are added to the gravity and magnetic data, respectively. These simulations, and the inversions we show subsequently, were run using the open-source SimPEG software (Cockett et al. 2015).
Individual inversions
As a first step, we invert each data set independently. We use a deterministic approach to solve the inverse problem. We formulate the inversion as an optimization problem where we minimize an objective function that is comprised of a data misfit term and a model norm term (Tikhonov and Arsenin 1977).
The data misfit term is a measure of how well the data that are predicted given a model m reproduce the observed data. The model norm (or regularization) term is included to overcome the non-uniqueness of the inverse problem and can be used to incorporate geologic knowledge or assumptions (see (Oldenburg and Li 2005) for a tutorial). A standard approach is to use a term of the form:
The first term is often referred to as the “smallness” term while the subsequent terms are referred to as first-order “smoothness” in the x, y, z directions, respectively. Generally p=q=2 is implemented but (p,q)<2 can be used to promote sparse or compact structures in the inversion (Fournier and Oldenburg 2019). Throughout, we will use l2 to indicate an inversion where p=q=2 , and lpq where p<2 and/or q<2.
Magnetics
In Figure 3, we show inversion results for the magnetics data using (a) a standard l2 inversion, and (c) a l01 inversion. To compensate for the inherent lack of depth resolution of potential field data, we apply depth weighting in both inversions (Li and Oldenburg 1996). An l2 inversion is the standard approach for inverting most data types and, in general, will favour models with smooth structures and variations in physical properties. There are some artefacts in the recovered model: there is some moderately susceptible material at depth and in the padding cells. We can also see the impacts of this regularization in the histogram of recovered susceptibility values in Figure 3 (c) where there is a smooth distribution of physical property values with the maximum at 0 SI (the background susceptibility. We also see that very few cells reach the true susceptibility of the serpentinized unit (0.15 SI).
To promote sharper contrasts and more compact structures, we can alter the norms. In an lpq inversion, l2 norms are used during the first stage of the optimization to produce a smooth initial model. Once the data are fit, the norms are then updated to allow for more compact structures and discontinuities in the final recovered model. Using p=0 and q=1, we recover the results shown in Figure 3c. In the histogram in Figure 3d, we see that there are many cells with the true background value, but this distribution is not “smeared’’ to other low susceptibility values. Additionally, we see that the l01 inversion recovers some moderate susceptibility values in the center, carbonated region, which in the true model has a susceptibility value of 0.05 SI.
Figure 3. Models recovered using (a) an l2 inversion of the magnetic data and (c) an l01 inversion of the magnetic data. The top plot in each is a depth slice at z=-950m and the bottom panel shows a cross section along the line y=0m. Plots (b) and (d) show a histogram of the recovered susceptibility values. The histograms are normalized so that the area integrates to 1. Note that the y-limits clip the true maximum value in each.
An important control on how much CO2 can be sequestered is the total volume of serpentinized rock, and thus it is of interest to use these results to estimate the volume of serpentinized rock. In our synthetic model, the serpentinized rocks are distinguished from the carbonated rocks and background by having a high susceptibility (0.15 SI). Ideally, we would like to invert data and delineate the region of the model having this value. However, inversions generally produce models which are smoothed in some sense and hence very few cells in either the l2 or l01 inversions reach the true value of 0.15 SI. In the l01 inversion, the maximum recovered susceptibility is 0.14 SI, so using the true value of 0.15 SI would result in a null volume. There is no clear way to choose a cut-off when looking at the physical property histograms in Figure 3. To examine how our choice of threshold impacts our estimate of volume of serpentinized rock, we plot estimated volume as a function of threshold for both inversion results in Figure 4. For each chosen threshold value, we compute the volume of material with a susceptibility larger than or equal to that value. We confine the region where we are developing an estimate to the x-y extent of the survey, and within a maximum depth of 3.5km.
Figure 4. Estimated volume of serpentinized rock for a range of susceptibility thresholds.
In Figure 4, we see that the threshold that gives us the correct volume for both inversions is 0.08 SI. However, without knowing the true volume, it is difficult to estimate a-priori what the correct threshold should be. In order to estimate this quantity we set up a proxy experiment in which the geology is a serpentinized block in a half-space. We choose a 4km by 4km block that extends from 500m to 1500m depth, and is positioned in the center of the survey area. We assign the same magnetic susceptibility as the true serpentinized unit, 0.15 SI. We use the same survey geometry as before and perform a forward simulation and add noise. We then perform l2 and l01 inversions of these data using identical model regularizations and parameter choices as in our initial inversions. The results are shown in Figure 5. We then ask the question: what threshold value gives us the correct estimated volume for the block? For both inversions, the correct threshold turns out to be ~0.07 SI; this value is in good agreement with the value of 0.08 SI for the three block model in Figure 4. Using a threshold of 0.07 SI for the three block model gives us estimates of 40 km3 and 43 km3 for the l2 and l01 inversions, respectively. The threshold values obtained by doing this proxy experiment can be expected to vary with the geometry of the model. We tested this by using a 6km by 6km block and it yielded the same result of 0.07 SI. The procedure has some degree of robustness but choosing a different geometry for the proxy block, changing its dimension, altering the depth, or assigning different values for the assumed physical properties would all have altered the outcome.
Figure 5. (a) Simulated magnetic data for a 4km by 4km block, (b) recovered l2 susceptibility model, (c) recovered l01 susceptibility model, (d) estimated volume of serpentinized rock as a function of threshold choice.
Gravity
Next, we perform inversions of the gravity data. Similar to the magnetics, we perform two inversions, first, an l2 inversion, shown in Figure 6(a) as well as an l01 inversion, shown in Figure 6(b). Visually, the gravity inversions show a better distinction between the 3 units than the magnetics. Relative to the background, the serpentinized units are less dense, whereas carbonated units are more dense. This is a more diagnostic contrast than what we observed with magnetics. The positive side-lobes we see in the l2 inversion result are common when working with blocky models (Oldenburg and Li 2005). The l01 produces more compact structures, as expected, and reduces the side-lobe effects. It also improves the estimate of the depth of the bottom of the serpentinized unit. Both inversion results show structure coming all of the way to the surface.
Figure 6. Models recovered using (a) an l2 inversion of the gravity data and (c) an l01 inversion of the gravity data. The top plot in each is a depth slice at z=-950m and the bottom panel shows a cross section along the line y=0m. Plots (b) and (d) show a histogram of the recovered density values.
In Figure 7, we plot the estimated volume of serpentinized rock over a range of density cut-offs. Instead of considering all material above a given threshold as we did in magnetics, we consider all material below a given threshold because the serpentinized rock is less dense than the background. As compared to the magnetics, we can see a greater distinction between the l2 and l01 inversion approaches. The l2 has a steeper slope as compared to the l01 inversion. This indicates that a threshold chosen from the l01 inversion has a less substantial impact on the estimated volume than the l2 inversion.
Figure 7. Estimated volume of serpentinized rock for a range of density thresholds.
To estimate the choice of an appropriate threshold, we again perform inversions for the simple 1-block models where the blocks have a density equal to that of the serpentinized unit: 2.7 g/cc. The results are shown in Figure 8. Similar to the 3 block model, we again see that the volume estimates from the l2 inversions have a much steeper slope than the l01. In this case, the threshold values for the two different regularizations are different: 2.83 g/cc for the l2 inversion and 2.79 g/cc for the l01. This agrees well with the best threshold choices for the three-block model of 2.84 g/cc and 2.81 g/cc for the l2 and l01 results in Figure 7. Applying thresholds of 2.83 g/cc and 2.79 g/cc to their respective inversions gives us the same volume estimate of 27 km3 for both. The good agreement with choice of threshold and resultant volume estimate supports our hypothesis that our proxy strategy is working so long as we keep all of the parameters in the inversion the same.
Figure 8. (a) Simulated density data for a 4km by 4km block, (b) recovered l2 density model, (c) recovered l01 density model, (d) estimated volume of serpentinized rock as a function of threshold choice.
Summary: Individual Inversions
Standard l2 inversions produce models with smooth structures. For both the gravity and magnetic inversions, we can see this tends to smear out structures in depth. The resultant distribution of physical properties also tends to be smoothly varying with a concentration of values near the background value. Replacing the l2 model norm with a l01 norm promotes more compact structures. In this example, it results in more easily interpretable structures. To estimate the appropriate physical property threshold, we set up a proxy experiment and inverted synthetic data over a simple one-block model. In Table 2, we summarize the estimated thresholds obtained from our proxy-experiment and the resultant estimates of the volume of serpentinized rock from each of the inverted models. The good agreement between both the thresholds and resultant volumes is encouraging for the use of a proxy model for estimating an appropriate physical property threshold to delineate a rock unit. In an analysis of field data, one might run several examples, varying the size and physical properties of a representative target to obtain an estimate and gauge the uncertainty of that estimate.
Table 2. Summary of the thresholds and volume estimates obtained using a proxy experiment of a one-block model. Note that the true volume of serpentinized rock is 35 km3.
In addition to estimating volumes, the structure and depths of each unit is of practical importance. If serpentinized rock is near the surface, then it may be a candidate for ex-situ carbon mineralization, whereas if it is deeper it is unlikely that it would be economic to mine, but may be a candidate for in-situ carbon mineralization. Mitchinson et al. 2020 uses a depth cutoff of 1km for rocks that may be feasible for ex-situ operations, and within 2km of the surface for in-situ operations. The magnetics inversion shows a good estimate to the depth of the top of the unit, but the gravity brings structure to the surface. So we are left with a question of how to reconcile these results. In the next section, we will examine the use of a Petrophysically and Geologically Guided (PGI) framework for jointly inverting these data.
Including petrophysical data in the inversion
The Petrophysically and Geologically guided Inversions (PGI) framework
The inverse problem is non-unique. Including knowledge of physical property and geologic information in an inversion can improve our ability to recover the structure of interest. Prior information is also a key component to design coupling terms that link together various physical property models, which enables joint inversion of multiple geophysical surveys. The Petrophysically and Geologically guided Inversion (PGI) framework gives us an avenue for including petrophysical information and jointly inverting multiple geophysical data sets.
Figure 9. The PGI framework (modified from Astic and Oldenburg (2019)) is composed of three interlocked inverse problems (diamond nodes) over the geophysical, petrophysical, and geological information (represented by the rectangular nodes).
The PGI framework, introduced in Astic and Oldenburg, 2019 and Astic, Heagy, and Oldenburg, 2020, formulates the inverse problem from a probabilistic perspective. Petrophysical and geological information are represented as probability distributions, which are then encoded in the regularization of the Tikhonov objective function (Equation [1]). To fit all the information at once, the framework cyclically solves three interlocking inverse problems over the geophysical, petrophysical, and geological data (Figure 9). After each geophysical inversion iteration (Figure 9, process (1)), the PGI framework learns a petrophysical representation of the subsurface from the petrophysical data (laboratory measurements, borehole logs, etc.) and current geophysical model (Figure 9, process (2)). The petrophysical model can either be held fixed if we have a strong confidence in the prior petrophysical data or it can be iteratively updated to benefit from the knowledge gained at each iteration from the geophysical data. Then PGI builds a quasi-geology model (Li et al. 2019), a representation of the geology based on the geophysics, by classifying the current geophysical model into rock units using the petrophysical model and geological information. This quasi-geology model serves to update the regularization weights and reference model for the next geophysical inversion iteration; the cycle continues until all geophysical, petrophysical, and geological data are fit. More details about these different processes are given in Astic and Oldenburg 2019.
Representation of petrophysical information: Gaussian Mixture Model (GMM)
Within the PGI framework, petrophysical information is represented as a probability distribution modelled by a Gaussian Mixture Model (GMM) (example in Figure 10). A GMM is composed of a weighted sum of normal distributions, also called Gaussians. Each Gaussian distribution represents a rock unit. The goal of a PGI is to find a physical properties model of the subsurface that both fits the geophysical data and reproduces the GMM distribution.
In the present case, we have three rock units: the background, the carbonated, and the serpentinized facies. Each formation is characterized by: (1) the mean of its physical properties (for example 3.0 g/cc and 0.05 SI for the carbonated; represented as triangles in Figure 10), (2) its covariance which quantifies how spread its petrophysical values are around the mean, and (3) the probability of each rock prior to geophysical observation. This last parameter weights each normal distribution and can encode local geological information, or just be set at a constant global value everywhere if no such information is available. The chosen global value has minimum impact on the inversion (Astic 2020). In Figure 10, the central panel presents the joint density and magnetic susceptibility contrasts probability distribution; this is what is used to perform multi-physics inversion with PGI, using the petrophysical distribution as a coupling term. The side panels represent the “marginal probability distribution” for each physical property individually; they will be used for the individual magnetic and gravity inversions.
Figure 10. Physical properties information represented as a GMM. The left and bottom panels show the magnetic susceptibility and density probability distribution, respectively. The middle panel represents their joint probability distribution.
Individual PGI: Magnetics
We start by performing a PGI inversion of each data set independently. In addition to requiring that the model fit the geophysical data through the use of a geophysical data misfit, we now also include a petrophysical misfit term. With the geophysical data misfit, we compare the data predicted from a forward simulation over the geophysical model with the true, observed data. With the petrophysical misfit term, we compare the physical property estimated from the geophysical model with the mean physical property of the rock unit it is classified as in the geologic identification. For example, a cell with a value of 0.13 SI would be classified as serpentinized and therefore compared with the mean value of the serpentinized facies (0.15 SI), whereas a cell with a value of 0.06 SI would be classified as carbonated and compared with that value (0.05 SI). To construct the misfit metric, standard deviations of the petrophysical data are also taken into account; these details are discussed in (Astic and Oldenburg 2019). PGI is an iterative procedure, so with each update to the physical property model in the geophysical inversion component, we re-evaluate the petrophysical misfit and update the associated quasi-geology model. In doing so, we obtain a model that respects both the geophysical and petrophysical data. Further, the quasi-geology model is constructed as a part of the inversion, meaning we do not need to select a threshold to identify serpentinized rock and estimate its volume, as we did in the previous inversions. Here, we can use the geologic classification produced by PGI and directly estimate the volume from the quasi-geology model.
In Figure 11, we show the outcome of PGI applied to the magnetic data. In this case, we see that, we successfully identify two blocks of higher susceptibility whose position and petrophysical signature correspond well to the serpentinized unit in Figure 11a). Contrary to the l2 and l01 inversions, the histogram in (c) shows distinct clusters that correspond to the petrophysical characteristics of each rock unit (background in white, serpentinized in green, and carbonated in blue). The solid line in (c) shows the GMM we chose for this inversion. The background has a mean susceptibility of 0 SI, the carbonated unit 0.05 SI, and the serpentinized unit 0.15 SI. This is what we compare the physical property values from the geophysical inversion to in the petrophysical misfit. A “perfect” inversion would have a histogram that aligns with the defined GMM. In this inversion, the distributions of the recovered susceptibilities for the carbonated and serpentinized unit are shifted to the left, having lower mean values. There are also more values attributed to the carbonated unit than the true distribution. This effect, as well as the “halo” structure we observe in the quasi-geology model (b), where the regions classified as serpentinized rock are surrounded by carbonated rock magnetic susceptibility, can occur in single physical property inversions. When two units have contrasts that are both positive or both negative compared to the starting value of the background unit, the inversion tries to explain the data first with smaller contrast (carbonated rock) and only puts larger values where it is strongly indicated in the data (see Astic (2020) for more information on “consecutive clusters”). As a result, we underestimate the volume of serpentinized rock, at 21 km3. With PGI, we are also able to estimate the volume of carbonated rock (116 km3). This would be challenging to do with the l2 and l01 inversions as it would require that we choose both upper and lower bounds to attribute to the carbonated rock.
Figure 11. Results from inverting the magnetics data with the PGI framework. (a) Magnetic susceptibility model, (b) associated quasi-geology model, and (c) histogram of the recovered magnetic susceptibility values.
Individual PGI: Gravity
Next, we perform PGI using the gravity data and show the results in Figure 12. We see that 3 distinct volumes are recovered, two with the petrophysical signature of the serpentinized unit (mean density of 2.7 g/cc), and one in the middle with the petrophysical signature of the carbonated unit (mean density of 3.0 g/cc). Since the physical property contrasts are opposite in sign (the density of the serpentinized unit is less than the background, whereas the density of the carbonated unit is larger), this is a simpler distribution to work with in the PGI framework than we encountered with magnetics. In the histogram in Figure 12 (c), we see that the mean contrasts are slightly underestimated, but adequately fit, and there is a clear separation between the physical properties of each unit. There are also improvements in the recovered model. The depth to the top of the unit identified to be serpentinized is well-recovered, in contrast to what we observed in the l2 and l01 inversions (Figure 6). The total volume of serpentinized rock is estimated from the quasi-geologic model, 53 km3 for the serpentinized unit and 13 km3 for the carbonated unit. The total volume of carbonated rock is quite close to the true volume (15 km3), and the volume of serpentinized rock is overestimated, which may be expected because the recovered densities of the serpentinized rock have a lower contrast than the true value.
Figure 12. Results from inverting the gravity data with the PGI framework. (a) Density model, (b) associated quasi-geology model, and (c) histogram of the recovered density values.
There are several benefits to the PGI framework for this application: (a) prior petrophysical knowledge can be included in the inversion (with uncertainties); (b) rather than a smooth range of physical property values recovered in the l2, l01 inversions, the recovered petrophysical signatures are more representative of the true petrophysical signature; and (c) the quasi-geology model, built by reproducing the petrophysical signature, gives a volume estimate that is independent of the choice of an ad hoc threshold.
Despite the listed benefits in using a PGI approach, we note that the quasi-geology models obtained from magnetics and gravity are different. For example, in the magnetics inversion, the serpentinized unit comes to the surface whereas in the gravity inversion, the depth to the top of this unit is 300m. Also the volume estimates of the serpentinized body are different and neither is in agreement with the true model. The differences between the two quasi-geology models can be reconciled by taking advantage of the complementary information in gravity and magnetic data. The next step is to use PGI to jointly invert these data sets.
Joint PGI
Now, we seek to find a model that fits both geophysical data sets and is consistent with our knowledge of the petrophysics. There are many avenues for approaching a joint inversion, and (Haber and Holtzman Gazit 2013; Moorkamp et al. 2016) provide an overview of some approaches. The two major categories use either structural or petrophysical characteristics to couple data that are sensitive to different physical properties. The PGI approach uses petrophysical information to couple multiple data sets. Now rather than the marginal distributions, we will work with the 2D GMM as shown in Figure 10. For the joint PGI, we are asking the inversion now to fit: a misfit associated with the magnetics data, a misfit associated with the gravity data, and a petrophysical misfit that is defined in 2D space. The petrophysical misfit is what couples the two data sets, and in looking at Figure 10, what is considered for that misfit is now distances in 2D space from the cluster centers. All of these misfit terms must be balanced appropriately, and we provide details on our approach in the Astic, Heagy, and Oldenburg 2020 paper.
Using the PGI framework, we perform a joint inversion with the gravity and magnetic surveys. The result of that joint inversion can be seen in Figure 13. First, examining the magnetic susceptibility model, we see that the depth to the highly susceptible serpentinized unit is in much better agreement with the true model, and we are no longer plagued by the “halo” effects we saw in Figure 11. Looking at the density model first (Figure 13a), we see that the contacts between the various units are much better recovered thanks to the addition of the magnetic information. Both models share the same underlying quasi-geology model (Figure 13c). The petrophysical information is fit (Figure 13d, where we show the coupling GMM and the crossplots of the physical properties models), it appears that the density signatures are better reproduced than in the individual PGI inversion. There is more spread in the distribution of magnetic susceptibility than the density values, but this is reflective of what is expected (e.g. see the GMM in Figure 13c), and by using multiple physical properties in the inversion, we can better handle this ambiguity. Using the joint PGI approach, we produce a single quasi-geology model that is consistent with the magnetics and gravity data, as well as our prior knowledge of the petrophysics. The total volume estimates of the serpentinized and carbonated rock units are 58 km3 and 21 km3, respectively.
Figure 13. Results from jointly inverting the magnetics and gravity data with the PGI framework.
Cumulative volume estimates
It is satisfying to obtain a result that is consistent between the magnetics and gravity data and overall does a good job delineating the three rock units as we did with the joint PGI inversion (Figure 13). The total estimated volume of serpentinized rock is overestimated, which is due to the extra material at depth, beneath the unit. Both gravity and magnetic methods have limited depth information. Sensitivity decreases with depth and at some point our recovered models reflect a-priori information and the nature of our algorithm. Thus, our estimates of total volume can be impacted by this material at depth which may not be substantially supported in the data. As mentioned previously, the total volume is of interest, but perhaps more important is the volume as a function of depth as this controls whether it is practical to perform ex-situ or in-situ carbon mineralization operations.
To compare the various approaches, we plot the cumulative volume of serpentinized rock estimated from each of the inversions as a function of depth in Figure 14. Depth extent is on the x-axis, and the volume of rock above a given depth that we attribute to serpentinized rock from each inversion is on the y-axis. For the l2 and l01 inversions, we use the threshold values from the 4km block as defined in Table 2. For the volume estimates as a function of depth from magnetics, we see that overall, there is good agreement between the inversion results, with the exception of the individual PGI. The poor estimation with depth from the PGI inversion of the magnetics data alone is consistent with the cross section observed in Figure 11. More of the susceptible material was attributed to the carbonated unit than the serpentinized unit. This highlights the value of running multiple styles of inversions. Where there is consistency between results, we can be more confident in the resultant model. Turning our attention to the gravity, we see that the l2 inversion puts material right at the surface and both the l2 and the l01 underestimate the volume as a function of depth at 1000m depth. The individual gravity PGI is an improvement. It slightly underestimates the volume between 500 and 1300m depth (the true depth of the bottom of the unit). The joint PGI models result in the best agreement between the surface and 1300 m depth. It also has the advantage that it is consistent for both susceptibility and density, because we require that the inversion produces a result that agrees with both data sets.
Figure 14. Estimated volume of serpentinized rock above a given depth for the (a) susceptibility and (b) density models recovered in each of the inversions.
Summary
Carbon mineralization could be an important mechanism for reducing atmospheric CO2 as we work to transition to a carbon neutral economy. When rocks with CO2 sequestration potential have physical properties that are distinct, geophysical data can be used to locate and delineate rocks with sequestration potential. Inversion is an essential tool for obtaining models of the subsurface that we can interpret geologically. The inverse problem is non-unique, and therefore assumptions and prior information need to be included. Defining a model norm is a standard approach. Typically a l2 norm is used, and this tends to produce models with smooth features. More compact structures can be promoted by using norms less than 2. Using l01 produced more compact structures. This can simplify the interpretation of the resultant model because boundaries between units are more easily defined, and some inversion artefacts were reduced. When estimating a volume of rock, which is an important factor for the amount of CO2 that can be sequestered, the use of an l01 norm produced estimates that were less sensitive to the exact choice of threshold than an l2 inversion.
Using the Petrophysically and Geologically Guided Inversion (PGI) framework, we illustrated how knowledge of physical properties can be brought into the inversion. There are several benefits that this inversion approach offers: (1) physical property information, with associated uncertainties, is fit as a part of the inversion, (2) a quasi-geology model is produced, and (3) joint inversions that use multiple data sets that are sensitive to different physical properties can be performed. By using the joint PGI framework, we produced a quasi-geology model that was consistent with both the gravity and magnetics data sets. The model that was produced was also the best for estimating the volume of rocks with sequestration potential as a function of depth.
The inversions employed, the l2, l01 and the PGI, increase in complexity and the user-knowledge that is required to obtain reasonable results. The more complex inversions can offer some benefits, as we have described. This is not to say that they should be the first inversion employed. In general, starting with the simpler l2 inversion is advisable for a first interpretation. It is much simpler to run l2 inversions to examine the impacts of choices of inversion parameters such as the noise characteristics and the use of depth or sensitivity weighting. Once these have been tuned on a simple inversion, they can be used for the more complex inversions. It is also advantageous to have multiple models of the subsurface to examine which features are consistent or not between them.
The model we chose to work with was seemingly quite simple, but shows that inversion and obtaining high quality models of the subsurface is an iterative process. We did not include the impacts of remnant magnetization in our analysis. This will be important to address and presents an opportunity for future research.
Acknowledgements
The authors are grateful to Greg Dipple, Peter Bradshaw, and Randy Enkin for multiple conversations about the geochemistry and important questions for carbon mineralization. We are grateful to Dianne Mitchinson for discussions as we designed a representative model and Jamie Cutts for making the physical property data available. Finally, we are grateful to the SimPEG team, and Dominique Fournier in particular for making his work on compact and sparse norms available through the project.
References
Ajo-Franklin, J. B., J. Peterson, J. Doetsch, and T. M. Daley. 2013. “High-resolution characterization of a CO2 plume using crosswell seismic tomography: Cranfield, MS, USA.” International Journal of Greenhouse Gas Control 18: 497–509. https://doi.org/10.1016/j.ijggc.2012.12.018.
Astic, Thibaut. 2020. “A framework for joint petrophysically and geologically guided geophysical inversion.” PhD thesis, University of British.
Astic, Thibaut, Lindsey J. Heagy, and Douglas W. Oldenburg. 2020. “Petrophysically and geologically guided multi-physics inversion using a dynamic Gaussian mixture model.” Geophysical Journal International 224 (1): 40–68. https://doi.org/10.1093/gji/ggaa378.
Astic, Thibaut, and Douglas W. Oldenburg. 2019. “A framework for petrophysically and geologically guided geophysical inversion using a dynamic Gaussian mixture model prior.” Geophysical Journal International 219 (3): 1989–2012. https://doi.org/10.1093/gji/ggz389.
Cockett, Rowan, Seogi Kang, Lindsey J. Heagy, Adam Pidlisecky, and Douglas W. Oldenburg. 2015. “SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications.” Computers & Geosciences 85 (December): 142–54. https://doi.org/10.1016/j.cageo.2015.09.015.
Cutts, J. A., K. Steinthorsdottir, C. Turvey, G. M. Dipple, R. J. Enkin, and S. M. Peacock. 2021. “Deducing Mineralogy of Serpentinized and Carbonated Ultramafic Rocks Using Physical Properties With Implications for Carbon Sequestration and Subduction Zone Dynamics.” Geochemistry, Geophysics, Geosystems 22 (9). https://doi.org/10.1029/2021GC009989.
Fournier, Dominique, Lindsey J. Heagy, and Douglas W. Oldenburg. 2020. “Sparse magnetic vector inversion in spherical coordinates.” GEOPHYSICS 85 (3): J33–49. https://doi.org/10.1190/geo2019-0244.1.
Fournier, Dominique, and Douglas W. Oldenburg. 2019. “Inversion using spatially variable mixed Lp norms.” Geophysical Journal International 218 (1): 268–82. https://doi.org/10.1093/gji/ggz156.
Gíslason, Sigurdur R., Hólmfrídur Sigurdardóttir, Edda Sif Aradóttir, and Eric H. Oelkers. 2018. “A brief history of CarbFix: Challenges and victories of the project’s pilot phase.” Energy Procedia 146: 103–14. https://doi.org/10.1016/j.egypro.2018.07.014.
Haber, Eldad, and Michal Holtzman Gazit. 2013. “Model Fusion and Joint Inversion.” Surveys in Geophysics 34 (5): 675–95. https://doi.org/10.1007/s10712-013-9232-4.
Hansen, Lyle D., Gregory M. Dipple, Terry M. Gordon, and Dawn A. Kellett. 2005. “Carbonated serpentinite (listwanite) at Atlin, British Columbia: A geological analogue to carbon dioxide sequestration.” Canadian Mineralogist 43 (1): 225–39. https://doi.org/10.2113/gscanmin.43.1.225.
Heagy, Lindsey J., and Douglas W. Oldenburg. 2019. “Direct current resistivity with steel-cased wells.” Geophysical Journal International 219 (1): 1–26. https://doi.org/10.1093/gji/ggz281.
International Energy Agency. 2021. “Global Energy Review 2021.” Global Energy Review 2020, 1–36. https://iea.blob.core.windows.net/assets/d0031107-401d-4a2f-a48b-9eed19457335/GlobalEnergyReview2021.pdf.
IPCC. 2021. “Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [Masson-Delmotte, V., P. Zhai, A. Pirani, S. L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen,” Cambridge University Press, no. In Press: In Press.
Kelemen, Peter, Sally M. Benson, Hélène Pilorgé, Peter Psarras, and Jennifer Wilcox. 2019. “An Overview of the Status and Challenges of CO2 Storage in Minerals and Geological Formations.” Frontiers in Climate 1 (November): 1–20. https://doi.org/10.3389/fclim.2019.00009.
Lackner, Klaus S., Christopher H. Wendt, Darryl P. Butt, Edward L. Joyce, and David H. Sharp. 1995. “Carbon dioxide disposal in carbonate minerals.” Energy 20 (11): 1153–70. https://doi.org/10.1016/0360-5442(95)00071-N.
Li, Yaoguo, Aline Melo, Cericia Martinez, and Jiajia Sun. 2019. “Geology differentiation: A new frontier in quantitative geophysical interpretation in mineral exploration.” Leading Edge 38 (1): 60–66. https://doi.org/10.1190/tle38010060.1.
Li, Yaoguo, and Douglas W. Oldenburg. 1996. “3-D inversion of magnetic data.” GEOPHYSICS 61 (2): 394–408. https://doi.org/10.1190/1.1443968.
McGrail, B. P., F. A. Spane, J. E. Amonette, C. R. Thompson, and C. F. Brown. 2014. “Injection and monitoring at the wallula basalt pilot project.” Energy Procedia 63: 2939–48. https://doi.org/10.1016/j.egypro.2014.11.316.
Mitchinson, Dianne, Jamie Cutts, Dominique Fournier, Annika Naylor, Greg Dipple, Craig J R Hart, Connor Turvey, Mana Rahimi, and Dejan Milidragovic. 2020. “The carbon mineralization potential of ultramafic rocks in British Columbia: a preliminary assessment.” Geoscience BC Report 2020-15/MDRU Publication 452, 25.
Moorkamp, Max, Peter G Lelièvre, Niklas Linde, and Amir Khan. 2016. Integrated imaging of the earth: Theory and applications. Vol. 218. John Wiley & Sons.
National Academies of Sciences, Engineering, and Medicine. 2019. Negative Emissions Technologies and Reliable Sequestration. Washington, D.C.: National Academies Press. https://doi.org/10.17226/25259.
Oldenburg, Douglas W. and Yaoguo Li. 2005. “Inversion for applied geophysics: A tutorial.” Edited by Dwain K. Butler. Investigations in Geophysics 13 (January): 1–85. https://doi.org/10.1190/1.9781560801719.ch5.
Seifritz, W. 1990. “CO2 disposal by means of silicates.” Nature 345 (6275): 486–86. https://doi.org/10.1038/345486b0.
Tikhonov, A. N., and V. Y. Arsenin. 1977. “Solutions of Ill-Posed Problems.” SIAM Review 32: 1320–22. http://www.getcited.org/pub/101760676.
Vasco, D. W., A. Rucci, A. Ferretti, F. Novali, R. C. Bissell, P. S. Ringrose, A. S. Mathieson, and I. W. Wright. 2010. “Satellite-based measurements of surface deformation reveal fluid flow associated with the geological storage of carbon dioxide.” Geophysical Research Letters 37 (3): 1–5. https://doi.org/10.1029/2009GL041544.
Wilt, Michael J., Evan Schankee Um, Edward Nichols, Chester J. Weiss, Gregory Nieuwenhuis, and Kris Maclennan. 2020. “Casing integrity mapping using top-casing electrodes and surface-based electromagnetic fields.” Geophysics 85 (1): E1–13. https://doi.org/10.1190/geo2018-0692.1.
Zhang, Shuo, and Donald J. DePaolo. 2017. “Rates of CO2 Mineralization in Geological Carbon Storage.” Accounts of Chemical Research 50 (9): 2075–84. https://doi.org/10.1021/acs.accounts.7b00334.
