Anisotropic viscoacoustic wave modelling in VTI media using frequency-dependent complex velocity

Under the conditions of acoustic approximation and isotropic attenuation, we derive the pseudo-and pure-viscoacoustic wave equations from the complex constitutive equation and the decoupled P-wave dispersion relation, respectively. Based on the equations, we investigate the viscoacoustic wave propagation in vertical transversely isotropic media. The favourable advantage of these formulas is that the phase dispersion and the amplitude dissipation terms are inherently separated. As a result, we can conveniently perform the decoupled viscoacoustic wavefield simulations by choosing different coefficients. In the computational process, a generalised pseudo-spectral method and a low-rank decomposition scheme are adopted to calculate the wavenumber-domain and mixed-domain propagators, respectively. Because low-rank decomposition plays an important role in the simulated procedure, we evaluate the approximation accuracy for different operators using a linear velocity model. To demonstrate the effectiveness and the accuracy of our method, several numerical examples are carried out based on the new pseudo- and pure-viscoacoustic wave equations. Both equations can effectively describe the viscoacoustic wave propagation characteristics in vertical transversely isotropic media. Unlike the pseudo-viscoacoustic wave equation, the pure-viscoacoustic wave equation can produce stable viscoacoustic wavefields without any SV-wave


Introduction
The ubiquitous seismic anisotropy is deemed as an important property that has an effect on wave propagation and reversetime migration. Generally, there are several numerical algorithms to achieve the wavefield simulation in anisotropic media. A typical approach is to directly solve the full elastic wave equations and then decouple the wavefield into P-and S-wave components. However, the wavefield separation in heterogeneous anisotropic media is very complicated and extremely computational inefficient (Yan & Sava 2009;Xu & Zhou 2014;Wang et al. 2016).
Given this problem, some scholars propose to simulate the P-wave propagation using a pseudo-acoustic (PSA) wave equation, which is derived from the acoustic approximation (Alkhalifah 1998(Alkhalifah , 2000Du et al. 2007;Fowler et al. 2010). Rao & Wang (2019) explicitly derive the numerical dispersion relation and stability condition for the PSA wave equation. Compared Fifth, three numerical examples are adopted to prove the applicability and the validity of our new approach during viscoacoustic wavefield simulation. Finally, conclusions are summarised based on analyses and experiments.

Pseudo-viscoacoustic wave equations
Based on acoustic approximation, the stress-strain relationship for VTI elastic media can be symbolically written as where is the density, ij and e ij are the components of stress and strain tensor, V P is the velocity along symmetry axis and and are the Thomsen (1986) anisotropy parameters. To simulate the attenuating property of viscoacoustic wave propagation, a frequency-dependent complex velocity is introduced in the frequency domain (Aki & Richards 2002;Yang & Zhu 2018a): Taking the isotropic quality factor into account, we get the following relation in attenuating media based on equations (1-4): in which Equation (5) has a similar structure to that of PSA wave equation in transversely isotropic (TI) media. After some transformations, we can obtain the equivalent 2D PSVA wave equations as follows: where H and V are the stresses in the horizontal and vertical directions, respectively. Based on equations (9) and (10), we can resolve the PSVA wavefield extrapolation in attenuating media. However, the modelling results of PSVA wave still contain some undesired SV-wave components. In addition, the stability of these two equations is seriously restricted by the condition of ≥ (Chu et al. 2011;Xu & Liu 2018). When ignoring the phase dispersion and the amplitude dissipation terms of the complex velocity, equations (9) and (10) then degenerate into the PSA wave equations.

Pure-viscoacoustic wave equation
Another effective scheme to implement the acoustic wave propagation is solving the decoupled P-wave dispersion relation (Tsvankin 1996). For the attenuating media, the complex-valued P-wave phase velocity can be briefly written as follows (Appendix A): 2 V 2 ( ) = D 11 sin 2 + D 33 cos 2 + [(D 11 sin 2 − D 33 cos 2 ) 2 + 4D 2 13 sin 2 cos 2 ] 1∕2 , where is the phase angle. After some simplifications, equation (11) can be reformulated as follows: where k x = k sin is the wavenumber in the x direction and k z = k cos is the wavenumber in the z direction. When applying the first-order Taylor series expansion to approximate the square-root in equation (12) and ignoring the influence of on denominator, we attain the following PUVA wave equation (Chu et al. 2011;Zhan et al. 2012;Yan & Liu 2016) After that, the PUVA wavefield extrapolation in VTI media can be effectively implemented by solving equation (13).

Numerical implementation
With respect to the amplitude dissipation term isgn( ) Q , we adopt the following approximation (Yang & Zhu 2018a) isgn( ) where k = (k x , k z ) is the wavenumber vector. Equation (14) can be conveniently computed by using a generalised pseudospectral method. With respect to the dispersion term 2 Q ln | 0 |, the wavenumber and the velocity are coupled in the natural logarithm, which cannot be straightforwardly computed by Fourier transform (FT). Given this issue, a low-rank decomposition method is introduced to approximate this operator. After that, the mixed-domain propagator can be decomposed into a separable representation (Engquist & Ying 2009;Fomel et al. 2013): where W(x, k m ) is a submatrix of W(x, k) that is composed of selected columns related to wavenumbers, W(x n , k) is another submatrix that comprises selected rows associated with spatial positions and a mn is the middle matrix coefficients. The lowrank approximation is highly accurate. Generally, the parameters M = N = 3 can satisfy the requirement of accuracy well.
To illustrate the mechanism of low-rank decomposition, we take D 11 2 H x 2 in equation (9) as an example: where t denotes the first-order time derivative, F and F −1 denote forward and inverse FTs, respectively. Using a low-rank decomposition method to approximate −k 2 Analogously, the other terms in equations (9), (10) and (13) also can be calculated by using a generalized pseudo-spectral method and a low-rank decomposition scheme (Appendix B). For each time circulation, solving equations (9) and (10) requires four forward FTs and (4 + 2N) inverse FTs, and solving equation (13) requires two forward FTs and (6 + 3N) inverse FTs. Note that N is a parameter of low-rank decomposition.

Accuracy analyses
(1) Accuracy analyses for the approximation of ≈ |k|V P During the numerical implementation, ≈ |k|V P is adopted to solve PSVA and PUVA wave equations. However, the approximation has some limitations in anisotropic media. Therefore, the accuracy analyses for this approximation are necessary.
In an anistropic medium, the relationship between angular frequency and wavenumber is written as: where v( ) is the phase velocity. Under the acoustic assumption, the phase velocity in VTI media can be expressed as (Thomsen 1986): and  In the complex stiffness coefficients D 11 , D 13 and D 33 , (1 + 2 Q ln | 0 |) is related to the phase dispersion, and isgn( ) Q associates with the amplitude dissipation. For the phase dispersion term, we have Using reference velocity V P to replace phase velocity v( ) in anisotropic media would produce some errors, which can be defined as follows: Based on equation (22), we calculate the errors with different reference velocities, phase angles and anisotropy parameters, shown in figure 1a (for convenience, the medium is elliptically anisotropic). It can be seen that the approximation error is very small. .
After importing the approximation ≈ |k|V P , the error of equation (23) can be defined as: Ignoring the effect of imaginary unit, the error becomes 1 ). Figure 1b displays the corresponding error variations with reference velocities, phase angles and anisotropy parameters in an elliptically anisotropic medium. Compared with the value of equation (23), the error is negligible.
Based on these analyses, we can see that the approximation of ≈ |k|V P is limited in anisotropic media, but the error is relatively small and thus acceptable under the assumption of weak anisotropy and attenuation.
(2) Accuracy analyses for low-rank decomposition To reduce the numerical error caused by computational algorithms, we jointly adopt the pseudo-spectral method and the low-rank decomposition scheme to solve equations (9), (10) and (13), respectively. For the sake of evaluating the accuracy of low-rank approximation, we design a simple model with linear velocity, which varies from 1500 to 3000 m·s −1 . Moreover, the wavenumbers of k x and k z are between 0 and 0.314 m −1 , and the reference frequency is set to 0 = 1000 rad s −1 . Figure 2 shows the approximations of k 2 x ln |

Numerical examples
In this section, three numerical examples are carried out using our proposed PSVA and PUVA wave equations. In each example, we set 0 = 1000 rad·s −1 . A Ricker signal with main frequency of 20 Hz is used to produce waves. The parameters  M = N = 3 are adopted in the low-rank approximation. To eliminate the artefact boundary reflections, a modified hybrid absorbing boundary condition is applied during the wavefield simulation (Liu & Sen 2010. For convenience, we first give some abbreviations of different numerical methods in Table 1.

A homogeneous model
Our first example is for a 2D homogeneous model with grids of 301 × 301. The grid intervals in the xand z-directions are 10 m. The velocity along the symmetry axis is 1500 m s −1 . The Thomsen anisotropy parameters are = 0.20 and = 0.15, respectively. The quality factor is Q = 20. The record length and the time step are 1.2 and 0.001 s, respectively. The source wavelet is placed at the centre of the model. The reference solution in the homogeneous model is directly computed by using the generalised pseudo-spectral method. Figure 3 shows several wavefield snapshots and their differences compared with the reference solutions. Unlike the acoustic modelling wavefields (figure 3a and e), the viscoacoustic modelling wavefields (figure 3b, c, f and g) have a large amount of energy loss. Besides, the wavefields computed by PSA (figure 3a) and PSVA (figure 3b and c) wave equations contain some strong SV-wave artefacts around the source, while they are effectively eliminated in the wavefields of PUA (figure 3e) and PUVA (figure 3f and g). In addition, the differences (figure 3d and h) between the reference solutions (figure 3b and f) and our new modelling results (figure 3c and g) are very small, which illustrates that our proposed algorithm is highly accurate. Because the phase dispersion and the amplitude dissipation are naturally decoupled in the complex velocity, we can conveniently implement different viscoacoustic behaviours by choosing different coefficients in equations (9), (10) and (13). For example, figure 4a and b show four wavefields at 0.7 s, calculated by the anisotropic acoustic, dispersion-effect, dissipationeffect and viscoacoustic wave equations, respectively. Compared with the acoustic wavefield, the dispersion-effect wavefield has little amplitude variation, but presents a delayed phase. In contrast, the dissipation-effect wavefield only occurs amplitude

A two-layered model
The second example is a two-layered model shown in figure 6. The computational area is divided into 301 × 301 grids with a grid spacing of 10 m. The simulated length is 1.2 s with the time step 0.001 s. The interface is located at the depth of 1600 m. Figure 7 displays several snapshots at 0.6 s calculated by using the PSA, PSVA, PUA and PUVA wave equations, respectively. Meanwhile, figure 8 plots the corresponding seismic records and extracted traces at the distance of 1000 m or at the time of 500 ms, respectively. Compared with acoustic modelling wavefields (figure 7a and c) and records (figure 8a and c), viscoacoustic simulated wavefields (figure 7b and d) and records (figure 8b and d) present obvious amplitude dissipation. Besides, seismograms from PSA and PSVA wave equations contain a large amount of SV-wave noise, while those from PUA and PUVA wave equations are completely free of SV-wave noise.

A heterogeneous model
To verify the accuracy and the robustness of our new schemes, a part of BP 2007 tilted TI (TTI) model with a constant dip angle = 0 • is used to perform the viscoacoustic wavefield simulation. The model is discretised by 501 × 361 with a grid spacing of 15 m. The record length is 1.5 s. The time step of 0.001 s is applied in the wave propagation. some model parameters, including the velocity v along the symmetry axis, the quality factor Q = 14 × (v∕1000) 2.2 , and the Thomsen anisotropy parameters and , respectively. The source signal, represented by a Ricker wavelet, located at (3750, 2175 m), is used to generate the waves. Figures 10 and 11 show several snapshots and seismograms for the modified BP 2007 TTI model computed by using the PSA, PSVA, PUA and PUVA wave equations, respectively. Compared with the acoustic modelling wavefield, the viscoacoustic modelling wavefield contains some amplitude dissipation effects along the wavefront. This phenomenon can be obviously observed from the extracted traces at the distance of 6000 m (figure 11e-h). In addition, we also compare the computational efficiency between PSVA and PUVA wave equations in a personal computer (Intel Core i7-8850H 2.60 GHz). The elapsed times of PSVA and PUVA wave equations are 511.66 and 713.11 s, respectively.

Conclusions
We have derived the new PSVA and PUVA wave equations with isotropic quality factor from the complex constitutive wave equation and the decoupled P-wave dispersion relation. Since phase dispersion and amplitude dissipation terms are inherently separated in these formulas, we can conveniently realise the decoupled viscoacoustic wavefield extrapolation. During the numerical implementation, the generalised pseudo-spectral method and the low-rank decomposition scheme were adopted to solve the wavenumber-domain and mixed-domain propagators, respectively. Meanwhile, we also examined the accuracy of low-rank approximations for different propagators. Several numerical experiments were performed to prove the validity and the robustness of our approach. Moreover, we compared the modelling results between PSVA and PUVA wave equations. Unlike the PSVA wave equation, our PUVA wave equation could obtain stable viscoacoustic wavefield without any SV-wave artefacts. mixed-domain propagators. As for equations (9) and (10), each term can be solved by the following equations: Based on the complex stiffness coefficients and the definition of quality factor, we get