CO2 subsurface mineral storage by its co-injection with recirculating water

Chemical composition of the target subsurface basalts

Cuttings collected from the injection and production wells were characterized by X-ray fluorescence (XRF) to estimate their chemical compositions. The XRF measurements were performed at the laboratories of Isotope Tracer Technologies Europe (IT²E) in Milan, Italy and at King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia.

At the IT²E laboratory, the samples were pulverized with an agate mortar and pestle and their volatile contents were determined by the loss on ignition method. The powders were then dried in an oven at 110 °C overnight, heated in a muffle oven to 1,000 °C, mixed with powdered boric acid and compressed into pellets. The pellets were loaded onto an ARL automatic ADVANT XP spectrometer equipped with a Rh front-window X-ray tube. Analyses were performed using an applied power of 3.0 kW. The count times on identified peaks were 10 s for major elements and 40 s for trace elements. Matrix correction and interelement effects were accounted for using the method of Lachance and Traill³⁰. Analytical uncertainties were determined by analysis of international standards with known compositions. For the major elements, the uncertainties in Si, Ti, Ca and K are less than 3%, whereas for Al, Mn, Mg, Na and P, the uncertainties are less than 7%.

The samples analysed at the KAUST laboratories were first ground to a powder. Approximately 1 g of each sample was mixed with 9 g of XRF Scientific X-ray flux powder composed of 66% lithium tetraborate and 34% lithium metaborate. This mixture was melted at 1,050 °C in an Eagon 2 fusion instrument and poured to form a homogenous pellet that was analysed on a Bruker S8 TIGER machine. The detection limit was <150 ppm for all elements and the uncertainty was <1% for the major elements.

The mineral phases present in the ground well cutting samples were characterized by X-ray diffraction (XRD) analyses conducted at the Exploration Core Labs Department (ECLD) of Saudi Aramco and at IT²E. The XRD analyses at ECLD were run on powders ground with an agate mortar and pestle using a Rigaku Ultima IV powder X-ray diffractometer with CuKα radiation (40 kV, 40 mA) over the 3°–70° (2θ) interval, with a step size of 0.02° increment and a scan speed of 1° s⁻¹. We interpreted the XRD patterns using X’pert HighScore software using specific crystallographic information files for Rietveld refinement and conducted cluster analysis using JADE Pro and its toolkit.

The XRD analyses at IT²E were performed on samples dried at 60 °C and then ground in an agate mortar and pestle. Each sample was then loaded on a onto a polymethylmethacrylate (PMMA) sample holder and placed in a Bruker D8 ADVANCE DaVinci automatic powder diffractometer, equipped with a LYNXEYE detector set to discriminate CuKα_1,2 radiation. The interpretation of the diffractogram for phase identification was carried out by comparison with crystalline phases of the PDF-2 International Centre for Diffraction Data (ICDD) and Crystallography Open Database (COD) databases. A preliminary semi-quantitative estimate of weight fractions was carried out by using the normalized reference intensity ratio method (also known as the Chung method³¹). This method uses scaling factors assigned according to the heights of the characteristic peaks of the different phases. No internal standard was added, so the presence of any amorphous phase was not checked. Quantitative phase analysis of each sample was achieved by the Rietveld profile fitting method as implemented in the Bruker TOPAS V.5 program. This is based on the fundamental parameters approach³². The crystal structure models of crystalline phases considered in the XRD profile fitting for andesine plagioclase, augite pyroxene, clinochlore chlorite, montmorillonite, quartz, laumontite zeolite, calcite and richterite amphibole were taken from the Inorganic Crystal Structure Database (ICSD) Release 2021-2 (ref. ³³) (Supplementary Table 4). Unit cell parameters, scale factors and crystal sizes were allowed to vary for all phases. Atomic coordinates and atomic displacement parameters were fixed while site occupancy factors of octahedral cations and extra-framework species were adjusted, with restraints, in richterite and laumontite, respectively, to account for the crystal-chemical variations of these phases compared with the model. Rietveld profile fitting allowed testing of the presence of mineral phases preliminarily identified and semi-quantified by the reference intensity ratio method³⁴. The resulting chemical and mineralogical compositions of the subsurface drill cuttings are provided in Extended Data Tables 1 and 2.

Chemical analysis of production well fluid samples

Fluid samples were regularly collected from a dedicated outlet port located at the production well. Fluid temperature, pH at the reservoir temperature of 45 ± 0.5 °C and total dissolved solids were measured directly at the fluid sampling port using a Myron L Ultrameter II 6PFC multimeter. The pH electrode was regularly calibrated using Mettler Toledo pH 4.01, 7.00 and 10.01 standard buffer solutions. The uncertainty of the pH measurements was ±0.02 based on replicate analyses of the standard buffer solutions. Measured on-site pH values are reported in Fig. 2. Samples for alkalinity measurement were collected in cleaned 500-ml polyethylene terephthalate (PET) bottles. Further samples were immediately acidified, after filtering using a 0.22-μm syringe filter, by adding 2–3 drops, or approximately 0.01 ml of double-distilled nitric acid, containing 67–69% by weight HNO₃, to 50-ml production well samples in cleaned PET bottles. The PET bottles originally contained drinking water but were rinsed several times first with the production well water before sampling. Both sets of samples were stored in an insulated cooler until transported for chemical analysis. The alkalinity of the first fluid sample was measured by acid titration using the Gran function plot method³⁵. The measured alkalinity was then used to calculate the DIC concentration using PHREEQC³⁶ together with the measured pH and fluid compositions of the major elements at the 45.5 °C subsurface temperature. PHREEQC is a geochemical modelling code designed to perform a variety of aqueous geochemical calculations, including calculations of saturation indices. The second fluid sample was analysed for major cation and Si concentrations by inductively coupled plasma optical emission (ICP-OES) spectroscopy using an Agilent 5110 ICP-OES at KAUST. This spectrometer analysed the compositions of Fe, K, Mg, Si, Ca and Na using wavelengths of 239.563, 769.879, 279.800, 251.611, 315.887 and 568.821 nm, respectively. This instrument was calibrated using Sigma-Aldrich ICP standard solutions of 0.1, 1, 10 and 100 ppm concentration. The limits of detection for Ca and Si were 0.01 and 0.02 ppm, respectively. The analytical uncertainty was <0.1 ppm for all measured elements. To assess the potential contamination from the use of the PET bottles used in fluid sampling, three blanks were prepared using distilled demineralized water. These blanks were prepared and analysed identically to those used for the production well sampling. In each case, the concentration of all measured elements in the prepared blanks was below the respective analytical detection limits. Measured concentrations over time of all sampled fluids are illustrated in Extended Data Fig. 2 and tabulated in Supplementary Table 1.

On the basis of the ICP-OES, pH and alkalinity measurements, the saturation state of the sampled fluids with respect to selected minerals was calculated using PHREEQC³⁶ with its Kinec_v3 database²³. The Kinec_v3 database is the most recently updated database for use with PHREEQC. As the concentration of Al was below the analytical detection limit of our measurements, the concentration of this element was set to be in equilibrium with diaspore in the calculations. This choice was made as diaspore is readily observed to precipitate during experimental studies of basalt dissolution and that chlorite and zeolite minerals, such as clinochlore and laumontite, are more soluble than diaspore at our field conditions. The resulting saturation indices are provided in Supplementary Tables 2 and 3 and Extended Data Fig. 3.

Carbon isotope measurements

Carbon isotope compositions in this manuscript are presented in the delta notation given by:

in which ^13/12C refers to the indicated molar ¹³C to ¹²C isotope ratio, δ¹³C provides the normalized value of this ratio and the subscripts Sample and V-PDB represent the sample of interest and the V-PDB international standard, respectively.

The oxygen isotope values of the solid carbonates are also reported according to the V-PDB standard as

Carbon and oxygen stable isotope compositions of solid carbonate samples and carbon isotope compositions of the CO₂-injected gas were determined at IT²E. The stable isotope analyses were conducted using a Finnigan MAT Delta^plus isotope mass spectrometer using the NBS-18, NBS-19, IT2-20 and IT2-21 standards. The standard deviation associated with the measurements is ±0.3‰ for δ¹³C and ±0.2‰ for δ¹⁸O. The δ¹³C and δ¹⁸O compositions of the carbonates recovered from the disabled production well pump are reported in Extended Data Table 3.

Identification of solids collected from the damaged submersible pump

Fourteen solid samples were removed from the damaged submersible pump and one sample of the dust collected from this pump while drying was analysed for their mineralogical and isotopic composition. The results are provided in Extended Data Table 3. All XRD patterns collected for this purpose are provided in Supplementary Fig. 1.

The δ¹³C composition of the injected CO₂ gas was −37.1‰ and that of the formation fluids before the injection of the CO₂ was −12.9‰. The δ¹⁸O of the circulating water was −3.3‰. Taking account of the equilibrium fractionation factors of these systems^25,26, the δ¹³C and δ¹⁸O values of calcite in isotopic equilibrium with the injected CO₂ and the local formation fluid at 46 °C would be −29.0 and −10.0‰, respectively. These values are close to the δ¹³C and δ¹⁸O values of the recovered disabled pump samples having the highest amounts of carbonate as determined by XRD (Extended Data Fig. 4). This indicates that the cements are fresh carbonate that precipitated from the injected CO₂ and local groundwater. The variability of the δ¹³C and δ¹⁸O values of the rest of the pump samples indicates that, as well as fresh carbonate, the samples also contain a detrital carbonate component (Extended Data Fig. 4). The δ¹³C and δ¹⁸O compositions of this detrital component fall within the ranges of the δ¹³C and δ¹⁸O compositions of carbonates extracted from quartz-carbonate veins associated with the hydrothermal alteration of the basalts, which vary from −3.2 to −6.7‰ and from −19.7 to −21.8‰, respectively (Extended Data Fig. 4).

We attribute the precipitation of carbonate cements on the submersible pump to the supersaturation of carbonate minerals through subsurface mineral reactions rather than degassing of the fluids in the subsurface. This is because the partial pressure of CO₂ in the production well fluids did not exceed 1 bar, whereas the pressure in the submersible pump was not less than 5 bar. Thus, there was no driving force for the degassing of CO₂ within the pump during our study. Furthermore, no CO₂ gas was observed in the recirculated fluid at the injection well, which included a gas trap and a CO₂ concentration sensor at the wellhead.

Calculation of the expected DIC concentration in the absence of chemical reaction

The percentages of injected CO₂ fixed into solid phases by chemical reactions in the subsurface during this study are determined by comparing DIC concentrations measured in the production well fluids to estimates of the DIC concentration expected in the absence of subsurface reactions. Two independent approaches were used to calculate the expected DIC concentrations in the absence of subsurface reactions. These two approaches are described below.

From NaF

NaF was added to the recirculating fluid as a single slug of 2.3 kg on 4 April 2023 to characterize the subsurface flow paths at the pilot test site. The concentration of NaF was then continuously monitored in the production well fluids. Samples for NaF determination were collected in clean PET bottles and stored in the dark until analysis. The concentrations of NaF were measured at KAUST by first adjusting all samples to a pH > 8.7. This was achieved by adding a pH 9 NH₄Cl–NH₄OH buffer solution to the original samples at a 1:1 ratio. The diluted samples were then analysed using a Cary Eclipse Fluorescence Spectrometer at a 512-nm wavelength. The spectrometer used a spectral bandwidth of 5 nm and an excitation of 5 nm. This system has a detection limit of 0.1 ppb and an uncertainty of approximately 1%, as determined by replicate analyses. The measured NaF concentrations at the production well over time are shown in Extended Data Fig. 5a. The NaF breakthrough was detected in the production well fluids 8 days after its addition to the injection well. The NaF concentration then increased to a maximum of 44 ppb on 6 May 2023. This concentration then decreased as the fluids were continuously recycled into the subsurface system. The DIC measured at the production well (Extended Data Fig. 5b) showed a trend similar to that of the NaF tracer, peaking at approximately 33.8 mmol kgw⁻¹. This peak value is 90% lower than the DIC injected into the injection well, which was around 350 mmol kgw⁻¹.

To analyse the NaF tracer and DIC production profiles, the first-moment analysis technique was applied. This method adjusts for the recycled chemicals within the production curves by incorporating the age distribution, also known as the residence time distribution, E(t) (refs. ^{37,38,39,40,41}), defined as

in which Q_p represents the fluid production rate (m³ day⁻¹), C_m(t) denotes the chemical concentration measured at the production well (kg m⁻³) and m_inj designates the total mass of the chemical injected (kg). The unit of E(t) is inverse time (day⁻¹). The first moment, τ, is used to calculate the mean residence time between the injection and production wells, defined by

Because the chemicals are reinjected, the measured concentration history reflects a combined effect of the initial injection and the continuing recycling of the produced injectate. Moment analysis, which addresses the response to a one-pass slug injection, requires the removal of the recycling effect to accurately determine chemical residence time and swept pore volumes. The convolution integral can be applied to isolate the response of the initial tracer injection, resulting in a single-pass tracer return profile³⁷. The adjusted concentration, C_adj(t), thus represents the concentration at which the chemical is removed from the production fluids before being reinjected, such that

To account for the effects of fluid recycling, the deconvolution operator was applied to both the NaF tracer and the DIC concentration profiles. This step was necessary to isolate the first-pass breakthrough of the injected fluid from subsequent recirculation. The resulting adjusted concentrations of the NaF tracer and DIC concentration are shown in Extended Data Fig. 5a,b, respectively. The difference between the measured concentrations and these adjusted profiles is attributable to the recirculation of the chemicals.

Because NaF was injected as a single pulse, whereas CO₂ was injected continuously, an extra adjustment was required to ensure that the two chemicals share the same injection profile for a direct comparison. To achieve this, the continuous-response NaF concentration, denoted as C_Cont,NaF(t), was calculated by convolving the pulse injection profile of NaF with the continuous CO₂ injection function, Inj_CO2(t). Therefore, the convolution operator was applied to the pulse NaF injection profile, C_Pulse,NaF(t), to ensure that the response of NaF mimics a continuous injection scenario, such that

The effective pore volume, PV, swept by these chemicals can be estimated from

Given that NaF is a conservative and non-reactive tracer, an estimate of the DIC concentration was made by assuming that carbon behaves similarly to NaF in the absence of subsurface reactions. This assumption requires that diffusion of CO₂ and the tracer if it diffuses into the rock matrix do so at the same rate. This would be likely, as they are both dissolved aqueous phases. The residence time distribution for NaF was then convolved to match the injection profile of CO₂. The calculated DIC concentration, DIC_cal,NaF (mol kgw⁻¹), in the production well fluids in the absence of subsurface reactions, estimated from the convolved NaF profile, E_NaF(t), is given by

in which M_CO2 is the total mass in mol of the injected CO₂.

The expected DIC concentrations in production well fluids, assuming no subsurface reactions, are compared with the measured DIC concentrations in Fig. 3a. This comparison is used to estimate the percentage of CO₂ mineralized. The data indicate progressive mineralization of the CO₂ up to 70 ± 5% within eight months after stopping CO₂ injection.

Several studies suggest that NaF can sorb onto negatively charged silicate mineral surfaces at neutral to acid pH, leading to the non-conservative behaviour of this tracer^42,43. Notably, one study⁴³ found that the concentration of NaF decreased substantially in formation waters when co-injected together with CO₂, NO₂, SO₄ and O₂ into the subsurface. This study concluded that NaF was not suitable as a tracer to quantify the fate of CO₂ injected into geologic formations and ascribed this decrease to the sorption of the tracer to mineral surfaces. Two lines of evidence indicate that such an effect did not occur in our study, First, most of the NaF injected during our study is accounted for in the recovered production well fluids. Second, the fraction of carbon mineralized determined from NaF and SF₆ are identical, suggesting that any effect of NaF sorption on subsurface mineral storage in our study region is negligible (Fig. 3). In this regard, it should be noted that, if some NaF had selectively sorbed onto subsurface minerals in our study, this process would have lowered the concentration of recovered NaF and thus the estimated DIC concentration in Fig. 3a. As a consequence, if selective NaF sorption had occurred, the percent carbon mineralized in this study would have been underestimated and the true percentage of CO₂ mineralized would have been higher than reported in the main text. The difference in NaF behaviour in the present study compared with the previous study⁴³ could potentially be attributed to the co-injection of oxic gases with the NaF tracer. NaF is commonly known to oxidize to a non-fluorescent product in the presence of a variety of oxidizing agents⁴⁴.

From SF₆

Following the approach of Matter et al.¹¹, SF₆ was used as an inert tracer to estimate the percentage of the CO₂ mineralized over time. To this end, SF₆ was co-injected with CO₂ at a constant SF₆ to CO₂ mass ratio of 1.05 × 10⁻⁶:1. Both the CO₂ and the SF₆ injection were stopped on 7 July 2023. The SF₆ was originally stored in gas cylinders. The gas was added to the recirculating water stream using a mechanical flow regulator connected to an AALBORG electronic mass gas flow meter recalibrated for SF₆. Fluid samples for SF₆ analysis were collected in pre-cleaned 1.0-l glass bottles provided by Spurenstofflabor. The sample bottles were completely submerged in a 25-l plastic container fed with water from the sampling port of the production well for sample collection. Concentrations of SF₆ were measured at Spurenstofflabor. Sample preparation included the extraction of a 40-ml aliquot injected into glass bottles filled with N₂ gas. After equilibration by shaking for 30 min, 10 ml of headspace gas was separated for further dilution with N₂ gas. Aliquots were analysed by gas chromatography using an electron capture detector. The detection limit is below 0.1 fmol l⁻⁶ SF₆. The reproducibility is <2% for air samples and <10% for water samples. The reliability of the measurements was determined by analysing two random duplicate samples (that is, DW-3: 6 July and DW-1: 9 August and 17 and 24 March) with and without the final dilution step. The differences between the results were always less than 10%.

The measured SF₆ concentrations are provided in Fig. 3b. The SF₆ concentration in the production well waters first increased, maximizing during 9–12 July, or shortly after the CO₂ and SF₆ injection was stopped. A slow decline in the concentration is observed after this time. This can be attributed to a continuous dilution of this tracer over time. Dilution of the fluid is favoured owing to the lowering of the water table at the production well. The water table at the production well decreased by approximately 150 m owing to pumping. This lower water table height drives groundwater from surrounding subsurface to the production well as well as the injection well, diluting the concentrations of both SF₆ and NaF.

The measured SF₆ concentrations are compared with the corresponding DIC concentrations in Fig. 3b. Similar to the analysis for NaF tracer, the SF₆ tracer is non-reactive, therefore an estimation of DIC concentration was made assuming that carbon behaves similarly to SF₆ in the absence of subsurface reactions. The calculated concentration, DIC_cal,SF6 (mol kgw⁻¹), in the production well fluids, in the absence of subsurface reactions, is estimated from the adjusted residence time distribution of SF₆, E_SF6(t), such that

The calculated DIC_cal,SF6 concentration matched closely with that of the observed DIC until September 2023. Afterwards, the observed DIC concentration decreased faster with time than that of SF₆. As the SF₆ is non-reactive, the difference between the two curves in Fig. 3b can be attributed to the precipitation of the injected carbon in the subsurface, indicating about 70 ± 5% mineralization by April 2024, as shown in Fig. 3b.

Reservoir characterization from the NaF tracer

The NaF tracer concentrations recovered from the production well revealed a complex distribution of fractures and matrix systems within the reservoir. This concentration profile suggests the presence of at least two distinct flow systems in the reservoir with distinct permeabilities. We conducted a detailed analysis to compare the observed tracer profiles against simulations. Initial attempts to match the tracer data using a single-channel model were unsuccessful⁴¹, as shown in Extended Data Fig. 6a. However, when a dual-channel model was applied, the tracer profile could be accurately matched (Extended Data Fig. 6b). This dual-porosity scenario implies that fluid flow within the reservoir occurred through two main pathways: (1) highly permeable conduits such as faults, which led to a rapid tracer breakthrough at the production well, and (2) through a more extensive network in the matrix composed of smaller, interconnected and less permeable fractures.

To quantify the effective swept pore volume, the method of moments was used³⁸. This approach uses moment analysis, requiring a full production profile of the tracer over time. The extended tracer concentration profile at the production well was estimated using an exponential decline model, which, when plotted on a semi-log scale, results in a straight line, as depicted in Extended Data Fig. 6c. The outcome of this analysis, coupled with an uncertainty analysis, revealed that the total effective pore volume ranged between 24,000 and 43,000 m³. Of this, approximately 10% is attributed to the fast-flow regions with a mean residence time of 50–65 days, whereas the remaining 90% is associated with the low-permeability medium, with residence times spanning 255–445 days.

In total, 1.55 kg of NaF was recovered from the production well fluids over the course of this study, up until 21 April 2024. After accounting for that reinjected through the circulation of the production well fluids, the adjusted NaF recovery was 1.06 kg. This compares with a total of 2.30 kg of NaF originally injected into the system. Thus, 46% of the originally injected NaF was recovered before 21 April 2024. It seems likely that the unrecovered tracer is primarily related to fracture–matrix diffusion and not sorption. The tracer initially diffuses from fractures into the low-permeability basalt matrix and then returns slowly from matrix to fractures over multi-month timescales. Two independent observations support this interpretation. First, the adjusted cumulative mass of NaF continues to increase for months after the breakthrough peak, indicating continuing release of stored tracer. Second, a semi-log linear tail of the recovery of this tracer is observed (Extended Data Fig. 6c). Such a single-timescale exponential tail is the classic signature of mobile (fracture)–immobile (matrix) mass transfer, not of equilibrium adsorption.

Estimating total CO₂ storage capacity of the field site and upscaling potential

The total mass of CO₂ that could be mineralized in this system is challenging to estimate. Field observations suggest that the total extent of carbonation of a basalt is limited by the availability of pore space rather than the mass of available reactive rocks. This is because water–CO₂–basalt interaction leads to an increase in the overall volume of solids, reducing the pore space and blocking flow paths over time. Also, this process provokes the formation of numerous hydrous silicate minerals, including clay and zeolite minerals, as well as CO₂-storing carbonate minerals^45,46. Overall, the volume of carbonate minerals produced by water–CO₂–basalt interaction may be only 20% or less of the total volume of precipitated secondary minerals⁴⁵. This further limits the space available for mineral carbonation reactions in basalts. The overall efficiency of mineral carbonation can vary by adjusting the fluid injection rates and the composition of this injected fluid, including CO₂ concentration⁴⁵. Each of these factors could be varied during a subsurface mineral carbon storage effort. If all of the effective flow volume was available for just calcite precipitation, the available 24,000–43,000 m³ of effective pore space in our pilot system would be sufficient to contain 22,000–40,000 tons of mineralized CO₂ in total. This value is largely an overestimate owing to the likely formation of hydrous silicate secondary phases, which will consume some of the available pore space. Also, it is unlikely that all of this pore space would be available for secondary mineral precipitation, as fluid pathways may be blocked. Some of these challenges and limitations might be overcome by fracking the subsurface rock formation⁴⁷, but this opinion has been relatively unexplored so far in basaltic systems.

One potential advantage of using the approach of this study for subsurface storage is that it limits the consumption of energy needed for subsurface storage. Energy is required to pump water from the subsurface, for which the amount of energy required is a function of the depth of the water table, flow rate and the efficiency of the pump. The injection of CO₂, using our approach, requires CO₂ to be delivered at a pressure that is high enough to exceed the hydrostatic pressure at its delivery depth in the wellbore. With further pressure that could be needed to enhance CO₂–water mixing by means of a static mixer or a nozzle, the total CO₂ delivery pressure at the surface would be around 12–14 bar. This range is 8–16 times lower than what is typically needed for conventional CCS. The generated water-charged CO₂ injection is essentially driven into the system by gravity. Thus, the overall energy penalty by making use of our CO₂ storage approach may be less than for conventional storage requiring injection of fluids into deeper and higher-pressure systems. Nevertheless, this proposed technology is seen as complementary to conventional CCS and not an alternative, as the geological conditions dominate.