Poynting vector minimal reflection boundary conditions

Abstract

A method for exploring for hydrocarbons, including: simulating a seismic waveform, using a computer, wherein computations are performed on a computational grid representing a subsurface region, said computational grid using perfectly matched layer (PML) boundary conditions that use an energy dissipation operator to minimize non-physical wave reflections at grid boundaries; wherein, in the simulation, the PML boundary conditions are defined to reduce computational instabilities at a boundary by steps including, representing direction of energy propagation by a Poynting vector, and dissipating energy, with the dissipation operator, in a direction of energy propagation instead of in a phase velocity direction; and using the simulated waveform in performing full waveform inversion or reverse time migration of seismic data, and using a physical property model from the inversion or a subsurface image from the migration to explore for hydrocarbons.

Description

FIELD OF THE INVENTION

Exemplary embodiments described herein generally relate to the field of geophysical prospecting for hydrocarbons and, more particularly, to seismic data processing. Specifically, the exemplary embodiments relate to the technical fields of seismic simulation, reverse time depth migration, and full waveform inversion.

BACKGROUND

This section is intended to introduce various aspects of the art, which may be associated with exemplary embodiments of the present technological advancement. This discussion is believed to assist in providing a framework to facilitate a better understanding of particular aspects of the present technological advancement. Accordingly, it should be understood that this section should be read in this light, and not necessarily as admissions of prior art.

The perfectly matched layers (PML) absorbing boundary condition by Berenger (1994) is commonly used to approximate the radiation boundary condition for the sides and bottom of an earth model where the earth model is assumed to have infinite extent but the computational model has finite extent. Waves should not reflect from external boundaries of the computational model that are designated to have the radiation boundary condition.

In the standard form for PML as described by Marcinkovich and Olsen (2003), every derivative normal to an exterior boundary has a wave field dissipation operator applied. Several issues arise with standard PML operators. For general anisotropy, if the group velocity and the phase velocity have different signs for the direction normal to the boundary, PML goes unstable and energy can be amplified rather than attenuated at the boundary. The conventional design of PML dissipates in the direction of phase velocity, not in the actual direction of energy propagation as will be discussed below in the detailed description section. In addition, for elastic or anisotropic elastic wave propagation, thin high-contrast shear-velocity velocity layers on the boundary can create instability due to boundary or interface waves not behaving like body waves and again having a phase velocity vector with a different sign compared to the group velocity vector component normal to the boundary.

An ad hoc fix to the thin high-contrast shear wave velocity anomaly in the boundary zone has been to smooth the shear velocity earth model in the PML boundary zone. The smoother on the shear wave velocity needs to honor the rock physics constraints for stability and can be sensitive to how a full waveform inversion updates parameters near the boundary.

The frequency-domain form of the PML operator (eqn. 1)

$\begin{array}{cc}{\left(\frac{\partial }{\partial x_l}\right)}_{\mathrm{PML}}F\left(\overrightarrow{x}\right)=\left(\frac{1}{1+\frac{{\omega }_l\left(\overrightarrow{x}\right)}{i\omega }}\right)\left(\frac{\partial }{\partial x_l}\right)F\left(\overrightarrow{x}\right)& \left(1\right)\end{array}$ and the non-split PML (NPML) operator (eqn. 2)

$\begin{array}{cc}{\left(\frac{\partial }{\partial x_l}\right)}_{\mathrm{NPML}}F\left(\overrightarrow{x}\right)=\left(\frac{\partial }{\partial x_l}\right)\left(\frac{1}{1+\frac{{\omega }_l\left(\overrightarrow{x}\right)}{i\omega }}\right)F\left(\overrightarrow{x}\right)& \left(2\right)\end{array}$

to be associated with spatial derivative terms are given above. The NPML can be easier to implement, but the results may not be as good as those achieved with the PML. To mitigate reflections from an external earth model boundary, derivative terms in the original set of wave propagation equations are replaced with either PML or NPML derivatives which damp the waves propagating to and from the boundary. In three-dimensional space, l=1, 2 or 3, and the above formulation allows the frequency in a direction normal to the boundary, which may be called the damping frequency parameter to be different from the damping frequency parameter in the other two directions. In the time domain, these operators correspond to a temporal convolution with a damped exponential in time. The difference between the PML operator and the NPML operator is an exchange of the order of the spatial derivative and the dissipation operator. These two operators are transposes of each other. They are not identical because the damping coefficients are spatially dependent. The dissipation operator D_i can accordingly be defined in the frequency domain as

$\begin{array}{cc}{D}_l\left(\omega ,\overrightarrow{x}\right)=\left(\frac{1}{1+\frac{{\omega }_l\left(\overrightarrow{x}\right)}{i\omega }}\right)& \left(3\right)\end{array}$ where ω_i may be called the damping frequency parameter and ω is the frequency of propagation.

SUMMARY

A method for exploring for hydrocarbons, including: simulating a seismic waveform, using a computer, wherein computations are performed on a computational grid representing a subsurface region, said computational grid using perfectly matched layer (PML) boundary conditions that use an energy dissipation operator to minimize non-physical wave reflections at grid boundaries; wherein, in the simulation, the PML boundary conditions are defined to reduce computational instabilities at a boundary by steps including, representing direction of energy propagation by a Poynting vector, and dissipating energy, with the dissipation operator, in a direction of energy propagation instead of in a phase velocity direction; and using the simulated waveform in performing full waveform inversion or reverse time migration of seismic data, and using a physical property model from the inversion or a subsurface image from the migration to explore for hydrocarbons.

BRIEF DESCRIPTION OF THE DRAWINGS

While the present disclosure is susceptible to various modifications and alternative forms, specific example embodiments thereof have been shown in the drawings and are herein described in detail. It should be understood, however, that the description herein of specific example embodiments is not intended to limit the disclosure to the particular forms disclosed herein, but on the contrary, this disclosure is to cover all modifications and equivalents as defined by the appended claims. It should also be understood that the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating principles of exemplary embodiments of the present technological advancement.

FIG. 1 is a flow chart of an exemplary embodiment of the present technological advancement.

FIG. 2 is a flow chart of exemplary method for modifying the PML to reduce computational instabilities at a boundary.

FIG. 3 is an exemplary computer system usable with the technological advancement.

DETAILED DESCRIPTION

Exemplary embodiments are described herein. However, to the extent that the following description is specific to a particular, this is intended to be for exemplary purposes only and simply provides a description of the exemplary embodiments. Accordingly, the present technological advancement is not limited to the specific embodiments described below, but rather, it includes all alternatives, modifications, and equivalents falling within the true spirit and scope of the appended claims.

Perfectly Matched Layer (PML) boundary conditions are used to simulate non-reflecting external boundaries of a computational earth model in applications such as seismic simulation, reverse time depth migration (RTM) and full waveform inversion (FWI). Unfortunately, these boundary conditions have instabilities for several situations of critical interest for these applications. Particularly for the full wave form inversion application, the earth model is regularly and automatically updated and the algorithm will fail if the boundary conditions go unstable due to earth model modifications near the boundary zone. Typically, these instabilities occur whenever the group velocity vector normal to a boundary is in the opposite direction to the phase velocity vector. Two causes of that condition are commonly encountered; seismic anisotropy for certain types, and levels of anisotropy. Thin high-contrast shear velocity zones extending into the boundary zone may also make this happen.

The present technological advancement can provide a stable solution to this problem enabling full wave form inversion to proceed efficiently. For the seismic simulation and RTM applications, the present technological advancement can simplify and stabilize the earth model building process by providing more stable boundary conditions that work properly for a wider range of earth models.

In an exemplary embodiment, the PML/NPML work flow is modified to dissipate the energy associated with wave propagation normal to the boundary. The direction of energy propagation for anisotropic visco-elastic wave propagation is given by the Poynting vector represented by P below, and defined in terms of the stress tensor σ_ij and the particle velocity vector ν_j.P_i=−σ_ijν_j  (4)

The Poynting vector describes energy flow for body waves, interface waves, guided waves and inhomogeneous waves in isotropic and anisotropic media. If a spatial derivative of the Poynting vector is taken normal to a boundary, then by the chain rule of differential calculus, this involves taking spatial derivatives of stress normal to the boundary scaled by particle velocity plus additional terms involving stress scaled by spatial derivatives of particle velocity normal to the boundary. These spatial derivatives normal to the boundary are the ones to be dissipated in the exemplary embodiment. Below is the directional spatial derivative of the i^th component of the Poynting vector in the k^th direction.

$\begin{array}{cc}\frac{\partial P_i}{\partial x_k}=-\left(\frac{\partial {\sigma }_{\mathrm{ij}}}{\partial x_k}{v}_j+{\sigma }_{\mathrm{ij}}\frac{\partial v_j}{\partial x_k}\right)& \left(5\right)\end{array}$

For simulation on a Cartesian grid, then the directions normal to a boundary are the same as the directions spatial derivatives are taken for computing strain rate or divergence of stress. Therefore, those spatial derivatives in directions normal to the boundary are the terms that should be dissipated to ensure small energy flow to or from the boundary. Therefore k will be equal to i for dissipation of the Poynting vector in the k^th direction. This is denoted using the Kronecker delta notation below. This is a frequency-domain equation and the dissipation operator D_k is in the frequency domain.

$\begin{array}{cc}{D}_k\frac{\partial P_i}{\partial x_k}{\delta }_{\mathrm{ik}}=-D_k\left(\frac{\partial {\sigma }_{\mathrm{ij}}}{\partial x_k}{v}_j+{\sigma }_{\mathrm{ij}}\frac{\partial v_j}{\partial x_k}\right){\delta }_{\mathrm{ik}}& \left(6\right)\end{array}$

Converting from particle velocity to particle displacement and applying the Kronecker delta on the right hand side yields the following.

$\begin{array}{cc}{D}_k\frac{\partial P_i}{\partial x_k}{\delta }_{\mathrm{ik}}=-D_k\left(\frac{\partial {\sigma }_{\mathrm{kj}}}{\partial x_k}{v}_j+{\sigma }_{\mathrm{kj}}\frac{\partial v_j}{\partial x_k}\right)=-i\omega D_k\left(\frac{\partial {\sigma }_{\mathrm{kj}}}{\partial x_k}{u}_j+{\sigma }_{\mathrm{kj}}{\varepsilon }_{\mathrm{kj}}\right)& \left(7\right)\end{array}$ where u_j is displacement and ε_kj is the strain tensor. In the frequency domain, these factors are all multiplicative. It is easy to see that applying the dissipation operator in the k^th direction on a stress component that includes index k will dissipate the Poynting vector in the k^th direction. Likewise, the same is true for such an operator applied to a strain component that includes index k.

Any boundary condition needs to be studied in terms of stability. The total strain energy E in a system can be written in terms of either (a) strain and stiffness c_ijkl or (b) stress and compliance s_ijkl.

$\begin{array}{cc}E=\frac{1}{2}{\varepsilon }_{\mathrm{ij}}{c}_{\mathrm{ijkl}}{\varepsilon }_{\mathrm{kl}}=\frac{1}{2}{\sigma }_{\mathrm{ij}}{s}_{\mathrm{ijkl}}{\sigma }_{\mathrm{kl}}& \left(8\right)\end{array}$

For stable rocks, the stiffness and compliance tensors are positive definite. Therefore any method that dissipates either strain or stress (or both) will reduce the total strain energy in the system and have stable characteristics.

Note that each Poynting vector boundary zone dissipation operator applied to dissipate the Poynting vector normal to the boundary can have the following characteristics. (1) The application of the boundary zone dissipation operator to the strain tensor keeps the strain tensor real and symmetric. (2) The application of the boundary zone dissipation operator to the stress tensor keeps the stress tensor real and symmetric. (3) The strain tensor mentioned in (1) can be strain with temporal derivatives or integrals of any order applied and the Poynting vector will still be dissipated. (4) The stress tensor mentioned in (2) can be stress with temporal derivatives or integrals of any order applied and the Poynting vector will still be dissipated.

The kinetic energy (per unit volume) in the system is the following.

$\begin{array}{cc}{E}_{\mathrm{kinetic}}=\frac{1}{2}\rho v_i^2& \left(9\right)\end{array}$

FIG. 1 is an exemplary method that embodies the present technological advancement. Step 102 includes simulating a seismic waveform, using a computer, wherein computations are performed on a computational grid representing a subsurface region, said computational grid using perfectly matched layer (PML) boundary conditions that use an energy dissipation operator to minimize non-physical wave reflections at grid boundaries (i.e., noise or those reflections at the boundaries that do not correspond to reflections from geological objects, as wave should not reflect from external boundaries of the computational model during the simulation). In the simulation, according to step 104, the PML can be modified to reduce computational instabilities at a boundary. Step 106 includes using the simulated waveform in performing full waveform inversion or reverse time migration of seismic data, and using a physical property model from the inversion or a subsurface image from the migration to explore for hydrocarbons. Step 108 includes using a physical property model from the inversion or a subsurface image from the migration to explore for or manage hydrocarbons. As used herein, hydrocarbon management includes hydrocarbon extraction, hydrocarbon production, hydrocarbon exploration, identifying potential hydrocarbon resources, identifying well locations, determining well injection and/or extraction rates, identifying reservoir connectivity, acquiring, disposing of and/or abandoning hydrocarbon resources, reviewing prior hydrocarbon management decisions, and any other hydrocarbon-related acts or activities.

FIG. 2 is an exemplary method for implementing step 104 (i.e., an exemplary method for modifying the PML to reduce computational instabilities at a boundary). Step 202 includes representing direction of energy propagation by a Poynting vector. Step 204 includes causing the dissipation operator to dissipate energy in the direction of energy propagation instead of in a phase velocity direction.

One embodiment of the present technological advancement may be described as follows. ω_i in the PML and NPML operators can be replaced by ω_jk, with j and k referring to the directions associated with the derivatives used to compute either (a) strain or (b) divergence of the stress. n_PML can be defined as the number of points in the PML boundary zone. First, the spatial derivatives are taken to create the strains and then dissipate the strains so that the order of the operators (e.g., whether the damping operator is performed first prior to the derivative (NPML) or the derivative is applied first and then damped (PML)) given matches that for PML. The damping frequency can be computed for derivatives in the j^th direction and the k^th direction as follows. (Here, the PML zone spatial length L_PML is defined in terms of the spatial increment h and the number of grid points in the PML zone as L_PML=h n_PML. The scale factor of 2.25 is a number that works well for many applications; however it is not unique to the technological advancement and other scale factors can be used.

$\begin{array}{cc}{\omega }_{\mathrm{jk}}=2.25v_{\max }\frac{\ln \left(n_{\mathrm{PML}}\right)}{L_{\mathrm{PML}}}\left[{\left(\frac{L_{\mathrm{PML}}-x_j}{L_{\mathrm{PML}}}\right)}^2+{\left(\frac{L_{\mathrm{PML}}-x_k}{L_{\mathrm{PML}}}\right)}^2\right]& \left(10\right)\end{array}$

The two indices jk replace the single index l shown in Eqn. (3) to deal with situations where the stress (or strain?) tensor has non-zero off diagonal elements. In the frequency domain, the dissipation operator D_jk for the present technological advancement can be written as follows.

$\begin{array}{cc}{D}_{\mathrm{jk}}\left(\omega ,\overrightarrow{x}\right)=\frac{1}{1+\frac{{\omega }_{\mathrm{jk}}\left(\overrightarrow{x}\right)}{i\omega }}& \left(11\right)\end{array}$

For a finite difference implementation, the D_jk operators should honor any grid staggering associated with a staggered grid (a grid where different stress components are defined on different grids that are staggered relative to each other) finite difference scheme.

For a frequency domain implementation, dissipation can be applied using the D_jk operator to each Voigt strain and stress component in the boundary zone. This is exactly PML for the normal strains or normal stresses when j equals k. This formulation is not strictly PML for the shear strains or shear stresses where j does not equal k. Instead this is dissipation applied in a way that improves stability for elastic propagation, as no unbalanced torque is ever applied and total strain energy is always decreased. In this formulation, Poynting vectors for energy propagation normal to the boundary are always dissipated for any type of anisotropy.

For a time domain implementation, the following work flow can be used. An additional advantage to the present technological advancement in the time domain compared to standard PML algorithms is that fewer memory variables are used, helping to improve the balance between the memory allocations in the boundary zone versus that in the interior. The temporal convolution of the computed strain with the damped exponential can be computed via the following steps using an associated memory variable.B=e^−ωjk|Δt|A=B−1ϕ_jk^nt=Bϕ_jk^nt−1+Aε_jk ε_jk^pml=ε_jk+ϕ_jk^nt  12

ε is a strain (or stress) and φ is a memory variable. The dissipation operator on strain components and strain memory variables can be written as one matrix operator step. This time domain operator corresponds to the PML algorithm.

$\begin{array}{cc}\left[\begin{array}{cc}B& 1-B\\ B& 2-B\end{array}\right]\left[\begin{array}{c}{\phi }_{\mathrm{jk}}^{n-1}\\ {\varepsilon }_{\mathrm{jk}}\end{array}\right]=\left[\begin{array}{c}{\phi }_{\mathrm{jk}}^n\\ {\varepsilon }_{\mathrm{jk}}^{\mathrm{pml}}\end{array}\right]& \left(13\right)\end{array}$

The matrix D_jk is defined to be the time domain dissipation operator.

$\begin{array}{cc}{D}_{\mathrm{jk}}=\left[\begin{array}{cc}B& 1-B\\ B& 2-B\end{array}\right]& \left(14\right)\end{array}$

The transpose dissipation operator on stress and stress memory variables can also be written as one matrix operator step. This operator corresponds to the NPML algorithm.

$\begin{array}{cc}\left[\begin{array}{cc}B& B\\ 1-B& 2-B\end{array}\right]\left[\begin{array}{c}{\phi }_{{\sigma }_{\mathrm{jk}}}^{n-1}\\ {\sigma }_{\mathrm{jk}}\end{array}\right]=\left[\begin{array}{c}{\phi }_{{\sigma }_{\mathrm{jk}}}^n\\ {\sigma }_{\mathrm{jk}}^{\mathrm{pml}}\end{array}\right]& \left(15\right)\end{array}$

The matrix D_jk^T is defined to be the transpose time domain dissipation operator.

Putting everything together provides a formula for computing density-weighted particle acceleration as a function of displacement and a system of operators that equal its own transpose. This type of operator design is useful for developing wave propagation operators with matching adjoint wave propagation operators for applications related to reverse time depth migration or full waveform inversion. The following equation (16) shows how this can be done, and is applicable to both the time domain and the frequency domain.

$\begin{array}{cc}\rho \left[\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right]=\left[\begin{array}{cccccc}\frac{\partial }{\partial x_1}& 0& 0& 0& \frac{\partial }{\partial x_3}& \frac{\partial }{\partial x_2}\\ 0& \frac{\partial }{\partial x_2}& 0& \frac{\partial }{\partial x_3}& 0& \frac{\partial }{\partial x_1}\\ 0& 0& \frac{\partial }{\partial x_3}& \frac{\partial }{\partial x_2}& \frac{\partial }{\partial x_1}& 0\end{array}\right]\hspace{1em}\left[\begin{array}{cccccc}{D}_{11}^T& 0& 0& 0& 0& 0\\ 0& D_{22}^T& 0& 0& 0& 0\\ 0& 0& D_{33}^T& 0& 0& 0\\ 0& 0& 0& D_{44}^T& 0& 0\\ 0& 0& 0& 0& D_{55}^T& 0\\ 0& 0& 0& 0& 0& D_{66}^T\end{array}\right]\times \left(\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& C_{14}& C_{15}& C_{16}\\ C_{12}& C_{22}& C_{23}& C_{24}& C_{25}& C_{26}\\ C_{13}& C_{23}& C_{33}& C_{34}& C_{35}& C_{36}\\ C_{14}& C_{24}& C_{34}& C_{44}& C_{45}& C_{46}\\ C_{15}& C_{25}& C_{35}& C_{45}& C_{55}& C_{56}\\ C_{16}& C_{26}& C_{36}& C_{46}& C_{56}& C_{66}\end{array}\right)\hspace{1em}\left[\begin{array}{cccccc}{D}_{11}& 0& 0& 0& 0& 0\\ 0& D_{22}& 0& 0& 0& 0\\ 0& 0& D_{33}& 0& 0& 0\\ 0& 0& 0& D_{44}& 0& 0\\ 0& 0& 0& 0& D_{55}& 0\\ 0& 0& 0& 0& 0& D_{66}\end{array}\right]\left[\begin{array}{ccc}\frac{\partial }{\partial x_1}& 0& 0\\ 0& \frac{\partial }{\partial x_2}& 0\\ 0& 0& \frac{\partial }{\partial x_3}\\ 0& \frac{\partial }{\partial x_3}& \frac{\partial }{\partial x_2}\\ \frac{\partial }{\partial x_3}& 0& \frac{\partial }{\partial x_1}\\ \frac{\partial }{\partial x_2}& \frac{\partial }{\partial x_1}& 0\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]& \left(16\right)\end{array}$

Variations on this algorithm can be done to choose wave propagation state vectors as strain and particle velocity, particle velocity and acceleration, particle velocity and stress, or many other combinations that fully describe the wave propagation initial conditions from a combination of seismic wavefields. Another variation can include a small background dissipation term for waves propagating in any direction within the boundary zone.

In all practical applications, the present technological advancement must be used in conjunction with a computer, programmed in accordance with the disclosures herein. FIG. 3 provides an exemplary computer system upon which the present technological advancement may be embodied.

FIG. 3 is a block diagram of a computer system 2400 that can be used to generate the ASO. A central processing unit (CPU) 2402 is coupled to system bus 2404. The CPU 2402 may be any general-purpose CPU, although other types of architectures of CPU 2402 (or other components of exemplary system 2400) may be used as long as CPU 2402 (and other components of system 2400) supports the operations as described herein. Those of ordinary skill in the art will appreciate that, while only a single CPU 2402 is shown in FIG. 3, additional CPUs may be present. Moreover, the computer system 2400 may comprise a networked, multi-processor computer system that may include a hybrid parallel CPU/GPU system. The CPU 402 may execute the various logical instructions according to various teachings disclosed herein. For example, the CPU 2402 may execute machine-level instructions for performing processing according to the operational flow described.

The computer system 2400 may also include computer components such as non-transitory, computer-readable media. Examples of computer-readable media include a random access memory (RAM) 2406, which may be SRAM, DRAM, SDRAM, or the like. The computer system 2400 may also include additional non-transitory, computer-readable media such as a read-only memory (ROM) 2408, which may be PROM, EPROM, EEPROM, or the like. RAM 2406 and ROM 2408 hold user and system data and programs, as is known in the art. The computer system 2400 may also include an input/output (I/O) adapter 2410, a communications adapter 2422, a user interface adapter 2424, and a display adapter 2418.

The I/O adapter 2410 may connect additional non-transitory, computer-readable media such as a storage device(s) 2412, including, for example, a hard drive, a compact disc (CD) drive, a floppy disk drive, a tape drive, and the like to computer system 2400. The storage device(s) may be used when RAM 2406 is insufficient for the memory requirements associated with storing data for operations of the present techniques. The data storage of the computer system 2400 may be used for storing information and/or other data used or generated as disclosed herein. For example, storage device(s) 2412 may be used to store configuration information or additional plug-ins in accordance with the present techniques. Further, user interface adapter 2424 couples user input devices, such as a keyboard 2428, a pointing device 2426 and/or output devices to the computer system 400. The display adapter 2418 is driven by the CPU 2402 to control the display on a display device 2420 to, for example, present information to the user regarding available plug-ins.

The architecture of system 2400 may be varied as desired. For example, any suitable processor-based device may be used, including without limitation personal computers, laptop computers, computer workstations, and multi-processor servers. Moreover, the present technological advancement may be implemented on application specific integrated circuits (ASICs) or very large scale integrated (VLSI) circuits. In fact, persons of ordinary skill in the art may use any number of suitable hardware structures capable of executing logical operations according to the present technological advancement. The term “processing circuit” encompasses a hardware processor (such as those found in the hardware devices noted above), ASICs, and VLSI circuits. Input data to the computer system 2400 may include various plug-ins and library files. Input data may additionally include configuration information.

The present techniques may be susceptible to various modifications and alternative forms, and the examples discussed above have been shown only by way of example. However, the present techniques are not intended to be limited to the particular examples disclosed herein. Indeed, the present techniques include all alternatives, modifications, and equivalents falling within the spirit and scope of the appended claims.

REFERENCES

• 1. Berenger, J., “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics 114, 185-200 (1994); doi:10.1006/jcph.1994.1159, which is incorporated by reference in its entirety.
• 2. Marcinkovich, C., K. Olsen, “On the implementation of perfectly matched layers in a three-dimensional fourth-order velocity-stress finite difference scheme,” Journal of Geophysical Research Solid Earth 108, 2276 (2003), which is incorporated by reference in its entirety.

Claims

1. A method for exploring for hydrocarbons, comprising: simulating a seismic waveform, using a computer, wherein computations are performed on a computational grid representing a subsurface region, said computational grid using perfectly matched layer (PML) boundary conditions;wherein, in the simulation, the PML boundary conditions are defined to reduce computational instabilities at a boundary by steps including, representing direction of energy propagation by a Poynting vector, anddissipating energy that represents noise or reflections at the boundary that do not correspond to reflections from geological objects, with a dissipation operator operating on the Poynting vector, in a direction of energy propagation instead of in a phase velocity direction; andusing the simulated waveform in performing full waveform inversion or reverse time migration of seismic data;obtaining a physical property model from the inversion or a subsurface image from the migration; andusing the physical property model or the subsurface image to explore for hydrocarbons.

2. The method of claim 1, wherein the dissipating energy in the direction of energy propagation comprises dissipating energy associated with terms in the Poynting vector that involve spatial derivatives of particle velocity or stress that are normal to the boundary.

3. The method of claim 2, wherein the Poynting vector Pi is given by Pi=−σijνj  (4)where σij is a stress tensor and νj is a particle velocity vector.

4. The method of claim 1, wherein the subsurface region comprises anisotropic viscoelastic media.