Unsplit complex frequency-shifted PML implementation using Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling auxiliary differential equations for seismic wave modeling

The complex-frequency-shifted perfectly matched layer (cid:1) CFS-PML (cid:2) technique can efﬁciently absorb near-grazing incident waves. In seismic wave modeling, CFS-PML has been implemented by the ﬁrst-order-accuracy convolutional PML technique or second-order-accuracy recursive convolution PML technique. Both use different algorithms than the numerical scheme for the interior domain to update auxiliary memory variables in the PMLand thus cannot be used directly with higher-order time-marching schemes. We work with an unsplit-ﬁeld CFS-PML implementation using auxiliary differential equations (cid:1) ADEs (cid:2) to update the auxiliary memory variables. This ADE CFS-PML results in complete ﬁrst-order differential equations. Thus, the numerical scheme for the interior domain can be used to solve ADE CFS-PML equations. We have implemented ADE CFS-PML in the ﬁnite-difference time-domain method and in a nonstaggered-gridﬁnite-differencemethodwiththefourth-order Runge-Kuttascheme,demonstratingitsstraightforwardimple-mentationindifferentnumericaltime-marchingschemes.We havealsotheoreticallyanalyzedtheroleofthescalingfactorof CFS-PML;ittransformsthePMLtoatransverselyisotropicma-terial,reducingtheeffectivewavespeednormaltothePMLlayer andbendingthewavefronttowardthenormaldirectionofthe PMLlayer.Ournumericaltestsindicatethattheoptimalvaluere-ducesthepointsperdominantwavelengthattheoutermost boundarytothree,abouthalfthevaluerequiredbythenumerical scheme.WealsohavefoundthatthePMLequationsshouldbe derivedtakingthefree-surfaceboundaryconditionintoaccount inﬁnite-differencemethods.Otherwise,thefreesurfaceinthe PMLlayercausesinstabilityorineffectiveabsorptionofsurface waves.


INTRODUCTION
Special treatments are needed at the computational boundaries to absorb waves propagating outward when we simulate seismic wave propagation in an unbounded space.Techniques based on two concepts are commonly used: absorbing boundary conditions ͑ABCs͒, which impose a proper boundary condition at the outermost boundaries to satisfy the condition that waves only propagate outward ͑Clayton and Engquist, 1977;Liao et al., 1984;Bayliss et al., 1986;Higdon, 1986Higdon, , 1990;;Randall, 1988͒, and absorbing boundary layers ͑ABLs, or sponge layers͒, which use finite layers to gradually damp wave amplitude, so using the Dirichlet boundary condition at the outermost boundary does not generate a strong reflection ͑Cerjan et al., 1985;Sochacki et al., 1987͒.Strengths and weaknesses of each group can be found in Festa and Vilotte ͑2005͒ and Komatitsch and Martin ͑2007͒.Perfectly matched layers ͑PMLs͒ ͑Bérenger, 1994͒, an ABL technique, effectively absorb waves at a wide range of incident angles with only several to tens of layers and has become widely used in elastic wave modeling ͑Chew and Liu, 1996;Hastings et al., 1996;Collino and Tsogka, 2001;Marcinkovich and Olsen, 2003;Wang and Tang, 2003͒.The original ͑or standard͒ PML in Bérenger ͑1994͒ adopts a nonphysical splitting of the wavefield components and equations, leading to two different sets of equations for the interior domain and the PML layer; these are not trivial to incorporate into an existing numerical modeling code.Furthermore, the split-field PML is mathematically weakly well-posed ͑Abarbanel and Gottlieb, 1997͒.Alternative intepretations of the PML as an artificial anisotropic medium ͑Sacks et al., 1995;Gedney, 1996͒ or complex coordinate stretching ͑Chew andWeedon, 1994;Teixeira and Chew, 2000͒ have led to unsplit-field PML implementations involving convolution terms ͑Wang and Tang, 2003; Komatitsch andMartin, 2007͒, integral terms ͑Zeng andLiu, 2004;Drossaert and Giannopoulos, 2007a͒, or auxiliary differential equations ͑Ramadan, 2003͒.Hagstrom ͑2003͒ introduces a modified modal solution to derive PML equations; it is used by Appelö and Kreiss ͑2006͒ in 2D elastic wave modeling.This technique can provide a stable PML implementation by choosing the parameters satisfying the geometric relations in anisotropic media ͑Bécache et al., 2003͒.Recently, Meza-Fajardo and Papageorgiou ͑2008͒ show that stable results can be obtained in anisotropic media using a 3D damping profile in the PML based on the concept of complex coordinate stretching.
Though the standard PML has been widely used in seismic wave modeling, it can generate large, spurious reflections for near-grazing incident waves, low-frequency waves, or evanescent waves ͑Festa and Vilotte, 2005; Komatitsch and Martin, 2007;Drossaert and Giannopoulos, 2007b͒.The complex-frequency-shifted PML ͑CFS-PML͒ ͑Kuzuoglu and Mittra, 1996͒ is more efficient in such circumstances using a frequency-dependent damping ͑Festa et al., 2005;Festa and Vilotte, 2005;Drossaert andGiannopoulos, 2007a, 2007b;Komatitsch and Martin, 2007͒.The CFS-PML originally was implemented in a split-field form; three auxiliary memory variables are needed for each derivative ͑Gedney, 1998͒.Roden and Gedney ͑2000͒ derive an unsplit-field CFS convolutional-PML ͑C-PML͒ equation involving a convolution term that can be calculated efficiently using a recursive convolution algorithm ͑Luebbers and Hunsberger, 1992͒.This technique has been used in seismic wave modeling in elastic media ͑Komatitsch and Martin, 2007;Drossaert andGiannopoulos, 2007b͒ andporoelastic media ͑Martin et al., 2008a͒.Drossaert and Giannopoulos ͑2007a͒ derive an alternative CFS-PML implementation that involves integral terms, which can be calculated by a recursive integration PML ͑RIPML͒ algorithm.The recursive integration in C-PML is of first-order accuracy ͑Giannopoulos, 2008͒, and the trapezoidal rule in RIPML is of second-order accuracy ͑Drossaert and Giannopoulos, 2007a͒.Both approaches need an algorithm different from the numerical scheme in the interior domain to update the auxiliary memory variables.Thus, they cannot be used directly with higher-order timemarching schemes, such as a fourth-order Runge-Kutta scheme in the nonstaggered finite-difference velocity-stress method ͑Zhang and Chen, 2006͒ on a boundary-conforming grid to simulate seismic wave propagation in the presence of surface topography, in which the higher-order multistage Runge-Kutta scheme is important to allow a larger time step ͑Hixon, 1997͒.
In the field of electromagnetic simulation, Ramadan ͑2003͒ proposes an unsplit-field implementation of the standard PML using auxiliary differential equations ͑ADEs͒.Wang and Liang ͑2006͒ extend this approach to CFS-PML with a 2D alternating-direction-implicit ͑ADI͒ finite-difference time-domain ͑FDTD͒ method.This ADE CFS-PML implementation results in complete first-order differential equations for the wavefield components and the auxiliary memory variables.The same numerical scheme is also used in the interior domain to solve the ADE CFS-PML equations, regardless of the time-stepping scheme.
In this paper, we implement the ADE CFS-PML in seismic wave simulation in the FDTD method and in the nonstaggered-grid finitedifference method ͑Zhang and Chen, 2006͒ using a fourth-order Runge-Kutta time-marching scheme.We derive the ADE CFS-PML equations for the velocity-stress equations.Then we analyze the free-surface boundary condition inside the PML layer and modify the PML equations based on the free-surface condition to avoid the instability caused by incompatibility between the PML equations and the free-surface boundary condition.We analyze the role of the scaling factor ␤ in the CFS stretching function, which is still unclear in the literature of seismic wave modeling.Next, we perform a series of numerical tests in a thin-slab model to validate the ADE CFS-PML implementation and to investigate the optimal values of the parameters in the CFS stretching function.We add the free surface to the thin-slab model to demonstrate the efficiency of the modified PML equations, taking into account the free-surface condition in absorbing surface waves.Finally, we simulate seismic wave propagation in a two-layer medium with surface topography in a narrow vertical slice to demonstrate the implementation of the ADE CFS-PML in the fourth-order Runge-Kutta scheme and its potential effectiveness in seismic exploration.

VELOCITY-STRESS EQUATION AND FINITE-DIFFERENCE NUMERICAL SCHEME
We consider the following velocity-stress formulation of the elastic wave equations for an isotropic medium in the Cartesian coordinate ͑x,y,z͒: where v ‫ס‬ ͑v x ,v y ,v z ͒ T is the particle velocity vector; T is the transpose operator; is the stress tensor; the components are ij ,i,j ͑x,y,z͒; is the mass density; and c is the stiffness tensor.In this paper, a comma followed by a subscript t, x, y, or z means a derivative with respect to t, x, y, or z.Sometimes it is more convenient to discuss the PML implementation in the component form.The expression of the v x component of equation 1 is We use two numerical schemes to solve the velocity-stress equations 1 and 2. One is the FDTD method, with a second-order-accuracy leapfrog time-marching scheme on a staggered grid ͑Levander, 1988; Graves, 1996͒.The other is the nonstaggered-grid finite-difference method with a fourth-order Runge-Kutta time-marching scheme on a boundary-conforming grid ͑Zhang and Chen, 2006͒.Because PML implementations do not rely on how spatial derivatives are calculated, we only list the time-marching schemes here.The second-order leapfrog scheme can be written as where ⌬t is the time step; superscripts n ‫מ‬ 1 / 2, n, n ‫ם‬ 1 / 2 and n ‫ם‬ 1 mean time level; L͑ n‫2/1ם‬ ͒ is the right-hand side of equation 2 evaluated at time level n ‫ם‬ 1 / 2 using an arbitrary-order-accuracy spatial staggered finite-difference operator; and L͑v n ͒ is the righthand side of equation 1 at time level n.
To describe the implementation of the fourth-order Runge-Kutta scheme, we write equations 1 and 2 in vector form: w ,t ‫ס‬ Aw ,x ‫ם‬ Bw ,y ‫ם‬ Cw ,z , ͑6͒ where w ‫ס‬ ͑v x ,v y ,v z , xx , yy , zz , yz , xz , xy ͒ is the velocity-stress vector and where A, B, and C are the coefficient matrices that can be easily derived.The fourth-order Runge-Kutta scheme can be expressed as where ␣ i and ␤ i ͑1 Յ i Յ 4͒ are the coefficients of the fourth-order Runge-Kutta scheme and L͑w n ͒ means the right-hand side of equation 6 calculated for variable w n using the spatial operator.Here, we use the DRP/opt MacCormack operator ͑Hixon, 1997͒, which utilizes dispersion-relation-preserving methodology ͑Tam and Webb, 1993͒ to optimize the dispersion and dissipation errors.

CFS-PML USING ADES
For simplicity of notation, we only discuss the PML at the positive x-axis face and also presume the PML starts from x ‫ס‬ 0.
Based on the concept of complex coordinate stretching ͑Chew and Weedon, 1994͒, the equations inside the PML layer have exactly the same form as in the physical domain ͑equations 1 and 2͒ except that the coordinate x is replaced with a complex stretched coordinate x ˜: x where s x is the complex stretching function that determines the characteristics of the PML model and is the variable of integration. From and the CFS-PML, in which In equations 11 and 12, d x Ն 0 is the attenuation factor that causes the amplitude of the wavefield to be reduced exponentially inside the PML layer, ␣ x Ն 0 is the frequency-shifted factor that makes the attenuation frequency dependent, and ␤ x Ն 1 is the scaling factor.The latter has been found in numerical tests to be important for absorption of evanescent waves ͑Liu, 1999; Drossaert and Giannopoulos, 2007b͒ and near-grazing incident waves ͑Drossaert and Giannopoulos, 2007a͒, but the reason has not been theoretically explained in the literature of seismic wave modeling.We discuss the role of ␤ x in later sections.The parameters d x ͑x͒, ␣ x ͑x͒, ␤ x ͑x͒, and therefore s x ͑x͒ are all functions of x.In the following, we omit the dependence on x for brevity.
The basic idea of ADE implementation of CFS-PML is to separate the derivative with respect to the complex coordinate into two parts and only keep i in the denominator in one of them: where the auxiliary memory variable Taking equation 14 into equation 10, we obtain the unsplit-field CFS-PML equation for the v x component in the frequency domain: Transforming equations 15 and 16 to the time domain, we get the final ADE CFS-PML equation of v x as It is worth noting that we write the time-domain PML equation as a correction to the original wave equation because it is easy to achieve parallel balance and to implement in an existing code ͑Drossaert and Giannopoulos, 2007b͒.The ADE CFS-PML equations for the y-and z-axes and other components can be obtained in the same manner.The complete ADE CFS-PML equations for the velocitystress equations are in Appendix A.
The velocity-stress equations 1 and 2 and the ADE CFS-PML equations 17 and 18 are all first-order partial-differential equations with respect to time t, so they can be updated by the same time-stepping scheme.In the staggered second-order leapfrog scheme, equation 18 can be discretized as

͑19͒
Then we can update T xx x and v x through The implementation in the fourth-order Runge-Kutta scheme is also straightforward.One can just include the memory variables in the solution vector w in equation 6, then use the Runge-Kutta scheme ͑equation 7͒ to update the PML equations, including the memoryvariable equations.
The computer memory requirement of ADE CFS-PML is the same as for C-PML ͑Komatitsch and Martin, 2007͒ or RIPML ͑Drossaert and Giannopoulos, 2007a͒ -one memory variable per normal derivative to the PML layer.

ADE CFS-PML BOUNDARY AND FREE-SURFACE BOUNDARY CONDITIONS
To our knowledge, there is no publication on how to implement the free-surface boundary condition in the PML region, which is crucial to absorb surface waves efficiently in the FDTD method or to avoid the stability problem in the nonstaggered finite-difference method.
At the flat surface in the physical region, the free-surface boundary condition requires zz ‫ס‬ 0, yz ‫ס‬ 0, xz ‫ס‬ 0. ͑22͒ In the PML layer, the zero-valued stress components should remain zero at the free surface because the effect of the PML is to exponentially damp the wavefield to zero.
Because the velocity derivatives are related to the stress components through the stress-strain relation ͑equation 2͒, the zero-stress condition leads to velocity derivative constraints at the free surface: Taking this equation into the stress-strain relation, the true equations to update xx and yy at the free surface are To absorb the surface wave effectively, the ADE CFS-PML equations of xx and yy at the intersection of the free surface and the PML should be derived from equations 24 and 25 as equation A-13 in Appendix A.
In the FDTD method, there are two possible locations for the free surface: coincident with the normal stress components or with the shear stress component xz ͑Gottschammer and Olsen, 2001; Kristek et al., 2002͒.As the numerical tests presented in the following section show, the modified ADE CFS-PML equation A-13 is crucial to absorb surface waves effectively if the free surface coincides with the normal stress components.If the free surface coincides with xz , the free-surface boundary condition is implemented as in the physical domain by setting xz and yz to zero at the free surface, and no special treatment is needed in the PML layer.In the nonstaggered finite-difference method, our numerical tests show that the modified PML equation A-13 is essential for the stability of the solution if the velocity-component free-surface boundary condition ͑equation 23͒ is applied.

ROLE OF THE SCALING FACTOR ␤
rection into the layer by increasing the corresponding component of the wave vector while leaving the tangential component unchanged; at the same time, it increases the damping in the normal direction into the layer" with a stretching function s x ‫ס‬ ␤ x ͑1 ‫ם‬ ͑d x / ͑␣ x ‫ם‬ i͒͒͒.In this section, we present a theoretical explanation of how ␤ affects the absorption of near-grazing incident waves.
We rewrite the CFS stretching function ͑equation 12͒ as The CFS-PML based on the velocity-stress equations 1 and 2 and the stretching-function equation 26 is equivalent to using a stretching function to stretch the following equations: The solution of equations 28 in a homogeneous isotropic medium is the same as the solution of the velocity-stress equations 1 and 2 in a transversely isotropic medium ͑x is the symmetry axis͒ that has the stiffness tensor

͑29͒
We can see that ␤ x introduces anisotropy and reduces the phase velocity normal to the PML layer to C / ␤ x .The plane-wave solution of equation 28 can be written as u ‫ס‬ u 0 e i͑t‫מ‬s•x͒ , ͑30͒ where u 0 is the polarization vector, s ‫ס‬ s C ˜is the slowness vector normal to the wavefront surface, and C ˜is the phase velocity depending on the propagation direction.For simplicity, we only consider waves propagating in the ͑x,z͒-plane.Equation 30 becomes where is the angle between the normal direction of the wavefront and x-axis.Taking equation 27 into equation 8, we can obtain the stretched x ˜: x By substituting equation 32 into equation 31, we get the plane-wave solution of the CFS-PML:

͑33͒
If the wave propagates normal to the PML layer ͑phase velocity C ‫ס‬ C / ␤ x and cos ‫ס‬ 1, assuming a constant ␤ x ͒, the damping term in which ␤ x is canceled out.This indicates that for normal incident waves, ␤ x only decreases the wave-propagation speed without increasing damping.We verify this by calculating the waveform in the PML with different ␤ x values for a normal incident plane wave using the FDTD method ͑␣ x ‫ס‬ 0 and polynomial scaled ␤ x and d x ͒.
Figure 1a shows the maximum amplitude for different ␤ values ͑␤ 0 ‫ס‬ 1 and 9͒ at each grid point inside a 25-cell-thick PML layer.The maximum amplitude distributions do not change with ␤ 0 , verifying that ␤ does not directly increase damping.Figure 1b compares the time series of the particle velocity for ␤ 0 ‫ס‬ 1 and ␤ 0 ‫ס‬ 9 at the point 15 cell into the PML.We can see that the wave using ␤ 0 ‫ס‬ 9 arrives later than the wave using ␤ 0 ‫ס‬ 1 because of the reduced wave speed by ␤ x .The stretching function in Petropoulos ͑2000͒ is of a form s x ‫ס‬ ␤ x ͑1 ‫ם‬ ͑d x / ͑␣ x ‫ם‬ i͒͒͒, whose damping term contains the product of ␤ x and d x in the exponent.Therefore, ␤ x in Petropoulos ͑2000͒ affects damping.
For grazing incident waves, cos is very small if ␤ x is not used.Consequently, D in equation 33 may not be small enough to efficiently damp the wave propagating subparallel to the PML layer, and strong spurious reflections can occur.If ␤ x is used, the normal direction of the wavefront is bent toward the normal direction of the PML layer as a result of the anisotropy caused by ␤ ; therefore, cos in the damping term is increased and the absorption is improved.We will demonstrate the effect of ␤ x using particle-velocity snapshots in the first numerical test in the next section.Note that ␤ x cannot be canceled out in the damping term for oblique incident waves because the phase velocity is a function depending on several components of the stiffness tensor.

NUMERICAL TESTS
To demonstrate the ADE CFS-PML, we present results of three numerical tests.One is performed in a thin, homogeneous 3D slab surrounded by PMLs on all surfaces to investigate the optimal values of the CFS parameters.In the second test, we add a free surface to the thin-slab model to validate the modified PML equations that take into account the free-surface boundary condition in the absorption of surface waves.We use the FDTD method in the first two tests.In the third test, we use the nonstaggered finite-difference method with the fourth-order Runge-Kutta scheme in a two-layer model with a 3D topography to validate the ADE CFS-PML implementation in a higher-order time-marching scheme.
Before showing the numerical test results, let us briefly review the choice of the parameter values in the CFS function.According to Roden and Gedney ͑2000͒, d, ␤ , and ␣ should be spatially scaled in the PML layer to reduce reflection errors in the discrete space.The value of d usually is zero at the PML-interior interface and maximum at the exterior boundary.The value of ␤ ‫ס‬ 1 at the PML-interior interface and maximum at the exterior boundary.The parameter ␣ is scaled in a reverse manner so the value is maximum at the PML-interior interface and gradually reduces to zero at the exterior boundary to absorb low-frequency waves.
In this paper, we choose the commonly used p-order polynomial scaling function in the PML layer ͑x Ն 0͒ ͑Gedney, 1998; Roden and Gedney, 2000͒:

͑36͒
where d 0 and ␤ 0 are the maximum values of d and ␤ at the exterior boundary ͑the optimal ␤ 0 is discussed in the next section͒, ␣ 0 is the maximum value of ␣ at the PML-interior interface, L is the width of the PML layer, and x is the distance to the PML-interior interface.
where c P is the compressional wave speed and R is the theoretical reflection coefficient for a normal-incident plane P-wave with a Dirichlet condition ͑v ‫ס‬ 0 and ‫ס‬ 0͒ at the exterior boundary of the PML layer.The trade-off between a small and a large d 0 is that a toosmall d 0 cannot damp the wavefield enough in the PML layer and strong spurious reflections from the Dirichlet boundary will propagate back to the physical domain.A too-large d 0 leads to spurious reflections caused by the dispersion ͑Collino and Tsogka, 2001͒.Thus, the optimal d 0 depends on L and R.
Collino and Tsogka ͑2001͒ propose R for different widths of the PML layer as R 5 ‫ס‬ 0.01, R 10 ‫ס‬ 0.001, and R 20 ‫ס‬ 0.0001, where the superscript denotes the width of the PML region in the grid spacing.Marcinkovich and Olsen ͑2003͒ approximate Collino and Tsogka's ͑2001͒ relationship between R and PML width by a third-order polynomial function to obtain the R value for an arbitrary PML width between 1 and 20 layers.We note that the relationship between R and PML width follows a logarithmic equation: Thus, R for an arbitrary PML layer width of N cell size ͑could be Ͼ20͒ can be expressed as log 10 ͑R͒ ‫מס‬ log 10 ͑N͒ ‫מ‬ 1 log 10 ͑2͒ ‫מ‬ 3. ͑39͒ Equation 37 is optimized for normal incident waves.For oblique incident waves, a larger d 0 than the value from equation 37 is needed to obtain optimal damping.From numerical tests, we find that increasing d 0 to two to three times the value from equation 37 yields better solutions for grazing incidence.Figure 1.͑a͒ Maximum amplitude distribution for different ␤ values ͑␤ 0 ‫ס‬ 1, solid line; ␤ 0 ‫ס‬ 9, dashed line͒ for a 25-cell-thick PML of a normal incident plane wave ͑␣ ‫ס‬ 0 and an optimal d are used͒.͑b͒ Velocity seismogram at the grid point 15 cells into the PML for ␤ 0 ‫ס‬ 1 ͑solid line͒ and ␤ 0 ‫ס‬ 9 ͑dotted line͒.The wave using ␤ 0 ‫ס‬ 9 arrives later than using ␤ 0 ‫ס‬ 1 because of the reduced wave speed caused by ␤ x .
Festa and Vilotte ͑2005͒ recommend ␣ 0 to be f c , where f c is the dominant frequency of the source-time function.In our numerical tests, different values of ␣ 0 ͑1.0 f c , 1.5 f c , and 2.0 f c ͒ did not make a significant difference.Thus, ␣ 0 ‫ס‬ f c seems a good choice in practice.

Thin, homogeneous 3D slab surrounded by PMLs
This test compares the performance difference between the ADE CFS-PML and the standard PML ͑by setting ␤ 0 ‫ס‬ 1 and ␣ 0 ‫ס‬ 0͒ for near-grazing incident body waves and demonstrates the effects of the parameters d 0 , ␤ 0 , and ␣ 0 .
Figure 2a shows the computational configurations.The medium has a compressional-wave speed c P ‫ס‬ 3200 m / s, shear-wave speed c S ‫ס‬ 1600 m / s, and density ‫ס‬ 2200 kg/ m 3 .The grid spacing is 10 m, the total grid number is 600ϫ 600ϫ 81, and the size of the model is 5990ϫ 5990ϫ 800 m.The vertical-force source is applied at ͑200, 300, 600͒ m, 10 cells away from the PML-interior interface at the left ͑negative x-axis͒ face and the top ͑positive z-axis͒ face.The receiver is located at ͑5790, 3000, 600͒ m, 10 cells away from the PML-interior interface at the left and top faces.͑Additional terms used are right ͓positive x-axis͔ face, front ͓negative y-axis͔ face, back ͓positive y-axis͔ face, and bottom ͓negative z-axis͔ face.͒ The source-time function is a Ricker wavelet with a center frequency of 5 Hz and a time delay of 0.6 s.The time step is 1.5 ms, and the total time step is 5000.The result is a time window of 7.5 s, in which the reflected S-wave can travel through the computational domain.The value of d 0 is obtained through equation 37. Different values of ␤ 0 and ␣ 0 ͑0,0.5 f c , 1.0 f c , 1.5 f c , and 2.0 f c ͒ are tested.
Figure 2b-f shows snapshots of the v z component at 3.75 s.The Pwave has passed, and the dominant phase is the S-phase.The result of a large model ͑two times larger along the x-and y-axes and 12 times larger along the z-axis with a 20-cell-thick PML͒ is used as the reference ͑Figure 2b͒.The standard PML shows spurious evanescent waves at the top face ͑Figure 2c͒.The parameter ␣ can eliminate such spurious evanescent waves, but there is a spurious reflected body wave immediately after the S-phase ͑Figure 2d͒.Increasing ␣ 0 to 2 f c does not significantly affect the result ͑not shown͒.Using both ␤ and ␣ reduces the spurious reflected body wave and achieves better results ͑Figure 2e͒.We also increase d 0 to three times larger, which can dramatically improve the absorption of grazing incident waves ͑Figure 2f͒, but the spurious reflection for normal incidence strengthens ͑at the back face, Figure 2f͒.
These results indicate that it is important to choose appropriate values for ␣ 0 , ␤ 0 , and d 0 for the absorption of grazing incident waves.We need a quantitative measurement of the results to determine which values of the parameters are optimal.For this purpose, we use the following integral as a local error at position x: where v PML ͑x,t͒ is the simulated wavefield in the PML model and v ref ͑x,t͒ is the counterpart in the reference model.
Figure 3 shows the 3D distribution of the local error and particlevelocity time series at the receiver in Figure 2a for different parameter values.These plots further confirm the different effects of ␣, ␤ , and d on absorption of the grazing incident wave, as shown in Figure 2b-f.There are clear, spurious reflected S-waves ͑larger amplitude of the S-wave in Figure 3a͒ and spurious evanescent waves ͑Figure 3b͒ in the standard PML, caused by near-grazing incident waves.The spurious reflected S-wave propagates in nearly the same direction as the near-grazing incident S-wave.This it is nearly invisible in the snapshot in Figure 2d but can be observed in the time series in Figure 3d for the CFS-PML using ␣ 0 ‫ס‬ f c .
Figure 3d suggests that although ␣ can eliminate the spurious evanescent waves, it does little in the absorption of near-grazing inci- ADE CFS-PML for seismic wave modeling T147 dent body waves.So we observe similar local-error distributions in Figure 3c and f.From the damping-term expression ͑without ␤ x ͒, we see that ␣ does not significantly affect the incident angle , and cos is very small for near-grazing incident waves.We already know that ␤ x bends the wavefront toward the normal direction of the PML, decreasing the slowness angle and increasing the damping of waves at near-grazing incidence.The results in Figure 3g-i for ␤ 0 ‫ס‬ 10 are as expected: The S-wave amplitude is close to the reference solution ͑Figure 3g͒, and the local error in Figure 3i is much smaller than the local error in Figure 3c and f.We can clearly observe the effect of ␤ x for grazing incident waves in this numerical test from comparing the snapshots in Figure 4a and  b.As shown in Figure 3j-l, a three-times-larger d 0 ͑with ␤ 0 ‫ס‬ 7, ␣ 0 ‫ס‬ f c ͒ has nearly perfect fits for P-and S-waves ͑Figure 3j-l͒.
To investigate the optimal value of ␤ , we use an averaged local error in a subgrid located 10 cells away from the PML-interior interfaces and in the xz-plane through y ‫ס‬ 3000 m ͑20Յ i Յ 580, j ‫ס‬ 301, 20Յ k Յ 61͒ as the global error.Figure 5a shows the contour plot of the global error as a function of ␤ 0 and ␣ 0 .It can be seen that ␣ 0 / ͑ f c ͒ between one and two does not affect the overall accuracy.The optimal value of ␤ 0 is about 10-12 in this case.The global error with respect to ␤ 0 in Figure 5b shows that the optimal ␤ 0 is insensitive to the PML thickness N but increases when shear-wave speed c S increases ͑equivalent to grid-spacing decreases or sourcefrequency decreases͒.Theoretically, a larger ␤ can bend waves more efficiently, but a very large ␤ can make the wavelength too short to be resolved by the numerical scheme, which can cause strong numerical spurious reflections.
Figure 5c and d shows the normalized points-per-dominantwavelength distribution inside the PML layer at the center frequency for the optimal ␤ 0 in Figure 5b   dominant wavelength at the exterior boundary is approximately 0.5 in difference cases.Thus, we can estimate the optimal ␤ 0 through Note that we use the same parameter values in all PMLs on the six faces in this numerical test.In practice, one can use larger ␤ 0 and d 0 only for a particular PML layer for which grazing incidence is important.

Thin 3D slab with the free surface
We change one face of the thin slab in the previous test to the free surface to demonstrate the absorption of surface waves of the ADE CFS-PML and the effect of the modified PML equations that take into account the free-surface boundary condition.
Figure 6a shows the computational configurations.The star symbol shows the location of the vertical force ͑200, 3000, 795͒ m, 10 cells away from the PML-interior interface in the negative x-direction and one-half cell below the free surface.All other parameters are the same as the model in Figure 2a.The free-surface boundary condition is implemented using the adjusted finite-difference approximations ͑AFDA͒ technique ͑Kristek et al., 2002͒, in which a compact finite-difference operator and biased finite-differ- The mimimal global error is around 11 for different PML thicknesses of 10, 12, and 20 ͑c S ‫ס‬ 1600 m / s͒, indicating the optimal ␤ 0 is insensitive to the thickness of the PML.For N ‫ס‬ 20 and c S ‫ס‬ 2200 m / s, the optimal ␤ 0 is near 15.͑c͒ Normalized points per dominant wavelength variation inside the PML layer for S-waves for different PML thickness N and S-wave speed c S .The normalized points per dominant wavelength is defined as PPW/ PPW 0 , where PPW is points per dominant wavelength and PPW 0 is the minimal PPW requirement of the numerical scheme.For the fourth-order staggered finite-difference method, PPW 0 ‫ס‬ 6. ͑d͒ Zoomed-in view of Figure 5c to highlight details.The normalized points per dominant wavelength is approximately 0.5 at the exterior boundary in difference cases.ence operators are used to calculate derivatives with respect to the z-axis at grid points near the free surface.Figure 6 is the snapshot of the v z component at 3.75 s.The free surface coincides with the location of the normal stress components.
Because waves do not impinge the PML at near-grazing incidence, the standard PML ͑Figure 6c͒ and the ADE CFS-PML ͑Figure 6d, ␤ 0 ‫ס‬ 1 and ␣ 0 ‫ס‬ f c ͒ both efficiently absorb surface waves.For the same reason, increasing ␤ 0 ͑not shown here͒ does not make a significant difference.If the modified ADE CFS-PML equation A-13, which incorporates the free-surface boundary, is not used, a relatively strong reflected-surface phase can be clearly seen ͑Figure 6e͒.If we set the location of xz as the free surface in the FDTD method, the amplitude of the reflected wave ͑Figure 6f͒ is smaller than Figure 6e but stronger than Figure 6c.Thus, the location of the free surface at normal stress position has better absorption of surface waves than at the xz position.

Two-layer medium with free-surface topography in a narrow vertical slice
In this test, we validate the ADE CFS-PML implementation in the fourth-order Runge-Kutta scheme and demonstrate the optimal performance of theADE CFS-PML in a narrow vertical slice to simulate wave propagation in a common situation involving seismic profiles.
The medium has two layers ͑Figure 7͒ separated at z ‫8מס‬ km.A boundary-conforming grid is used to conform the grid with the surface topography.The horizontal grid spacing is 50 m; a varying grid spacing ͑ϳ33-60 m͒ is used vertically.The source-time function is a Ricker wavelet with a center frequency of 1.5 Hz and a time delay of 1 s.The value of d 0 is calculated from equation 37, where the c P of the first layer ͑3000 m / s͒ is used for the PML layers at the x-and y-axes faces and where c P of the lower layer ͑5000 m / s͒ is used for the bottom PML layer.The value of ␣ 0 is f c .The optimal ␤ 0 is estimated to be eight from equation 41, using the surface-wave speed of the first layer.We also simulate with ␤ 0 ‫ס‬ 1 for comparison with the optimal ␤ 0 value.Numerical solution of a larger model ͑double the length and depth, five times the width, a 25-cell-thick PML͒ is used as reference.
The difference between the ADE CFS-PML ͑Figure 8b, ␣ 0 ‫ס‬ f c and ␤ 0 ‫ס‬ 1͒ solution and the reference ͑Figure 8a͒ is so small that it is hard to observe from the snapshots when the color is scaled to 1% of the maximum amplitude.We cannot tell the difference between the v z time series of the reference ͑solid line͒, PML with ␤ 0 ‫ס‬ 1 ͑dashed-dotted line͒, and PML with ␤ 0 ‫ס‬ 8 ͑dotted line͒ in Figure 9a for the receiver in Figure 7.In Figure 9b, we expand Figure 9a by limiting the time axis to approximately 13-21 s and reducing the amplitude axis to 0.2%.Even at this scale, the solution of ␤ 0 ‫ס‬ 8 ͑dotted line͒ still overlaps with the reference.Figure 9c  It is clear that the optimal ␤ 0 yields a better result than ␤ 0 ‫ס‬ 1.The amplitude of the spurious waves is below 0.1% of the maximum amplitude.It can be seen that the ADE CFS-PML can efficiently absorb waves in this narrow, two-layer topographic model.The spurious waves do not significantly affect the reflections from the inner interface.This numerical test indicates that we can use a narrow-slice mesh with the ADE CFS-PML to simulate full-wave propagation in complex structured models in a common situation involving seismic profiles.

CONCLUSION
We have introduced an efficient unsplit-field CFS-PML implementation using ADEs in seismic wave modeling.We derive the complete ADE CFS-PML equations for the velocity-stress equations.Because the ADE CFS-PML equations are first-order partialdifferential equations, they can be solved by the same numerical scheme used in the inner domain no matter what time-marching scheme is used.We have demonstrated this by implementing the ADE CFS-PML in the second-order leapfrog scheme and the fourthorder Runge-Kutta scheme.
Also, we have analyzed the role of the scaling factor ␤ in CFS-PML, demonstrating that it introduces an intermediate transformation of the velocity-stress equations into an anisotropic medium that has a lower wave speed normal to the PML layer.For normal incident waves, it only delays the waves propagating into the PML layer without increased damping.For oblique incident waves, it can bend the waves toward the normal direction, increasing cos in the damping term and improving absorption.Numerical tests indicate that the frequency-shifted factor ␣ only eliminates the evanescent waves caused by near-grazing incident waves but does not significantly improve the absorption of the reflected body waves.The scaling factor ␤ is more efficient and important in absorbing grazing incident waves.
Our numerical tests indicate that the incompatibility between PML equations and the free-surface boundary-condition implementation can cause an unstable wave originating from the free surface in the PML layer in the nonstaggered-grid finite-difference method.Thus, we derive the PML equations, taking the free-surface boundary condition into account.Numerical tests validate that these modified PML equations solve the instability problem in the nonstaggered finite-difference method and improve the efficiency of absorbing surface waves in the staggered finite-difference method when the free surface coincides with the normal stress components.
Three numerical tests validate the ADE CFS-PML implementation in the FDTD method and the nonstaggered-grid finite-difference method using the fourth-order Runge-Kutta scheme.We have investigated the different values of the CFS parameters on the efficiency of absorption based on quantitative error criteria.Our numerical tests indicate that the optimal ␤ reduces the points per dominant wavelength at the exterior boundary to three for the fourth-order finite-difference operator, about half the value required by the numerical scheme.Different values of ␣ 0 ͑1.0 f c , 1.5 f c , and 2.0 f c ͒ do not make a significant difference.Thus, ␣ 0 ‫ס‬ f c seems a good choice in practice.For grazing incidence, adjusting only ␣ and ␤ is not enough; one also needs to increase the maximum value d 0 of the damping profile d for optimal absorption.
A numerical test of a two-layer model with a surface topography in a narrow vertical slice indicates that ADE CFS-PML yields accurate results compared with the reference solution in a much larger and thicker model, making it possible to carry out efficient 3D forward simulation in a common situation involving seismic profiles in seismic exploration.

ACKNOWLEDGMENTS
This research was supported by the U. S. National Science Foundation ͑grant NSF-0551117͒ and the U. S. Air Force Research Laboratory ͑contract fa8718-06-c-0014͒.We thank reviewer Antonios Giannopoulos for his reminder of the first-order accuracy of the C-PML implementation and comments on ADE CFS-PML implementation in the FDTD method.Constructive comments by editors Jeff Shragge and Evert Slob, reviewer Valeriy Philip Chkuaseli, and an anonymous reviewer improved the manuscript.

APPENDIX A ADE CFS-PML EQUATIONS FOR ELASTIC WAVE EQUATIONS
ADE CFS-PML equations for the velocity-stress formulation of the elastic wave equations are The amplitude is zoomed in to 0.2% of that in Figure 9a.The difference can be seen between the reference and ␤ 0 ‫ס‬ 1.The solution with ␤ 0 ‫ס‬ 8 overlaps with the reference.͑c͒ Difference between solutions from the CFS ADE-PML and the reference.The solid line shows the difference for ␤ 0 ‫ס‬ 1, and the dotted line shows the difference for the optimal ␤ 0 ‫ס‬ 8.The maximum y-axis is about 0.5% of Figure 9a.

Figure 3 .Figure 4 .
Figure 3.The top row shows v z at the receiver in Figure 2a for different PML parameters.The amplitude axis of the middle row is scaled to 1% of that of the left column to highlight the difference.The bottom row shows the local-error distribution.The first arrival in the middle row is the P-wave, which is nearly invisible in the left row.͑a-c͒ Standard PML; ͑d-f͒ ␤ 0 ‫ס‬ 1 and ␣ 0 ‫ס‬ f c ; ͑g-i͒ ␤ 0 ‫ס‬ 10 and ␣ 0 ‫ס‬ f c ; ͑j-l͒ ␤ 0 ‫ס‬ 7, ␣ 0 ‫ס‬ f c and a three-times-larger d 0 .

Figure 5 .
Figure 5. ͑a͒ Contours of the global error as a function of ␣ 0 and ␤ 0 .͑b͒ Global-error variation with respect to ␤ 0 , PML thickness N, and S-wave velocities c S ͑␣ 0 is fixed to f c ͒.The mimimal global error is around 11 for different PML thicknesses of 10, 12, and 20 ͑c S ‫ס‬ 1600 m / s͒, indicating the optimal ␤ 0 is insensitive to the thickness of the PML.For N ‫ס‬ 20 and c S ‫ס‬ 2200 m / s, the optimal ␤ 0 is near 15.͑c͒ Normalized points per dominant wavelength variation inside the PML layer for S-waves for different PML thickness N and S-wave speed c S .The normalized points per dominant wavelength is defined as PPW/ PPW 0 , where PPW is points per dominant wavelength and PPW 0 is the minimal PPW requirement of the numerical scheme.For the fourth-order staggered finite-difference method, PPW 0 ‫ס‬ 6. ͑d͒ Zoomed-in view of Figure5cto highlight details.The normalized points per dominant wavelength is approximately 0.5 at the exterior boundary in difference cases.

Figure 6 .
Figure 6.Snapshots of v z of the free-surface slab model at 3:75 s.The color is scaled to 0.1% of the maximum absolute amplitude.White lines show the PML-interior interfaces.͑a͒ Thin slab model with the free surface.A 10-cell-thick PML ͑dashed lines͒ is applied on all the faces except the free surface.͑b͒ Reference solution; ͑c͒ standard PML; ͑d͒ ␤ 0 ‫ס‬ 1 and ␣ 0 ‫ס‬ f c ; ͑e͒ standard PML without using the modified PML equation A-13.Strong spurious reflected surface waves can be clearly seen.͑f͒ Standard PML with the free surface coincident with the location of the xz component.
plots the v z difference between the PML with ␤ 0 ‫ס‬ 1 and reference ͑v z ␤ 0 ‫1ס‬ ‫מ‬ v z ref ͒ and between the PML with ␤ 0 ‫ס‬ 8 and reference ͑v z ␤ 0 ‫8ס‬ ‫מ‬ v z ref ͒.

Figure 8 .
Figure 7.The two-layer model with the free-surface topography, defined by z ‫ס‬ 1000 exp͑‫מ‬r 2 / 1000 2 ͒, where r ‫ס‬ ͱ x 2 ‫ם‬ y 2 .When r Ͼ 2000 m, z 0 is damped to zero by a Gaussian function.The model is 19.95 km long, 6.15 km wide, and 10.7 km deep ͑including PML layers͒.The interface is at z ‫8מס‬ km, and the medium parameters are c P ‫ס‬ 3000 m / s, c S ‫ס‬ 2000 m / s, and ‫ס‬ 1200 kg/ m 3 for the upper layer and c P ‫ס‬ 5000 m / s, c S ‫ס‬ 2500 m / s, and ‫ס‬ 1800 kg/ m 3 for the lower layer.A12-cell-thick PML is applied on all faces except the free surface ͑white lines͒.The star on the left side shows the location of the vertical force ‫,0008מ͑‬ 0, ‫מ‬ 500͒ m, and the triangle ͑right side͒ indicates the receiver's location ͑9000, 0, 0͒ m, seven cells'distance to the PML-interior interface at the right face.

Figure 9 .
Figure9.͑a͒ Comparison of seismograms v z for the reference ͑solid line͒, ␤ 0 ‫ס‬ 1 ͑dashed line͒, and ␤ 0 ‫ס‬ 8 ͑dotted line͒ at the receiver shown in Figure7.͑b͒ Detailed comparison within the time window of approximately 13-21 s.The amplitude is zoomed in to 0.2% of that in Figure9a.The difference can be seen between the reference and ␤ 0 ‫ס‬ 1.The solution with ␤ 0 ‫ס‬ 8 overlaps with the reference.͑c͒ Difference between solutions from the CFS ADE-PML and the reference.The solid line shows the difference for ␤ 0 ‫ס‬ 1, and the dotted line shows the difference for the optimal ␤ 0 ‫ס‬ 8.The maximum y-axis is about 0.5% of Figure9a.
Martin et al. ͑2008b͒p ␤ , and p ␣ typically range from one to four, and two is commonly used.We use p d ‫ס‬ 2 and p ␤ ‫ס‬ 2, and we choose p ␣ ‫ס‬ 1 becauseMartin et al. ͑2008b͒show the linear variation of ␣ gets a more pronounced decay of energy.The maximum value of d is obtained through a theoretical reflection coefficient relation ͑Collino and Tsogka, 2001͒: