View License

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Highlights from
Okada: Surface deformation due to a finite rectangular source

Join the 15-year community celebration.

Play games and win prizes!

» Learn more

5.0
5.0 | 8 ratings Rate this file 21 Downloads (last 30 days) File Size: 8.25 KB File ID: #25982 Version: 1.7
image thumbnail

Okada: Surface deformation due to a finite rectangular source

by

 

30 Nov 2009 (Updated )

Computes Okada's 1985 solution for displacements, tilts and strains due to fault dislocation.

| Watch this File

File Information
Description

The Okada [1985] model calculates analytical solution for surface deformation due to shear and tensile faults in an elastic half-space. This model is widely used to simulate ground deformation produced by local perturbation like tectonic faults (earthquakes) or volcanic dykes (magmatic intrusion). Given rectangular fault geometry (length, width, depth, strike, dip) and 3-component dislocation amplitude (rake, slip and open), it computes the displacements, tilts and strains at the free-surface.
The proposed Matlab script is a literal transcription of the Okada's equations, except that it is transposed in a geographical referential (East, North, Up), where the fault is defined by a strike angle relative to the North, and dislocation parameters are given by: rake, slip and opening (instead of U1, U2, U3), following Aki & Richards [1980] definition. All coordinates and depth are relative to fault centroid. Lamé's constants λ and μ are replaced by Poisson's ratio ν (with a default value of 0.25 for isotropic medium), since the equations are independent of other elastic parameters. The equations are also vectorized for (x,y) coordinates and all input parameters except dip angle.

To check the consistency of numerical calculations, run the script okada85_checklist.m, a transcription of table 2 cases 2, 3, and 4 checklist from [Okada, 1985] paper (needs also the roundsd.m function).

See help for further details, syntax, example, and script comments for technical details.

MATLAB release MATLAB 6.5 (R13)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (11)
26 Oct 2016 Raquel Felipe

no entiendo por que me marca un error cuando hago [UE, ONU, uZ] = okada85 (E, N,14000,193,24,600000,200000,81,11.0586, 0);

Comment only
12 Jul 2016 FK

FK (view profile)

a

Comment only
03 Jun 2016 he

he (view profile)

After making comparison with Okada's code revised in 2002 by himself, I find that this MATLAB version works well as Okada's ! Thank you!

03 Jun 2016 he

he (view profile)

Great code! For deformation on surface, there are 9 independent variables to give for users. But maybe one should note that: uZN=-uNZ, uZE=-uEZ, (uNN+uEE)/u(ZZ)=(1-nu)/nu on free surface without stress. And these equations can be easily given by constitutive relation on free surface.

(nu is Poisson's ratio, uZZ is strain in z direction)

20 Apr 2016 Noverina Alfiany  
19 Dec 2014 Ding

Ding (view profile)

 
10 Dec 2013 chen

chen (view profile)

 
27 Mar 2013 hui zhou

Hello,I'am puzzled about how to get E,N when the code is used for InSAR measurements:
[uE,uN,uZ] = okada85(E,N,14000,193,24,600000,200000,81,11.0586,0) ;
Thank you

Comment only
17 Aug 2012 Matang

Matang (view profile)

Great code!
Try this:

Case of the 2011, Mw 9.0 Tohoku earthquake

Slip (D):
Mo = µ D S
µ = 30 GPa (elastic shear modulus)
S = 600000 m * 200000 m (fault surface)
log10(M0) = 1.5 * Mw + 9.1
Mw = 9.0
obtained: D = 11.0586 m

%----------------------
%Epicenter
LatEpic=38.297;
LonEpic=142.372;
% Observation coordinate at GEONET GPS station (#0175)
LatGPS_sta=38.68269882693;
LonGPS_sta=141.4493698139;

E=(LonGPS_sta-LonEpic)*111*1000;
N=(LatEpic-LatGPS_sta)*111*1000;

[uE,uN,uZ] = okada85(E,N,14000,193,24,600000,200000,81,11.0586,0)
%----------------------

got the the solution:
uE = 4.0607
uN = -1.8354
uZ = -0.4523

The GPS data from the ARIA team at JPL and Caltech (version 2.0) for the station:
uE = 4.0394
uN = -1.612
uZ = -0.6559

--Muzli
GFZ-Potsdam

26 Jun 2011 reza zohd  
05 Oct 2010 Forrest Brett

very clean coding.

Updates
19 Mar 2010 1.1

Corrects some output arguments description in the help (inversion of uN and uE), and better explanation of input arguments.
Adds a 3-D graphical output the of fault geometry when using without output argument.

24 Sep 2010 1.2

ATTENTION: correction of 3 errors in some equations (functions I1, K2 and uyy_tf) affecting some component values. Now the function fits all the numerical values of Okada's checklist cases 2, 3 and 4. Thanks to Dmitry Nikolski for his contribution.

05 Dec 2010 1.3

Coordinates are now relative to fault centroid (instead of middle top edge). This is more convenient for earthquake hypocenter use.

05 Mar 2011 1.4

Bug correction on fault centroid exact position. Add optional plot and example.

06 Mar 2011 1.5

Corrects another bug in the plot figure.
Minor modifications in the code and help text.

30 Aug 2012 1.6

- allows vectorization of RAKE, SLIP and OPEN
- checklist function to validate numerical consistancy

28 May 2014 1.7

fixes a problem with K1 function (tilt) when DIP=90. Thanks to Halldor Geirsson.

Contact us