Matlab Users
Layered elastic Models
Email from Tom Fournier (June 17, 2008)
Jeff recently got some layered half-space codes for an isotropic point source and a planer dislocation that run in matlab. I have been playing around with them, so if you are interested here are some details that might make your life easier. The functions are very easy to use, so really hats off to Kaj Johnson who developed them.
There are 4 m-files (that I assume are in /user/local/matlab, if not I'll put them there)
LayeredGreens.m
MogiLayers.m
get_scaleN.m
get_scaleN_mogi.m
The real work horses are LayeredGreens and MogiLayers. The inputs for these two functions are the same except for the model parameter vector. The functions require as inputs; a model vector, station locations, depth of layers, relative Lame constants for each layer, and an optional scale ( default scaleN=1). The output of these functions is a familiar 3xN matrix of site displacments.
I wrote a simple script that calculates the Lame constants from a velocity profile, vel2elastic, (now in usr/local/matlab). This takes as inputs the p-wave velocity, density, poisson's ratio, and optionally the p-wave/s-wave ratio. The density can be constant or vary by layer. Poisson's ratio is assumed to be constant unless you give the p/s ratio, then poisson's ratio is ignored and the elastic parameters are calculated using Vs. The outputs are Lame's parameters; Lame1 (lambda) and shear modulus.
ScaleN controls how many terms are used in the Hankel transform, which you can determine using get_scaleN and get_scaleN_mogi. These two functions simply compare the layered half-space solution to the homogeneous half-space solution. The output are the deformations from each case and a figure ploting the displacements for comparison. These functions call mogi.m and disloc3d.m, so make sure you have a path to a working version of these defined. The inputs are; a model vector, a maximim station distance, and scaleN. The maximum distance is so that you can compare results at the distance of your farthest station (or the farthest you want the solution to be accurate). ScaleN controls the accuracy of the solution but also the time it takes to calculate it, so depending on your problem, you want to pick scaleN large enough to give good results, but also small enough that you don't have to wait forever to get some answers. Ideally an optimization function should be written that will return a scale value optimized for time with a given tolerance between the homogeneous and layered half-spaces solutions. This type of function wouldn't be hard to write, nevertheless it does not exist in our lab.
A couple of caveats about get_scaleN functions. For a given scaleN, the accuracy of the layered model is dependent on depth, dip and other parameters that affect the depth to any portion of the fault. Keep this in mind when choosing a scale value.
The other warning I should bring up is in regard to the LayeredGreens.m function. This function will crash in cases where the depth of the plane is too deep. I simply changed the length scales from 1km to 10km, it worked sometimes. If you are interested in this problem, particularly if you want to try and fix it, let me know and I can tell you more about it.
Example: Creating a layered model
This example uses the velocity model at Veniaminof volcano [J.J. Sanchez, 2005, UAF thesis].
The information from this table is imported into MATLAB and converted into an elastic model using VEL2ELASTIC.m.
The variable "data" is the same as the table.
VEL2ELASTIC takes as inputs the p-wave velocity (km/s), density in each layer(kg/m3), and Poisson's ratio. Alternatively VEL2ELASTIC can also use the Vp/Vs (p-wave/s-wave velocity ratio), in which case it ignores the value of Poisson's ratio. The output are the Lame parameters for each layer.
The scripts that calculate deformation in a layered half-space, LayeredGreens.m (Okada dislocation), and MogiLayers (isotropic point source), use the Lame constants and the depth to the bottom of each layer. The bottom layer extends to infinity, so the depth vector has a length n-1, while the Lame parameter vectors have a length of n, for a model with n layers.
