The Evolution of Paleo-Porosity in Basalts: Reversing Pore-Filling Mechanisms Using X-Ray Computed Tomography

Often carrying a high-volume fraction of vesicles, basaltic rocks can be an important reservoir horizon in petroleum systems, and are considered an excellent candidate for CO2 storage by in situ mineral trapping. The frequency of amygdaloidal basalts in many sequences highlights the prevalence of mineralisation, but when the vesicle network has been filled, the basalts can act as impermeable seals and traps. Characterising the spatial and temporal evolution of the porosity and permeability is critical to understanding the petro-physical properties and CO2 storage potential of basalts. We exploit X-ray computed tomography (XCT) to investigate the precipitation history of an amygdaloidal basalt containing a pore-connecting micro fracture network now partially filled by calcite as an analogue for CO2 mineral trapping in a vesicular basalt. The fracture network likely represents a preferential pathway for CO2-rich fluids during mineralisation. We investigate and quantify the evolution of basalt porosity and permeability during pore-filling calcite precipitation by applying novel numerical erosion techniques to “back-strip” the calcite from the amygdales and fracture networks. We provide a semi-quantitative technique for defining reservoir potential and quality through time and understanding sub-surface flow and storage. We found that permeability evolution is dependent on the precipitation mechanism and rates, as well as on the presence of micro fracture networks, and that once the precipitation is sufficient to close off all pores, permeability reaches values that are controlled by the micro fracture network. These results prompt further studies to determine CO2 mineral trapping mechanisms in amygdaloidal basalts as analogues for CO2 injections in basalt formations. Different mechanisms of mineral precipitation affect pore connectivity and permeability evolution; When precipitation instantaneously fills all pores, porosity is disconnected quite sharply. When precipitation fills the porosity with a constant layer of calcite at each time step, porosity is disconnected more slowly. Once precipitation closes all pores, permeability evolution is controlled by the micro fracture network. Different mechanisms of mineral precipitation affect pore connectivity and permeability evolution; When precipitation instantaneously fills all pores, porosity is disconnected quite sharply. When precipitation fills the porosity with a constant layer of calcite at each time step, porosity is disconnected more slowly. Once precipitation closes all pores, permeability evolution is controlled by the micro fracture network.


Introduction
Vesicular basalts are common in many basaltic sequences. The vesicles generated by gas exsolution from the melt are often filled by the precipitation of secondary minerals such as quartz or calcite during hydrothermal or later ground water circulation, depending on the fluid chemistry percolating through the lava. Basalts are therefore considered to be excellent candidates for CO 2 storage by mineral trapping because they make up ~ 10% of Earth's continental surface and the vast majority of seafloors (Goldberg et al. 2008;Matter et al. 2016;Oelkers et al. 2008), and are often capped by impermeable sediments, offering additional advantage towards preventing post-injection leakage (Goldberg et al. 2009). These abundant and accessible units contain up to 25% by weight of calcium, magnesium and iron, causing high reactivity with CO 2 and facilitating the production of stable and nontoxic pore-filling minerals (Goldberg et al. 2009(Goldberg et al. , 2008Luhmann et al. 2017;Matter et al. 2016;Wu et al. 2021). This reactivity translates into carbonate mineralisation on the timescale of weeks to years, compared to hundreds or thousands of years for sedimentary rocks (Matter et al. 2016;Menefee et al. 2018). However, understanding the petrophysics, and characterising the porosity, pore connectivity and permeability during mineral precipitation are critical if we are to fully exploit the storage potential of basalts.. Cementation has been widely investigated using numerical simulations in porous network models. Initially applied on packed beds of sphere as simplest pore geometry to model (Bakke and Øren 1997), increasing computational power has now seen these methods applied to complex porous media such as tight sandstones, and to microporosity in carbonates (Mehmani and Prodanović 2014;Prodanović et al. 2015). Other studies have used XCT images as input for level set methods (Prodanović et al. 2013) or lattice-Boltzmann methods to similar geological processes such as chemical dissolution and precipitation (Kang et al. 2014;Miller et al. 2017;Zambrano et al. 2018). While they show porosity-permeability trends similar to those found in nature (Bernabé et al. 2003), these studies do not consider how precipitation may occur heterogeneously within it. Although several studies have characterised vesicles morphology and connectivity in basalt using X-ray computed tomography (XCT) (e.g. Couves et al. 2016;Degruyter et al. 2010;Giachetti et al. 2011;Polacci et al. 2006;Song et al. 2001), exploiting the method to calculate porosities, perform flow simulations and estimate CO 2 storage through mineral precipitation has been limited (Callow et al. 2018). Low-resolution (~ 52.4 µm) XCT has showed that olivine-tholeiite Icelandic basalts with 2.05-8.76 vol.% connected porosity had a reactive surface area of 0.10-0.33 m −1 (~ 0.33 Gt of CO 2 ) for a 1.8 km 3 basalt volume, and found an XCT derived permeability that was two orders of magnitude higher than the one calculated from hydromechanical tests (Callow et al. 2018). By artificially precipitating a ~ 50-µm-thick film of calcite to investigated the role of pore diameter reduction and pore clogging during in situ precipitation, the authors suggest that these mechanisms may be insufficient to account for the observed drop in permeability (Callow et al. 2018). Mass balance models, which preserve the mass balance of the available precipitate, are not uncommon but they have been used mostly to investigate the effects of CO 2 -rich brines into cement fractures (Rod et al. 2020), the role of fluid transport regimes on the precipitation of secondary carbonates in fractured basalt as well as in the fundamental understanding of mineral carbonation (Adeoye et al. 2017;Haque et al. 2019).
Here, we explore these mechanisms more fully using an amygdaloidal basalt where the vesicles have been infilled with calcite precipitated from CO 2 -rich fluids. This is an analogue for the end-product of CO 2 injection into basalt for the purposes of CO 2 storage. Using novel "back-stripping" methods of 3D image analysis at a higher resolution we progressively dilate the paleo-porosity, replicating a reverse of the calcite precipitation process as analogue to mass balance models, providing insights into the temporal evolution of pore connectivity and permeability during calcite precipitation.

Methods and Analysis
The sample was obtained from a road cutting exposure on the Isle of Mull, UK (Ordnance Survey NM 51196 53869, 56.610441-6.056248), and comes from the olivine basalt Mull Plateau Lava Formation dated at ca. 60.5 Ma (Emeleus et al. 2005). The host basalts are heavily altered adjacent to the calcite vein, although plagioclase appears relatively fresh. Clumped isotope analysis of a calcite-filled fracture from this locality records a calcite precipitation temperature of 83-85 °C and an isotopically depleted parent fluid δ 18 O of ca. −10‰, indicating calcite precipitation from meteoric water (MacDonald et al. 2019).

XCT Data Analysis
X-ray computed tomography (XCT) data were acquired on a Nikon XT H 320 LC laboratory scanner (University of Strathclyde) using a transmission X-ray source, with an X-ray energy of 110 kV and a beam current of 109 µA. 3141 projections were acquired (no frame averaging) with an exposure time of 1.42 s per projection. The reconstructed voxel size was 12.58 µm. Data analysis was performed using a combination of Avizo® 9.3 and ImageJ.

Image Processing
A volume of interest (VOI) (1142 × 863 × 1072 voxels) was extracted from the XCT dataset to remove edge effects, and a median filter (kernel size 2) applied for noise reduction. The Trainable Weka Segmentation (Arganda-carreras et al. 2016) (ImageJ®) was applied using the random forest algorithm to segment and classify the dataset into four phases: the volume and connectivity of each was then analysed in Avizo® using labelling (neighbourhood 26) and label analysis (neighbourhood 26) operators. The volume fractions of vesicles (i.e. maximum pore volume) and fractures were calculated from the filled and unfilled volumes, 1 3 as was the maximum paleo-porosity. We did not define or separate individual amygdales based on their pore throats, as pore geometries are complex and the throats are large, making arithmetic-based separation a poor fit to the geological understanding. The individual measures presented are for each unseparated cluster of amygdales.

Back-Stripping Method
Back-stripping-method 1: instantaneous pore-scale dissolution where we assume any filled pore or fracture is dissolved within one time step, as soon as it is connected to the pore network. To understand how the calcite filled the paleo-porosity, we imposed a workflow based on morphological dilation (Fig. 1a). We define t 0 to be the present, and select all voxels of the porosity that would be accessible to a fluid flow parallel to the Z orientation of the VOI (Zϕ c ), defining this as the surface connected porosity, t 0 Zϕ c , prior to backstripping (e.g. current connected porosity). The percentage connectivity of the porosity that is connected (or accessible) is defined as % conn = t 0 Zϕ c tot × 100 , where ϕ tot represents the total volume of paleo-porosity. t 0 Zϕ c is then dilated by a unit voxel (5 in this case) using a spherical dilation operator, within ϕ tot (basalt volume remains unchanged). After the dilation, the t 1 Zϕ c is calculated. This will include the t 0 porosity, the voxels added by the (a) (b) (c) Fig. 1 Sketch illustrating the different methods of back-stripping (a, b) and the simulated precipitation (c), for two amygdales partially filled by calcite (C), and including a central pore filled by air (P). The air-filled paleo-porosity (t 0 ) is indicated by a solid black line dilation and any disconnected air-filled t 0 porosity that is reconnected by the dilation. This process is then repeated n times until all the calcite is removed. This back-stripping was also performed for t 0 Xϕ c and t 0 Yϕ c (for X and Y parallel flow, respectively). For each step of back-stripping, the specific surface area (SSA, surface area/volume) of the newly connected network is calculated from the measured surface area (SA) and volume (t n Xϕ c , t n Yϕ c , t n Zϕ c ). This enables us to understand how calcite precipitation progressed through the sample. Assuming precipitation is isotropic, a dilation step is the reverse of a period precipitation (hence the term back-stripping). While pore volume could be reconnected after a spherical dilation with r = 1 voxel, back-stripping with r = 5 (~ 63.7 µm) was sufficient to track the process and is equivalent to assuming a longer time step or faster rates of precipitation per unit time. We assume that fractures are open throughout. Permeability here is calculated using LBflow lattice-Boltzmann fluid flow simulator (Llewellin 2010a, b), which has been calibrated over a wide range of porosity against explicit models for the permeability of body-centred cubic (Llewellin 2010b) and simple cubic packs of spheres (Wadsworth et al. 2017). The LBflow model takes the segmented pore phase in the total domain as an input and discretizes the pore phase into a lattice of fluid nodes (one fluid node per voxel). Fluid mass parcels are numerically propagated through the lattice positions with time and the result of collisions are computed on a D 3 Q 15 lattice. In this framework, we apply a numerical pressure gradient of ∇p = 0.01Pam −1 using a pore fluid with viscosity f = 1.8205 × 10 −5 Pas and density f = 1.2047kgm −3 , corresponding to air at ambient conditions. The simulation converges at steady-state flow determined via a criterion that the average fluid speed across the lattice does not vary by more than a 10 −5 factor over 50 iterations twice. The back-stripping processes are performed on the full volume, but because of computational limitations we run all LBFlow simulation on a cubic subdomain with edge length L (300 and 400 voxels) at a fixed position in XYZ. The two subvolumes are centred within the same XYZ coordinate system. Permeability was calculated for each pore configuration and mutually perpendicular directions of simulated fluid flow. The LBflow algorithm produces a spatial distribution of fluid velocities at steady state, u . The average fluid velocity in the direction of the applied pressure gradient is given by ⟨u⟩ = tot u which can be used in Darcy's law to determine the permeability k in a given direction via k = −∇p f ∕⟨u⟩ . The input properties f and f and condition ∇p are checked to ensure flow in the creeping low-Reynolds number regime and low-Mach number regime. The Reynolds number is given by Re = f L⟨u⟩∕ f and the Mach number is given by Ma = ⟨u⟩∕c where L is the domain edge length and c is the speed of sound in the fluid medium. We find that 10 −8 ≤ Re ≤ 10 −5 and 10 −12 ≤ Ma ≤ 10 −9 ; regimes in which Darcy's law is valid and in which k does not depend on the fluid properties chosen.
Back-stripping Method 2: constant dissolution rate where we assume that a constant thickness of calcite is dissolved from the fluid accessible surface at any time step (Fig. 1b). We assume fractures are open throughout the back-stripping. Here, t 0 Zϕ c is dilated by n voxels (here n = 50 to allow complete back-stripping in a computationally feasible number of steps) using a spherical dilation operator, within ϕ tot (basalt volume remains unchanged). After the dilation, the t 1 Zϕ c is calculated and the process is repeated n times until all the calcite is removed, as the other back-stripping approach. Calcite was stripped out in 10 dilation steps. The connected porosity did not reach 100% because of the presence of disconnected pores which did not take part to the flow. Because of the larger number of steps required to strip back out the calcite, LBflow was considered too computationally extensive: Avizo has a comparable simplified permeability version, which solves Darcy's flow. The two methods have proven to yield comparable results, and therefore, permeability for method 2 and 3 was calculated using Avizo permeability simulator for computational efficiency.
Method 3-Simulated precipitation of paleo-porosity where we assume that, for a given time, we have a constant thickness of precipitation across the entire surface accessible to the fluid (Fig. 1c). We assume fractures remain open throughout. Precipitation is performed by applying a spherical erosion operator (n = 1) to the paleo-porosity, with ϕ c and permeability recalculated every fifth erosion for X and Y parallel flow, and every second erosion for Z parallel flow. Permeability was reduced to 0 mD by 35 erosion time steps in all directions.

Microstructure
The XCT analysis shows that the amygdales are preferentially elongated parallel to the Y direction of the XCT data and are almost entirely filled with calcite. The amygdales (primary porosity) are not all fully filled with calcite ( Fig. 2a, b), with some pores within the amygdales representing incomplete calcite precipitation. The sample also has secondary porosity in the form of lacy porosity (regions of calcite with complex fracture networks) and fractures ( Fig. 2a, b, e). We can identify two generations of cross-cutting fractures, some of which are partially filled with calcite ( Fig. 2c, d). The amygdales form disconnected clusters (Fig. 2e), with almost all (99.8 vol.%) clusters having volumes greater than 1 mm 3 (largest has a volume of 54.8 mm 3 , equivalent diameter of 4.7 mm) (Fig. 2f). The scanned palaeo-porosity of the basalt is 18.49 vol.%: 18.19 vol.% amygdales (0.19% vol.% still air-filled), and 0.40 vol.% introduced by the fractures (0.24 vol% are partially air-filled).

Back-Stripping the Calcite: Reversing the Precipitation in Time
Method 1. Results from the first back-stripping approach indicate that the percentage connectivity at t 0 is equal to ~ 76.5%, along the Z direction (Table 1, Figs. 3,4). The 3D rendering shows that quite a lot of paleo-porosity is missing as it is not connected (Figs. 3,5). At t 1 , the percentage connectivity of the paleo-porosity is ~ 99%, and by t 2 is 99.9%, indicating that almost the entire paleo-porosity is connected (Table 1, Figs. 4,5). In the 3D renderings (Fig. 5), it is possible to note that this is caused by new amygdales becoming progressively connected through the back-stripping process.
When applying the back-stripping method along the other two directions, we found that the paleo-porosity is more connected initially, with ~ 85.7 and 93.3% connectivity in X and Y, respectively, where Y represents the elongation direction of the amygdales (Figs. 3, 5). After t 1 , the data follow the same trend as in Z, and the paleo-porosity becomes connected after two steps of dilation in each direction (Figs. 4,5). As compared to the sample volume, paleo-porosity was decreased by ~ 4% in connectivity with increase in time of precipitation, decreasing from 18.4 to 14.1 vol.%, along the Z direction (Table 1, from t 3 to t 0 ). Along X and Y, paleo-porosity decreased by 2.6 and 1.2 vol.%, respectively.
For each step of dilation, we also calculated the surface area and the specific surface area of the connected paleo-porosity. The specific surface area provides an estimate for calcite precipitation. Values of specific surface area range from 90.66 to 144.96 m −1 . 2 a, b XCT vertical slices of XZ and YZ planes, of the scanned samples: vesicles (V), lacy porosity (L), inter-calcite porosity (IC) and fractures (F) are present within the basaltic matrix (BM). c 3D render of the partially filled fractures (blue) and air-filled fractures (red), cross-cutting each other. d 3D render of residual inter-clast pore indicating incomplete precipitation intersected by fracture (both rendered in red) within a calcite-filled amygdale (transparent grey colour). e 3D render of the individual unconnected amygdale clusters. f 3D render of the largest amygdale cluster (0.09 to 0.14 mm −1 ). Back-stripping the paleo-porosity highlights how pore filling by ~ 188 µm (from t 3 to t 0 in Table 1) along Z causes the available specific surface area to decrease by 0.05 mm −1 , equivalent to nearly 36% of the original specific surface area. Similar values are found in the other two directions (Table 1). Method 2. We originally imposed a morphological dilation of 5 voxels in each direction, equals to ~ 63 µm, considering instantaneous pore-scale dissolution. However, precipitation of calcite may have occurred with constant precipitation rate where a constant thickness of calcite is deposited to the fluid accessible surface at any time step. The results from this method (method 2) indicate that the X direction is the most connected at time 0, with a percentage connectivity of 2.5% (Fig. 6a, b). In Y and Z direction, porosity is connected only after 3 and 4 events of dilation (hence dissolution), respectively (~ 1885-2514 µm), testifying basalt heterogeneity. The connectivity along the X direction starts increasing once the back-stripping is sufficient to dissolve the calcite in the paleo-amygdales and connect them to the pore network (Fig. 7). Similarly, when Table 1 Tabulation of ϕ c as percentage connectivity (ϕ conn % ), its volume in mm 3 , the specific surface area (SSA), and as against the volume of the sample (ϕ conn in the sample), for each step of back-stripping (t 0 -t 3 ) and for the three directions simulating calcite precipitation (method 3), we found that precipitation takes longer along the X direction to disconnect the connected porosity, as compared to the other two directions, testifying the larger degree of connectivity along the X direction provided by the fracture network (Figs. 6c, 8).

Permeability Results
Back-stripping Method 1. Results indicate that there is no permeability along the Z direction, considering the Z direction of back-stripping, at the scale of the subvolumes. However, this pore network is still permeable in the other two directions of fluid flow. When calculating the permeability along the X and Y direction, for the relative back-stripping data, permeability stays constant with time, increasing only after the first back-stripping  (Fig. 9).
To corroborate the results from the LBflow model, we also skeletonized the pore network for each back-stripping datasets along the Z and X directions, to understand how the pore coordination number and the number of skeletons evolve with precipitation (see Suppl. Material Section S2). We found that there is no significant change in the frequency of pore coordination numbers with precipitation, nor in the evolution of the number of segments per skeletons (Figs. S2, S3).
Back-stripping Method 2. We calculated the permeability for this method to check how a different precipitation mechanism might affect permeability and CO 2 storage capacity (Fig. 10a). We found that the permeability on the subvolume along the X direction increases after the first back-stripping event when the pores become open, and then it remains constant through relative time (Figs. 10a, 11a). This means that as soon as calcite precipitation is sufficient to close the fractures, the sample becomes less permeable. The permeability along the Y direction in the subvolume is equal to zero for the current time, as porosity is not connected, but it then increases and remains constant with the back-stripping and increasing connected porosity. The permeability along the Z direction is very low as compared to the other two directions, indicating that the subvolume remains impermeable throughout the back-stripping process. Dilating the remaining open-air porosity does not change the permeability for any direction. We also calculated the permeability evolution along the X direction for 4 relative time steps during the simulated event of precipitation (Fig. 10b). We found that, after the first precipitation, permeability is drastically reduced, in agreement with findings from the permeability evolution back through time during the back-stripping process (cfr. 10a). 3D visualisations of the subvolume used for the permeability for different steps of the back-stripping process and the simulated precipitation enabled us to identify that the critical change in permeability during the two processes Fig. 4 The graph shows calcite precipitation through relative time. t n ϕ c is reported (in mm 3 and as percentage connectivity (ϕ c /ϕ tot ) against time of precipitation for each direction. The precipitation time, right to left (t 3 to t 0 ), is represented by different stages of morphological dilation, from left to right. 3D renderings for each direction are represented in blue models is caused by pore openings (during the back-stripping) and pore closing and eventually snapping-off (during the precipitation) (Fig. 11), indicating that permeability during the final stages of precipitation and initial stages of back-stripping is controlled by the micro fracture network, as testified by the results of modelling permeability along the fracture network only.

Porosity and Permeability Results From XCT Morphological Back-Stripping
XCT analysis reveals that this amygdaloidal basalt is characterised by a heterogeneous distribution of amygdales, elongated in one direction. The XCT data also revealed the presence of two sets of fractures, one of which seems to be only partially filled with calcite. We found that these amygdales are connected along the Y direction, without taking into consideration the fracture network. This means that precipitation of calcite could have occurred also prior to the fracture network.  5 Temporal evolution of ϕ c (blue) with calcite precipitation for each direction of flow and for each step of back-stripping. As we reverse calcite precipitation (from left to right), pore connectivity increases (refer to Fig. 3). Vice versa, as calcite precipitation occurs, pore connectivity decreases (right to left). Numbers refer to percentage connectivity (ϕ c /ϕ tot ). The remaining space of the sample is filled by the basaltic matrix, that is not represented here. Red circles identify new connected pores as we reverse calcite precipitation. Vertical side is 13.6 mm long In this study, we are investigating how precipitation affected pore connectivity and permeability in this sample, to understand what factors may be involved in the process and may need to be taken into consideration for upscaling models of CO 2 sequestration in basalts. We have simulated different mechanisms of back-stripping and precipitation through time. Our XCT analyses show that the total paleo-porosity, combining the amygdales and the Fig. 6 a The graph reports the evolution through relative time of ϕ c /ϕ tot during the back-stripping for each direction. The precipitation time, right to left, is represented by different stages of morphological dilation, from left to right. 3D renderings for the X direction are represented in blue models for 3 different time steps. b The graph reports the evolution through relative time of the volume (in mm 3 ) of ϕ c during the backstripping for each direction. Note the along Y and Z direction porosity is not connected despite being dilated (refer to a). c The graph represents the evolution of the simulated precipitation over time. The evolution of ϕ c /ϕ tot is reported through relative precipitation time for each direction. From right to left, the volume percentage decreases with precipitation fracture space, is ~ 18.4 vol% (as the volume percentage within the sample). This is in good agreement with values found in other basaltic lava flows, ranging between 5 and 40 vol% (Callow et al. 2018;Franzson et al. 2008;Kanakiya et al. 2017). In the first method, we considered the pores to be dissolved instantaneously by the back-stripping process as soon as they are open. We found that current connected porosities, reflecting the final stage of Fig. 7 Temporal evolution of connected porosity (blue) during the back-stripping method 2 for each direction and for four steps of back-stripping (as indicated by the values of ϕ c /ϕ tot ). The porosity that is not connected has been rendered in grey. Note that there is no connected porosity along Y and Z directions in the current time considering this method: this is because only the fractures connect from side to side. Note that, as we dissolve calcite (from left to right), pore connectivity increases. The remaining space of the sample is filled by the basaltic matrix, that is not represented here. Vertical side is 13.6 mm long Fig. 8 Temporal evolution of connected porosity (ϕ c ) in blue with calcite precipitation along the X direction for four steps of the simulated precipitation. The percentages indicate the volume of surface connected porosity against total porosity (ϕ c /ϕ tot , cfr. to Fig. 6c). The remaining space of the sample is filled by the basaltic matrix, that is not represented here. Vertical side is 13.6 mm long precipitation, vary between 14 and 17.1 vol.%, along different directions, reflecting ~ 4 to 1.3 vol% decrease from fully connected porosity due to in situ precipitation and indicative of basalt heterogeneity. This is in agreement with studies from the literature proposing that geochemical alteration and in situ mineral precipitation can cause up to 10 vol.% decrease in porosity content (Callow et al. 2018;Luhmann et al. 2017). Translating these values as percentage of total porosity gives percentage connectivity between 76.6 and 93.3%, values which are considerably higher than that found at the CarbFix site by Callow et al. (2018). This is, however, in good agreement with the lower porosity content found at the CarbFix site, and overall higher in situ precipitation. However, if we do not consider that pores are dissolved as soon as they are open, but rather than a constant thickness of calcite is dissolved at each time step as in method 2, we found that the current open-air porosity is only connected along the X direction, giving 2.5% of the total paleo-porosity. This value is more consistent with what found at the CarbFix site.
It has been shown experimentally that carbonate and secondary minerals are likely to precipitate in diffusion-dominated areas of fractures (Luhmann et al. 2017;Wu et al. 2021). When considering that in situ precipitation may have localised within fracture-like pores and initiated at fracture-wall interfaces, our results show that mineralisation causes porosity to become disconnected more slowly (Figs. 6a, b, 7). Precipitation along the fracture network would take longer and more flux input would be needed to disconnect porosity, particularly along the X direction, along which the fracture network is distributed. This confirms some other experimental results according to which precipitation of secondary minerals on the surface would reduce the available reactive surface area, thus slowing down dissolution rates of basalts with time and hence CO 2 mineralisation Wu et al. 2021). When we consider constant dissolution rate, where a constant thickness of calcite is dissolved from the fluid accessible surface at any time step, and fractures remain open, our results indicate that along the X direction dissolution of calcite back Fig. 9 Permeability-porosity relationship for four back-stripping steps as calculated by the LBflow model applied on two subvolumes (300 voxels, 400 voxels), as they have been dilated to mimic calcite precipitation in reverse. The absence of a critical decrease in permeability with decrease in paleo-porosity content, indicates that pore reduction is the not the only mechanism responsible for permeability reduction during in situ calcite precipitation. Note that k tot and ϕ tot refers to the permeability and porosity of the maximum paleo-porosity found in the subvolumes. Also note that the last three steps (increased connected porosity) have similar porosities and equal permeabilities, so they plot on top of each other through relative time would take longer to occur (10 steps to reconnect 95% surface connected porosity) (Fig. 6b). This indicates that only new interfaces between fracture-fracture or fracture-amygdale open up during precipitation, indicating that those could have been preferential sites for precipitation and that, overall, they would have contributed more slowly to decrease porosity connectivity, as evidenced by Fig. 10a. This implies a decrease in dissolution and CO 2 mineralization rates along those preferential areas due to the lower reactive areas available. Our simulated precipitation event also shows that precipitation is slower along the X direction than in the other two directions (Figs. 6c, 10b). Looking closely at the XCT data provides evidence that this is caused by pore closing and snappingoff over time (Fig. 10b).
In all data, it is possible to note that the surface connected porosity trend yields a tipping point (Figs. 3, 4, 6a). Porosity-permeability trends have been found to describe geological processes (Bernabé et al. 2003;Van der Land et al. 2013) and tipping points become critical to assess one process versus another. During CO 2 -water-rock interactions, two main Fig. 10 Permeability-porosity evolution through relative time during the back-stripping method at variable precipitation rate on the chosen subvolume, for each direction. The porosity fraction of the surface connected porosity is defined here as the ratio between ϕ c (connected porosity) and the total paleo-porosity in the subvolume (ϕ max paleo ). The permeability is calculated as the ratio between the permeability of the connected porosity (k c ) and the total permeability (k tot ) processes drive changes in porosity and permeability, and related pore connectivity. On the one hand, dissolution of the silicate minerals in basalts can increase the porosity and release the divalent cations necessary for mineral carbonation (Kanakija et al. 2017;Sigurdur Reynir et al. 2010). On the other hand, mineral carbonation requires CO 2 to combine with the divalent cations to form carbonate minerals and therefore reduces the porosity of the rock (Kanakija et al. 2017;Sigurdur Reynir et al. 2010). Precipitation may be affected by various mechanisms and factors. It is well known that temperature, pressure, pH, saturation indexes and the presence of impurities affect reaction rates and therefore precipitation rates (Ahkami et al. 2020;Noiriel et al. 2016). At the microscale, factors such as particle Fig. 11. 3D renderings of the evolution of the surface connected porosity along the X direction in the subvolume used for the permeability calculations during the back-stripping process method 2 (a) and the simulated precipitation event (b). a Note how pores are opening and being dissolved (indicated by the red arrows) immediately after one back-stripping step through relative time, causing the permeability to increase (cfr. Figure 9a). b Similarly, note how there are pore throats drastically closing (indicated by the red arrows) after one precipitation step through relative time, causing the permeability to drop, which will eventually cause pores to snap-off through time. The side of the subvolume is 5.03 mm wide. Porosity fraction is reported as the ratio between the connected porosity in the subvolume (ϕ c ), and the total paleoporosity in the subvolume (ϕ max paleo ) shape and size may also impact dissolution and precipitation rates (Pluymakers and Spiers 2015). It has been suggested that, in addition to secondary mineral precipitation, which may alter fluid pathways and pore throats hence affecting overall precipitation, clay swelling may play a significant role in controlling precipitation, as, while absorbing water, their mineral structure increases as well as the volume, hence modifying the porosity and pore connectivity, and lastly permeability. However, this depends also on the porosity content and type of matrix grain (Aksu et al. 2015). The presence of clays was not found in the geological formation from which the sample comes from (MacDonald et al. 2019). Pressuresolution creep can also play a major role in controlling and driving precipitation across porous media (Gratier et al. 2013). Changes in porosity, permeability and pore connectivity will therefore depend on which one of these processes is predominant. Porosity and permeability evolution during CO 2 sequestration in basalts can also depend on the nature of the mineral filling the amygdales. According to Aksu et al. (2015) and Franzson et al. (2008), the amount of clay infilling the amygdales is inversely proportional to the porosity/permeability. Aksu et al. (2015) observed up to 40% decrease in permeability in swelling clay coated samples as compared to non-clay coated samples. Pore clogging by precipitation products and alteration minerals was also observed in experimental studies (e.g. Oelkers and Sigurdur Reynir 2001) and in similar natural outcrops (e.g. Alfredsson et al. 2013;Kelemen et al. 2019). The presence of roughness on the reactive surface areas as well as the shape of the pores may also play an important role in determining dissolution versus CO 2 mineralisation processes (Kanakiya et al. 2017). The values of permeability that we calculated along the X and Y directions in the first back-stripping method stay constant with precipitation ad porosity reduction, and no tipping point is observed. The magnitude of absolute permeability along these directions ranges between ~ 200 and 1000 mD (10 -12 and 10 -13 m 2 ), values which are in agreement with published values (Saar and Manga 1999) for vesicular basalts, and bulk permeabilities in deep sea basalts (Callow et al. 2018;Fisher 1998;Goldberg et al. 2008). However, results from permeability analyses using the second method of back-stripping highlight that permeability increases as soon as pores become open (Fig. 10b); similarly, permeability drops in the simulated precipitation event as soon as pores are being closed and eventually snapped off, to reach a value of permeability that is determined by the fracture network (Fig. 10b). It follows that permeability reduction during CO 2 mineral trapping is highly dependent on the precipitation mechanisms and rate, as well as on the present of directional micro fracture networks. Callow et al. (2018) results showed a two order of magnitude difference between the permeabilities calculated with the hydromechanical tests and the permeability derived from the XCT data. Because a 10 voxels pore reduction did not lead to a drop in the measured permeability from the XCT data to explain the results of the hydromechanical tests, they concluded that other factors, in addition to pore diameter reduction/pore collapse, cause a significant decrease in permeability. While we agree that other factors may play a role in determining porosity-permeability trend, we speculate that their sample might have also been controlled by fracture networks present at the sub-voxel scale, which may have been missed by the XCT imaging and that may be responsible for the difference in permeability in the hydromechanical tests. These results prompt further studies to determine CO 2 mineral trapping mechanisms in amygdaloidal basalts as analogues for CO 2 injections in basalt formations. The relatively simple back-stripping and simulated precipitation here can be replaced by more complex processes if necessary. It is possible to define the location and rate of any precipitation or dissolution process from the flow or pressure field generated by the permeability simulations, or to use more complex geometric dependencies based on the observed porosity and/ or fracture networks.
Lastly, limitations in XCT resolution may have missed the presence of microcracks or micropore network that could have contributed to the overall pore connectivity and permeability results. However, as demonstrated by Peng et al. (2014) this is not always the case. In their study, they found that sub-voxel porosity does not affect permeability in a Berea sandstone. Arguably, in XCT data imaging there is always a trade-off between representative sample size and optimal voxel size. Because of the highly heterogenous nature of these basalts, a smaller sample would have resulted in a less representative volume of the overall geological formation. Multi-scale techniques are also encouraged to address the problem of limitations in XCT imaging and the possible presence of microcracks and pores which could have affected permeability evolution.

CO 2 Sequestration Potential of Basalt
Our values of calculated specific surface area calculated in the first back-stripping method are considerably lower than what found by other studies on vesicular and young basalts, but they are overall of the same order of magnitude as the one found at the CarbFix site in Iceland reflecting basalt heterogeneity (Callow et al. 2018;Kanakiya et al. 2017;Saar and Manga 1999). We estimated the CO 2 potential in this sample of amygdaloidal basalt and we found that we can sequestrate 0.08 tonnes of CO 2 per tonne of basalt. Callow et al. (2018) reported a mean value of 0.33 Gt of CO 2 in 1.8 km 3 of basalt formation. Converting these values using the dry density of their samples (2.4 g/cm 3 ) resulted in 0.069 tonnes of CO 2 per tons of basalt, at the CarbFix site, which is the same order of magnitude for our basalt formation. However, our values of connected porosities are considerably higher than the ones reported for the CarbFix site, thus testifying basalts heterogeneity.

Conclusions
In this study, we used XCT to quantify how natural CO 2 mineral trapping occurred in an amygdaloidal basalt, as an analogue to future engineered CO 2 sequestration in basalts, to provide further insights into the evolution of porosity during CO 2 infiltration. We used novel back-stripping imaging methods to understand how porosity and pore connectivity may have evolved through time during calcite mineral precipitation in this sample, and to gain insights into CO 2 sequestration potential of basalts. Our back-stripping simulations have shown that when precipitation instantaneously fills all pores, porosity is disconnected quite sharply, in only 4 time-steps. However, when simulating precipitation at a constant rate, porosity is disconnected more slowly, particularly along the X direction, along which the fracture network develops. In this scenario, permeability simulations show that permeability increases within one time step as soon as the pores are open (during the backstripping), and drops, also within one time step, as soon as pores close and snap-off (during the precipitation event), to reach a value controlled by the fracture network. These findings indicate that during the final stages of precipitation permeability is controlled by the dominant fractures in the system, and they highlight the importance of further studies into the CO 2 potential of basalt formations.