File Exchange

image thumbnail

Okada: Surface deformation due to a finite rectangular source

version 1.7.0.0 (8.25 KB) by François Beauducel
Computes Okada's 1985 solution for displacements, tilts and strains due to fault dislocation.

14 Downloads

Updated 28 May 2014

View Version History

View License

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.

Cite As

François Beauducel (2020). Okada: Surface deformation due to a finite rectangular source (https://www.mathworks.com/matlabcentral/fileexchange/25982-okada-surface-deformation-due-to-a-finite-rectangular-source), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (24)

he

Whyjay Zheng

Been using this code for years and it is very reliable.

Wthee

Would it be that difficult to add a quiver plot, like in this plot from Segall, https://ibb.co/V2qhKLK

Shamim Aziz

Great code. Thanks

François Beauducel

Dear Shamim Aziz,

The function works in relative plan geometry only, not in lat/lon. You must define a rectangular area in m or km around the fault rupture, for example :

[E,N] = meshgrid(linspace(-1000,1000,500)); % this defines a square area grid of 2 x 2 km around the fault with 2000/500 = 4-m step.
[uE,uN,uZ] = okada85(E,N,13.9,261,60,69,16,215,3.96,3.96,'plot'); % if you just need to plot displacements no need to output other parameters
figure, surf(E,N,uN)

Shamim Aziz

I want to plot the fault displacement plot using the fault parameter as below :
Longitude= 134.4138
Latitude =35.7569,
Rupture length =69
Rupture width =16,
Slip lenth=3.96
strike angle=261
dip angle=60
rake angle=215
fault depth=13.9(top=15,bottom=1.1)
can anyone please help me .
i have tried with

[E,N] = meshgrid(linspace(35.7569,134.4138,15));
[uE,uN,uZ,uZE,uZN,uNN,uNE,uEN,uEE] = okada85(E,N,13.9,261,60,69,16,215,3.96,3.96,'plot');
figure, surf(E,N,uN)

ben_jiang

Amri Rasyidi

can someone help me ?
i tried input as the example in the script, but i have error

>> okada85
Error: File: okada85.m Line: 465 Column: 1
Unexpected MATLAB operator.

Seok Goo Song

I want to compute the 6 independent stress components using the okada85, but a bit confused because of the tilt components in radian and compression is positive, and so on. Does any one give me a clear guideline for the calculation of the 6 stress comps?

Thanks!

Elizabeth Hearn

Ardiansyah Fauzi

cheng zeng

Hello,great code for deformation on surface.And in my case the rupture is on a finite fault which has some rectangular subfaults, The rake and slippage for each subfault is different.I wonder how to deal with this problem in my case, To get the displacement on the whole fault through which way? Please advise me on this.very thank you.

run li

Great code! thank

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);

FK

a

he

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!

he

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)

Noverina Alfiany

Ding

chen

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

Matang

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

reza zohd

Forrest Brett

very clean coding.

MATLAB Release Compatibility
Created with R13
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!