Capillary pressure correction in irregularly shaped pore ... - ISSMGE
Capillary pressure correction in irregularly shaped pore channel under fully wetting condition Correction de la pression capillaire dans un canal poreux de forme irrégulière sous un état de mouillage complet Hyoung Suk Suh, Dong Hun Kang, Tae Sup Yun School of Civil and Environmental Engineering, Yonsei University, Korea, [email protected]
ABSTRACT: Pore network modeling (PNM) has been useful to simulate water retention process of the capillary dominated systems. Randomly structured 3D pore geometry can be simplified as the network of pore channels of cylindrical shape. Radii of cylindrical pore channels are determined by minimum inscribing spheres which lead to the overestimation of capillary pressures of irregularly shaped pore throats. The correct capillary pressure can be obtained by bubble analysis using lattice Boltzmann method (LBM). Pore domain larger than its representative elementary volume, however, contains more than hundreds of pore throats. LB bubble simulations for estimating corrected capillary pressures of every pore throats are computationally demanding. Therefore, this study presents a simple shape analysis technique by Euclidean distance transform (EDT) to identify the effective cross-sectional area which is defined as the expected location of the nonwetting fluid during the invasion into a capillary tube filled with wetting fluid. By extending MS-P theory, we calculate capillary pressures of tubes whose effective cross-sections are taken from 3D X-ray CT image compared with LB bubble simulation. We also investigate the applicability of the shape analysis technique by calibrating the pore network drainage simulation results. RÉSUMÉ: La modélisation du réseau de pores (PNM) a été utile pour simuler le processus de rétention d'eau des systèmes capillaires dominés. La géométrie des pores 3D à structure aléatoire peut être simplifiée en tant que réseau de canaux poreux de forme cylindrique. Les rayons des canaux de pores cylindriques sont déterminés par des sphères d'inscription minimales qui conduisent à une surestimation des pressions capillaires de gorges de pores de forme irrégulière. La pression capillaire correcte peut être obtenue par analyse de bulles en utilisant la méthode de Boltzmann en treillis (LBM). Un domaine de pores plus grand que son volume élémentaire représentatif, cependant, contient plus de centaines de gorges de pores. Les simulations de bulles LB pour estimer les pressions capillaires corrigées de chaque gorge de pore sont très exigeantes en termes de calcul. Par conséquent, cette étude présente une technique d'analyse de forme simple par la transformée de distance euclidienne (EDT) pour identifier la section transversale efficace qui est définie comme l'emplacement attendu du fluide non mouillant pendant l'invasion dans un tube capillaire rempli de fluide mouillant. En prolongeant la théorie MS-P, nous calculons les pressions capillaires de tubes dont les coupes efficaces sont prises à partir de l'image de tomodensitométrie 3D par rapport à la simulation de bulles LB. Nous étudions également l 'applicabilité de la technique d' analyse de forme en calibrant les résultats de simulation de drainage de réseau de pores. KEYWORDS: Capillary pressure, Pore throat geometry, Pore network model, Lattice Boltzmann method 1 INTRODUCTION Carbon dioxide storage (Al-Menhali and Krevor 2016) or gas production from hydrate-bearing sediments (Dai and Seol 2014, Jang and Santamarina 2014) involves capillary-driven two phase fluid flow. Such systems include functional expressions for capillary pressure pc and wetting fluid saturation Sw as the governing equations (Brooks and Corey 1964, van Genuchten 1980). Efforts has been made to develop numerical models to capture pc-Sw relationship as an alternative of time consuming laboratory experiments (Frette and Helland 2010). It is well known that pore-scale simulators such as lattice Boltzmann models (LBM), (Liu et al. 2012) are computationally expensive compared to pore network models (PNM) due to the discretization methods (Joekar-Niasar and Hassanizadeh 2012). Advances in 3D imaging technology allows the direct extraction of pore network from 3D image, which is composed of pore chambers and pore throats, however, pore throats have been usually modeled as cylindrical capillary tubes with radius of minimum inscribing sphere rins (Silin and Patzek 2006) which leads to an overestimation of capillary entry pressure by adopting the Young-Laplace equation as follows: pc ,ins
2Ts cos rins
(1)
where Ts and θ indicate the interfacial tension and contact angle, respectively. Entry pressure for fluid invasion in tube having constant crosssection along the longitudinal direction is usually calculated by adopting Mayer and Stowe-Princen theory (MS-P theory) which is based on the energy balance (Mayer and Stowe 1965, Princen 1969a, 1969b, 1970). Analytical solutions for particular geometries have been reported by previous studies (Mason and Morrow 1991, Dong and Chatzis 1995, Lago and Araujo 2001), and semi-analytical solutions for cross-sections from 2D rock images have also been reported (Frette and Helland 2010). Otherwise, accurate measurement of pc can be made by LB bubble simulation (Ramstad et al. 2009). In this study, we introduce a shape analysis technique using Euclidean distance transform (EDT) to identify the effective pore cross-sectional area which is defined as the expected location of nonwetting fluid bubble during the invasion into capillary tube filled with wetting fluid. By extending MS-P theory, we calculate capillary pressures of tubes whose crosssections are acquired from 3D X-ray CT image of Berea sandstone compared with lattice Boltzmann simulation results under fully wetting condition (θ = 0) to represent a water-wet mineral surface (Jang and Santamarina 2011). Additionally, its application to pore network drainage simulation is made and its results before and after capillary pressure calibration is discussed.
- 855 -
Proceedings of the 19th International Conference on Soil Mechanics and Geotechnical Engineering, Seoul 2017
2 MATERIALS AND METHODS
2 .3
2 .1
As a pre-processing stage, synthetic tubes with constant crosssections, which are acquired during pore network extraction, are constructed as a simulation domain. Then, measurement of accurate pc of acquired pore throat cross-section is made by LB bubble simulation (Ramstad et al. 2009). Modified RothmanKeller model for D3Q19 lattice proposed by Liu et al. (2012) is used. Simulation starts with locating the nonwetting phase bubble at the center of the tube while the rest of the tube is filled with wetting phase fluid. Without any external forces, two phase fluid system would reach the equilibrium state after sufficiently long time steps. Then pc is determined by measuring the pressure difference between inside and outside of the equilibriated bubble as shown in Figure 3.
3D X-ray CT imaging
Berea sandstone sample (Cleveland Ltd., permeability: 100-200 mD, porosity: ~0.2) is subjected to 3D X-ray computed tomography for the quantitative evaluation of its pore structure. The resolution of reconstructed 3D image is 1.9 μm/pixel, and has 400×400×400 voxels. The 16-bit grayscaled images are binarized to extract the pore space by Otsu’s method (Otsu 1979). The porosity of binarized pore domain is 0.2193 which is calculated as the number of pore voxels over the total number of voxels. 2 .2
Pore network extraction
Lattice Boltzmann bubble simulation
Pore network extraction from 3D image aims to replicate the geometrical and topological features of real porous media. Pore network used in this study are composed of spherical pore chambers connected to cylindrical pore throats. Pore throats are assumed to have negligible volume compared with pore chambers. Figure 1 shows the schematic of pore network extraction process which will be discussed below.
Figure 3. Schematic of capillary pressure measurement during LB bubble simulation Figure 1. Schematic of pore network extraction process
Discretized pore space is initially subjected to medial axis transform to obtain pore skeleton. Junctions and dead ends of pore skeleton is defined as potential pore chambers. Then, potential pore chambers are dilated until their radii reach their Euclidean distance values. Remaining segments of pore skeleton which connect two different potential chambers are assigned as pore channels. Single element voxel among the pore skeleton which has minimum Euclidean distance (i.e., center of minimum inscribing sphere) is defined as a center of pore throat. Throat radius rt (which will be calibrated) is initially defined as the radius of minimum inscribing sphere rins. For the calibration, 2D cross-sectional images of pore throats are acquired by inserting the plane through the center of pore throat which is perpendicular to the pore channel skeleton. Plane insertion is performed by cubic interpolation to have 20 pixels per rins to achieve sufficient resolution for LB bubble simulation and pore shape analysis. The size of the plane is set to be 200×200 to capture the irregularity of the pore cross-section. Finally, potential pore throats are again dilated until they touch the inserted planes. The number of voxels consisting pore chamber are assigned as the volume of pore chamber. Figure 2 illustrates the binary 3D X-ray CT image of Berea sandstone, segmented pore space and extracted equivalent pore network used in this study. Number of pore chambers is 1193 and number of pore throats is 1719, and their cross-sectional images are subjected to shape analysis and LB simulation.
Figure 2. Binary 3D X-ray CT image of Berea sandstone, its segmented pore space and the equivalent pore network.
As a validation process, measured pc of regular polygon tubes are compared with analytical solutions by MS-P theory. For the further calibration, we introduce a concept of the normalized capillary pressure () which is defined as the ratio between measured or analytically solved pc and pc,ins which is the upper bound of the capillary pressure. Based on the analytical expression of pc for n-sided regular polygon by Lago and Araujo (2001), normalized pc based on MS-P theory under fully wetting condition can be expressed as follows:
pc pc ,ins
Ts rins
cot 1 n n 1 cot 1 2Ts 2 n n rins
(2)
Figure 4 shows the computed normalized capillary pressure by LB bubble simulation compared with analytical solutions expressed in Eq. (2). Good agreement between two confirms the validity of the measured pc of constructed synthetic tubes which will be discussed.
Figure 4. Comparison between measured by LB bubble simulation and analytically driven by MS-P theory.
- 856 -
Technical Committee 103 / Comité technique 103
2 .4
Pore throat shape analysis
Calibration of capillary pressure can be done by identifying the effective cross-section of the pore throat. The procedure starts with identifying the lower bound of the radius of arc meniscus (AM). General form of Young-Laplace equation under fully wetting condition gives: 1 1 pc Ts r r MTM AM
relationship between theoretical values of Geff and highlight that the pore shape analysis technique is valid and can be used for the correction of capillary pressure.
(3)
where rAM is the radius of AM and rMTM is the radius of main terminal meniscus (MTM). Note that pc could be written only in terms of rAM when sufficiently far away from the fluid interface (Lago and Araujo 2001, Frette and Helland 2010). Then, we can approxiate the effective radius of curvature reff which is the lower bound of rAM as follows: 1 1 Ts 2T lim Ts pc ,ins s r r r rins MTM AM eff
rMTM
(4)
By approximating reff as the half of rins, identification of effective cross-sectional area can be made. Let p = (p1, p2, …, pn) the pixels of the pore region of the pore throat crosssectional image and as s = (s1, s2, …, sm) the solid region. Obtained pore throat cross-sectional image I could be written as I = p∪s. Then, determine the region R1 where Euclidean distance values from the solid region are larger than the lower bound of rAM. Determine another region R2, the pore region where Euclidean distance values from R1 are smaller than the lower bound of rAM. Then, the identified region of effective pore throat cross-section could be written as R1∪R2. The schematic process is shown in Figure 5.
Figure 5. Schematic process of image-based effective pore cross-section identification
Obtaining area Aeff and perimeter Peff of effective cross-section allows the computation of pc thorugh MS-P theory. We define the effective shape factor Geff, which shares the same concept with . Ts Peff Geff
Aeff rins Peff pc 2Ts pc ,ins 2 Aeff rins
Figure 6. Correlation between normalized capillary pressure () and effective shape factor (Geff)
Note that pc values of pore throats are deterministically calculated by Eq. (1) when equivalent pore network is extracted. The correlation shown in Figure 6 and identified Geff values of the extracted pore network allows the correction of capillary pressures as follows:
pc Geff pc ,ins
(6)
By adopting Eq. (6), quasi-static drainage simulation is conducted to obtain water retention curve of Berea sandstone. We adopted the simulation method proposed by Joekar-Niasar et al. (2008). Initially, pore network is assumed to be saturated with wetting phase fluid. Pore chambers connected to the upper plane of the cubic domain are set to be wetting phase reservoir and its pressure is assumed to be zero and not changing. The pressure of wetting phase in all pores connected to the wetting phase reservoir is zero. A plane parallel to the wetting phase reservoir is set to be nonwetting phase reservoir and four perpendicular planes are set to be impermeable. Drainage simulation starts by increasing the pressure of nonwetting phase reservoir (dP). Nonwetting fluid invades when pressure difference at two pore chambers of different phase is higher than the imposed pc. Thus, we could obtain the static fluid configurations in the network at every dP stage. Figure 7 shows the changes in pore throat distribution before and after calibration.
(5)
3 RESULTS AND DISCUSSIONS Figure 6 shows the correlation between Geff and with superimposed effective cross-sections of acquired pore throats. Note that total 1719 pore throats are both subjected to LB bubble simulation and shape analysis. White region of the superimposed pore throat cross-section indicates the estimated area occupied by nonwetting fluid and blue region indicates the area occupied by wetting fluid. High R2 values and perfect 1:1
Figure 7. Pore throat distribution of Berea sandstone before and after capillary pressure calibration
- 857 -
Proceedings of the 19th International Conference on Soil Mechanics and Geotechnical Engineering, Seoul 2017
Both pore throat distributions follow log-normal distributions corroborated by previous study (Phandis and Santamarina 2011). Calibrated pore throat radius (rt,c) is determined by rt,c = G-1eff ·rt. Note that rt is initially defined by the radius of minimum inscribing sphere rins. Figure 8 shows the results of pore network drainage simulations before and after capillary pressure correction. Solid symbols indicate the simulation results fitted by van Genuchten equation (van Genuchten 1980). Result highlight that capillary pressure correction helps to overcome the limitation of pore network extraction process, which overestimates the capillary pressures of the pore throats.
Figure 8. Water retention curves of Berea sandstone before and after capillary pressure correction
4 CONCLUSIONS 3D X-ray CT provides sufficient spatial resolution to capture the pore structure of Berea sandstone. To overcome the capillary-pressure-overestimating limitation of the pore network extraction procedure, a simple shape analysis technique for computation of capillary pressures in irregularly shaped pore geometry is developed. The validity of the technique is supported by lattice Boltzmann bubble simulation results which is again validated by analytical solutions obtained by MS-P theory. Result highlight that simple pore shape analysis helps to determine more realistic capillary pressure of irregularly shaped pore channel and can be applied to pore network simulation to obtain realistic water retention curve.
5 ACKNOWLEDGMENTS This work was supported by the Korea CCS R&D Center (KCRC) grant and the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No. 2012-0008929, 2011-0030040, 2016R1A2B4011292).
6
REFERENCES
Al-Menhali A. S. and Krevor S. 2016. Capillary Trapping of CO2 in Oil Reservoirs: Observations in a Mixed-Wet Carbonate Rock. Environmental Science and Technology, 50 (5), 2727–2734. Brooks R. H. and Corey A. T. 1964. Hydraulic properties of porous media and their relation to drainage design. Transactions of the ASAE, 7 (1), 26–0228. Chen L. Kang Q. Robinson B. A. He Y. L. and Tao W. Q. 2013. Porescale modeling of multiphase reactive transport with phase
transitions and dissolution-precipitation processes in closed systems. Physical Review E, 87 (4), 1–16. Dai S. and Seol Y. 2014. Water permeability in hydrate-bearing sediments: A pore-scale study. Geophysical Research Letters, 41 (12), 4176–4184. Dong H. and Blunt M. J. 2009. Pore-network extraction from microcomputerized-tomography images. Physical Review E, 80 (3), 36307. Dong M. and Chatzis I. 1995. The imbibition and flow of a wetting liquid along the corners of a square capillary tube. Journal of Colloid and Interface Science, 172 (2), 278–288. Frette O. I. and Helland J. O. 2010. A semi-analytical model for computation of capillary entry pressures and fluid configurations in uniformly-wet pore spaces from 2D rock images. Advances in Water Resources, 33 (8), 846–866. Jang J. and Santamarina J. C. 2011. Recoverable gas from hydratebearing sediments: Pore network model simulation and macroscale analyses. Journal of Geophysical Research: Solid Earth, 116 (8), 1–12. Jang J. and Santamarina J. C. 2014. Evolution of gas saturation and relative permeability during gas production from hydrate-bearing sediments: Gas invasion vs. gas nucleation. Journal of Geophysical Research: Solid Earth, 119 (1), 116–126. Joekar-Niasar V. and Hassanizadeh S. M. 2012. Analysis of Fundamentals of Two-Phase Flow in Porous Media Using Dynamic Pore-Network Models: A Review. Critical Reviews in Environmental Science and Technology, 42 (18), 1895–1976. Joekar-Niasar V. Hassanizadeh S. M. and Leijnse A. 2008. Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using pore-network modeling. Transport in Porous Media, 74 (2), 201–219. Krevor S. C. M. Pini R. Li B. and Benson S. M. 2011. Capillary heterogeneity trapping of CO2 in a sandstone rock at reservoir conditions. Geophysical Research Letters, 38 (15), 1–5. Lago M. and Araujo M. 2001. Threshold Pressure in Capillaries with Polygonal Cross Section. Journal of Colloid and Interface Science, 243 (1), 219–226. Liu H. Valocchi A. J. and Kang Q. 2012. Three-dimensional lattice Boltzmann model for immiscible two-phase flow simulations. Physical Review E, 85 (4), 1–14. Liu M. Meakin P. and Huang H. 2006. Dissipative particle dynamics with attractive and repulsive particle-particle interactions. Physics of Fluids, 18 (1). Mason G. and Morrow N. R. 1991. Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. Journal of Colloid and Interface Science, 141 (1), 262–274. Mayer R. P. and Stowe R. 1965. Mercury porosimetry—breakthrough pressure for penetration between packed spheres. Journal of Colloid Science, 20 (8), 893–911. Phadnis H. S. and Santamarina J. C. 2011. Bacteria in sediments: pore size effects. Géotechnique Letters, 1 (4), 91–93. Princen H. M. 1969. Capillary phenomena in assemblies of parallel cylinders: II. Capillary rise between two cylinders. Journal of Colloid and Interface Science, 30 (3), 359–371. Princen H. M. 1969. Capillary phenomena in assemblies of parallel cylinders. I. Capillary rise between two cylinders. Journal of Colloid And Interface Science, 30 (1), 69–75. Princen H. M. 1970. Capillary phenomena in assemblies of parallel cylinders. III. Liquid Columns between Horizontal Parallel Cylinders. Journal of Colloid And Interface Science, 34 (2), 171– 184. Ramstad T. Oren P. E. and Bakke S. 2009. Simulation of Two Phase Flow in Reservoir Rocks Using a Lattice Boltzmann Method. SPE Annual Technical Conference and Exhibition. Silin D. and Patzek T. 2006. Pore space morphology analysis using maximal inscribed spheres. Physica A, 371 (2), 336–360. van Genuchten M. T. 1980. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1. Soil Science Society of America Journal.
- 858 -
ABSTRACT: Pore network modeling (PNM) has been useful to simulate water retention process of the capillary dominated systems. Randomly structured 3D pore geometry can be simplified as the network of pore channels of cylindrical shape. Radii of cylindrical pore channels are determined by minimum inscribing spheres which lead to the overestimation of capillary pressures of irregularly shaped pore throats. The correct capillary pressure can be obtained by bubble analysis using lattice Boltzmann method (LBM). Pore domain larger than its representative elementary volume, however, contains more than hundreds of pore throats. LB bubble simulations for estimating corrected capillary pressures of every pore throats are computationally demanding. Therefore, this study presents a simple shape analysis technique by Euclidean distance transform (EDT) to identify the effective cross-sectional area which is defined as the expected location of the nonwetting fluid during the invasion into a capillary tube filled with wetting fluid. By extending MS-P theory, we calculate capillary pressures of tubes whose effective cross-sections are taken from 3D X-ray CT image compared with LB bubble simulation. We also investigate the applicability of the shape analysis technique by calibrating the pore network drainage simulation results. RÉSUMÉ: La modélisation du réseau de pores (PNM) a été utile pour simuler le processus de rétention d'eau des systèmes capillaires dominés. La géométrie des pores 3D à structure aléatoire peut être simplifiée en tant que réseau de canaux poreux de forme cylindrique. Les rayons des canaux de pores cylindriques sont déterminés par des sphères d'inscription minimales qui conduisent à une surestimation des pressions capillaires de gorges de pores de forme irrégulière. La pression capillaire correcte peut être obtenue par analyse de bulles en utilisant la méthode de Boltzmann en treillis (LBM). Un domaine de pores plus grand que son volume élémentaire représentatif, cependant, contient plus de centaines de gorges de pores. Les simulations de bulles LB pour estimer les pressions capillaires corrigées de chaque gorge de pore sont très exigeantes en termes de calcul. Par conséquent, cette étude présente une technique d'analyse de forme simple par la transformée de distance euclidienne (EDT) pour identifier la section transversale efficace qui est définie comme l'emplacement attendu du fluide non mouillant pendant l'invasion dans un tube capillaire rempli de fluide mouillant. En prolongeant la théorie MS-P, nous calculons les pressions capillaires de tubes dont les coupes efficaces sont prises à partir de l'image de tomodensitométrie 3D par rapport à la simulation de bulles LB. Nous étudions également l 'applicabilité de la technique d' analyse de forme en calibrant les résultats de simulation de drainage de réseau de pores. KEYWORDS: Capillary pressure, Pore throat geometry, Pore network model, Lattice Boltzmann method 1 INTRODUCTION Carbon dioxide storage (Al-Menhali and Krevor 2016) or gas production from hydrate-bearing sediments (Dai and Seol 2014, Jang and Santamarina 2014) involves capillary-driven two phase fluid flow. Such systems include functional expressions for capillary pressure pc and wetting fluid saturation Sw as the governing equations (Brooks and Corey 1964, van Genuchten 1980). Efforts has been made to develop numerical models to capture pc-Sw relationship as an alternative of time consuming laboratory experiments (Frette and Helland 2010). It is well known that pore-scale simulators such as lattice Boltzmann models (LBM), (Liu et al. 2012) are computationally expensive compared to pore network models (PNM) due to the discretization methods (Joekar-Niasar and Hassanizadeh 2012). Advances in 3D imaging technology allows the direct extraction of pore network from 3D image, which is composed of pore chambers and pore throats, however, pore throats have been usually modeled as cylindrical capillary tubes with radius of minimum inscribing sphere rins (Silin and Patzek 2006) which leads to an overestimation of capillary entry pressure by adopting the Young-Laplace equation as follows: pc ,ins
2Ts cos rins
(1)
where Ts and θ indicate the interfacial tension and contact angle, respectively. Entry pressure for fluid invasion in tube having constant crosssection along the longitudinal direction is usually calculated by adopting Mayer and Stowe-Princen theory (MS-P theory) which is based on the energy balance (Mayer and Stowe 1965, Princen 1969a, 1969b, 1970). Analytical solutions for particular geometries have been reported by previous studies (Mason and Morrow 1991, Dong and Chatzis 1995, Lago and Araujo 2001), and semi-analytical solutions for cross-sections from 2D rock images have also been reported (Frette and Helland 2010). Otherwise, accurate measurement of pc can be made by LB bubble simulation (Ramstad et al. 2009). In this study, we introduce a shape analysis technique using Euclidean distance transform (EDT) to identify the effective pore cross-sectional area which is defined as the expected location of nonwetting fluid bubble during the invasion into capillary tube filled with wetting fluid. By extending MS-P theory, we calculate capillary pressures of tubes whose crosssections are acquired from 3D X-ray CT image of Berea sandstone compared with lattice Boltzmann simulation results under fully wetting condition (θ = 0) to represent a water-wet mineral surface (Jang and Santamarina 2011). Additionally, its application to pore network drainage simulation is made and its results before and after capillary pressure calibration is discussed.
- 855 -
Proceedings of the 19th International Conference on Soil Mechanics and Geotechnical Engineering, Seoul 2017
2 MATERIALS AND METHODS
2 .3
2 .1
As a pre-processing stage, synthetic tubes with constant crosssections, which are acquired during pore network extraction, are constructed as a simulation domain. Then, measurement of accurate pc of acquired pore throat cross-section is made by LB bubble simulation (Ramstad et al. 2009). Modified RothmanKeller model for D3Q19 lattice proposed by Liu et al. (2012) is used. Simulation starts with locating the nonwetting phase bubble at the center of the tube while the rest of the tube is filled with wetting phase fluid. Without any external forces, two phase fluid system would reach the equilibrium state after sufficiently long time steps. Then pc is determined by measuring the pressure difference between inside and outside of the equilibriated bubble as shown in Figure 3.
3D X-ray CT imaging
Berea sandstone sample (Cleveland Ltd., permeability: 100-200 mD, porosity: ~0.2) is subjected to 3D X-ray computed tomography for the quantitative evaluation of its pore structure. The resolution of reconstructed 3D image is 1.9 μm/pixel, and has 400×400×400 voxels. The 16-bit grayscaled images are binarized to extract the pore space by Otsu’s method (Otsu 1979). The porosity of binarized pore domain is 0.2193 which is calculated as the number of pore voxels over the total number of voxels. 2 .2
Pore network extraction
Lattice Boltzmann bubble simulation
Pore network extraction from 3D image aims to replicate the geometrical and topological features of real porous media. Pore network used in this study are composed of spherical pore chambers connected to cylindrical pore throats. Pore throats are assumed to have negligible volume compared with pore chambers. Figure 1 shows the schematic of pore network extraction process which will be discussed below.
Figure 3. Schematic of capillary pressure measurement during LB bubble simulation Figure 1. Schematic of pore network extraction process
Discretized pore space is initially subjected to medial axis transform to obtain pore skeleton. Junctions and dead ends of pore skeleton is defined as potential pore chambers. Then, potential pore chambers are dilated until their radii reach their Euclidean distance values. Remaining segments of pore skeleton which connect two different potential chambers are assigned as pore channels. Single element voxel among the pore skeleton which has minimum Euclidean distance (i.e., center of minimum inscribing sphere) is defined as a center of pore throat. Throat radius rt (which will be calibrated) is initially defined as the radius of minimum inscribing sphere rins. For the calibration, 2D cross-sectional images of pore throats are acquired by inserting the plane through the center of pore throat which is perpendicular to the pore channel skeleton. Plane insertion is performed by cubic interpolation to have 20 pixels per rins to achieve sufficient resolution for LB bubble simulation and pore shape analysis. The size of the plane is set to be 200×200 to capture the irregularity of the pore cross-section. Finally, potential pore throats are again dilated until they touch the inserted planes. The number of voxels consisting pore chamber are assigned as the volume of pore chamber. Figure 2 illustrates the binary 3D X-ray CT image of Berea sandstone, segmented pore space and extracted equivalent pore network used in this study. Number of pore chambers is 1193 and number of pore throats is 1719, and their cross-sectional images are subjected to shape analysis and LB simulation.
Figure 2. Binary 3D X-ray CT image of Berea sandstone, its segmented pore space and the equivalent pore network.
As a validation process, measured pc of regular polygon tubes are compared with analytical solutions by MS-P theory. For the further calibration, we introduce a concept of the normalized capillary pressure () which is defined as the ratio between measured or analytically solved pc and pc,ins which is the upper bound of the capillary pressure. Based on the analytical expression of pc for n-sided regular polygon by Lago and Araujo (2001), normalized pc based on MS-P theory under fully wetting condition can be expressed as follows:
pc pc ,ins
Ts rins
cot 1 n n 1 cot 1 2Ts 2 n n rins
(2)
Figure 4 shows the computed normalized capillary pressure by LB bubble simulation compared with analytical solutions expressed in Eq. (2). Good agreement between two confirms the validity of the measured pc of constructed synthetic tubes which will be discussed.
Figure 4. Comparison between measured by LB bubble simulation and analytically driven by MS-P theory.
- 856 -
Technical Committee 103 / Comité technique 103
2 .4
Pore throat shape analysis
Calibration of capillary pressure can be done by identifying the effective cross-section of the pore throat. The procedure starts with identifying the lower bound of the radius of arc meniscus (AM). General form of Young-Laplace equation under fully wetting condition gives: 1 1 pc Ts r r MTM AM
relationship between theoretical values of Geff and highlight that the pore shape analysis technique is valid and can be used for the correction of capillary pressure.
(3)
where rAM is the radius of AM and rMTM is the radius of main terminal meniscus (MTM). Note that pc could be written only in terms of rAM when sufficiently far away from the fluid interface (Lago and Araujo 2001, Frette and Helland 2010). Then, we can approxiate the effective radius of curvature reff which is the lower bound of rAM as follows: 1 1 Ts 2T lim Ts pc ,ins s r r r rins MTM AM eff
rMTM
(4)
By approximating reff as the half of rins, identification of effective cross-sectional area can be made. Let p = (p1, p2, …, pn) the pixels of the pore region of the pore throat crosssectional image and as s = (s1, s2, …, sm) the solid region. Obtained pore throat cross-sectional image I could be written as I = p∪s. Then, determine the region R1 where Euclidean distance values from the solid region are larger than the lower bound of rAM. Determine another region R2, the pore region where Euclidean distance values from R1 are smaller than the lower bound of rAM. Then, the identified region of effective pore throat cross-section could be written as R1∪R2. The schematic process is shown in Figure 5.
Figure 5. Schematic process of image-based effective pore cross-section identification
Obtaining area Aeff and perimeter Peff of effective cross-section allows the computation of pc thorugh MS-P theory. We define the effective shape factor Geff, which shares the same concept with . Ts Peff Geff
Aeff rins Peff pc 2Ts pc ,ins 2 Aeff rins
Figure 6. Correlation between normalized capillary pressure () and effective shape factor (Geff)
Note that pc values of pore throats are deterministically calculated by Eq. (1) when equivalent pore network is extracted. The correlation shown in Figure 6 and identified Geff values of the extracted pore network allows the correction of capillary pressures as follows:
pc Geff pc ,ins
(6)
By adopting Eq. (6), quasi-static drainage simulation is conducted to obtain water retention curve of Berea sandstone. We adopted the simulation method proposed by Joekar-Niasar et al. (2008). Initially, pore network is assumed to be saturated with wetting phase fluid. Pore chambers connected to the upper plane of the cubic domain are set to be wetting phase reservoir and its pressure is assumed to be zero and not changing. The pressure of wetting phase in all pores connected to the wetting phase reservoir is zero. A plane parallel to the wetting phase reservoir is set to be nonwetting phase reservoir and four perpendicular planes are set to be impermeable. Drainage simulation starts by increasing the pressure of nonwetting phase reservoir (dP). Nonwetting fluid invades when pressure difference at two pore chambers of different phase is higher than the imposed pc. Thus, we could obtain the static fluid configurations in the network at every dP stage. Figure 7 shows the changes in pore throat distribution before and after calibration.
(5)
3 RESULTS AND DISCUSSIONS Figure 6 shows the correlation between Geff and with superimposed effective cross-sections of acquired pore throats. Note that total 1719 pore throats are both subjected to LB bubble simulation and shape analysis. White region of the superimposed pore throat cross-section indicates the estimated area occupied by nonwetting fluid and blue region indicates the area occupied by wetting fluid. High R2 values and perfect 1:1
Figure 7. Pore throat distribution of Berea sandstone before and after capillary pressure calibration
- 857 -
Proceedings of the 19th International Conference on Soil Mechanics and Geotechnical Engineering, Seoul 2017
Both pore throat distributions follow log-normal distributions corroborated by previous study (Phandis and Santamarina 2011). Calibrated pore throat radius (rt,c) is determined by rt,c = G-1eff ·rt. Note that rt is initially defined by the radius of minimum inscribing sphere rins. Figure 8 shows the results of pore network drainage simulations before and after capillary pressure correction. Solid symbols indicate the simulation results fitted by van Genuchten equation (van Genuchten 1980). Result highlight that capillary pressure correction helps to overcome the limitation of pore network extraction process, which overestimates the capillary pressures of the pore throats.
Figure 8. Water retention curves of Berea sandstone before and after capillary pressure correction
4 CONCLUSIONS 3D X-ray CT provides sufficient spatial resolution to capture the pore structure of Berea sandstone. To overcome the capillary-pressure-overestimating limitation of the pore network extraction procedure, a simple shape analysis technique for computation of capillary pressures in irregularly shaped pore geometry is developed. The validity of the technique is supported by lattice Boltzmann bubble simulation results which is again validated by analytical solutions obtained by MS-P theory. Result highlight that simple pore shape analysis helps to determine more realistic capillary pressure of irregularly shaped pore channel and can be applied to pore network simulation to obtain realistic water retention curve.
5 ACKNOWLEDGMENTS This work was supported by the Korea CCS R&D Center (KCRC) grant and the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No. 2012-0008929, 2011-0030040, 2016R1A2B4011292).
6
REFERENCES
Al-Menhali A. S. and Krevor S. 2016. Capillary Trapping of CO2 in Oil Reservoirs: Observations in a Mixed-Wet Carbonate Rock. Environmental Science and Technology, 50 (5), 2727–2734. Brooks R. H. and Corey A. T. 1964. Hydraulic properties of porous media and their relation to drainage design. Transactions of the ASAE, 7 (1), 26–0228. Chen L. Kang Q. Robinson B. A. He Y. L. and Tao W. Q. 2013. Porescale modeling of multiphase reactive transport with phase
transitions and dissolution-precipitation processes in closed systems. Physical Review E, 87 (4), 1–16. Dai S. and Seol Y. 2014. Water permeability in hydrate-bearing sediments: A pore-scale study. Geophysical Research Letters, 41 (12), 4176–4184. Dong H. and Blunt M. J. 2009. Pore-network extraction from microcomputerized-tomography images. Physical Review E, 80 (3), 36307. Dong M. and Chatzis I. 1995. The imbibition and flow of a wetting liquid along the corners of a square capillary tube. Journal of Colloid and Interface Science, 172 (2), 278–288. Frette O. I. and Helland J. O. 2010. A semi-analytical model for computation of capillary entry pressures and fluid configurations in uniformly-wet pore spaces from 2D rock images. Advances in Water Resources, 33 (8), 846–866. Jang J. and Santamarina J. C. 2011. Recoverable gas from hydratebearing sediments: Pore network model simulation and macroscale analyses. Journal of Geophysical Research: Solid Earth, 116 (8), 1–12. Jang J. and Santamarina J. C. 2014. Evolution of gas saturation and relative permeability during gas production from hydrate-bearing sediments: Gas invasion vs. gas nucleation. Journal of Geophysical Research: Solid Earth, 119 (1), 116–126. Joekar-Niasar V. and Hassanizadeh S. M. 2012. Analysis of Fundamentals of Two-Phase Flow in Porous Media Using Dynamic Pore-Network Models: A Review. Critical Reviews in Environmental Science and Technology, 42 (18), 1895–1976. Joekar-Niasar V. Hassanizadeh S. M. and Leijnse A. 2008. Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using pore-network modeling. Transport in Porous Media, 74 (2), 201–219. Krevor S. C. M. Pini R. Li B. and Benson S. M. 2011. Capillary heterogeneity trapping of CO2 in a sandstone rock at reservoir conditions. Geophysical Research Letters, 38 (15), 1–5. Lago M. and Araujo M. 2001. Threshold Pressure in Capillaries with Polygonal Cross Section. Journal of Colloid and Interface Science, 243 (1), 219–226. Liu H. Valocchi A. J. and Kang Q. 2012. Three-dimensional lattice Boltzmann model for immiscible two-phase flow simulations. Physical Review E, 85 (4), 1–14. Liu M. Meakin P. and Huang H. 2006. Dissipative particle dynamics with attractive and repulsive particle-particle interactions. Physics of Fluids, 18 (1). Mason G. and Morrow N. R. 1991. Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. Journal of Colloid and Interface Science, 141 (1), 262–274. Mayer R. P. and Stowe R. 1965. Mercury porosimetry—breakthrough pressure for penetration between packed spheres. Journal of Colloid Science, 20 (8), 893–911. Phadnis H. S. and Santamarina J. C. 2011. Bacteria in sediments: pore size effects. Géotechnique Letters, 1 (4), 91–93. Princen H. M. 1969. Capillary phenomena in assemblies of parallel cylinders: II. Capillary rise between two cylinders. Journal of Colloid and Interface Science, 30 (3), 359–371. Princen H. M. 1969. Capillary phenomena in assemblies of parallel cylinders. I. Capillary rise between two cylinders. Journal of Colloid And Interface Science, 30 (1), 69–75. Princen H. M. 1970. Capillary phenomena in assemblies of parallel cylinders. III. Liquid Columns between Horizontal Parallel Cylinders. Journal of Colloid And Interface Science, 34 (2), 171– 184. Ramstad T. Oren P. E. and Bakke S. 2009. Simulation of Two Phase Flow in Reservoir Rocks Using a Lattice Boltzmann Method. SPE Annual Technical Conference and Exhibition. Silin D. and Patzek T. 2006. Pore space morphology analysis using maximal inscribed spheres. Physica A, 371 (2), 336–360. van Genuchten M. T. 1980. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1. Soil Science Society of America Journal.
- 858 -
