cake

Classical seismic ray theory for layered earth models (layer cake models).

This module can be used to e.g. calculate arrival times, ray paths, reflection and transmission coefficients, take-off and incidence angles and geometrical spreading factors for arbitrary seismic phases. Computations are done for a spherical earth, even though the module name may suggests something flat.

The main classes defined in this module are:

• Material - Defines an isotropic elastic material.

• PhaseDef - Defines a seismic phase arrival / wave propagation

  history.

• Leg - Continuous propagation in a PhaseDef.

• Knee - Conversion/reflection in a PhaseDef.

• LayeredModel - Representation of a layer cake model.

• Layer - A layer in a LayeredModel.

  • HomogeneousLayer - A homogeneous Layer.

  • GradientLayer - A gradient Layer.

• Discontinuity - A discontinuity in a LayeredModel.

  • Interface - A Discontinuity between two Layer instances.

  • Surface - The surface Discontinuity on top of a LayeredModel.

• RayPath - A fan of rays running through a common sequence of layers / interfaces.

• Ray - A specific ray with a specific (ray parameter, distance, arrival time) choice.

• RayElement - An element of a RayPath.

  • Straight - A ray segment representing propagation through one Layer.

  • Kink - An interaction of a ray with a Discontinuity.

exception CakeError

exception InvalidArguments

class Material(vp=None, vs=None, rho=2600.0, qp=None, qs=None, poisson=None, lame=None, qk=None, qmu=None, burgers=None)

Isotropic elastic material.

Parameters:

• vp – P-wave velocity [m/s]

• vs – S-wave velocity [m/s]

• rho – density [kg/m^3]

• qp – P-wave attenuation Qp

• qs – S-wave attenuation Qs

• poisson – Poisson ratio

• lame – tuple with Lame parameter lambda and shear modulus [Pa]

• qk – bulk attenuation Qk

• qmu – shear attenuation Qmu

• burgers (tuple) – Burgers rheology paramerters as tuple. transient viscosity [Pa], <= 0 means infinite value, steady-state viscosity [Pa] and alpha, the ratio between the effective and unreleaxed shear modulus, mu1/(mu1 + mu2).

If no velocities and no lame parameters are given, standard crustal values of vp = 5800 m/s and vs = 3200 m/s are used. If no Q values are given, standard crustal values of qp = 1456 and qs = 600 are used. If no Burgers material parameters are given, transient and steady-state viscosities are 0 and alpha=1.

Everything is in SI units (m/s, Pa, kg/m^3) unless explicitly stated.

The main material properties are considered independant and are accessible as attributes (it is allowed to assign to these):

vp, vs, rho, qp, qs

Other material properties are considered dependant and can be queried by instance methods.

astuple()

Get independant material properties as a tuple.

Returns a tuple with (vp, vs, rho, qp, qs).

lame()

Get Lame’s parameter lambda and shear modulus.

lame_lambda()

Get Lame’s parameter lambda.

Returned units are [Pa].

shear_modulus()

Get shear modulus.

Returned units are [Pa].

poisson()

Get Poisson’s ratio.

bulk()

Get bulk modulus.

youngs()

Get Young’s modulus.

vp_vs_ratio()

Get vp/vs ratio.

qmu()

Get shear attenuation coefficient Qmu.

qk()

Get bulk attenuation coefficient Qk.

burgers()

Get Burger parameters.

rayleigh()

Get Rayleigh velocity assuming a homogenous halfspace.

Returned units are [m/s].

describe()

Get a readable listing of the material properties.

class Leg(departure=None, mode=None)

Represents a continuous piece of wave propagation in a phase definition.

Attributes:

To be considered as read-only.

departure

One of the constants UP or DOWN indicating upward or downward departure.

mode

One of the constants P or S, indicating the propagation mode.

depthmin

None, a number (a depth in [m]) or a string (an interface name), minimum depth.

depthmax

None, a number (a depth in [m]) or a string (an interface name), maximum depth.

exception InvalidKneeDef

class Knee(*args)

Represents a change in wave propagation within a PhaseDef.

Attributes:

To be considered as read-only.

depth

Depth at which the conversion/reflection should happen. this can be a string or a number.

direction

One of the constants UP or DOWN to indicate the incoming direction.

in_mode

One of the constants P or S to indicate the type of mode of the incoming wave.

out_mode

One of the constants P or S to indicate the type of mode of the outgoing wave.

conversion

Boolean, whether there is a mode conversion involved.

reflection

Boolean, whether there is a reflection involved.

headwave

Boolean, whether there is headwave propagation involved.

matches(discontinuity, mode, direction)

Check whether it is relevant to a given combination of interface, propagation mode, and direction.

out_direction()

Get outgoing direction.

Returns one of the constants UP or DOWN.

class Head(*args)

exception UnknownClassicPhase(phasename)

exception PhaseDefParseError(definition, position, exception)

Exception raised when an error occures during parsing of a phase definition string.

class PhaseDef(definition=None, classicname=None)

Definition of a seismic phase arrival, based on ray propagation path.

Parameters:

definition – string representation of the phase in Cake’s phase syntax

Seismic phases are conventionally named e.g. P, Pn, PP, PcP, etc. In Cake, a slightly different terminology is adapted, which allows to specify arbitrary conversion/reflection histories for seismic ray paths. The conventions used here are inspired by those used in the TauP toolkit, but are not completely compatible with those.

The definition of a seismic ray propagation path in Cake’s phase syntax is a string consisting of an alternating sequence of legs and knees.

A leg represents seismic wave propagation without any conversions, encountering only super-critical reflections. Legs are denoted by P, p, S, or s. The capital letters are used when the take-off of the leg is in downward direction, while the lower case letters indicate a take-off in upward direction.

A knee is an interaction with an interface. It can be a mode conversion, a reflection, or propagation as a headwave or diffracted wave.

• conversion is simply denoted as: (INTERFACE) or DEPTH

• upperside reflection: v(INTERFACE) or vDEPTH

• underside reflection: ^(INTERFACE) or ^DEPTH

• normal kind headwave or diffracted wave: v_(INTERFACE) or v_DEPTH

The interface may be given by name or by depth: INTERFACE is the name of an interface defined in the model, DEPTH is the depth of an interface in [km] (the interface closest to that depth is chosen). If two legs appear consecutively without an explicit knee, surface interaction is assumed.

The phase definition may end with a backslash \, to indicate that the ray should arrive at the receiver from above instead of from below. It is possible to restrict the maximum and minimum depth of a leg by appending <(INTERFACE) or <DEPTH or >(INTERFACE) or >DEPTH after the leg character, respectively.

Examples:

• P - like the classical P, but includes PKP, PKIKP, Pg

• P<(moho) - like classical Pg, but must leave source downwards

• pP - leaves source upward, reflects at surface, then travels as P

• P(moho)s - conversion from P to S at the Moho on upgoing path

• P(moho)S - conversion from P to S at the Moho on downgoing path

• Pv12p - P with reflection at 12 km deep interface (or the interface closest to that)

• Pv_(moho)p - classical Pn

• Pv_(cmb)p - classical Pdiff

• P^(conrad)P - underside reflection of P at the Conrad discontinuity

Usage:

>>> from pyrocko.cake import PhaseDef
# must escape the backslash
>>> my_crazy_phase = PhaseDef('pPv(moho)sP\\')
>>> print my_crazy_phase
Phase definition "pPv(moho)sP":
 - P mode propagation, departing upward
 - surface reflection
 - P mode propagation, departing downward
 - upperside reflection with conversion from P to S at moho
 - S mode propagation, departing upward
 - surface reflection with conversion from S to P
 - P mode propagation, departing downward
 - arriving at target from above

Note

1. These conventions might be extended in a way to allow to fix wave propagation to SH mode, possibly by specifying SH, or a single character (e.g. H) instead of S. This would be benificial for the selection of conversion and reflection coefficients, which currently only deal with the P-SV case.

static classic(phasename)

Get phase definitions based on classic phase name.

Parameters:

phasename – classic name of a phase

Returns:

list of PhaseDef objects

This returns a list of PhaseDef objects, because some classic phases (like e.g. Pg) can only be represented by two Cake style PhaseDef objects (one with downgoing and one with upgoing first leg).

first_leg()

Get the first leg in phase definition.

last_leg()

Get the last leg in phase definition.

legs()

Iterate over the continuous pieces of wave propagation (legs) defined within this phase definition.

knees()

Iterate over conversions and reflections (knees) defined within this phase definition.

definition()

Get original definition of the phase.

given_name()

Get entered classic name if any, or original definition of the phase.

used_repr()

Translate into textual representation (cake phase syntax).

copy()

Get a deep copy of it.

psv_surface_ind(in_mode, out_mode)

Get indices to select the appropriate element from scatter matrix for free surface.

psv_surface(material, p, energy=False)

Scatter matrix for free surface reflection/conversions.

Parameters:

• material – material, object of type Material

• p – flat ray parameter [s/m]

• energy – bool, when True energy normalized coefficients are returned

Returns:

Scatter matrix

The scatter matrix is ordered as follows:

[[PP, PS],
 [SP, SS]]

The formulas given in Aki & Richards are used.

psv_solid_ind(in_direction, out_direction, in_mode, out_mode)

Get indices to select the appropriate element from scatter matrix for solid-solid interface.

psv_solid(material1, material2, p, energy=False)

Scatter matrix for solid-solid interface.

Parameters:

• material1 – material above, object of type Material

• material2 – material below, object of type Material

• p – flat ray parameter [s/m]

• energy – bool, when True energy normalized coefficients are returned

Returns:

Scatter matrix

The scatter matrix is ordered as follows:

[[P1P1, S1P1, P2P1, S2P1],
 [P1S1, S1S1, P2S1, S2S1],
 [P1P2, S1P2, P2P2, S2P2],
 [P1S2, S1S2, P2S2, S2S2]]

The formulas given in Aki & Richards are used.

exception BadPotIntCoefs

exception PathFailed

exception SurfaceReached

exception BottomReached

exception MaxDepthReached

exception MinDepthReached

exception Trapped

exception NotPhaseConform

exception CannotPropagate(direction, ilayer)

class Layer(ztop, zbot, name=None)

Representation of a layer in a layered earth model.

Parameters:

• ztop – depth of top of layer

• zbot – depth of bottom of layer

• name – name of layer (optional)

Subclasses are: HomogeneousLayer and GradientLayer.

potint_coefs(mode)

Get coefficients for potential interpolation.

Parameters:

mode – mode of wave propagation, P or S

Returns:

coefficients (a, b)

contains(z)

Tolerantly check if a given depth is within the layer (including boundaries).

inner(z)

Tolerantly check if a given depth is within the layer (not including boundaries).

at_bottom(z)

Tolerantly check if given depth is at the bottom of the layer.

at_top(z)

Tolerantly check if given depth is at the top of the layer.

pflat_top(p)

Convert spherical ray parameter to local flat ray parameter for top of layer.

pflat_bottom(p)

Convert spherical ray parameter to local flat ray parameter for bottom of layer.

pflat(p, z)

Convert spherical ray parameter to local flat ray parameter for given depth.

xt_potint(p, mode, zpart=None)

Get travel time and distance for for traversal with given mode and ray parameter.

Parameters:

• p – ray parameter (spherical)

• mode – mode of propagation (P or S)

• zpart – if given, tuple with two depths to restrict computation to a part of the layer

This implementation uses analytic formulas valid for a spherical earth in the case where the velocity c within the layer is given by potential interpolation of the form

c(z) = a*z^b

test(p, mode, z)

Check if wave mode can exist for given ray parameter at given depth within the layer.

zturn_potint(p, mode)

Get turning depth for given ray parameter and propagation mode.

propagate(p, mode, direction)

Propagate ray through layer.

Parameters:

• p – ray parameter

• mode – propagation mode

• direction – in direction (UP or DOWN

resize(depth_min=None, depth_max=None)

Change layer thinkness and interpolate material if required.

exception DoesNotTurn

class HomogeneousLayer(ztop, zbot, m, name=None)

Representation of a homogeneous layer in a layered earth model.

Base class: Layer.

class GradientLayer(ztop, zbot, mtop, mbot, name=None)

Representation of a gradient layer in a layered earth model.

Base class: Layer.

class Discontinuity(z, name=None)

Base class for discontinuities in layered earth model.

Subclasses are: Interface and Surface.

class Interface(z, mabove, mbelow, name=None)

Representation of an interface in a layered earth model.

Base class: Discontinuity.

class Surface(z, mbelow)

Representation of the surface discontinuity in a layered earth model.

Base class: Discontinuity.

class RayElement

An element of a RayPath.

class Straight(direction_in, direction_out, mode, layer)

A ray segment representing wave propagation through one Layer.

class HeadwaveStraight(direction_in, direction_out, mode, interface)

class Kink(in_direction, out_direction, in_mode, out_mode, discontinuity)

An interaction of a ray with a Discontinuity.

exception PRangeNotSet

class RayPath(phase)

Representation of a fan of rays running through a common sequence of layers / interfaces.

copy()

Get a copy of it.

endgaps(zstart, zstop)

Get information needed for end point adjustments.

used_phase(p=None, eps=1.0)

Calculate phase definition from ray path.

pmax()

Get maximum valid ray parameter.

pmin()

Get minimum valid ray parameter.

xmin()

Get minimal distance.

xmax()

Get maximal distance.

kinks()

Iterate over propagation mode changes (reflections/transmissions).

straights()

Iterate over ray segments.

first_straight()

Get first ray segment.

last_straight()

Get last ray segment.

efficiency(p)

Get product of all conversion/reflection coefficients encountered on path.

spreading(p, endgaps)

Get geometrical spreading factor.

xt_endgaps(p, endgaps, which='both')

Get amount of distance/traveltime to be subtracted at the generic ray path’s ends.

xt_endgaps_ptest(p, endgaps)

Check if ray parameter is valid at source and receiver.

xt(p, endgaps)

Calculate distance and traveltime for given ray parameter.

xt_limits(p)

Calculate limits of distance and traveltime for given ray parameter.

iter_zxt(p)

Iterate over (depth, distance, traveltime) at each layer interface on ray path.

zxt_path_subdivided(p, endgaps, points_per_straight=20, x_for_headwave=None)

Get geometrical representation of ray path.

interpolate_x2pt_linear(x, endgaps)

Get approximate ray parameter and traveltime for distance.

critical_pstart(endgaps)

Get critical ray parameter for source depth choice.

critical_pstop(endgaps)

Get critical ray parameter for receiver depth choice.

ranges(endgaps)

Get valid ranges of ray parameter, distance, and traveltime.

describe(endgaps=None, as_degrees=False)

Get textual representation.

exception RefineFailed

class Ray(path, p, x, t, endgaps, draft_pxt)

Representation of a ray with a specific (path, ray parameter, distance, arrival time) choice.

Attributes:

path

RayPath object containing complete propagation history.

p

Ray parameter (spherical) [s/rad]

x

Radial distance [deg]

t

Traveltime [s]

endgaps

Needed for source/receiver depth adjustments in many RayPath methods.

given_phase()

Get phase definition which was used to create the ray.

Returns:

PhaseDef object

used_phase()

Compute phase definition from propagation path.

Returns:

PhaseDef object

takeoff_angle()

Get takeoff angle of ray.

The angle is returned in [degrees].

incidence_angle()

Get incidence angle of ray.

The angle is returned in [degrees].

efficiency()

Get conversion/reflection efficiency of the ray.

A value between 0 and 1 is returned, reflecting the relative amount of energy which is transmitted along the ray and not lost by reflections or conversions.

spreading()

Get geometrical spreading factor.

zxt_path_subdivided(points_per_straight=20)

Get geometrical representation of ray path.

Three arrays (depth, distance, time) with points on the ray’s path of propagation are returned. The number of points which are used in each ray segment (passage through one layer) may be controlled by the points_per_straight parameter.

exception DiscontinuityNotFound(depth_or_name)

exception LayeredModelError

class LayeredModel

Representation of a layer cake model.

There are several ways to initialize an instance of this class.

1. Use the module function load_model() to read a model from a file.

2. Create an empty model with the default constructor and append layers and discontinuities with the append() method (from top to bottom).

3. Use the constructor LayeredModel.from_scanlines(), to automatically create the Layer and Discontinuity objects from a given velocity profile.

An earth model is represented by as stack of Layer and Discontinuity objects. The method arrivals() returns Ray objects which may be e.g. queried for arrival times of specific phases. Each ray is associated with a RayPath object. Ray objects share common ray paths if they have the same conversion/reflection/propagation history. Creating the ray path objects is relatively expensive (this is done in gather_paths()), but they are cached for reuse in successive invocations.

copy_with_elevation(elevation)

Get copy of the model with surface layer stretched to given elevation.

Parameters:

elevation – new surface elevation in [m]

Elevation is positiv upward, contrary to the layered models downward z axis.

append(element)

Add a layer or discontinuity at bottom of model.

Parameters:

element – object of subclass of Layer or Discontinuity.

elements(direction=4)

Iterate over all elements of the model.

Parameters:

direction – direction of traversal DOWN or UP.

Objects derived from the Discontinuity and Layer classes are yielded.

layers(direction=4)

Iterate over all layers of model.

Parameters:

direction – direction of traversal DOWN or UP.

Objects derived from the Layer class are yielded.

layer(z, direction=4)

Get layer for given depth.

Parameters:

• z – depth [m]

• direction – direction of traversal DOWN or UP.

Returns first layer which touches depth z (tolerant at boundaries).

material(z, direction=4)

Get material at given depth.

Parameters:

• z – depth [m]

• direction – direction of traversal DOWN or UP

Returns:

object of type Material

If given depth z happens to be at an interface, the material of the first layer with respect to the the traversal ordering is returned.

discontinuities()

Iterate over all discontinuities of the model.

discontinuity(name_or_z)

Get discontinuity by name or depth.

Parameters:

name_or_z – name of discontinuity or depth [m] as float value

adapt_phase(phase)

Adapt a phase definition for use with this model.

This returns a copy of the phase definition, where named discontinuities are replaced with the actual depth of these, as defined in the model.

path(p, phase, layer_start, layer_stop)

Get ray path for given combination of ray parameter, phase definition, source and receiver layers.

Parameters:

• p – ray parameter (spherical) [s/rad]

• phase – phase definition (PhaseDef object)

• layer_start – layer with source

• layer_stop – layer with receiver

Returns:

RayPath object

If it is not possible to find a solution, an exception of type NotPhaseConform, MinDepthReached, MaxDepthReached, CannotPropagate, BottomReached or SurfaceReached is raised.

gather_paths(phases=PhaseDef('P'), zstart=0.0, zstop=0.0)

Get all possible ray paths for given source and receiver depths for one or more phase definitions.

Parameters:

• phases – a PhaseDef object or a list of such objects. Comma-separated strings and lists of such strings are also accepted and are converted to PhaseDef objects for convenience.

• zstart – source depth [m]

• zstop – receiver depth [m]

Returns:

a list of RayPath objects

Results of this method are cached internally. Cached results are returned, when a given combination of source layer, receiver layer and phase definition has been used before.

arrivals(distances=[], phases=PhaseDef('P'), zstart=0.0, zstop=0.0, refine=True)

Compute rays and traveltimes for given distances.

Parameters:

• distances – list or array of distances [deg]

• phases – a PhaseDef object or a list of such objects. Comma-separated strings and lists of such strings are also accepted and are converted to PhaseDef objects for convenience.

• zstart – source depth [m]

• zstop – receiver depth [m]

• refine – bool flag, whether to use bisectioning to improve (p, x, t) estimated from interpolation

Returns:

a list of Ray objects, sorted by (distance, arrival time)

classmethod from_scanlines(producer)

Create layer cake model from sequence of materials at depths.

Parameters:

producer – iterable yielding (depth, material, name) tuples

Creates a new LayeredModel object and uses its append() method to add layers and discontinuities as needed.

profile(get)

Get parameter profile along depth of the earthmodel.

Parameters:

get (string) – property to be queried ( 'vp', 'vs', 'rho', 'qp', or 'qs', or 'z')

min(get='vp')

Find minimum value of a material property or depth.

Parameters:

get – property to be queried ( 'vp', 'vs', 'rho', 'qp', or 'qs', or 'z')

max(get='vp')

Find maximum value of a material property or depth.

Parameters:

get – property to be queried ( 'vp', 'vs', 'rho', 'qp', 'qs', or 'z')

simplify(max_rel_error=0.001)

Get representation of model with lower resolution.

Returns an approximation of the model. All discontinuities are kept, but layer stacks with continuous model parameters are represented, if possible, by a lower number of layers. Piecewise linear functions are fitted against the original model parameter’s piecewise linear functions. Successively larger numbers of layers are tried, until the difference to the original model is below max_rel_error. The difference is measured as the RMS error of the fit normalized by the mean of the input (i.e. the fitted curves should deviate, on average, less than 0.1% from the input curves if max_rel_error = 0.001).

extract(depth_min=None, depth_max=None)

Extract LayeredModel from LayeredModel.

Parameters:

• depth_min – depth of upper cut or name of Interface

• depth_max – depth of lower cut or name of Interface

Interpolates a GradientLayer at depth_min and/or depth_max.

perturb(rstate=None, keep_vp_vs=False, zmax=None, **kwargs)

Create a perturbed variant of the earth model.

Randomly change the thickness and material parameters of the earth model from a uniform distribution.

Parameters:

• kwargs (dict) – Maximum change fraction (e.g. 0.1) of the parameters. Name the parameter, prefixed by p. Supported parameters are ph, pvp, pvs, prho, pqs, pqp.

• rstate (numpy.random.RandomState, optional) – Random state to draw from, defaults to None

• keep_vp_vs (bool, optional) – Keep the Vp/Vs ratio, defaults to False

Returns:

A new, perturbed earth model

Return type:

LayeredModel

perturbed_model = model.perturb(ph=.1, pvp=.05, prho=.1)

read_hyposat_model(fn)

Reader for HYPOSAT earth model files.

To be used as producer in LayeredModel.from_scanlines().

Interface names are translated as follows: 'MOHO' -> 'moho', 'CONR' -> 'conrad'

read_nd_model(fn)

Reader for TauP style ‘.nd’ (named discontinuity) files.

To be used as producer in LayeredModel.from_scanlines().

Interface names are translated as follows: 'mantle' -> 'moho', 'outer-core' -> 'cmb', 'inner-core' -> 'icb'.

The format has been modified to include Burgers materials parameters in columns 7 (burger_eta1), 8 (burger_eta2) and 9. eta(3).

load_model(fn='ak135-f-continental.m', format='nd', crust2_profile=None)

Load layered earth model from file.

Parameters:

• fn – filename

• format – format

• crust2_profile – (lat, lon) or pyrocko.crust2x2.Crust2Profile object, merge model with crustal profile. If fn is forced to be None only the converted CRUST2.0 profile is returned.

Returns:

object of type LayeredModel

The following formats are currently supported:

format	description
'nd'	‘named discontinuity’ format used by the TauP programs
'hyposat'	format used by the HYPOSAT location program

The naming of interfaces is translated from the file format’s native naming to Cake’s own convention (See read_nd_model() and read_hyposat_model() for details). Cake likes the following internal names: 'conrad', 'moho', 'cmb' (core-mantle boundary), 'icb' (inner core boundary).

castagna_vs_to_vp(vs)

Calculate vp from vs using Castagna’s relation.

Castagna’s relation (the mudrock line) is an empirical relation for vp/vs for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al., 1985]

vp = 1.16 * vs + 1360 [m/s]

Parameters:

vs – S-wave velocity [m/s]

Returns:

P-wave velocity [m/s]

castagna_vp_to_vs(vp)

Calculate vp from vs using Castagna’s relation.

Castagna’s relation (the mudrock line) is an empirical relation for vp/vs for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al., 1985]

vp = 1.16 * vs + 1360 [m/s]

Parameters:

vp – P-wave velocity [m/s]

Returns:

S-wave velocity [m/s]

filled(v, *args, **kwargs)

Create NumPy array filled with given value.

This works like numpy.ones() but initializes the array with v instead of ones.