Optimum boundaries for finite-difference modelling of waves

Optimum boundaries for finite-difference modelling of waves

Peter M. Manning* and Gary F. Margrave

Dept. of Geology and Geophysics, University of Calgary,

2500 University Dr. N.W., Calgary, AB, T2N 1N4

manning@geo.ucalgary.ca

ABSTRACT

A technique is described to generate an optimum set of points beyond the boundary of a finite-difference model; points which can be used for finite-

difference operations, and which eliminate reflections from that boundary. The rationale is explained, the mathematics are developed, and the final formulae are given. Examples are shown for the case of a two-dimensional elastic pressure wave. The limitations of the technique are shown to be consistent with the general limitations of simulating continuous waves by finite-difference

techniques.

Introduction

There have been many techniques developed to reduce the reflections from the boundaries of finite-difference models, and simulate the infinite real earth. These techniques are useful for economic reasons, so that a model size can remain small and yet simulate the effects of specified internal boundaries without interference from the model edges.

The basic technique is to provide extra rows and columns of points around the edges of a model. The amplitudes at these points are needed to allow the finite-difference operations to be executed within the model, but cannot themselves be generated by the same techniques because of their edge position. Unique

algorithms, or in some cases unique conditions, must be used to calculate these amplitudes.

The earlier techniques used to reduce boundary effects were called absorbing boundaries, and simulated the effects of having a highly attenuating material around the model. This technique is very practical where the modelling already accounts for viscous effects on the particle motion (Kelly and Marfurt, 1990). The viscosity is simply made very high for several rows and columns around the model’s area of interest.

Another absorbing technique that can be used is to taper, at each time step, the amplitudes toward the model edge by a minimal amount. Cerjan et al. (1985) got very successful results by tapering to a maximum of 0.92 across a boundary

zone of 20 points in width. With absorbing boundaries, the edge-point amplitudes are calculated by an approximate algorithm, but any errors that this introduces is shielded by the attenuating zone. The increased overhead caused by providing the attenuating zone is usually not a major barrier with modern computers.

An alternative to absorbing boundary conditions can be called transmitting boundary conditions. Reynolds (1978) called his boundaries transmitting, and although Clayton and Engquist (1977) called their boundaries absorbing, they used algorithms similar to Reynolds. These algorithms project amplitudes into the boundary zones from the values already calculated for the zone of interest.

Clayton and Engquist adapted a migration algorithm to project boundary values. Reynolds factored the wave equation and then used approximations for finite-differencing. A requirement of these techniques is to select only those solutions that advance into the boundary, and suppress solutions that advance out of the boundary (the reflections). They are found to work very well with waves moving directly toward the boundary, but not so well with waves approaching the boundary at an acute angle.

This paper describes a transmitting boundary solution for the second-order finite-difference elastic wave equation. In this space of digital values, two unknowns must be found. The first unknown is the extra boundary value amplitude, and the second unknown is the advanced time-step amplitude that is calculated using the extra boundary value. The first of the two equations that is required for a solution is, of course, the time stepping equation. We have found that the second

required equation is the one that relates all the first derivatives of an unimpeded advancing wave (the eikonal equation). Any solution that does not satisfy this equation must involve some reflected energy.

The above simultaneous solution takes the form of a quadratic. The root of the quadratic must be chosen so that the slope of the wave toward the boundary is consistent with the slope in time of an advancing wave. In particular, a slope down toward the boundary must accompany more positive amplitudes with time, and vice-versa.

Theory

The development of the theory starts with the definition of a scalar plane-wave, which may be chosen to advance with time

P=F((zcosθ+xsinθ)k?ωt). (1) Then an equation relating the derivatives of the function may be shown to be k2??P?1??P???P???P???+??=2??=2??. v??t?ω??t???x???z?2222(2) This is the well-known eikonal equation. Note that the squaring of the derivatives destroys the sign of an inward or outward wave, so that the selection of an outward advancing wave must be made by choosing the correct root.

The equation of the scalar function P may be translated into a finite-difference version using central differences. If m, n, and k represent the function at x, z, and t respectively, then

1?P(m+1)?P(m?1)??P(n+1)?P(n?1)??=2?+??2?x2?zv????22?P(k+1)?P(k?1)?? (3) ?2?t??2

where P(m,n,k) has the independent variables omitted unless they have been incremented or decremented. At the x border, where m now represents the x edge, the decremented variables can be assumed to represent interior spaces and older times, and are therefore known. Also all the n’s are known in the z direction (n+1 and n-1). The two unknown amplitudes are P(m+1) and P(k+1). If R is defined by

and C=?R2=(P(n+1)?P(n?1)), 2(4) ?t (where ?x = ?z, and generally C≥1 for stability), then

(P(m+1)?P(m?1))2+R2=C(P(k+1)?P(k?1))2. 2(5) The wave equation, ?2P?2P1?2P, +=?x2?z2v2?t2(6)

has a finite-difference representation of the form

P(m+1)?2P+P(m?1)

?x2+P(n+1)?2P+P(n?1)?z2=1P(k+1)?2P+P(k?1),(7) 22v?twhich may be converted into the form

P(k+1)=2P?P(k?1)+1(P(m+1)+P(m?1)+P(n+1)+P(n?1)?4P). (8) 2C

Again, the two unknown amplitudes are P(m+1) and P(k+1).

If D is defined by the equation,

D=2C2P?2CP(k?1)+P(m?1)+P(n+1)+P(n?1)?4P, 2(9) then, when equation (8) is substituted into equation (5) the quadratic equation in P(m+1) which can be derived is

?C2P(m?1)+D?C2P2(m?1)+C2R2?D2

?=0, (10) P(m+1)?2?P(m+1)+22??C?1C?1??2

which can be solved for P(m+1), and using (8), for P(k+1). The root that must be chosen for the right (x) boundary is the one for which

(P(m+1)?P(m?1))?(P(k+1)?P(k?1))≤0, (11) ensuring that when the slope in the x direction is positive, the slope in the time direction is negative, and vice-versa. These are the conditions for a wave advancing into the boundary.

For the special case C = 1, the solution for P(m+1) is

R2

P(m+1)=2P+2P(k?1)?P(n+1)?P(n?1)+. 2Pm?1+D(12)

Fig. 1 shows the relation between the unknowns and knowns for the closely

related finite-difference time-stepping case with only one spatial dimension. Each time step adds a trace on the the lower right edge of the surface. The wave is moving into the boundary on the lower left edge of the surface.

The five points of the one spatial dimension finite-difference equation, equivalent to equation (7), can be found at the ends and centre of the black cross. The time-stepped value at the right limb of the cross is determined from equation (7) and the known amplitudes at the other four points.

The five points of the red cross include two unknowns: the point for the new time as above, and the point at the lower edge from beyond the boundary. The equivalent of equation (5) is added to the usual time-stepping equation for a simultaneous solution for the two output points From The Three Known Points. Application To The Elastic Wave Equation

The boundary conditions specified above have been adapted to the elastic wave equation by converting the elastic wave displacements into two sets of scalar amplitudes near the boundaries. Fig. 2 shows the relative positions of the x and z displacements within the staggered grid at the right side (x) boundary. At the positions between the displacements where the displacements converge (or diverge), a pressure may be found by the formula,

P(m,j)=(Ux(m,j)?Ux(m?1,j))+(Uz(m,j+1)?Uz(m,j)). (13) Similarly, a twist may be fond at the intermediate positions staggered from the pressure positions using the formula,

T(m,j)=(Uz(m,j+1)?Uz(m?1,j+1))?(Ux(m?1,j+1)?Ux(m?1,j)). (14) Note that the relative values of the indices for Ux, Uz, P and T depend on how a particular staggered grid system is defined.

The positions of Ux, Uz, P and T are related graphically in Fig. 2. The

displacements calculated within the normal (interior) steps of the model are shown by the arrows in black. The blue Ps and Ts are calculated from these displacements. The Ps and Ts are estimated at the red positions by the scalar

projection described above. The projected Ts allow the red z displacements to be calculated, and then the Ps may be used to calculate the x displacements Application Example

Figure 3 shows an application of the transmitting boundary for a pressure wave within an elastic model. In the left panel, the directly arriving P-wave has

generated minimal reflections from the boundary. From a more shallow source, the interaction of the transmitting boundary and the free surface boundary has generated an unacceptable point source.

Corners involving transmitting boundaries require special coding to reduce undesirable artefacts (see Clayton and Engquist, 1977). This has not yet been attempted for this modelling code.

Limitations Of The Transmitting Boundaries Code

The techniques developed here were designed for plane waves, and may be found to be less than optimum for non-planar waves and combinations of plane waves with differing basic frequencies. Testing will show if these other waves and combinations cause significant problems.

Application of the transmitting boundary code can also be expected to suppress reflections only to the extent that the finite-difference code simulates the

continuous conditions. In other words, reflections will be suppressed to the same extent that the sample rates chosen can represent the continuous equations (see Manning and Margrave, 2000). Tests done with one-dimensional models confirm this expectation.

Fig. 4 shows a series of tests for application of transmitting boundaries on one-dimensional data. The boundary condition is the one specified in Reynolds, and developed for central differences by a method parallel to that specified above for two-dimensional data. The bottom graph in the figure shows the perfect sampling case where dt = dx/v, and here the reflection is completely suppressed. The

upper three graphs show cases where dt = dx/2*v, or half the dt value needed for the perfect sampling case. These results all show noticeable reflections, but finer sample rates reduce the undesirable effect.

Conclusions

Theory shows that a finite-difference transmitting boundary condition can be specified that will eliminate reflections from a boundary under certain conditions. The conditions are: the incident energy is an isolated plane-wave, the finite-difference scheme is second order, and the elimination is within the limits expected for the spatial and time sample-rates of the finite-difference scheme used.

Practical tests done to date give reason to hope that the proposed method will be useful for most forms of wave energy.