Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional Laplacians

We evaluated a time-domain wave equation for modeling acoustic wave propagation in attenuating media. The wave equation was derived from Kjartansson’s constant-Q constitutive stress-strain relation in combination with the mass and momentum conservation equations. Our wave equation, expressed by a second-order temporal derivative and two fractional Laplacian operators, described very nearly constant-Q attenuation and dispersion effects. The advantage of using our formulation of two fractional Laplacians over the traditional fractional time derivative approach was the avoidance of time history memory variables and thus it offered more economic computations. In numerical simulations, we formulated the first-order constitutive equations with the perfectly matched layer absorbing boundaries. The temporal derivative was calculated with a staggered-grid finite-difference approach. The fractional Laplacians are calculated in the spatial frequency domain using a Fourier pseudospectral implementation. We validated our numerical results through comparisons with theoretical constant-Q attenuation and dispersion solutions, field measurements from the Pierre Shale, and results from 2D viscoacoustic analytical modeling for the homogeneous Pierre Shale. We also evaluated different formulations to show separated amplitude loss and dispersion effects on wavefields. Furthermore, we generalized our rigorous formulation for homogeneous media to an approximate equation for viscoacoustic waves in heterogeneous media. We then investigated the accuracy of numerical modeling in attenuating media with differentQ-values and its stability in largecontrast heterogeneous media. Finally, we tested the applicability of our time-domain formulation in a heterogeneous medium with high attenuation. INTRODUCTION Conventional seismic modeling approaches usually ignore attenuation effects on wavefields. In earth media, however, attenuation causes reduced energy and a distorted phase of seismic waves, especially when the attenuation is high due to the presence of, for example, gas pockets (Dvorkin and Mavko, 2006). To accurately characterize wave propagation in real media, attenuation and associated dispersion effects should be taken into account in seismic wavefield modeling. Over the past three decades, much work has been done to model acoustic attenuation effects during wave propagation. Basically, these forward-modeling approaches are split into two categories: One category implements attenuation in the frequency domain by allowing the velocity to be complex (Aki and Richards, 1980; Liao and McMechan, 1996; Štekl and Pratt, 1998). The other introduces Q in a time-domain wave equation. Time-domain approaches commonly use superposition of mechanical elements (e.g., Maxwell and standard linear solid [SLS] elements) to describe Q behavior, which are known as approximate constant-Q (ACQ) models (Liu et al., 1976; Emmerich and Korn, 1987; Carcione, 2007; Zhu et al., 2013). Nevertheless, the memory and computation time requirements, especially for, say, 3D applications, may limit these approaches from being widely applied for seismic migration and inversion. Instead of such ACQ models, use of a constant-Q model (Kjartansson, 1979) in the time-domain wave equation is attractive because it is accurate in producing desirable constant-Q behavior at all frequencies in the band. However, rigorous use of such a constant-Qmodel introduces a fractional time derivative, i.e., noninteger (irrational or fractional) power of the time derivative (Caputo and Mainardi, 1971; Mainardi and Tomirotti, 1997). Solving the fractional wave equation has been quite common in mathematics, acoustics, and bioengineering (e.g., Mainardi, 2010; Manuscript received by the Editor 3 July 2013; revised manuscript received 14 December 2013; published online 12 March 2014. Stanford University, Department of Geophysics, Stanford, California, USA. E-mail: tyzhu@stanford.edu; jerry.harris@stanford.edu. © 2014 Society of Exploration Geophysicists. All rights reserved. T105 GEOPHYSICS, VOL. 79, NO. 3 (MAY-JUNE 2014); P. T105–T116, 12 FIGS. 10.1190/GEO2013-0245.1 D ow nl oa de d 03 /1 2/ 14 to 1 71 .6 5. 23 8. 13 0. 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 / Treeby and Cox, 2010; Caputo et al., 2011). Recently, this fractional calculus was introduced to model seismic wave propagation (Carcione et al., 2002; Carcione, 2009, 2010). Carcione et al. (2002) and Carcione (2009) simulate acoustic and elastic waves using the constant-Q wave equation with the fractional time derivative that was computed by the Grünwald-Letnikov (GL) approximation (Podlubny, 1999). Nevertheless, it is notable that the fractional time derivative of a single variable depends on all the previous values of this variable. This property requires a large memory of stressstrain history from time t 1⁄4 0 (Carcione, 2010). Therefore, in numerical simulations we have the expense of memory resources even though in some situations it is possible to truncate the operator (Podlubny, 1999; Carcione et al., 2002). To overcome this memory expense, Lu and Hanyga (2004) introduce a set of secondary internal variables in addition to the primary internal variables to solve fractional differential equations of arbitrary order. Nevertheless, additional computational costs are still paid. Without introducing extra variables, Chen and Holm (2004) propose to use the fractional Laplacian operator to model anomalous attenuation behavior in fractional wave equations for homogeneous media. The fractional Laplacian is easily evaluated in the spatial frequency domain using the generalized Fourier approach (Carcione, 2010). As a result, the fractional Laplacian successfully avoids additional memory required by the fractional time derivative. Carcione (2010) gives a homogeneous constant-Q wave equation using the fractional Laplacian. His formulation captures amplitude loss and velocity dispersion in a single term. Separation of amplitude loss and dispersion operators in the fractional Laplacian wave equation may be preferable because separated forms are more useful in compensating for attenuation loss in inverse problems, e.g., reverse time imaging by only reversing sign of the attenuation operator and leaving the sign of the dispersion operator unchanged (Treeby et al., 2010; Zhu, 2014; Zhu et al., 2014). Zhang et al. (2010) also derive the approximate constant-Q wave equation with decoupled amplitude loss and dispersion effects, but their derivation using the normalization transform for decoupling amplitude loss and dispersion is not clearly described in the abstract. Following Treeby and Cox’s (2010) approach, here, we derive a time-domain wave equation using the fractional Laplacian to model constant-Q behavior. Starting fromKjartansson’s constant-Q constitutive stress-strain relation, we present the formulation for homogeneous media and then generalize to smoothly heterogeneous media. This fractional Laplacian based nearly constant-Q (NCQ) wave equation introduces separated amplitude loss and phase velocity dispersion operators based on the fractional Laplacian. For numerical simulations, perfectly matched layer (PML) absorbing boundaries are incorporated with the first-order constant-Q conservation equations. The spatial derivative is computed by a staggeredgrid pseudospectral method while a staggered-grid finite-difference approach is used for temporal discretization. The paper is organized as follows: The beginning is a review of the constant-Q stress-strain relation and a derivation of the fractional time derivative wave equation. Next, we provide detailed derivations of the fractional NCQ wave equation for homogeneous and heterogeneous media. Following that, we describe numerical experiments for a model of the homogeneous Pierre Shale to investigate the accuracy of our modeling by comparison with theoretical constant-Q solutions and for a two-layer model of heterogeneous medium to investigate stability for large-contrast heterogeneous media. Finally, we test the capability and efficiency of our formulation for heterogeneous media. CONSTANT-Q AND FRACTIONAL LAPLACIAN Constant-Q stress-strain relation and fractional time derivative Even though there is evidence for frequency-dependent Q from field and laboratory data (Kan et al., 1983; Sams et al., 1997; Batzle et al., 2006), in this paper we focus on a mathematical and computational frequency-independent Q model, which is considered to be a practical approximate Q model for seismology problems (McDonal et al., 1958; Aki and Richards, 1980). Kjartansson (1979) explicitly gives a linear description of attenuation that exhibits the exact constant-Q characteristic. His formulation uses fractional derivatives that are based on a relaxation function of the form t−2γ , instead of integer power of t. The relaxation function is given as ψ1−2γðtÞ 1⁄4 M0 Γð1 − 2γÞ t t0 −2γ HðtÞ; t > 0; (1) where M0 1⁄4 ρ0c0 cos2ðπγ∕2Þ is a bulk modulus, ρ0 is the acoustic density, Γ is Euler’s gamma function, and t0 is a reference time (t0 1⁄4 1∕ω0). The parameter γ 1⁄4 arctanð1∕QÞ∕π is dimensionless, and we know that 0 < γ < 0.5 for any positive value of Q. HðtÞ is the Heaviside step function. Note that only three parameters, the phase velocity c0 at a reference frequency ω0 and the Q, are needed to describe loss mechanisms, rather than several relaxation times and the complex modulus in the SLS models (Carcione, 2007). In terms of the Caputo’s fractional derivative (Caputo, 1967), the isotropic stress-strain (σ-ε) relation may be deduced in the following form: σ 1⁄4 ψ1−2γðtÞ dε dt 1⁄4 M0 t 0 ∂2γε ∂t2γ ; (2) where the symbol * denotes the convolution operator. The strain ε has the fractional order 2γ of the time derivative. According to equation 2, the behavior of the fractional time derivative ∂2γε∕∂t2γ depends on the characteristics


INTRODUCTION
Conventional seismic modeling approaches usually ignore attenuation effects on wavefields.In earth media, however, attenuation causes reduced energy and a distorted phase of seismic waves, especially when the attenuation is high due to the presence of, for example, gas pockets (Dvorkin and Mavko, 2006).To accurately characterize wave propagation in real media, attenuation and associated dispersion effects should be taken into account in seismic wavefield modeling.
Over the past three decades, much work has been done to model acoustic attenuation effects during wave propagation.Basically, these forward-modeling approaches are split into two categories: One category implements attenuation in the frequency domain by allowing the velocity to be complex (Aki and Richards, 1980;Liao and McMechan, 1996;Štekl and Pratt, 1998).The other introduces Q in a time-domain wave equation.Time-domain ap-proaches commonly use superposition of mechanical elements (e.g., Maxwell and standard linear solid [SLS] elements) to describe Q behavior, which are known as approximate constant-Q (ACQ) models (Liu et al., 1976;Emmerich and Korn, 1987;Carcione, 2007;Zhu et al., 2013).Nevertheless, the memory and computation time requirements, especially for, say, 3D applications, may limit these approaches from being widely applied for seismic migration and inversion.Instead of such ACQ models, use of a constant-Q model (Kjartansson, 1979) in the time-domain wave equation is attractive because it is accurate in producing desirable constant-Q behavior at all frequencies in the band.However, rigorous use of such a constant-Q model introduces a fractional time derivative, i.e., noninteger (irrational or fractional) power of the time derivative (Caputo and Mainardi, 1971;Mainardi and Tomirotti, 1997).
Solving the fractional wave equation has been quite common in mathematics, acoustics, and bioengineering (e.g., Mainardi, 2010;Treeby and Cox, 2010;Caputo et al., 2011).Recently, this fractional calculus was introduced to model seismic wave propagation (Carcione et al., 2002;Carcione, 2009Carcione, , 2010)).Carcione et al. (2002) and Carcione (2009) simulate acoustic and elastic waves using the constant-Q wave equation with the fractional time derivative that was computed by the Grünwald-Letnikov (GL) approximation (Podlubny, 1999).Nevertheless, it is notable that the fractional time derivative of a single variable depends on all the previous values of this variable.This property requires a large memory of stressstrain history from time t ¼ 0 (Carcione, 2010).Therefore, in numerical simulations we have the expense of memory resources even though in some situations it is possible to truncate the operator (Podlubny, 1999;Carcione et al., 2002).
To overcome this memory expense, Lu and Hanyga (2004) introduce a set of secondary internal variables in addition to the primary internal variables to solve fractional differential equations of arbitrary order.Nevertheless, additional computational costs are still paid.Without introducing extra variables, Chen and Holm (2004) propose to use the fractional Laplacian operator to model anomalous attenuation behavior in fractional wave equations for homogeneous media.The fractional Laplacian is easily evaluated in the spatial frequency domain using the generalized Fourier approach (Carcione, 2010).As a result, the fractional Laplacian successfully avoids additional memory required by the fractional time derivative.Carcione (2010) gives a homogeneous constant-Q wave equation using the fractional Laplacian.His formulation captures amplitude loss and velocity dispersion in a single term.Separation of amplitude loss and dispersion operators in the fractional Laplacian wave equation may be preferable because separated forms are more useful in compensating for attenuation loss in inverse problems, e.g., reverse time imaging by only reversing sign of the attenuation operator and leaving the sign of the dispersion operator unchanged (Treeby et al., 2010;Zhu, 2014;Zhu et al., 2014).Zhang et al. (2010) also derive the approximate constant-Q wave equation with decoupled amplitude loss and dispersion effects, but their derivation using the normalization transform for decoupling amplitude loss and dispersion is not clearly described in the abstract.
Following Treeby and Cox's (2010) approach, here, we derive a time-domain wave equation using the fractional Laplacian to model constant-Q behavior.Starting from Kjartansson's constant-Q constitutive stress-strain relation, we present the formulation for homogeneous media and then generalize to smoothly heterogeneous media.This fractional Laplacian based nearly constant-Q (NCQ) wave equation introduces separated amplitude loss and phase velocity dispersion operators based on the fractional Laplacian.For numerical simulations, perfectly matched layer (PML) absorbing boundaries are incorporated with the first-order constant-Q conservation equations.The spatial derivative is computed by a staggeredgrid pseudospectral method while a staggered-grid finite-difference approach is used for temporal discretization.
The paper is organized as follows: The beginning is a review of the constant-Q stress-strain relation and a derivation of the fractional time derivative wave equation.Next, we provide detailed derivations of the fractional NCQ wave equation for homogeneous and heterogeneous media.Following that, we describe numerical experiments for a model of the homogeneous Pierre Shale to investigate the accuracy of our modeling by comparison with theoretical constant-Q solutions and for a two-layer model of heterogeneous medium to investigate stability for large-contrast heterogeneous media.Finally, we test the capability and efficiency of our formulation for heterogeneous media.

Constant-Q stress-strain relation and fractional time derivative
Even though there is evidence for frequency-dependent Q from field and laboratory data (Kan et al., 1983;Sams et al., 1997;Batzle et al., 2006), in this paper we focus on a mathematical and computational frequency-independent Q model, which is considered to be a practical approximate Q model for seismology problems (McDonal et al., 1958;Aki and Richards, 1980).Kjartansson (1979) explicitly gives a linear description of attenuation that exhibits the exact constant-Q characteristic.His formulation uses fractional derivatives that are based on a relaxation function of the form t −2γ , instead of integer power of t.The relaxation function is given as where M 0 ¼ ρ 0 c 2 0 cos 2 ðπγ∕2Þ is a bulk modulus, ρ 0 is the acoustic density, Γ is Euler's gamma function, and t 0 is a reference time (t 0 ¼ 1∕ω 0 ).The parameter γ ¼ arctanð1∕QÞ∕π is dimensionless, and we know that 0 < γ < 0.5 for any positive value of Q. HðtÞ is the Heaviside step function.Note that only three parameters, the phase velocity c 0 at a reference frequency ω 0 and the Q, are needed to describe loss mechanisms, rather than several relaxation times and the complex modulus in the SLS models (Carcione, 2007).
In terms of the Caputo's fractional derivative (Caputo, 1967), the isotropic stress-strain (σ-ε) relation may be deduced in the following form: where the symbol * denotes the convolution operator.The strain ε has the fractional order 2γ of the time derivative.According to equation 2, the behavior of the fractional time derivative ∂ 2γ ε∕∂t 2γ depends on the characteristics of the relaxation function ψ 1−2γ ðtÞ.
Figure 1 shows how ψ 1−2γ ðtÞ varies with Q-values (0.1, 1, 10, and 100).When Q is large, the relaxation function ψ 1−2γ ðtÞ decays slowly and tends toward the step function.Consequently, the fractional time derivative ∂ 2γ ε∕∂t 2γ approaches ε.When Q is infinite (the elastic case), the current stress will be linear with the current strain σ ¼ M 0 ε.When Q tends toward zero (infinite attenuation), the function ψ 1−2γ ðtÞ approaches a one-sided delta function and thus the fractional time derivative ∂ 2γ ε∕∂t 2γ reduces to ∂ε∕∂t.When Q lies between zero and infinity, the calculation of the fractional time derivative depends on the state of previous time steps; i.e., the current stress must be calculated from the time history of the strain.As a result, numerical evaluation of a fractional time derivative operator requires large memory and computation time to store and use the time history (Carcione et al., 2002;Carcione, 2009;Caputo et al., 2011), even if it is possible to truncate its operator based on some short memory principle (Podlubny, 1999;Carcione et al., 2002).

The fractional Laplacian
To overcome this memory requirement and the associated numerical difficulties, the fractional Laplacian operator ð−∇ 2 Þ α∕2 is proposed by Chen and Holm (2004) to model power law attenuation behavior instead of the fractional time derivative (∂ α ∕∂t α , where α is real).In this case, the fractional Laplacian, calculated in the spatial domain, avoids storage of the wavefields from previous time steps.When incorporating this with a Fourier pseudospectral numerical scheme, the fractional Laplacian can be implemented without any particular numerical difficulties.
A constant-Q wave equation using the fractional Laplacian For homogeneous acoustic media, the linear first-order momentum conservation equation can be written as and the strain-velocity equation is where σ is the stress field, ε is the strain field, and υ ¼ ðυ x ; υ y ; υ z Þ is the acoustic particle velocity vector.
Combining equation 2, the stress-strain relation, with equations 5 and 6, the first-order conservation equations, the fractional wave equation with constant density ρ 0 is obtained: where ∇ 2 is the Laplacian operator and c 2 ¼ M 0 ∕ρ 0 ¼ c 2 0 cos 2 ðπγ∕2Þ.The exponent 2 − 2γ is the fractional order of the time derivative, the so-called fractional time derivative.It is evident that equation 7 reduces to the classical acoustic wave equation as γ → 0 (Q → ∞), and the diffusion equation as γ → 1∕2 (Q → 0).This fractional wave equation is also well-known as the Caputo's wave equation.Carcione et al. (2002) successfully solve the fractional time wave equation using Grünwald-Letnikow and central-difference approximations for the time discretization and Fourier method to compute the spatial derivative.Recently, Carcione ( 2009) extends its derivation to build the elastic wave equation with constant-Q.In the following, we will replace the fractional time derivative with the fractional Laplacian operator.
The dispersion relation for fractional wave equation 7 can be obtained by substituting the plane wave solution expðiωt − i k • rÞ, where k is the complex wavenumber vector, ω is the angular frequency, and r is the spatial coordinate vector: The attenuation and phase dispersion are derived in Appendix A.
After some calculations (see Appendix B), we obtain a new dispersion relation as Equation 9 approximates equation 8; thus, we refer to it as a NCQ dispersion relation.Figure 2 shows how well the attenuation and dispersion in equation 9 (dashed line) approximate to the exact constant-Q attenuation and dispersion (see Appendix A) in equation 8 (solid line) using different Q values (2, 10, and 100).The other media parameters used in Figure 2 are listed in the first example of the numerical experiments section.One can see that the differences increase with decreasing Q value, even though they match very well for relatively high Q values (low attenuation).
In homogeneous media, velocity c 0 and Q are independent of the space variable.We can apply 2D inverse spatial and temporal Fourier transform (see the notation in equation 4) to equation 9, resulting in the fractional Laplacian k2γþ2 → ð−∇ 2 Þ γþ1 and k2γþ1 → ð−∇ 2 Þ γþ1∕2 ; thus, we obtain the NCQ wave equation for homogeneous media: Note that when Q is small, the curve approximates a one-sided delta function.In contrast, when Q is large, the dependence becomes more flat (toward a step function).
A nearly constant-Q wave equation where the coefficients are Note that equation 10 essentially serves the same purpose as equation 7, the fractional time wave equation.When Q is infinite (no attenuation), i.e., γ ¼ 0, our fractional Laplacian equation 10 reduces exactly to the familiar acoustic wave equation.Numerical stability analysis of equation 10 is shown in Appendix C.

Freezing-unfreezing theory for heterogeneous media
In smoothly heterogeneous media c 0 ðrÞ and QðrÞ are assumed to vary smoothly in space.Using the locality principle (freezing-unfreezing) of the pseudodifferential operator (Stein [1993], p. 230), we now adapt the fractional NCQ wave equation for homogeneous media (equation 10) to describe propagation in heterogeneous media.The locality principle allows us to replace the constant parameters appearing in the pseudodifferential operator in equation 10 by the variable parameters.The detailed derivation is shown in Appendix D. The fractional NCQ wave equation for smoothly heterogeneous media is then obtained: where the coefficients vary in space as follows:

Numerical implementation
Modeling wave propagation using the first-order conservation equations has previously been implemented in the finite-difference time-domain approach (Zhu et al., 2013).It is easy to implement the PML boundary to absorb waves leaving the domain of interest (Berenger, 1994).The temporal derivative is solved with the staggeredgrid finite-difference approach.The first-order spatial derivatives are calculated by the staggered-grid pseudospectral method that is known to reduce spatial numerical dispersion and nonphysical ringing (Özdenvar and McMechan, 1996;Carcione, 1999).Here, the fractional Laplacian is calculated by taking the spatial Fourier transform, as defined in equation 3. Thus, in homogeneous media, the computational form of the fractional Laplacians in equation 17 is written as ð−∇ 2 Þ γ εðr; tÞ ¼ F −1 fk 2γ F½εðr; tÞg; (18) where the symbols F and F −1 are the Fourier transform and inverse Fourier transform operators, respectively.This result highlights the advantage of the derived fractional Laplacian operators; in the spatial frequency domain, the fractional Laplacian is easily calculated.Whereas in smoothly heterogeneous media (γðrÞ is slowly varying in space), equations 18 and 19 will be approximately satisfied.In our implementation, we use the average value of γðrÞ to approximate the fractional power terms γðrÞ in equations 18 and 19.

NUMERICAL EXPERIMENTS Attenuation and dispersion
To illustrate the utility of the proposed viscoacoustic wave equation to correctly model constant-Q attenuation and dispersion, we compare in the first example the numerical results  From Carcione (2009), the quality factor of Pierre Shale rock is Q ¼ 32 and the P-wave acoustic velocity is about 2164 m∕s at the reference frequency 100 Hz.For a homogeneous lossy medium, we choose a typical shale density of 2.2 g∕cm 3 .The simulations are performed in 2D with a 512 × 512 grid, 1.0-m spacing in the xand z-directions, and a time step of 0.23 × 10 −3 s.A Ricker wavelet source with 100-Hz center frequencies is located at (256 m, 256 m) in simulations.
The attenuation coefficient α and the phase velocity are calculated by the amplitude spectral ratio in decibels/meter and the phase change at each frequency (Buttkus, 2000;Treeby and Cox, 2010): where A 1;2 are the amplitude spectra at two receivers, ϕ 1;2 are phase spectra, and d is the propagation distance in meters.To keep the simulation free from numerical artifacts (e.g., numerical dispersion), we placed two receivers at 20 and 70 m from the source.The geometric spreading for 2D problem is corrected by the factor ffiffi ffi r p . Figure 3 shows the accuracy of calculated attenuation and dispersion for the Pierre Shale model.The experimental data from McDonal et al. (1958) are represented by diamond points.The solid lines in Figure 3 give the theoretical attenuation and dispersion based on equations A-4 and A-5.The open circles denote the calculated attenuation and dispersion values extracted from our numerical simulations at two receivers.They appear to fit well with the theoretical values as well as the experimental data.In Figure 3b, the calculated phase velocity drops rapidly below 40 Hz and cannot fit this portion of the theoretical results well.The difference is due to the underestimation of small attenuation values at lower frequencies by our NCQ model.Similar observations are made by Wuenschel (1965).
To illustrate the effects of separated amplitude loss and phase dispersion operators on wavefields, we conduct simulations with different partial formulations of equation 10.For example, if one only considers the phase dispersion effects and ignores the amplitude loss effects, the dispersion-dominated wave equation will be Similarly, the loss-dominated wave equation that contains the amplitude loss effects will be The second term in equation 23 describes attenuation of the constant-Q model.We assume the reference velocity c 0 ¼ 2164 m∕s is defined at high frequency, say, 1500 Hz.We simulate wavefields using the model with Q ¼ 10.A nearly constant-Q wave equation To test the accuracy of our numerical results of this NCQ wave equation over long distances (more attenuation), we increase our model size to a 2048 × 3584 grid to capture wavefields over the longer distance.The mesh spacing in the xand z-directions is 1.0 m, respectively.Other parameters are the same as the previous experiment.We compare numerical results at receivers located at 500, 1500, and 3000 m from the source with the corresponding analytical solutions for a 2D homogeneous Pierre Shale medium (Figure 5).The 2D analytical solution is obtained by convolving the constant-Q Green's function with the source function (Carcione, 2007).We evaluated the accuracy of numerical solutions using the root-mean-square (rms) errors, which is defined by where nt is the number of time samples of the seismic trace, d n j is the calculated value of the numerically simulated trace at sample j, and the superscript a denotes the corresponding analytical value.We found that the rms error deceases with offset, 1.3 × 10 −3 at 500 m (Figure 5a), 3.2 × 10 −3 at 1500 m (Figure 5b), and 7.1 × 10 −3 at 3000 m (Figure 5c).Numerical solutions agree with the analytical solutions very well at these distances.This section investigates the accuracy of the wave equation for different Q media.The velocity model is homogeneous as above.The source is located at the center of the model (0.5, 0.5) km.We place a receiver at the location of (0.5, 0.7) km. Figure 6 shows snapshots taken at 200 ms in the NCQ wave equation with Q = (a) 100, (b) 30, (c) 10, and (d) 4. We can see decreased amplitudes and delayed phases with decreasing Q values.We do not simulate for Q < 4 because the simulation is unstable at this distance.Figure 7  corresponding to the Q values in Figure 6.We found that the smaller Q value results in the larger error, which is consistent with the observation in Figure 2. The corresponding rms error values are 2.3 × 10 −4 , 1.3 × 10 −3 , 4.7 × 10 −3 , and 1.1 × 10 −1 , respectively.

Stability in large-contrast heterogeneous media
To demonstrate the stability of simulations in media with large contrasts, we design a two-layer model with P-wave velocity of 1800 m∕s in the top layer and 3600 m∕s in the bottom layer and a Q value of 30 in the top and 100 in the bottom.The interface is at a depth of 800 m.The source function is a Ricker wavelet with 25-Hz center frequency located at (968 m, 520 m) in the simulation.We discretize model in a 200 × 200 grid with 8.0 m spacings in the xand z-directions.The time step is 0.8 × 10 −3 s.The accuracy of the solutions is compared with SLS solutions, which are roughly accurate to model constant-Q behavior in simulation over a short distance, as Zhu et al. (2013) show.
We test three scenarios: (1) two-layer velocity model with homogeneous Q ¼ 30, (2) two-layer Q model with homogeneous P-wave velocity ¼ 1800 m∕s, and (3) two-layer velocity and Q models.Figure 8 shows the snapshots at 280 ms using the proposed NCQ wave equation and the viscoacoustic wave equation based on the SLS model (Zhu et al., 2013).The first row shows the effects of velocity discontinuities on the wavefields.The two numerical results agree with each other, as shown in Figure 8c.The second row shows the effects of Q discontinuities on the wavefields.One can see that the weak reflection is due to Q contrast, as shown in the second row of Figure 8.In the third row, we show simulations using velocity and Q contrasts.Overall, the solution by NCQ in Figure 8g matches that by the SLS in Figure 8h very well.These comparisons demonstrate that numerical simulations using the proposed NCQ wave equation in such large contrast velocity and Q models are stable and practically accurate.

SYNTHETIC MODEL EXAMPLE
Figure 9 shows the P-wave velocity and Q models.The density model is constant of 2.2 g∕cm 3 .Zone #4 in the model represents a high attenuation reservoir (Q ¼ 30).The model is discretized in a mesh with 200 × 600 grid points.The space is equally sampled in x and z at 0.5 m.A PML absorbing boundary with 20 grid points is implemented at four sides.An explosive source with a 200-Hz center-frequency Ricker wavelet is chosen for simulation.The time step is 0.014 ms.The source is located at a 90-m horizontal distance and at 110-m depth.Receivers are in the left well of 0.5 m horizontal distance.They are distributed from depths of 5-292.5 m with space 2.5 m.
Figure 10 shows the snapshots at 20 and 45 ms, respectively.To illustrate the attenuation effects on the wavefield, we show the acoustic cases (Figure 10a and 10b) as reference.The arrow labeled 3-4 indicates reflections from the interface of zones #3 and #4. Figure 10c and 10d show the constant-Q viscoacoustic cases.
Figure 10c shows overall reduced reflection amplitude compared with acoustic reflections.Reflections from the interface at the bottom of the high attenuation zone (arrow 4-5) suffer from significantly reduced energy because the wave propagates through zone #4.These observations are more visible in the seismograms shown in Figure 11.
To check the efficiency and accuracy of the NCQ wave equation, we also ran SLS modeling using a single SLS, which is considered to be reasonably accurate in practical applications (Zhu et al., 2013).All simulations were run on an Intel workstation (Xeon X5650 2.66 GHz).The acoustic simulation takes 226 s.The running time for the proposed NCQ algorithm and the SLS algorithm is 303 and 368 s, respectively.It may appear surprising that the present NCQ algorithm is more computationally efficient than the single SLS modeling algorithm that is considered as the most efficient viscoelastic modeling approach (Zhu et al., 2013).The reason is that the  Figure 12 shows a comparison of simulated traces between our proposed NCQ wave equation and the SLS wave equation.The two agree with each other very well in details, even later-arriving reflections.

DISCUSSION
The fractional Laplacian NCQ wave equation is an alternative to the fractional time wave equation.Theoretically, the NCQ wave equation is able to model approximately constant-Q behavior when Q is not very small, a feature that is further verified by numerical simulations.For smaller Q values, the rms errors of simulations compared with the analytical solutions become large.We found that numerical simulation becomes unstable when Q < 4. We conclude that this NCQ wave equation is highly accurate in Q environments for most seismic applications.On the other hand, the advantage of solving the fractional Laplacian wave equation is computational efficiency in comparison with the fractional time wave equation.Computing the fractional Laplacian is accomplished without using and storing the previous time wavefields required by the fractional time derivative via GL approximation (e.g., Carcione et al., 2002).A synthetic example also indicates that the proposed NCQ wave equation is more efficient than other SLS modeling approaches (Zhu et al., 2013).
To obtain the NCQ wave equation for heterogeneous media, we used the locality principle for pseudodifferential operators, which is valid for smoothly heterogeneous media.However, it may not hold for high-contrast heterogeneous media, of which the gradient in velocity and Q is large.Fortunately, from limited numerical tests, we found that numerical simulations are stable for reasonably large contrast heterogeneous velocity variations and Q variations.
We only consider the average γðrÞ to approximate the spatial varying γðrÞ in the fractional Laplacian terms when using the Fourier method.To the best of the authors' knowledge, numerical methods for spatially varying γðrÞ using the Fourier method are still not reported in the literature.Alternatively, finite-difference approximations for the spatial varying γðrÞ in the fractional Laplacian terms may be used but are still at an early stage of development (Lin et al., 2009).The finite-difference implementation is untested on practical problems.Future research is needed to investigate this topic.
Compared to the constant-Q wave equation by Carcione (2010), our NCQ wave equation with two fractional Laplacian operators can model amplitude loss and phase dispersion effects separately.Treeby et al. (2010) report that the inherent separation of loss and dispersion effects in their lossy wave equation is beneficial for attenuation compensation in medical imaging.Similarly, in our wave equation, the back-propagation in seismic inversion/migration can be done by reversing the sign of the loss operator in a way that A nearly constant-Q wave equation  (Zhu, 2014;Zhu et al., 2014).It appears that the present work has potential applications to seismic inversion/migration with attenuation compensation.

CONCLUSIONS
We have derived a NCQ viscoacoustic wave equation in the time domain using the fractional Laplacian.The advantage of using our fractional Laplacian formulation over the traditional fractional time derivative approach is avoidance of storing the time history of variables, and thus it is more economic in computational costs.In addition, our viscoacoustic wave equation has two fractional Laplacian operators to describe separated amplitude loss and phase dispersion effects, respectively.We provided a detailed formulation of this viscoacoustic wave equation for homogeneous media and smoothly heterogeneous media.We discussed the stability condition of the viscoacoustic wave equation.We also gave the first-order conservation equations for numerical simulation.
Numerical results demonstrate that the proposed viscoacoustic wave equation models the approximate constant-Q attenuation and dispersion.Modeling with the proposed wave equation is accurate in typical Q environments and stable in large contrast heterogeneous velocity and Q media.We believe that this fractional viscoacoustic wave equation is promising for use in accurate wavefield modeling for attenuating media.
Then, the equation is expressed in a matrix form where a ¼ c 2 0 ηk 2γþ2 and b ¼ c 2 0 τk 2γþ1 .We analyze the stability condition of the NCQ wave equation using the eigenvalue method, in which the eigenvalues of the matrix must be less than or equal to 1 in magnitude for stable simulation (Gazdag, 1981).Solving the eigenvalues λ of matrix A and assuming λ ≤ 1, the stability condition of NCQ wave equation is ðc 0 kÞ γþ1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðπγÞ p : (C-3) To guarantee that the solution is stable for all waves, c 0 ¼ c max is the maximum velocity, k ¼ k Ny ¼ π∕Δx is the Nyquist spatial wavenumber, and γ ¼ 1∕π tan −1 ð1∕Q min Þ.As γ ¼ 0, it becomes the stability condition of the acoustic wave equation (Gazdag, 1981): (C-4) To our best knowledge, the stability condition in equation C-3 is slightly stricter than equation C-4.In other words, attenuation (1∕Q min ) has a limited effect on the stability condition.

APPENDIX D CONSTANT-Q WAVE EQUATION IN HETEROGENEOUS MEDIA
In homogeneous media (c 0 and γ are constant with position), the NCQ wave equation (equation 10) is written as ∇ 2 σðr; tÞ: (D-1) If we write equation D-1 in the first-order mass-momentum conservation equations 15 and 16, we get the constitutive equation as follows: σðr; tÞ ¼ Lεðr; tÞ; (D-2) In smoothly heterogeneous media, we localize the pseudodifferential operator L at local point r 0 .All parameters related to velocity c 0 and γ are constant at fixed point r 0 .We have the same pseudodifferential operator L as equation D-3.It is reasonable to suppose that the L we want should be well approximated by Lðr 0 Þ when r is near r 0 ; thus, we unfreeze the localized point r 0 to general spatial r and define the pseudodifferential operator LðrÞ, LðrÞ ¼ M 0 ðrÞ ηðrÞð−∇ 2 Þ γðrÞ þ τðrÞ ∂ ∂t ð−∇ 2 Þ γðrÞ−1∕2 : (D-4) This generalized pseudodifferential operator from homogeneous to heterogeneous parameters is also called the freezing-unfreezing principle that appears in the literature in the context of hyperbolic and elliptic partial differential equations (Stein [1993] 12/14 to 171.65.238.130.Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Figure 1 .
Figure1.The relaxation function ψ 1−2γ ðtÞ varying with the previous time t in the constant-Q model.The curves denote values of Q ¼ 0.1 (solid), 1 (dash), and 10 (triangle), and 100 (circle).Note that when Q is small, the curve approximates a one-sided delta function.In contrast, when Q is large, the dependence becomes more flat (toward a step function).
Figure 4 shows four wavefields at 130 ms: (a) acoustic (lossless), (b) loss-dominated (equation 23), (c) dispersion-dominated (equation 22), and (d) viscoacoustic (lossy), respectively.Compared to the acoustic wavefield, Figure4bshows significantly reduced amplitude due to the loss while weak energy before the reference wavefront indicated by the red dots in Figure4is the result of the acausal attenuation term in equation 23.The dispersive acoustic wavefield in Figure4chas delayed phase but comparable amplitude; the viscoacoustic wavefield in Figure4dhas reduced amplitude and delayed phase.

Figure 3
Figure 3. (a) Attenuation and (b) phase velocity with frequency for Pierre Shale with Q ¼ 32.The numerical attenuation coefficients and phase velocity (circles) are computed using equations 20 and 21.The theoretical constant-Q values (solid lines) are computed using equations A-4 and A-5.The experimental results (diamond points) are from McDonal et al. (1958).

Figure 4 .
Figure 4. Four wavefield parts split by the dashed lines are generated by: acoustic (a), loss dominated (b), dispersion dominated (c), and viscoacoustic (d), respectively.The red dots indicate the reference acoustic wavefront (no attenuation).The simulation configuration is the same as that in the first example.The medium parameters are c 0 ¼ 2164 m∕s and Q ¼ 10.

Figure 8 .
Figure 8. Snapshots at 224 ms.A two-layer model with velocities of 1800 m∕s in the top layer and 3600 m∕s in the bottom layer, and a Q value of 30 in the top and 100 in the bottom is used.The interface is at a depth of 800 m.First row: homogeneous Q model with two-layer velocity: (a) computed by solving the NCQ wave equation, (b) computed by solving the viscoacoustic (SLS) wave equation, and (c) comparison of traces at horizontal 500 and 1000 m between (a) NCQ (circle line) and (b) SLS (solid line).Second row: homogeneous velocity model with two-layer Q: (d) computed by solving the NCQ wave equation, (e) computed by solving the SLS wave equation, and (f) comparison of traces at horizontal 500 and 1000 m between (d) NCQ (circle line) and (e) SLS (solid line).Third row: two-layer velocity and Q models: (g) computed by solving the NCQ wave equation, (h) computed by solving the SLS wave equation, and (j) comparison of traces at horizontal 500 m and 1000 m between (g) NCQ (circle line) and (h) SLS (solid line).The dotted line indicates the difference.

Figure 11 .Figure 12 .
Figure 11.Synthetic data calculated by (a) acoustic and (b) NCQ.Due to the different attenuation area, reflections at interface 3-4 are less attenuated and reflections at interface 4-5 experience significant loss.
To avoid calculating the spatial derivatives of the medium parameters in the wave equation 12, we solve instead the equivalent firstorder conservation equations with the additional of an external body force f: tÞ is the stress field; υðr; tÞ is the particle velocity vector; εðr; tÞ is the strain field at position r at time t; and cðrÞ and ρ 0 ðrÞ are the spatial acoustic velocity and density, respectively.Again, equation 17 reduces to the lossless constitutive stress-strain relation when γ → 0 (Q → ∞).
Downloaded 03/12/14 to 171.65.238.130.Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/ with laboratory measured data and Kjartanssan constant-Q theoretical data in a homogeneous medium.McDonal et al. (1958) use measured data of Pierre Shale in Colorado to find an approximate constant-Q behavior with the attenuation coefficient α ¼ 0.12f dB∕kft ≈ 0.0453f neper∕km, where f is the frequency in hertz.
Comparisons of the approximate dispersion in equation 9 (dashed line) with the exact constant-Q dispersion (see Appendix A) in equation 8 (solid line).We calculate these curves using different values of Q; i.e., Q ¼ 2, 10, and 100.The smaller the Q values, the larger the errors.T108 Zhu and Harris Downloaded 03/12/14 to 171.65.238.130.Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/SLS modeling algorithm requires more computer time in solving memory-variable equations at each time step.

T113
Downloaded 03/12/14 to 171.65.238.130.Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/naturally compensates for the amplitude loss.The back-propagation has been successfully demonstrated in applications of Q reversetime migration , p. 230).Plugging equation D-4 into equation D-2, and the massmomentum conservation equations 15 and 16, we have the approximate heterogeneous second-order wave equation: