Geol 456/656
Velocity Structure of the Earth

Spheres of the Earth

The phenomenon of plate tectonics gives evidence of sources of energy within the Earth. Although the Earth's heat output each day is much smaller than the amount of solar energy it receives, it is enough to build mountains that resist the solar-driven processes of erosion. To understand how plate tectonics works, we must understand the nature of the Earth's interior.
Question: what are the sources of heat energy within the earth?
earth layers diagram
(J. Louie) Above is a simplified view of the composition and state of regions within the Earth. Note how the regions are arranged in radial shells. Compositions and properties vary vetween shells a great deal, while relatively little within any one shell. While planetary physics, the properties of the geomagnetic field, Earth's gravity field, tides, geochemistry, and geology all contribute to knowledge of the Earth's interior, defining and characterizing the layers of the Earth has basically been the business of earthquake seismologists.

earth velocity plot
(from Kearey & Vine, copyright Blackwell Sci. Publ.) Timing of the arrivals of seismic waves is, of course, most sensitive to the velocity property of the materials encountered. While other properties, such as incompressibility, rigidity, and density may be inferred from velocities, only velocities can be measured accurately and directly.

Note how the transitions between layers show both gradients and discontinuities:

Earth Discontinuities
Depth, km	Name

5144 Lehman - Fe solid against FeO, FeS fluid 2885 Gutenberg - fluid FeO, FeS against (Mg, Fe) silicates, velocity decrease, density increase 2870 D'' - thin, mixing of mantle and core material? 670 "670 km" - worldwide, no earthquakes deeper, debates over whether a composition, phase, or viscosity change 400 "400 km" - worldwide, structure more variable above, phase change to spinels 50-200 LVZ - really a couple of gradients, regionally variable 4-55 Moho (mo-ho-RHO-vi-chich) - sharp compositional change to crust, tectonically active? 5-30 Conrad - mafic to felsic crust, often absent

prism photo (original image from the Exploratorium; used by permission) The above discontinuities were discovered by Lehman, Gutenberg, Mohorovicic, and Conrad by virtue of their ability to refract or bend seismic waves, just as a prism bends light waves.

Question: state Snell's Law of refraction in terms of refractive index, and in terms of velocity.

earth ray diagram
(J. Louie) Beno Guterberg located the core-mantle boundary from the P-wave shadow zone, where P-waves are bent away from the boundary, and from the larger S-wave shadow zone, showing that the core is fluid.

core diagram

(J. Louie) Just by knowing the delta degree angle of the onset of the S-wave shadow zone, and the radius of the Earth, you can make a simple estimate of the radius of the core.

Question: knowing the true core radius, how much do velocity changes in the mantle bend S waves up?

Seismic Refraction

According to Snell's Law, a seismic wave hitting any velocity discontinuity will partition into different wave phases (just as melting rock will partition into different chemical phases). Each wave phase has a characteristic pattern of arrival time versus distance from the earthquake, which we represent as an equation on a time-distance (t-x) plot:

constant velocity phase sketch
In a constant velocity medium, there are no refractions, and the t-x plot shows a straight line through the origin, with a slope of the inverse of the velocity.

discontinuity phase sketch
A discontinuity gives you the original direct phase as well as a refraction (an line offset from the origin) and a reflection (a hyperbola). There is also a wave transmitted through the discontinuity that appears elsewhere ... Pg P* Pn phase diagram
(from Kearey & Vine, copyright Blackwell Sci. Publ.) such as where a deeper discontinuity interacts with it further.

Question: How many reflections are not shown above?
gradient phase sketch
A simple velocity gradient produces a refraction that takes a circular path, having a hyperbolic sine (or Error Function) shape in the t-x plot.

gradient discontinuity phase sketch
A gradient discontinuity produces hybrid reflection/refractions, forming a triplication pattern in t-x, since some observation distances will note 3 arrivals from the discontinuity: a direct wave, a refraction, and a reflection.

The cut off, of the triplication's back branches at B and C, are useful in estimating the steepness of the gradient, which together with experimental data can determine whether the discontinuity is a compositional or a phase change.

M. Walck data and synthetics figure
In her thesis Maryanne Walck used these data and synthetics to control the steepness of the gradients at the 400 and 670 km discontinuities below the Gulf of California. The t-x diagrams have been skewed by a reduced time t - t * 10 (km/s), for presentation, so velocities less than 10 km/s tilt left, and more than 10 km/s tilt right. Walck estimated the 400 km discontinuity to be a shallow-gradient phase change, while the 670 km discontinuity had the steep-gradient character of a compositional change.

Question: Identify a 400 km back branch and a 670 km back branch in the data above.
crustal transition sketch
The above methods assume a purely radial or layered symmetry, with no lateral discontinuities. Repeated, localized experiments can show general changes, as in the thickness of the crust, as long as the experiments did not intersect any lateral transition.

Pn velocity map of US
(copyrighted figure) This map of velocity just below the Moho was contoured from many localized experiments around the 48 states. The higher velocities below the Sierras are an artifact of the experiments' spanning a lateral discontinuity.

Question: what is the significance of the low Pn velocities below the Great Basin?

Seismic Tomography

To properly measure lateral discontinuities, we need to make lots of measurements and use tomography.
Tomography is the reconstruction of an image from its projections, or shadows.
Doctors in the analog age used to x-ray a patient from several directions, arrange the films around a disk of ground glass, and shine lights through them from the outside to back-project the shadows into the glass.
Note: you won't need to know any of the fancy matrix manipulations below to pass this class. However, you will be tested later on vector algebra, after we go through it.
slowness block sketch A tomographic study starts by dividing the region (here a map) into blocks. Number each block with a block number b. Each block has a slowness s: s(b) = 1/v(b).

With many sources and many receivers, the number of rays will be (#sources)(#receivers) = nr . Each ray, numbered r, has a travel time t: t(r) = sum(b) l(rb) s(b). lrb is the length of ray r in block b.

Of course the length of any one ray in most blocks is zero. For a couple of rays, the time summation looks like:
t1 = l11s1+l12s2+...
We can write the summation above in matrix form:
(t1,t2,...)=(l11,l12,...;l21,l22,...;...)(s1,s2,...)
or t=Ls. The ray length matrix L has dimensions nr by nb, and is very sparse.

The matrix equation t=Ls is a statement of the forward problem, or how to model synthetic travel times when given a velocity model input. (Arbitrarily numbering our time data and slowness model values and packing them into a vector takes some getting used to, but makes the math a lot simpler.) In geophysics we are given data input as a set of travel times, and we wish to solve for an estimated velocity model.

Suppose we could somehow invert L:
s = L^(-1) t We can likely at least estimate the ray paths, so we have no problem finding the ray length matrix L. Then we could determine the slowness model from the time data.

Question: what's in the Identity matrix I?
Since our number of observations never matches the number of blocks, and L is sparse, we use least squares by pre-multiplying the forward equation by L-transpose to form the normal equations
L^T t = L^T L s
Then the problem reduces to finding the inverse of the square matrix L^T L to get the slowness solution
(L^T L)^-1 L^T t = s
Question: why are these called the normal equations?
However, a typical tomography problem will be characterizing a region 100 by 100 blocks in size, or having 10,000 total blocks. This means that L^T L will be a 10,000 by 10,000 matrix, with 100 million elements. Unless you have a very large supercomputer at your disposal, it is extremely difficult to invert such a large matrix with gaussian elimination or singular-value decomposition.
Question: even if your computer has a large amount of virtual memory, it won't help you here. Why not?
The secret of tomography is to make the initial assumption that only the diagonal of L^T L matters, and approximate it with a diagonal matrix D, which is trivial to invert:
L^T L ~ D , D(b)=sum(r)(l(rb))^2
Inverting D instead of L^T L is a first-order linear approximation to the full least-squares inverse.
Question: what in general would be a way to improve upon a first-order solution? Think, for example about the steps in a long division.
Substituting D for L^T L in the least-squares slowness solution above, and re-formulating as summations, we get the tomographic solution:
s(b) = sum(r)t(r)l(rb)/sum(r)l(rb)^2 ray blocks sketch
The solution shows that we are just back-projecting each ray onto the model, one at a time, and summing the slownesses into each block. With a total ray length L, we are giving each block hit the average slowness of the ray t/L. Each block accumulates the slownesses of all rays that hit it, and scales them at the end.

We can solve large tomography problems on any computer, since we only need two bins of storage for each block (for the 100x100 block model, maybe just 80 kbytes) and we process the rays one at a time. We never have to keep a huge L^T L matrix in memory.

Tomographic solutions do have some serious problems. Most important is a tendency for anomalies to be drawn out in the direction of most of the rays passing through a region. In geophysics it is very difficult to achieve an even and isotropic (omni-directional) ray coverage.

mantle velocity section
Above is a tomographic inversion result for a section of the Earth's mantle, at the equator (Original 0.3 Mbyte PostScript file from Harvard, used by permission). Cool colors represent represent positive deviations of velocity from the radially-symmetric average at that depth, and warm colors represent negative velocity deviations. Since the mantle has a nearly constant composition, velocity deviations are thought to be due to differences in temperature, with cool colors for fast cold mantle and warm colors for slow warm mantle. The dashed circle is the 670 km discontinuity; plate boundaries are yellow on the index map at center.

Warm, viscous parts of the mantle are less dense than their surroundings and will rise buoyantly over geologic time (as fast as your fingernails grow). Relatively cold mantle will sink. Driven by heat escaping the core, and cooled in the oceans, the mantle will circulate or convect, just like a boiling pot. A complete turnover takes hundreds of millions of years. The interiors of all planet-sized bodies must be actively convecting, to release their heat of formation. Any planet with a radius over ~1500 km cannot conduct its internal heat away within the age of the universe, so it must convect viscously to release its heat, or it would melt and then convect as a fluid.

3-d mantle plume image
(3.3 Mbyte printable PostScript file from Harvard, used by permission) A 3-d view of a mercator projection of the mantle, with orange surfaces surrounding warm blobs of mantle, which should be rising plumes.

Question: what areas are above rising mantle plumes?
3-d mantle slab image
(3.3 Mbyte printable PostScript file from Harvard, used by permission) A 3-d view of a mercator projection of the mantle, with blue surfaces surrounding cold blobs of mantle, which should be sinking slabs.
Question: what areas are producing sinking mantle slabs?
Question: do you see any evidence in the above results for or against two layers of convection in the mantle?

Previous: Seismology Fundamentals -- Up to Syllabus -- Next: Earth Compositions