Three-Dimensional Imaging and Quantiﬁcation of Gas Storativity in Nanoporous Media via X-rays Computed Tomography

: This study provides the engineering science underpinnings for improved characterization and quantiﬁcation of the interplay of gases with kerogen and minerals in shale. Natural nanoporous media such as shale (i.e., mudstone) often present with low permeability and dual porosity, making them difﬁcult to characterize given the complex structural and chemical features across multiple scales. These structures give nanoporous solids a large surface area for gas to sorb. In oil and gas applications, full understanding of these media and their sorption characteristics are critical for evaluating gas reserves, ﬂow, and storage for enhanced recovery and CO 2 sequestration potential. Other applications include CO 2 capture from industrial plants, hydrogen storage on sorbent surfaces, and heterogeneous catalysis in ammonia synthesis. Therefore, high-resolution experimental procedures are demanded to better understand the gas–solid behavior. In this study, CT imaging was applied on the sub-millimeter scale to shale samples (Eagle Ford and Wolfcamp) to improve quantitative agreement between CT-derived and pulse decay (mass balance) derived results. Improved CT imaging formulations are presented that better match mass balance results, highlighting the signiﬁcance of gas sorption in complex nanoporous media. The proposed CT routine implemented on the Eagle Ford sample demonstrated a 17% error reduction (22% to 5%) when compared to the conventional CT procedure. These observations are consistent in the Wolfcamp sample, emphasizing the reliability of this technique for broader implementation of digital adsorption studies in nanoporous geomaterials. the porosity and storativity results for the Eagle Ford and Wolfcamp samples obtained via CT. Comparisons are drawn against the mass balance results as well as results obtained using the conventional CT method. method (pulse decay) and the imaging approach that takes into account Kr density more accurately. The proposed CT routine and conventional CT method yield storativity values within 5% and 22%, respectively, of the mass balance derived results The CT curves are corrected to account for an additional 0.4 cm 3 of pore volume captured by the pulse decay material balance but not visible to the CT scanner.


Introduction
Nanoporous materials are ubiquitous throughout industrial applications, such as separations and catalysis, and in natural settings including energy resources. In unconventional shale resources, pore sizes may be on the order of a nm but properties vary over greater length scales (micrometer to millimeter). Characterization of static and dynamic properties over multiple spatial scales is frustrated by the inherent trade off in imaging between image resolution and sample size. X-rays computed tomography (CT) has the potential to provide quantitative details of in situ properties, such as porosity (void fraction), with fine spatial resolution over macroscopic length scales; however, processing of CT images to obtain such details is frustrated by gas sorption. Improved methods of laboratory scale characterization of shale samples with 1 mm distance between each slice. The large X-ray energy of 140 kV was chosen given that metals and rocks are difficult to penetrate [25].
A CT number is computed by normalizing the measured attenuation coefficient, µ x , so that water and air have CT numbers of 0 and −1000 in Hounsfield Units (HU), respectively, where CT is dimensionless HU (0.1% change in density with respect to calibration density scale), µ w is the linear attenuation coefficient of water, and µ a is the linear attenuation coefficient of air. The linear attenuation coefficient of a material depends on the incident X-ray intensity, the remaining detected intensity, and the thickness of the object [26,27]. X-ray sources are limited to large energies (140 kV) to minimize beam hardening that affects the linear attenuation coefficient as the incident beam intensity depends on the voltage, current, and scan time [28]. Beam hardening is an error that results from preferential attenuation of less energetic X-rays at sample interfaces. The result is a shift in detected attenuation to higher (harder) values. Typically, beam hardening appears as an apparent ring of large density material around a core image due to the abnormally large CT number measured. Therefore, 5% of the core periphery is cropped typically before any CT analysis to further mitigate this issue. Cylindrical geometry is used to avoid x-shaped artifacts [25].
Typically, water is used as a penetrant when measuring porosity using the fluid substitution method, but the extreme difficulty of injecting water through a shale sample due to its small permeability suggests the use of smaller-viscosity gases that attenuate X-rays. Kr attenuates X-rays strongly in the energy range typical of clinical CT scanners, compared to air or vacuum, and is relatively easy to use to saturate shale samples. The procedure starts by scanning a dry sample to obtain the voxel by voxel CT number of the vacuum-evacuated rock, CT r,ar . Then, the sample is saturated with Kr and scanned once again to get the voxel by voxel CT number of the Kr-saturated rock, CT r,Kr . It is necessary however, to measure the CT number of gases saturating the pore network at different pressures because the CT number is typically a function of fluid density.
Shale apparent porosity is obtained conventionally as [25]: where CT Kr is the CT number of pure Kr and CT ar is the CT number of air. Alnoaimi (2016) [29] [30] was used to construct Figure 1, that shows the pure Kr CT number versus pressure at 20 • C . The CT number of pure Kr, CT Kr , is usually measured by pressurizing Kr gas in a hollow aluminum core holder while recording the corresponding CT numbers. Pure Kr CT numbers acquired using this technique account for the free gas phase density only. This may lead to inaccurate values of porosity when Equation (2) is applied to adsorptive samples (e.g., high kerogen content) as the in situ Kr CT number measured captures both the free and adsorbed phases. Therefore, a new methodology is needed to capture the pure Kr CT number in an adsorptive rock that accounts for the effective attenuation from the adsorbed and free gas phases.
Converting Kr gas pressure values at 20 • C, as shown in Figure 1(left), to their corresponding density values, Kr density can be related to pure Kr CT number as shown in Figure 1(right). The curve depicted in Figure 1 is then translated to compute the density change associated with the measured change in the CT number after Kr saturation of the rock, as shown in Figure 2. The slope, β, of 5.9 HUm 3 kg −1 , as shown in Figure 2, is used to capture the in situ effective (free and adsorbed) density of pure Kr at a particular pore pressure given the assumption that the density-CT correlation remains linear and constant with pressure. This assumption is validated given that µ x is approximated as [31] where Z e is the effective atomic number, K c is a constant that depends on the electron shell, ρ is the fluid density, A r is the atomic weight, and E p is the energy of the incident photon beam. Notice that the linear attenuation coefficient is most sensitive to the effective atomic number and the energy of the incident photon beam, and it is linearly proportional to fluid density. Clearly, β depends on temperature. The voxel-by-voxel Kr density, ρ Kr−v , is computed using Equation (4)  Knowing the CT scanner voxel dimensions, the mass of Kr in each voxel, m Kr−v , and the total mass of Kr in the imaged rock, m Kr−t , are computed as shown below: where V v is the voxel volume, m Kr−s is Kr mass in each imaged slice (cross section), n vs is the number of voxels per slice based on the sample's diameter, and n s is the number of acquired slices based on the sample's length. The voxels cropped to eliminate artifacts are accounted for during these computations as their exclusion result in a significant loss in accuracy of mass. After reaching pressure equilibrium, the computed total Kr mass is converted to apparent Kr pore volume, V Kr−app , using the molecular weight of Kr, MW Kr , and the real gas law: where Z is the real gas compressibility factor, T is the temperature, P p is the pore pressure, and R is the universal gas constant. Finally, the apparent Kr porosity, φ Kr−app , computed via CT is determined using where n v is the total number of voxels in the image. We use the imaged rock dimensions instead of the true rock dimensions as a small portion of the rock periphery is cropped during analysis of porosity to mitigate beam-hardening artifacts. A more common method to quantify and compare rock storativity, S R , is by presenting it in standard cubic feet (SCF) of gas stored in a ton of rock. To achieve this, the computed Kr apparent pore volume is expanded to atmospheric pressure, P atm , and divided by the imaged rock mass, , as shown in Equation (10), where ρ B is the rock's bulk density.

Sample Preparation
Experiments were conducted on an Eagle Ford core and a Wolfcamp core. Both cores are 2.54 cm in diameter and the lengths of the Eagle Ford and Wolfcamp samples are 4.24 and 6.35 cm, respectively. The dry weights and depths of the Eagle Ford and Wolfcamp samples are 49.8 g and 11,184 ft and 76.3 g and 9708 ft, respectively. The mineralogy of both samples is captured in Table 1.
Core samples were sleeved between two 316 L grade stainless-steel end caps with 3-4 Teflon heat-shrink tube layers. Aluminum foil was applied in between the Teflon tube layers to minimize gas diffusion between the pore space and the confining fluid. Deionized water was used as the confining fluid by injecting it in between the Teflon tube and the outer aluminum body of the core holder using an ISCO pump. The aluminum body allows for sufficient CT transparency to image the in situ rock and saturating gas. Core preparations for drying and sleeving are identical to the procedure detailed by Elkady and Kovscek (2020) [24].

Experimental Procedure
We characterized the rock samples porosity and storativity using the pulse decay technique detailed in Elkady and Kovscek (2020) [32]. The accuracy of the pulse decay storativity measurements depends on the pressure transducers' precision, the accuracy of the pre-measurement of the upstream and downstream reservoir volumes, and the pore volume to reservoir volume ratio. The pressure transducers' precision is within 1 psi and the accuracy of the upstream and downstream reservoir volumes are within 2%. The reservoir volumes were chosen to be comparable to the sample's pore volume to achieve reasonable pressure sensitivity. Reservoir volumes that are too large result in minimal pressure change during a pressure pulse experiment, while reservoir volumes that are too small introduce very little gas into the rock resulting in small pressure buildups. The experimental setup shown in Figure 3 was designed to obtain the mass balance derived storativity. Therefore, the setup was placed on a CT table for image acquisition to compare CT-derived storativity results with the pulse decay results. Once the core was sleeved and vacuum evacuated, a full CT scan of the rock was acquired. The upstream reservoir was then pressurized by opening the first (leftmost) valve to a Kr source while keeping the second valve closed. The first valve was then closed and the upstream pressure was monitored until the reading stabilizes. The second valve was then opened to allow Kr to pass through the core to the downstream end. The last (right-most) valve was always kept closed. The pressure upstream and downstream of the sample was recorded using pre-calibrated Precise Sensors transducers rated up to 2000 psig. A second full scan was obtained after the pore pressure remains constant with time. The equilibrium condition was identified when the upstream and downstream pressures were equal and had not changed for 6-12 h.
The confining pressure was initially set to 500 psia above the average pressure of the charged upstream and downstream reservoirs. Once a better estimate of the equilibrium pressure became apparent, the confining pressure was adjusted again during the saturation process to ensure the 500 psia effective stress at equilibrium. Successive pulses were initiated by repeating the same steps while cautiously maintaining a 500 psia effective stress. Here, instead of starting at vacuum conditions, each additional pulse added more Kr moles to the closed system, thus increasing equilibrium pore pressure, as illustrated in Figure 4 for the Eagle Ford sample. Kr successive pressure pulses were performed under the CT scanner so that the vacuum evacuated rock image and Kr-saturated rock images were acquired and used at the end of each pulse to compute the CT-derived porosity and storativity values. The CT-derived results were compared against the results obtained via mass balance. Note that the procedure and results summarized in Figure 4 are also useful to measure the matrix to fracture transfer function [33].

Results
This section presents the porosity and storativity results for the Eagle Ford and Wolfcamp samples obtained via CT. Comparisons are drawn against the mass balance results as well as results obtained using the conventional CT method.

Krypton Storativity via CT
We compare the pulse decay results with the proposed CT methodology for computing Kr storage capacity. This is essential as CT computed Kr apparent porosity is often referenced against helium (He) pulse decay porosity to highlight adsorption. However, there is a lack of testing to prove that both the pulse decay and CT methods are capable of producing similar results for Kr.
We initially used Equation (2) to compute Kr storativity via CT and realized that the both the shape and values of the curve do not match the mass balance results. Fundamentally, Equation (2) normalizes the subtracted rock images with the CT number of Kr in free phase that does not take into account the denser adsorbed phase. The proposed CT method, however, does produce a storativity curve that better matches the material balance results (within 5%) for the Eagle Ford sample, as shown in Figure 5. The proposed method captures the effect of adsorption by assuming a linear relationship between Kr density and CT number that results in effective CT numbers that lay on the curve shown in Figure 2. Therefore, an effective Kr density in each voxel is obtained directly. Knowing the voxel dimensions, the Kr mass in each voxel is visualized via CT to show areas with low and high storativity, as seen in Figure 6.
Both CT methods are corrected for 0.4 cm 3 dead space captured by the material balance method but is invisible to the CT. This includes small volumes occurring from the slightly uneven rock surface that creates small dead pockets between the rock and end caps as well as any minuscule crevices between the rock and sleeving material. This is especially relevant given that 5% of the rock periphery is cropped out before CT analysis due to artifacts. An uneven rock surface produced during rock cutting was confirmed via CT showing a small pocket between one of the end caps and the rock (see Figure 7). In future experiments, this correction can be eliminated or reduced by precise cutting, tight sleeving of the core, and upstream and downstream reservoir volume measurements after rock sleeving. Note the excellent agreement between the mass balance method (pulse decay) and the imaging approach that takes into account Kr density more accurately. The proposed CT routine and conventional CT method yield storativity values within 5% and 22%, respectively, of the mass balance derived results The CT curves are corrected to account for an additional 0.4 cm 3 of pore volume captured by the pulse decay material balance but not visible to the CT scanner.  The lighter (2000 HU range) areas reflect etching in the end caps that allows for a more distributed inflow of gas into the core. The changing CT numbers or density in these etched areas reflect the uneven contact surface between the rock and end cap. The low CT number (less dense) blue zone suggests that a void exists between end cap and core in that region.
The Wolfcamp core Kr CT porosity results are compared in Table 2 for three different pressure pulses. The Wolfcamp sample permeability was 2-3 times smaller than the Eagle Ford sample, therefore pressure pulses took longer to equilibrate. Given that the pulse decay experiments were conducted under CT, only three pressure pulses were attainable in a reasonable time frame on this multiuser scanner. Nevertheless, each data point (whether Eagle Ford or Wolfcamp) was used independently to test the new CT processing approach.
Apparent porosity is used as an alternative to storativity to illustrate that these quantities may be used interchangeably to quantify the storage capacity of the rock. The results presented in Table 2 show that the new CT formulation is in better agreement with the pulse decay derived Kr porosity. Moreover, the decreasing trend captured by the pulse decay method is also reflected by the new CT method that is a result of the increasing free gas phase contribution to total porosity. The conventional method predicts a smaller Kr porosity that is even less than the pulse decay measured He porosity (11.0%) [32]. Additionally, Table 2 shows an increasing Kr porosity with increasing pore pressure. This contradicts observations made from the pulse decay results. The increasing porosity captured by the conventional method is due to the increasing contribution of the free phase gas to the total porosity as pore pressure increases. As the conventional method considers the free phase gas only, the increased free phase contribution with pressure is reflected as greater porosity. Figure 8 depicts 3D CT images showing the apparent porosity distribution using both CT methods for the Wolfcamp sample. The proposed method suggests that adsorption results in local Kr apparent porosity close to 25%, whereas other areas have almost no gas storativity.
Observing Table 2, the new CT formulation Kr porosity is still smaller than the pulse decay results. The difference is about 20% for the measurement of 556 psia. A 0.5 cm 3 void space in-between the core, sleeving Teflon, and ends caps would account for this discrepancy. Nonetheless, the new CT formulation agrees substantially with the pulse decay Kr porosity compared to the conventional method. In most, if not all, cases, the new CT method will produce slightly lower porosity values compared to the pulse decay method. One could argue that the new CT method is more accurate as the pulse decay method slightly overestimates porosity due to very small dead spaces around the core that cannot be avoided and are virtually unquantifiable. Further experimentation and comparisons are needed before a definitive conclusion can be made with regard to the method that is more accurate. In conclusion, a better CT processing method is proposed that is more consistent with pulse decay results.

Summary
Both pulse decay and X-ray CT imaging were used independently to establish consistency between storativity results derived from each method. New image processing routines for CT data were developed that better match mass balance derived porosity and storativity results compared to conventional CT methods. This consistency between both techniques was not previously attainable in shale using the conventional CT method.
Measurement of in situ porosity via CT scanning is improved by incorporating the contribution from both the free and adsorbed gas phases in our formulation. In shale characterization work found in the literature, little cross-examination is made between He pulse decay porosity and storage capacity determined using Kr in shale because of the difficulty of correcting for the sorption of gases. This motivated initiating comparisons between both characterization methods to establish consistency. The newly proposed methodology was tested on two shale samples (Eagle Ford and Wolfcamp) at various equilibrium pore pressures. In all cases tested, the proposed method had better agreement with the pulse decay results than the conventional method. After dead volume correction, the proposed CT routine and conventional CT method yield storativity values within 5% and 22%, respectively, of the mass balance derived results. Moreover, the reduction in apparent porosity with increasing pore pressure was also consistent between pulse decay and the proposed CT method. The reduction in Kr apparent porosity with increasing pore pressure is a consequence of the increasing contribution of the free gas phase. On the other hand, the conventional method occasionally exhibited unrealistic apparent porosity numbers where Kr CT porosity was lower than the He pulse decay porosity. We strongly encourage the use of the newly proposed method for future CT porosity computations in rocks that strongly adsorb gas, such as shale.
Further testing can be made to compare results for strongly adsorbing (shale) and weakly adsorbing rocks (sandstones and carbonates) to test the performance of both CT methods. Additionally, upstream and downstream reservoir volumes should be measured while the core is sleeved into the setup to account for void spaces that may lead potentially to slight discrepancies between CT and mass balance methods.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: