Minor Squirt in Unconsolidated Sands versus Strong Squirt in Compressed Glass Beads

Squirt driven by local pressure imbalance between contact of grains (or throat) and the main pore space is a mechanism of P-wave attenuation in consolidated rocks. In this paper, we investigate squirt in unconsolidated and consolidated porous media (represented by Toyoura sands and compressed glass beads, respectively). The former sample has very small bulk modulus and shear modulus, manifested by relatively free/mobile grains. As such, solid stress on the surface tends to be uniform and squirt should be minor. Biot’s theory improved with dynamic permeability successfully predicts the ultrasonic velocity and quality factor of P wave measured in the unconsolidated sands, confirming the aforementioned judgment. Dynamic permeability inverted at a high frequency is far lower than Darcy permeability. However, the improved model remains to be incapable of predicting velocity and attenuation measured in the latter sample; the reason is that the compressed beads allows a high pressure at compliant throat which drives strong squirt between throat and the main pore space.


Introduction
Biot [1,2] proposed a theory for waves in fluid saturated rocks, based on single porosity, elastic solid, and viscous fluid. The dilatational wave of the first kind in the Biot theory was fast P wave, while the dilatational wave of the second kind was slow P wave. Fast P wave has a high speed well observed in seismology [3], while slow P wave is diffusive in which fluid pressure drives unsteady flow in porous media [4,5].
The Biot [1,2] theory was found to yield inaccurate prediction of ultrasonic attenuation in consolidated rocks [6,7]. Therefore, Mavko and Nur [8], Budiansky and O'Connell [9], and Murphy et al. [10] proposed a mechanism of squirt due to local pressure difference. Dvorkin et al. [11,12] developed a BISQ model considering the squirt of fluid between compliant microcracks and the main pore space. For porous media, such microcracks are contact of grains or throat.
The initial BISQ model [11] was no longer considered to be valid as it at the low frequency limit was inconsistent with the Gassmann equation [13]; i.e., it violated fluid mass conservation. The BISQ model [12] was based on the idea that a cylindrical volume could represent the squirt region (called modified solid) and the skeleton was modified accordingly (called as modified frame). Change of pore pressure with hydrostatic confining pressure on the cylinder boundary (the interface between the modified solid and main pore space) was simply quantified using the Gassmann equation, rather than the Biot theory [1,2]. However, the BISQ model had a major disadvantage; i.e., it simulated phase velocity well but could not well predict the quality factor [12].
Considering (intergranular) squirt from contact of grains to the main pore space, the squirt mechanism has recently been reformulated and incorporated into the Biot theory [14,15]. The work showed that the effect of squirt does not only change the fluid momentum equation but also affects the fluid mass equation.
Double porosity models [16][17][18][19][20][21][22][23] were a popular way to parameterize squirt and to simulate P wave attenuation. Pride and Berryman [19,20] constructed their model to a linear and shift-invariable system [24], and the spectrum method was extensively used. Nevertheless, the number of parameters was escalated and how to accurately quantify the detailed mechanism of squirt remains unclear in view of physics.
The quality factor is a dimensionless quantity describing the ratio of energy loss to the total energy in one cycle/period [25]. For a fixed frequency, the factor is low when attenuation is large and vice versa. Toksöz et al. [26] documented the ultrasonic velocity and quality factor of P and S waves in brine-saturated Berea sandstone. Lucet [27] measured the quality factor of P wave in water-saturated limestone, and the ultrasonic value was as low as 5-6. Sams et al. [28] measured the quality factor of P wave at the Imperial College borehole site in UK, with the result that the minimum quality factor occurred as low as 10 in the sonic band (8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24).
Local pressure difference is a problem of scale due to pore heterogeneity [29,30]. The fluid pressure in the vicinity of throat varies around the average pore pressure defined in traditional poroelasticity, which drives squirt and generates dissipation. An interesting question is whether such a mechanism is important in unconsolidated porous media. Sands are relatively homogenous and isotropic, without any microcracks or fissures. Moreover, they have very small bulk modulus and shear modulus. The skeleton is very easy to deform. If fluid pressure varies significantly on the surface of a free grain, it will move quickly to decrease the pressure difference. Therefore, via Newton's third law, stress exerted by the grains onto pore fluid should tend to be uniform and the resulting squirt should be minor.
Besides squirt, another mechanism of large P wave attenuation is dynamic permeability [31,32]. As frequency increases from zero, the velocity profile in pores will change from parabolic distribution to Stokes boundary layer [33,34] and accordingly, permeability will decrease from Darcy permeability to zero.
Matsushima et al. [35] conducted laboratory experiments in which the ultrasonic velocity and quality factor of P wave in unconsolidated Toyoura sands were measured accurately. If squirt in such sands is indeed minor as discussed above, the Biot theory [1,2] with dynamic permeability should be capable to model the measured velocity and attenuation simultaneously and vice versa. In the P wave direction, the REV (representative element volume) is two main pores connected by a throat. Squirt in this paper is defined using the pressure difference between the throat and the main pore versus the pressure difference between the two main pores. If the former pressure difference is larger than the latter pressure difference, squirt is strong, and vice versa.
Following the same idea that Biot [1,2] formulated, the kinetic energy of the system (using the coupling density ρ 12 ). Morozov and Deng [36,37] used coupling matrix, attempting to formulate the Lagrange function and dissipation function. They coupled the volume strains of the skeleton/frame and fluid with a new independent variable θ, via coupling matrix (Q and Q ' ). By setting the terms inside the coupling matrix, they controlled how the coupling of θ was with the skeleton, fluid, or both. However, the model suffered from two disadvantages: (1) variable θ had no unit, which was not quite convincing in physics; and (2) it was still a model with frequencyindependent parameters trying to summarize all specific mechanisms, but such a model does not exist. Consequently, despite P wave velocity tending to Gassmann [13] velocity, the quality factor yielding from the model was far lower than measurement. In this regard, we stand with most researchers that frequency-dependent permeability (rather than frequency-independent parameters) should be used.
In Section 2, the Biot theory [1] is improved by generalizing Darcy permeability to dynamic permeability. In Section 3, dynamic permeability at the low frequency limit is shown to be Darcy permeability. In Section 4, both the Biot theory [1] and the improved model are used to predict velocity and attenuation in unconsolidated Toyoura sands. Laboratory experiment conducted by Li [38] on compressed glass beads is used to show that the improved model remains unsuccessful due to strong squirt. The final are discussion and conclusions.

The Biot Theory Improved with Dynamic Permeability
In this section, we suppose a simple 1D scenario in which P wave propagates in the x direction and t denotes time.

Fluid Momentum Equation.
The lower part of (6.7) [1] stated the fluid momentum equation as follows: where ϕ is porosity, μ and ρ f are fluid viscosity and density, respectively, k D is Darcy permeability, ρ 12 is the coupling factor, P p is pore/fluid pressure, and v and q are solid velocity and Darcy flux rate in the wave direction, respectively. According to Newton's second law, the first and second terms on the right hand side of equation (1) are pressure gradient (the force on pore fluid by the neighboring fluid on its two ends) and fluid acceleration, respectively. The left hand side of (1) is force on pore fluid by solid. If the second terms on both sides are ignored, equation (1) will degenerate to Darcy law [4].
At a lower frequency, the term of relative acceleration, i.e., the second term on the left hand side of (1) can be neglected. Fourier Transform of equation (1) and generalizing k D to frequency-dependent permeability, k p ðωÞ where ω is angular frequency, yields in which we construct a linear and shift-invariable system [24] as Pride and Berryman [19,20] did. Equation (2) is a monochromatic equation and is stationary (both sides are independent of time).

Geofluids
Inverse Fourier Transform of equation (2) yields in which causality condition, k p ðτ < 0Þ = 0 is used. k p ðtÞ is the permeability kernel (with the unit of Darcy/s) defined as inverse Fourier Transform of k p ðωÞ. Because (2) is a monochromatic equation (k p ðωÞ is a function of argument ω), the resulting equation (3) in time domain has convolution on its right hand side. At the low frequency limit (ω → 0), k p ðωÞ → k D .

Fluid Constitutional Relation.
In the lowest part of (2.11) [1], fluid pressure is a linear combination of the volumetric strains of solid and fluid. Reversely, the volumetric strain of fluid is linearly related to the confining pressure (the mean of three normal stresses) and fluid pressure, which is the fluid constitutional equation. The equation is identical to fluid mass conservation as follows, except that the latter is the time derivative of the former [39,40].
where P C is the confining pressure and β f , β eff , and β s are compressibility coefficients of fluid, skeleton, and solid material, respectively.

Momentum Equation of
Rock System. Rock system is composed of solid and fluid. Adding the upper part with the lower part of equations (6) and (7) [1] cancels the force between solid and fluid (the internal force), yielding the momentum equation of rock system as follows: where σ 11 is the normal stress in the x direction and ρ s is skeleton density. The first and second terms on the right hand side of (5) are solid acceleration and fluid acceleration, respectively. The equation is precisely Newton's second law for rock system.

Solid Constitutional
Relation. Equation (2.12) [1] is the constitutional equation of rock system. In 1D case, it simplifies to the following equation [41,42]: where α B is the Biot-Willis coefficient, K eff is bulk modulus of the skeleton (the reciprocal of bulk compressibility β eff ), G is shear modulus of the skeleton, and u is solid displacement.
Please note that the left hand side of (6) is called effective stress [42,43]. Solving the above four equations, namely, equation (2), (4), (5), and (6) in terms of planar wave yields the wavenumber equation, an algebra equation of complex number to the fourth power. The wavenumber equation has two branches of solution. In this paper, only the branch corresponding to fast P wave is used. Finally, solving the wavenumber equation yields phase velocity (V p ) and quality factor (Q p ) as functions of frequency.

Dynamic Permeability at Low Frequency Limit
Solid motion associated with P wave is often much slower than fluid motion, such that equation (2) is approximately As long as ððωρ f k p ðωÞÞ/ðϕμÞÞ ≪ 1, the second term inside the bracket of (7) will vanish and the equation will degenerate to Darcy law. As k p ðω → 0Þ = k D , we get the condition for Darcy law to be valid: As a reminder, macroscopic equation (7) cannot yield vanishing dynamic permeability at the high frequency limit. k p ðω → +∞Þ → 0 requires investigating the microscopic profile of fluid velocity in pores [33].

Examples
4.1. Unconsolidated Toyoura Sands. The results of laboratory experiments on unconsolidated Toyoura sands [35,44] are used to investigate whether the Biot theory [1] and the model with dynamic permeability are capable of predicting the measured ultrasonic velocity and attenuation.
Parameters of Toyoura sands are listed in Table 1. The sands have grain diameter ranging from 0.1 to 0.3 mm (with an average of 0.225 mm). The shape of Toyoura sands is elongated, with high angularity, because the sands in Japan did not experience long-distance erosion. The high Darcy permeability (23 Darcy) of Toyoura sands arises from very good hydraulic connectivity between pores (inherent in the unconsolidated feature). Due to the unconsolidated feature, there was no S wave observed in the sands and so shear modulus of the skeleton is set to zero. Bulk modulus of the skeleton is also set to be far smaller than the bulk modulus of fluid. The brine salinity was 2% [35]. At temperature of 0°C, the brine was not frozen due to salt. The salinity of the brine was very low. However, water at temperature of 0°C has a viscosity significantly higher than water at a temperature of 20°C [42].

Geofluids
The coupling factor ρ 12 [1] has been estimated by Li [34]. Using the parameters in Table 1 as model input, the Biot theory [1] yields V p and Q p that are plotted in Figures 1(a)  and 1(b), respectively. The measured velocity and quality factor of P wave in the sands at temperature of 0°C was 1600 m/s and 25-33, respectively. The centroid frequency of measured P wave ranged from 350 to 600 kHz [35]. It is evident that the Biot theory [1] cannot simultaneously predict the measured velocity and quality factor of ultrasonic P wave. We can use frequency (500 kHz) and velocity (1600 m/s) to calculate the wavelength of ultrasonic P wave, which is 3.2 mm and one order of magnitude larger than the grain diameter.
As depicted in Figures 2(a) and 2(b), at a frequency of 350 kHz, the model improved with dynamic permeability of 0.074 Darcy yields V p of 1582 m/s and Q p of 25 at frequency 350 kHz, surprisingly close to the measured velocity (1600 m/s) and quality factor (25). It appears that as frequency increases from 0 to 350 kHz, dynamic permeability will decrease from 23 to 0.074 Darcy, which well confirms the theoretical trend of dynamic permeability in Section 3.

Compressed Glass Beads.
Water-saturated glass beads compressed in an apparatus [38] are used to examine whether the improved model can simulate the ultrasonically measured velocity and attenuation. Parameters of the beads are listed in Table 2. Glass beads are round in shape.
The permeability of glass beads in our experiment was not measured. Nevertheless, it can be estimated by using a Theoretical quality factor (Q p ) versus measurement (in which the minimum Q p will appear at lower frequency with increased value, as coupling factor ρ 12 decreases from zero). Note that Darcy permeability is 23 Darcy.  [45] Shear modulus of skeleton, G 0 P a Porosity, ϕ 0.41 [35] Darcy permeability, k D 23 × 10 −12 m 2 [44] Density of brine, ρ f 1020 kg m -3 [35] Viscosity of brine, μ 0.0018 Pa·s [46] Compressibility of brine, β f 4:6 × 10 −10 Pa -1 [47] 4 Geofluids comparative approach. The diameter of our glass beads was 45 μm, and the one in [48] had the bead diameter of 80 μm. The ratio of pore scale between our glass beads and that in [48] is the same as the bead diameter ratio of 0.562, because of the geometrical similarity. Using the Poiseuille flow analogy, the permeability is proportional to the square of pore diameter for a given porosity. Therefore, the permeability ratio between our experiment and [48] is the pore-scale ratio to the second power, which is 0.316. The permeability of glass beads in [48] is 1.9 Darcy, and therefore, the permeability of our glass beads is 0.60 Darcy approximately. In our experiment, the apparatus was cylindrical, placed in the vertical direction. Two transducers were placed on its top and bottom. The procedure of the experiment was this way. At first, the apparatus was filled with water and then with glass beads little by little, to ensure that no air bubbles be trapped inside the beads. Then the beads in the apparatus were compressed by a compressor, with some water outlet from four small holes on the top of the apparatus. At last, transducers emitted and received P wave in the vertical direction. During the experiment, the beads were kept compressed by the compressor. P wave in dry glass beads was measured in a similar condition, except for no any water in the apparatus.
The velocity of S wave in the compressed dry beads was measured to be 583 m/s [38], yielding shear modulus of 0.53 GPa. The velocity of P wave in the dry beads was also measured to be 1417 m/s, yielding bulk modulus of 2.42 GPa. Based on the parameters listed in Table 2, the velocity and quality factor of P wave in compressed beads are calculated and plotted in Figures 3(a) and 3(b), respectively.
With frequency of 540 kHz, the measured velocity of P wave in saturated beads was 2130 m/s, which is appreciably larger than the predicted velocity in Figure 3   Bulk compressibility of skeleton, β eff 0:413 × 10 −9 Pa -1 [38] Compressibility of solid material, β s 2:70 × 10 −11 Pa -1 [45] Shear modulus of skeleton, G 0:530 × 10 9 Pa [38] Porosity, ϕ 0.38 [38] Darcy permeability, k D 0:60 × 10 −12 m 2 Density of water, ρ f 1000 kg m -3 [46] Viscosity of water, μ 0.001 Pa·s [46] Compressibility of water, β f 4:6 × 10 −10 Pa -1 [47] 5 Geofluids method of the spectral amplitude ratio to be 21 [38]. However, the quality factor yielding from the model with dynamic permeability of 0.176 Darcy is 112 (at frequency of 540 kHz). In short, both the modeled velocity and attenuation are distinctively different from the measured counterparts, suggesting that squirt plays a significant role in P wave. We use the above frequency (540 kHz) and velocity (2130 m/s) to calculate the wavelength of ultrasonic P wave, which is 3.9 mm and two orders of magnitude larger than the grain diameter (45 μm).

Discussion
In the P wave direction, there is a throat between two main pores. Berryman (2003a, 2003b) defined squirt using throat pressure minus the main pore pressure, which is not sufficiently accurate. In the following, squirt is resolved carefully.
Without squirt, fluid pressure in the throat is the average between those in the two main pores. In this case (Figure 2), flow in the throat is one way from one main pore to the other main pore, and the pressure difference between the two main pores is very small. Consequently, the friction/viscous force associated with the slow throat flow is very small. In contrast, if the pressure difference between the throat and the main pore is much larger than the pressure difference between the two main pores, squirt in the throat is very strong and towards both pores. The large friction/viscous force associated with the fast throat squirt will significantly attenuate wave energy ( Figure 3).
As shown in Figures 2(a) and 2(b), the Biot theory [1] improved with dynamic permeability is very successful in modeling velocity and attenuation of P wave in unconsolidated sands. By nature, permeability is the square of the length to upscale viscous stress between solid and fluid [4]. k p ðωÞ has been shown to be frequency dependent by the microscopic wave solution of the Navier-Stokes equation in a pipe laterally bounded by solid [33].
In Figure 3(a), the measured velocity is appreciably larger than the model prediction, which is called stiffening of the major pores [49]. In the following, we use equation (6) to explain the phenomena.
the Biot theory is based on single porosity. Alternatively, the theory is largely dependent on the compressibility of the main pore space, rather than on contact of grains (throat). At the first order of approximation, the main pore alone (without squirt) can represent the single porosity. In this regard, equation (6) is applicable to the main pore space. At the second order of approximation, we need to consider the effect of squirt. When P wave compresses throat off the pressure equilibrium state, fluid will squirt from it to the main pore space and fluid mass in the main pore space will significantly increase, which adds a positive fluid pressure to the second term (P p ) on the left hand side of (6). Eventually, the magnitude of solid strain (∂u/∂x) is decreased, or equivalently, the modulus term inside the bracket of (6) is increased. When P wave stretches throat off the pressure equilibrium state, fluid will squirt/suck from the main pore space to it and fluid mass in the main pore space will significantly decrease, which adds a negative fluid pressure into the second term (P p ) on the left hand side of (6), eventually decreasing the magnitude of solid strain, or equivalently, increasing the modulus term inside the bracket of (6). In both phases, squirt invariably increases the modulus term inside the bracket of (6) and thus increases P wave velocity as observed in Figure 3(a).
As depicted in Figure 3(b), the theoretical attenuation is much smaller than the measured attenuation, suggesting that dynamic permeability [31,32] and the associated Stokes boundary layer [33,34] in the main pore space is not dominant in P wave attenuation. Instead, the large attenuation  6 Geofluids measured is mainly induced by squirt in the narrow space such as contact of grains or throat. This is very reasonable from the perspective of fluid mechanics, as fast flow inside narrow conduits has larger friction and consumes more mechanical energy than slow flow inside large voids does [42]. With very small bulk and shear modulus, unconsolidated sands have relative free grains. Fluid pressure onto their surface tends to be uniform (otherwise, a free grain can move to decrease pressure difference onto its surface), and therefore, the associated squirt is minor. In contrast, for compressed beads, the space at throat changes with a larger ratio than the main pore space does, resulting significant pressure difference between them which drives squirt. Consolidated rocks (represented by our compressed beads) allow fluid pressure at throat to be significantly higher than fluid pressure in the main pore space.
Using three phases composed of ice, unfrozen brine, and sands, Zhan and Matsushima [50] modeled the velocity and attenuation observed in partially frozen systems by Matsushima et al. [35]. The study suggested that squirt between ice grains is a dominant factor for P-wave attenuation around the freezing point. In contrast, our study only uses the data of the experiment with unfrozen brine in Toyoura sands, showing that squirt is minor in such cases.

Conclusions
Squirt is a well-acknowledged mechanism of P-wave attenuation in consolidated rocks. Therefore, previous studies often questioned the Biot theory for its validness in modeling attenuation.
This paper is remarkable in three aspects: (1) the Biot theory with dynamic permeability (without the consideration of squirt) is used to successfully predict both the ultrasonic velocity and quality factor measured in unconsolidated sands; (2) why squirt is minor in unconsolidated sands is well-explained, i.e., such sands have relatively free/mobile grains and thus cannot maintain large pressure difference on their surface required for strong squirt; and (3) squirt is delicately redefined (than double porosity models), which helps to better describe the detailed mechanism in the near future. The reason why the improved model remains to be incapable of predicting the ultrasonic velocity and quality factor in compressed glass beads is that the beads allow fluid pressure at throat to be significantly higher than that in the main pores. The large pressure difference arises from compliant throat, drives strong squirt (which significantly dissipates wave energy), and causes stiffening of the main pores (which appreciably increases P wave velocity).