Geol 455/655 - Lab 4: Huygens, Fermat, Snell

Instructions: Turn in your answers to the questions in italics below.

1. We examine the input file carefully to understand the geological model we are computing waves through.
   1. Download and examine the e3d input file test4.in. Note that this is a two-dimensional grid that is a cross section, not a map. (Up to here we were looking at the surfaces of 3-d grids in map view.) The block parameter lines set the velocity model; we won't need ModelAssembler for this lab. The model in section is a simple crust over a mantle. How deep in kilometers is the top of the mantle? What is the total amount of time that will be computed?
   2. Run e3d with test4.in and observe the wave propagation in cross section. What is the Courant condition? Given the source frequency, how many grid points per wavelength are there, at least? Do you think there will be significant grid dispersion artifacts?

   Huygen's Principle

   Huygen's Principle states that each point on a propagating wave front acts as a new source of a spherical wave. (In this 2-d simulation, we have circular wavefronts.)

   

   1. Note that test4.in has seven additional source lines that are commented out. Remove the ``#'' comment characters and run e3d with this input. The visual starts out like the picture above. After a while the wavefronts from the 8 sources should coalesce into a fairly straight wavefront, moving straight down (and then reflecting off the crust-mantle boundary and moving straight up- the Moho reflects 10-20% of the incident amplitude, but the surface reflects 100% of the incident amplitude). The straight wavefront is called a plane wave (again, a line wavefront in 2-d).
   2. Construct an input file like test4.in, but one that will use 8 or more sources to emulate a plane wave that will propagate at 15 degrees away from vertical (left or right, your choice). Test the input and observe the computation. Turn in your altered input file, by email, so I can test it.
   3. Construct an input file like test4.in, but one that will use 8 or more sources to emulate a plane wave that will propagate horizontally from left to right across the section. Test the input and observe the computation. The plane wave breaks up. Describe why, or what about the geological situation descibed in the input file that causes the plane wave to separate into two (or more) plane waves. Turn in your altered input file, by email, so I can test it.

   Fermat's Principle

   Your run for 2.c. above is a good demonstration of what is called Fermat's Principle, which states that a wave will take the path of least time between its source and the observation point. This may not be the path that has the shortest distance. If it can, and still obey Fermat's Principle, a wave will take a longer, out-of-the-way route that finds the fastest velocities around.

   

   1. Given the length of a path d, and a constant seismic velocity along the path v, derive the wave travel time t. Turn in your equation for t, with a note on the units of each variable.
   2. Get the original test4.in back, with all the extra Huygen's sources commented out. Retain the one source on the left, and uncomment the third block line labeled ``batholith''. Considering that a true fluid has no rigidity and thus zero shear-wave velocity, is this model ``batholith'' entirely molten, or only partially?
   3. Run e3d with the batholith model and observe the wave propagation around the anomaly. Note that waves first get into the far end of the ``batholith'' from the far end, away from the source! This kind of wave propagation near strong lateral velocity anomalies can only be modeled with a finite-difference method such as e3d. Do waves bend toward or away from a low-velocity anomaly?

   Snell's Law

   The bending of waves to obey Fermat's Principle is described mathematically by Snell's Law, for a wave at an interface between two media of different velocity:

   sin(i₁)/v₁ = sin(i₂)/v₂
   where i₁ is the angle (from the vertical) of the propagation direction of the incident wave, in medium one with velocity v₁, and i₂ is the angle (from the vertical) of propagation of the bent wave, in medium two with velocity v₂.

   Note that these angles of wave propagation direction are perpendicular to the wavefronts. The path of propagation is called a ray. You can see the wavefronts easily on the sections below, but not the rays, so I've drawn a few in for you. There is a separate ray to every little Huygens source on a wavefront. Each of those rays started at the source at a slightly different angle of incidence, and then bends according to Snell's Law as it propagates. (The incident angles are really from the normal to the interface. Since we are talking here only about flat-laying interfaces, the normals are always vertical.)

   If velocity increases at depth, as it does at the Moho in our models, you can compute from Snell's Law that the ray's angle from vertical must also increase. So rays flatten and shallow out as they propagate through the Moho, for instance. You can see from the rays annotated below that this flattening is necessary for the wave to find the high-velocity mantle, shortening the travel time, and then curve back toward the crust.

   
   
   
   1. From top to bottom, the snapshots of wave propagation above are at times t=200, t=400, and t=700. At what time, in seconds, is each picture?
   2. How much do the rays bend if the angle of incidence is zero?
   3. Think for a moment about the bending at the Moho, where velocity increases. As you increase the incidence angle (say from the yellow P ray to the green Pn ray) eventually i₂ becomes 90 degrees. Since the sine of 90 degrees is one, the incident angle i₁ at which this occurs is called the critical angle i_c, and sin(i_c) = v₁/v₂ What is the critical angle for P waves i_c at the Moho, with the model in test4.in. Does this match what I've drawn for Pn above? Note that the Pn wavefront visible in the lower picture should be at angle i_c from the horizontal.
   4. Explain what would happen to the critical angle if velocity were to decrease with depth: with v₁ > v₂ . Is there still a refracted wave Pn? Convince yourself by altering test4.in to have lower velocities in the mantle, and observing a run.
   5. Snell's Law holds for reflections too. Of course, for the purple PmP reflection off the Moho shown above, both the incident and reflected wave are in the same medium with the same velocity. So i₁=i_r because: sin(i₁)/v₁ = sin(i_r)/v₁ So incident and reflected waves are at the same angle. Not so fast! Often, a wave will hit a boundary and reflect both a P wave and an S wave. The S wave is traveling at the shear velocity, slower than the P wave. Use Snell's Law to figure the reflected angle off the Moho of a converted PmS wave, from a P wave incident at 60 degrees. Run the original model and look for this converted PmS reflection in the visual. You can enhance it by raising, or drastically lowering, the shear velocity below the Moho.

Resources:

• E3D modeling code manual, defining inputs and outputs
• test4.in E3D input file, defining how E3D will run.