Seismic tomography: the easy way

When I was an undergrad, I had to take a year of geophysics, and on the first day, the professor explained that geophysics meant finding things underground without digging.  The rest of my lab group here at the UW are all marine geophysicists, and they routinely have to find things underground without digging, but with the added complication of having to get through a couple of kilometers of water, too.

You may remember my interview with Rob Weekly a few weeks back. This week’s interview is with my office-mate, Dax Soule. Dax studies marine geology and geophysics in the School of Oceanography at the University of Washington. Like Rob, he focuses on active source seismology. In fact, he’s looking at some of the same data as Rob, just different parts. And since we already went over a lot of the basic background in Rob’s post, I’m going to jump right into some fun tomography stuff. In a nutshell, tomography means building up a 3D picture of what’s underground by measuring how long it takes for seismic waves to travel through. In tomography, you’re not sampling the material directly, but you’re measuring how the seismic velocity varies. Seismic velocity is a measure of how quickly a seismic wave can travel through a certain type of rock, so if you know the velocity structure you can make inferences about the geological structure.

Imaging the earth’s crust

When you’re doing seismic tomography, it’s a bit like peeling back the layers of an onion: you have to figure out the shallow stuff before the deep. Rob’s work focused on obtaining a three-dimensional seismic velocity structure from the seafloor down to about 2.5-3 km depth. Dax is taking Rob’s results and extending downward to include the entire depth of the crust. He’s looking between the seafloor and the Mohorovicic discontinuity, or the “Moho”, which is the boundary between the earth’s crust and the mantle.

In seafloor maps of the region, you can see that there’s a bathymetric high – a plateau where the water depth is shallower than the surrounding seafloor. A seismic reflection study a few years back [1] showed that this plateau likely corresponds to a thickening of the oceanic crust at this location. Dax’s tomography work will help to clarify exactly what’s going on, and will give a more detailed picture of how the crustal thickness varies in the area.

Dax explained that measuring crustal thickness variations near a mid-ocean spreading center can tell you about how that crust was produced – was a episodic or constant? And where was the source of new crust – was it on the axis of the ridge, or off-axis?

How the experiment works

Here’s a little cartoon showing the basic geometry of the seismic experiment.SeismicRays_600px

A seismic source is generated near the sea surface, and the energy travels through the water column and into the sea floor. Once in the earth’s crust, the energy is converted into to types of waves: primary waves (p-waves) and secondary waves (s-waves). There are many different types of paths that these waves can take as they travel through the earth, and each type of path is known as a “phase”. In the above figure, the two phases that are shown are Pg and PmP. Rob used the Pg phase arrivals to image the upper portion of the crust. Dax looks at the PmP phases, which are the ones that bounce off the Moho.

Dax goes through the data and manually picks out the PmP phase arrival times. This part of the job is not the most exciting, but someone’s got to do it! Fortunately, all that work does pay off, and once Dax has his picks, he can start digging into the tomography part of his work.

Five easy steps
Here’s a handy-dandy summary of the steps that Dax goes through to build up a tomographic inversion.


So that’s what my office-mate does… good to know!  Dax is still working through the data and the very complex inversion code, so stay tuned for a future post on what he finds deep in the crust, and what the crustal thickness can tell us about seafloor production at the Endeavour Ridge.

[1] Carbotte, Suzanne M., Mladen R. Nedimović, Juan Pablo Canales, Graham M. Kent, Alistair J. Harding, and Milena Marjanović. “Variable crustal structure along the Juan de Fuca Ridge: Influence of on‐axis hot spots and absolute plate motions.” Geochemistry, Geophysics, Geosystems 9, no. 8 (2008).

Interpolation using Scipy/Numpy

As part of a little script I’m writing, I need to do some simple linear interpolation. The Matlab equivalent of what I’m doing is called ‘interp1’. So far I’ve come across two ways to do 1-d, linear interpolation. One is ‘interp’ in numpy. The other is ‘interp1d’ in scipy.interpolate. I’m not completely sure of the difference. Is one faster? Is one older? I’ve googled around a bit, but still haven’t found a clear answer. For now, I’ve implemented ‘interp1d’, which is less similar to Matlab’s ‘interp1’. You first define an interpolated object given your x and y vectors. Then you call that object with your new x-values to generate the new y-values.

And here is an example, chopped out of my code:

[sourcecode language=”python”]
depthvector1 = r_[0:nadirdepth:DepthIncrement]
interpfun1 = interp1d(depth1,ssp1) # my x and y vectors
sspvector1 = interpfun1(depthvector1) # find the new y’s for this x

Python – matrices or arrays?

Sometimes I think I have my biggest breakthroughs while biking home. Tonight, for instance, it dawned on me (on the stretch between the beach and Patterson) just how awful the Matlab version of my ray tracing code really is. And I realised that I could vastly improve its simplicity and efficiency by using matrices (Yikes! Why on earth did I implement all of those nested loops?). I pictured the steps and the math involved, the Matlab commands that I would need, and got really excited – until I realized that I don’t know how to do any of the same stuff in Python yet. First of all, I need to figure out how to work with arrays in Matlab. Thankfully, there’s a special page on Matlab users who are trying to migrate to Numpy/Scipy. Yes! And here it is: Numpy for Matlab Users. My first question was, “How do I make it do element-wise operations?” In Matlab, you need to modify the operator, for example using “.*” to multiply by elements, and just “*” to do matrix multiplication. Matlab thinks in linear algebra for 2d matrices. But, as I quickly learned, the default array in Numpy does element by element operations, unless you create a special matrix object. Awesome.

Tkinter and GUI building

I’m back to ray tracing again. And since I’m trying to re-create and expand on the Matlab version of this code, I am finding myself having to learn lots of little details – and big details. One of these things is GUI building. The standard GUI builder that comes with Python is TKinter. I’ve used this a bit before, but now I’m having to think about doing more complicated things with it. The Matlab version of my code used an interactive GUI. In Matlab, there’s a utility called GUIDE that gives you a GUI to build a GUI 🙂 You can put buttons, sliders, graphs, etc, into your GUI window, and then it will automatically generate a skeleton script that you can fill with commands telling Matlab what to do with it. So now I’m wondering if there’s something I can use with TKinter that does a similar thing. I’m sure there is, and I just need to find it…

Python for ray tracing

I just found some python scripts that I started writing about a year ago, including one in which I was trying to read in a sound speed profile through the water column, and using Snell’s law to compute a *very* simple piecewise ray trace.  So I create several beams that launch from one location (ie. a multibeam head), all departing a different launch angle between -75:75 degrees. The purpose of it is to see how much effect an incorrect sound speed profile will have on the estimate of the seafloor depth for each beam.  I already wrote it in Matlab, so re-writing it in Python was really just an exercise to learn Python.  I’d only ever gotten as far as asking the user which sound speed files to use.  I think it’s a good project for me to work through. 🙂