forked from IPGP/deformation-lib
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathokubo92.m
276 lines (244 loc) · 10.2 KB
/
okubo92.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
function varargout=okubo92(varargin)
%OKUBO92 Surface gravity change due to a finite rectangular source.
% [dG,dH] = OKUBO92(E,N,DEPTH,STRIKE,DIP,LENGTH,WIDTH,RAKE,SLIP,OPEN,RHO,RHOP)
% computes total gravity change at the free surface of an elastic
% half-space, due to a dislocation defined by RAKE, SLIP, and OPEN on a
% rectangular fault defined by orientation STRIKE and DIP, and size LENGTH
% and WIDTH. The fault centroid is located (0,0,-DEPTH).
%
% E,N : vector or matrix of coordinates of observation points in a
% geographic referential (East,North,Up) relative to fault
% centroid (units are described below)
% DEPTH : depth of the fault centroid (DEPTH > 0)
% STRIKE : fault trace direction (0 to 360° relative to North), defined
% so that the fault dips to the right side of the trace
% DIP : angle between the fault and a horizontal plane (0 to 90°)
% LENGTH : fault length in the STRIKE direction (LENGTH > 0)
% WIDTH : fault width in the DIP direction (WIDTH > 0)
% RAKE : direction the hanging wall moves during rupture, measured
% relative to the fault STRIKE (-180 to 180°).
% SLIP : dislocation in RAKE direction (in m)
% OPEN : dislocation in tensile component (in m)
% RHO : density of the medium (in kg/m^3)
% RHOP : optional density of the cavity-filling matter (in kg/m^3);
% takes RHO value if not specified.
%
% returns the following variables (same matrix size as E and N):
% dG : total gravity change (in m/s^2)
% dH : elevation changes as given by Okada [1985] model (in m)
%
% Length unit consistency: E, N, DEPTH, LENGTH, and WIDTH must have the
% same unit (e.g. km) which can be different from that of SLIP and OPEN
% (e.g. m to to produce dG in m/s^2 and dH in m).
%
% [...] = OKUBO92(...,RHOP,BETA,NU) specifies gravity vertical gradient
% BETA (default is 0.309e-5 s^-2) and Poisson's ratio NU (default is 0.25
% for an isotropic medium).
%
% Formulas and notations from Okubo [1992] solution excepted for the fault
% geometry parameters after Aki & Richards [1980], e.g.:
% DIP=90, RAKE=0 : left lateral (senestral) strike slip
% DIP=90, RAKE=180 : right lateral (dextral) strike slip
% DIP=70, RAKE=90 : reverse fault
% DIP=70, RAKE=-90 : normal fault
%
% Equations are all vectorized excepted for argument DIP which must be
% a scalar; all other arguments can be scalar or matrix of the same size.
%
% Example:
%
% [E,N] = meshgrid(linspace(-10,10,50));
% dG = okubo92(E,N,6,90,90,10,10,0,5,0,2670);
% figure, pcolor(E,N,1e8*dG), shading interp
% hold on, [c,h]=contour(E,N,1e8*dG,-50:10:50,'k'); clabel(c,h), hold off
% colorbar, polarmap
%
% reproduces the figure 4a in Okubo [1992] paper, considering a 10x10km
% fault at 6km depth, oriented N90°-strike, 90°-dip, and 5m senestral
% dislocation. Gravity changes are computed on a regular grid from -10 to
% 10km, and are plotted as a surface in 1e-8 m/s^2 unit. POLARMAP is an
% author's colormap.
%
% See also OKADA85 author's function to compute complete displacements,
% tilt and strain at free surface.
%
% Author: François Beauducel <[email protected]>
% Institut de Physique du Globe de Paris
% Created: 2012-06-11
% Updated: 2013-07-26
%
% References:
% Aki K., and P. G. Richards, Quantitative seismology, Freemann & Co,
% New York, 1980.
% Okubo S., Gravity and Potential Changes due to Shear and Tensile Faults
% in a Half-Space, J. Geophys. Res., 97:B5, 7137-7144, 1992.
%
% Acknowledgments: Bernard F. Whiting, Eric Chassande-Mottin, Martin
% Fuchs
% Copyright (c) 2013, François Beauducel, covered by BSD License.
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
if nargin < 11 || nargin > 14
error('Wrong number of input arguments.')
end
for i = 1:nargin
if ~isnumeric(varargin{i})
error('Input arguments E,N,DEPTH,STRIKE,DIP,LENGTH,WIDTH,RAKE,SLIP,OPEN,RHO,RHOP must be all numeric.')
end
end
if numel(varargin{5}) ~= 1
error('DIP argument must be scalar.')
end
% Newton's gravitational constant
G = 6.67384e-11; % (m^3/kg/s^2)
% Default values for optional input arguments
nu = 0.25; % isotropic Poisson's ratio
beta = 0.309e-5; % free-air gravity gradient (m.s^-2/m)
% Assigns input arguments
e = varargin{1};
n = varargin{2};
depth = varargin{3};
strike = varargin{4};
dip = varargin{5}; % 'delta' in Okubo's equations
L = varargin{6};
W = varargin{7};
rake = varargin{8};
slip = varargin{9};
U3 = varargin{10};
rho = varargin{11};
if nargin > 11
rhop = varargin{12};
else
rhop = rho;
end
if nargin > 13
beta = varargin{13};
nu = varargin{14};
end
% Defines dislocation in the fault plane system
U1 = cosd(rake).*slip;
U2 = sind(rake).*slip;
% Converts fault coordinates (E,N,DEPTH) relative to centroid
% into Okada's reference system (X,Y,D)
d = depth + sind(dip).*W/2; % d is fault's base edge as Okada's convention
ec = e + cosd(strike)*cosd(dip).*W/2;
nc = n - sind(strike)*cosd(dip).*W/2;
x = cosd(strike).*nc + sind(strike).*ec + L/2;
y = sind(strike).*nc - cosd(strike).*ec + cosd(dip).*W;
% Variable substitution (independent from xi and eta)
p = y.*cosd(dip) + d.*sind(dip);
q = y.*sind(dip) - d.*cosd(dip);
% Elevation changes dH (must be computed first) [equation (57) p. 7139]
dH = 1/(2*pi) * (U1*chinnery(@Sh,x,p,L,W,q,dip,nu) ... % strike-slip
+ U2 * chinnery(@Dh,x,p,L,W,q,dip,nu) ... % dip-slip
+ U3 * chinnery(@Th,x,p,L,W,q,dip,nu) ... % tensile fault
);
% Total gravity changes dG [equation (49) p. 7139]
dG = rho*G * (U1*chinnery(@Sg,x,p,L,W,q,dip,nu) ... % strike-slip
+ U2 * chinnery(@Dg,x,p,L,W,q,dip,nu) ... % dip-slip
+ U3 * chinnery(@Tg,x,p,L,W,q,dip,nu)) ... % tensile fault
+ (rhop-rho)*G*U3*chinnery(@Cg,x,p,L,W,q,dip,nu) ... % filling cavity
- beta*dH; % free-air effect
% Assigns output arguments
switch nargout
case 1
varargout = {dG};
case 2
varargout = {dG,dH};
otherwise
disp('Unvalid number of output arguments.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =================================================================
% --- Chinnery's notation [equation (23) p. 7138]
function u=chinnery(f,x,p,L,W,q,dip,nu)
u = feval(f,x,p,q,dip,nu) ...
- feval(f,x,p-W,q,dip,nu) ...
- feval(f,x-L,p,q,dip,nu) ...
+ feval(f,x-L,p-W,q,dip,nu);
% =================================================================
% Displacement subfunctions (Okada model regiven by Okubo)
% --- strike-slip displacement subfunction [equation (58) p. 7139]
function u=Sh(xi,eta,q,dip,nu)
R = sqrt(xi.^2 + eta.^2 + q.^2);
db = eta*sind(dip) - q*cosd(dip);
u = - db.*q./(R.*(R + eta)) ...
- q.*sind(dip)./(R + eta) ...
- I4(db,eta,q,dip,nu,R).*sind(dip);
% --- dip-slip displacement subfunction [equation (59) p. 7139]
function u=Dh(xi,eta,q,dip,nu)
R = sqrt(xi.^2 + eta.^2 + q.^2);
db = eta*sind(dip) - q*cosd(dip);
u = - db.*q./(R.*(R + xi)) ...
- sind(dip).*atan(xi.*eta./(q.*R)) ...
+ I5(xi,eta,q,dip,nu,R,db).*sind(dip).*cosd(dip);
% --- tensile fault displacement subfunction [equation (60) p. 7139]
function u=Th(xi,eta,q,dip,nu)
R = sqrt(xi.^2 + eta.^2 + q.^2);
db = eta*sind(dip) - q*cosd(dip);
u = (eta*cosd(dip) + q*sind(dip)).*q./(R.*(R + xi)) ...
+ cosd(dip).*(xi.*q./(R.*(R + eta)) ...
- atan(xi.*eta./(q.*R))) ...
- I5(xi,eta,q,dip,nu,R,db).*sind(dip).^2;
% --- I4 displacement subfunction [equations (61) and (63) p. 7139-7140]
function I=I4(db,eta,q,dip,nu,R)
if cosd(dip) == 0
I = -(1 - 2*nu) * q./(R + db);
else
I = (1 - 2*nu) * 1/cosd(dip) * (log(R + db) - sind(dip)*log(R + eta));
end
% --- I5 displacement subfunction [equations (62) and (64) p. 7140]
function I=I5(xi,eta,q,dip,nu,R,db)
%X = sqrt(xi.^2 + q.^2);
if cosd(dip) == 0
I = -(1 - 2*nu) * xi.*sind(dip)./(R + db);
else
I = (1 - 2*nu) * 2./cosd(dip) ...
.* atan((-q.*cosd(dip) + (1 + sind(dip)).*(R + eta))./(xi.*cosd(dip)));
% .* atan((eta.*(X + q.*cosd(dip)) + X.*(R + X).*sind(dip))./(xi.*(R + X).*cosd(dip)));
end
% =================================================================
% Gravity subfunctions
% --- strike-slip gravity subfunction [equation (52) p. 7139]
function u=Sg(xi,eta,q,dip,nu)
R = sqrt(xi.^2 + eta.^2 + q.^2);
u = - q*sind(dip)./R + (q.^2 * cosd(dip))./(R.*(R + eta));
% --- dip-slip gravity subfunction [equation (53) p. 7139]
function u=Dg(xi,eta,q,dip,nu)
R = sqrt(xi.^2 + eta.^2 + q.^2);
db = eta*sind(dip) - q*cosd(dip);
u = 2*I2(xi,eta,q,R)*sind(dip) - q.*db./(R.*(R + xi));
% --- tensile fault gravity subfunction [equation (54) p. 7139]
function u=Tg(xi,eta,q,dip,nu)
R = sqrt(xi.^2 + eta.^2 + q.^2);
yb = eta*cosd(dip) + q*sind(dip);
u = 2*I2(xi,eta,q,R)*cosd(dip) + q.*yb./(R.*(R+xi)) + q.*xi*cosd(dip)./(R.*(R+eta));
% --- cavity-filling gravity subfunction [equation (55) p. 7139]
function u=Cg(xi,eta,q,dip,nu)
R = sqrt(xi.^2 + eta.^2 + q.^2);
u = 2*I2(xi,eta,q,R)*cosd(dip) - sind(dip)*log(R + xi);
% --- I2 subfunction [equation (32) p. 7139]
function I=I2(xi,eta,q,R)
I = atan((R + xi + eta)./q);