Abstract
Elastic wave attenuation in partially saturated porous rock is primarily due to wave-induced fluid flow, which arises from the contrast in compressibility between air and water and is influenced by the water distribution within the rock. We propose a method for constructing a numerical model that predicts mesoscopic dispersion and attenuation. Initially, we use fluid distribution data sourced from 3D X-ray Computed Tomography images to construct the numerical model, utilizing Biot’s poroelastic equations as the governing equations. Subsequently, we implement the finite element method to derive solutions for the numerical model. Our focus is centered on two key challenges: 1) reducing memory cost, and 2) efficiently handling element intersection during the meshing process. The solutions illustrate the evolution of fluid pressure distribution and the frequency-dependent advancement of the elastic moduli, coupled with their corresponding attenuation. Ultimately, we compare these numerical predictions with previously published experimental data from a study on partially saturated Indiana limestone. The considerable agreement between our numerical results and the experimental data confirms the validity of our method, which crucially incorporates the actual fluid distribution (captured from 3D CT images) as a vital input.
1 Introduction
Characterization of fluid distribution in a reservoir is essential in several scenarios, such as monitoring CO2 geological storage and gas and oil production exploration (; ). Seismic waves are known to be affected by fluid; therefore, it is a valuable tool for detecting in situ fluid properties (; ; ). At the mesoscopic scale, for a biphasic saturated rock, like gas and water, seismic waves induce a pore pressure gradient due to the difference in the fluid bulk moduli, causing diffusion between the different fluid phases and, thus, energy transfer (; ). This diffusion process causes attenuation and dispersion of seismic waves, known as patchy-flow or mesoscopic-wave-induced fluid flow (). The mesoscopic scale refers to heterogeneities in the fluid distribution and rock fabric (e.g., ; ; ; Zhao Luanxiao et al., 2021b) greater than the pore size but smaller than the wavelength. The mesoscopic scale serves as a crucial bridge between the microscopic and macroscopic levels, enabling the significant upscaling of properties from the pore level to a broader, macroscopic perspective. The effect of mesoscopic flow on the dispersion and attenuation has been reported in a lot of experiments (e.g., ; ; ; ; ; ; Zhao Liming et al., 2021a; ). If we focus on fluid heterogeneity at the mesoscale, many analytical and numerical models can quantitatively assess its effect. One classical analytical model is the White model (e.g., ; White et al., 1975; ; ). It assumed that fluid patches are composed of periodic layers or spheres. The layer’s thickness or the sphere’s radius, i.e., the so-called patchy size, determines the critical frequency for dispersion and attenuation. A second kind of analytical model is to assume a random distribution of the fluid (e.g., ; ; ; ; ; ; Zhang et al., 2022). It assumed that the fluid distribution is stochastic and characterized by a correlation length, which can be used to predict the critical frequency of dispersion and attenuation. However, the prediction of the correlation length, the key parameter in this model, is not straightforward. Another way to predict the effect of mesoscopic flow on dispersion/attenuation is to use numerical models. It usually takes the fluid distribution as an input, and uses the finite element method to obtain the solution (e.g., ; ; ; ; ). The numerical model is computationally expensive compared to the analytical model. However, there is no assumption regarding the fluid distribution, making it more widely applicable. Fluid distribution can be obtained, for instance, from CT scan techniques (e.g., ; ; Zhu et al., 2017; 2023; ; ). In a recent study, measured the velocity dispersion and attenuation in a biphasic saturated sandstone (water and CO2 gas) and obtained the 3D fluid distribution using CT images. Their results indicate that the majority of the gas is situated towards the end of the sample, resembling a two-layer fluid distribution. Thus, they used an effective 1D numerical model and did not have to consider the cost of a 3D numerical simulation. More recently, measured the velocity dispersion and attenuation in a partially saturated (air/water) Indiana limestone. also obtained the 3D fluid distribution using the micro-CT image, and used the finite-element method to predict dispersion and attenuation. However, their numerical simulation was conducted in a 2D space due to unresolved memory consumption issues within the 3D numerical simulation. In addition, they observed a discrepancy between the 2D simulations and the experimental data, which they attributed to the difference between a 2D and 3D numerical simulation. A method for computing dispersion and attenuation in fully saturated rocks was presented by to predict squirt flow using a 3D CT image as input. As focused on squirt flow, the simulation was done on a cube containing several cracks leading to a cube size of (∼300 μm3), using 0.8 TB of RAM. This approach works perfectly for a local prediction, as for squirt flow; however, it is infeasible for mesoscopic flow: i) a very fine mesh would be needed to represent the distribution and geometry of the two fluid phases, ii) the simulation should be done at a larger scale (∼cm3).
The study describes a new and detailed method for numerically predicting dispersion and attenuation due to mesoscopic flow using a 3D fluid distribution obtained by a micro-CT image as an input. The finite element method solves the frequency-domain Biot’s equations to predict the fluid diffusion process. We present a method to overcome the problems of the element intersections in meshing and memory cost in solving Biot’s equations. Finally, the 3D numerical predictions are compared and discussed with experimental data published by .
2 Methods
Our proposed method consists of five steps: 1) reconstruct the fluid distribution to make the numerical model; 2) mesh the numerical model; 3) apply Biot’s equations as the governing equations; 4) set the boundary condition; 5) solve the numerical solution using the finite element method. These steps are tested on an Indiana specimen (Figure 1A). This carbonate rock has a porosity of 10.8% and a permeability of 2x10−17 m2. The dispersion of elastic wave velocity under confining pressure was investigated under dry and water saturation by and under partial saturation (air/water) by . Additionally, X-ray images under dry, fully water-saturated, and partially saturated conditions were obtained by .
FIGURE 1
2.1 Fluid distribution reconstruction
The estimation of fluid distribution is the first step of the method. For an homogenous dry sample, the fluid distribution dominates the heterogeneity of the partial saturated sample. Following
, recent studies like
,
, and
, the 3D fluid distribution can be estimated using the CT gray image of dry, partially and fully saturated sample. The distribution of the fluid is obtained following two steps.
(i) The gray image of partially and fully saturated samples should be normalized, referring to the gray value of two reference materials, for example, aluminum and sleeve. The rescaled image for the fully water-saturated sample can be obtained following Eq. 1, which is adapted from and :
where
is the normalized gray value,
is the gray value of the raw image,
and
are the average gray value of the two reference materials measured during the scan of the water-saturated sample.
and
are the average gray value of the two reference materials measured during the scan of the dry sample.
(ii) The images for the water-saturated sample and partially saturated sample are rescaled referring to the dry sample using Eq. 1. Then, the air saturation is calculated as:
We use the CT data from obtained on an Indiana limestone partially saturated by the drainage method to calculate the fluid saturation distribution according to Eq. 1 and Eq. 2. The air saturation is shown in Figure 1B. In this sample, the global water saturation obtained by drainage is 88%. In Figure 1B, white zones correspond to full air saturation while blue zones to pure water saturation. Air patches are distributed over the entire sample with sizes in the range of 0.5 mm–10 mm.
2.2 Numerical model meshing
The second step of the method is to mesh the fluid heterogeneities (Figure 1B). converted the CT images into a surface format in AVIZO to create triangular elements on every surface between solids and pores and on the boundaries of the investigated volume. Then, they imported the mesh ‘*.stl’ in COMSOL Multiphysics. However, this procedure cannot be used for partial saturation, as shown in Figure 1B. Indeed, Figure 1B shows that the volume contains many air patches with complex geometries; in particular, the meshing process in AVIZO leads to too many intersections or overlap elements, which are difficult to remove.
We use a method presented by to overcome the limitation. This method was developed first for medical CT scan images and allows incorporating complex geometries with non-uniform material properties in COMSOL Multiphysics. It is a practical alternative, as no intersections or overlapping elements occur, and thus, it avoids the need for critical geometry simplifications that may compromise the accuracy of the simulation. In our case (Figure 1B), the non-uniform property is the fluid saturation . The mesh must be refined at the water-air interfaces. We thus define an adaptive mesh using the following steps: first, a uniform 3D cylinder is constructed according to the size of the CT image and then divided with a coarse mesh; afterward, we refine the mesh at the water-air interfaces using a function :where, and are the space coordinates, is the fine mesh size, that is fixed. denotes the air saturation at a given spatial coordinate . For example, signifies the air saturation at the location , while denotes the air saturation at the location , and so forth. The function is 1 for the air-saturated zones and 0 for the water-saturated zones and varies in the partially saturated zones. Finally, considering the 5% uncertainty in the fluid distribution, the mesh size is defined as:where d and are the fine and coarse mesh sizes, respectively. Figure 1E shows the mesh size in the YZ section (Figure 1D) for the 3D fluid distribution given in Figure 1B according to Eq. 4, with C=3.5 mm and d=0.35 mm. Finally, a tetrahedral mesh is created and shown in Figure 1F. As expected, the mesh is coarse in the pure water saturation zone and refined in the partially saturated zone.
2.3 Governing equations
We use Biot’s equations (; ; ; ; ) in the frequency-space domain:where is the Hamiltonian operator, is the angle frequency, is Biot-Willis coefficient, is the drained bulk modulus, is the bulk modulus of the grain. The density of the saturated sample is:where the and are the densities of fluid and grain, respectively. is the porosity. The complex density is:where is the tortuosity of the pore and can be estimated roughly by according to and . is the fluid viscosity, is the permeability, and is the imaginary unit.
The displacement vector of the rock matrix is , and the corresponding strain tensor is defined as , where are Euclidean space dimensions.
The stress tensor is related to the displacement of the matrix and fluid pressure , and:where the , and the is the shear modulus, is Kronecker (delta) tensor. The so-called pore-space modulus () is defined as:
For a biphasic saturated sample (air/water), the effective fluid bulk modulus , density and viscosity are defined respectively as:where and are the bulk modulus of water and air, respectively, is the air saturation. and are the densities of air and water, respectively. Eq. 13 for the mixed-fluid viscosity follows the work of , where and are the viscosities of air and water, respectively.
2.4 Oscillatory relaxation test
The third step is the oscillatory relaxation test (e.g., ; ; ; ).
For computing the P-wave modulus, the boundary conditions are defined as follows: i) an axial oscillation stress is loaded at the top boundary of the sample, with an amplitude of 0.1 MPa; ii) the vertical displacement at the bottom boundary of the sample is set to zero; iii) for the side boundaries, the normal displacement is set to zero, i.e., ; iv) all the boundaries are impermeable for the fluid (no flow across the boundaries). The initial conditions for displacements are set to zero. The P-wave modulus and P-wave attenuation are obtained using the following:where and are the axial strain and stress, respectively.
To compute the shear modulus (directly related to the S-velocity), the boundary conditions are defined as follows: i) an oscillation stress is loaded at the top boundary of the sample, with an amplitude of 0.1 MPa; ii) the horizontal displacement at the bottom boundary of the sample is set to zero; iii) for the side boundaries, the axial displacement is set to zero, i.e., ; iv) all the boundaries are impermeable for the fluid (no flow across the boundaries). The initial conditions for displacements are set to zero. The shear modulus and attenuation are obtained using:where and are the shear strain and stress along the y direction, respectively. For both oscillatory relaxation tests, the physical properties of the rock sample and fluids used are deduced from and and shown in Table 1 and Table 2.
TABLE 1
| Properties | Indiana |
|---|---|
| Porosity (%) | 10.8 |
| Permeability-(m2) | 2x10−17 |
| Drained Bulk modulus- (GPa) | 24 ( = 5 MPa) |
| Bulk modulus of grain- (GPa) | 77 |
| Undrained Bulk modulus - (GPa) | 32.5 ( = 5 MPa) |
| Shear modulus- (GPa) | 15.2 ( = 5 MPa) |
| Biot-Willis coefficient- | 0.688 |
| Skempton’s coefficient-B | 0.38 |
| Density (kg/m3) | 2,369.2 |
Rock properties and the elastic parameters for the numerical prediction. is the confining pressure.
TABLE 2
| Properties | Water | Air |
|---|---|---|
| Bulk modulus- (GPa) | 2.25 | 1x10−4 |
| Density (kg/m3) | 1,000 | 1 |
| Viscosity- (Pa*s) | 10–3 | 2x10−5 |
Fluid properties for the numerical prediction.
Finally, the complex bulk modulus () can be deduced from the complex P-wave modulus and shear modulus by:
2.5 Numerical solution
The method’s fourth step is to solve Biot’s equations numerically (Eq. 5 and Eq. 6). We adopt a hybrid method (), i.e., Newton iteration method (NIM) and LU matrix factorization method (LUM). Eq. 5 and Eq. 6 are rewritten in the following form:where the variable , is composed of the solid displacement vector [ and fluid pressure . Here we drop the superscript s for the solid displacement vector to leave a space for a new superscript counting the iteration number. We use a hybrid method to solve Equation 20: the displacement vector [ is solved using the NIM method, and the fluid pressure is solved by the LUM method. The detailed steps for the hybrid method are as follows:
Step A: Set the initial condition, ,
Do loop on iteration number:
Step B: Update the calculated variable ( =0 for i=0). is achieved using the NIM method: The Jacobian of the linear equation with respect to the independent variable is ; then , where the updated term is ().
Step C: Update the calculated variable . is estimated using the NIM method: Take the updated into , then calculate Jacobian ; and finally where updated term .
Step D: Update the calculated variable . is also estimated by the NIM method: Take the updated and into , the corresponding Jacobian is ; then , where the updated term .
Step E: Update the calculated variable . Take , and into , then is obtained using the LU matrix decomposition method (e.g., ; ).
Step F:, if , end the loop; Otherwise, go back to step B and . is a relative error, defined as 10–3, which is a measure of the error relative to the size of each solution component.
Figure 2 shows the evolution of the relative error for displacements and fluid pressure as a function of the iteration number during the solving process. After ten iterations, Eq. 20 is solved with a relative error below 10–3. In the case of the hybrid method, only one independent variable is considered in every step, thus reducing the memory cost and calculation time. For comparison, we solved Eq. 20 using the hybrid method and the classical LU decomposition method and show the results in Figure 3: With the hybrid method (see the black line in Figure 3A), the memory cost is divided by a factor of 3 in comparison with the conventional LU method (see the red line in Figure 3A), and the computation time is reduced by a factor of 5 (Figure 3B).
FIGURE 2
FIGURE 3
3 Results and discussion
3.1 Axial strain , fluid pressure and local bulk modulus evolution
Using the physical properties (Table 1; Table 2) and the patchy air-water distribution (Figure 1B), we conducted an oscillatory-compressibility test (Section 2.4) to calculate the strain and fluid pressure as a function of the frequency oscillation. The distribution of the i) pore fluid pressure normalized to axial stress and ii) axial strain are shown in Figure 4 and Figure 5, respectively, for different frequencies. It can be observed that pore pressure gradients take their highest values at the air-water interfaces with higher values for frequencies above 10 Hz. However, at the low frequency of 1Hz, air and water pressures are equilibrated to a very low value (Figure 4A). We can refer to this state as a ‘relaxed state’ under undrained conditions. Indeed, during an axial oscillation of 1Hz, the pressure of the water increases due to the Skempton effect (), but the frequency is sufficiently low to give time for water to diffuse in the air-saturated zone, as air is much more compressible than water. At the highest frequency of 1 kHz (Figure 4D), water is pressurized and has no time to flow in the air-saturated zone, i.e., the distribution of the overpressure (red color in Figure 4D) is close to the distribution of the water saturation zones. We refer to this state as an “unrelaxed state” under undrained conditions. To estimate the increase of water pressure at 1 kHz, we recall that under the P-wave boundary condition (), the ratio in a representative elementary volume (REV) fully saturated with water is deduced aswhere is the Skempton’s coefficient (), is the shear modulus. is the undrained bulk modulus obtained by Biot-Gassmann’s equation ():
FIGURE 4
FIGURE 5
Using the parameters given in Table 1, the Skempton’s coefficient B =0.38 and , which is consistent with values shown in the water-saturated zones (red color in Figure 4D). For the air-saturated zones, as the compressibility of air is large, (Eq. 22 using air bulk modulus for ) and no pressurization is expected in agreement with the blue color in Figure 4D. Figure 4B and Figure 4C give the results of the water pressurization at the intermediate frequencies of 10 Hz and 100 Hz and illustrate the evolution of the pore pressure gradient in the sample as the frequency increases.
The distribution of axial strain at different frequencies is given in Figure 5. At the low frequency of 1 Hz (Figure 5A), there is no pressurization of the pore fluid (Figure 5A), and the axial deformation is homogeneous and can be approximately estimated () as:which is consistent with the value predicted by the numerical simulation in Figure 5A. On the other hand, in the case of an unrelaxed state at a REV scale, the axial strain is approximately given by:
As a result, the axial strain varies from 2.3×10−6 to 1.8×10−6 with an increase in frequency. This signifies that the behavior of porous rock shifts from a relaxed to an unrelaxed regime at the Representative Elementary Volume (REV) scale. With the rising frequency, the spatial distribution of the axial strain undergoes changes, which align with the evolution of the pore pressure (see Figure 4; Figure 5). Specifically, areas saturated with water exhibit less deformation compared to those saturated by air. Additionally, as the frequency increases, the count of less deformable patches escalates.
Finally, we performed an oscillatory-shear test to assess the influence of frequency on pore pressure, shear strain, and shear attenuation. As anticipated, the numerical simulations indicate that oscillatory-shear stress does not induce fluid pressurization. Furthermore, shear strain is found to be independent of frequency, and there is no observable shear attenuation. This corroborates the foundational assumption in poroelasticity theory: the fluid has no effect on the shear modulus.
3.2 Global P-wave, bulk, and shear moduli
We determined the global P-wave modulus and its corresponding attenuation (represented by black curves in Figure 6A; Figure 6B), shear modulus and its corresponding attenuation (black curves in Figures 6C,D), and bulk modulus with its associated attenuation (black curves in Figure 6E; Figure 6F) by the entire specimen. These are collectively referred to as the global modulus. The four frequencies (1, 10, 100, and 1,000 Hz) highlighted in Figure 4 and Figure 5 are illustrated as red lines in Figure 6A and Figure 6B.
FIGURE 6
Our initial observation revealed that the shear modulus (see Figure 6C) is independent of frequency, and there is no associated attenuation (indicating no water effect) (refer to Figure 6D). On the other hand, the bulk modulus (represented by the black curves in Figure 6E) ranges from 24.2 GPa to 31.3 GPa. This span corresponds to the bulk modulus in both the relaxed and unrelaxed states under the undrained boundary condition (see Figure 6).
At lower frequencies, the fluid pressure has ample time to equilibrate, yielding a homogeneously mixed fluid. Consequently, the bulk modulus of this mixed fluid can be treated as a single-phase effective fluid bulk modulus by applying Wood’s law (1946):
Here, S represents the water saturation. S=0.88, =10–4 GPa, and =2.25 GPa. This results in 8x10−4 GPa, which is approximately equal to . By extending Gassmann’s theory (as per Eq. 22) with , we infer a bulk modulus (24.2 GPa) that closely mirrors the drained bulk modulus (24 GPa). This represents the minimum value for the bulk modulus, often referred to as the low-frequency limit or the Gassmann-Wood limit.
At higher frequencies, there is insufficient time for fluid flow and pressure equalization. Under these conditions, individual fluid phases are effectively isolated, allowing for the use of Eq. 22 to define an undrained bulk modulus for each region saturated by its respective fluid. Following this, can be used to define an effective bulk modulus for the entire sample:
Here, represents the drained bulk modulus and is the undrained bulk modulus fully saturated with water (). This upper limit, referred to as the Gassmann-Hill limit, results in 31.3 GPa in our case, which is less than GPa. Both the drained bulk modulus and the Gassmann-Hill limit are depicted as dashed grey lines in Figure 6E.
In a short summary, the bulk modulus of the entire sample escalates from the Gassmann-Wood limit to the Gassmann-Hill limit with increasing frequency. The dispersion is associated with an attenuation (represented by the black curve in Figure 6F) peaking at 0.075 at 20 Hz. Lastly, we plot the P-wave modulus and its corresponding attenuation as functions of frequency in Figures 6A,B. As the P-wave modulus is a linear amalgamation of the shear and bulk moduli, its behavior closely mirrors that of the bulk modulus. Specifically, the P-wave modulus increases from 44 GPa to 52 GPa, accompanied by a peak attenuation of 0.045 at 20 Hz.
3.3 Numerical prediction vs experimental measurements
In numerous laboratory experiments (e.g., ; ; ; ), researchers measure the strain’s evolution with frequency under oscillatory stress using local strain gauges. This approach yields a locally measured bulk modulus. To emulate such experiments, we average strain over a span of 6 mm—the typical length of a strain gauge—and simulate four strain gauges that are averaged. These gauges are situated at the sample’s half-length (represented by red lines in Figure 1C), a quarter-length from the bottom (blue lines in Figure 1C), and a quarter-length from the top (green lines in Figure 1C). The results are presented in Figure 6E: In general, the local bulk modulus and attenuation measured at the midpoint (50%, red curve) and three-quarters (75%, green curve) of the sample length closely align with the bulk properties of the entire sample (black curve). However, the frequency-dependent evolution of the local bulk modulus at one-quarter (25%, blue curve) of the sample deviates significantly from the evolution of the global bulk modulus (black curve). Furthermore, the high-frequency limit of local measurements differs from the Gassmann-Hill limit, as the strain gauges only capture the effects of local saturation. Consequently, high-frequency results may approach the undrained bulk modulus fully saturated with water, which surpasses the Gassmann-Hill limit. Interestingly, the disparity between local and global measurements in our study is not as marked as in . This reduced difference can likely be attributed to i) our approach of consolidating local results from four strain gauges, which serves to lessen the discrepancy between local and global responses, and ii) a comparatively uniform air/water distribution in our experiment, as opposed to the more varied distribution seen in Chapman and Quintal’s work (2018).
The progression of the moduli with frequency for the Indiana sample, saturated to 88% (Figures 1A,B), was studied by . We chose the Indiana sample for this investigation due to its lack of micro-cracks, thus eliminating the squirt-flow mechanism (). In these tests, strain measurements were taken by averaging readings from four strain gauges situated in the half-length of the specimen. We have represented the experimental data in Figure 6 with black square dots. Figures 6C,D juxtapose the projected and observed values of the shear modulus and attenuation. The frequency-independent shear modulus (depicted by the black curve in Figure 6C) aligns with the measurements (square dots in Figure 6C) spanning the seismic bands (0.1–100 Hz) and ultrasonic frequency (1 MHz). Additionally, the predicted shear attenuation (black curve in Figure 6D) aligns with the measurement (square dots in Figure 6D), given that the measurement error range for attenuation is within 0.02.
Figures 6E,F juxtapose the anticipated and actual measurements of the bulk modulus and attenuation. Generally, there’s an excellent agreement between the measurements and the numerical predictions, as illustrated by the black square dots and red curves. The numerical simulation, based on the CT-scan images, accurately replicates i) the dispersion of the modulus and ii) both the low and high-frequency limit. There’s also a strong correspondence between the measured attenuation and the numerical simulation. This comparison between the numerical simulation and experimental data confirms the efficacy of the method detailed in this paper in predicting mesoscopic dispersion and attenuation.
3.4 2D model vs 3D model
We now turn to comparing the results derived from 2D and 3D numerical simulations. Utilizing the YZ section of the CT image (Figure 1D), the bulk modulus and attenuation for the 2D model were computed. The global response of the 2D model is determined by averaging the strain across the entire section. As shown by the solid gray line in Figure 6F, the prediction using the 2D model exhibits a higher critical frequency (100 Hz) and peak attenuation (0.08) compared to the 3D model (black curve in Figure 6F). Furthermore, the bulk modulus derived from the 2D model (solid gray line in Figure 6E) deviates from the one obtained from the 3D model (black curve in Figure 6E). The discrepancy between the 2D and 3D results is attributed to the fluid flow in the ZX and ZY direction, which is not accounted for in the 2D numerical simulation. This finding elucidates the mismatch observed by between experimental data and numerical simulation, which can be attributed to their use of a 2D model.
4 Conclusion
In this paper, we introduce a novel method aimed at predicting velocity dispersion and attenuation attributable to mesoscopic flow, leveraging actual fluid distribution data derived from CT images. The numerical model is initially established by meshing CT images through a technique adept at handling intricate geometries and non-uniform material properties, thus effectively bypassing element intersections, a common issue associated with AVIZO. To find the solution for the numerical model governed by Biot’s equations, we use a hybrid method that significantly curtails memory cost compared to the LU matrix factorization method.
The solution from the numerical model forecasts the evolution of pore pressure distribution with frequency, thereby anticipating the advancement of the elastic moduli and their attenuation. We also model the development of the moduli for the entire sample (global moduli) as well as those measured by a strain gauge (local moduli). The discrepancies observed between local and global responses can be attributed to the heterogeneity in fluid distribution. Importantly, the 3D model’s predictions are validated by experimental data collected from Indiana limestone.
The presented method successfully addresses issues pertaining to memory consumption and calculation time, thereby setting the stage for quantifying the relationship between fluid distribution and seismic attenuation. This innovative approach holds the potential to serve as a robust tool for upscaling at the reservoir scale.
Statements
Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.
Author contributions
CS: Methodology, Writing–original draft, Writing–review and editing. JF: Methodology, Writing–original draft, Writing–review and editing. GT: Supervision, Writing–review and editing. SW: Supervision, Writing–review and editing.
Funding
The author(s) declare financial supportwas received for the research, authorship, and/or publication of this article. This work is supported by the National Natural Science Foundation of China (42104111, 42274142, 41930425, 41774143, 41804104) and also supported by Open Fund (WX-KFJJ-2022-08) of SINOPEC Key Laboratory of Geophysics, and State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, and Xuzhou Science and Technology Bureau Young Talents Project (No. KC22018). The authors declare that this study received funding from SINOPEC Key Laboratory of Geophysics. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.
Acknowledgments
We acknowledge the help of Dr. Jan V. M. Borgomano at Ecole Normale Supérieure for scanning the sample.
Conflict of interest
Authors GT and SW were employed by China National Petroleum Corporation.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
AbbasbandyS.EzzatiR.JafarianA. (2006). LU decomposition method for solving fuzzy system of linear equations. Appl. Math. Comput.172 (1), 633–643. 10.1016/j.amc.2005.02.018
2
AdelinetM.DorbathC.RavalecM. L.FortinJ.GuéguenY. (2011). Deriving microstructure and fluid state within the Icelandic crust from the inversion of tomography data. Geophys. Res. Lett.38 (3). 10.1029/2010GL046304
3
AdelinetM.FortinJ.GuéguenY.SchubnelA.GeoffroyL. (2010). Frequency and fluid effects on elastic properties of basalt: experimental investigations. Geophys. Res. Lett.37 (2). 10.1029/2009GL041660
4
AnwerH. M.AliA.AlvesT. M. (2017). Bayesian inversion of synthetic AVO data to assess fluid and shale content in sand-shale media. J. Earth Syst. Sci.126 (3), 42. 10.1007/s12040-017-0818-y
5
BaJ.CarcioneJ. M.SunW. (2015). Seismic attenuation due to heterogeneities of rock fabric and fluid distribution. Geophys. J. Int.202 (3), 1843–1847. 10.1093/gji/ggv255
6
BaJ.XuW.FuL.-Y.CarcioneJ. M.ZhangL. (2017). Rock anelasticity due to patchy saturation and fabric heterogeneity: A double double-porosity model of wave propagation. J. Geophys. Res. Solid Earth122 (3), 1949–1976. 10.1002/2016JB013882
7
BartelsR. H.GolubG. H. (1969). The simplex method of linear programming using LU decomposition. Commun. ACM12 (5), 266–268. 10.1145/362946.362974
8
BatzleM. L.HanD.-H.HofmannR. (2006). Fluid mobility and frequency-dependent seismic velocity — direct measurements. GEOPHYSICS71 (1), N1–N9. 10.1190/1.2159053
9
Ben-IsraelA. (1966). A Newton-Raphson method for the solution of systems of equations. J. Math. Analysis Appl.15 (2), 243–252. 10.1016/0022-247X(66)90115-6
10
Berryman (1982). “Elastic waves in fluid-saturated porous media (Vol 154),” in Macroscopic Properties of Disordered Media, Rome, June 1–3, 1981.
11
BiotM. A. (1962). Mechanics of deformation and acoustic propagation in porous media. J. Appl. Phys.33 (4), 1482–1498. 10.1063/1.1728759
12
BiotM. A. (1956b). Theory of propagation of elastic waves in a fluid‐saturated porous solid. II. Higher frequency range. J. Acoust. Soc. Am.28 (2), 179–191. 10.1121/1.1908241
13
BiotM. A. (1956a). Theory of propagation of elastic waves in a fluid‐saturated porous solid. I. Low‐Frequency range. J. Acoust. Soc. Am.28 (2), 168–178. 10.1121/1.1908239
14
BorgomanoJ. V. M.PimientaL. X.FortinJ.GuéguenY. (2019). Seismic dispersion and attenuation in fluid-saturated carbonate rocks: effect of microstructure and pressure. J. Geophys. Res. Solid Earth124 (12), 12498–12522. 10.1029/2019JB018434
15
ChapmanS.QuintalB. (2018). “Numerical assessment of local versus bulk strain measurements to quantify seismic attenuation in partially saturated rocks,” in SEG technical program expanded abstracts 2018 (Texas, United States: Society of Exploration Geophysicists), 1–0, 3547–3551. 10.1190/segam2018-2992202.1
16
CadoretT.MavkoG.ZinsznerB. (1998). Fluid distribution effect on sonic attenuation in partially saturated limestones. GEOPHYSICS63 (1), 154–160. 10.1190/1.1444308
17
CadoretT.MarionD.ZinsznerB. (1995). Influence of frequency and fluid distribution on elastic wave velocities in partially saturated limestones. J. Geophys. Res. Solid Earth100 (B6), 9789–9803. 10.1029/95JB00757
18
CavalliniF.CarcioneJ. M.Vidal de VentósD.Engell-SørensenL. (2017). Low-frequency dispersion and attenuation in anisotropic partially saturated rocks. Geophys. J. Int.209 (3), 1572–1584. 10.1093/gji/ggx107
19
CepedaJ. F.BirlaS.SubbiahJ.ThippareddiH. (2013). A practical method to model complex three-dimensional geometries with non-uniform material properties using image-based design and COMSOL multiphysics.
20
ChapmanS.BorgomanoJ. V. M.QuintalB.BensonS. M.FortinJ. (2021). Seismic wave attenuation and dispersion due to partial fluid saturation: direct measurements and numerical simulations based on x‐ray ct. J. Geophys. Res. Solid Earth126 (4). 10.1029/2021JB021643
21
ChapmanS.TisatoN.QuintalB.HolligerK. (2016). Seismic attenuation in partially saturated berea sandstone submitted to a range of confining pressures: seismic attenuation in berea sandstone. J. Geophys. Res. Solid Earth121 (3), 1664–1676. 10.1002/2015JB012575
22
DuttaN. C.OdéH. (1979). Attenuation and dispersion of compressional waves in fluid-filled porous rocks with partial gas saturation (White model); Part II, Results. Geophysics44 (11), 1789–1805. 10.1190/1.1440939
23
GassmannF. (1951). Über die elastizität poröser medien: vier. der natur. Ges. Zürich96, 1–23.
24
GurevichB.BrajanovskiM.GalvinR. J.MüllerT. M.Toms‐StewartJ. (2009). P-wave dispersion and attenuation in fractured and porous reservoirs – poroelasticity approach. Geophys. Prospect.57 (2), 225–237. 10.1111/j.1365-2478.2009.00785.x
25
Halimi Bin IbrahimI.SkoteM. (2013). Effects of the scalar parameters in the Suzen‐Huang model on plasma actuator characteristics. Int. J. Numer. Methods Heat Fluid Flow23 (6), 1076–1103. 10.1108/HFF-05-2011-0108
26
HeY.-X.WangS.YuanS.TangG.WuX. (2020). An improved approach for hydrocarbon detection using Bayesian inversion of frequency- and angle-dependent seismic signatures of highly attenuative reservoirs. IEEE Geoscience Remote Sens. Lett.19, 1–5. 10.1109/LGRS.2020.3017627
27
HillR. (1963). Elastic properties of reinforced solids: some theoretical principles. J. Mech. Phys. Solids11 (5), 357–372. 10.1016/0022-5096(63)90036-X
28
KlimentosT. (1995). Attenuation of P‐ and S‐waves as a method of distinguishing gas and condensate from oil and water. GEOPHYSICS60 (2), 447–458. 10.1190/1.1443782
29
KümpelH.-J. (1991). Poroelasticity: parameters reviewed. Geophys. J. Int.105 (3), 783–799. 10.1111/j.1365-246X.1991.tb00813.x
30
LinQ.BijeljicB.RaeiniA. Q.RiekeH.BluntM. J. (2021). Drainage capillary pressure distribution and fluid displacement in a heterogeneous laminated sandstone. Geophys. Res. Lett.48 (14), e2021GL093604. 10.1029/2021GL093604
31
LinQ.BijeljicB.RiekeH.BluntM. J. (2017). Visualization and quantification of capillary drainage in the pore space of laminated sandstone by a porous plate method using differential imaging X-ray microtomography: imaging of capillary drainage using dipp. Water Resour. Res.53 (8), 7457–7468. 10.1002/2017WR021083
32
LissaS.RufM.SteebH.QuintalB. (2021). Digital rock physics applied to squirt flow. GEOPHYSICS86 (4), MR235–MR245. 10.1190/geo2020-0731.1
33
MavkoG.MukerjiT. (1998). Bounds on low‐frequency seismic velocities in partially saturated rocks. GEOPHYSICS63 (3), 918–924. 10.1190/1.1444402
34
MikhaltsevitchV.LebedevM.GurevichB. (2016). Laboratory measurements of the effect of fluid saturation on elastic properties of carbonates at seismic frequencies: effect of fluid saturation on carbonates. Geophys. Prospect.64 (4), 799–809. 10.1111/1365-2478.12404
35
MikhaltsevitchV.LebedevM.GurevichB. (2015). A laboratory study of attenuation and dispersion effects in glycerol-saturated Berea sandstone at seismic frequencies. Texas, United States: Society of Exploration Geophysicists, 3085–3089. 10.1190/segam2015-5898429.1
36
MonachesiL. B.WollnerU.DvorkinJ. (2020). Effective pore fluid bulk modulus at patchy saturation: an analytic study. J. Geophys. Res. Solid Earth125 (1), e2019JB018267. 10.1029/2019JB018267
37
MüllerT. M.Toms‐StewartJ.WenzlauF. (2008). Velocity-saturation relation for partially saturated rocks with fractal pore fluid distribution. Geophys. Res. Lett.35 (9), L09306. 10.1029/2007GL033074
38
MüllerT. M.GurevichB.LebedevM. (2010). Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks A review. Geophysics75 (5), 75A147–75A164. 10.1190/1.3463417
39
MüllerT. M.GurevichB. (2004). One-dimensional random patchy saturation model for velocity and attenuation in porous rocks. Geophysics69 (5), 1166–1172. 10.1190/1.1801934
40
MüllerT. M.GurevichB. (2005). Wave-induced fluid flow in random porous media: attenuation and dispersion of elastic waves. J. Acoust. Soc. Am.117 (5), 2732–2741. 10.1121/1.1894792
41
PrideS. R.BerrymanJ. G.HarrisJ. M. (2004). Seismic attenuation due to wave-induced flow. J. Geophys. Res. Solid Earth109 (B1), B01201. 10.1029/2003JB002639
42
QiQ.MüllerT. M.GurevichB.LopesS.LebedevM.CaspariE. (2014). Quantifying the effect of capillarity on attenuation and dispersion in patchy-saturated rocks. GEOPHYSICS79 (5), WB35–WB50. 10.1190/geo2013-0425.1
43
QuintalB.FrehnerM.MadonnaC.TisatoN.KuteynikovaM.SaengerE. H. (2011). Integrated numerical and laboratory rock physics applied to seismic characterization of reservoir rocks. Lead. Edge30 (12), 1360–1367. 10.1190/1.3672480
44
RubinoJ. G.CaspariE.MüllerT. M.MilaniM.BarbosaN. D.HolligerK. (2016). Numerical upscaling in 2-D heterogeneous poroelastic rocks: anisotropic attenuation and dispersion of seismic waves. J. Geophys. Res. Solid Earth121 (9), 6698–6721. 10.1002/2016JB013165
45
RubinoJ. G.RavazzoliC. L.SantosJ. E. (2009). Equivalent viscoelastic solids for heterogeneous fluid-saturated porous rocks. GEOPHYSICS74 (1), N1–N13. 10.1190/1.3008544
46
SantosJ.CarcioneJ.BaJ. (2021). Two-phase flow effects on seismic wave anelasticity in anisotropic poroelastic media. Energies14 (20), 6528. 10.3390/en14206528
47
SantosJ. E.RavazzoliC. L.GauzellinoP. M.CarcioneJ. M. (2005). Numerical simulation of ultrasonic waves in reservoir rocks with patchy saturation and fractal petrophysical properties. Comput. Geosci.9 (1), 1–27. 10.1007/s10596-005-2848-9
48
SunC.FortinJ.BorgomanoJ. V. M.WangS.TangG.BultreysT.et al (2022). Influence of fluid distribution on seismic dispersion and attenuation in partially saturated limestone. J. Geophys. Res. Solid Earth127 (5). 10.1029/2021JB023867
49
SunC.TangG.ZhaoJ.ZhaoL.WangS. (2018). An enhanced broad-frequency-band apparatus for dynamic measurement of elastic moduli and Poisson’s ratio of rock samples. Rev. Sci. Instrum.89 (6), 064503. 10.1063/1.5018152
50
SunW. (2017). Determination of elastic moduli of composite medium containing bimaterial matrix and non-uniform inclusion concentrations. Appl. Math. Mech.38 (1), 15–28. 10.1007/s10483-017-2157-6
51
TejaA. S.RiceP. (1981). Generalized corresponding states method for the viscosities of liquid mixtures. Industrial Eng. Chem. Fundam.20 (1), 77–81. 10.1021/i100001a015
52
TesterJ. W.AndersonB. J.BatchelorA. S.BlackwellD. D.DiPippoR.DrakeE. M.et al (2007). Impact of enhanced geothermal systems on US energy supply in the twenty-first century. Philosophical Trans. R. Soc. A Math. Phys. Eng. Sci.365 (1853), 1057–1094. 10.1098/rsta.2006.1964
53
TisatoN.QuintalB. (2013). Measurements of seismic attenuation and transient fluid pressure in partially saturated berea sandstone: evidence of fluid flow on the mesoscopic scale. Geophys. J. Int.195 (1), 342–351. 10.1093/gji/ggt259
54
TomsJ.MüllerT. M.GurevichB. (2007). Seismic attenuation in porous rocks with random patchy saturation. Geophys. Prospect.55 (5), 671–678. 10.1111/j.1365-2478.2007.00644.x
55
Toms-StewartJ.MüllerT. M.GurevichB.PatersonL. (2009). Statistical characterization of gas-patch distributions in partially saturated rocks. GEOPHYSICS74 (2), WA51–WA64. 10.1190/1.3073007
56
WangS.RuspiniL. C.ØrenP.-E.Van OffenwertS.BultreysT. (2022a). Anchoring multi-scale models to micron-scale imaging of multiphase flow in rocks. Water Resour. Res.58 (1), e2021WR030870. 10.1029/2021WR030870
57
WangY.ZhaoL.CaoC.YaoQ.YangZ.CaoH.et al (2022b). Wave-induced fluid pressure diffusion and anelasticity in partially saturated rocks: the influences of boundary conditions. Geophysics87 (5), MR247–MR263. 10.1190/geo2021-0809.1
58
WhiteJ. E. (1975). Computed seismic speeds and attenuation in rocks with partial gas saturation. Geophysics40 (2), 224–232. 10.1190/1.1440520
59
WhiteJ. E.MihailovaN.LyakhovitskyF. (1975). Low‐frequency seismic waves in fluid‐saturated layered rocks. J. Acoust. Soc. Am.57 (S1), S30. 10.1121/1.1995164
60
WoodA. B. (1946). A textbook of sound. 2. New York: Macmillan. Available at: http://archive.org/details/in.ernet.dli.2015.15768.
61
ZhangL.BaJ.CarcioneJ. M.WuC. (2022). Seismic wave propagation in partially saturated rocks with a fractal distribution of fluid‐patch size. J. Geophys. Res. Solid Earth127 (2). 10.1029/2021JB023809
62
ZhaoL.TangG.SunC.ZhaoJ.WangS. (2021a). Dual attenuation peaks revealing mesoscopic and microscopic fluid flow in partially oil-saturated Fontainebleau sandstones. Geophys. J. Int.224 (3), 1670–1683. 10.1093/gji/ggaa551
63
ZhaoL.WangY.YaoQ.GengJ.LiH.YuanH.et al (2021b). Extended Gassmann equation with dynamic volumetric strain: modeling wave dispersion and attenuation of heterogeneous porous rocks. Geophysics86 (3), MR149–MR164. 10.1190/geo2020-0395.1
64
ZhuW.ZhaoL.ShanR. (2017). Modeling effective elastic properties of digital rocks using a new dynamic stress-strain simulation method. Geophysics82 (6), MR163–MR174. 10.1190/geo2016-0556.1
65
ZhuW.ZhaoL.YangZ.CaoH.WangY.ChenW.et al (2023). Stress relaxing simulation on digital rock: characterize attenuation due to wave‐induced fluid flow and scattering. J. Geophys. Res. Solid Earth128 (2), e2022JB024850. 10.1029/2022jb024850
Summary
Keywords
3D CT image, dispersion, attenuation, mesoscopic-flow, numerical modeling
Citation
Sun C, Fortin J, Tang G and Wang S (2023) Prediction of dispersion and attenuation on elastic wave velocities in partially saturated rock based on the fluid distribution obtained from three-dimensional (3D) micro-CT images. Front. Earth Sci. 11:1267522. doi: 10.3389/feart.2023.1267522
Received
26 July 2023
Accepted
14 August 2023
Published
24 August 2023
Volume
11 - 2023
Edited by
Qiaomu Qi, Chengdu University of Technology, China
Updates
Copyright
© 2023 Sun, Fortin, Tang and Wang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chao Sun, kang2008ping2008@163.com
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.