- Email: [email protected]
Improved Vinegar & Wellington calibration for estimation of fluid saturation and porosity from CT images for a core flooding test under geologic carbon storage conditions
Micron 124 (2019) 102703
Contents lists available at ScienceDirect
Micron journal homepage: www.elsevier.com/locate/micron
Improved Vinegar & Wellington calibration for estimation of fluid saturation and porosity from CT images for a core flooding test under geologic carbon storage conditions
T
⁎
Xiuxiu Miaoa,b,1, Yan Wanga,1, Liwei Zhanga,b, , Ning Weia,b, Xiaochun Lia,b a State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan, 430071, Hubei Province, China b University of Chinese Academy of Sciences, Beijing, 100049, China
A R T I C LE I N FO
A B S T R A C T
Keywords: X-ray CT Image processing Two-phase displacement flooding Saturation mapping Geologic carbon sequestration
X-ray computed tomography (CT) of fluid flow in formation rocks is an important characterization technique in geologic carbon sequestration research to provide insight into the migration and capillary trapping of CO2 under reservoir conditions. An improved calibration method adapted from traditional Vinegar & Wellington calibration is proposed to map the 3D pore and fluid distributions from the CT images of CO2/brine displacement flooding. Similar to Vinegar & Wellington calibration, the proposed method adopts the linear scaling law of CT number transformation to mass density. However, different from Vinegar & Wellington calibration that uses a 100% brine-saturated core image and a 100% CO2-saturated core image as references to calculate CO2 and brine saturations at all time steps, the proposed method uses the CT numbers of CO2 and brine to calculate the incremental of CO2 and brine saturations from time step i to time step i +1. The method is intended for cases in which the two 100% brine saturation and 100% CO2 saturation images can not be successfully obtained. Overall, the improved calibration proposed by this study presents more reasonable results of CO2 and brine distribution in a Berea sandstone core, as compared to traditional Vinegar & Wellington calibration. The reconstructed porosity image agrees with the laminated structure of the Berea sandstone core, and the average porosity evaluated over the entire core (0.176) is comparable to the physical porosity (0.165). Furthermore, the reconstructed saturation images using the improved calibration reveal a flat piston-like flooding front from a homogeneous longitudinal-section of the 3D orthogonal view and preferential fingerings from another nonhomogeneous longitudinal-section, which are not present in the reconstructed saturation images using traditional Vinegar & Wellington calibration. Concerns and causes with respect to the uncertainty of linear CT number calibration are also explained, and approaches to alleviate the uncertainty are suggested.
1. Introduction Successful geologic carbon sequestration requires long-term trapping of CO2 in deep formations without significant CO2 leakage back to the atmosphere. The residual trapping mechanism that depends on the capillary equilibrium between brine and injected CO2 in the porous rock formation is believed to be crucial to determine the storage capacity of CO2 and the long-term fate of injected CO2 (U.S. Geological Survey Geologic Carbon Dioxide Storage Resources Assessment Team, 2013). Also, it is necessary to understand migration behavior of both CO2 and water in subsurface CO2 storage formations. In laboratory studies, core samples with sizes in centimeter range are used in CO2 and
brine flooding experiments under simulated deep reservoir conditions, so as to obtain valuable information about the physical mechanisms that control pressure-driven CO2 and brine migration in porous rocks. X-ray CT technology combined with quantitative image analysis has made visualization of the spatial distribution of CO2 and brine in heterogeneous rock media possible. Moreover, the movement and flow characteristics of CO2 and brine can be revealed from CT image analysis. Mapping saturation levels and porosity within core samples by Xray CT has been practiced since the 1980s. Some of the earliest works were conducted in 1987 by Vinegar and Wellington for two-phase (oil/ water, Wellington and Vinegar, 1987) and three-phase (oil/water/air, Vinegar and Wellington, 1987) displacement flooding. In recent studies,
⁎
Corresponding author. E-mail address: [email protected] (L. Zhang). 1 Co-1st author: Dr. Xiuxiu Miao and Prof. Yan Wang have made equal contribution to this paper. https://doi.org/10.1016/j.micron.2019.102703 Received 7 March 2019; Received in revised form 7 June 2019; Accepted 7 June 2019 Available online 29 June 2019 0968-4328/ © 2019 Elsevier Ltd. All rights reserved.
Micron 124 (2019) 102703
X. Miao, et al.
Zhang et al. (2014) imaged the CO2 and brine flow during drainage and imbibition in a Berea sandstone core using a medical X-ray CT scanner; Menke et al. (2018) used 4D X-ray micro‒computed tomography (μCT) to study the impact of flow heterogeneity on CO2-induced dissolution of limestone cores at both the pore and core scales under reservoir conditions; Cao et al. (2013) quantified changes in wellbore cement integrity by X-ray μCT analysis of pore volume and structure changes of wellbore cement during continuous flooding of CO2-saturated brine. With CT imaging of CO2 flooding of core samples saturated with water, studies have shown that the heterogeneity of pore distribution significantly influences brine and CO2 distributions in rocks (Perrin and Benson, 2010; Shi et al., 2011). The feasibility of using core CT images to compute porosity and fluid saturation, as well as for many other quantification applications, is based on a linear relationship between CT image grayscale intensity (CT number) and material density. For medical CT scans, CT number is expressed in Hounsfield Units (HU) as in Eq. (1)(Akin and Kovscek, 2003; Bryant et al., 2012):
CT = 1000
μ − μW μW
SiA =
(1)
(2)
where ρ is density, σ(E) is the Klein-Nishina coefficient that is presumed “nearly energy-independent” for CT X-ray with energy higher than 100 keV; b is a constant that is approximately 9.8 × 10−24; Σ(fkZk3.8) represents the effective atomic number for a mixture of atomic species with fk and Zk being the fraction of electron and the atomic number of the kth species, respectively; E denotes the X-ray beam energy. As shown in Eq. (1), the CT number of vacuum (zero density) is -1000 HU and the CT number of water is 0 HU by definition. For any other substance, the CT number is energy-dependent. Moreover, for a given substance at a given energy, the CT number is linearly correlated with density. This is the basis for Vinegar & Wellington’s linear calibration for porosity and fluid saturation from CT images, which is achieved by replacing the density of all involved phases with their attenuation coefficients or their CT numbers in the mass and density relationship. For a pixel saturated with fluids A and B in a core image, the relationship is re-written as:
2. Core-flooding CT experiment description A core flooding test that displaced brine in core samples with CO2 was conducted at National Energy Technology Laboratory (NETL), USA using a medical CT scanner with a Berea sandstone core. The Berea core was approximately 5 cm in diameter and 16 cm in length, and the porosity was about 0.165. The temperature of the flooding test was maintained at 67 °C, a confining pressure of 25.6 MPa was applied to the core, and the outlet fluid pressure was maintained at 21.4 MPa. Prior to the flooding cycles, the air-saturated core was scanned as the first reference. After that, supercritical CO2 was injected till the core was fully saturated with CO2, and then the core was scanned as the second reference. Following the supercritical CO2 saturation, the displacement flooding was initiated with a brine injection cycle, and alternated between a CO2 cycle and a brine cycle at constant flow rates of 3.43 ml/min, 1.72 ml/min and 5.13 ml/min, subsequently. Within each flooding cycle, several flooding stages with different amounts of injected fluid were scanned. Each flooding cycle stopped at a stage when no more displacement was detected. In-situ scanning of the core at different CO2/brine saturation stages was simultaneously implemented. To enhance brine-to-CO2 contrast, the brine solution was conditioned with 5 wt% KI to enhance the attenuation coefficient of brine phase, taking advantage of the Kα absorption of X-ray beam by element iodine. The medical CT machine, flooding schematic and the Berea core are shown in Fig. 1. The scanning dataset obtained at different flooding
j μCore = ϕAj μA + ϕBj μ B + ϕRj μRj = ϕ jSAj μA + ϕ j (1 − SAj ) μ B + (1 − ϕ j ) μRj or j CTCore = ϕ jSAj CTA + ϕ j (1 − SAj ) CTB + (1 − ϕ j ) CTRj
(3) where, the superscript j denotes the pixel-wise value; ΦA, ΦB and ΦR are the volume fractions of fluid A, fluid B and rock matrix, respectively; SA is the saturation of fluid A, and the saturation of fluid B equals to 1 minus the saturation of A (i.e. SB = 1 - SA); the attenuation coefficients (CT numbers) of fluid A, B and the rock matrix are denoted as μA, μB and μR (CTA, CTB and CTR), respectively. By imaging a 100% A-saturated core and a 100% B-saturated core, the 3D porosity image can be solved by Eq. (4), and then the saturation image of fluid A or B can be calculated using Eq. (5), given a core image with pores occupied by both A and B.
ϕ=
Bsat Asat CTCore − CTCore CTB − CTA
(5)
The aforementioned calibration (Eq. 4–5), which is referred to as Vinegar & Wellington calibration, has been widely adopted for core porosity and fluid saturation estimation in various CT scanning cases (Akin and Kovscek, 2003; Alemu et al., 2011; Crandall and Bromhal, 2013, 2014; Crandall et al., 2014; Ning et al., 2014; Niu et al., 2015; Ott et al., 2012; Perrin and Benson, 2010; Pini et al., 2012; Zhang et al., 2014) involving different core types (e.g., sandstone, limestone, mudstone, shale, etc.), and different mixtures of miscible or immiscible fluids (e.g., CO2 and water, oil and water, CO2 and H2S, CO2 and N2, etc.) in displacement or co-injection flooding tests. In this study, Vinegar & Wellington calibration was initially employed to estimate porosity and saturations of CO2 and brine from medical CT images obtained from alternating CO2 and brine displacement flooding cycles under laboratory-simulated reservoir conditions. The calibration was based on CT images of Berea sandstone samples with pores filled with CO2 and brine in a CO2 flooding experiment. Step-by-step procedures of CT image analysis were presented. However, Vinegar & Wellington calibration gave unintelligible results of fluid saturation. A possible reason for Vinegar & Wellington calibration to give unsatisfactory results was that certain presumptions for Vinegar & Wellington calibration were not met in our experiment. Therefore, we modified Vinegar & Wellington calibration to analyze our images, and established comprehensive procedures for the reconstruction of porosity and saturation images in hope that our approach can be applied in future studies experiencing similar situations. Our approach uses the CT numbers of both CO2 and brine, and fluid content images of CO2 and brine at time step i to map the fluid content images of CO2 and brine at time step i + 1. The porosity image is obtained by combining the maximum fluid content at each pixel from all time steps. The saturation images of CO2 and brine are then derived by dividing the fluid content images with the porosity image. All raw CT images were processed with the use of ImageJ, a powerful Java-based image processing program (Smith et al., 2019). CO2 and brine saturation maps based on the calculation of the modified Vinegar & Wellington calibration were generated by ImageJ as well. The slice-average porosity, fluid saturation variation along the core and the average fluid saturation for the entire core at each flooding stage were computed and plotted using MATLAB.
where μ is the attenuation coefficient of substance for X-ray beam, and μW is the attenuation coefficient of water. The attenuation coefficient of substance is a function of electron density (or bulk density), effective atomic number, and beam energy. One formula suggested by Vinegar and Wellington (Vinegar and Wellington, 1987; Wellington and Vinegar, 1987) is:
⎡ ⎤ μ = ρ ⎢σ (E ) + b ∑ (fk Zk3.8) E 3.2⎥ k ⎣ ⎦
Bsat Asat − CT iCore CTCore CT iCore − CTCore or SiB = Bsat Asat Bsat Asat CTCore − CTCore CTCore − CTCore
(4) 2
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 1. Core flooding experiment: (a) medical CT scanner; (b) flooding schematic and reconstructed and original images of the core (Φ 5 cm × 16 cm).
slice-average grayscale intensity of the core is taken. Fig. 3 showed the average grayscale intensity at each CT slice for selected flooding stages (each note, i.e., Inj1-B1, represents a flooding stage), and the average grayscale intensities of each flooding stage were plotted along the axis of the core (the z-direction, representing the locations of CT slices from one end to the other end of the core). As shown in Fig. 3, the intensity fluctuates from slice to slice along the Z-direction. Continuous intensity increase (from curves Inj3-B1, B2 to B6) during brine invasion and gradual intensity decrease (from curves Inj4-C2, C3 to C6) during CO2 injection were noticed. Observation of these same curves also suggests that brine injection is much easier than CO2, as the change from Inj4-C6 to Inj5-B1 is much more significant compared to the change from Inj3B6 to Inj4-C2. Another important phenomenon is that all brine invading curves terminate somewhere close to Inj3-B6 (considered as the brinesaturated curve). However, the CO2 invading curves stop at a point much higher than Inj0-CO2 (taken as the CO2-saturated curve). This gap between the CO2-saturated curve and the CO2 invasion terminating curves could be related to residual brine. Further discussion on the displacement pattern is available in the Methodology of porosity and fluid saturation calculation section. In general, this intensity change is consistent with fluid content change. CO2 has higher CT number than air, therefore, when air is displaced by CO2, the average grayscale intensity increases. Similarly, as KI dosed brine has higher CT number than CO2, the core region intensity increases when brine displaces CO2, and decreases when CO2 displaces brine.
Table 1 CT scan dataset obtained at different flooding stages. Flooding cycle
Flow rate
Number of stages
Data label
air saturation CO2 saturation 1st brine flooding 1st CO2 flooding 2nd brine flooding 2nd CO2 flooding 3rd brine flooding 3rd CO2 flooding
/ / 3.43 ml/min
1 1 6 8 6 6 3 5
Inj0-Air Inj0-CO2 Inj1-B1 to Inj2-C1 to Inj3-B1 to Inj4-C1 to Inj5-B1 to Inj6-C1 to
1.72 ml/min 5.13 ml/min
B6 C8 B6 C6 B3 C5
stages is enclosed in Table 1. Each data label in Table 1 relates to a CT scan of a flooding stage. For instance, during the second flooding cycle of brine (Inj3), a total of 6 flooding stages were imaged, and the 3D image data obtained from the first to the sixth flooding stages within this flooding cycle were labelled from Inj3-B1 to Inj3-B6 sequentially. Description of the CT machine and flooding apparatus can be found in Crandall’s previous publications (Crandall and Bromhal, 2013, 2014; Crandall et al., 2014), and the flooding test protocol is also similar to those described in previous papers. 3. CT image description Porosity and saturation calculations were performed on reconstructed CT image sequences. An image sequence consisted of 337 contiguous 2D image slices at a size of 512-pixel by 512-pixel. The pixel resolution was 0.43 mm/pixel and the depth resolution was 0.5 mm/ slice. The 3D view of scanned object could be reconstructed from the image sequence by sequentially stacking contiguous 2D slices. Presented in Fig. 2 are the longitudinal-sectional and cross-sectional views of the reconstructed flooding apparatus and the core. Though the scanning resolution was not adequate to discriminate pore, rock, CO2 and brine in the core region, it was sufficient to reveal the artificial channel at the flooding entrance, different core textures between the flooding entrance and the main core-body. The components that made up the flooding apparatus, including rubber sleeve, confining fluid, and holder, were also revealed. Slices at the beginning and the end of an image sequence were typically covered with heavy noises, and must be excluded from porosity and saturation calculation. For in-situ CT imaging of porous media using constant scanning parameters, we assume that if pore space occupation at a given location does not change, then the “grayscale intensity” (an image term, the same as CT number) at this location remains the same in different “frames” (an image term, means CT images taken at different times, the same meaning as “time steps”). Therefore, any variation of grayscale intensity indicates changes of pore space occupation. In the case of flooding experiment, the intensity change is related to displacement of fluid. To further explore the core region information, a close look at the
4. Methodology of porosity and fluid saturation calculation 4.1. Image pre-processing To compute porosity and fluid saturation is to find a direct conversion of grayscale intensity change to the change of fluid content, porosity and fluid saturation are then derived from estimated fluid content. A difference image representing grayscale intensity change can be generated by subtracting 2 frames of in-situ CT images, for instance, a brine/CO2 flooding frame and an air-saturated frame, as shown in Fig. 4. However, due to long staging time and flooding pressure disturbance, slight drift of the sample occurred, which had to be addressed before image subtraction would work. Apart from sample drift, inevitable CT noises also needed to be evaluated and mitigated when necessary. Due to sample drift, several image frames are not aligned. 3D rigid registration was applied to align all image frames to the first frame (the air-saturated frame), so as to save the trouble of manual selection of problematic frames. Rigid registration is a linear transform of image without scaling; it is one of the most sophisticated practices to align different time-frame images obtained from the same source. We used the descriptor-based image registration plugin (Preibisch et al., 2010) 3
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 2. Anatomy of core flooding CT image, revealing the core, the sleeve, confining fluid, the holder and air: (a) longitudinalsection view; (b) slice with artificial channel from flooding entrance, and typical heavy-noise slice that may present at the front and the end of an image sequence; (c) transitional core slice where artificial channel ends; (d) slice from core-body region, showing texture different from flooding entrance.
Fig. 3. Change of average grayscale intensity along the core during flooding experiment (partial dataset; for full dataset see Fig. SI5).
window sizes were tested on the difference images to remove noises. It is believed that the radius 1 (pixel) Gaussian filter is adequate after comparing the Gaussian filtered images with different window sizes. The comparison is elaborated in Section 2 of the Supplementary information.
from Fiji (an ImageJ distribution) (Schneider et al., 2012) to implement the 3D rigid registration. Details on how the 3D rigid registration was implemented and how the registration improved image alignment are explained in Section 1 of the Supplementary information. Since all CT images are liable to noise, and all noise filtering algorithms will inevitably altering true intensity signals to some degree, we suggest a conservative noise-filtering practice, which is to apply filtering as rare and as mild as possible. In this study, no filtering was implemented before image registration, so that the difference signal is maximally preserved. After registration, Gaussian filters of different
4.2. Fluid saturation by vinegar & Wellington calibration According to Vinegar & Wellington linear calibration (Eq. 4–5), we could calculate core porosity and saturations of brine and CO2 in the 4
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 4. Application of the subtraction method to create difference image: (a) brine/CO2 invaded core image; (b) air saturated core image; (c) difference image.
core image obtained at an arbitrary stage i (SiCO2 and SiBrine) by the following equation: Brine-sat Air-sat − CTCore ⎧ ϕ = CTCore ⎪ BA CTBrine − CTAir ⎪ Brine-sat 2 -sat ⎪ − CTCO CTCore Core ϕBC = CTBrine − CTCO2 ⎨ ⎪ CO2 -sat Air-sat − CTCore CT ⎪ϕ = Core ⎪ CA CTCO2 − CTAir ⎩
calibration does not need the 100% CO2 and brine saturated images. If the two 100% saturation images for reference are correctly acquired, then ΔθiBrine = -ΔθiCO2 and θ0CO2 = φ, and this method is consistent with Vinegar & Wellington calibration.
i Brine-sat − CT Core CTCore ⎫ Brine-sat 2 -sat ⎪ − CTCO CTCore ⎪ Core CO2 -sat ⎬ i CT Core − CTCore ⎪ = Brine-sat 2 -sat ⎪ − CTCO CTCore Core ⎭
i = SCO 2
SiBrine
Inj0-CO2 Air-sat − CTCore CTCore ⎧ 0 ⎪ θCO2 = CTCO2 − CTAir ⎪ +1 = θiCO θiCO2 + ΔθiCO2 2 ⎨ ⎪ θiCO2 SiCO2 = ⎪ ϕ ⎩
(6)
where, φBA is the porosity obtained at the stage of brine displacing air; φBC is the porosity obtained at the stage of brine displacing CO2; φCA is the porosity obtained at the stage of CO2 displacing air. Theoretically φBA = φBC = φCA. The CT numbers of KI-dosed brine (CTBrine), CO2 (CT CO2) and air (CTAir) are 735, -400, -1000, respectively. In Fig. 5a, the porosity calculated using different linear calibration formula did not conform to each other, and the CO2 saturation significantly fluctuated along the core (Fig. 5b). The reason that the saturation values calibrated by Vinegar and Wellington method have significant variation is explained in the chessboard simulation section (Section 3 of the SI). Also, frame-average fluid saturation within one CO2 or brine injection cycle did not progressively increase or decrease as expected (Fig. 5c). For Vinegar & Wellington calibration to work, the assumptions are that there are two 100% saturation states established for brine and CO2, respectively. All effective pores in the single-phase referential saturation image must be occupied by only one fluid phase. However, in our case, those presumptions are not met. It is possible that in our experiment, not all effective pores were occupied by CO2 when the CO2 saturated image (Inj0-CO2) was obtained, and that the brine saturated image (Inj3-B6) was much likely acquired at a co-saturation state of both brine and CO2. The fact that no CO2 was detected at the flow outlet during a brine displacing CO2 stage does not ensure no residual CO2 in the brine-saturated core. The fact is also the case for CO2 displacing air or CO2 displacing brine scenarios. We simulated the core matrix with a 6-by-6 chessboard to demonstrate the possibility for the predicted values of porosity and saturation to go beyond the 0–1 limit using the Vinegar & Wellington calibration. The description of the chessboard simulation is provided in Section 3 of the Supplementary information.
0 =0 θBrine ⎧ i ⎪ θi + 1 = θi Brine Brine + ΔθBrine i ⎨ θ SiBrine = Brine ⎪ ϕ ⎩
(7)
The implementation of our calibration is described in this paragraph. Based on flooding sequence (Fig. 1b & Table 1), the initial brine and CO2 contents (θ0Brine and θ0CO2) were computed for the first CO2 displacing air frame (Inj0-CO2). Then incremental of brine or CO2 content at imaging stage i+1 against its previous stage i was evaluated from Eq. (8). By simultaneously comparing pixel j at stage i+1 with its previous stage i, and pixel j at the air-saturated stage, three fluid displacement patterns from stage i to stage i+1 at pixel j, including CO2brine displacement, brine-air displacement, and CO2-air displacement, were recognized and computed. Due to system noises, even air intensity fluctuated from frame to frame (Figure SI7), and the relative mild Gaussian filter could not completely remove this fluctuation (Figure SI7). Therefore, we introduced a tolerance (Tol = 0.05) to Eq. (8) to address the fluctuation. (8) The porosity image was calculated from the maximum fluid content at each pixel j from all frames, as expressed in Eq. (9). The formula meant that the estimation was based on a record of all pores once occupied by fluid during the whole displacement flooding process.
ϕ = max(θ j ) ∀i
(9)
5. Results and discussion Following the procedures of registration, filtering and reconstruction of fluid content images and the porosity image, the saturation images were then computed by dividing fluid content images with the porosity image. Presented in Fig. 6 are the porosity, brine and CO2 saturation images at one CO2 flooding stage (Inj4-C4) processed by the improved calibration described in this paper. The orthogonal view of the porosity image (Fig. 6a) revealed the intrinsic heterogeneity of the core. The longitudinal section was marked by a strip-like pattern that was parallel to the main flow direction, which agreed with the laminated structure of Berea sandstone (Fig. 1b). It was also obvious that most pixels concurrently hosted pore space and the rock material, though a few pixels had porosity values very close to 0 or very close to 1. Due to heterogeneous porosity distribution, brine and CO2 flow behaviors shown in different longitudinal-sections (Fig. 6b-c) varied significantly from each other. A flat piston-like flooding front was observed on the left
4.3. Description of the new approach to reconstruct porosity and fluid content images To handle the failure patterns when applying Vinegar & Wellington calibration to process our images, the authors suggested a new calibration method to reconstruct porosity and fluid content images. Given the fact that the 100% brine saturation frame and the 100% CO2 saturation frame could not be acquired in our flooding experiment, and there was possible air-phase interference, the authors suggested using an incremental calibration (Eq. 7) to estimate brine and CO2 saturation. The advantage of this new incremental calibration is that the 5
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 5. Porosity and fluid saturation converted from CT number by traditional Vinegar & Wellington calibration: (a) porosity; (b) CO2 saturation (partial dataset; for full dataset see Fig. SI6); (c) frame-average CO2 and brine saturation. Poor average CO2 saturation results and frame-average saturation results from direct linear calibration can be clearly seen.
longitudinal-sections in Fig. 6b-c, while preferential fingering dominated the right longitudinal-sections in Fig. 6b-c. The CO2 fingering was able to penetrate through the brine phase, and much less brine was displaced on the longitudinal-section with preferential fingering compared to the longitudinal-section with piston-like flooding front. Fig. 7 serves as a comparison between results of Vinegar & Wellington calibration and the calibration method proposed in this paper. The CO2 saturation profile for the frame (Inj6-C2) computed by our proposed method exhibited pixels showing clear and distinguishable flooding signals (Fig. 7b, saturation value is cut off to the range of 0 to 1). The actual saturation value distribution from -1.1 to 1 is shown in
Fig. 7c, demonstrating that more than 99% pixels have saturation values from 0 to 1. In contrast, the saturation images obtained by Vinegar & Wellington calibration were spoiled by extremes, and barely showed any distinguishable signals (Fig. 7a, saturation value is cut off to the range of −100 to 100). The actual saturation value distribution from −695 to 613 is shown in Fig. 7c, demonstrating that less than 45% pixels have saturation values from 0 to 1. These extreme values mainly originated from the residual CO2 pixels in the referential brine saturation image, described as the first failure pattern of Vinegar & Wellington calibration in the chessboard simulation (Section 3 of the Supplementary information). The slice-average saturation in Fig. 7d also 6
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 6. Orthogonal view of brine and CO2 saturation in a CO2 flooding stage (Inj4-C4) processed by the improved calibration: (a) the porosity image; (b) brine saturation image; (c) CO2 saturation image. In each image, two longitudinal sections and two cross sections are shown.
in the core at different flooding stages and injection rates. It should be noted that the images were not strictly acquired at the same interval of time, and the injection rate within each injection cycle (corresponding to a row of images) also differed from each other due to injection pressure variation. For instance, during the fifth flooding cycle of brine (Inj5-B1 to B3), brine swiftly took up most of the pores at the first injection (that was the reason why only 3 images were obtained), while at lower injection rate 6 images were acquired for the third flooding cycle of brine (Inj3-B1 to B6). Therefore, we can observe that the brine saturation on the Inj5-B1 frame was much higher than that on Inj3-B1. Nevertheless, the cyclic fluid intrusion/recession during the flooding cycle was clear. Barely any brine was injected at the Inj1-B1 flooding stage, and CO2 saturation was initially high, yet did not reach 100% saturation. Following the flooding cycles, the amount of CO2 in the core
shows larger variation produced by Vinegar & Wellington calibration, as compared to our method. Meanwhile, the slice-average porosity of the Berea core obtained from our calibration approach (about 0.176) was closer to the physical porosity (0.165) than the average porosities from VW method (the average porosities of VW-BA (brine displacing air), VW-BC (brine displacing CO2) and VW-CA (CO2 displacing air) were 0.145, 0.118 and 0.178, respectively). Neither the Vinegar & Wellington calibration nor our approach has considered CO2 dissolution in brine, which could enhance the CT number of brine and lower the CT number of CO2, rendering higher CT number difference for brine and CO2 than the expected value; hence, both calibration methods tend to overestimate fluid content if CO2 dissolution is significant. Fig. 8 is a collection of all brine and CO2 saturation images reconstructed by our method, exhibiting the distribution of brine and CO2
Fig. 7. Comparison between traditional Vinegar & Wellington calibration (referred as VW in the figure) and the improved calibration (referred as New in the figure), taking CO2 saturation profile of Inj6-C2 stage for example: (a) CO2 saturation map by Vinegar & Wellington calibration; (b) CO2 saturation map by the improved method; (c) comparison of the distribution of saturation values computed by the two methods; (d) comparison of average porosity and CO2 saturation along Zdirection computed by the two methods.
7
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 8. Visualization of CO2 (magenta color) and brine (cyan color) saturation processed by the improved calibration for all flooding stages.
periodically decreased and increased as expected. The preferential flow behaviour was present in almost all image frames. Based on information from the frames upon the end of CO2 flooding (Inj2-C8, Inj4-C6 and Inj6-C5 in Fig. 8b), residual brine occupied lots of pore space across the entire core. Take another look at the frames upon the end of brine flooding (Inj1-B6, Inj3-B6 and Inj5-B3 in Fig. 8a), residual CO2 did not occupy as much pore space as brine. CO2 trap could still be observed in the Inj1-B5, Inj3-B5 and Inj5-B2 frames. This could be the proof of capillary trapping of CO2 that is commonly described as residual CO2 bubbles or patches surrounded by formation brine (Juanes et al., 2010; Krevor et al., 2011; Saadatpoor et al., 2010). Aside from visualization of fluid distribution, average fluid content and saturation values were evaluated for every slice at each frame with the use of our calibration method. Illustrated in Fig. 9 is the sliceaverage brine and CO2 content for some selected frames as compared with the slice-average porosity, and in Fig. 10a-b is the slice-average brine and CO2 saturation for these same frames. The frame-average fluid saturation is given in Fig. 10c as well. The average fluid content and average fluid saturation at each CT slice for selected flooding stages (each note, i.e., Inj1-B1, represents a flooding stage) are demonstrated, and the average fluid contents and average fluid saturations of each flooding stage are plotted along the axis of the core (the z-direction, representing the locations of CT slices from one end to the other end of the core). The porosity curve implied a uniform pore distribution along the longitudinal axis. However, the pores were not evenly distributed at cross-section (X–Y plane), as shown in Fig. 6a. The two segments of the core could be easily recognized from these porosity, fluid content and saturation curves, where two distinct flow patterns were observed at the two sides of the transitional slice. From the main core-body side, the fluid content curves (Fig. 9 and Figure SI8) show that after several flooding cycles, the upper limit of CO2 content dropped from 0.16 (Inj0CO2) to slightly more than 0.1 (Inj6-C5), and the lower limit of CO2 content remained at around 0.05 (Inj3-B6). Meanwhile, the upper limit of brine content stayed somewhere close to 0.11 (Inj3-B6), and the
lower limit was around 0.05 (Inj6-C5). Those curves indicated a substantial residual fluid in the core at the end of each flooding cycle, suggesting that the 100% saturation of a single phase was never achieved after the “Inj0-CO2” stage, although according to our calculation, even at the Inj0-CO2 stage, CO2 did not occupy all the effective pores. In summary, by comparison between Figs. 5 and 10, and based on the information in Fig. 7, our method of fluid saturation computation is more reasonable than the fluid saturation computation based on Vinegar & Wellington calibration. The slice-average and frame-average saturation curves obtained by our approach are in line with displacement physics, the transition between the two segments of the core is clear by looking at the slice-average curves, and the frame-average curves have gradual increasing or decreasing trend within one flooding cycle, which makes sense. 6. Other issues Issues related to the reliability of using CT number to estimate material density have been raised, especially for cone-beam CT in medical applications (Bryant et al., 2012; Coolens and Childs, 2003; Lagravère et al., 2006; Loubele et al., 2006; Sarkis et al., 2007; Yamashina et al., 2008). Due to poor quality and inaccurate Hounsfield units, cone-beam CT images cannot be directly used for adaptive radiation therapy and especially for dose calculation (Kidar et al., 2017; Richter et al., 2008). Efforts to convert the pixel values of cone-beam CT images to Hounsfield units of medical CT images for improved calibration are yet on going, and the accuracy of these modifications depends on the conversion algorithm (Almatani et al., 2016; Kidar et al., 2017; Ming et al., 2008; Richter et al., 2008; Yusuke et al., 2014). The uncertainties associated with the correlation of CT number with calculated porosity and saturation were also addressed in Niu et al. (2015) and Pini et al. (2012). The cause to this uncertainty is a complex issue that requires a thorough inspection of the X-ray attenuation mechanism and CT image 8
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 9. Slice-average fluid content (for full dataset see Fig. SI8) and porosity from the improved calibration: (a) CO2 content; (b) brine content. i μCore (x , y, E ) = μCore =
reconstruction process. The process of reconstructing CT image is actually a process of solving the attenuation coefficient matrix (μ(x,y) for medical fan-beam CT and μ(x,y,z) for cone-beam CT) of a sample by measuring X-ray intensity change after X-ray passing the sample. The Xray beam is polychromatic for all currently available types of CT scanners, and the detectors have an energy-dependent efficiency. The true received intensity by a detector according to Beer’s linear attenuation law can be expressed by Eq. (10).
I=
∫E
Eh
l
L ⎡ dI0 ⋅ε (E )⋅e− ∫0 μ (x , y, E ) dl⎤ dE ⎣ dE ⎦
⎡
(ρA ϕAi
+
ρB ϕBi
i ρA ϕA
+
MA
+
⎢ ρR ϕRi ) ⎢σ (E ) ⎢ ⎣
∑l (fl,A Zl3.8 ,A ) +
ρB ϕBi MB
∑m (fm,B Zm3.8,B) + i
E3.2 ρA ϕA ⎛ b MA
⎝
(10)
+
ρB ϕBi MB
+
ρR ϕRi MR
ρR ϕRi MR
⎞ ⎠
⎤ ∑n (fn,R Zn3.8 ,R ) ⎥ ⎥ ⎥ ⎦
(11)
where the superscript i indicates the parameters are pixel-wise, associated to position (x, y) and local beam energy (E); A, B and R are the subscripts for the two fluids and the rock; Φ is the volumetric fraction; M represents the molar mass. CT images are usually reconstructed with variants of filtered backprojection algorithm. The implication is that the X-ray source provides a narrow beam and an effective energy. As a result, the beam can be regarded as monochromatic, and detectors have constant efficiency. Thus, the true situation (Eq. 10) is idealized into Eq. (12).
where I0 is the incident beam intensity, I is the beam intensity received by detector after penetrating a sample, L is the path length from beam source to detector, ε(E) stands for the energy-dependent detector efficiency, μ(x,y,E) or μ(x,y,z,E) is the energy-dependent attenuation coefficient at any given position of the sample, and E relates to the polychromatic spectrum energy ranging from El to Eh. Taking two-phase immiscible flow as an example, when a core is saturated with two immiscible fluids A and B, the attenuation coefficient of the core according to Eq. (2) is (based on the hypothesis that the fluids and rock matrix are pure substances):
sP =
9
ln (I0 I ) L ∫0 dl
=
L ∫0 μ (x , y ) dl L ∫0 dl
(12)
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 10. Average fluid saturation (for full dataset see Fig. SI9) results from the improved calibration presented in this study: (a) slice-average CO2 saturation; (b) slice-average brine saturation; (c) frame-average saturation. Compared with results in Fig. 5, the average CO2 saturation results and frame-average saturation results from improved calibration make more sense.
images, which require some linear and nonlinear filters to correct. Moreover, even if there exists a characteristic “effective energy”, it takes bold approximation to establish the following linear relation that is equivalent to Eq. (3).
According to Eq. (12), the attenuation coefficient matrix can be derived by mapping the signal from CT projections (sP) discretely at different angular views (of abundant number). However, such idealisation and discrete imaging and mapping will give rise to defected 10
Micron 124 (2019) 102703
X. Miao, et al.
A,B,R
μCore (x , y ) ≈
∑ ph
⎡ ∑l (fl, ph Zl3.8 , ph ) ⎤ ⎫ i i ϕph ρph ⎢σ + ⎥ = E3.2 ⎨ ⎢ ⎥⎬ b ⎣ ⎦⎭ ⎩ ⎧
References
A,B,R
∑ ph
i i ϕph μph
Akin, S., Kovscek, A.R., 2003. Computed Tomography in Petroleum Engineering Research, vol.215. pp. 23–38. https://doi.org/10.1144/GSL.SP.2003.215.01.03. (1). Alemu, B.L., Aker, E., Soldal, M., Johnsen, Ø., Aagaard, P., 2011. Influence of CO2 on rock physics properties in typical reservoir rock: a CO2 flooding experiment of brine saturated sandstone in a CT-scanner. Energy Procedia 4 (1), 4379–4386. Almatani, T., Hugtenburg, R.P., Lewis, R.D., Barley, S.E., Edwards, M.A., 2016. Automated algorithm for CBCT-based dose calculations of prostate radiotherapy with bilateral hip prostheses. Br. J. Radiol. 89 (1066), 20160443. Bryant, J.A., Drage, N.A., Richmond, S., 2012. CT number definition. Radiat. Phys. Chem. 81 (4), 358–361. Cao, P., Karpyn, Z.T., Li, L., 2013. Dynamic alterations in wellbore cement integrity due to geochemical reactions in CO2 -rich environments. Water Resour. Res. 49 (7), 4465–4475. Coolens, C., Childs, P.J., 2003. Calibration of CT Hounsfield units for radiotherapy treatment planning of patients with metallic hip prostheses: the use of the extended CT-scale. Phys. Med. Biol. 48 (11), 1591–1603. Crandall, D., Bromhal, G., 2013. Experimental examination of fluid flow in fractured carbon storage sealing formations. Int. J. Geosci. 04, 1175–1185. https://doi.org/10. 4236/ijg.2013.48111. Crandall, D., Bromhal, G., 2014. Characterization of Experimental Fracture Alteration and Fluid Flow in Fractured Natural Seals. Department of Energy, National Energy Technology Laboratory, USA NRAP-TRS-III-003-2014, NRAP Technical Report Series. Crandall, D., Mcintyre, D., Jarvis, K., Lapeer, R., Tennant, B., 2014. Dynamic Imaging of Multiphase Flows in Rock Using Computed Tomography. American Society of Mechanical Engineers, Fluids Engineering Division. FEDSM 2014https://doi.org/10. 1115/FEDSM2014-21575. Juanes, R., Macminn, C.W., Szulczewski, M.L., 2010. The footprint of the CO2 plume during carbon dioxide storage in saline aquifers: storage efficiency for capillary trapping at the basin scale. Transp. Porous Media 82 (1), 19–30. Kidar, H.S., Hacene, A., Boukellouz, W., Academy, I., 2017. Evaluation of CT to CBCT deformable registration algorithms for adaptive radiation therapy. 5th International Conference on Control Engineering&Information Technology (CEIT-2017), Proceeding of Engineering and Technology 33, 40–44. Krevor, S.C.M., Pini, R., Li, B., Benson, S.M., 2011. Capillary heterogeneity trapping of CO2 in a sandstone rock at reservoir conditions. Geophys. Res. Lett. 38 (15), 98–106. Lagravère, M.O., Fang, Y., Carey, J., Toogood, R.W., Packota, G.V., Major, P.W., 2006. Density conversion factor determined using a cone-beam computed tomography unit NewTom QR-DVT 9000. Radiol. 35 (6), 407–409. Loubele, M., Maes, F., Schutyser, F., Marchal, G., Jacobs, R., Suetens, P., 2006. Assessment of bone segmentation quality of cone-beam CT versus multislice spiral CT: a pilot study. Oral Surg. Oral Med. Oral Pathol. Oral Radiol. Endod. 102 (2), 225–234. Menke, H.P., Reynolds, C.A., Andrew, M.G., Nunes, J.P.P., Bijeljic, B., Blunt, M.J., 2018. 4D multi-scale imaging of reactive flow in carbonates: assessing the impact of heterogeneity on dissolution regimes using streamlines at multiple length scales. Chem. Geol. 481, 27–37. Ming, C., Yaoqin, X., Lei, X., 2008. Auto-propagation of contours for adaptive prostate radiation therapy. Phys. Med. Biol. 53 (17), 4533. Ning, W., Gill, M., Crandall, D., Mcintyre, D., Yan, W., Bruner, K., Li, X., Bromhal, G., 2014. CO2 flooding properties of Liujiagou sandstone: influence of sub-core scale structure heterogeneity. Greenh. Gases Sci. Technol. 4 (3), 400–418. Niu, B., Al‐Menhali, A., Krevor, S.C., 2015. The impact of reservoir conditions on the residual trapping of carbon dioxide in Berea sandstone. Water Resour. Res. 51 (4), 2009–2029. Ott, H., De, K.K., Van, B.M., Vos, F., Van, P.A., Legerstee, P., Bauer, A., Eide, K., Van, d.L.A., Berg, S., 2012. Core-flood experiment for transport of reactive fluids in rocks. Rev. Sci. Instrum. 83 (8), 269. Perrin, J.C., Benson, S., 2010. An experimental study on the influence of sub-core scale heterogeneities on CO2 distribution in reservoir rocks. Transp. Porous Media 82 (1), 93–109. Pini, R., Krevor, S.C.M., Benson, S.M., 2012. Capillary pressure and heterogeneity for the CO2/water system in sandstone rocks at reservoir conditions. Adv. Water Resour. 38 (2), 48–59. Preibisch, S., Saalfeld, S., Schindelin, J., Tomancak, P., 2010. Software for bead-based registration of selective plane illumination microscopy data. Nat. Methods 7 (6), 418. Richter, A., Hu, Q., Steglich, D., Baier, K., Wilbert, J., Guckenberger, M., Flentje, M., 2008. Investigation of the usability of conebeam CT data sets for dose calculation. Radiat. Oncol. 3 (1) (2008-12-16), 3(1), 42-42. Saadatpoor, E., Bryant, S.L., Sepehrnoori, K., 2010. New trapping mechanism in carbon sequestration. Transp. Porous Media 82 (1), 3–17. Sarkis, A., Noujeim, M., Nummikoski, P., Sarkis, A., Noujeim, M., Nummikoski, P., 2007. Bone density measurements in cone-beam computed tomography. Oral Surg. Oral Med. Oral Pathol. Oral Radiol. Endodontol. 103 (2) e53-e53. Schneider, C.A., Rasband, W.S., Eliceiri, K.W., 2012. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9 (7), 671–675. Shi, J.Q., Xue, Z., Durucan, S., 2011. Supercritical CO2 core flooding and imbibition in Tako sandstone—influence of sub-core scale heterogeneity. Int. J. Greenh. Gas Control. 5 (1), 75–87. Smith, C.A., Keiser, D.D., Miller, B.D., Aitkaliyeva, A., 2019. Comparison of manual and automated image analysis techniques for characterization of fission gas pores in irradiated U-Mo fuels. Micron 119, 98–108. U.S. Geological Survey Geologic Carbon Dioxide Storage Resources Assessment Team, 2013. National Assessment of Geologic Carbon Dioxide Storage Resources—Results
(13)
In conclusion, many factors contribute to the uncertainty when correlating CT number with porosity and density, which can hardly be cleared by post-scan image processing. Wellington & Vinegar (Vinegar and Wellington, 1987; Wellington and Vinegar, 1987) paid great attention to pre-calibrating the CT scanner, using fused quartz standard next to the core sample to check instrument drift, and selecting appropriate source energy and dopant dosage to ensure the validity of Eq. (13). Other than these, attentions to the stability of CT machine and selection of filter boards are also needed. In summary, when linear calibration of CT number is implemented, it is necessary to be cautious to ensure the accuracy of calibration results.
7. Conclusions X-ray CT imaging of an alternating CO2/brine displacement flooding experiment was carried out to characterize the two-phase fluid flow and capillary trapping behaviour. The initial trial to map the 3D porosity and saturation profiles by means of Vinegar and Wellington calibration was not successful due to failure to acquire the two 100% saturation images of CO2 and brine. Therefore, an improved calibration method adapted from Vinegar & Wellington calibration was proposed, and comprehensive image processing procedures were described in details. Comparison of the results computed by the two methods was made, and our method produced more reasonable results that were consistent with rock structure and fluid flow physics than Vinegar & Wellington calibration. From the slice-average porosity and saturation curves, the two different segments of the core were distinct. Gradual increasing or decreasing of the saturation for CO2 phase was consistent with the injecting or expelling operations, which was also the case for brine phase. The 3D porosity profile agrees with the laminated structure of Berea sandstone. Meanwhile, from the vertical profiles of reconstructed saturation images, a flat piston-like flooding front was revealed from a homogeneous longitudinal-section, and preferential fingerings were unveiled from the other non-homogeneous longitudinal-section. In the end, issues related to the uncertainty of the linear calibration were raised, causes of the uncertainty were explored, and approaches to alleviate this uncertainty were suggested. In future research works, we plan to conduct similar flooding experiments using a cone-beam microCT at the Institute of Rock and Soil Mechanics, Chinese Academy of Sciences with different core samples and different fluid combinations to test the versatility of our method to process different CT images.
Acknowledgements This work is supported by the Thousand Talent Program for Outstanding Young Scientists (Y731101B01), the National Natural Science Foundation of China project (41807275), CAS-ITRI collaborative research funding (CAS-ITRI2019011), and the China Postdoctoral Science Foundation project (2018M632948). The experiment of X-ray CT imaging of fluid displacement in Berea core samples is supported by an appointment to Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, under the memorandum of CAS-NETL-PNNL. The authors would like to thank Dustin Crandall, Jonathan Moore and Magdalena Gill at National Energy Technology Laboratory (NETL) for their generous help during the CT scanning experiment.
Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.micron.2019.102703. 11
Micron 124 (2019) 102703
X. Miao, et al.
Facial Radiol. 37 (5), 245. Yusuke, O., Noriyuki, K., Yukio, F., Kazuhiro, A., Suguru, D., Ken, T., Kazuma, K., Rei, U., Haruo, M., Keiichi, J., 2014. Evaluation of on-board kV cone beam computed tomography-based dose calculation with deformable image registration using Hounsfield unit modifications. Int. J. Radiat. Oncol. Biol. Phys. 89 (2), 416–423. Zhang, Y., Nishizawa, O., Kiyama, T., Chiyonobu, S., Xue, Z., 2014. Flow behaviour of supercritical CO2 and brine in Berea sandstone during drainage and imbibition revealed by medical X-ray CT images. Geophys. J. Int. 197 (3), 1789–1807.
(ver. 1.1) U.S. Geological Survey Circular. 1386. pp. 41. Vinegar, H.J., Wellington, S.L., 1987. Tomographic imaging of three‐phase flow experiments. Rev. Sci. Instrum. 58 (1), 96–107. Wellington, S.L., Vinegar, H.J., 1987. X-Ray computerized tomography. J. Peterol. Technol. 39 (8), 885–898. https://doi.org/10.2118/16983-PA. Yamashina, A., Tanimoto, K., Sutthiprapaporn, P., Hayakawa, Y., 2008. The reliability of computed tomography (CT) values and dimensional measurements of the oropharyngeal region using cone beam CT: Comparison with multidetector CT. Dento Maxillo
12
Contents lists available at ScienceDirect
Micron journal homepage: www.elsevier.com/locate/micron
Improved Vinegar & Wellington calibration for estimation of fluid saturation and porosity from CT images for a core flooding test under geologic carbon storage conditions
T
⁎
Xiuxiu Miaoa,b,1, Yan Wanga,1, Liwei Zhanga,b, , Ning Weia,b, Xiaochun Lia,b a State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan, 430071, Hubei Province, China b University of Chinese Academy of Sciences, Beijing, 100049, China
A R T I C LE I N FO
A B S T R A C T
Keywords: X-ray CT Image processing Two-phase displacement flooding Saturation mapping Geologic carbon sequestration
X-ray computed tomography (CT) of fluid flow in formation rocks is an important characterization technique in geologic carbon sequestration research to provide insight into the migration and capillary trapping of CO2 under reservoir conditions. An improved calibration method adapted from traditional Vinegar & Wellington calibration is proposed to map the 3D pore and fluid distributions from the CT images of CO2/brine displacement flooding. Similar to Vinegar & Wellington calibration, the proposed method adopts the linear scaling law of CT number transformation to mass density. However, different from Vinegar & Wellington calibration that uses a 100% brine-saturated core image and a 100% CO2-saturated core image as references to calculate CO2 and brine saturations at all time steps, the proposed method uses the CT numbers of CO2 and brine to calculate the incremental of CO2 and brine saturations from time step i to time step i +1. The method is intended for cases in which the two 100% brine saturation and 100% CO2 saturation images can not be successfully obtained. Overall, the improved calibration proposed by this study presents more reasonable results of CO2 and brine distribution in a Berea sandstone core, as compared to traditional Vinegar & Wellington calibration. The reconstructed porosity image agrees with the laminated structure of the Berea sandstone core, and the average porosity evaluated over the entire core (0.176) is comparable to the physical porosity (0.165). Furthermore, the reconstructed saturation images using the improved calibration reveal a flat piston-like flooding front from a homogeneous longitudinal-section of the 3D orthogonal view and preferential fingerings from another nonhomogeneous longitudinal-section, which are not present in the reconstructed saturation images using traditional Vinegar & Wellington calibration. Concerns and causes with respect to the uncertainty of linear CT number calibration are also explained, and approaches to alleviate the uncertainty are suggested.
1. Introduction Successful geologic carbon sequestration requires long-term trapping of CO2 in deep formations without significant CO2 leakage back to the atmosphere. The residual trapping mechanism that depends on the capillary equilibrium between brine and injected CO2 in the porous rock formation is believed to be crucial to determine the storage capacity of CO2 and the long-term fate of injected CO2 (U.S. Geological Survey Geologic Carbon Dioxide Storage Resources Assessment Team, 2013). Also, it is necessary to understand migration behavior of both CO2 and water in subsurface CO2 storage formations. In laboratory studies, core samples with sizes in centimeter range are used in CO2 and
brine flooding experiments under simulated deep reservoir conditions, so as to obtain valuable information about the physical mechanisms that control pressure-driven CO2 and brine migration in porous rocks. X-ray CT technology combined with quantitative image analysis has made visualization of the spatial distribution of CO2 and brine in heterogeneous rock media possible. Moreover, the movement and flow characteristics of CO2 and brine can be revealed from CT image analysis. Mapping saturation levels and porosity within core samples by Xray CT has been practiced since the 1980s. Some of the earliest works were conducted in 1987 by Vinegar and Wellington for two-phase (oil/ water, Wellington and Vinegar, 1987) and three-phase (oil/water/air, Vinegar and Wellington, 1987) displacement flooding. In recent studies,
⁎
Corresponding author. E-mail address: [email protected] (L. Zhang). 1 Co-1st author: Dr. Xiuxiu Miao and Prof. Yan Wang have made equal contribution to this paper. https://doi.org/10.1016/j.micron.2019.102703 Received 7 March 2019; Received in revised form 7 June 2019; Accepted 7 June 2019 Available online 29 June 2019 0968-4328/ © 2019 Elsevier Ltd. All rights reserved.
Micron 124 (2019) 102703
X. Miao, et al.
Zhang et al. (2014) imaged the CO2 and brine flow during drainage and imbibition in a Berea sandstone core using a medical X-ray CT scanner; Menke et al. (2018) used 4D X-ray micro‒computed tomography (μCT) to study the impact of flow heterogeneity on CO2-induced dissolution of limestone cores at both the pore and core scales under reservoir conditions; Cao et al. (2013) quantified changes in wellbore cement integrity by X-ray μCT analysis of pore volume and structure changes of wellbore cement during continuous flooding of CO2-saturated brine. With CT imaging of CO2 flooding of core samples saturated with water, studies have shown that the heterogeneity of pore distribution significantly influences brine and CO2 distributions in rocks (Perrin and Benson, 2010; Shi et al., 2011). The feasibility of using core CT images to compute porosity and fluid saturation, as well as for many other quantification applications, is based on a linear relationship between CT image grayscale intensity (CT number) and material density. For medical CT scans, CT number is expressed in Hounsfield Units (HU) as in Eq. (1)(Akin and Kovscek, 2003; Bryant et al., 2012):
CT = 1000
μ − μW μW
SiA =
(1)
(2)
where ρ is density, σ(E) is the Klein-Nishina coefficient that is presumed “nearly energy-independent” for CT X-ray with energy higher than 100 keV; b is a constant that is approximately 9.8 × 10−24; Σ(fkZk3.8) represents the effective atomic number for a mixture of atomic species with fk and Zk being the fraction of electron and the atomic number of the kth species, respectively; E denotes the X-ray beam energy. As shown in Eq. (1), the CT number of vacuum (zero density) is -1000 HU and the CT number of water is 0 HU by definition. For any other substance, the CT number is energy-dependent. Moreover, for a given substance at a given energy, the CT number is linearly correlated with density. This is the basis for Vinegar & Wellington’s linear calibration for porosity and fluid saturation from CT images, which is achieved by replacing the density of all involved phases with their attenuation coefficients or their CT numbers in the mass and density relationship. For a pixel saturated with fluids A and B in a core image, the relationship is re-written as:
2. Core-flooding CT experiment description A core flooding test that displaced brine in core samples with CO2 was conducted at National Energy Technology Laboratory (NETL), USA using a medical CT scanner with a Berea sandstone core. The Berea core was approximately 5 cm in diameter and 16 cm in length, and the porosity was about 0.165. The temperature of the flooding test was maintained at 67 °C, a confining pressure of 25.6 MPa was applied to the core, and the outlet fluid pressure was maintained at 21.4 MPa. Prior to the flooding cycles, the air-saturated core was scanned as the first reference. After that, supercritical CO2 was injected till the core was fully saturated with CO2, and then the core was scanned as the second reference. Following the supercritical CO2 saturation, the displacement flooding was initiated with a brine injection cycle, and alternated between a CO2 cycle and a brine cycle at constant flow rates of 3.43 ml/min, 1.72 ml/min and 5.13 ml/min, subsequently. Within each flooding cycle, several flooding stages with different amounts of injected fluid were scanned. Each flooding cycle stopped at a stage when no more displacement was detected. In-situ scanning of the core at different CO2/brine saturation stages was simultaneously implemented. To enhance brine-to-CO2 contrast, the brine solution was conditioned with 5 wt% KI to enhance the attenuation coefficient of brine phase, taking advantage of the Kα absorption of X-ray beam by element iodine. The medical CT machine, flooding schematic and the Berea core are shown in Fig. 1. The scanning dataset obtained at different flooding
j μCore = ϕAj μA + ϕBj μ B + ϕRj μRj = ϕ jSAj μA + ϕ j (1 − SAj ) μ B + (1 − ϕ j ) μRj or j CTCore = ϕ jSAj CTA + ϕ j (1 − SAj ) CTB + (1 − ϕ j ) CTRj
(3) where, the superscript j denotes the pixel-wise value; ΦA, ΦB and ΦR are the volume fractions of fluid A, fluid B and rock matrix, respectively; SA is the saturation of fluid A, and the saturation of fluid B equals to 1 minus the saturation of A (i.e. SB = 1 - SA); the attenuation coefficients (CT numbers) of fluid A, B and the rock matrix are denoted as μA, μB and μR (CTA, CTB and CTR), respectively. By imaging a 100% A-saturated core and a 100% B-saturated core, the 3D porosity image can be solved by Eq. (4), and then the saturation image of fluid A or B can be calculated using Eq. (5), given a core image with pores occupied by both A and B.
ϕ=
Bsat Asat CTCore − CTCore CTB − CTA
(5)
The aforementioned calibration (Eq. 4–5), which is referred to as Vinegar & Wellington calibration, has been widely adopted for core porosity and fluid saturation estimation in various CT scanning cases (Akin and Kovscek, 2003; Alemu et al., 2011; Crandall and Bromhal, 2013, 2014; Crandall et al., 2014; Ning et al., 2014; Niu et al., 2015; Ott et al., 2012; Perrin and Benson, 2010; Pini et al., 2012; Zhang et al., 2014) involving different core types (e.g., sandstone, limestone, mudstone, shale, etc.), and different mixtures of miscible or immiscible fluids (e.g., CO2 and water, oil and water, CO2 and H2S, CO2 and N2, etc.) in displacement or co-injection flooding tests. In this study, Vinegar & Wellington calibration was initially employed to estimate porosity and saturations of CO2 and brine from medical CT images obtained from alternating CO2 and brine displacement flooding cycles under laboratory-simulated reservoir conditions. The calibration was based on CT images of Berea sandstone samples with pores filled with CO2 and brine in a CO2 flooding experiment. Step-by-step procedures of CT image analysis were presented. However, Vinegar & Wellington calibration gave unintelligible results of fluid saturation. A possible reason for Vinegar & Wellington calibration to give unsatisfactory results was that certain presumptions for Vinegar & Wellington calibration were not met in our experiment. Therefore, we modified Vinegar & Wellington calibration to analyze our images, and established comprehensive procedures for the reconstruction of porosity and saturation images in hope that our approach can be applied in future studies experiencing similar situations. Our approach uses the CT numbers of both CO2 and brine, and fluid content images of CO2 and brine at time step i to map the fluid content images of CO2 and brine at time step i + 1. The porosity image is obtained by combining the maximum fluid content at each pixel from all time steps. The saturation images of CO2 and brine are then derived by dividing the fluid content images with the porosity image. All raw CT images were processed with the use of ImageJ, a powerful Java-based image processing program (Smith et al., 2019). CO2 and brine saturation maps based on the calculation of the modified Vinegar & Wellington calibration were generated by ImageJ as well. The slice-average porosity, fluid saturation variation along the core and the average fluid saturation for the entire core at each flooding stage were computed and plotted using MATLAB.
where μ is the attenuation coefficient of substance for X-ray beam, and μW is the attenuation coefficient of water. The attenuation coefficient of substance is a function of electron density (or bulk density), effective atomic number, and beam energy. One formula suggested by Vinegar and Wellington (Vinegar and Wellington, 1987; Wellington and Vinegar, 1987) is:
⎡ ⎤ μ = ρ ⎢σ (E ) + b ∑ (fk Zk3.8) E 3.2⎥ k ⎣ ⎦
Bsat Asat − CT iCore CTCore CT iCore − CTCore or SiB = Bsat Asat Bsat Asat CTCore − CTCore CTCore − CTCore
(4) 2
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 1. Core flooding experiment: (a) medical CT scanner; (b) flooding schematic and reconstructed and original images of the core (Φ 5 cm × 16 cm).
slice-average grayscale intensity of the core is taken. Fig. 3 showed the average grayscale intensity at each CT slice for selected flooding stages (each note, i.e., Inj1-B1, represents a flooding stage), and the average grayscale intensities of each flooding stage were plotted along the axis of the core (the z-direction, representing the locations of CT slices from one end to the other end of the core). As shown in Fig. 3, the intensity fluctuates from slice to slice along the Z-direction. Continuous intensity increase (from curves Inj3-B1, B2 to B6) during brine invasion and gradual intensity decrease (from curves Inj4-C2, C3 to C6) during CO2 injection were noticed. Observation of these same curves also suggests that brine injection is much easier than CO2, as the change from Inj4-C6 to Inj5-B1 is much more significant compared to the change from Inj3B6 to Inj4-C2. Another important phenomenon is that all brine invading curves terminate somewhere close to Inj3-B6 (considered as the brinesaturated curve). However, the CO2 invading curves stop at a point much higher than Inj0-CO2 (taken as the CO2-saturated curve). This gap between the CO2-saturated curve and the CO2 invasion terminating curves could be related to residual brine. Further discussion on the displacement pattern is available in the Methodology of porosity and fluid saturation calculation section. In general, this intensity change is consistent with fluid content change. CO2 has higher CT number than air, therefore, when air is displaced by CO2, the average grayscale intensity increases. Similarly, as KI dosed brine has higher CT number than CO2, the core region intensity increases when brine displaces CO2, and decreases when CO2 displaces brine.
Table 1 CT scan dataset obtained at different flooding stages. Flooding cycle
Flow rate
Number of stages
Data label
air saturation CO2 saturation 1st brine flooding 1st CO2 flooding 2nd brine flooding 2nd CO2 flooding 3rd brine flooding 3rd CO2 flooding
/ / 3.43 ml/min
1 1 6 8 6 6 3 5
Inj0-Air Inj0-CO2 Inj1-B1 to Inj2-C1 to Inj3-B1 to Inj4-C1 to Inj5-B1 to Inj6-C1 to
1.72 ml/min 5.13 ml/min
B6 C8 B6 C6 B3 C5
stages is enclosed in Table 1. Each data label in Table 1 relates to a CT scan of a flooding stage. For instance, during the second flooding cycle of brine (Inj3), a total of 6 flooding stages were imaged, and the 3D image data obtained from the first to the sixth flooding stages within this flooding cycle were labelled from Inj3-B1 to Inj3-B6 sequentially. Description of the CT machine and flooding apparatus can be found in Crandall’s previous publications (Crandall and Bromhal, 2013, 2014; Crandall et al., 2014), and the flooding test protocol is also similar to those described in previous papers. 3. CT image description Porosity and saturation calculations were performed on reconstructed CT image sequences. An image sequence consisted of 337 contiguous 2D image slices at a size of 512-pixel by 512-pixel. The pixel resolution was 0.43 mm/pixel and the depth resolution was 0.5 mm/ slice. The 3D view of scanned object could be reconstructed from the image sequence by sequentially stacking contiguous 2D slices. Presented in Fig. 2 are the longitudinal-sectional and cross-sectional views of the reconstructed flooding apparatus and the core. Though the scanning resolution was not adequate to discriminate pore, rock, CO2 and brine in the core region, it was sufficient to reveal the artificial channel at the flooding entrance, different core textures between the flooding entrance and the main core-body. The components that made up the flooding apparatus, including rubber sleeve, confining fluid, and holder, were also revealed. Slices at the beginning and the end of an image sequence were typically covered with heavy noises, and must be excluded from porosity and saturation calculation. For in-situ CT imaging of porous media using constant scanning parameters, we assume that if pore space occupation at a given location does not change, then the “grayscale intensity” (an image term, the same as CT number) at this location remains the same in different “frames” (an image term, means CT images taken at different times, the same meaning as “time steps”). Therefore, any variation of grayscale intensity indicates changes of pore space occupation. In the case of flooding experiment, the intensity change is related to displacement of fluid. To further explore the core region information, a close look at the
4. Methodology of porosity and fluid saturation calculation 4.1. Image pre-processing To compute porosity and fluid saturation is to find a direct conversion of grayscale intensity change to the change of fluid content, porosity and fluid saturation are then derived from estimated fluid content. A difference image representing grayscale intensity change can be generated by subtracting 2 frames of in-situ CT images, for instance, a brine/CO2 flooding frame and an air-saturated frame, as shown in Fig. 4. However, due to long staging time and flooding pressure disturbance, slight drift of the sample occurred, which had to be addressed before image subtraction would work. Apart from sample drift, inevitable CT noises also needed to be evaluated and mitigated when necessary. Due to sample drift, several image frames are not aligned. 3D rigid registration was applied to align all image frames to the first frame (the air-saturated frame), so as to save the trouble of manual selection of problematic frames. Rigid registration is a linear transform of image without scaling; it is one of the most sophisticated practices to align different time-frame images obtained from the same source. We used the descriptor-based image registration plugin (Preibisch et al., 2010) 3
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 2. Anatomy of core flooding CT image, revealing the core, the sleeve, confining fluid, the holder and air: (a) longitudinalsection view; (b) slice with artificial channel from flooding entrance, and typical heavy-noise slice that may present at the front and the end of an image sequence; (c) transitional core slice where artificial channel ends; (d) slice from core-body region, showing texture different from flooding entrance.
Fig. 3. Change of average grayscale intensity along the core during flooding experiment (partial dataset; for full dataset see Fig. SI5).
window sizes were tested on the difference images to remove noises. It is believed that the radius 1 (pixel) Gaussian filter is adequate after comparing the Gaussian filtered images with different window sizes. The comparison is elaborated in Section 2 of the Supplementary information.
from Fiji (an ImageJ distribution) (Schneider et al., 2012) to implement the 3D rigid registration. Details on how the 3D rigid registration was implemented and how the registration improved image alignment are explained in Section 1 of the Supplementary information. Since all CT images are liable to noise, and all noise filtering algorithms will inevitably altering true intensity signals to some degree, we suggest a conservative noise-filtering practice, which is to apply filtering as rare and as mild as possible. In this study, no filtering was implemented before image registration, so that the difference signal is maximally preserved. After registration, Gaussian filters of different
4.2. Fluid saturation by vinegar & Wellington calibration According to Vinegar & Wellington linear calibration (Eq. 4–5), we could calculate core porosity and saturations of brine and CO2 in the 4
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 4. Application of the subtraction method to create difference image: (a) brine/CO2 invaded core image; (b) air saturated core image; (c) difference image.
core image obtained at an arbitrary stage i (SiCO2 and SiBrine) by the following equation: Brine-sat Air-sat − CTCore ⎧ ϕ = CTCore ⎪ BA CTBrine − CTAir ⎪ Brine-sat 2 -sat ⎪ − CTCO CTCore Core ϕBC = CTBrine − CTCO2 ⎨ ⎪ CO2 -sat Air-sat − CTCore CT ⎪ϕ = Core ⎪ CA CTCO2 − CTAir ⎩
calibration does not need the 100% CO2 and brine saturated images. If the two 100% saturation images for reference are correctly acquired, then ΔθiBrine = -ΔθiCO2 and θ0CO2 = φ, and this method is consistent with Vinegar & Wellington calibration.
i Brine-sat − CT Core CTCore ⎫ Brine-sat 2 -sat ⎪ − CTCO CTCore ⎪ Core CO2 -sat ⎬ i CT Core − CTCore ⎪ = Brine-sat 2 -sat ⎪ − CTCO CTCore Core ⎭
i = SCO 2
SiBrine
Inj0-CO2 Air-sat − CTCore CTCore ⎧ 0 ⎪ θCO2 = CTCO2 − CTAir ⎪ +1 = θiCO θiCO2 + ΔθiCO2 2 ⎨ ⎪ θiCO2 SiCO2 = ⎪ ϕ ⎩
(6)
where, φBA is the porosity obtained at the stage of brine displacing air; φBC is the porosity obtained at the stage of brine displacing CO2; φCA is the porosity obtained at the stage of CO2 displacing air. Theoretically φBA = φBC = φCA. The CT numbers of KI-dosed brine (CTBrine), CO2 (CT CO2) and air (CTAir) are 735, -400, -1000, respectively. In Fig. 5a, the porosity calculated using different linear calibration formula did not conform to each other, and the CO2 saturation significantly fluctuated along the core (Fig. 5b). The reason that the saturation values calibrated by Vinegar and Wellington method have significant variation is explained in the chessboard simulation section (Section 3 of the SI). Also, frame-average fluid saturation within one CO2 or brine injection cycle did not progressively increase or decrease as expected (Fig. 5c). For Vinegar & Wellington calibration to work, the assumptions are that there are two 100% saturation states established for brine and CO2, respectively. All effective pores in the single-phase referential saturation image must be occupied by only one fluid phase. However, in our case, those presumptions are not met. It is possible that in our experiment, not all effective pores were occupied by CO2 when the CO2 saturated image (Inj0-CO2) was obtained, and that the brine saturated image (Inj3-B6) was much likely acquired at a co-saturation state of both brine and CO2. The fact that no CO2 was detected at the flow outlet during a brine displacing CO2 stage does not ensure no residual CO2 in the brine-saturated core. The fact is also the case for CO2 displacing air or CO2 displacing brine scenarios. We simulated the core matrix with a 6-by-6 chessboard to demonstrate the possibility for the predicted values of porosity and saturation to go beyond the 0–1 limit using the Vinegar & Wellington calibration. The description of the chessboard simulation is provided in Section 3 of the Supplementary information.
0 =0 θBrine ⎧ i ⎪ θi + 1 = θi Brine Brine + ΔθBrine i ⎨ θ SiBrine = Brine ⎪ ϕ ⎩
(7)
The implementation of our calibration is described in this paragraph. Based on flooding sequence (Fig. 1b & Table 1), the initial brine and CO2 contents (θ0Brine and θ0CO2) were computed for the first CO2 displacing air frame (Inj0-CO2). Then incremental of brine or CO2 content at imaging stage i+1 against its previous stage i was evaluated from Eq. (8). By simultaneously comparing pixel j at stage i+1 with its previous stage i, and pixel j at the air-saturated stage, three fluid displacement patterns from stage i to stage i+1 at pixel j, including CO2brine displacement, brine-air displacement, and CO2-air displacement, were recognized and computed. Due to system noises, even air intensity fluctuated from frame to frame (Figure SI7), and the relative mild Gaussian filter could not completely remove this fluctuation (Figure SI7). Therefore, we introduced a tolerance (Tol = 0.05) to Eq. (8) to address the fluctuation. (8) The porosity image was calculated from the maximum fluid content at each pixel j from all frames, as expressed in Eq. (9). The formula meant that the estimation was based on a record of all pores once occupied by fluid during the whole displacement flooding process.
ϕ = max(θ j ) ∀i
(9)
5. Results and discussion Following the procedures of registration, filtering and reconstruction of fluid content images and the porosity image, the saturation images were then computed by dividing fluid content images with the porosity image. Presented in Fig. 6 are the porosity, brine and CO2 saturation images at one CO2 flooding stage (Inj4-C4) processed by the improved calibration described in this paper. The orthogonal view of the porosity image (Fig. 6a) revealed the intrinsic heterogeneity of the core. The longitudinal section was marked by a strip-like pattern that was parallel to the main flow direction, which agreed with the laminated structure of Berea sandstone (Fig. 1b). It was also obvious that most pixels concurrently hosted pore space and the rock material, though a few pixels had porosity values very close to 0 or very close to 1. Due to heterogeneous porosity distribution, brine and CO2 flow behaviors shown in different longitudinal-sections (Fig. 6b-c) varied significantly from each other. A flat piston-like flooding front was observed on the left
4.3. Description of the new approach to reconstruct porosity and fluid content images To handle the failure patterns when applying Vinegar & Wellington calibration to process our images, the authors suggested a new calibration method to reconstruct porosity and fluid content images. Given the fact that the 100% brine saturation frame and the 100% CO2 saturation frame could not be acquired in our flooding experiment, and there was possible air-phase interference, the authors suggested using an incremental calibration (Eq. 7) to estimate brine and CO2 saturation. The advantage of this new incremental calibration is that the 5
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 5. Porosity and fluid saturation converted from CT number by traditional Vinegar & Wellington calibration: (a) porosity; (b) CO2 saturation (partial dataset; for full dataset see Fig. SI6); (c) frame-average CO2 and brine saturation. Poor average CO2 saturation results and frame-average saturation results from direct linear calibration can be clearly seen.
longitudinal-sections in Fig. 6b-c, while preferential fingering dominated the right longitudinal-sections in Fig. 6b-c. The CO2 fingering was able to penetrate through the brine phase, and much less brine was displaced on the longitudinal-section with preferential fingering compared to the longitudinal-section with piston-like flooding front. Fig. 7 serves as a comparison between results of Vinegar & Wellington calibration and the calibration method proposed in this paper. The CO2 saturation profile for the frame (Inj6-C2) computed by our proposed method exhibited pixels showing clear and distinguishable flooding signals (Fig. 7b, saturation value is cut off to the range of 0 to 1). The actual saturation value distribution from -1.1 to 1 is shown in
Fig. 7c, demonstrating that more than 99% pixels have saturation values from 0 to 1. In contrast, the saturation images obtained by Vinegar & Wellington calibration were spoiled by extremes, and barely showed any distinguishable signals (Fig. 7a, saturation value is cut off to the range of −100 to 100). The actual saturation value distribution from −695 to 613 is shown in Fig. 7c, demonstrating that less than 45% pixels have saturation values from 0 to 1. These extreme values mainly originated from the residual CO2 pixels in the referential brine saturation image, described as the first failure pattern of Vinegar & Wellington calibration in the chessboard simulation (Section 3 of the Supplementary information). The slice-average saturation in Fig. 7d also 6
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 6. Orthogonal view of brine and CO2 saturation in a CO2 flooding stage (Inj4-C4) processed by the improved calibration: (a) the porosity image; (b) brine saturation image; (c) CO2 saturation image. In each image, two longitudinal sections and two cross sections are shown.
in the core at different flooding stages and injection rates. It should be noted that the images were not strictly acquired at the same interval of time, and the injection rate within each injection cycle (corresponding to a row of images) also differed from each other due to injection pressure variation. For instance, during the fifth flooding cycle of brine (Inj5-B1 to B3), brine swiftly took up most of the pores at the first injection (that was the reason why only 3 images were obtained), while at lower injection rate 6 images were acquired for the third flooding cycle of brine (Inj3-B1 to B6). Therefore, we can observe that the brine saturation on the Inj5-B1 frame was much higher than that on Inj3-B1. Nevertheless, the cyclic fluid intrusion/recession during the flooding cycle was clear. Barely any brine was injected at the Inj1-B1 flooding stage, and CO2 saturation was initially high, yet did not reach 100% saturation. Following the flooding cycles, the amount of CO2 in the core
shows larger variation produced by Vinegar & Wellington calibration, as compared to our method. Meanwhile, the slice-average porosity of the Berea core obtained from our calibration approach (about 0.176) was closer to the physical porosity (0.165) than the average porosities from VW method (the average porosities of VW-BA (brine displacing air), VW-BC (brine displacing CO2) and VW-CA (CO2 displacing air) were 0.145, 0.118 and 0.178, respectively). Neither the Vinegar & Wellington calibration nor our approach has considered CO2 dissolution in brine, which could enhance the CT number of brine and lower the CT number of CO2, rendering higher CT number difference for brine and CO2 than the expected value; hence, both calibration methods tend to overestimate fluid content if CO2 dissolution is significant. Fig. 8 is a collection of all brine and CO2 saturation images reconstructed by our method, exhibiting the distribution of brine and CO2
Fig. 7. Comparison between traditional Vinegar & Wellington calibration (referred as VW in the figure) and the improved calibration (referred as New in the figure), taking CO2 saturation profile of Inj6-C2 stage for example: (a) CO2 saturation map by Vinegar & Wellington calibration; (b) CO2 saturation map by the improved method; (c) comparison of the distribution of saturation values computed by the two methods; (d) comparison of average porosity and CO2 saturation along Zdirection computed by the two methods.
7
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 8. Visualization of CO2 (magenta color) and brine (cyan color) saturation processed by the improved calibration for all flooding stages.
periodically decreased and increased as expected. The preferential flow behaviour was present in almost all image frames. Based on information from the frames upon the end of CO2 flooding (Inj2-C8, Inj4-C6 and Inj6-C5 in Fig. 8b), residual brine occupied lots of pore space across the entire core. Take another look at the frames upon the end of brine flooding (Inj1-B6, Inj3-B6 and Inj5-B3 in Fig. 8a), residual CO2 did not occupy as much pore space as brine. CO2 trap could still be observed in the Inj1-B5, Inj3-B5 and Inj5-B2 frames. This could be the proof of capillary trapping of CO2 that is commonly described as residual CO2 bubbles or patches surrounded by formation brine (Juanes et al., 2010; Krevor et al., 2011; Saadatpoor et al., 2010). Aside from visualization of fluid distribution, average fluid content and saturation values were evaluated for every slice at each frame with the use of our calibration method. Illustrated in Fig. 9 is the sliceaverage brine and CO2 content for some selected frames as compared with the slice-average porosity, and in Fig. 10a-b is the slice-average brine and CO2 saturation for these same frames. The frame-average fluid saturation is given in Fig. 10c as well. The average fluid content and average fluid saturation at each CT slice for selected flooding stages (each note, i.e., Inj1-B1, represents a flooding stage) are demonstrated, and the average fluid contents and average fluid saturations of each flooding stage are plotted along the axis of the core (the z-direction, representing the locations of CT slices from one end to the other end of the core). The porosity curve implied a uniform pore distribution along the longitudinal axis. However, the pores were not evenly distributed at cross-section (X–Y plane), as shown in Fig. 6a. The two segments of the core could be easily recognized from these porosity, fluid content and saturation curves, where two distinct flow patterns were observed at the two sides of the transitional slice. From the main core-body side, the fluid content curves (Fig. 9 and Figure SI8) show that after several flooding cycles, the upper limit of CO2 content dropped from 0.16 (Inj0CO2) to slightly more than 0.1 (Inj6-C5), and the lower limit of CO2 content remained at around 0.05 (Inj3-B6). Meanwhile, the upper limit of brine content stayed somewhere close to 0.11 (Inj3-B6), and the
lower limit was around 0.05 (Inj6-C5). Those curves indicated a substantial residual fluid in the core at the end of each flooding cycle, suggesting that the 100% saturation of a single phase was never achieved after the “Inj0-CO2” stage, although according to our calculation, even at the Inj0-CO2 stage, CO2 did not occupy all the effective pores. In summary, by comparison between Figs. 5 and 10, and based on the information in Fig. 7, our method of fluid saturation computation is more reasonable than the fluid saturation computation based on Vinegar & Wellington calibration. The slice-average and frame-average saturation curves obtained by our approach are in line with displacement physics, the transition between the two segments of the core is clear by looking at the slice-average curves, and the frame-average curves have gradual increasing or decreasing trend within one flooding cycle, which makes sense. 6. Other issues Issues related to the reliability of using CT number to estimate material density have been raised, especially for cone-beam CT in medical applications (Bryant et al., 2012; Coolens and Childs, 2003; Lagravère et al., 2006; Loubele et al., 2006; Sarkis et al., 2007; Yamashina et al., 2008). Due to poor quality and inaccurate Hounsfield units, cone-beam CT images cannot be directly used for adaptive radiation therapy and especially for dose calculation (Kidar et al., 2017; Richter et al., 2008). Efforts to convert the pixel values of cone-beam CT images to Hounsfield units of medical CT images for improved calibration are yet on going, and the accuracy of these modifications depends on the conversion algorithm (Almatani et al., 2016; Kidar et al., 2017; Ming et al., 2008; Richter et al., 2008; Yusuke et al., 2014). The uncertainties associated with the correlation of CT number with calculated porosity and saturation were also addressed in Niu et al. (2015) and Pini et al. (2012). The cause to this uncertainty is a complex issue that requires a thorough inspection of the X-ray attenuation mechanism and CT image 8
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 9. Slice-average fluid content (for full dataset see Fig. SI8) and porosity from the improved calibration: (a) CO2 content; (b) brine content. i μCore (x , y, E ) = μCore =
reconstruction process. The process of reconstructing CT image is actually a process of solving the attenuation coefficient matrix (μ(x,y) for medical fan-beam CT and μ(x,y,z) for cone-beam CT) of a sample by measuring X-ray intensity change after X-ray passing the sample. The Xray beam is polychromatic for all currently available types of CT scanners, and the detectors have an energy-dependent efficiency. The true received intensity by a detector according to Beer’s linear attenuation law can be expressed by Eq. (10).
I=
∫E
Eh
l
L ⎡ dI0 ⋅ε (E )⋅e− ∫0 μ (x , y, E ) dl⎤ dE ⎣ dE ⎦
⎡
(ρA ϕAi
+
ρB ϕBi
i ρA ϕA
+
MA
+
⎢ ρR ϕRi ) ⎢σ (E ) ⎢ ⎣
∑l (fl,A Zl3.8 ,A ) +
ρB ϕBi MB
∑m (fm,B Zm3.8,B) + i
E3.2 ρA ϕA ⎛ b MA
⎝
(10)
+
ρB ϕBi MB
+
ρR ϕRi MR
ρR ϕRi MR
⎞ ⎠
⎤ ∑n (fn,R Zn3.8 ,R ) ⎥ ⎥ ⎥ ⎦
(11)
where the superscript i indicates the parameters are pixel-wise, associated to position (x, y) and local beam energy (E); A, B and R are the subscripts for the two fluids and the rock; Φ is the volumetric fraction; M represents the molar mass. CT images are usually reconstructed with variants of filtered backprojection algorithm. The implication is that the X-ray source provides a narrow beam and an effective energy. As a result, the beam can be regarded as monochromatic, and detectors have constant efficiency. Thus, the true situation (Eq. 10) is idealized into Eq. (12).
where I0 is the incident beam intensity, I is the beam intensity received by detector after penetrating a sample, L is the path length from beam source to detector, ε(E) stands for the energy-dependent detector efficiency, μ(x,y,E) or μ(x,y,z,E) is the energy-dependent attenuation coefficient at any given position of the sample, and E relates to the polychromatic spectrum energy ranging from El to Eh. Taking two-phase immiscible flow as an example, when a core is saturated with two immiscible fluids A and B, the attenuation coefficient of the core according to Eq. (2) is (based on the hypothesis that the fluids and rock matrix are pure substances):
sP =
9
ln (I0 I ) L ∫0 dl
=
L ∫0 μ (x , y ) dl L ∫0 dl
(12)
Micron 124 (2019) 102703
X. Miao, et al.
Fig. 10. Average fluid saturation (for full dataset see Fig. SI9) results from the improved calibration presented in this study: (a) slice-average CO2 saturation; (b) slice-average brine saturation; (c) frame-average saturation. Compared with results in Fig. 5, the average CO2 saturation results and frame-average saturation results from improved calibration make more sense.
images, which require some linear and nonlinear filters to correct. Moreover, even if there exists a characteristic “effective energy”, it takes bold approximation to establish the following linear relation that is equivalent to Eq. (3).
According to Eq. (12), the attenuation coefficient matrix can be derived by mapping the signal from CT projections (sP) discretely at different angular views (of abundant number). However, such idealisation and discrete imaging and mapping will give rise to defected 10
Micron 124 (2019) 102703
X. Miao, et al.
A,B,R
μCore (x , y ) ≈
∑ ph
⎡ ∑l (fl, ph Zl3.8 , ph ) ⎤ ⎫ i i ϕph ρph ⎢σ + ⎥ = E3.2 ⎨ ⎢ ⎥⎬ b ⎣ ⎦⎭ ⎩ ⎧
References
A,B,R
∑ ph
i i ϕph μph
Akin, S., Kovscek, A.R., 2003. Computed Tomography in Petroleum Engineering Research, vol.215. pp. 23–38. https://doi.org/10.1144/GSL.SP.2003.215.01.03. (1). Alemu, B.L., Aker, E., Soldal, M., Johnsen, Ø., Aagaard, P., 2011. Influence of CO2 on rock physics properties in typical reservoir rock: a CO2 flooding experiment of brine saturated sandstone in a CT-scanner. Energy Procedia 4 (1), 4379–4386. Almatani, T., Hugtenburg, R.P., Lewis, R.D., Barley, S.E., Edwards, M.A., 2016. Automated algorithm for CBCT-based dose calculations of prostate radiotherapy with bilateral hip prostheses. Br. J. Radiol. 89 (1066), 20160443. Bryant, J.A., Drage, N.A., Richmond, S., 2012. CT number definition. Radiat. Phys. Chem. 81 (4), 358–361. Cao, P., Karpyn, Z.T., Li, L., 2013. Dynamic alterations in wellbore cement integrity due to geochemical reactions in CO2 -rich environments. Water Resour. Res. 49 (7), 4465–4475. Coolens, C., Childs, P.J., 2003. Calibration of CT Hounsfield units for radiotherapy treatment planning of patients with metallic hip prostheses: the use of the extended CT-scale. Phys. Med. Biol. 48 (11), 1591–1603. Crandall, D., Bromhal, G., 2013. Experimental examination of fluid flow in fractured carbon storage sealing formations. Int. J. Geosci. 04, 1175–1185. https://doi.org/10. 4236/ijg.2013.48111. Crandall, D., Bromhal, G., 2014. Characterization of Experimental Fracture Alteration and Fluid Flow in Fractured Natural Seals. Department of Energy, National Energy Technology Laboratory, USA NRAP-TRS-III-003-2014, NRAP Technical Report Series. Crandall, D., Mcintyre, D., Jarvis, K., Lapeer, R., Tennant, B., 2014. Dynamic Imaging of Multiphase Flows in Rock Using Computed Tomography. American Society of Mechanical Engineers, Fluids Engineering Division. FEDSM 2014https://doi.org/10. 1115/FEDSM2014-21575. Juanes, R., Macminn, C.W., Szulczewski, M.L., 2010. The footprint of the CO2 plume during carbon dioxide storage in saline aquifers: storage efficiency for capillary trapping at the basin scale. Transp. Porous Media 82 (1), 19–30. Kidar, H.S., Hacene, A., Boukellouz, W., Academy, I., 2017. Evaluation of CT to CBCT deformable registration algorithms for adaptive radiation therapy. 5th International Conference on Control Engineering&Information Technology (CEIT-2017), Proceeding of Engineering and Technology 33, 40–44. Krevor, S.C.M., Pini, R., Li, B., Benson, S.M., 2011. Capillary heterogeneity trapping of CO2 in a sandstone rock at reservoir conditions. Geophys. Res. Lett. 38 (15), 98–106. Lagravère, M.O., Fang, Y., Carey, J., Toogood, R.W., Packota, G.V., Major, P.W., 2006. Density conversion factor determined using a cone-beam computed tomography unit NewTom QR-DVT 9000. Radiol. 35 (6), 407–409. Loubele, M., Maes, F., Schutyser, F., Marchal, G., Jacobs, R., Suetens, P., 2006. Assessment of bone segmentation quality of cone-beam CT versus multislice spiral CT: a pilot study. Oral Surg. Oral Med. Oral Pathol. Oral Radiol. Endod. 102 (2), 225–234. Menke, H.P., Reynolds, C.A., Andrew, M.G., Nunes, J.P.P., Bijeljic, B., Blunt, M.J., 2018. 4D multi-scale imaging of reactive flow in carbonates: assessing the impact of heterogeneity on dissolution regimes using streamlines at multiple length scales. Chem. Geol. 481, 27–37. Ming, C., Yaoqin, X., Lei, X., 2008. Auto-propagation of contours for adaptive prostate radiation therapy. Phys. Med. Biol. 53 (17), 4533. Ning, W., Gill, M., Crandall, D., Mcintyre, D., Yan, W., Bruner, K., Li, X., Bromhal, G., 2014. CO2 flooding properties of Liujiagou sandstone: influence of sub-core scale structure heterogeneity. Greenh. Gases Sci. Technol. 4 (3), 400–418. Niu, B., Al‐Menhali, A., Krevor, S.C., 2015. The impact of reservoir conditions on the residual trapping of carbon dioxide in Berea sandstone. Water Resour. Res. 51 (4), 2009–2029. Ott, H., De, K.K., Van, B.M., Vos, F., Van, P.A., Legerstee, P., Bauer, A., Eide, K., Van, d.L.A., Berg, S., 2012. Core-flood experiment for transport of reactive fluids in rocks. Rev. Sci. Instrum. 83 (8), 269. Perrin, J.C., Benson, S., 2010. An experimental study on the influence of sub-core scale heterogeneities on CO2 distribution in reservoir rocks. Transp. Porous Media 82 (1), 93–109. Pini, R., Krevor, S.C.M., Benson, S.M., 2012. Capillary pressure and heterogeneity for the CO2/water system in sandstone rocks at reservoir conditions. Adv. Water Resour. 38 (2), 48–59. Preibisch, S., Saalfeld, S., Schindelin, J., Tomancak, P., 2010. Software for bead-based registration of selective plane illumination microscopy data. Nat. Methods 7 (6), 418. Richter, A., Hu, Q., Steglich, D., Baier, K., Wilbert, J., Guckenberger, M., Flentje, M., 2008. Investigation of the usability of conebeam CT data sets for dose calculation. Radiat. Oncol. 3 (1) (2008-12-16), 3(1), 42-42. Saadatpoor, E., Bryant, S.L., Sepehrnoori, K., 2010. New trapping mechanism in carbon sequestration. Transp. Porous Media 82 (1), 3–17. Sarkis, A., Noujeim, M., Nummikoski, P., Sarkis, A., Noujeim, M., Nummikoski, P., 2007. Bone density measurements in cone-beam computed tomography. Oral Surg. Oral Med. Oral Pathol. Oral Radiol. Endodontol. 103 (2) e53-e53. Schneider, C.A., Rasband, W.S., Eliceiri, K.W., 2012. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9 (7), 671–675. Shi, J.Q., Xue, Z., Durucan, S., 2011. Supercritical CO2 core flooding and imbibition in Tako sandstone—influence of sub-core scale heterogeneity. Int. J. Greenh. Gas Control. 5 (1), 75–87. Smith, C.A., Keiser, D.D., Miller, B.D., Aitkaliyeva, A., 2019. Comparison of manual and automated image analysis techniques for characterization of fission gas pores in irradiated U-Mo fuels. Micron 119, 98–108. U.S. Geological Survey Geologic Carbon Dioxide Storage Resources Assessment Team, 2013. National Assessment of Geologic Carbon Dioxide Storage Resources—Results
(13)
In conclusion, many factors contribute to the uncertainty when correlating CT number with porosity and density, which can hardly be cleared by post-scan image processing. Wellington & Vinegar (Vinegar and Wellington, 1987; Wellington and Vinegar, 1987) paid great attention to pre-calibrating the CT scanner, using fused quartz standard next to the core sample to check instrument drift, and selecting appropriate source energy and dopant dosage to ensure the validity of Eq. (13). Other than these, attentions to the stability of CT machine and selection of filter boards are also needed. In summary, when linear calibration of CT number is implemented, it is necessary to be cautious to ensure the accuracy of calibration results.
7. Conclusions X-ray CT imaging of an alternating CO2/brine displacement flooding experiment was carried out to characterize the two-phase fluid flow and capillary trapping behaviour. The initial trial to map the 3D porosity and saturation profiles by means of Vinegar and Wellington calibration was not successful due to failure to acquire the two 100% saturation images of CO2 and brine. Therefore, an improved calibration method adapted from Vinegar & Wellington calibration was proposed, and comprehensive image processing procedures were described in details. Comparison of the results computed by the two methods was made, and our method produced more reasonable results that were consistent with rock structure and fluid flow physics than Vinegar & Wellington calibration. From the slice-average porosity and saturation curves, the two different segments of the core were distinct. Gradual increasing or decreasing of the saturation for CO2 phase was consistent with the injecting or expelling operations, which was also the case for brine phase. The 3D porosity profile agrees with the laminated structure of Berea sandstone. Meanwhile, from the vertical profiles of reconstructed saturation images, a flat piston-like flooding front was revealed from a homogeneous longitudinal-section, and preferential fingerings were unveiled from the other non-homogeneous longitudinal-section. In the end, issues related to the uncertainty of the linear calibration were raised, causes of the uncertainty were explored, and approaches to alleviate this uncertainty were suggested. In future research works, we plan to conduct similar flooding experiments using a cone-beam microCT at the Institute of Rock and Soil Mechanics, Chinese Academy of Sciences with different core samples and different fluid combinations to test the versatility of our method to process different CT images.
Acknowledgements This work is supported by the Thousand Talent Program for Outstanding Young Scientists (Y731101B01), the National Natural Science Foundation of China project (41807275), CAS-ITRI collaborative research funding (CAS-ITRI2019011), and the China Postdoctoral Science Foundation project (2018M632948). The experiment of X-ray CT imaging of fluid displacement in Berea core samples is supported by an appointment to Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, under the memorandum of CAS-NETL-PNNL. The authors would like to thank Dustin Crandall, Jonathan Moore and Magdalena Gill at National Energy Technology Laboratory (NETL) for their generous help during the CT scanning experiment.
Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.micron.2019.102703. 11
Micron 124 (2019) 102703
X. Miao, et al.
Facial Radiol. 37 (5), 245. Yusuke, O., Noriyuki, K., Yukio, F., Kazuhiro, A., Suguru, D., Ken, T., Kazuma, K., Rei, U., Haruo, M., Keiichi, J., 2014. Evaluation of on-board kV cone beam computed tomography-based dose calculation with deformable image registration using Hounsfield unit modifications. Int. J. Radiat. Oncol. Biol. Phys. 89 (2), 416–423. Zhang, Y., Nishizawa, O., Kiyama, T., Chiyonobu, S., Xue, Z., 2014. Flow behaviour of supercritical CO2 and brine in Berea sandstone during drainage and imbibition revealed by medical X-ray CT images. Geophys. J. Int. 197 (3), 1789–1807.
(ver. 1.1) U.S. Geological Survey Circular. 1386. pp. 41. Vinegar, H.J., Wellington, S.L., 1987. Tomographic imaging of three‐phase flow experiments. Rev. Sci. Instrum. 58 (1), 96–107. Wellington, S.L., Vinegar, H.J., 1987. X-Ray computerized tomography. J. Peterol. Technol. 39 (8), 885–898. https://doi.org/10.2118/16983-PA. Yamashina, A., Tanimoto, K., Sutthiprapaporn, P., Hayakawa, Y., 2008. The reliability of computed tomography (CT) values and dimensional measurements of the oropharyngeal region using cone beam CT: Comparison with multidetector CT. Dento Maxillo
12