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 Cadoret et al. (1995), recent studies like Chapman et al. (2021), Sun et al. (2022), and Wang S. et al. (2022a), 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 Lin et al. (2017) and Wang S. et al. (2022a):

We use the CT data from Sun et al. (2022) 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 S A i r 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.

The second step of the method is to mesh the fluid heterogeneities (Figure 1B). Lissa et al. (2021) 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 Cepeda et al. (2013) 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 S A i r . 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 A i r a r e a : A i r a r e a = S A i r x − d , y , z + S A i r x + d , y , z + S A i r x , y − d , z + S A i r x , y + d , z + S A i r x , y , z − d + S A i r x , y , z + d 6 (3) where, y , and z are the space coordinates, d is the fine mesh size, that is fixed. S A i r x , y , z denotes the air saturation at a given spatial coordinate x , y , z . For example, S A i r x − d , y , z signifies the air saturation at the location x − d , y , z , while S A i r x + d , y , z denotes the air saturation at the location x + d , y , z , and so forth. The function A i r a r e a 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: M e s h   s i z e = d , 0.05 ≤ A i r a r e a ≤ 0.95 C , A i r a r e a > 0.95   a n d   A i r a r e a < 0.05 (4) where d and C 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.

We use Biot’s equations (Biot, 1956a; 1956b; 1962; Rubino et al., 2009; 2016) in the frequency-space domain: − ω 2 ρ b − ρ f 2 ρ c ω u s − ∇ ∙ σ = ρ f ρ c ω ∇ P f (5) ∇ ∙ − 1 ρ c ω ∇ P f − ω 2 ρ f u s − ω 2 P f M = α ω 2 ∇ ∙ u s (6) where ∇ is the Hamiltonian operator, ω is the angle frequency, α = 1 − K d K g is Biot-Willis coefficient, K d is the drained bulk modulus, K g is the bulk modulus of the grain. The density of the saturated sample is: ρ b = 1 − ϕ ρ s + ϕ ρ f (7) where the ρ f and ρ s are the densities of fluid and grain, respectively. ϕ is the porosity. The complex density is: ρ c ω = τ ρ f ϕ + η i ^ ω κ (8) where τ is the tortuosity of the pore and can be estimated roughly by τ = 1 2 1 + 1 ϕ according to Berryman (1982) and Rubino et al. (2009). η is the fluid viscosity, κ is the permeability, and i ^ is the imaginary unit.

The displacement vector of the rock matrix is u s = u i s , and the corresponding strain tensor is defined as ε i j = 1 2 u i , j s + u j , i s , where i , j = 1,2,3 are Euclidean space dimensions.

The stress tensor σ is related to the displacement of the matrix and fluid pressure P f , and: σ i j = 2 μ ε i j + δ i j λ m ∇ ∙ u s − α P f (9) where the λ m = K d − 2 3 μ , and the μ is the shear modulus, δ i j = 1   i = j 0   i ≠ j is Kronecker (delta) tensor. The so-called pore-space modulus (Gurevich et al., 2009) is defined as: M = ϕ K f + α − ϕ K g − 1 (10)

For a biphasic saturated sample (air/water), the effective fluid bulk modulus K f , density ρ f and viscosity η are defined respectively as: K f = S A i r K A i r + 1 − S A i r K w − 1 (11) ρ f = ρ A i r S A i r + ρ w 1 − S A i r (12) η = η A i r η w η A i r 1 − S A i r (13) where K A i r and K w are the bulk modulus of water and air, respectively, S A i r is the air saturation. ρ A i r and ρ w are the densities of air and water, respectively. Eq. 13 for the mixed-fluid viscosity η follows the work of Teja and Rice (1981), where η A i r and η w are the viscosities of air and water, respectively.

The third step is the oscillatory relaxation test (e.g., Rubino et al., 2009; 2016; Chapman and Quintal, 2018; Santos et al., 2021).

For computing the P-wave modulus, the boundary conditions are defined as follows: i) an axial oscillation stress σ 33 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., u 1 = u 2 = 0 ; 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 M P and P-wave attenuation Q P − 1 are obtained using the following: M P = σ 33 / ε 33 (14) Q P − 1 = I m a g M P / R e a l M P (15) where ε 33 and σ 33 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 σ 32 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., u 3 = 0 ; 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 Q S − 1 are obtained using: μ = σ 32 / 2 ε 32 (16) Q S − 1 = I m a g M S / R e a l M S (17) where ε 32 and σ 32 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 Borgomano et al. (2019) and Sun et al. (2022) and shown in Table 1 and Table 2.

Finally, the complex bulk modulus ( K ) can be deduced from the complex P-wave modulus M P and shear modulus μ by: K = M P − 4 3 μ (18) Q K − 1 = I m a g K / R e a l K (19)

The method’s fourth step is to solve Biot’s equations numerically (Eq. 5 and Eq. 6). We adopt a hybrid method (Halimi Bin Ibrahim and Skote, 2013), i.e., Newton iteration method (NIM) and LU matrix factorization method (LUM). Eq. 5 and Eq. 6 are rewritten in the following form: L X = L 1 u 1 , u 2 , u 3 , P f L 2 u 1 , u 2 , u 3 , P f L 3 u 1 , u 2 , u 3 , P f L 4 u 1 , u 2 , u 3 , P f = 0 (20) where the variable X = u 1 u 2 u 3 P f , is composed of the solid displacement vector u s = [ u 1 , u 2 , u 3 and fluid pressure P f . Here we drop the superscript s for the solid displacement vector to leave a space for a new superscript i counting the iteration number. We use a hybrid method to solve Equation 20: the displacement vector [ u 1 , u 2 , u 3 is solved using the NIM method, and the fluid pressure P f is solved by the LUM method. The detailed steps for the hybrid method are as follows:

Step A: Set the initial condition, X i = 0 = u 1 i = 0 u 2 i = 0 u 3 i = 0 P f i = 0 = 0 ,

Do loop on iteration number i :

Step B: Update the calculated variable u 1 i ( u 1 0 =0 for i=0). u 1 i + 1 is achieved using the NIM method: The Jacobian of the linear equation L 1 with respect to the independent variable u 1 is J 1 = ∂ L 1 ∂ u 1 ; then u 1 i + 1 = u 1 i + ∆ u 1 , where the updated term is ∆ u 1 = − L 1 J 1 (Ben-Israel, 1966).

Step C: Update the calculated variable u 2 i . u 2 i + 1 is estimated using the NIM method: Take the updated u 1 i + 1 into L 2 , then calculate Jacobian J 2 = ∂ L 2 ∂ u 2 ; and finally u 2 i + 1 = u 2 i + ∆ u 2 where updated term ∆ u 2 = − L 2 J 2 .

Step D: Update the calculated variable u 3 i . u 3 i + 1 is also estimated by the NIM method: Take the updated u 1 i + 1 and u 2 i + 1 into L 3 , the corresponding Jacobian is J 3 = ∂ L 3 ∂ u 3 ; then u 3 i + 1 = u 3 i + ∆ u 3 , where the updated term ∆ u 3 = − L 3 J 3 .

Step E: Update the calculated variable P f i . Take u 1 i + 1 , u 2 i + 1 and u 3 i + 1 into L 4 , then P f i + 1 is obtained using the LU matrix decomposition method (e.g., Bartels and Golub, 1969; Abbasbandy et al., 2006).

Step F: X i + 1 = u 1 i + 1 u 2 i + 1 u 3 i + 1 P f i + 1 , if X i + 1 − X i < ϵ , end the loop; Otherwise, go back to step B and i = i + 1 . ϵ 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).

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 σ 33 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 (Kümpel, 1991), 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 ( ε 11 = ε 22 = 0 ), the ratio P f σ 33 in a representative elementary volume (REV) fully saturated with water is deduced as P f σ 33 = B K u K u + 4 μ / 3 (21) where B = 1 − K d K u α is the Skempton’s coefficient (Kümpel, 1991), μ is the shear modulus. K u is the undrained bulk modulus obtained by Biot-Gassmann’s equation (Gassmann, 1951): K u = K d + α 2 ϕ K f + α − ϕ K g − 1 (22)

Using the parameters given in Table 1, the Skempton’s coefficient B =0.38 and P f σ 33 = 0.38 ∗ 32.5 32.5 + 4 ∗ 15.2 / 3 = 0.23 , 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, K u ≈ K d (Eq. 22 using air bulk modulus for K f ) 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 ( K u ≈ K d ) as: ε 33 = σ 33 K d + 4 μ / 3 ≈ 2.3 × 10 − 6 (23) 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: ε 33 = σ 33 K u + 4 μ / 3 ≈ 1.8 × 10 − 6 (24)

As a result, the axial strain varies from 2.3×10⁻⁶ to 1.8×10⁻⁶ 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.

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.

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): K f L F = 1 − S K a i r + S K w a t − 1 (25)

Here, S represents the water saturation. S=0.88, K a i r =10^–4 GPa, and K w a t =2.25 GPa. This results in K f L F = 8x10⁻⁴ GPa, which is approximately equal to K a i r . By extending Gassmann’s theory (as per Eq. 22) with K f L F , 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, Hill’s law (1963) can be used to define an effective bulk modulus for the entire sample: K H F = S K u + 4 / 3 μ + 1 − S K d + 4 / 3 μ − 1 − 4 / 3 μ (26)

Here, K d represents the drained bulk modulus and K u is the undrained bulk modulus fully saturated with water (Mavko and Mukerji, 1998). This upper limit, referred to as the Gassmann-Hill limit, results in K H F = 31.3 GPa in our case, which is less than K u = 32.5 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.

In numerous laboratory experiments (e.g., Batzle et al., 2006; Adelinet et al., 2010; Mikhaltsevitch et al., 2015; Sun et al., 2018), 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 Chapman and Quintal (2018). 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 Sun et al., 2022. We chose the Indiana sample for this investigation due to its lack of micro-cracks, thus eliminating the squirt-flow mechanism (Borgomano et al., 2019). 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.

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 Sun et al. (2022) between experimental data and numerical simulation, which can be attributed to their use of a 2D model.