Nonlinear Acoustic Waves in Fluid-Saturated Porous Rocks – Poro-Acoustoelasticity Theory

1. Introduction

The two elastic constants of bulk and shear moduli are insufficient to describe the nonlinear acoustic nature of solid materials under higher confining pressure. The theory of third-order elastic constants (nonlinear acoustics) has long been developed to analyze the velocity vs stress relationships in theoretical and experimental research for hyperelastic solid materials [1-7].

The 3rd order constants theory in solid material was completed with the papers by Toupin and Bernstein, Jones and Kobett and Truesdell [8-10]. Truesdell used four 3rd-order elastic constants in his general theory for isotropic solid. Brugger [11] gave the thermodynamic definition of higher order elastic constants. In 1973, Green reviewed 3rd-order constants measurements of various crystals, and gave the relations between the 3rd-order constant notations in isotropic solids by different authors [12].

From late 1980s to the present, 3rd-order nonlinear elasticity theory has been applied to rock experiments [13, 14]. For the special case of porous solid materials, Winkler et al. [15] measured third order elastic constants based on pure solid's acoustoelasticity theory in a variety of dry rocks and found that pure solid's third-order elasticity could successfully describe dry rock's velocity-stress relationships. A similar approach performed on water-saturated rocks showed that traditional 3rd-order elasticity theory for the isotropic solid could not fully describe the stress dependence of velocities in water-saturated rocks [16]. For fluid-saturated porous material, as the confining pressure increases, the unrelaxed fluid's effect must be taken into account when probing into the quantitative relationships between two-phase rock's velocities and confining pressures.

For the fluid/solid composite, Biot derived the two-phase wave equations on the basis of the linear elasticity, in which the coupling motion of solid and fluid was first analyzed. Based on the same linear assumption, some new expansions on Biot's theory have included local fluid flow, dynamic permeability and multi-scale heterogeneity in last two decades [17-23]. Around linear poroelasticity, the research interests of recent years are focused on the frequency-dependent P- and S- waves' velocity and attenuation features which are influenced by patchy saturation, pore distribution and rock microstructure.

According to Biot theory, an isotropic poroelastic medium has four independent static moduli (2nd-order elastic constants). Biot [24] developed semilinear mechanics of porous solids, in which total seven physical constants describes the semilinear properties of an isotropic media, four to characterize the linear behavior and three to characterize the nonlinear behavior. The theory is adapted to solid/fluid composite systems by Norris, Sinha and Kosteck in 1994 [25]. Biot [26] also gave eleven-constant elastic potential function for fluid-saturated porous media, in which seven 3rd-order elastic constants describe the nonlinear features of the two-phase system. Drumheller and Bedford chose a Eulerian reference frame to derive nonlinear equations for wave propagation through porous elastic media in 1980 [27], while Berryman and Thigpen chose a Lagrange reference frame in their theory [28].

Donskoy, Khashanah and McKee [29] derived the nonlinear acoustic wave equations for porous media which established a correlation between the measurable effective nonlinear parameter and structural parameters of the porous medium. Their work is based on the semilinear approximation of Biot's poroelasticity theory. The assumption of semilinearity implies a linear relation between the volume change of the solid matrix and the effective stress, therefore a modified strain (Equ. 7 in [29], Equ. 54 in [24]) is used and the effective stress is included, so that the nonlinear acoustoelasticity equations can be simplified and only three independent 3rd-order elastic constants are needed to be considered in application. Dazel and Tournat [30] considered one-dimensional problems following this theoretical approach and derived the solutions for the second harmonic Biot waves. The wave velocity dispersion and dissipation were analyzed in a semi-infinite medium. Zaitsev, Kolpakov and Nazarov [31, 32] theoretically analyzed the experimental results of the propagation of low-frequency video-pulse signals in dry and water-saturated river sand and demonstrated the nonlinearity of such a loose granular media can be non-quadratic in amplitude.

Grinfeld and Norris [33] derived wave speed formulas in closed-pore jacketed and open-pore jacketed configurations to find a complete set of seven 3rd-order elastic moduli for poroelastic medium. Nevertheless, solid and fluid's finite strain is not included in Grinfeld and Norris's work, so that 2nd-order elastic constants will not appear in the velocity square vs compression ratio slope relations for both transverse and longitudinal waves. This means there will be no compatibility between their poro-acoustoelasticity theory and the pure solid material's acoustoelasticity theory, within which 2nd-order elastic constants do affect the velocity square vs compression ratio slopes (i.e. the slope in equation (1a) is$$7\lambda +10\mu +6l+4m$$, where $$\lambda$$ and $$\mu$$ are 2nd-order elastic constants, and $$l$$ are $$m$$ the 3rd-order elastic constants.).

In this paper, we derive the wave propagation equations by substituting the 11-term potential function and the Biot kinetic energy function into Lagrange equations. Finite strain is included, and the velocities and quality factors of P- and S-waves are analytically expressed as functions of components' strains, confining pressure, pore fluid pressure and rock's 2nd- and 3rd-order moduli in cases of hydrostatic and uniaxial loading conditions. The reduction of this theory to a pure solid's acoustoelasticity theory is discussed. Numerical examples for velocity and attenuation are given which are corresponding to the open-pore jacketed and closed-pore jacketed rock tests, respectively. In the last section, we perform ultrasonic P-velocity measurements in the "open-pore jacketed" and "close-pore jacketed" tests with hydrostatic loading. The comparisons between theory and rock experimental data are then given.

2. Nonlinear wave equations for poroelasticity

2.1... Review on nonlinear acoustoelasticity for pure solid material

3rd-order elastic constants theory has long been developed to describe the wave propagation and mechanical deformation in pure solid material. The strain energy function for isotropic solid was given with three independent 3rd-order elastic constants (see Table1).

For nonlinear waves propagating in an isotropic solid media, relations between velocity square and confining pressure was derived as [3, 5, 12, 34].

$$\begin{array}{l}{\rho }_0{V}_{LP}^2=(\lambda +2\mu )-\frac{P}{3\lambda +2\mu }(7\lambda +10\mu +6l+4m),\text{(a)}\\ {\rho }_0{V}_{TP}^2=\mu -\frac{P}{3\lambda +2\mu }(3\lambda +6\mu +3m-\frac{1}{2}n),\text{(b)}\\ {\rho }_0{V}_{La}^2=(\lambda +2\mu )-\frac{P}{3\lambda +2\mu }[\frac{\lambda +\mu }{\mu }(4\lambda +10\mu +4m)+\lambda +2l],\text{(c)}\\ {\rho }_0{V}_{Lb}^2=(\lambda +2\mu )-\frac{P}{3\lambda +2\mu }[2l-\frac{2\lambda }{\mu }(\lambda +2\mu +m)],\text{(d)}\\ {\rho }_0{V}_{Ta}^2=\mu -\frac{P}{3\lambda +2\mu }\left(4\lambda +4\mu +m+\frac{\lambda n}{4\mu }\right),\text{(e)}\\ {\rho }_0{V}_{Tb}^2=\mu -\frac{P}{3\lambda +2\mu }(m+\frac{\lambda n}{4\mu }+\lambda +2\mu ),\text{(f)}\\ {\rho }_0{V}_{Tc}^2=\mu -\frac{P}{3\lambda +2\mu }(m-\frac{\lambda +\mu }{2\mu }n-2\lambda ).\text{(g)}\end{array}$$E1

Here, $$\lambda$$and $$\mu$$ are 2nd-order elastic constants.$$l$$, $$m$$and $$n$$ are 3rd-order elastic constants. $$P$$is confining pressure. $${\rho }_0$$is density. The subscripts $$L$$ and $$T$$ denote longitudinal and transverse waves. The subscript $$P$$ denotes hydrostatic loading. The subscript $$a$$ denotes waves propagating parallel to uniaxial pressure. In waves propagating perpendicular to uniaxial pressure, subscript $$b$$ and $$c$$ respectively denote waves' vibrating parallel and perpendicular to uniaxial pressure.

Brugger [11]	Truesdell [10]	Toupin & Bernstein [8]	Landau & Lifshitz [6], Goldberg [7]	Murnaghan [2, 3], Hughs & Kelly [5]
$$4C_{456}$$	$$-4\mu +4{\alpha }_6$$	$$4{\nu }_3$$	$$A$$	$$n$$
$$C_{144}$$	$$-\frac{\mu }{2}{\alpha }_4$$	$${\nu }_2$$	$$B$$	$$m-\frac{1}{2}n$$
$$\frac{1}{2}{C}_{123}$$	$$\lambda +\mu {\alpha }_3+\frac{\mu }{2}{\alpha }_4$$	$$\frac{1}{2}{\nu }_1$$	$$C$$	$$l-m+\frac{1}{2}n$$

Table 1.

Brief relations between 3rd-order elastic constants by different authors [12]

2.2. Review on linear, semilinear and nonlinear mechanics for poroelasticity

In an isotropic medium, the strain energy is a function of four variables (see equ. 3.1 in [35]), the three invariants$$I_1$$,$$I_2$$ ,$$I_3$$ of solid strain components and the fluid displacement divergence$$\zeta$$.

$$W=W(I_1,I_2,I_3,\zeta ),$$E2

Biot demonstrated that the expansion to the third degree of equation (2) leads to eleven elastic constants, and he only considered the linear relations to derive the quadratic form for $$W$$ in [35]. The semilinear mechanics for porous media was developed by Biot in [24].

In nonlinear investigations for poroelasticity, perform Taylor expansion on equation (2) to the 3rd order so that we derive the power series of four independent variables$$I_1$$,$$I_2$$ $$I_3$$and $$\zeta$$(see equ. 5.9 in [26]).

$$\begin{array}{l}2W=M_1{I}_1^2+M_2{I}_2+M_3{\zeta }^2+M_4\zeta I_1+M_5{I}_1^3+M_6{I}_3+M_7{\zeta }^3\\ +M_8{I}_1{I}_2+M_9\zeta I_2+M_{10}{I}_1^2\zeta +M_{11}{I}_1{\zeta }^2.\end{array}$$E3

Equation (3) is an elastic potential expression with four 2nd-order elastic constants and seven 3rd-order elastic constants.

We list in table 2 the different notations for 2nd- and 3rd-order elastic constants which have been used in literatures of poroelasticity.

There are 4 2nd-order terms in equation (3). $$M_1{I}_1^2$$(the dilatation term) and $$M_2{I}_2$$ (the shear term) are mainly dependent on the elastic characteristics of the solid skeleton, $$M_3{\zeta }^2$$is related to the pore fluid, and $$M_4\zeta I_1$$ is associated with the coupling effect between solid and fluid. If pore fluid has a extremely high bulk modulus (fluid's 2nd-order elastic constants) while solid grain has relatively lower bulk modulus for a specific porous material, $$M_3{\zeta }^2$$may become the most significant term. But as to the natural rocks, solid grain's bulk modulus is much higher than pore fluid (generally solid modulus is one magnitude higher than fluid), therefore, $$M_1{I}_1^2$$and $$M_2{I}_2$$ are the most significant, $$M_4\zeta I_1$$is the moderate, and $$M_3{\zeta }^2$$ is the least significant. For a gas-saturated rock (in this case pore fluid has a very low elastic constant), $$M_4\zeta I_1$$and $$M_3{\zeta }^2$$ tend to zero, and $$M_1{I}_1^2$$ and $$M_2{I}_2$$ dominate the effective elastic nature of the fluid/solid composite. For the rocks saturated with water or light oil, $$M_4\zeta I_1$$is generally more significant than$$M_3{\zeta }^2$$, nevertheless, $$M_4\zeta I_1$$'s effects will get weaker with a high-porosity and softer solid skeleton. For the high-porosity poor-consolidated sand/sandstone, the term of $$M_3{\zeta }^2$$ may become more significant than$$M_4\zeta I_1$$.

	Biot [17, 18]; Ba; [23, 36]	Biot [26]	Biot [24]; Donskoy, Khashanah & McKee [29]; Dazel & Tournat [30]	Pierce [37]; Norris, Sinha & Kostek; [25, 38]	Berryman & Thigpen; [28]	Grinfeld & Norris [33]	Bemer, Boutéca, Vincké, Hoteit & Ozanam [39]
2nd	$$A$$,$$N$$ ,$$Q$$ ,$$R$$	$$C_1$$,$$C_2$$ ,$$C_6$$ ,$$C_8$$	$$\lambda$$,$$\mu$$ ,$$\alpha$$ ,$$M$$	$$A$$,$$N_1$$	$$a$$,$$b$$ ,$$c$$ ,$$d$$	$${\Gamma }_{11}$$,$${\Gamma }_2$$ ,$${\Gamma }_{1\zeta }$$ ,$${\Gamma }_{\zeta \zeta }$$	$$K_0$$,$$K_s$$ ,$$K_{fl}$$ ,$${\varphi }_0$$
3rd	-	$$C_3$$,$$C_4$$ ,$$C_5$$ ,$$C_7$$ ,; $$C_9$$,$$C_{10}$$ ,$$C_{11}$$	$$D$$,$$F$$ ,$$G$$	$$B$$	$$e$$,$$f$$ ,$$g$$ ,$$m$$ ,$$n$$	$${\Gamma }_{111}$$,$${\Gamma }_{12}$$ ,$${\Gamma }_3$$ ,$${\Gamma }_{1\zeta \zeta }$$ ,; $${\Gamma }_{11\zeta }$$,$${\Gamma }_{2\zeta }$$ ,$${\Gamma }_{\zeta \zeta \zeta }$$	$$D$$,$$F$$

Table 2.

Different notations of 2nd and 3rd order elastic constants for poroelasticity

	1st	2nd	3rd	4 th	5th
Soft solid skeleton saturated with extremely high-modulus fluid	$$M_7{\zeta }^3$$	$$M_{11}{I}_1{\zeta }^2$$	$$M_{10}{I}_1^2\zeta$$	$$M_9\zeta I_2$$	$$M_5{I}_1^3$$,$$M_6{I}_3$$ ,$$M_8{I}_1{I}_2$$
Rock saturated with water/oil	$$M_5{I}_1^3$$,$$M_6{I}_3$$ ,$$M_8{I}_1{I}_2$$	$$M_{10}{I}_1^2\zeta$$	$$M_9\zeta I_2$$	$$M_{11}{I}_1{\zeta }^2$$	$$M_7{\zeta }^3$$
Rock saturated with gas	$$M_5{I}_1^3$$,$$M_6{I}_3$$ ,$$M_8{I}_1{I}_2$$	-	-	-	$$M_7{\zeta }^3$$,$$M_9\zeta I_2$$ ,$$M_{10}{I}_1^2\zeta$$ ,$$M_{11}{I}_1{\zeta }^2$$

Table 3.

The significance of relative-magnitude of the 7 3rd-order terms in the nonlinear elastic potential expression

As to the 7 3rd-order terms in equation (3), $$M_5{I}_1^3$$, $$M_6{I}_3$$and $$M_8{I}_1{I}_2$$ are related to solid skeleton, $$M_7{\zeta }^3$$is related to pore fluid, while$$M_9\zeta I_2$$, $$M_{10}{I}_1^2\zeta$$and $$M_{11}{I}_1{\zeta }^2$$ are related to the fluid-solid coupling effect. The significance of the relative magnitude of these terms are listed in Table 3. As to the actual rocks in nature, the terms of$$M_5{I}_1^3$$, $$M_6{I}_3$$and $$M_8{I}_1{I}_2$$ are in the most significant level. In gas-saturated rocks, $$M_7{\zeta }^3$$, $$M_9\zeta I_2$$, $$M_{10}{I}_1^2\zeta$$and $$M_{11}{I}_1{\zeta }^2$$ can be neglected in comparison with$$M_5{I}_1^3$$, $$M_6{I}_3$$and$$M_8{I}_1{I}_2$$, which means fluid's effect can be completely neglected. Therefore, the acoustoelasticity theory for a pure solid can be sufficient to describe the nonlinear wave phenomena in gas-saturated or dry rock samples. But as the fluid's elastic modulus increases, the terms of $$M_9\zeta I_2$$ and $$M_{10}{I}_1^2\zeta$$ will play an important role in its nonlinear elastic effects on the whole solid/fluid system. In this case, the pure solid's three 3rd-order elastic constants may be not enough to give a full description for nonlinear phenomena in a water/oil-saturated rock, which has been discussed in the experimental studies by Winkler et al [15, 16].

2.3. Finite strain in solid/fluid system

If the fluid-saturated rock is loaded under high confining pressure, the infinitesimal strain expression is insufficient to describe solid’s and fluid’s microscale finite deformation. On this occasion, the Lagrangian strain tensor has much higher precision than infinitesimal strain for solid material. It can be written

$${\epsilon }_{ij}=\frac{1}{2}(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}+\frac{\partial u_l}{\partial x_i}\frac{\partial u_l}{\partial x_j}),i,j,l=1,2,3$$E4

where $$u_i$$ denotes solid displacement in $$i$$ direction. The convention of summation over repeated indices is adopted.

Different from Lagrangian description in solid material, the nonlinear problems of the acoustics of fluids are usually formulated in terms of an Eulerian description of wave motion [38, 40, 41]. Kostek, Sinha and Norris [38] gave the explicit relations for the two descriptions in inviscid fluid. In dealing with nonlinear problems involving both fluid and solid, the unified treatment of Lagrangian variables is used in this paper. If an assumption is reasonably made that fluid viscosity's effect on rock's shear deformation can be neglected so that fluid's shear deformation is completely ignored, the fluid's finite strain under high confining pressure can be approximately written

$${\zeta }_{ii}=\frac{\partial U_i}{\partial x_i}+\frac{1}{2}\frac{\partial U_i}{\partial x_i}\frac{\partial U_i}{\partial x_i},i=1,2,3$$E5

where $$U_i$$ denotes fluid displacement in $$i$$ direction. The convention of summation over repeated indices is not adopted in equation (5). Here we only take into account the longitudinal finite deformation for those low viscous fluid like water and gas. For some non-Newtonian fluid such as bitumen and heavy oil, finite shear deformation need be considered.

In this paper, two types of rock experimental configurations, "open-pore jacketed" and "closed-pore jacketed", are considered. The "closed-pore jacketed" system [33] corresponds to constancy of pore fluid mass, which means the fluid/solid composite are closed-pore jacketed with a impervious deformable jacket under confining pressure so that pore fluid cannot flow out solid matrix in experiments. The "closed-pore jacketed" system is the simplest and rather difficult to realize in practice. The "open-pore jacketed" system [33, 42] corresponds to the constancy of fluid pressure, which means the fluid/solid composite are closed-pore jacketed with a thin impermeable jacket and the inside of the jacket is made to communicate with the atmosphere through a tube, so that in experiment pore fluid will flow out under confining pressure. The "open-pore jacketed" configuration is very common in actual rock tests.

Pore fluid will not be finitely deformed like solid if the porous solid skeleton is "open-pore jacketed", especially in case of highly permeable rocks being saturated with low viscous fluid. In the open-pore jacketed porous media, fluid particles would rather flow out the matrix through connected pores and throats when being subject to pressure gradient, rather than get compressed into a finite deformation state. On the other hand, if the porous structure is "closed-pore jacketed" or the whole fluid system could be regarded as a closed system, finite deformation of the inclosed pore fluid will happen.

Several wave speed formulas for specific experimental configurations are derived by Grinfeld and Norris [33]. However, finite strain (equations 4~5) has not been considered in their derivation, therefore their poro-acoustoelasticity velocity expressions can not be reduced to pure solid's acoustoelasticity expressions (equations 1a~1g) if we perform a gedanken experiment by replacing pore fluid with solid in a solid/fluid composite.

2.4. Nonlinear wave equations in fluid-saturated solid structure

In this paper, we use a theoretical approach similar to Biot's [35] and Norris's [33], but in our development finite strain is used instead of infinitesimal strain, so that this newly developed theory will be more appropriate to describe wave phenomena in the finite-deformed solid/fluid composed system.

The dissipation function and the kinetic energy of a unit volume for the isotropic fluid-solid composed system [35] is given by

$$2D={\varphi }^2\frac{\eta }{k}({({\dot{u}}_1-{\dot{U}}_1)}^2+{({\dot{u}}_2-{\dot{U}}_2)}^2+{({\dot{u}}_3-{\dot{U}}_3)}^2)$$E6

$$2T={\rho }_{11}({\dot{u}}_1^2+{\dot{u}}_2^2+{\dot{u}}_3^2)+2{\rho }_{12}({\dot{u}}_1{\dot{U}}_1+{\dot{u}}_2{\dot{U}}_2+{\dot{u}}_3{\dot{U}}_3)+{\rho }_{22}({\dot{U}}_1^2+{\dot{U}}_2^2+{\dot{U}}_3^2)$$E7

Mass parameters$${\rho }_{11}$$, $${\rho }_{12}$$and $${\rho }_{22}$$ have been defined by [17,18].$$\varphi$$, $$\eta$$and $$k$$ are the constants of porosity, fluid viscosity and permeability.

By applying Lagrange's equations and take $$u_i$$ and $$U_i$$ as generalized coordinates, generalized forces of solid and fluid phase are derived.

$$\begin{array}{l}{f}_i=\frac{d}{dt}(\frac{\partial T}{\partial {\dot{u}}_i})+\frac{\partial D}{\partial {\dot{u}}_i},\text{(a)}\\ F_i=\frac{d}{dt}(\frac{\partial T}{\partial {\dot{U}}_i})+\frac{\partial D}{\partial {\dot{U}}_i}\text{(b)}\end{array}$$E8

where $$f_i$$ and $$F_i$$ satisfy

$$f_i=\frac{d}{d{x}_j}(\frac{\partial W}{\partial (\frac{\partial u_i}{\partial x_j})})$$E9

,

$$F_i=\frac{d}{d{x}_j}(\frac{\partial W}{\partial (\frac{\partial U_i}{\partial x_j})})$$E10

.

Moreover,

$$a_{ij}=\frac{\partial W}{\partial {\epsilon }_{ij}},b_{ij}=\frac{\partial W}{\partial {\zeta }_{ij}},i,j=1,2,3$$E11

and

$$J={\left[\begin{array}{l}1+u_{1,1}{u}_{1,2}{u}_{1,3}\\ u_{2,1}1+u_{2,2}{u}_{2,3}\\ u_{3,1}{u}_{3,2}1+u_{3,3}\end{array}\right]}^T,K={\left[\begin{array}{l}1+U_{1,1}00\\ 01+U_{2,2}0\\ 001+U_{3,3}\end{array}\right]}^T$$E12

where $$u_{i,j}$$ designates$$\frac{\partial u_i}{\partial x_j}$$.

By substituting equations (9~10) into equations (8a~b), nonlinear wave equations in 2-phase medium are written as

$$\begin{array}{l}{a}_{ik,j}{J}_{kj}+a_{ik}{J}_{kj,j}={\rho }_{11}{\ddot{u}}_i+{\rho }_{12}{\ddot{U}}_i+b({\dot{u}}_i-{\dot{U}}_i)\text{(a)}\\ b_{ik,j}{K}_{kj}+b_{ik}{K}_{kj,j}={\rho }_{12}{\ddot{u}}_i+{\rho }_{22}{\ddot{U}}_i-b({\dot{u}}_i-{\dot{U}}_i),\text{(b)}\end{array}$$E13

where

$$b={\varphi }^2\frac{\eta }{k}$$E14

.

By substituting equations (4~5) into equation (3), substituting equation (3) into equation (9), then solving $$a_{ij}$$ and$$b_{ij}$$, and finally substituting into (11a~b), we derive the nonlinear acoustic wave propagation equations in 2-phase medium in which both solid and fluid's finite strain are considered. The full expansion scheme for equations (11a~b) is very complicated. In this paper, we only give simplified schemes for some particular experimental configurations.

3. Wave speeds for specific experimental configurations

3.1. The 3rd-order elastic constants for solid/fluid composite

There are four 2nd-order and seven 3rd-order elastic constants in poro-acoustoelasticity theory. The determination of four 2nd-order elastic constants has been discussed by [42]. In this subsection methods of measurement are described for the determination of the seven 3rd-order elastic constants. With the determination of all the 11 elastic constants, the poro-acoustoelasticity theory will be directly applicable to nonlinear solid/fluid systems.

Studies of Winkler et al. [16] showed that traditional 3rd-order elastic constants theory of pure solid material gives much better description for observed velocity-pressure relations in dry rocks than in saturated ones. Therefore we assume the three 3rd-order elastic constants of$$l_b$$, $$m_b$$and $$n_b$$ are sufficient for describing nonlinear acoustic features in dry rocks.

Three gedanken experiments have been discussed to determine the four 2nd-order elastic constants for poroelasticity. Similar approach was also adopted on the determination of the six 2nd-order elastic constants in double-porosity models [43]. For the case of poro- acoustoelasticity, the problem comes to much more complicated for determining the total seven 3rd-order elastic constants.

Five gedanken experiments are designed.

1. Because the fluid viscosity effect has been neglected and it will not directly contribute to the whole structure's higher-order distortion, the 3rd-order elastic coefficient of $$I_3$$ in equation (3) will satisfy$$M_6=2n_b$$.

2. If the fluid-saturated rock sample surrounded by a flexible rubber is subject to a hydrostatic pressure $$P_h$$ and the fluid is allowed to squirt out, all the outside pressure will be confined to the solid skeleton because fluid will squirt out when rock sample is squeezed. For the generalized stress tensors for solid and fluid phase $$a_{ij}$$ and$$b_{ij}$$, because the rock sample only responds to frame stiffness and fluid pressure is nearly zero,

$$\begin{array}{l}{a}_{ii}=P_h=K_{b}e+l_b{e}^2+\frac{1}{9}{n}_b{e}^2,\text{(a)}\\ b_{ii}=0,\text{(b)}\end{array}$$E15

where$$e={\epsilon }_{11}+{\epsilon }_{22}+{\epsilon }_{33}$$. $$K_b$$is the bulk modulus of solid matrix.

For isotropic rocks under hydrostatic loading,$${\epsilon }_{11}={\epsilon }_{22}={\epsilon }_{33}$$. By substituting equations (12a~b) into equations (9), we get one relation between$$M_5$$, $$M_7$$, $$M_8$$, $$M_9$$, $$M_{10}$$, $$M_{11}$$and$$l_b$$, $$m_b$$,$$n_b$$.

1. The fluid-saturated rock sample surrounded by a flexible rubber is subject to a uniaxial pressure $$P_u$$ along axis 1, fluid being allowed to squirt out.

The solid and fluid pressures are expressed as

$$\begin{array}{l}{a}_{11}=P_u=(K_b+\frac{4}{3}{\mu }_b)e-2{\mu }_b({\epsilon }_{22}+{\epsilon }_{33})+(l_b+2m_b)e^2\\ -2m_b[e({\epsilon }_{22}+{\epsilon }_{33})+{\epsilon }_{22}{\epsilon }_{33}+{\epsilon }_{11}{\epsilon }_{33}+{\epsilon }_{11}{\epsilon }_{22}]+n_b{\epsilon }_{22}{\epsilon }_{33},\text{(a)}\\ a_{22}=a_{33}=0=(K_b+\frac{4}{3}{\mu }_b)e-2{\mu }_b({\epsilon }_{11}+{\epsilon }_{22})+(l_b+2m_b)e^2\\ -2m_b[e({\epsilon }_{11}+{\epsilon }_{33})+{\epsilon }_{22}{\epsilon }_{33}+{\epsilon }_{11}{\epsilon }_{33}+{\epsilon }_{11}{\epsilon }_{22}]+n_b{\epsilon }_{11}{\epsilon }_{33},(\text{b)}\\ b_{ii}=0(\text{c})\end{array}$$E16

where solid strain in two transverse directions are identical as$${\epsilon }_{22}={\epsilon }_{33}$$.

By substituting equations (13a~c) into equations (9), we get one more relation.

1. Let the open-pore jacketed porous sample be subjected to a uniform hydrostatic pressure$$P_f$$.

In the acoustoelasticity of solid grains and fluid, strains can be solved through

$$\begin{array}{l}{P}_f=K_{s}e+K_{s3}{e}^2,\text{(a)}\\ P_f=K_f\zeta +K_{f3}{\zeta }^2\text{(b)}\\ a_{ii}=(1-\varphi )P_f,\text{(c)}\\ b_{ii}=\varphi P_f\text{(d)}\end{array}$$E17

where $$K_s$$ and $$K_f$$ are the solid grain and fluid's bulk modulus. $$K_{s3}$$and $$K_{f3}$$ indicate the solid grain and fluid's 3rd-order moduli [38].

$$\varphi$$is the porosity.

By substituting equations (14a~d) into equations (9), we get another two relations.

1. The closed-pore jacketed porous rock is subjected to a uniaxial pressure $$P_u$$ along axis 1, while at the same time the pore pressure is kept with$$P_p$$. Let$$P_u=P_p$$, we derive the last two relations as below.

$$\begin{array}{l}{P}_p=K_s({{\epsilon }^{\prime }}_{11}+{{\epsilon }^{\prime }}_{22}+{{\epsilon }^{\prime }}_{33})+K_{s3}{(}^{{{\epsilon }^{\prime }}_{11}},\text{(a)}\\ P_p=K_f\zeta +K_{f3}{\zeta }^2(\text{b)}\\ 0=(K_b+\frac{4}{3}{\mu }_b)e^{″}-2{\mu }_b({{\epsilon }^{″}}_{22}+{{\epsilon }^{″}}_{33})+(l_b+2m_b){e^{″}}^2\\ -2m_b[e^{″}({{\epsilon }^{″}}_{22}+{{\epsilon }^{″}}_{33})+{{\epsilon }^{″}}_{22}{{\epsilon }^{″}}_{33}+{{\epsilon }^{″}}_{11}{{\epsilon }^{″}}_{33}+{{\epsilon }^{″}}_{11}{{\epsilon }^{″}}_{22}]+n_b{{\epsilon }^{″}}_{22}{{\epsilon }^{″}}_{33},\text{(c)}\\ -P_p=(K_b+\frac{4}{3}{\mu }_b)e^{″}-2{\mu }_b({{\epsilon }^{″}}_{11}+{{\epsilon }^{″}}_{22})+(l_b+2m_b){e^{″}}^2\\ -2m_b[e^{″}({{\epsilon }^{″}}_{11}+{{\epsilon }^{″}}_{33})+{{\epsilon }^{″}}_{22}{{\epsilon }^{″}}_{33}+{{\epsilon }^{″}}_{11}{{\epsilon }^{″}}_{33}+{{\epsilon }^{″}}_{11}{{\epsilon }^{″}}_{22}]+n_b{{\epsilon }^{″}}_{11}{{\epsilon }^{″}}_{33}.\text{(d)}\\ a_{11}=(1-\varphi )P_p(\text{e)}\\ a_{22}=a_{33}=-\varphi P_p\text{(f)}\\ b_{ii}=\varphi P_p(\text{g})\end{array}$$E18

The total strain in solid frame is then calculated with$${\epsilon }_{11}={{\epsilon }^{\prime }}_{11}+{{\epsilon }^{″}}_{11}$$, $${\epsilon }_{22}={{\epsilon }^{\prime }}_{22}+{{\epsilon }^{″}}_{22}$$, $${\epsilon }_{33}={{\epsilon }^{\prime }}_{33}+{{\epsilon }^{″}}_{33}$$, where $${{\epsilon }^{\prime }}_{ij}$$ and $${{\epsilon }^{″}}_{ij}$$ denote solid strains induced by pore pressure and effective pressure respectively.

By substituting equations (15a~g) into equations (9), we get the last two relations.

The total seven relations in gedanken experiments (1~5) are well-posed to solve out the seven unknown 3rd-order elastic constants.

3.2. Wave velocity expressions under hydrostatic confining pressure

Based on poro-acoustoelasitcity theory, relations between acoustic wave velocity and loading pressure will be derived, which can be directly applicable to the actual measurements for fluid-saturated porous solids, especially in rock experiments.

In this paper we only consider the case of rocks being subject to hydrostatic or uniaxial loading in experiments. No transverse stress is loaded on samples therefore both shear deformation of solid and fluid can be neglected.

For simplicity let us consider the plane P-waves propagating along the axis ''1". Wave equations (11a~b) are reduced to

$$\begin{array}{l}{a}_{11,1}+a_{11,1}{u}_{1,1}+a_{11}{u}_{1,11}={\rho }_{11}{\ddot{u}}_1+{\rho }_{12}{\ddot{U}}_1+b({\dot{u}}_1-{\dot{U}}_1),\text{(a)}\\ b_{11,1}+b_{11,1}{U}_{1,1}+b_{11}{U}_{1,11}={\rho }_{12}{\ddot{u}}_1+{\rho }_{22}{\ddot{U}}_1-b({\dot{u}}_1-{\dot{U}}_1)\text{(b)}\end{array}$$E19

For hydrostatic loading, the "large" part of the static strain in "small-on-large" constitutive relations induced by static confining pressure is expressed as

$${\epsilon }_{ij}=\left(\begin{array}{l}\alpha 00\\ 0\alpha 0\\ 00\alpha \end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}\beta 00\\ 0\beta 0\\ 00\beta \end{array}\right)$$E20

where $$\alpha$$ is a component of the solid's deformation in the 3D space under hydrostatic loading, while $$\beta$$ is a component of the fluid's deformation.

The small dynamic strain induced by the plane P-waves' vibration along axis 1 is

$${\epsilon }_{ij}=\left(\begin{array}{l}S{e}^{i[\omega t-k{x}_1(1+\alpha )]}00\\ 000\\ 000\end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}{S}^{\prime }{e}^{i[\omega t-k{x}_1(1+\beta )]}00\\ 000\\ 000\end{array}\right)$$E21

where $$S$$ and $$S^{\prime }$$ are wave amplitudes in solid and fluid phases. $$\omega$$and $$k$$ denote angle frequency and wave number.

By substituting equations (17~18) into equations (16a~b), neglecting all terms whose power orders are higher than 1 and eliminating $$S$$ and$$S^{\prime }$$, equations (16a~b) will be simplified to

$$\left|\begin{array}{l}{\Upsilon }_1{k}^2-{\rho }_{11}{\omega }^2+ib\omega ({\Upsilon }_2+M_4\beta )k^2-{\rho }_{12}{\omega }^2-ib\omega \\ ({\Upsilon }_2+M_4\alpha )k^2-{\rho }_{12}{\omega }^2-ib\omega {\Upsilon }_3{k}^2-{\rho }_{22}{\omega }^2+ib\omega \end{array}\right|=0,$$E22

where

$${\Upsilon }_1=M_1+(7M_1+M_2+9M_5+2M_8)\alpha +(\frac{3}{2}{M}_4+3M_{10})\beta ,$$E23

$${\Upsilon }_2=\frac{1}{2}{M}_4+(\frac{1}{2}{M}_4+M_9+3M_{10})\alpha +(3M_{11}+\frac{1}{2}{M}_4)\beta ,$$E24

and

$${\Upsilon }_3=M_3+(\frac{3}{2}{M}_4+3M_{11})\alpha +(9M_7+7M_3)\beta .$$E25

Frequency-dependent P velocity can be solved through$$V=k/\omega$$. The fast and slow P-wave velocities respectively correspond to the two solutions of the quadratic equation.

Similar to the reduction process from Biot equations to Gassmann equations [44] in zero frequency limit, if we neglect the relative motion between solid matrix and pore fluid, which means fluid particles have the same vibration amplitudes as solid in wave propagation (that is,$$S=S^{\prime }$$), equation (19) will be reduced to

$$\rho V^2={\Upsilon }_1+2{\Upsilon }_2+{\Upsilon }_3+M_4(\alpha +\beta ).$$E26

where $$\rho ={\rho }_{11}+2{\rho }_{12}+{\rho }_{22}$$ denotes the average density [35].

Equation (20) is also a nonlinear extension of Gassmann theory. The difference between solid and fluid's "small" dynamic strain induced by wave vibration is neglected, while solid and fluid's "large" static strain induced by loading are still different with$$\alpha \ne \beta$$.

Another necessary condition must be satisfied to reduce equation (20) to equation (1a). By assuming pore fluid have the same static deformation as solid skeleton under hydrostatic loading, which means solid and fluid share the same constitutive relation, both in "small" dynamic and "large" static part, and then substituting $$\alpha =\beta$$ in equation (20),

$$\begin{array}{l}\rho V^2=M_1+M_3+M_4\\ +[7(M_1+M_3+M_4)+M_2+9\left(M_5+M_7+M_{10}+M_{11})+2(M_8+M_9\right))]\alpha \end{array}$$E27

Equation (21) is identical to equation (1a) if we take the average body deformation of $$-\frac{P}{3\lambda +2\mu }$$ for$$\alpha$$. The constant relations between different studies are listed in table 4 (As to the elastic potential $$W$$ in equation 3, a similar notation in equation 5.9 in [26] has been defined as a function of the relative variation of fluid content$$m$$,$$m=\varphi (\zeta -I_1)$$. Since [26] has used a different invariant for fluid phase, the constant relations between [26] and this paper are very complicated. If we assume $$\zeta$$ is also used in equation 5.9 in [26] by replacing $$m$$ with$$\zeta$$, the constants can be listed in the column 4 of table 4.)

This paper	Biot [17, 18]	Murnaghan [2]; Hughs & Kelly [5]	Biot [26]	Grinfeld &; Norris [33]
$$M_1$$	$$A+2N$$, $$P$$	-	$$2C_1$$	$$2{\Gamma }_{11}$$
$$M_3$$	$$R$$	-	$$2C_6$$	$$2{\Gamma }_{\zeta \zeta }$$
$$M_4$$	$$2Q$$	-	$$2C_8$$	$$2{\Gamma }_{1\zeta }$$
$$M_1+M_3+M_4$$	$$A+2N+2Q+R$$	$$\lambda +2\mu$$	-	-
$$\frac{M_2}{4}$$	$$-N$$	$$-\mu$$	$$\frac{C_2}{2}$$	$$\frac{{\Gamma }_2}{2}$$
$$M_5$$	-	-	$$2C_3$$	$$2{\Gamma }_{111}$$
$$M_6$$	-	$$2n$$	$$2C_5$$	$$2{\Gamma }_3$$
$$M_7$$	-	-	$$2C_7$$	$$2{\Gamma }_{\zeta \zeta \zeta }$$
$$M_8$$	-	-	$$2C_4$$	$$2{\Gamma }_{12}$$
$$M_9$$	-	-	$$2C_{11}$$	$$2{\Gamma }_{2\zeta }$$
$$M_{10}$$	-	-	$$2C_9$$	$$2{\Gamma }_{11\zeta }$$
$$M_{11}$$	-	-	$$2C_{10}$$	$$2{\Gamma }_{1\zeta \zeta }$$
$$\frac{M_5+M_7+M_{10}+M_{11}}{2}$$	-	$$\frac{l+2m}{3}$$	-	-
$$\frac{M_8+M_9}{2}$$	-	$$-2m$$	-	-

Table 4.

able 4.The relations between elastic constants as used in different studies. (The relation between column 1 and column 3 holds when poro-acoustoelasticity theory is reduced to acoustoelasticity theory.)

For plane S-waves propagating along 2 axis and vibrating along 1 axis,

$$\begin{array}{l}{a}_{12,2}+a_{12,2}{u}_{1,1}+a_{22}{u}_{1,22}={\rho }_{11}{\ddot{u}}_1+{\rho }_{12}{\ddot{U}}_1+b({\dot{u}}_1-{\dot{U}}_1)\text{(a)}\\ 0={\rho }_{12}{\ddot{u}}_1+{\rho }_{22}{\ddot{U}}_1-b({\dot{u}}_1-{\dot{U}}_1)(\text{b)}\end{array}$$E28

The small dynamic strain induced by S-waves is

$${\epsilon }_{ij}=\left(\begin{array}{l}S{e}^{i[\omega t-k{x}_2(1+\alpha )]}00\\ 000\\ 000\end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}{S}^{\prime }{e}^{i[\omega t-k{x}_2(1+\beta )]}00\\ 000\\ 000\end{array}\right)$$E29

By substituting equation (17) and equation (23) into equations (22a~b),

$$\left|\begin{array}{l}\Psi k^2-{\rho }_{11}{\omega }^2+ib\omega -{\rho }_{12}{\omega }^2-ib\omega \\ -{\rho }_{12}{\omega }^2-ib\omega -{\rho }_{22}{\omega }^2+ib\omega \end{array}\right|=0$$E30

where

$$\Psi =-\frac{M_2}{4}+(3M_1-\frac{M_6}{4}-\frac{3}{4}{M}_8)\alpha +(\frac{3}{2}{M}_4-\frac{3}{4}{M}_9)\beta$$E31

.

Frequency-dependent S-velocity can be solved through equation (24).

By assuming pore fluid has the same static deformation as solid skeleton, equation (24) comes to

$$\rho V^2=-\frac{M_2}{4}+(3M_1+\frac{3}{2}{M}_4-\frac{M_6}{4}-\frac{3}{4}{M}_8-\frac{3}{4}{M}_9)\alpha .$$E32

Equation (25) can not be directly reduced to equation (1b), because low viscous fluid's shear constitutive functions has been neglected in Biot theory. Return to pure solid's nonlinear acoustic theory needs to regard pore fluid as solid grains with the same shear constitutive behaviors. Therefore, neglecting solid/fluid relative deformation and adding the neglected term of fluid's shear deformation, S-velocity expression is reduced to

$$\rho V^2=-\frac{M_2}{4}+(3M_1+3M_4+3M_3-\frac{M_6}{4}-\frac{3}{4}{M}_8-\frac{3}{4}{M}_9)\alpha ,$$E33

which is identical to equation (1b).

3.3. Velocity expressions under Uniaxial confining pressure

For uniaxial loading along axis 1, the "large" static strain in "small-on-large" constitutive relations induced by confining pressure is expressed as

$${\epsilon }_{ij}=\left(\begin{array}{l}\alpha 00\\ 0{\alpha }^{\prime }0\\ 00{\alpha }^{\prime }\end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}\beta 00\\ 0\beta 0\\ 00\beta \end{array}\right)$$E34

The small dynamic strain induced by plane P-waves vibration along axis 1 is

$${\epsilon }_{ij}=\left(\begin{array}{l}S{e}^{i[\omega t-k{x}_1(1+\alpha )]}00\\ 000\\ 000\end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}{S}^{\prime }{e}^{i[\omega t-k{x}_1(1+\beta )]}00\\ 000\\ 000\end{array}\right)$$E35

By substituting equations (27~28) into equations (16a~b),

$$\left|\begin{array}{l}{\Upsilon }_1{k}^2-{\rho }_{11}{\omega }^2+ib\omega ({\Upsilon }_2+M_4\beta )k^2-{\rho }_{12}{\omega }^2-ib\omega \\ ({\Upsilon }_2+M_4\alpha )k^2-{\rho }_{12}{\omega }^2-ib\omega {\Upsilon }_3{k}^2-{\rho }_{22}{\omega }^2+ib\omega \end{array}\right|=0,$$E36

where

$${\Upsilon }_1=M_1+(5M_1+3M_5)\alpha +(2M_1+M_2+6M_5+2M_8){\alpha }^{\prime }+(\frac{3}{2}{M}_4+3M_{10})\beta$$E37

,

$${\Upsilon }_2=\frac{1}{2}{M}_4+(\frac{1}{2}{M}_4+M_{10})\alpha +(M_9+2M_{10}){\alpha }^{\prime }+(3M_{11}+\frac{1}{2}{M}_4)\beta$$E38

and

$${\Upsilon }_3=M_3+(\frac{1}{2}{M}_4+M_{11})\alpha +(M_4+2M_{11}){\alpha }^{\prime }+(9M_7+7M_3)\beta$$E39

.

The two frequency-dependent P velocities are solved out from equation (29). The equation will be reduced to equation (19) in hydrostatic loading configuration if $$\alpha ={\alpha }^{\prime }$$ is satisfied.

If the difference between solid and fluid deformation is neglected for uniaxial loading,

$$\rho V^2={{\Upsilon }^{\prime }}_1+2{{\Upsilon }^{\prime }}_2+{{\Upsilon }^{\prime }}_3+M_4(\alpha +\beta ).$$E40

Equation (30) can not be directly reduced to equation (1c) like the reduction process from (19) to (1a) in hydrostatic case. In equation (30), if we replace the fluid strain $$\beta$$ with the solid strain $$\alpha$$ along axis 1 and replace $$\beta$$ with $${\alpha }^{\prime }$$ along axis 2 and axis 3 (the two orthogonal directions perpendicular to axis 1), we can derive an analytical scheme which is completely compatible with equation (1c).

$$\begin{array}{l}\rho V^2=M_1+M_3+M_4\\ +[5(M_1+M_3+M_4)+3(M_5+M_7+M_{10}+M_{11})]\alpha \\ +[2(M_1+M_3+M_4)+M_2+6(M_5+M_7+M_{10}+M_{11})+2(M_8+M_9)]{\alpha }^{\prime }\end{array}$$E41

where

$$\alpha \approx -\frac{(\lambda +\mu )P_u}{\mu (3\lambda +2\mu )},{\alpha }^{\prime }\approx \frac{\lambda P_u}{2\mu (3\lambda +2\mu )}$$E42

.

With the same "large" static strain, for transverse waves propagating perpendicular to uniaxial loading and vibrating along loading,

$${\epsilon }_{ij}=\left(\begin{array}{l}S{e}^{i[\omega t-k{x}_2(1+{\alpha }^{\prime })]}00\\ 000\\ 000\end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}{S}^{\prime }{e}^{i[\omega t-k{x}_2(1+\beta )]}00\\ 000\\ 000\end{array}\right),$$E43

Then

$$\left|\begin{array}{l}\Psi k^2-{\rho }_{11}{\omega }^2+ib\omega -{\rho }_{12}{\omega }^2-ib\omega \\ -{\rho }_{12}{\omega }^2-ib\omega -{\rho }_{22}{\omega }^2+ib\omega \end{array}\right|=0,$$E44

where

$$\Psi =-\frac{M_2}{4}+(M_1-\frac{1}{4}{M}_8)\alpha +(2M_1-\frac{1}{2}{M}_8-\frac{M_6}{4}){\alpha }^{\prime }+(\frac{3}{2}{M}_4-\frac{3}{4}{M}_9)\beta$$E45

.

By replacing $$\beta$$ with $$\alpha$$ and $${\alpha }^{\prime }$$ and include fluid's shear constitutive relations, equation (33) will be reduced to an identical scheme of equation (1f).

$$\rho V^2=-\frac{M_2}{4}+[(M_1+M_3+M_4)(\alpha +2{\alpha }^{\prime })-\frac{M_6}{4}{\alpha }^{\prime }-\frac{1}{4}(M_8+M_9)(\alpha +2{\alpha }^{\prime })]$$E46

For longitudinal waves propagating perpendicular to uniaxial loading, assume rock sample is subject to a uniaxial pressure along axis 2 and the P-waves transmit along axis 1. The "large" strains are expressed as

$${\epsilon }_{ij}=\left(\begin{array}{l}{\alpha }^{\prime }00\\ 0\alpha 0\\ 00{\alpha }^{\prime }\end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}\beta 00\\ 0\beta 0\\ 00\beta \end{array}\right)$$E47

The small dynamic strain induced by plane P-waves propagating along axis 1 is

$${\epsilon }_{ij}=\left(\begin{array}{l}S{e}^{i[\omega t-k{x}_1(1+{\alpha }^{\prime })]}00\\ 000\\ 000\end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}{S}^{\prime }{e}^{i[\omega t-k{x}_1(1+\beta )]}00\\ 000\\ 000\end{array}\right)$$E48

Therefore,

$$\left|\begin{array}{l}{\Upsilon }_1{k}^2-{\rho }_{11}{\omega }^2+ib\omega ({\Upsilon }_2+M_4\beta )k^2-{\rho }_{12}{\omega }^2-ib\omega \\ ({\Upsilon }_2+M_4{\alpha }^{\prime })k^2-{\rho }_{12}{\omega }^2-ib\omega {\Upsilon }_3{k}^2-{\rho }_{22}{\omega }^2+ib\omega \end{array}\right|=0$$E49

where

$${\Upsilon }_1=M_1+(M_1+\frac{1}{2}{M}_2+3M_5+M_8)\alpha +(6M_1+\frac{1}{2}{M}_2+6M_5+M_8){\alpha }^{\prime }+(\frac{3}{2}{M}_4+3M_{10})\beta$$E50

,

$${\Upsilon }_2=\frac{1}{2}{M}_4+(\frac{1}{2}{M}_9+M_{10})\alpha +(\frac{1}{2}{M}_4+\frac{1}{2}{M}_9+2M_{10}){\alpha }^{\prime }+(3M_{11}+\frac{1}{2}{M}_4)\beta$$E51

and

$${\Upsilon }_3=M_3+(\frac{1}{2}{M}_4+M_{11})\alpha +(M_4+2M_{11}){\alpha }^{\prime }+(9M_7+7M_3)\beta$$E52

.

The two frequency-dependent P velocities can be solved from equation (37).

If solid/fluid relative motion is neglected, replace fluid strain $$\beta$$ with $$\alpha$$ and$${\alpha }^{\prime }$$,

$$\begin{array}{l}\rho V^2=M_1+M_3+M_4\\ +[(M_1+M_3+M_4)+\frac{1}{2}{M}_2+3(M_5+M_7+M_{10}+M_{11})+(M_8+M_9)]\alpha \\ +[6(M_1+M_3+M_4)+\frac{1}{2}{M}_2+6(M_5+M_7+M_{10}+M_{11})+(M_8+M_9)]{\alpha }^{\prime }.\end{array}$$E53

Equation (38) is identical to equation (1d).

For transverse waves propagating along uniaxial loading and vibrating perpendicular to loading, the static strain is expressed as equation (35) and the dynamic strain is

$${\epsilon }_{ij}=\left(\begin{array}{l}S{e}^{i[\omega t-k{x}_2(1+\alpha )]}00\\ 000\\ 000\end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}{S}^{\prime }{e}^{i[\omega t-k{x}_2(1+\beta )]}00\\ 000\\ 000\end{array}\right)$$E54

Then

$$\left|\begin{array}{l}\Psi k^2-{\rho }_{11}{\omega }^2+ib\omega -{\rho }_{12}{\omega }^2-ib\omega \\ -{\rho }_{12}{\omega }^2-ib\omega -{\rho }_{22}{\omega }^2+ib\omega \end{array}\right|=0,$$E55

where

$$\Psi =-\frac{M_2}{4}+(M_1-\frac{M_2}{2}-\frac{1}{4}{M}_8)\alpha +(2M_1+\frac{M_2}{2}-\frac{1}{2}{M}_8-\frac{M_6}{4}){\alpha }^{\prime }+(\frac{3}{2}{M}_4-\frac{3}{4}{M}_9)\beta$$E56

.

A similar operation for the reduction to pure solid's acoustoelasticity leads to

$$\rho V^2=-\frac{M_2}{4}+[(M_1+M_3+M_4)(\alpha +2{\alpha }^{\prime })-\frac{M_2}{2}(\alpha -{\alpha }^{\prime })-\frac{M_6}{4}{\alpha }^{\prime }-\frac{1}{4}(M_8+M_9)(\alpha +2{\alpha }^{\prime })]$$E57

which is identical to equation (1e).

For transverse waves propagating perpendicular to uniaxial loading and vibrating perpendicular to uniaxial loading, the dynamic strains are expressed as

$${\epsilon }_{ij}=\left(\begin{array}{l}S{e}^{i[\omega t-k{x}_3(1+{\alpha }^{\prime })]}00\\ 000\\ 000\end{array}\right),{\zeta }_{ij}=\left(\begin{array}{l}{S}^{\prime }{e}^{i[\omega t-k{x}_3(1+\beta )]}00\\ 000\\ 000\end{array}\right)$$E58

Therefore

$$\left|\begin{array}{l}\Psi k^2-{\rho }_{11}{\omega }^2+ib\omega -{\rho }_{12}{\omega }^2-ib\omega \\ -{\rho }_{12}{\omega }^2-ib\omega -{\rho }_{22}{\omega }^2+ib\omega \end{array}\right|=0,$$E59

where

$$\Psi =-\frac{M_2}{4}+(M_1+\frac{M_2}{2}-\frac{1}{4}{M}_8-\frac{M_6}{4})\alpha +(2M_1-\frac{M_2}{2}-\frac{1}{2}{M}_8){\alpha }^{\prime }+(\frac{3}{2}{M}_4-\frac{3}{4}{M}_9)\beta$$E60

.

The reduction operation leads to

$$\rho V^2=-\frac{M_2}{4}+[(M_1+M_3+M_4)(\alpha +2{\alpha }^{\prime })+\frac{M_2}{2}(\alpha -{\alpha }^{\prime })-\frac{M_6}{4}{\alpha }^{\prime }-\frac{1}{4}(M_8+M_9)(\alpha +2{\alpha }^{\prime })]$$E61

which is identical to equation (1g).

Equations (19, 24, 29, 33, 37, 40, 43) are the central result of this paper, which are applicable in rock tests.

4. Numerical tests for poro-acoustoelasticity

The elastic and density coefficients of a virtual water-saturated porous rock sample are listed in table 5. The 2nd-order elastic constants and the density coefficients come from the actual sample of water-saturated Coldlake sandstone [45]. The seven 3rd-order elastic constants are assumed for numerical tests.

4.1. Wave speeds for sample under hydrostatic loading

If the rock sample is open-pore jacketed and subject to a hydrostatic pressure$$P_h$$, and the fluid pressure is hold on one atmosphere (can be zero in approximation), we have the Biot's constitutive relations as $$P_h=Ae+\frac{2}{3}Ne+Q\zeta$$ and$$0=Qe+R\zeta$$. In the virtual tests, the hydrostatic pressure changes from 0 to 30MPa. The solid and fluid's strains are directly solved. By substituting $$\alpha =e/3$$ and $$\beta =\zeta /3$$ into equation (19) and equation (24), open-pore jacketed wave speeds are calculated.

For closed-pore jacketed rock sample in hydrostatic loading, the Biot's constitutive relations come to$$P_s=Ae+\frac{2}{3}Ne+Q\zeta$$, $$P_f=Qe+R\zeta$$, $$P_h=P_s+P_f$$and$$P_f/\varphi =K_f\zeta$$. Hydrostatic loading changes from 0 to 30MPa. By solving solid and fluid's strain and substituting into equation (19) and equation (24), "closed pore" wave speeds are calculated.

The waves center frequency ranges from 100 to 1000 KHz. The numerical results for open- and closed-pore jacketed rock sample under hydrostatic loading are shown in figures 1~2.

Figure 1 shows the relations between wave velocity and attenuation and confining pressure. Fast P wave velocity, Fast P wave inverse quality factor and S wave velocity are sensitive to the change in hydrostatic pressure. Fast P wave attenuation is significantly increased in loading process. Slow P wave's velocity in open-pore jacketed test reaches a maximum value around 13MPa. Moreover, S wave attenuation is entirely not affected by hydrostatic pressure (it holds on $$2.1569\times {10}^{-4}$$ at 1MHz), and slow P attenuation has very small change in relation to loading, therefore both have not been drawn in figure 1.

$$M_1$$	8.566 GPa	$${\rho }_{11}$$	2111 Kg/ m^3
$$M_2$$	-11.704 GPa	$${\rho }_{12}$$	-348 Kg/ m^3
$$M_3$$	0.7285 GPa	$${\rho }_{22}$$	697 Kg/ m^3
$$M_4$$	2.7 GPa	$$\varphi$$	0.335
$$M_5$$	244.3 GPa	$$k$$	1×10^-12 m^2
$$M_6$$	131.2 GPa	$$K_b$$	2.5 GPa
$$M_7$$	7.0 GPa
$$M_8$$	-449.8 GPa
$$M_9$$	-149.9 GPa
$$M_{10}$$	65.1 GPa
$$M_{11}$$	24.1 GPa

Table 5.

able 5.The coefficients of the virtual rock sample.

Comparing with open-pore jacketed test, closed-pore jacketed test has lower fast P velocity (2380~2580m/s) and higher S and slow P velocity (1177~1284m/s and 52.5~61.7m/s respectively). Fast P inverse quality factor changes from 0.5$$\times {10}^{-6}$$ to 4.6$$\times {10}^{-6}$$, which are much lower than open-pore jacketed tests (0.5$$\times {10}^{-6}$$~4.5$$\times {10}^{-5}$$). The sign of the predicted $$1/Q$$ is positive. The $$1/Q$$ value is below$${\text{10}}^{-4}$$, which agrees with the description of dissipation of traditional Biot theory [17, 18] and is obviously lower than the level of the attenuation which can be caused by local fluid flow mechanism. Local fluid flow is not considered in this study, and according to Dvorkin et al. 1994 [46], the magnitude of $$1/Q$$ produced by local fluid flow may be 1~2 orders magnitude higher than the one produced by Biot friction. Pressure is mainly confined on solid skeleton in open-pore jacketed configuration, so solid's nonlinear acoustic wave feature dominates the whole rock's speeds. Pore water is more likely to be compressed in closed-pore jacketed test and its finite deformation undertakes part of pressure, therefore water's acoustoelasticity should be considered for closed-pore jacketed case. Actually, water's velocity- pressure slope is much lower than solid, therefore, closed-pore jacketed test will have lower fast P velocity than open-pore jacketed test since part of pressure is consumed by pore water instead of by solid matrix. The higher slow P velocity and lower inverse quality factors physically means solid/fluid coupling effect is strengthened, while Biot dissipation is weakened.

Velocities and inverse quality factors vs pressure and attenuation are drawn in figure 2. For both cases, it is obvious that fast P wave and S wave are more sensitive to the changes in confining pressure than the changes in wave frequencies, while slow P waves are opposite. As to inverse quality factors, only fast P wave is sensitive to both the confining pressure and the wave frequency. It reaches peak value at the highest center frequency and the largest pressure. S wave's inverse quality factor slightly increases when center frequency increases and shows no response to the loading variations.

Figure 1.

The results in hydrostatic loading with central frequency at 1M Hz. (a) Fast P wave velocity (dashed line for open-pore jacketed, solid line for closed-pore jacketed), (b) Fast P wave inverse quality factor, (c) S wave velocity, (d) Slow P wave velocity.

Comparing open-pore jacketed results with closed-pore jacketed results shows that S wave inverse quality factors of the two cases are entirely equal and not affected by loading. Slow P inverse quality factors of the two cases have very small difference. S and slow P velocities of closed-pore jacketed test are always slightly higher than open-pore jacketed test. Fast P velocity of open-pore jacketed test is significantly higher than closed-pore jacketed test as confining pressure increases. Fast P inverse quality factor in open-pore jacketed test is lower than closed-pore jacketed test when confining pressure is below 5MPa, but gets much higher than closed-pore jacketed test when pressure reaches 30MPa.

Figure 2.

The velocities and inverse quality factors vs pressure and centre frequency in hydrostatic loading (red faces for open-pore jacketed, blue faces for closed-pore jacketed). (a) Fast P wave velocity, (b) Fast P wave inverse quality factor, (c) S wave velocity, (d) S wave inverse quality factor, (e) Slow P wave velocity, (f) Slow P wave inverse quality factor.

Mechanisms of "local fluid flow" have not been taken into account in this work. The pores in the rock are assumed to be equant when we derive the nonlinear acoustical wave equations for solid/fluid composite, which is also a basic assumption in Biot's theory [35]. In real rocks which contain both equant pores and soft pores, the contributions of wave-induced local fluid flow may be important for the explanation of the observed attenuation in rock tests. The produced attenuation in this study is mainly associated to the traditional Biot friction and may be underestimated in comparison with an actual laboratory measurement. In the limitation of Biot's theory, the prediction of dissipation can be appropriate only in the case that local fluid flow can be neglected.

4.2. Wave speeds for open-pore jacketed sample under uniaxial loading

In this subsection a virtual test is analyzed with uniaxial loading configuration. The uniaxial pressure is$$P_u$$. The fluid pressure is hold around zero. Biot's constitutive relations are written as$$P_u=Ae+2N{e}_{11}+Q\zeta$$, $$0=Ae+2N{e}_{22}+Q\zeta$$, and$$0=Qe+R\zeta$$. Let uniaxial pressure change from 0 to 30MPa, $$e_{11}$$, $$e_{22}$$and $$\zeta$$ are solved out. Substitute$$\alpha =e_{11}$$, $${\alpha }^{\prime }=e_{22}$$and $$\beta =\zeta /3$$ into equations (29, 33, 37, 40, 43). Wave center frequency changes from 100 to 1000 KHz.

Figure 3a shows that fast P velocity changes from 2380m/s to 2690m/s along loading direction and changes from 2383m/s to 2389m/s perpendicular to loading direction, while figure 3c shows slow P wave velocity increases slightly in transverse direction and decreases slightly in parallel direction in the loading process. On the same loading level, fast P inverse quality factor is higher in parallel direction than in transverse direction when loading pressure exceeds 5MPa, and both are much lower than in hydrostatic test. The $$1/Q$$ vs pressure curve seems to have a trough around 6MPa along loading and around 8MPa perpendicular to loading. The sign of the predicted $$1/Q$$ is positive and the value is below$${\text{10}}^{-4}$$.

Numerical results for S waves are shown in figure 3d. All S velocities in three directions are sensitive to the changes of confining pressure. S wave transmitting along loading direction has the highest speed comparing with two other directions, changing from 1177m/s to 1228m/s. S wave whose propagating and vibrating direction are both perpendicular to loading has relatively lower speed changing from 1177m/s to 1224m/s. S wave vibrating along uniaxial loading has the lowest speed changing from 1177m/s to 1184m/s. All three S velocities in uniaxial configuration are lower than the open-pore jacketed hydrostatic test, where S speed changes from 1177m/s to 1280m/s. For all the three directions in uniaxial test and for the unique direction in hydrostatic test, S waves share the same inverse quality factor independent to the pressure. It only changes with wave centre frequency.

As has been shown in figure 1b and figure 3b, the predicted attenuation increases with confining pressure in average, which disagrees with the experimental results in [47]. Local fluid flow, which dominates the wave attenuation and dissipation phenomena in the rocks under low confining pressure, is depressed since the flat throats and soft pores tend to close under higher confining pressure, so that the attenuation in actual experiment generally decreases with effective stress. The mechanism of local fluid flow and its relation to the confining pressure has not been analyzed in this study. The numerical results are reasonable in condition that the rock pores are assumed to be equant and the local fluid flow mechanism can be neglected.

Figure 3.

The open-pore jacketed results under uniaxial loading with central frequency at 1M Hz. (a) Fast P wave velocity (dashed line for along loading, solid line for perpendicular to loading), (b) Fast P wave inverse quality factor, (c) Slow P wave velocity, (d) S wave velocity (direction 1: $$\parallel P_u$$propagating, $$\perp P_u$$vibrating; direction 2: $$\perp P_u$$propagating, $$\parallel P_u$$vibrating; direction 3: $$\perp P_u$$propagating, $$\perp P_u$$vibrating).

4.3. Wave speeds for closed-pore jacketed sample under uniaxial loading

If water-saturated rock sample is closed-pore jacketed and subject to uniaxial pressure. In process of loading, because compressed pore water is not allowed to flow out from rubber jacket, the fluid phase will impose an extra stress increment on solid matrix in transverse directions. The six poroelastic constitutive relations are$$P_{su}=Ae+2N{e}_{11}+Q\zeta$$, $$P_f=Qe+R\zeta$$, $$P_u=P_{su}+P_f$$, $$P_{s\text{v}}=Ae+2N{e}_{22}+Q\zeta$$, $$0=P_{sv}+P_f$$and$$P_f/\varphi =K_f\zeta$$, where $$P_{su}$$ and $$P_{s\text{v}}$$ respectively denote solid stress components along the loading direction and perpendicular to loading.

Numerical results for P waves are figured in figures 4a~4c. First particular feature in this configuration is that fast P speed in vertical direction decrease slightly ( from 2385 to 2345 m/s ) when confining pressure increases, while all results in former three configurations show opposite trends. Fast P velocity changes from 2380 to 2660m/s along loading, slightly lower than the open-pore jacketed uniaxial test. As is shown in figure 4b, another particular feature is that fast P inverse quality factors along uniaxial loading decreases with rise of pressure. Fast P inverse quality factor in vertical direction increases with loading and is much higher than in parallel direction ( 0.5~5.5 $$\times {10}^{-6}$$ vs 0.5~3.5 $$\times {10}^{-6}$$ ), while in open-pore jacketed test fast P inverse quality factor in vertical direction is lower than in parallel direction ( 5.2~2.4 $$\times {10}^{-7}$$ vs 0.5~3.2 $$\times {10}^{-6}$$ ).

Figure 4.

The closed-pore jacketed results under uniaxial loading with central frequency at 1M Hz. (a) Fast P wave velocity (dashed line for along loading, solid line for perpendicular to loading), (b) Fast P wave inverse quality factor, (c) Slow P wave velocity, (d) S wave velocity (direction 1: $$\parallel P_u$$propagating, $$\perp P_u$$vibrating; direction 2: $$\perp P_u$$propagating, $$\parallel P_u$$vibrating; direction 3: $$\perp P_u$$propagating, $$\perp P_u$$vibrating).

Numerical results for S waves in closed-pore jacketed sample under uniaxial loading are shown in figure 4d. S wave transmitting along loading has the highest speed ranging from 1177m/s to 1231m/s. S velocities in two vertical directions respectively range from 1177m/s to 1225m/s and from 1177m/s to 1186m/s. S-velocities in all three directions increase as loading rises. S velocities of closed-pore jacketed uniaxial test are very close to open-pore jacketed uniaxial tests. In all test configurations of this paper, S inverse quality factors seem to share the same value which is only slightly dependent to wave centre frequency.

5. Experimental data

5.1. Physical setup

Two hydrostatic loading tests (one for "open-pore jacketed" configuration, another for "closed-pore jacketed" configuration) are performed on a sandstone sample with moderate porosity (13.26 percent) and low permeability (1.21mD). The sample is collected from a gas reservoir in southwest China. It is from the depth of around 2000 meters from surface. The sandstone sample is mainly constructed by quartz and feldspar. Minor clays and rock fragments reside inside pores and grains. It is moderately sorted. The grain size ranges from below 0.1mm to around 1 mm. The pore size ranges from below 0.1 mm to around 0.4 mm. Most grains contact well, so as to form a rigid solid skeleton. The average grain density is 2.659g/cm³. The grain bulk modulus is 39.0 GPa.

The experimental setup consists of a digital oscilloscope and a pulse generator. In the test, the rock sample is jacketed with a rubber tubing to isolate it from the confining pressure. The receiving transducer is connected to the digitizing board in the PC through a signal amplifier. A pore fluid inlet in the endplate allows passage of pore fluid through the sample and can help to control the pore pressure inside rocks in experiments. We give a small modification on the original inlet instrument by adding a valve (it is closely connected to the inlet), so that the "closed-pore jacketed" configuration can be realized. The open-pore and closed-pore tests are performed respectively. In each test, the rock is full saturated with water, and both confining pressure and pore pressure is raised to 10 MPa before the P-wave speed measurements. In the "closed-pore jacketed" test, we keep the valve closed and raises the confining pressure from 10 to 62 MPa with an interval of around 4 MPa. The P-wave speed is measured and recorded in each pressure level. In the "open-pore jacketed" test, the valve is kept open in measuring process. The pore pressure is close to atmospheric pressure when confining pressure is not so high (it can be neglected comparing with the confining loading).

5.2. Discussions on the theory and experimental data

The relationship between the measured velocities of P waves and the confining pressure in the "open-pore jacketed" and "closed-pore jacketed" tests are shown in figure 5. In the both cases, measured P velocity increases as the confining pressure gets higher. Comparing with the closed-pore test, the velocity-pressure relationship in the open-pore test has a higher slope in a range of 10-35 MPa if it is fitted with a liner correlation, and the difference of measured velocities between open-pore and closed-pore tests increases as the confining pressure rises. However, when confining pressure increases in a higher range of 40-62 MPa, this difference decreases as the confining pressure rises. The observed velocity-pressure results of open-pore test will approach to the results of closed-pore test.

Based on the poro-acoustoelasticity theory, 11 elastic constants have to be determined if we want to predict the velocity-pressure relationships in actual rock. The four 2nd order elastic constants can be theoretically calculated with a Biot-Gassmann theory [48]. The dry bulk modulus can be estimated through Pride's frame moduli expression for consolidated sandstone [49] as$$K_d=K_s(1-\varphi )/(1+c\varphi )$$, where$$K_d$$, $$K_s$$and $$\varphi$$ are dry rock modulus, grain bulk modulus and porosity. $$c$$are the consolidation parameter. We choose $$c=8$$ since it is a moderate porosity and good consolidated sandstone. According to Table 3, the 4 2nd order elastic constants can be calculated as$$M_1=37.\text{4}GPa$$, $$M_{\text{2}}=\text{-53}\text{.6}GPa$$, $$M_3=0.273GPa$$and$$M_{\text{4}}=\text{1}\text{.84}GPa$$. As to the 7 3rd-order elastic constants, more prior information about each component's higher-order elastic modulus should be fully realized if we want to give an applicable theoretical estimate. Moreover, some uniaxial loading tests and S velocity measurements should also be performed in the "open/closed-pore jacketed" configurations, so that all the 7 3rd-order elastic constants can be precisely determined and discussed. Nevertheless, these tests also require further improvements on our current experimental instruments. We will leave the issue for a future consideration.

Figure 5.

Comparison between experimental data and theoretical curves.

In this paper, we use 4 3rd order elastic constants ($$M_{\text{5}}$$, $$M_{\text{6}}$$, $$M_{\text{8}}$$and$$M_{\text{10}}$$) to fit the measured results in the open-pore and closed-pore jacketed tests, and neglect other 3 terms, since these four terms are considered as the most important for actual rocks as has been discussed in section 2.

Based on the numerical approach in Section 4, we perform numerical calculation and draw the velocity-pressure curves in figure 5. They are both in good agreement with open-pore and closed-pore measurements in a low effective pressure range. We choose$$M_{\text{5}}=1800GPa$$, $$M_{\text{6}}\text{=}-2500GPa$$, $$M_{\text{8}}\text{=2300}GPa$$and $$M_{\text{10}}\text{=100}GPa$$ in calculation.$$M_{\text{5}}$$, $$M_{\text{6}}$$and $$M_{\text{8}}$$ are mostly associated with the solid skeleton's higher order elastic characteristics. The order of magnitude of these constants (KGPa) agrees well with the measurement of 3rd-order elastic constants on rocks performed by Winkler et. al. [15, 16]. As a complement to the pure solid's 3rd-order elastic constant theory, the new 3rd elastic constant $$M_{\text{10}}$$ is related to the higher-order coupling effect between solid and fluid. Its order of magnitude is around 100GPa. It effects on the wave speed expressions for fluid-saturated rocks, and holds a lower magnitude than$$M_{\text{5}}$$, $$M_{\text{6}}$$and$$M_{\text{8}}$$, which has been discussed in Section 2. A higher value of $$M_{\text{10}}$$ will produce unreasonable results for wave speeds.

The observed differences between "open-pore jacketed" and "closed-pore jacketed" tests decrease as we still increase the confining pressure in a high pressure range (higher than 40 MPa), which do not agree with our theoretical analysis. Actually, the "open-pore jacketed" measured speeds approach to the "closed-pore jacketed" results in high pressure end. In high confining pressure range, the theoretical results fit well with the "closed-pore" measurement, while give an estimate higher than the actually measured speeds in the "open-pore" tests. As confining pressure is getting higher for an open-pore system, the stress difference between solid skeleton and pore fluid is consequently increased. The effective pressure is also increased. Pore fluid continues to flow out and the soft pores such as grain cracks will tend to close. In a high effective pressure range, the pore system of rock in a "open-pore" test will tend to act like a theoretical "closed-pore" system, where soft pores are closed and fluid stops flowing out. The basic assumption of a theoretical "open-pore" configuration breaks down in this case.

In a real sandstone, at higher pressure (above 50 MPa), the increase in the values of wave velocities essentially saturates [50, 51]. This saturation effect is caused exclusively by gradual closing of the soft crack-like fraction of the rock-frame porosity. This crack-like fraction strongly dominates the nonlinearity of the material in the pressure range of several tens MPa, whereas visible pores are less affected. The compressibility of crack-like pores strongly depends on the fact whether the filling fluid is confined inside or can be squeezed out. Thus below 50MPa the difference between "open-pore jacketed" and "closed-pore jacketed" tests is obvious until the most part of the cracks become closed. The application of this theory is limited to the pressure range lower than 50MPa. Especially for the "open-pore jacketed" case, the theoretical predictions will be unreliable when the effective pressure exceeds the limitation.

6. Conclusions

The nonlinear acoustic wave propagation equations are derived in this paper. For hydrostatic and uniaxial loading configurations of rock tests, the analytical deformation-dependent P- and S- speed expressions are presented. We discuss the reductions from current theoretical approach to the traditional works of acoustoelasticity in pure solid material. Neglecting solid/fluid differences both in a large static part of strain induced by confining pressure and in a small dynamic part of strain induced by wave's propagation and vibration, the velocity equations in this paper (equations 21, 26, 31, 34, 38, 41 & 44) are completely compatible with the expressions of acoustoelasticity theory in pure solid material (equations 1a~g).

We design a virtual rock sample with seven assumed 3rd-order elastic constants. We perform four numerical tests on this virtual rock sample with different considerations: open-pore jacketed sample under hydrostatic loading, open-pore jacketed sample under uniaxial loading, closed-pore jacketed sample under hydrostatic loading and closed-pore jacketed sample under uniaxial loading. Numerical results show that fast P wave's inverse quality factor is both sensitive to the confining pressure and wave centre frequency, which agrees with the laboratory observations of our former studies [47] in the aspect that P wave's inverse quality factors are more sensitive to pressure than S wave's. For the virtual rock sample, velocities of fast P waves and S waves seem to be more sensitive to pressure change than to changes in wave centre frequency, while slow P wave shows an opposite feature. In closed-pore jacketed uniaxial loading configuration, fast P wave's velocity in transverse direction is in decreasing trend with the rise of loading pressure, which is completely contrary to other three configurations. In all test configurations of this paper, S inverse quality factors share the same value which is independent to confining loading. S inverse quality factor can be influenced by wave centre frequency.

We perform the "open-pore jacketed" and "closed-pore jacketed" P-velocity measurements on a sandstone sample under hydrostatic loading. The observed P-velocity differences between the "open-pore jacketed" and "closed-pore jacketed" tests are getting larger as the confining pressure is increased in the range of 10~35 MPa, which agrees well with theoretical prediction, while in a higher confining loading range (>40MPa), the tested "open-pore jacketed" velocity approaches to the "closed-pore jacketed" velocity. When hydrostatic confining pressure is increased and pore pressure keeps constant in an "open-pore" experiment, fluid will continue to flow out and microscopic connections between pores tend to close. When effective pressure reaches a high level, the rock in an "open-pore" test will approach to the theoretical configuration of the "closed-pore" system, where the pore fluid will actually stop flowing out under loading. The basic assumptions of the theoretical "open-pore jacketed" configuration is violated. A rock under high effective pressure in an "open-pore jacketed" test is actually "closed". In a high confining loading, the crack-like fraction will become gradually closed while the visible pores remain almost affected. Mostly the equant pores remain survived for which the compressibility weakly depend on the presence of the fluid (does not matter "closed-" or "open-pore"). The theory presented in this paper can give reasonable predictions for wave speeds in a confining pressure range lower than 50 MPa.

Moreover, in this work, the mechanism of "local fluid flow" have been neglected when we theoretically derive the nonlinear wave equations. The heterogeneity in pore structures and fluid distributions will induce local fluid flow if elastic waves squeeze the porous rock, which surely will effect on the waves' propagation velocities and intrinsic attenuation [23, 36, 50, 52]. However, to include the mechanism of local fluid flow in a poro-acoustoelasticity theory will produce too complicated expressions to be handled on the present stage, and there probably will be several additional 2nd order terms [52, 53] and much more 3rd order terms in the strain energy expression. The present paper is limited to the area of applicability of the Biot approach, whereas rocks in nature usually contain the stiff and the soft pores. It may be important and essential that the local fluid flow mechanism should be analyzed in a poro-acoustoelasticity context, since the effects of confining pressure on pore shapes are also believed to be relevant to the local fluid flow and will consequently change the observed waves' velocities and attenuation in laboratory, which can be described in detail by introducing additional 2nd and 3rd order terms. This extension will be a future study.

Acknowledgement

Discussions with Winkler K. W. were helpful. Wu X. Y. and Hao Z. B. helped in designing the experimental setup for a "closed-pore jacket" rock test and performing the measurements. This research was sponsored by the Major State Basic Research Development Program of China (973 Program, 2007CB209505), the CNPC 12-5 Basic Research Plan (2011A-3601) the Natural Science Foundation of China (41074087, 41104066) and the RIPED Youth Innovation Foundation (2010-A-26-01).