modelling

pyrocko.modelling.okada

class AnalyticalSource(**kwargs)[source]

Base class for analytical source models.

name

str, optional, default: ''

time

time_float, optional, default: 0.0

Source origin time

vr

float, optional, default: 0.0

Rupture velocity [m/s]

class AnalyticalRectangularSource(**kwargs)[source]

Rectangular analytical source model.

Coordinates on the source plane are with respect to the origin point given by (lat, lon, east_shift, north_shift, depth).

strike

float, default: 0.0

Strike direction in [deg], measured clockwise from north

dip

float, default: 90.0

Dip angle in [deg], measured downward from horizontal

rake

float, default: 0.0

Rake angle in [deg], measured counter-clockwise from right-horizontal in on-plane view

al1

float, default: 0.0

Left edge source plane coordinate [m]

al2

float, default: 0.0

Right edge source plane coordinate [m]

aw1

float, default: 0.0

Lower edge source plane coordinate [m]

aw2

float, default: 0.0

Upper edge source plane coordinate [m]

slip

float, optional, default: 0.0

Slip on the rectangular source area [m]

class OkadaSource(**kwargs)[source]

Rectangular Okada source model.

opening

float, optional, default: 0.0

Opening of the plane in [m]

poisson

float, optional, default: 0.25

Poisson’s ratio \nu

lamb

float, optional

First Lame’ s parameter \lambda [Pa]

shearmod

float, optional, default: 32000000000.0

Shear modulus along the plane \mu [Pa]

poisson

Calculation of Poisson’ s ratio \nu (if not given).

The Poisson’ s ratio \nu can be calculated from the Lame parameters \lambda and \mu using \nu =
\frac{\lambda}{2(\lambda + \mu)} (e.g. Mueller 2007).

lamb

Calculation of first Lame’ s parameter (if not given).

Poisson’ s ratio \nu and shear modulus \mu must be available.

shearmod

Calculation of shear modulus \mu (if not given).

Poisson ratio’ s \nu must be available.

Important

We assume a perfect elastic solid with K=\frac{5}{3}\mu

Through \mu = \frac{3K(1-2\nu)}{2(1+\nu)} this leads to \mu = \frac{8(1+\nu)}{1-2\nu}

seismic_moment

Scalar Seismic moment M_0.

Code copied from Kite. It disregards the opening (as for now). We assume M_0 = mu A D

Important

We assume a perfect elastic solid with K=\frac{5}{3}\mu

Through \mu = \frac{3K(1-2\nu)}{2(1+\nu)} this leads to \mu = \frac{8(1+\nu)}{1-2\nu}

Returns:Seismic moment release
Return type:float
moment_magnitude

Moment magnitude M_\mathrm{w} from Seismic moment.

We assume M_\mathrm{w} = {\frac{2}{3}}\log_{10}(M_0) - 10.7

Returns:Moment magnitude
Return type:float
source_patch()[source]

Build source information array for okada_ext.okada input.

Returns:array of the source data as input for okada_ext.okada
Return type:numpy.ndarray, (1, 9)
source_disloc()[source]

Build source dislocation for okada_ext.okada input.

Returns:array of the source dislocation data as input for okada_ext.okada
Return type:numpy.ndarray, (1, 3)
discretize(nlength, nwidth, *args, **kwargs)[source]

Discretize fault into rectilinear grid of fault patches.

Fault orientation, slip and elastic parameters are passed to the sub-faults unchanged.

Parameters:
  • nlength (int) – Number of patches in strike direction
  • nwidth (int) – Number of patches in down-dip direction
Returns:

Discrete fault patches

Return type:

list of OkadaPatch objects

class OkadaPatch(parent=None, *args, **kwargs)[source]

Okada source with additional 2D indexes for bookkeeping.

ix

int

Relative index of the patch in x

iy

int

Relative index of the patch in y

make_okada_coefficient_matrix(source_patches_list, pure_shear=False, rotate_sdn=True, nthreads=1, variant='normal')[source]

Build coefficient matrix for given fault patches.

The boundary element method (BEM) for a discretized fault and the determination of the slip distribution from stress drop is based on the relation stress = coefmat \cdot displ. Here the coefficient matrix is built, based on the displacements from Okada’s solution and their partial derivatives.

Parameters:
  • source_patches_list (list of OkadaSource objects) – Source patches, to be used in BEM.
  • pure_shear (optional, bool) – If True, only shear forces are taken into account.
  • rotate_sdn (optional, bool) – If True, rotate to strike, dip, normal.
  • nthreads (optional, int) – Number of threads.
Returns:

Coefficient matrix for all source combinations.

Return type:

numpy.ndarray, (len(source_patches_list) * 3, len(source_patches_list) * 3)

invert_fault_dislocations_bem(stress_field, coef_mat=None, source_list=None, pure_shear=False, epsilon=None, nthreads=1, **kwargs)[source]

BEM least squares inversion to get fault dislocations given stress field.

Follows least squares inversion approach by Menke (1989) to calculate dislocations on a fault with several segments from a given stress field. The coefficient matrix connecting stresses and displacements of the fault patches can either be specified by the user (coef_mat) or it is calculated using the solution of Okada (1992) for a rectangular fault in a homogeneous half space (source_list).

Parameters:
  • stress_field (numpy.ndarray, shape (nsources, 3)) – Stress change [Pa] for each source patch (as stress_field[isource, icomponent] where isource indexes the source patch and icomponent indexes component, ordered (strike, dip, tensile).
  • coef_mat (optional, numpy.ndarray, shape (len(source_list) * 3, len(source_list) * 3)) – Coefficient matrix connecting source patch dislocations and the stress field.
  • source_list (optional, list of OkadaSource objects) – list of all source patches to be used for BEM.
  • epsilon (optional, float) – If given, values in coef_mat smaller than epsilon are set to zero.
  • nthreads (int) – Number of threads allowed.
Returns:

Inverted displacements as displacements[isource, icomponent] where isource indexes the source patch and icomponent indexes component, ordered (strike, dip, tensile).

Return type:

numpy.ndarray, (nsources, 3)

pyrocko.modelling.cracksol

class GriffithCrack(**kwargs)[source]

Analytical Griffith crack model

width

float, default: 1.0

Width equals to 2 \cdot a

poisson

float, default: 0.25

Poisson ratio \nu

shearmod

float, default: 1000000000.0

Shear modulus \mu [Pa]

stressdrop

numpy.ndarray (pyrocko.guts_array.Array), default: array([0., 0., 0.])

Stress drop array:[\sigma_{3,r} - \sigma_{3,c}, \sigma_{2,r} - \sigma_{2,c}, \sigma_{1,r} - \sigma_{1,c}] = [\Delta stress_{strike}, \Delta stress_{dip}, \Delta stress_{tensile}]

a

Half width of the crack

disloc_infinite2d(x_obs)[source]

Calculation of dislocation at crack surface along x2 axis

Follows equations by Pollard and Segall (1987) to calculate dislocations for an infinite 2D crack extended in x3 direction, opening in x1 direction and the crack extending in x2 direction.

Parameters:x_obs (numpy.ndarray, (N,)) – Observation point coordinates along x2-axis. If x_{obs} < -a or x_{obs} > a, output dislocations are zero
Returns:dislocations at each observation point in strike, dip and tensile direction.
Return type:numpy.ndarray, (N, 3)
disloc_circular(x_obs)[source]

Calculation of dislocation at crack surface along x2 axis

Follows equations by Pollard and Segall (1987) to calculate displacements for a circulat crack extended in x2 and x3 direction and opening in x1 direction.

Parameters:x_obs (numpy.ndarray, (N,)) – Observation point coordinates along axis through crack centre. If x_{obs} < -a or x_{obs} > a, output dislocations are zero
Returns:dislocations at each observation point in strike, dip and tensile direction.
Return type:numpy.ndarray, (N, 3)
displ_infinite2d(x1_obs, x2_obs)[source]

Calculation of displacement at crack surface along different axis

Follows equations by Pollard and Segall (1987) to calculate displacements for an infinite 2D crack extended in x3 direction, opening in x1 direction and the crack tip in x2 direction.

Parameters:
  • x1_obs (numpy.ndarray, (M,)) – Observation point coordinates along x1-axis. If x1_obs = 0., displacment is calculated along x2-axis
  • x2_obs (numpy.ndarray, (N,)) – Observation point coordinates along x2-axis. If x2_obs = 0., displacment is calculated along x1-axis
Returns:

displacements at each observation point in strike, dip and tensile direction.

Return type:

numpy.ndarray, (M, 3) or (N, 3)