Revision as of 00:34, 2 October 2008

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 layer depths are 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 calculates 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. **The layered model scripts require the length scale to be km.**

First make a vector with the depths to the bottom of each layer. Create a model vector, calculate the elastic parameters, and choose a scale term. **Note the scaleN is not required the default is 1.** The elastic parameter inputs need to be relative to the first layer, so divide by the first value in each vector. Notice that the program crashes, even though all of the inputs seem to be valid. You can get around this problem by adjusting the length scale, instead of using 1km as the length unit, use 10km. Now everything works just fine! Notice that the fault slip/opening does not have to be rescaled. Just remember that the units of displacement are the same as the units of slip/opening, whatever unit you choose.

Choosing a proper scaleN value is important, particularly if you are going to make many function calls to the layer model scripts. Use get_scaleN.m and get_scaleN_mogi.m to find the best scale value. The scaleN scripts require a model vector, a maximum distance, and a scale value.

   >>m = [3 3 5 -45];   % For the Okada model only specify the length, width depth and dip
   >>rmax = 50;
   >>scaleN = 0.3;
   >> [Uprop,Uokada]=get_scaleN(m,rmax,scaleN);

The figure is created automatically, it is useful for double checking that the geometry of your model is what you think it is, and to be sure that the layered model solution matches the homogeneous solution. The best approach is to make scaleN as small as possible while still matching the analytical solution to within some tolerance, say 0.1mm or 0.01mm. Remember that the accuracy of the result, with a given scaleN value, may depend on the depth and dip of the dislocation, so if you are inverting for source geometry you should test several model geometries to ensure that you are picking a scale value that gives accurate results for all conditions.

Simulated Annealing

There are a handful of Matlab functions that are designed for the easy implementation of the simulated annealing algorithm (Cervelli et al., 2001).