For wave propagation in prestressed media, there exist three configurations involved: the stress-free natural configuration as a reference state, the initial configuration of finite static deformations as a prestressed state, and the final configuration associated with wave-induced dynamic deformations. The theory of acoustoelasticity superposes small-amplitude dynamic deformations of a wave onto a prestressed finite deformation. The relevant acoustoelastic equation of motion can be derived either in the natural or in the initial states.

As defined by Pao et al. (1984) for a predeformed medium, the dynamic displacement u ( ξ , t ) from the initial to the final state satisfies the following acoustoelastic equation of motion in the natural frame of reference, A α β γ δ ∂ 2 u γ ∂ ξ β ∂ ξ δ = ρ ∂ 2 u α ∂ t 2   (1) with A α β γ δ = T β δ i δ α γ + Γ α β γ δ (2) where   δ α γ represent the Kronecker symbol,   Γ α β γ δ   and   T β δ i   denote the fourth-order stiffness tensor and the second Piola-Kirchhoff stress tensor caused by a finite static deformation, respectively. The mass density   ρ   refers to the natural state. The fourth-order stiffness tensor   A α β γ δ   can be expanded about the material’s elastic constants in the initial state, A α β γ δ = c α β γ δ ( 1 − e η η ) + c α β γ δ ϵ η e ϵ η + c ϵ β γ δ ∂ u α i ∂ ξ ϵ + c α ϵ γ δ ∂ u β i ∂ ξ ϵ + c α β ϵ δ ∂ u γ i ∂ ξ ϵ + c α β γ ϵ ∂ u δ i ∂ ξ ϵ (3) where   e ϵ η = 1 2 ( ∂ u ϵ i ∂ ξ η + ∂ u η i ∂ ξ ϵ )   is the infinitesimal strain, u i ( ξ ) denotes the static displacement caused by the applied static loading of material points for the natural state, and   c α β γ δ   and c α β γ δ ϵ η   are the second-order and third-order elastic constants of the medium in the prestressed configuration, respectively. The coefficients   A α β γ δ   depend on both the material constants and the prestress strains. These coefficients possess the same symmetry (Pao et al., 1984) as the Hookean stiffness tensor   c α β γ δ , namely, A α β γ δ = A β α γ δ = A α β δ γ = A γ δ α β (4)

For convenience, we use the following Voigt’s compressed notation for the tensorial indices to contract the indices of   c α β γ δ (α, β, γ, δ = 1, 2, 3) to   c p q   (p, q = 1, 2, …, 6) whereby.

11→1, 22→2, 33→3, 23→4, 13→5, and 12→6. (5)

The second-order tensor   c α β γ δ   is replaced hereafter with a 6 by 6 matrix   c p q   as follows (Pao et al., 1984), [ c 11 c 12 c 13 0 0 0 c 21 c 22 c 23 0 0 0 c 31 c 32 c 33 0 0 0 0 0 0 c 44 0 0 0 0 0 0 c 55 0 0 0 0 0 0 c 66 ] (6)

The medium has 9 independent second-order elastic constants because of the symmetry of   c p q   and 20 third-order elastic constants with the following symmetry, { c 111 = c 222 = c 333 , c 144 = c 255 = c 366 , c 112 = c 223 = c 133 = c 113 = c 122 = c 233 c 155 = c 244 = c 344 = c 166 = c 266 = c 355 c 123 ,   a n d   c 456 (7)

Compared to the linear elasticity with two elastic constants for an isotropic medium, the third-order nonlinear elasticity invokes three additional elastic constants (sometimes called A, B, and C) (Green, 1973). Following Winkler and Liu (1996), the static strain tensor caused by static loading with respect to the natural state can be expressed as a function of u i ( ξ ) , ε α β = 1 2 ( u α , β i + u β , α i + u δ , α i u δ , β i ) (8)

The elastic strain energy density is expressed as W = μ ε α β 2 + ( K 2 − μ 3 ) ε δ δ 2 + A 3 ε α β ε α δ ε β δ + B ε α β 2 ε δ δ + C 3 ε δ δ     3 (9) where μ is the shear modulus, K   is the bulk modulus, and (A, B, C) are the third-order elastic constants. We see that there are two sources of nonlinearity in Eqs. 8, 9, as described by Pau and Vestroni (2019), one from the geometrical nonlinearity term in the strain tensor and the other from the cubic terms referred to physical nonlinearity through the potential energy.

For infinitesimally small strains, the geometrical nonlinearity term in Eq. 8 and the cubic terms in Eq. 9 are negligible, leading to the standard equation of linear elasticity. The (A, B, C) nomenclature used above is common, but it is not universal. Green (1973) presents a table that correlates five definitions of the third-order elastic constants used by different authors. With the Goldberg, 1961, the third-order elastic constants   c p q r   can be represented by { c 111 = c 333 = 2 C + 6 B + 2 A , c 133 = c 113 = 2 C + 2 B , c 155 = c 355 = c 111 − c 113 4 = B + A / 2 (10)

Based on the above third-order elastic constants and the Voigt abbreviated symbol in Eq. 5, the second-order effective elastic constants   A α β γ δ   for a 2D homogeneous medium can be reduced to { A 11 = ( λ + 2 μ ) ( 1 + 3 e 11 − e 33 ) + ( 6 B + 2 C + 2 A ) e 11 + ( 2 B + 2 C ) e 33 , A 13 = λ ( 1 + e 11 + e 33 ) + ( 2 B + 2 C ) ( e 11 + e 33 ) , A 33 = ( λ + 2 μ ) ( 1 − e 11 + 3 e 33 ) + ( 6 B + 2 C + 2 A ) e 33 + ( 2 B + 2 C ) e 11 , A 51 = 2 ( B + A 2 ) e 13 + 2 ( λ + 2 μ ) e 13 , A 53 = 2 ( B + A 2 ) e 13 + 2 ( λ + 2 μ ) e 13 , A 55 = μ ( 1 + e 11 + e 33 ) + ( B + A 2 ) ( e 11 + e 33 ) (11) where   e α β     are the components of the total prestrain tensor, and   λ   and   μ   are the Lamé constants in the initial configuration.

The plane-wave analysis of acoustoelastic equations provides a simple way to understand the physics of wave propagation in acoustoelastic media. In this section, we follow the plane-wave analysis in Section 3.2 of Pao et al. (1984), with attempts to compare with experimental measurements (Winkler and Liu, 1996). The resulting analytical solution will be used to validate our numerical results.

Let us assume that the displacement components in Eq. 1 can be described by a time-harmonic plane wave, u α = U α ⁡ exp [ i ( k l β ξ β − ω t ) ] (12) where   ω   is the angular frequency, l β   denotes the component of the wave normal (a unit vector),   k is the wave number with the wave velocity v = ω / k , U α is the amplitude constant, and i = − 1 . SubstitutingEq. 12 into (1), we obtain a system of equations for the amplitude vector U, [ A α β γ δ l β l δ − ρ v 2 δ α γ ] U γ = 0 (13)

Considering the case of orthotropic predeformations, the orthotropic body possesses three axes of symmetry, called the principal axes of orthotropy. For a plane wave propagating in the direction of the ξ 3 axis, Eq. 13 can be simplified to [ A α 3 γ 3 − ρ v 2 δ α γ ] U γ = 0 (14)

Equation 14 has a solution under the condition that the determinant of coefficients is zero, namely, det [ A α 3 γ 3 − ρ v 2 δ α γ ] = 0 (15)

The corresponding characteristic equation in the contracted notation is, det [ A 55 − ρ v 2 A 54 A 53 A 45 A 44 − ρ v 2 A 43 A 35 A 34 A 33 − ρ v 2 ] = 0 (16) where the second-order effective elastic constants are defined by Pao et al. (1984) as { A 33 = T 33 i + c 33 ( 1 + 2 e 33 ) + c 331 e 11 + c 332 e 22 + c 333 e 33 , A 44 = T 33 i + c 44 ( 1 + 2 e 22 ) + c 441 e 11 + c 442 e 22 + c 443 e 33 , A 55 = T 33 i + c 55 ( 1 + 2 e 11 ) + c 551 e 11 + c 552 e 22 + c 553 e 33 , A 34 = A 43 = c 33 ( ∂ u 2 i / ∂ ξ 3 ) + c 44 ( ∂ u 3 i / ∂ ξ 2 ) + 2 c 344 e 23 , A 35 = A 53 = c 33 ( ∂ u 1 i / ∂ ξ 3 ) + c 55 ( ∂ u 3 i / ∂ ξ 1 ) + 2 c 355 e 31 , A 45 = A 54 = c 44 ( ∂ u 1 i / ∂ ξ 2 ) + c 55 ( ∂ u 2 i / ∂ ξ 1 ) + 2 c 456 e 12 (17)

Because ξ 3 is in the direction of one of the strains, we have A 43 =   A 53 = 0 (due to   e 31 = e 32 = 0 ), Eq. 16 reduces to   ( A 33 − ρ v 2 ) [ ( ρ v 2 ) 2 − ( A 55 + A 44 ) ρ v 2 + A 55 A 44 − A 45 2 ] = 0 (18)

The solution to this equation for three eigenvalues can be obtained in terms of the following two cases. For the cubic-symmetric medium under the confining prestress, the resulting three eigenvalues lead to one longitudinal mode and two transverse modes as follows, { v P 2 = A 11 / ρ , v S V 2 = A 44 / ρ , v S H 2 = A 44 / ρ (19)

For the symmetric plane, for example when θ ₂ = 0°, the velocities can be expressed as (Zong, 2014) { v P 2 = ( A 11 + A 44 + K ) / ρ , v S V 2 = ( A 11 + A 44 − K ) / ρ , v S H 2 = A 44 / ρ (20) where K = 4 A 11 2 ⁡ sin 4 θ 1 − 8 A 11 A 44 ⁡ sin 4 θ 1 − 4 A 12 2 ⁡ sin 4 θ 1 − 8 A 12 A 44 ⁡ sin 4 θ 1 − 4 A 11 2 ⁡ sin 2 θ 1 + 8 A 11 A 44 ⁡ sin 2 θ 1 + 4 A 12 2 ⁡ sin 2 θ 1 + 8 A 12 c 44 ⁡ sin 2 θ 1 + ( A 11 − A 44 ) 2 (21)

Figure 1 shows the definition of angles θ ₁, θ ₂, and θ ₃, where θ ₁ is the angle between the z-axis and the line linking a point in the space and the origin of coordinates, θ ₂ is the angle between the x-axis and the projection line on the xoy plane, and θ ₃ is the angle between the y-axis and the projection line on the xoy plane. These angles define the relationship between the reference-coordinate and orthotropy axes which should be coincident for the cases of confining stress, uniaxial stress, and pure shear stress, but with a difference of 45° for the case of simple shear stress. We see that these modes are strongly related to stress changes. The stress-induced velocity anisotropies depend on the orientation of external stresses.

Based on the experimental measurements (Winkler and Liu, 1996) under the hydrostatic pressures for Portland sandstone, we calculate velocities for the longitudinal mode against pressure. The physical constants for the simulation are picked from Tables III-IV of (Winkler and Liu, 1996). The prestress condition in the simulation is regarded as the confining stress with its implementation described in detail in Section 4.1. The result is shown in Figure 2, where the circles (experimental measurements) are obtained by picking the pixel coordinates in Figure 4 of (Winkler and Liu, 1996). We see a weak nonlinearity for velocity variations with increasing pressures, resulting from the third-order elasticity with a cubic term for the strain energy function. Some departures from experimental measurements can be observed during the pressure of 20–60 MPa possibly because of the deformation of pores in the sandstone which the third-order elasticity for solids cannot account for. In general, for most consolidated rocks, the third-order elasticity could be enough to describe the nonlinear stress dependence of velocity variations.

Equation 1 can be written as the first-order velocity-stress equation, { ρ v β , t = τ α β , α , τ α β , t = A α β γ δ v γ , δ (22) where   v   and   τ   are velocity and stress, respectively. With Eq. 11 and the Cartesian tensor notation with the x- and z-coordinates for convenience, we rewrite Eq. 22 as its velocity-stress FD formulation, { ρ v x , t = τ x x , z + τ x z , z , ρ v z , t = τ x z , x + τ z z , z , τ x x , t = A 11 v x , x + A 13 v z , z + A 15 ( v x , z + v z , x ) , τ z z , t = A 13 v x , x + A 33 v z , z + A 35 ( v x , z + v z , x ) , τ x z , t = A 15 v x , x + A 35 v z , z + A 55 ( v x , z + v z , x ) (23)

The nonlinear effect of predeformation is introduced by the effective elastic constants. Equation 23 can be reduced to the classical elastic wave equation when no static stress is applied.

The SSG-FD method is widely used in seismology to simulate elastic wave propagation. The method defines velocity components, stress components, and physical properties in four grids crosswise shown in Figure 3A, whereas the RSG-FD method defines velocity components in one grid and stress components and elastic parameters in another grid shown in Figure 3B. The RSG-FD method is employed in this study to solve the first-order velocity-stress acoustoelastic equations.

As shown in Figure 3B, the RSG technique calculates the spatial derivative along the diagonal direction of grids, and then interpolates the result along the normal coordinate axis to obtain the spatial derivative along the horizontal and vertical directions. We rotate the direction of spatial derivatives from the horizontal ( x ) and vertical ( z ) directions to the diagonal directions ( x ¯ and z ¯ ) as { x ¯ = Δ x Δ r x − Δ z Δ r z , z ¯ = Δ x Δ r x + Δ z Δ r z , (24) where   Δ r = Δ x 2 + Δ z 2 . The first-order spatial derivatives along the horizontal and vertical directions become, { ∂ ∂ x = Δ r 2 Δ x ( ∂ ∂ z ¯ + ∂ ∂ x ¯ ) , ∂ ∂ z = Δ r 2 Δ z ( ∂ ∂ z ¯ − ∂ ∂ x ¯ ) . (25) With Eq. 25, it is easy to define the differentiation operators   D x ¯     and   D z ¯   , which implement the spatial derivatives along diagonal directions ( x ¯ and z ¯ ) in the time domain (Saenger et al., 2000), { D x ¯ u ( x , z , t ) = 1 Δ r [ u ( x − Δ x 2 , z + Δ z 2 , t ) − u ( x + Δ x 2 , z − Δ z 2 , t ) ] , D z ¯ u ( x , z , t ) = 1 Δ r [ u ( x + Δ x 2 , z + Δ z 2 , t ) − u ( x − Δ x 2 , z − Δ z 2 , t ) ] . (26) With Eqs. 25, 26, the numerical differentiation operators are obtained that perform the spatial derivatives along the x and z directions in the rotated staggered grid by a linear combination of the derivatives along the x ¯ and z ¯ directions, { ∂ ∂ x u ( x , z , t ) ≈ Δ r 2 Δ z ( D x ¯ u ( x , z , t ) + D z ¯ u ( x , z , t ) ) , ∂ ∂ z u ( x , z , t ) ≈ Δ r 2 Δ z ( D x ¯ u ( x , z , t ) − D z ¯ u ( x , z , t ) ) . (27)

A detailed process of the implementation of the RGS-FD scheme is addressed as follows.

By that analogy, we express the L-order RSG of each physical quantity in the discrete form as v x , x ( x , z ) ≈ ∑ l = 1 L / 2 c n 2 Δ x { v x ( x + ( l − 1 / 2 ) Δ x , z + ( l − 1 / 2 ) Δ z ) + v x ( x + ( l − 1 / 2 ) Δ x , z − ( l − 1 / 2 ) Δ z ) − v x ( x − ( l − 1 / 2 ) Δ x , z − ( l − 1 / 2 ) Δ z ) − v x ( x − ( l − 1 / 2 ) Δ x , z + ( l − 1 / 2 ) Δ z ) } v z , x ( x , z ) ≈ ∑ l = 1 L / 2 c n 2 Δ x { v z ( x + ( n − 1 / 2 ) Δ x , z + ( n − 1 / 2 ) Δ z ) + v z ( x + ( n − 1 / 2 ) Δ x , z − ( n − 1 / 2 ) Δ z ) − v z ( x − ( n − 1 / 2 ) Δ x , z − ( n − 1 / 2 ) Δ z ) − v z ( x − ( n − 1 / 2 ) Δ x , z + ( n − 1 / 2 ) Δ z ) } v x , z ( x , z ) ≈ ∑ l = 1 L / 2 c n 2 Δ z { v x ( x + ( l − 1 / 2 ) Δ x , z + ( l − 1 / 2 ) Δ z ) − v x ( x + ( l − 1 / 2 ) Δ x , z − ( l − 1 / 2 ) Δ z ) − v x ( x − ( l − 1 / 2 ) Δ x , z − ( l − 1 / 2 ) Δ z ) + v x ( x − ( l − 1 / 2 ) Δ x , z + ( l − 1 / 2 ) Δ z ) } v z , z ( x , z ) ≈ ∑ l = 1 L / 2 c n 2 Δ z { v z ( x + ( l − 1 / 2 ) Δ x , z + ( l − 1 / 2 ) Δ z ) − v z ( x + ( l − 1 / 2 ) Δ x , z − ( l − 1 / 2 ) Δ z ) − v z ( x − ( l − 1 / 2 ) Δ x , z − ( l − 1 / 2 ) Δ z ) + v z ( x − ( l − 1 / 2 ) Δ x , z + ( l − 1 / 2 ) Δ z ) } τ x x , x ( x + Δ x / 2 , z + Δ z / 2 ) ≈ ∑ l = 1 L / 2 c n 2 Δ x { τ x x ( x + l Δ x , z + l Δ z ) + τ x x ( x + l Δ x , z − ( l − 1 ) Δ z )   − τ x x ( x − l Δ x , z − ( l − 1 ) Δ z ) − τ x x ( x − l Δ x , z + l Δ z ) } τ x z , x ( x + Δ x / 2 , z + Δ z / 2 ) ≈ ∑ l = 1 L / 2 c n 2 Δ x { τ x z ( x + l Δ x , z + l Δ z ) + τ x z ( x + l Δ x , z − ( l − 1 ) Δ z )   − τ x z ( x − l Δ x , z − ( l − 1 ) Δ z ) − τ x z ( x − l Δ x , z + l Δ z ) } τ z z , x ( x + Δ x / 2 , z + Δ z / 2 ) ≈ ∑ n = 1 L / 2 c n 2 Δ x { τ z z ( x + n Δ x , z + n Δ z ) + τ z z ( x + n Δ x , z − ( n − 1 ) Δ z )   − τ z z ( x − n Δ x , z − ( n − 1 ) Δ z ) − τ z z ( x − n Δ x , z + n Δ z ) } τ x x , z ( x + Δ x / 2 , z + Δ z / 2 ) ≈ ∑ n = 1 L / 2 c n 2 Δ z { τ x x ( x + n Δ x , z + n Δ z ) − τ x x ( x + n Δ x , z − ( n − 1 ) Δ z )   − τ x x ( x − n Δ x , z − ( n − 1 ) Δ z ) + τ x x ( x − n Δ x , z + n Δ z ) } τ x z , z ( x + Δ x / 2 , z + Δ z / 2 ) ≈ ∑ l = 1 L / 2 c n 2 Δ z { τ x z ( x + l Δ x , z + l Δ z ) − τ x z ( x + l Δ x , z − ( l − 1 ) Δ z )   − τ x z ( x − l Δ x , z − ( l − 1 ) Δ z ) + τ x z ( x − l Δ x , z + l Δ z ) } τ z z , z ( x + Δ x / 2 , z + Δ z / 2 ) ≈ ∑ l = 1 L / 2 c n 2 Δ z { τ z z ( x + l Δ x , z + l Δ z ) − τ z z ( x + l Δ x , z − ( l − 1 ) Δ z )   − τ z z ( x − l Δ x , z − ( l − 1 ) Δ z ) + τ z z ( x − l Δ x , z + l Δ z ) } (28) where the differential coefficients   c l   will be specified subsequently. For the time derivatives, we use the second-order FD approximation with   v i ( ( l + 1 / 2 ) Δ t ) = v i l + 1 / 2   for velocity, where Δ t is the time step.

Substituting Eq. 28 into (23) and adding a source field in the simulation, we obtain the discrete form of the first-order velocity-stress acoustoelastic equations, { v ˙ x = ( ∂ x τ x x + ∂ z τ x z − f x ) / ρ , v ˙ z = ( ∂ x τ x z + ∂ z τ z z − f z ) / ρ , τ ˙ x x = A 11 ∂ x v x + A 13 ∂ z v z + A 15 ( ∂ z v x + ∂ x v z ) + f ˙ x x , τ ˙ z z = A 13 ∂ x v x + A 33 ∂ z v z + A 35 ( ∂ z v x + ∂ x v z ) + f ˙ z z , τ ˙ x z = A 15 ∂ x v x + A 35 ∂ z v z + A 55 ( ∂ z v x + ∂ x v z ) + f ˙ x z , (29) where the cap dot represents the first-order partial derivative of time, and f   denotes the components of the external body force. The set of equations can be combined into a matrix system as v ˙ + s = Mv (30) where the product Mv is given by Mv = [ 0 0 A 11 ∂ x + A 15 ∂ z A 13 ∂ x + A 35 ∂ z A 15 ∂ x + A 55 ∂ z 0 0 A 13 ∂ z + A 15 ∂ x A 33 ∂ z + A 35 ∂ x A 35 ∂ z + A 55 ∂ x ∂ x / ρ 0 0 0 0 0 ∂ z / ρ 0 0 0 ∂ z / ρ ∂ x / ρ 0 0 0 ] [ v x v z τ x x τ z z τ x z ] , (31) and the source vector is s = [ f x / ρ f z / ρ f ˙ x x f ˙ z z f ˙ x z ] T (32)

For the time integration from t to t + Δ t , with the discretizing time as t = n Δ t , we have, { v x n + 1 = v x n + Δ t ( ∂ x τ x x n + ∂ z τ x z n − f x n ) / ρ , v z n + 1 = v z n + Δ t ( ∂ x τ x z n + ∂ z τ z z n − f z n ) / ρ , τ x x n + 1 = τ x x n + Δ t ( A 11 ∂ x v x n + A 13 ∂ z v z n + A 15 ( ∂ z v x n + ∂ x v z n ) + f x x n ) , τ z z n + 1 = τ z z n + Δ t ( A 13 ∂ x v x n + A 33 ∂ z v z n + A 35 ( ∂ z v x n + ∂ x v z n ) + f z z n ) , τ x z n + 1 = τ x z n + Δ t ( A 15 ∂ x v x n + A 35 ∂ z v z n + A 55 ( ∂ z v x n + ∂ x v z n ) + f x z n ) . (33)

The set of Eq. 33 can be combined into a matrix system as. v n + 1 = M n v n + s n .

The conventional PML absorption boundary is unable to handle the reflections (evanescent waves) at grazing incidences, especially for low-frequency components. The CPML absorption boundary improves the absorption effect of evanescent waves, with efficient computations by introducing auxiliary variables to avoid convolutional calculations that need to store past-time wavefields. The implementation of the CPML boundary does not need to split velocity and stress fields, and is easily incorporated into the RSG-FD program.

The CPML absorption boundary with the selection of parameters is detailed in Appendix. The algorithm is developed by modifying the complex coefficient S k and introducing the auxiliary variables d k , α k , and χ k , where d k denotes the damping profile, and two real variables χ k ≥ 1 and α x ≥ 0 . Substituting Eq. A4 into Eq. A2 yields, ∂ x ˜ k = 1 S k ∂ x k = [ 1 χ k − d k χ k 2 ( i ω + a k + d k / χ k ) ] ∂ x k = 1 χ k ∂ x k + ψ ˜ k (34) where ψ ˜ k is a memory variable that can be rewritten as i ω ψ ˜ k + ( a k + d k χ k ) ψ ˜ k = − d k χ k 2 ∂ x k (35)

Using an inverse Fourier transform to Eq. 35, we have, ∂ ψ k ∂ t + ( a k + d k χ k ) ψ k = − d k χ k 2 ∂ x k (36) where ψ k is the inverse of ψ ˜ k . This equation has an iterative solution of the form, ψ k n = ψ k n − 1 e − ( a k + d k / χ k ) Δ t + d k ∂ x k χ k ( a k χ k + d k ) ( e − ( a k + d k / χ k ) Δ t − 1 ) (37)

Equation 37 is derived simply by solving a first-order differential equation, leading to the same result in Komatitsch and Martin (2007) obtained by a recursive convolution method.

The CPML absorbing boundary for the first-order velocity-stress formulation of acoustoelastic equations can be implemented by 1) applying Eqs. 27–38 for the estimation of various memory variables in the rotated operators, 2) substituting these memory variables into Eq. 34 to obtain the spatial derivatives of all the field components in the stretched coordinates, and 3) substituting the resulting spatial derivatives into the velocity–stress formulas (Eq. 22) to yield the following C-PML formulation of acoustoelastic equations, { ρ v ˙ x = 1 χ x ∂ x τ x x + ψ x , τ x x + 1 χ z ∂ z τ x z + ψ z , τ x z , ρ v ˙ z = 1 χ x ∂ x τ x z + ψ x , τ x z + 1 χ z ∂ z τ z z + ψ z , τ z z , τ ˙ x x = A 11 ( 1 χ x ∂ x v x + ψ x , v x ) + A 13 ( 1 χ x ∂ z v z + ψ z , v z ) + A 15 ( 1 χ z ∂ z v x + ψ z , v x + 1 χ x ∂ x v z + ψ x , v z ) , τ ˙ z z = A 13 ( 1 χ x ∂ x v x + ψ x , v x ) + A 33 ( 1 χ z ∂ z v z + ψ z , v z ) + A 35 ( 1 χ z ∂ z v x + ψ z , v x + 1 χ x ∂ x v z + ψ x , v z ) , τ ˙ x z = A 15 ( 1 χ x ∂ x v x + ψ x , v x ) + A 35 ( 1 χ z ∂ z v z + ψ z , v z ) + A 55 ( 1 χ z ∂ z v x + ψ z , v x + 1 χ x ∂ x v z + ψ x , v z ) . (38)

The stability criterion for the velocity–stress RSG-FD operator with equal grid spacings can be generally expressed as the following inequality (Masson et al., 2006), Δ t V max Δ r ≤ C     , (39) with C = 1 D ∑ k = 1 n | c n |     , (40) where   V max   represents the maximum phase velocity in the medium, D represents the spatial dimension, and   c n   are the differential coefficients depending on the order of the spatial operator. For an eighth-order spatial and second-order temporal 2D FD operators with coefficients c 1 = 1225 1024 ,   c 2 = − 245 3072 ,   c 3 = 49 5120 , and   c 4 = − 5 7168 , we can calculate C = 0.5497.

Numerical dispersion is one of the main factors affecting the accuracy of the finite-difference method. To reduce the level of numerical dispersions in the RSG-FD method, a discretization interval should be selected to ensure a good number of grid points n λ per minimum wavelength, which is generally calculated by n λ = V min Δ r f max , (41) where   V min is the minimum phase velocity and   f max     is the maximum frequency, usually taken as four times the center frequency of the source. Dispersion errors accumulate with increasing propagation distances, which can be reduced by refining grids. Numerical experiments by Chen et al. (2006) show that no less than 3 grid points per wavelength are required for the RSG eighth-order operator.

For this case, the stress field is isotropic under the hydrostatic stress P , leading to the same stress and strain in all the directions. The principal strain components can be expressed as, { e 11 = e 33 = − P 3 K , e 13 = 0. (43)

Substituting into Eq. 11 leads to a stiffness matrix { A 11 = λ + 2 μ + ( 2 λ + 4 μ + 4 C + 8 B + 2 A ) ( − P / 3 K ) , A 33 = λ + 2 μ + ( 2 λ + 4 μ + 4 C + 8 B + 2 A ) ( − P / 3 K ) , A 55 = μ + ( 2 μ + 2 B + A ) ( − P / 3 K ) , A 13 = λ + ( 2 λ + 4 C + 4 B ) ( − P / 3 K ) , A 51 = 0 , A 53 = 0. (44)

Wavefield snapshots are calculated with increasing hydrostatic stresses, as shown in Figure 4 for the x-component of the particle velocity at t = 0.12 ms. We see that the stress-induced velocity variations are isotropic, with the amplitude along the wavefront changing in accordance with the characteristics of a vertical force source. We calculate the theoretical wave velocity by Eqs. 20, 21. Figure 5 compares the theoretical and simulated wave velocities at the same angle but with different prestresses. Figure 6 compares the theoretical and calculated wave velocities under the same pressure (10 MPa) but with different directions. Both the theoretical and simulated wave velocities agree well under the different prestress conditions. These comparisons validate numerical simulations at least in travel time.

For this case, the stress field is anisotropic under the uniaxial stress P. The principal strain components can be expressed as { e 11 = P ( λ + μ ) μ ( 3 λ + 2 μ ) , e 33 = − P λ 2 μ ( 3 λ + 2 μ ) , e 13 =   0. (45)

The corresponding stiffness matrix becomes, { A 11 = λ + 2 μ + ( 3 λ + 6 μ + 2 C + 6 B + 2 A ) e 11 + ( 2 B + 2 C − λ − 2 μ ) e 33 , A 33 = λ + 2 μ + ( 2 B + 2 C − λ − 2 μ ) e 11 + ( 6 B + 2 C + 2 A + 3 λ + 6 μ ) e 33 , A 55 = μ + ( μ + B + A 2 ) e 11 + ( μ + B + A 2 ) e 33 , A 13 = λ + ( λ + 2 C + 2 B ) e 11 + ( 2 B + 2 C + λ ) e 33 , A 51 = 0 , A 53 = 0 , (46)

Figure 7 shows the wavefield snapshots with increasing uniaxial stresses for the x-component of the particle velocity at t = 0.12 ms. With increasing uniaxial stresses, the P-wave and S-wave are gradually coupled together, with the wavefronts of qP-wave and qS-wave becoming more elliptical due to the stress-induced anisotropy of velocities. For p = 10 MPa, the qP-wave becomes elliptical slightly, but the qS-wave is almost a circle. For p = 30 MPa, we see a trigeminal region occurring symmetrically to connect the qS-wave and the qP-wave. The oval aspect ratio of the qP-wave becomes larger, implying an increasing velocity anisotropy on the longitudinal and transverse symmetry axes. For p = 50 MPa, we see more pronounced trigeminal region, approximately rhombic qS-wave, and elliptical qP-wave. Under the anisotropic predeformed condition, the polarization direction changes with the propagation direction, but neither parallel nor perpendicular to the propagation direction. These polarized waves are usually denoted by qP and qS, with their wavefronts deforming along the x-axis as the uniaxial stress increases. The azimuthal anisotropy of qP-wave is more and more obvious than that of the qS-wave, consistent with the mechanical mechanism of uniaxial stress. We calculate the principal strain direction according to Eq. 45 and mark it in the figure. It coincides with the directions of maximum/minimum velocities, which validates our numerical simulations. Thus, the directions of the velocity principal axis can indicate the directions of the principal strain for this case.

As shown in Figure 8, the shear prestress can be classified into two cases: pure-shear and simple-shear deformations (Means and Williams, 1976). Different pre-deformations correspond to the different geological features. In order to identify these geological features in the seismic wave field, the simulations of pure shear and simple shear wavefield are carried out. For the former (see Figure 8A), the principal strain direction does not rotate with the progression of deformations, which is also called the non-rotational deformation or the coaxial deformation. In this case, the horizontal elongation and vertical contraction of deformations occur simultaneously under the uniaxial loading. The geological tensile and compressive effects produce a pure shear deformation.

The principal strain components for the pure-shear prestress become, { e 11 = − e 33 = P λ + 2 μ , e 13 = 0 , (47) yielding the following stiffness matrix, { A 11 = λ + 2 μ + ( 4 λ + 8 μ + 4 B − 2 A ) e 11   , A 33 = λ + 2 μ − ( 4 λ + 8 μ + 4 B − 2 A ) e 11   , A 55 = μ , A 13 = λ , A 53 = 0 , A 51 = 0. (48)

Wavefield snapshots are calculated with increasing pure-shear stresses and shown in Figure 9 for the x-component of the particle velocity at t = 0.12 ms. Similar to the case of uniaxial stresses, the circular wavefront at 0 MPa becomes more and more elliptical with increasing pure-shear stresses. The same changes in azimuthal anisotropy and amplitude can be seen as the uniaxial case. Unlike the uniaxial pressure condition, the absolute values of strains along the two axes of a pure-shear stress field are the same, implying that the stress-induced velocity anisotropy is more sensitive under the pure-shear prestress. Therefore, the velocities of qP-wave and qS-wave increase much more in the z-axis direction compared to the uniaxial case, as shown in Figure 7. We also calculate the principal strain direction according to Eq. 47 and mark it in the figure. It coincides with the directions of maximum/minimum velocities. Similarly, the directions of the velocity principal axis indicate the directions of the principal strain.