Reflection and transmission (scattering) of plane waves at a planar boundary between two elastic half-spaces are important fundamental processes in seismology. Such plane waves may be compressional (P) or shear (S) in an elastic medium. In this article we apply Mathematica to computing the complex algebraic expressions describing the reflection and transmission amplitudes, phases, and angles of propagation. I also illustrate the energy flux quantities. Brewster’s angle, which represents where zero reflected energy occurs, is an important variable for SH waves, which are polarized out of the cross section, and is dependent on both the velocity and density ratio of the two half-spaces. For the P-SV system, the P and SV waves are polarized in the plane of the cross section. Poisson’s ratio (a function of the P and S velocities) affects the energy flux behavior in the case of incident SV but not incident P. Also for the P-SV system, there are four critical angles that affect the energy flux behavior. I also study the limiting P-SV case where the top medium becomes a vacuum, thus creating a “free surface.” For the P-SV system, energy flux of the four scattered waves for any incident wave can be partitioned into top and bottom layer contributions and into total P and total SV contributions, providing further insight into the nature of P-SV scattering. The single-boundary formulation can easily be extended to a stack of layers giving the amplitude and phase at receivers offset along the surface from the source.

Introduction

The reflection and transmission of waves across boundaries in elastic media is an interesting phenomenon with applications in both global Earth structure and exploration seismology. Waves generally lose energy as they pass through a system of layers, with some energy trapped in the layers. The amount of trapped energy is dependent on the angle of incidence on the stack of layers. The earliest formulas for the reflection and transmission coefficients at a plane boundary between two half-spaces were given by Knott [1] and Zoeppritz [2]. Additional papers added to the understanding of plane waves in layered media, as summarized in [3]. Unfortunately, all of the early papers apparently had one or more algebraic or typographical errors that made their implementation problematic. The work of Nafe [4], however, provided a set of equations with an internal check on the accuracy and consistency; and implementations based on his work will provide correct coefficients. Aki and Richards [5] essentially repeated his formulation with additional clarifications, and many modern applications are based on their published equations. Published examples of modern computed results based on this formulation can be found in [6], for instance. Although some authors have looked at various departures from homogeneous layers and planar boundaries, this article will deal only with plane waves in media of two semi-infinite half-spaces separated by a planar boundary and having homogeneous properties. This article also treats the energy partition between adjacent half-spaces, showing the partition of P and SV energy across a boundary and investigating the zeros for reflected energy (known as Brewster’s angle for compressional waves in liquid media). Also, current Mathematica computational capability allows the amplitudes or energy to be visualized not with just two-dimensional line graphs showing results for one set of media parameters, but with three-dimensional density plots showing results over a range of media parameters. These results are newly presented using the existing theory, with no new theoretical results. Due to the complexity of the equations, these newly presented results have not been sufficiently explored before. Mathematica has been the enabling factor. Several .m files of functions for computing reflection and transmission responses are available at content.wolfram.com/sites/19/2012/01/von_Seggern.zip (see Erratum). These files are used to develop this article throughout. Readers of this notebook should install these files and augment their $Path variable if necessary so that they can be read in when evaluating the computation cells of this article. The Mathematica code contained in this article was developed using Mathematica 8.0.

Erratum

Prior to May 8, 2012, the downloadable .m file AVOfunctions.m had an algorithm flaw that may cause earth models other than the one used in the article to fail in the raytracing function rayparms when rays become evanescent; that is, go beyond incidence angles at which there is complete reflection. A corrected AVOfunctions.m file has been provided, and any readers who downloaded the original version should download the current version.

General Approach

I describe the boundary problem as that of two elastic half-spaces, top and bottom (also called layer 1 and layer 2, or top layer and bottom layer), each of homogeneous material and separated by a planar boundary. The top and bottom layers are described with velocity and density parameters: for compressional velocity, for shear velocity, and for density. An important parameter, called Poisson’s ratio, relates the compressional and shear velocities thus:

(1)

A special case, often assumed in practice, is that , thus giving . A medium having this property is called a Poisson solid. When , the medium is a fluid with . I will allow the constant relating and to vary in each half-space, thus allowing the Poisson ratio to vary.

Consider an incident wave whose phase angle is constant along a plane (a “plane wave”) propagating in these media at an arbitrary angle. Let the vector normal to the wavefront define a vertical plane that forms a cross section of the media. Let the normal to the boundary be the axis, positive upwards, and the horizontal along the boundary be the axis, positive in the rightward direction. We no longer consider the dimension of the media, as the action can be described in the – plane. Compressional waves (P) have their displacement in this vertical plane. Shear waves have their displacement either in this plane (vertical shear, or SV) or transverse to it (horizontal shear, or SH). I treat only isotropic media, in which a wave propagates with constant velocity in any direction. For an isotropic medium, the SV and SH waves are decoupled while the P and SV waves have a complex interaction when they are scattered at the boundary. Thus, incident SH only scatters into SH on both sides of the boundary while incident P or SV scatters into both P and SV on both sides.

I adopt the equations and notation of [5] for most of this article, but equations 5.35 through 5.40 in [5] for the reflection and transmission coefficients of P-SV plane waves. These equations can be utilized to compute the coefficients at arbitrary angles of incidence. They describe the 16 possible interactions at an interface between P and SV waves. The incident wave may be taken as one of four cases: downgoing P (Pd), downgoing S (Sd), upgoing P (Pu), or upgoing S (Su). Each incident wave case gives rise to two reflections (P and S) and to two transmissions (P and S), giving four scattered waves. A similar treatment is given in [5] for horizontally polarized S waves (SH), and I use equation 5.33 from [5] for the four possible reflection/transmission coefficients. In the SH case, no conversions take place, and the results are fairly simple compared to the P-SV case and akin to compressional waves in fluid-fluid media. In general, computed coefficients are complex quantities. In all cases, one can obtain the amplitude by taking the absolute value of the computed coefficient and the phase by taking the argument of the computed coefficient.

A visual rendering of the reflection/transmission coefficients is obtained by putting the expressions in the Plot function of Mathematica. Prior to doing that, the expressions can be greatly simplified by realizing that the four velocities involved for P-SV can be reduced to three simple ratios and that the two densities involved can be reduced to one ratio. Mathematica enables one to easily make the substitutions and cancellations in the algebraic expressions to achieve simpler expressions. These simpler expressions, even though containing many terms, can still be quickly evaluated and plotted by the function Plot. Similarly, the SH expressions can be simplified using a velocity ratio and a density ratio; the resulting expressions are very simple compared to the P-SV ones.

I define vector wave slowness as the inverse of velocity, , where is any of the four velocities involved. The vertical component of the slowness is termed and is given by

(2)

where is the angle measured between the axis and the vector normal to the wavefront. The parameter is used throughout in the expressions of [5].

An internal check exists for the correctness of the reflection/transmission coefficients for any particular incident wave, based on the requirement that the sum of the vertical energy fluxes on a unit area of the horizontal boundary of the four scattered P-SV waves must equal the vertical energy flux of the incident wave on a unit area at that same boundary. In [5] this relation is given for downgoing P in equation 5.42. By evaluating equations of this type for any incident wave and for all possible horizontal slownesses, one can verify that the computed reflection and transmission coefficients are correct. A similar energy flux equation pertains to the SH case. Expressions for these energy quantities can also be plotted directly using Plot.

The propagation direction of the scattered waves relative to that of the incident wave is governed by Snell’s law. Let be the velocity of the incident wave and be the velocity of the scattered wave. Then Snell’s law is

(3)

A fundamental limit in all the computations involves the maximum of the vertical slowness of a wave, defined as where , or . For a scattered wave, depending on its exact type, the limit of will be reached at a smaller angle of the incident wave when of the scattered wave is larger than , or when . At that value of , the scattered wave becomes “evanescent” and has zero theoretical amplitude while the incident wave still exists. This occurs for 0, 1, 2, or 3 of the four possible scattered waves, with consequent effects on all the reflection/transmission coefficients, often quite abrupt. In order to get complete results, scattering computations need to be made up to and including the largest possible value of , which is the inverse of the maximum of the four velocities. Results past the evanescent limit for a particular scattered wave return Indeterminate values from Mathematica. For energy computations, the evanescent waves have zero vertical energy flux across the boundary and no longer contribute to the energy flux equations.

I start by examining SH phenomena because the equations are fairly simple and the results are nearly intuitive. Then I move to the P-SV case, where the scattering phenomena become much more complicated. Lastly I consider the limiting P-SV case when the upper medium becomes a vacuum, thus creating a “stress-free boundary” above which no waves exist.

To use the functions in this article, make sure the $Path variable will point to the files RTfunctions.m and AVOfunctions.m and load them.

SH Waves in Elastic Media

SH Amplitude and Phase

Consider now the amplitude and phase of the two scattered waves at the boundary relative to that of a downward incident SH wave. The downward incident wave is labeled Sd for S down. The scattered waves are either Su (for S up) or Sd. Similarly, an upward incident wave is designated Su, again having scattered waves Su or Sd. Then, for instance, an interaction is labeled SdSu for upward reflection of a downward incident wave. The algebraic results for the four scattered waves are given as equation 5.33 in [5] and repeated here:

(4)

(5)

(6)

(7)

There are four layer parameters here: , , , . For the sake of compactness and generality, I make these into ratios and . The expressions (equations 4–7) are then converted to ones using only ratios and in the Mathematica code. Also the cosine functions are changed to sine functions. The denominator is identical for all four coefficients. The angles and refer to the top scattered-wave angle and the bottom scattered-wave angle, respectively. Figures 1 and 2 show the graphs of the results of evaluating these expressions in the function plotSHcoef for an incident Sd wave and an incident Su wave for given and . This function takes the following arguments:

1. incident wave ("Sd" or "Su")
2. density ratio r (bottom over top)
3. velocity ratio s (bottom over top)

The plotting requires a careful handling of the critical angles, as indicated by the Mathematica code in the file RTfunctions.m, in order to give a neat presentation. Recall that the coefficients are computed up to the maximum incident angle of for all scattered waves and that some of the results are therefore indeterminate, depending on how the two velocities relate to one another.

Figure 1. Amplitude and phase results for scattered waves for an SH wave incident downward on the boundary between elastic half-spaces.

Figure 2. Amplitude and phase results for scattered waves for an SH wave incident upward on the boundary between elastic half-spaces.

There are two important angles in the SH results. Consider the case where the top velocity is less than the bottom velocity. (The opposite case can easily be generated as well.) In Figure 1, as the incident wave increases its angle from zero (straight downward), the first angle of interest is Brewster’s angle [6], where no reflected wave exists (), and the energy of the incident wave is entirely transmitted. Brewster’s angle can be computed given the media parameters, as shown below. The second important angle is the critical angle where the transmitted wave becomes evanescent (). This angle is always greater than Brewster’s angle and is the incident angle beyond which no energy is transmitted to the layer of higher velocity. An SdSu phase shift of occurs at Brewster’s angle, and then a gradual phase shift begins at the critical angle, eventually coming to a full change from the initial phase angle at 90° incidence angle. The transmitted wave undergoes no phase shift, though. Figure 2 shows reciprocal results for an Su incident wave, now from the higher-velocity medium below. There is no critical angle for the transmitted wave, but Brewster’s angle still remains in the reflected wave (), with its accompanying phase shift. The sum of the Brewster’s angles for SdSu and SuSd is equal to . There is no phase shift in the transmitted wave for any angle of the incident Su wave. Similar plots can be made for any specified ratios of layer parameters by changing and in the input cells.

Brewster’s Angle

The angle at which there is no reflected energy is called Brewster’s angle in electromagnetic theory [7]; it is normally given in terms of the refraction indices by:

(8)

where and are the refractive indices for layers 1 and 2. However, density is not a factor in the electromagnetic formulation, and this angle has a more complicated formula for elastic media with variable density. It can, in fact, be derived from the reflection coefficient expression for SdSu in equation (4). Setting the numerator to zero (Brewster’s condition of zero reflection), we get

(9)

Square the above and apply the trigonometric identity to get

(10)

Multiplying out the denominator on the RHS and rearranging, we get

(11)

Now apply Snell’s law (equation 3) to the first term on the LHS of equation (11) to obtain:

(12)

Now let and rearrange to get

(13)

and finally, letting this be given as , we get

(14)

Letting simplifies the expression to

(15)

For the special case of , we get

(16)

Factoring out , we obtain

(17)

and, making the obvious change to the arc tangent and substituting back in , we obtain finally

(18)

for a medium with no density contrast. This is the same as the electromagnetic solution (8). An interesting result, true only for , is that for a given is complementary to that for , such that , as shown in Figures 1 and 2.

Using equation (3) (Snell’s law), the critical angle , the angle at which the refracted wave becomes evanescent and the normal to the wavefront becomes theoretically horizontal, is

(19)

We can now plot the critical angle and Brewster angle versus velocity ratio as in Figure 3. The input function plotBrewster allows one to modify this plot as desired, with the following arguments:

1. minimum velocity ratio (bottom over top)
2. maximum velocity ratio (bottom over top)

Results for three density contrasts of [0.9, 1.0, 1.1] are shown in Figure 3, but these could also be modified in the RTfunctions.m file. Brewster’s angle follows a connected line, as given by equation (18), only when there is no density contrast (). At , the angle is exactly 45°. Brewster’s angle for cases where there is a density contrast is more complex, as given by equation (15), and it is discontinuous at . Figure 3 also shows the critical angle (19), which reaches 90° when there is no velocity contrast (). The critical angle and Brewster’s angle converge for higher velocity ratios and theoretically approach, but do not equal, zero.

Figure 3. Brewster’s angle (dashed) and the critical angle (solid) for SH waves in two half-spaces with a variable velocity ratio. Brewster’s angle is shown for three density contrasts of 0.9, 1.0, and 1.1 in red, green, and blue.

SH Energy

We can take the pattern of equation 5.42 in [5] for the P-SV system and apply it to SH waves to compute energy flux terms, as defined above. In the SH case, only two terms remain on the RHS: one for the energy flux of the scattered SH wave above the boundary and the other for the energy flux of the scattered SH wave below the boundary. For SH incident from the top and bottom, the energy flux equations are:

(20)

(21)

where the LHS represents the incident wave. Dividing by the LHS variables, we obtain

(22)

(23)

I call these the normalized energy fluxes because the LHS incident wave term is unity. Hereafter, I use the term “energy flux” for normalized energy flux. Thus, the combined energy flux of the scattered waves, in this case that of the reflected wave plus that of the transmitted wave, should equal unity, regardless of angle of incidence. In [6], for instance, one finds computed graphs exemplifying these relations versus incidence angle of the incident wave for given medium parameters. I treat the cases of an incident wave in the top and bottom layers here.

The input cells below transform these expressions to ones using the ratios and . The function plotSHflux allows one to vary the plot with the following arguments:

1. incident wave ("Sd" or "Su")
2. density ratio r (bottom over top)
3. velocity ratio s (bottom over top)

This approach transforms the expressions for the coefficients before plotting them. The absolute values of the coefficients must be squared, and an additional factor must be multiplied into the first terms on the RHS of equations (22) and (23) to get the correct energy flux terms for SdSd. Figure 4 shows plots of the energy flux of the reflected and transmitted SH waves for an incident Sd wave. Figure 5 shows the results for the same media, but now having an incident Su wave. To simplify the calculations, I have assumed the density is equal in the two layers (); there is, however, no computational penalty for using . In both cases, the two plotted curves are complementary and therefore sum to unity at any point in order to satisfy the energy flux equations. Figures 4 and 5 show the location of the critical angle and of Brewster’s angle. Again, Brewster’s angle is where the reflected energy in the incident layer goes to zero. The critical angle bounds a sudden increase or decrease in energy flux for the scattered waves, occurring after the zero in the reflected energy.

Figure 4. Energy partition into the top and bottom layers for an SH wave incident downward on the boundary between two elastic half-spaces.

Figure 5. Energy partition into the top and bottom layers for an SH wave incident upward on the boundary between two elastic half-spaces.

SH Energy for Variable Media Parameters

Here we take advantage of the computing capability in Mathematica to show the energy relations as plots over two variables: incidence angle and ratio of velocity in the bottom layer to that in the top layer. I have shown that the terms in the RHS of equation (22) for incident Sd can be parameterized in and . I again assume that the density is equal in both the top and bottom half-spaces, meaning ; thus the propagation medium is characterized by just one parameter (). Then the two terms that sum to the unit energy flux in the RHS of equation (22) can each be plotted against the velocity ratio and the incidence angle with the Mathematica function ContourPlot in order to see how the energy flux behaves as the medium is varied. This can be done in the following section, and Figure 6 is a representative example of the type of plot that it produces. The function plotSHenergy has the following arguments:

1. incident wave type ("Sd" or "Su")
2. density ratio r (bottom over top)
3. minimum velocity ratio (bottom over top)
4. maximum velocity ratio (bottom over top)

The two plots are mutually complementary, summing to unity at every point as required by equation (22). (This evaluation can take minutes.)

Figure 6. Energy partition into the top and bottom layers for an SH wave incident downward on the boundary between two elastic half-spaces. The velocity ratio of the two layers is the vertical axis variable, and the angle of incidence for the SH wave is the horizontal axis.

Figure 6 could be overlaid with the Brewster’s angle curve for and the critical angle curve of Figure 3. The critical angle curve bounds a sudden increase or decrease in energy flux for the scattered waves. The locus of the Brewster’s angle would track the zero contour of the reflected SH energy flux, although it is not visually clear exactly where this lies, due to the low resolution of the color coding scheme. Conversely, the transmitted wave’s energy flux is unity along this locus. The horizontal slice at gives the result for a medium with no velocity contrast, for which the energy flux of the reflected wave is always zero and that of the transmitted wave always unity. Figure 6 shows that for a medium of low velocity contrast, the reflected SH wave generally has less energy flux than the transmitted SH, except past the critical angle when . Thus, SH waves propagating though layered media with interfaces of fairly low contrast will be mostly transmitted until the incidence angle becomes large; then large-angle reflections with significant amplitude abruptly appear.

A horizontal slice through both plots in Figure 6 at any specific velocity ratio would give a pair of two-dimensional line graphs like Figure 4, which showed the complementary nature of the top and bottom energy flux. Until now only these two-dimensional line graphs have been shown in the literature, but Mathematica easily allows for more informative plots such as Figure 6.

P-SV Waves in Elastic Media

P-SV Amplitude and Phase

The case of P and SV waves in elastic media is much more complex than the case of SH waves. There are now four scattered waves for each of the four possible incident waves, a total of 16 possible interactions. The downgoing waves are now Pd and Sd, and the upgoing ones are Pu and Su. There are now four critical angles that influence reflection and transmission. Consider first the case of incident Pd waves in the top layer. Two critical angles will be expressed in terms of the incidence angle of the Pd wave. The first is the PdPd critical angle and is only relevant when the P velocity of the bottom layer exceeds that of the top layer:

(24)

The second is the PdSd critical angle and is only relevant when the S velocity of the bottom layer exceeds the P velocity of the top layer:

(25)

Now consider the incident Sd wave in the top layer. The third critical angle is the SdSd critical angle and is only relevant when the S velocity of the bottom layer exceeds the S velocity of the top layer:

(26)

The fourth is the SdPd critical angle and is only relevant when the P velocity of the bottom layer exceeds the S velocity of the top layer:

(27)

Similar relations hold if we consider incident waves in the bottom layer. When Poisson’s ratio for both layers is equal, equations (24) and (26) give the same curve. The function plotPSVcoef will plot the reflection and transmission coefficients for the P-SV system. The arguments to this functions are:

1. incident wave ("Pd", "Pu", "Sd", or "Su")
2. density ratio r (bottom over top)
3. velocity ratio s (bottom over top)
4. velocity ratio v (compressional to shear in top layer)
5. velocity ratio w (compressional to shear in bottom layer)

Typical results for reflection/transmission coefficients for the P-SV system are shown in Figure 7 for an incident Pd wave in a medium with high velocity contrast and in Figure 8 for a incident Sd wave in the same medium. The graphs of amplitude are much more complex than for the typical SH examples of Figures 1 and 2. The critical angle for P or S at approximately 30° is evident, however. For the P-SV case, there is nothing equivalent to Brewster’s angle for SH. Plots similar to Figures 7 and 8 can be made for upward incident Pu or Su waves using the same input cells.

Figure 7. Amplitude and phase results for scattered waves for a P wave incident downward on the boundary between two elastic half-spaces.

Figure 8. Amplitude and phase results for scattered waves for an SV wave incident downward on the boundary between two elastic half-spaces.

P-SV Energy

The pattern of equation 5.42 in [5] for the downgoing P for a P-SV case can be extended easily to the remaining three possible incident waves, and the energy fluxes in the top and bottom layers can be computed as for the SH wave above. In this case, four terms exist on the RHS, two for the P and SV in the top layer and two for the P and SV in the bottom layer. The four energy flux equations for cases of P and SV incident waves in the top layer and then P and SV incident waves in the bottom layer are:

(28)

(29)

(30)

(31)

Line graphs of the energy relationships for given medium parameters are shown in [6]. As for the reflection/transmission coefficients, I convert these expressions, using Mathematica manipulations, to expressions involving , , , and , all previously defined. The function plotPSVflux computes and plots the energy fluxes; it has the following arguments:

1. incident wave ("Pd", "Pu", "Sd", or "Su")
2. density ratio r (bottom over top)
3. velocity ratio s (bottom over top)
4. velocity ratio v (compressional to shear in top layer)
5. velocity ratio w (compressional to shear in bottom layer)

Again, the sum of all the energy fluxes must equal unity, regardless of the layer parameters and the angle of incidence of the incident wave. Here I show typical results in Figures 9 and 10 for the cases of an incident Pd wave and an incident Sd wave. In Figures 11 and 12, I have made partitions of the energy flux into top and bottom layers and then again into P and S for the Pd and Sd incident waves using the function plotPSVfluxp, which has the same arguments as plotPSVflux. In Figures 11 and 12, the energy fluxes in the top and bottom layers are complementary, summing to unity at any angle of incidence, as required by equations (28–31). Similarly, in the same figures, the partition of the four energy fluxes into P and SV, rather than top and bottom, gives curves that are complementary.

Figure 9. Energy flux in scattered waves for a P wave incident downward on the boundary between two elastic half-spaces.

Figure 10. Energy flux in scattered waves for an SV wave incident downward on the boundary between two elastic half-spaces.

Figure 11. Energy flux partition in top and bottom scattered waves and into P and SV scattered waves for a P wave incident downward on the boundary between two elastic half-spaces.

Figure 12. Energy flux partition in top and bottom scattered waves and into P and SV scattered waves for an SV wave incident downward on the boundary between two elastic half-spaces.

P-SV Energy for Variable Media Parameters

Here again, as for SH above, I take advantage of Mathematica computing capability to show these relations as energy flux plots over two variables: incidence angle and ratio of velocity in the bottom layer to that in the top layer. Taking the expressions formed in the input cells for Figures 9 and 10 but now setting only while allowing to vary, Figures 13 and 14 show the partition of energy flux terms for the four scattered waves, assuming and (Poisson solids). The function plotPSVenergy has the following arguments:

1. incident wave type ("Pd", "Sd", "Pu", or "Su")
2. density ratio r (bottom over top)
3. velocity ratio v (compressional to shear in top layer)
4. velocity ratio w (compressional to shear in bottom layer)
5. minimum velocity ratio smin (bottom over top)
6. maximum velocity ratio smax (bottom over top)

Again, as for SH, the assumption that is not a computational requirement and neither is the assumption . The correctness of the results can be partially checked by summing all the terms of equation (28) or (29) and plotting the total flux (not shown). For both the incident Pd and incident Sd, such plots appear as solid color, indicating unity over all incidence angles and velocity ratios. (The evaluations in this subsection can take minutes.)

Figure 13. Energy flux for the four scattered waves of a P wave incident downward on the boundary between two elastic half-spaces. The velocity ratio of the two layers is the vertical axis variable, and the angle of incidence for the P wave is the horizontal axis.

Figure 14. Energy flux for the four scattered waves of a SV wave incident downward on the boundary between two elastic half-spaces. The velocity ratio of the two layers is the vertical axis variable, and the angle of incidence for the SV wave is the horizontal axis.

Figures 13 and 14 for incident Pd and incident Sd show the energy flux for the four scattered waves from each incident wave. Horizontal slices through any of these plots at any specific velocity ratio would give a line graph to compare with the simpler calculations shown in Figures 9 and 10. This is another check on the correctness of the algebraic expressions for PdPu etc. or SdPu etc.

In Figures 13 and 14, three curves are shown in gray. These are related to the four critical angles (24–27). Due to the fact that the ratio of compressional-to-shear velocity was assumed to be identical for the two layers (), the two curves from equations (24) and (26) overlap exactly; they separate off the area for which there is no transmission to the bottom layer for a wave of the same type as the incident wave. The curve to the left corresponds to equation (27), and it has an obvious effect on the contours for incident Sd only. The curve to the right corresponds to equation (25), and this indicates where the S wave in the bottom layer becomes evanescent for an incident P wave in the top layer. The dotted vertical line is at the incidence angle given by , or equivalently in this case. This angle is invariant with the ratio of bottom to top velocity if as in this case. For a Poisson solid (), this angle is , as shown here. The energy flux magnitude is suddenly decreased or increased near this line, but only for incident Sd; there is no effect for incident P. The reflected energy flux is always finite along this line except, of course, at the one point where the medium has no velocity contrast. The effect associated with Poisson’s ratio remains for cases of as the vertical line and its associated effect on the plotted contours move left or right with the exact value of .

Similarly to the graphs shown in Figures 11 and 12, partitions of the P-SV energy flux can be made into total reflected and total transmitted flux and into total P and total S energy flux. Each partition involves summing two particular terms of equations (28) or (29) and then summing the remaining two terms. The function plotPSVenergyp does this, and its arguments are the same as for plotPSVenergy. Figure 15 shows the reflected , transmitted , total S, and total P energy for an incident Pd wave, and Figure 16 shows the same for an incident Sd wave. Because the total energy flux must equal unity, any two horizontally adjacent plots in either Figure 15 or Figure 16 are complementary to one another. This provides a further check on the correctness of the extremely complex algebraic expressions.

Figure 15. Energy flux partitions for the four scattered waves of a P wave incident downward on the boundary between two elastic half-spaces. The velocity ratio of the two layers is the vertical axis variable, and the angle of incidence of the P wave is the horizontal axis.

Figure 16. Energy flux partitions for the four scattered waves of an SV wave incident downward on the boundary between two elastic half-spaces. The velocity ratio of the two layers is the vertical axis variable, and the angle of incidence of the SV wave is the horizontal axis.

P-SV Amplitude and Phase at a Stress-Free Boundary

It is worthwhile to investigate one special case for the P-SV interaction at a stress-free boundary (also called a free surface). This occurs when one of the half-spaces has elastic parameters equal to zero. This serves as a good approximation to conditions at the surface of the Earth and is of broad importance in seismology. Stress, but not displacement, must go to zero at such a boundary, and this greatly simplifies the interactions. Letting the top half-space be the vacuum, I will then consider only the upcoming P and SV waves at this boundary. Each scatters into P and SV waves, so that a matrix of coefficients is sufficient to describe the interactions. For an SH wave, the interaction becomes even simpler, with only the reflected SH wave to consider. The reflected SH from a stress-free boundary has an amplitude equal to the incident wave.

The required expressions are given in [5] as equations 5.27, 5.28, 5.31, and 5.32. Using Mathematica‘s algebraic manipulation functions, one can recast those expressions into expressions that depend only on the ratio of the P to S velocities of the elastic medium, or parameter . The density does not enter into the expressions at all. The function plotPSVfs plots the coefficients; the arguments to this function are:

1. incident wave ("Pu" or "Su")
2. velocity ratio u (compressional-to-shear velocity in the bottom elastic layer)

Figures 17 and 18 show typical results for the incident P wave and incident S wave. The SuPd amplitude becomes relatively large just before the Pd becomes evanescent. The maximum value increases inversely with .

Figure 17. Amplitude and phase results for scattered waves for a P wave upward incident on the boundary between a vacuum and an elastic medium.

Figure 18. Amplitude and phase results for scattered waves for an S wave upward incident on the boundary between a vacuum and an elastic medium.

Amplitude-Versus-Offset (AVO) in Layered Media

Introduction to AVO

Although the study of reflection and transmission at planar boundaries between half-spaces is interesting and can have fairly complex formulations, as shown in the previous sections, seismology tends to look at more realistic media with multiple such boundaries representing the real Earth. The total response of such layered media to waves originating at or below the surface and recorded at or below the surface is of widespread interest. Much effort has been expended to make robust computer programs to handle the total response for not only layered media but also media with properties varying irregularly in three dimensions. However, here we merely exploit the computational efficiency for the two-layer problem developed in previous sections to a multilayer medium with a free surface and a specified reflector boundary at depth. Solving this problem has been necessary in the analysis of amplitude versus offset (AVO) in the seismological literature. “Offset” is defined as the horizontal distance along the free surface between the source and the receiver. For a review of AVO, consult [8] and the accompanying articles. It is well studied for structures over a broad scale, from meters in shallow seismic prospecting to hundreds of kilometers in deep-Earth seismology. It is a subject with many complexities which cannot be handled within the formulations of this article; for instance, anisotropy of real Earth layers, inelasticity of real Earth layers, and spherical spreading of waves emanating from a point. Yet the approach using just the basic formulations of this article gives a good first-order approximation and is therefore useful.

Therefore, I use a layered medium where all the boundaries are parallel and where the top boundary is a free surface; the layers can have variable thickness. The expressions to compute the reflection and transmission coefficients have already been developed, and it remains only for the serial interactions to be accumulated between the surface and a reflecting boundary at depth and then back to the surface for a series of specified offsets. If one defines a “ray” as being normal to the wavefront, then a ray which starts out at a particular angle with respect to the axis can be traced down to the reflector and back up via Snell’s law (3). Snell’s law can be extended, layer after layer, by a series of equalities. Each side in these equalities can be expressed as the unvarying ray parameter (equation 2).

AVO in the SH Case

Consider first the SH-wave problem. A particular raypath is developed along a constant . The reflection at the specified layer boundary is symmetrical about the axis, and the ray returns to the surface along an upward path symmetrical to the downward path. The amplitude of the arrival at the free surface at a particular offset is governed by a series of downward transmission coefficients, say , and a similar series of upward transmission coefficients. In between the two transmission series is the reflection at the specified boundary, say . Expressions for these coefficients have already been developed in a previous section. Each interaction contributes multiplicatively in the final result, which is complex, giving amplitude and phase. Numbering the boundaries with index such that the first boundary lies at the bottom of the first layer etc., and defining the reflecting boundary index as , the final result will be termed SS for an S wave reflected as an S wave:

(32)

In applying equation (32), one must use one set of angles on the way down for incident waves and a different set on the way up for incident waves, both given by Snell’s law. The final amplitude is , given by in Mathematica, and the final phase is , given by in Mathematica. The geometry of a specific raypath can be traced from boundary to boundary using Snell’s law to bend the raypath and using the given layer thicknesses. The difficult part of the computation is finding the raypaths that end at predetermined offsets. The computation makes use of Mathematica‘s built-in functions Interpolation and FindRoot. Then the AVO amplitudes and phases are computed at each boundary along these raypaths. Individual boundary-interaction amplitudes are multiplied to get the final amplitude, and boundary-interaction phases are summed to get the final phases.

The following input cells implement these AVO calculations and plot the results in Mathematica. This notebook also contains several auxiliary steps, such as reading an Earth model, plotting the Earth model, plotting the raypaths, and even constructing a “section” of time series at regular offsets wherein the amplitudes and phases of the AVO results can be seen. Typical AVO results are shown here for an Earth model developed for the Great Basin of Nevada [9]. The model data is retrieved with the function getmodel, which takes the model file name as its one argument. The model parameters are returned as a list. Then the model is plotted with plotmodel, which has the following arguments:

1. model parameters returned by getmodel
2. "pvel" or "svel" for compressional or shear velocity
3. deepest layer to plot (numbered from the top layer)
4. a label for the plot
5. the maximum offset horizontally, from source to receiver
6. the vertical exaggeration to apply to the plot

This Earth model goes to approximately 400 km depth, but here we only look at the reflection from the Mohorovic (“Moho”) discontinuity at 35 km depth, the bottom of the third layer. Figure 19 illustrates the model for the three layers down to the Moho. This figure shows compressional velocities—shear velocities are approximately a factor of less for each layer.

Figure 19. An Earth model, comprised of flat layers, for which AVO calculations will be done. Velocity versus depth is color coded as shown.

It is informative to plot the rays, which are the normals to the wavefront from source to receiver. The function rayparms first determines the ray parameters; its arguments are:

1. model parameters returned by getmodel
2. wave ("PP", "SS", "PS", or "SP"), which is composed of two letters, the first for the downgoing ray and the second for the upcoming ray
3. number of rays desired, regularly spaced from zero to maximum offset
4. the maximum offset horizontally, from source to receiver (km)
5. deepest layer through which rays propagate (the reflection occurs at the bottom of this layer)

The function plotrays then plots the rays; its arguments are:

1. model parameters returned by getmodel
2. wave ("PP", "SS", "PS", or "SP"), which is composed of two letters, the first for the downgoing ray and the second for the upcoming ray
3. deepest layer through which rays propagate (the reflection occurs at the bottom of this layer)
4. a label for the plot
5. the maximum offset horizontally, from source to receiver (km)
6. the vertical exaggeration to apply to the plot
7. pvso, the list returned by rayparms

The exact SS ray geometry corresponding to the model in Figure 19 is shown in Figure 20. The source lies at the left of the model, and rays propagate rightward to receivers at the specified offsets. The vertical scale has been expanded by the same factor as the model plot in order to show the rays clearly; as a consequence, the ray angles are no longer displayed correctly. When seismometers are set out to record the motions in the real Earth, they usually have three components of motion in a Cartesian system. The vertical () motion is registered on the vertical component. The horizontal () motion is registered on a horizontal component aligned along the line of receivers; this is often called the radial component. The horizontal motion () normal to the – plane is recorded on the second horizontal component, often called the tangential component. Only the SH source produces waves which have tangential motion.

This computes and plots the rays. The input cell to plot the model (see above) must be evaluated first.

Figure 20. Rays for the SS reflection in the Earth model of Figure 19. Vertical exaggeration distorts the true angles.

The function plotSHAVO then graphs the amplitude versus offset of the arrivals, taking into account all the transmission coefficients and the one reflection coefficient at the bottom of the specified reflection layer. The arguments for plotSHAVO are:

1. model parameters returned by getmodel
2. deepest layer through which rays propagate (the reflection occurs at the bottom of this layer)
3. a label for the plot
4. the maximum offset horizontally, from source to receiver (km)
5. pvso, the list returned by rayparms

Figure 21 shows the computed AVO results evaluated at discrete increments of 10 km to 500 km (ListPlot joins the points, though). The top two graphs show the amplitude and phase of SS, the bottom-left graph shows the travel time of the rays versus offset, and the bottom-right graph shows the takeoff angle of the rays at the source with respect to vertical. The amplitude and phase results are for a plane wave and so do not include geometrical spreading effects that would accompany a point source. The functions getmodel (above) and rayparms (above, with ) must be run before calling plotSHAVO.

Figure 21. AVO calculations to a maximum offset of 500 km for the Earth model shown in Figure 19. The reflection is at 35 km depth.

In order for the reader to fully visualize the results, I construct wavelets arriving at each receiver and place them in a “section” showing what would be recorded at each offset from the single source location. Realize again that this approach leaves out several important effects which are present with real wave propagation in the Earth: anelasticity, anisotropy, and spherical spreading being the most important. The wavelet construction begins by taking a so-called “Ricker” wavelet as a basis:

(33)

where is a time scaling factor that gives the exact distance to the first zero crossing. Figure 22 shows a Ricker wavelet for seconds. The basic Ricker wavelet is zero phase. The input cell allows one to vary the time scaling factor. The Ricker wavelet is non-causal (starts before time zero) and is therefore not possible in the real Earth.

Figure 22. Ricker wavelet. Note the zero-phase (symmetrical) character of the basic wavelet.

The spectrum of the Ricker wavelet is:

(34)

The Fourier transform functions of Mathematica readily show that equations (33) and (34) form a Fourier transform pair. To apply the Ricker wavelet to varying offset, I use equation (34) for the spectrum, add the phase shown by the upper-right graph in Figure 21, inverse transform to the time domain, and then apply the amplitude factor shown in the upper-left graph of Figure 21. Lastly, the time delay shown in the bottom-left graph of Figure 21 is applied to make the wavelet appear on the section at the correct arrival time. This is all done in the function plotarrivals. The section of time series for regular-spaced offsets formed in this way is shown in Figure 23. Trace amplitude and spacing can be varied in the input cell. The function plotarrivals has the following arguments:

1. model = the output of getmodel
2. system = 1 for SH or 2 for P-SV
3. wave = "PP" or "SS" or "PS" or "SP" for P-SV system or "SH" for SH system
4. component = "R" for radial or "Z" for vertical (applies only to system = 2); "T" is appropriate for system = 1
5. lastlayer = layer number at which wave will reflect off bottom
6. pvso = a list of lists, each level 1 element being the pair defining a ray from source to receiver
7. ascale = scaling factor on wavelet height
8. plotmode = "relative" for viewing amplitudes in correct relation to one another or “uniform” for constant amplitude
9. t0 = time scaling constant of Ricker wavelet (s)
10. twave = Ricker wavelet span (s) from negative to positive cutoff
11. dt = sampling interval for time series (s)
12. plotlabel = label for the plot
13. aspect = ratio of plot’s vertical axis (offset) over horizontal axis (time)

The next command plots a section of traces showing the relative amplitudes that would be recorded versus offset. (This function takes some time to evaluate.)

Figure 23. A section of waveforms that would be recorded for the raypaths shown in Figure 20 using the amplitude and phase information of equation (32) applied to a Ricker wavelet.

AVO in the P-SV Case

Consider now the P-SV wave system. A formulation similar to SS in the SH system works for PP, the reflection of a P wave into a P wave in the P-SV system, or for SS in the P-SV system. Again, all the necessary reflection and transmission coefficient expressions have been developed in a previous section. The AVO results for PP and SS results can be expressed as:

(35)

(36)

Often the results for asymmetrical raypaths consisting of downgoing P and upgoing S, and vice versa, are of interest. In one case the reflection at the specified boundary is PdSu and in the other it is SdPu; I use the terms PS and SP to designate these AVO types. The expressions developed for the P-SV case are used to get the PS and SP results. The geometry of the PS or SP raypaths is not as simple as for the symmetrical SS or PP cases because of different downward and upward paths. However, both parts of the raypath are still governed by a single, constant ; and so the raypath is easily traced for the asymmetrical cases. The AVO results for PS and SP can be expressed as:

(37)

(38)

The cells below do the calculations for the symmetrical PP and SS rays and also the asymmetrical PS and SP waves. The function that plots the AVO results for the P-SV system is plotPSVAVO, and its arguments are:

1. model parameters returned by getmodel
2. wave = "PP" or "SS" or "PS" or "SP" for P-SV system
3. deepest layer through which rays propagate (the reflection occurs at the bottom of this layer)
4. a label for the plot
5. the maximum offset horizontally, from source to receiver (km)
6. pvso, the list returned by rayparms

The rays and AVO results for PP are shown in Figures 24 and 25, and those for SS are shown in Figures 26 and 27, again for the Moho reflection. The same model as for the SH-wave case above is assumed. Because the ray geometry is only affected by the layer velocities and the compressional and shear velocities in each layer are in proportion as , the raypaths for PP and SS are identical.

This command computes and plots the rays. The input cell to plot the model (see above) must be evaluated first.

Figure 24. Rays for PP in the Earth model shown in Figure 19. Vertical exaggeration distorts the true angles.

The functions getmodel (above) and rayparms (above, with ) must be run before calling plotPSVAVO.

Figure 25. AVO calculations for PP in the Earth model shown in Figure 19.

This computes and plots the rays. The input cell to plot the model (see above) must be evaluated first.

Figure 26. Rays for SS in the Earth model shown in Figure 19. Vertical exaggeration distorts the true angles.

The functions getmodel (above) and rayparms (above, with ) must be run before calling plotPSVAVO.

Figure 27. AVO calculations for SS in the Earth model shown in Figure 19.

The asymmetrical raypaths are shown in Figure 28 for PS, again for the Moho reflection; the corresponding AVO results for PS are shown in Figure 29. Raypaths for SP are shown in Figure 30, and the corresponding AVO results in Figure 31. The AVO amplitudes for PS and SP in Figures 29 and 31 differ only by a constant, and the phases are identical.

Compute and plot the rays. The input cell to plot the model (see above) must be evaluated first.

Figure 28. Rays for PS in the Earth model shown in Figure 19. Vertical exaggeration distorts the true angles.

The functions getmodel (above) and rayparms (above, with ) must be run before calling plotPSVAVO.

Figure 29. AVO calculations for PS in the Earth model shown in Figure 19.

Compute and plot the rays. The input cell to plot the model (see above) must be evaluated first.

Figure 30. Rays for SP in the Earth model shown in Figure 19. Vertical exaggeration distorts the true angles.

The functions getmodel (above) and rayparms (above, with ) must be run before calling plotPSVAVO.

Figure 31. AVO calculations for SP in the Earth model shown in Figure 19.

Record sections for the PP, SS, PS, or SP arrivals can be created in the same manner as for the SS reflection in the SH-wave case, but are not shown here. The function plotarrivals will take into account the emergent direction of the P or SV waves and compute the motion on the radial (R) or tangential (T) component, as prescribed in the function argument for the component.

Summary

I have written the formulas for reflection and transmission coefficients for plane elastic waves scattered at a plane boundary between two half-spaces as Mathematica expressions. Plotting these coefficients gives results similar to those previously published as line graphs with angle of incidence as the independent variable. I have also written the formulas for vertical energy flux, based on these coefficients, as Mathematica expressions. Again, plots of the energy flux are similar to those previously published.

For the case of SH waves in elastic media, I have used Mathematica‘s algebraic manipulation functions to change the expressions into ones dependent on the ratio of the velocities of the two half-spaces and on the ratio of the densities. This simplifies the computations without loss of generality. For the SH case, I have investigated Brewster’s angle (angle at which there is no reflected wave). This has a single-valued curve versus angle of incidence of the SH wave when there is no density contrast between the half-spaces, but it becomes fairly complex for an interface with a density contrast. By contour-plotting the SH reflection-transmission energy fluxes in terms of the ratio of the bottom and top velocities and the incidence angle, one sees the complete behavior of the energy flux relations through a color scheme applied to the magnitude of the energy flux. Previously, the energy fluxes were only shown as line graphs for specified media parameters, due to computational limitations.

In the case of P-SV waves in elastic media, I have parameterized the reflection-transmission coefficients in terms of four dimensionless parameters: density ratio of the two half-spaces, velocity ratio of the two half-spaces, and compressional-to-shear velocity ratios of each half-space, thereby reducing the number of variables from six to four without loss of generality. The P-SV case is much more complex than the SH case. By contour-plotting the P-SV reflection-transmission energy fluxes in terms of the ratio of the bottom and top shear velocities and the incidence angle, one sees the complete behavior of the energy flux relations through a color scheme applied to the magnitude of the energy flux, as for SH waves. As for SH, the P-SV energy fluxes were only shown as line graphs for specified media parameters, due to computational limitations. The three-dimensional visualization allows for more complete understanding of the complex scattering behavior at an interface. The critical angle relations have a conspicuous role in the shape of the contours in these plots and usually bound a strong increase or decrease in energy flux. We also see the effect of Poisson’s ratio in the case of SV incidence. In all cases, the energy flux in the top layer and that in the bottom layer sum to unity, such that the plotted energy fluxes are reciprocal relations; this is also true when the energy fluxes are repartitioned into total P and total S.

I have also investigated the amplitudes of waves reflected from a deep reflector within a stack of layers as the distance between source and receiver along the surface varies. This involves applying the reflection/transmission coefficients to several layers. Results were interpolated from equal increments of ray parameter to equal increments of source-to-receiver distance (offset). Results were then applied to a Ricker wavelet and put in record sections showing the arrivals in time versus offset. Termed “amplitude versus offset” or “AVO” in the literature, I have shown that this behavior can be investigated efficiently and easily using Mathematica.

Conclusion

The scattering of plane waves at the boundary between two elastic half-spaces is a problem which can be readily formulated and presented in Mathematica. Mathematica handles all the algebraic manipulations, the calculations using these formulas, and also draws all the plots necessary to understand the results. No additional tools were required in composing this article.

These plots shown in this article provide a theoretical explanation of why small-angle (relative to vertical) reflections in layered media of small velocity contrasts give small reflected amplitudes for either P or S incident from either above or below, but then can give relatively large reflected amplitudes for incidence angles near grazing at a boundary. The phase of the scattered waves is as important a feature as the amplitude and should always be considered in the analysis of real data.

Acknowledgments

The author wishes to thank the many previous authors who have treated the complex interactions at the boundary between elastic half-spaces. Without tools such as Mathematica, these authors often laboriously performed the calculations by hand or with first-generation computers. The text of Keiiti Aki and Paul Richards [5] has been invaluable in allowing the author to evaluate and understand the complex interactions.

References

[1]	C. G. Knott, “Reflection and Refraction of Elastic Waves, with Seismological Applications,” Philosophical Magazine, 5th series, 48, 1899 pp. 64-97.
[2]	K. Zoeppritz, “Erdbebenwellen VII B, Über Reflexion and Durchgang seismischer Wellen durch Unstetigkeitsfläschen,” Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, 1, 1919 pp. 66-84. bssa.geoscienceworld.org/cgi/content/abstract/66/6/1881.
[3]	G. B. Young and L. W. Braille, “A Computer Program for the Application of Zoeppritz’s Amplitude Equations and Knott’s Energy Equations,” Bulletin of the Seismological Society of America, 66(6), 1976 pp.1881-1885. bssa.geoscienceworld.org/cgi/content/abstract/66/6/1881.
[4]	J. E. Nafe, “Reflection and Transmission Coefficients at a Solid-Solid Interface of High Velocity Contrast,” Bulletin of the Seismological Society of America, 47(3), 1957 pp. 205-219. bssa.geoscienceworld.org/content/vol47/issue3.
[5]	K. Aki and P. G. Richards, Quantitative Seismology, 2nd ed., paperback edition, Sausalito, CA: University Science Books, 2009.
[6]	E. D. Guy, S. J. Radzevicius, and J. P. Conroy, “Computer Programs for Application of Equations Describing Elastic and Electromagnetic Wave Scattering from Planar Interfaces,” Computers and Geosciences, 29(5), 2003 pp. 569-575. www.ohio.edu/people/guye/EDGuy_Papers/EDGuy-Computers-and-Geosciences.pdf.
[7]	W. C. Elmore and M. A. Heald, Physics of Waves, Mineola, NY: Dover Publications, 1985.
[8]	S. Chopra and J. P. Castagna, “Introduction to This Special Section—AVO,” The Leading Edge, 26(12), 2007 pp. 1506-1507. tle.geoscienceworld.org/cgi/content/abstract/26/12/1506.
[9]	K. Priestley and J. Brune, “Surface Waves and the Structure of the Great Basin of Nevada and Western Utah,” Journal of Geophysical Research, 83(B5), 1978 pp. 2265-2272. europa.agu.org/?view=article&uri=/journals/jb/JB083iB05p02265.xml.
D. H. von Seggern, “Exploring Reflection and Transmission Coefficients in Elastic Media,” The Mathematica Journal, 2012. dx.doi.org/doi:10.3888/tmj.14-2.

List of Additional Material

Additional electronic files:

1. AVOfunctions.m
2. RTfunctions.m
3. p-b.dat

About the Author

David von Seggern began his career in 1967 working in the field of detecting, locating, and characterizing underground nuclear explosions. In 1982 he joined Phillips Petroleum to apply seismology to oil and gas exploration. His interests there were in wave propagation and in 2D and 3D seismic imaging. He relocated to the University of Nevada, Reno in 1992 to become the manager of the Yucca Mountain Seismic Network, and retired from that role in 2005. Von Seggern continues as emeritus faculty at UNR.

David H. von Seggern
Nevada Seismological Laboratory, MS 174
University of Nevada
Reno, NV 89557
vonseg@seismo.unr.edu