Efficient modeling of wave propagation in a vertical transversely isotropic attenuative medium based on fractional Laplacian

Abstract

To efficiently simulate wave propagation in a vertical transversely isotropic (VTI) attenuative medium, we have developed a viscoelastic VTI wave equation based on fractional Laplacian operators under the assumption of weak attenuation (Q ≫ 1), where the frequency-independent Q model is used to mathematically represent seismic attenuation. These operators that are nonlocal in space can be efficiently computed using the Fourier pseudospectral method. We evaluated the accuracy of numerical solutions in a homogeneous transversely isotropic medium by comparing with theoretical predictions and numerical solutions by an existing viscoelasticanisotropic wave equation based on fractional time derivatives. To accurately handle heterogeneous Q, we select several Q values to compute their corresponding fractional Laplacians in the wavenumber domain and interpolate other fractional Laplacians in space. We determined its feasibility by modeling wave propagation in a 2D heterogeneous attenuative VTI medium. We concluded that the new wave equation is able to improve the efficiency of wave simulation in viscoelasticVTI media by several orders and still maintain accuracy. INTRODUCTION The derivation of suitable wave equations to model seismic wave propagation in the earth media is a key subject in the development of seismic imaging/inversion for studying subsurface geologic properties and exploring earth resources. In past decades, attention has been paid to quantifying seismic attenuation effects during wave propagation in isotropic media (e.g., Carcione et al., 1988; Zhu and Harris, 2014; Shen et al., 2018). Like velocity anisotropy, attenuation anisotropy is recognized in laboratory and field measurements. For instance, aligned fractures often cause measurable velocity and attenuation anisotropy (Liu et al., 1993; Lynn et al., 1999; Carter and Kendall, 2006; Usher et al., 2017). Fluid-filled fractures also cause significant attenuation anisotropy (Liu et al., 2007). Guo and McMechan (2017) conduct full-waveform modeling tests to study the effects of fracture parameters on anisotropic attenuation and suggest that attenuation anisotropy is a potential tool for reservoir characterization. Therefore, the ability to model the wave behavior in the presence of anisotropic attenuation is important for further quantitative reservoir studies. An early theoretical approach to study the effective anisotropic attenuation used mechanical models (i.e., standard linear solid or Maxwell bodies) to simulate waves in viscoelastic finely layered media using Backus averaging (Carcione, 1992). By extending the viscoelastic relaxation times to the anisotropic relaxation times τij (Carcione, 1995; Komatitsch and Tromp, 1999), we can model, in theory, arbitrary attenuation anisotropy in seismic modeling. Following this, Ruud and Hestholm (2005) implement a 2D time-domain finite-difference modeling of orthorhombic attenuation, and Bai and Tsvankin (2016) present a detailed numerical demonstration and analysis of vertical transversely isotropic (VTI) attenuation by 2D time-domain finite-difference modeling. Zhu (2017) derives a viscoelastic-anisotropic wave equation based on frequency-independent Q (Kjartansson, 1979). Although his equation is able to parameterize general Q anisotropy and can be used to simulate wave simulation with arbitrarily anisotropic Q, the efficiency of his numerical computation may hamper large-scale applications. The reason is that these fractional time-derivative operators are nonlocal in time, which means that their numerical evaluation requires storing the time history of the field variables in memory. Hence, the memory required to evaluate the fractional time derivative term can be very significant, particularly in cases of a long recording length and/or 3D problems. Carcione et al. (2002) reformulate the wave equation by adding an integer to the fractional order. They show that this reformulation Manuscript received by the Editor 18 July 2018; revised manuscript received 5 November 2018; published ahead of production 21 January 2019; published online 11 March 2019. The Pennsylvania State University, Department of Geosciences, EMS Energy Institute, and Institute of Natural Gas Research, University Park, Pennsylvania 16802, USA. E-mail: [email protected] (corresponding author). Colorado School of Mines, Center for Wave Phenomena, Golden, Colorado, USA. E-mail: [email protected]. © 2019 Society of Exploration Geophysicists. All rights reserved. T121 GEOPHYSICS, VOL. 84, NO. 3 (MAY-JUNE 2019); P. T121–T131, 13 FIGS., 3 TABLES. 10.1190/GEO2018-0538.1 D ow nl oa de d 03 /1 8/ 19 to 6 6. 71 .1 10 .1 77 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / can reduce the memory storage by one order, but not degrade the accuracy. The same idea was used for accelerating the solving of the viscoelastic case (Carcione, 2009) and the Kevin-Voigt relation for medical dissipation media (Caputo et al., 2011). One approach proposed by MacDonald et al. (2015) is choosing “optimal” time steps to compute temporally weighted history that includes contributions from the entire past history. It seems that fewer past wavefields than the truncated strategy actually are stored. They claim that this optimal storage not only improves computational efficiency, but it also maintains the same or higher accuracy than the truncation method. All of the above strategies partially alleviate the computational problem with less storage. Partially or optimally storing the past history of wavefields in a realistic 3D dimension problem still poses a big challenge for seismic problems. Another interesting approach may overcome the storage of the past wavefields using the diffusive representation of fractional time derivatives (Jazia et al., 2014). So far, in our opinion, the most promising approach is to replace fractional time derivatives with fractional Laplacians, which are nonlocal in space rather than time, first suggested by Chen and Holm (2004). Computing these operators can be efficient without needing to store the time history of the field variables (Carcione, 2010; Treeby and Cox, 2010). Zhu and Harris (2014) develop a fractional Laplacian viscoacoustic wave equation for seismic modeling and demonstrate that the fractional Laplacian formulation can approximate the frequency-independent Q model very well. This strategy has been implemented for solving viscoelastic wave equations (Zhu and Carcione, 2014; Wang et al., 2018) and further used for (least-squares) reverse time imaging problems (Zhu, 2014; Zhu et al., 2014; Guo et al., 2016; Sun et al., 2016; Zhu and Sun, 2017), wave-equation migration Q analysis (Shen and Zhu, 2015), and accelerating fullwaveform inversion (Chen and Zhou, 2017; Xue et al., 2018). In this paper, we show how this idea can be extended to modeling the propagation of seismic waves in viscoelastic VTI media. First, the fractional time constitutive equations of the frequency-independent Q model are described. Second, a fractional Laplacian constitutive equation and its corresponding wave equation are derived. The discretization and implementation of the proposed wave equation using the Fourier pseudospectral method is discussed. The accuracy of the numerical solutions of the fractional Laplacian wave equation is evaluated. Finally, two synthetic models are tested to illustrate the efficiency of the developed modeling. REVIEW OF FRACTIONAL TIME CONSTITUTIVE EQUATIONS Before proceeding to a discussion of fractional Laplacian constitutive equations, it will be beneficial to briefly revisit the fractional time formulation of the anisotropic-viscoelastic wave equation (Zhu, 2017). This section will serve as the basis for later derivation. The generalized time-domain constitutive equation for the anisotropic-viscoelastic medium is expressed as (Carcione, 2014) σij 1⁄4 ψ ijkl ∂tekl; i; j; k; l 1⁄4 1; 2; 3; (1) where σij is the stress tensor, ekl is the strain tensor, ∂t is the time derivative, denotes the time convolution, and ψ ijkl are the relaxation function stiffness tensor of the medium properties. In terms of the frequency-independent Q (or constant-Q) model (Kjartansson, 1979), the anisotropic relaxation function of the frequency-independent Q model can be defined as ψ ijklðtÞ 1⁄4 Cijkl Γð1 − tÞ t t0 −2γijkl HðtÞ; (2) where Cijkl is the real (elastic) stiffness tensor, and the anisotropic γijkl can be calculated by γijkl 1⁄4 arctanð1∕QijklÞ∕π. The anisotropic Qijkl is the ratio of the real and imaginary parts of the corresponding stiffness coefficient (Carcione, 1992; Zhu and Tsvankin, 2006), Qijkl 1⁄4 Reð ~ CijklÞ∕Imð ~ CijklÞ. Assume that the symmetry of the attenuationQijkl in anisotropic media is the same as that of the stiffness tensor Cijkl. Because of the symmetry in the stiffness tensor Cijkl and Qijkl, in Voigt notation, we use a symmetric stiffness matrix CIJ and QIJ hereafter (I, J 1⁄4 1, 2, 3, 4, 5, 6). The most general anisotropic linear elastic material has 21 components. In VTI media, the elastic model is fully characterized by five elastic parameters and density. Similarly, five anelastic parameters are needed to characterize such a VTI attenuation behavior. In orthorhombic media, nine anelastic parameters are required to describe the viscoelastic behavior. By substituting the relaxation function in equation 2 into equation 1, the stress-strain equation in the Voigt’s matrix notation becomes 0 BBBBBBB@ σ11 σ22 σ33 σ23 σ13 σ12 1 CCCCCCCA 1⁄4 0 BBBBBBB@ D11 D12 D13 D14 D15 D16 D12 D22 D23 D24 D25 D26 D13 D23 D33 D34 D35 D36 D14 D24 D34 D44 D45 D46 D15 D25 D35 D45 D55 D56 D16 D26 D36 D46 D56 D66 1 CCCCCCCA × 0 BBBBBBB@ e11 e22 e33 2e23 2e13 2e12 1 CCCCCCCA ; (3) where DIJ 1⁄4 CIJωIJ 0 ∂IJ t , γIJ 1⁄4 arctanð1∕QIJÞ∕π, and ∂IJ t is the fractional order time derivative (γIJ ≪ 0.5). Stress-strain relation in 2D VTI media We consider 2D modeling of P-